Abstract
A new approach is presented for the numerical evaluation of arbitrary singular domain integrals. In this method, singular domain integrals are transformed into a boundary integral and a radial integral which contains singularities by using the radial integration method. The analytical elimination of singularities condensed in the radial integral formulas can be accomplished by expressing the nonsingular part of the integration kernels as a series of cubic B-spline basis functions of the distance r and using the intrinsic features of the radial integral. In the proposed method, singularities involved in the domain integrals are explicitly transformed to the boundary integrals, so no singularities exist at internal points. A few numerical examples are provided to verify the correctness and robustness of the presented method.
1. Introduction
When the boundary element method (BEM) [1, 2] is used to solve complicated engineering problems [3, 4], such as compressible potential flows, heat conduction with heat generation, seepage problems with sources, electromagnetic field problems with electric charge, and elastostatics with body forces, a large number of regular and/or singular domain integrals may appear in the basic integral equations [2]. The accurate evaluation of general regular and/or singular domain integrals is very important and is still a crucial area of research in the BEM [5–11].
To compute the domain integrals, the frequently used approach is to apply some special schemes to transform domain integrals into boundary integrals, thus avoiding the discretization of the internal domain. The most extensively used technique is the dual reciprocity method (DRM) [12], which transforms the domain integrals to the boundary by approximating the given body force effect quantities as a series of prescribed basis functions and then employing particular solutions corresponding to these basis functions. As an extension of the idea of DRM, the multiple reciprocity method (MRM) [13, 14] was presented by Nowak and Brebbia [13] and the triple-reciprocity method was developed by Ochiai and Sekiya [15]. The radial integration method (RIM) presented by Gao [16–18] is another powerful method for transforming domain integrals into boundary integrals based on pure mathematical treatments, without using the Laplace operator and particular solutions of the problem.
Recently, an efficient approach was presented by Gao and Peng [7] for evaluating arbitrarily high-order singular domain integrals in a unified way based on the radial integration method. In the method, the radial integral can be integrated analytically by expressing the nonsingular part of integrand as polynomials of the global distance . Therefore, this approach can deal with any type of regular, weakly, and high-order singular domain integrals. However, the convergence performance of the computational result is not very good and is volatility with respect to the order of polynomials for some complicated function over complicated geometry due to the characteristics of high-order polynomial.
In this paper, the nonsingular part of the integrand involved in the domain integral is expressed as a series of third-degree B-spline basis functions [19] of the distance . This approach has two advantages. One is that the coefficients calculation becomes much simpler than that in [7] because of the local support and the endpoint interpolatory properties of B-spline basis functions with open knot vectors [20]. The other is that the integration results are very stable, less dependent on the number of order of polynomials in [7]. A number of singular domain integrals are given to verify the correctness and validity of the presented method.
2. Singular Domain Integrals
In this paper, the following singular domain integral is considered: where is bounded everywhere, being the nonsingular part of the integrand, is the domain of the problem to be considered, and is the distance between the source point and the field point ; that is,
The level of the singularity of domain integrals in (1) depends on the value of . For a finite integration domain , the regular and weakly singular domain integrals [7] always exist (i.e., results are finite). However, for strongly, hyper- and supersingular domain integrals [7], integrals only exist under some conditions, depending heavily on the characteristics of the function [18]. Usually, integrals stemming from real physical problems should exist. In this paper, the singular integrals are evaluated in the Cauchy principal value sense; that is, where is a small circular domain for 2D problems or a small spherical domain for 3D problems with the radius around the source point .
3. The Radial Integration Method [19, 20]
According to the radial integration method, the singular domain integral represented by (1) can be transformed into the equivalent boundary integral as follows: where In (4) and (5), for two-dimensional (2D) problems and for three-dimensional (3D) problems; denotes the coordinates at the boundary point ; is the distance from the source point to the boundary point ; is the boundary of the integration domain ; and is the unit outward normal vector to the boundary .
In order to evaluate the radial integral (5), the coordinate in the integrand needs to be expressed in terms of the integration variable ; that is, where represents the th component of Cartesian coordinates at the field point and is defined as follows:
From (4) and (5), we can see that the key task for the transformation of singular domain integral to the boundary integral is to accurately evaluate the radial integral (5). For simple functions, the radial integral (5) can be integrated analytically. However, for some complicated functions, analytically evaluating the radial integral is not available. In particular, for high-order singular integrals, the numerical evaluation of radial integrals may have infinite value problems. In this paper, a new approach will be presented for evaluating the radial integral (5) analytically in a unified way. This will be described in detail in the next section.
4. Evaluation of Singular Domain Integrals
In the section, the analytical evaluation of radial integral (5) is carried out by expressing the nonsingular part as a series of third-degree B-spline basis functions [19] of the distance .
4.1. Introduction of B-Spline Basis Functions
Let be a nondecreasing sequence of real numbers; that is, , . are called knots, and is the knot vector. The th B-spline basis function of -degree, denoted by , is defined as [19, 21] where ; . If knots are equally spaced in the parametric space, they are said to be uniform [19, 20]. A knot vector is said to be open if its first and last knots appear times. An initial example of the results of applying (8a) and (8b) to a uniform knot vector [22] is presented in Figure 1.
(a)
(b)
(c)
(d)
Important properties of B-spline basis functions are as follows [19, 20].(1)They constitute a partition of unity; that is, , .(2)The support of each is compact and contained in the interval .(3)Each basis function is nonnegative; that is, , .(4)For an open, nonuniform knot vector, the basis functions are interpolatory at the ends of the interval.
In Sections 4.2 and 4.3, it will be seen that properties (2) and (3) of B-spline basis functions can bring simplicity to our calculation.
4.2. Express Nonsingular Part of the Integration Kernel as a Series of Third-Degree B-Spline Basis Functions
It is assumed that the radial integration equations (4) and (5) are in the Cauchy principle value sense, such that (4) and (5) are written as where
In order to analytically evaluate the singular radial integral (10), the nonsingular part of the integrand is expressed as a series of third-degree B-spline basis functions of the distance ; that is, where are the third-degree B-spline basis functions and are the coefficients to be determined. Once knot vector is defined, the third-degree B-spline basis functions can be determined. In this paper, the following thirteen types of knot vectors are used, corresponding to , respectively; where .
The formulation of third-degree B-spline basis functions can be deduced from (8a) and (8b) as follows:
It should be noted that the open knot vectors listed in (12) are -continuous [19, 20] except at the endpoints of the interval, which makes the calculation more accurate. Moreover, the basis functions formed from open knot vectors are interpolatory at the ends of the parametric space interval ; that is,
The other coefficients are determined by collocating distances () over the integration region , in which where , .
So the other coefficients can be solved using where is a square matrix of order ; that is, and and are the following vectors:
Due to the local support property [19, 20] of B-spline basis functions, the coefficient matrix is banded and (16) is easy to be solved.
4.3. Evaluation of Singular Domain Integrals
Once the coefficients are obtained from the above equations, substituting (11) into (10) yields where
Substituting (13) into (20) yields From (13), it can be seen that is cubic polynomial functions of independent variable . Therefore, can be expressed as where can be determined by the polynomial arithmetic of (13).
Accordingly, appearing in (21) can be expressed as where
For a physical problem, the infinite terms in (24) should be cancelled out by free terms. That is, if the integral exists, after considering the contributions of all elements around points with the same size of , the last term in (24) must be zero [2, 18]. Based on this assumption, the final integral results can be expressed as
Substituting (25) into (23), the result into (21), and finally applying (19) yield
Substituting (26) into (9), it follows that
Now the singular domain integral has been transformed into boundary integral. As a result, no singularity exists for internal points and the boundary integral in (27) can be evaluated using Gaussian quadrature [2, 18]. For each Gauss point used in the Gaussian quadrature, the coefficient and quantity are computed using (16) and (25), respectively. However, for boundary points, singularity still exists and more treatment should be performed (see [22, 23]).
5. Numerical Examples
In order to verify the presented approach in this paper, some 2D and 3D singular domain integrals are analyzed in this section.
5.1. 2D Highly Singular Domain Integrals over a Rectangular Region
The first example is aimed at verifying the correctness of the presented formulations. The following domain integral is analyzed:
The integration domain considered is a rectangular region with and (Figure 2). To evaluate the domain integral, the boundary of the rectangle is discretized into 12 equally spaced three-noded quadratic elements with 24 boundary nodes. Two source points (located at nodes 25 and 26 in Figure 2) are considered, the coordinates of which are shown in Table 1. Table 1 lists the computed results for singularity order from 1 to 6. The computational results for the selected points using are denoted by RIM_B. For the integral, the analytical solution is easy to be obtained and also listed in Table 1.
(a)
(b)
From Table 1 it can be seen that the results (RIM_B) obtained by the derived formulations in the paper are in excellent agreement with the analytical solutions. It demonstrates that the method described in this paper can indeed evaluate high-order singular domain integrals accurately.
It is noted that the current computed results for , which are corresponding to knot vectors defined as (12) for , are the same. That is because the nonsingular part of the integrand is simple. It should also be noted that the more the Gauss points (here the number is 20) used are, the more accurate the results are. And the finer the mesh is, the more accurate the results are.
5.2. 2D Domain Integrals with Logarithm and Exponent over a Multiconnected Region
The second example is to demonstrate the capability of the method presented in this paper to evaluate complicated domain integrals over multiconnected regions. The domain considered is a quarter of fan-shaped sector with the inner and outer radiuses being 3 and 12, respectively, which is cut by two cavities (with inner radius of 6 and outer radius of 9, and with a 90/8-degree angle), as shown in Figure 3.
Two domain integrals are examined as follows:
To compute the singular domain integrals, the boundary of the domain is discretized into 114 three-noded quadratic elements with 228 boundary nodes (Figure 3). Three internal points A, B, and C are selected as the source points (see Figure 3), the coordinates of which are shown in Table 2. Computational results for these points using are listed in Table 3 and denoted by RIM_B. To verify the correctness of the computational results, the two domain integrals are also computed using the method presented in the literature [7]. The corresponding computational results are listed in the columns denoted by [7] in Table 3.
Table 3 shows that when , 1, or 2, the two sets of computational results are in good agreement, while, for high-order singular domain integrals, they are moderately close.
To examine the convergence performance of the computational result with respect to the number of B-spline basis functions, the domain integral was evaluated for different values of . Table 4 gives the computational results of the source point A for . It can be seen that good results can be obtained using .
5.3. 3D Singular Domain Integrals over a Hexagonal Prism
The third example is used to examine the performance of the proposed method in solving 3D high-order singular domain integrals. The following two singular domain integrals are considered:
The integration domain is a hexagonal prism with the height and side length being 1 and , respectively (Figure 4). To evaluate the integrals and , the surfaces of the hexagonal prism are discretized into 90 eight-noded quadratic elements with a total of 272 boundary nodes (Figure 4). Three internal points (273, 274, and 275) with the coordinates of , , and are selected as the source points. The computational results using are listed in Table 5 and denoted by RIM_B. The two domain integrals are also analyzed by the proposed method in [7] and the computational results are shown in Table 5 in the columns denoted by [7].
Comparing the numerical results in Table 5, it reveals that the two sets of computational results are in excellent agreement.
5.4. 3D Strongly Singular Domain Integral over a Quarter Circular Tube
The fourth example is concerned with the following strongly singular domain integral over a 3D curved circular tube with the radius (Figure 5). The curvature of the pipe is determined by a radius of . The following singular domain integral is considered: where in which the subscripts are for the dimensions in the 3D problem, is the Poisson ratio (here assumed to be 0.3), is determined by (7), and . can be found in the evaluation of the internal stresses in elastoplastic mechanics [2, 24]. Here only , , and are computed.
(a)
(b)
To evaluate these three integrals, the surfaces of the tube are discretized into 282 eight-noded quadratic elements with a total of 848 boundary nodes (Figure 5).
The problem was analyzed by Gao and Peng [7] and the results are used here as comparison. Table 7 gives the computed results using for three source points (849, 850, and 851) with coordinates listed in Table 6. It can be found that they are in excellent agreement from Table 7.
6. Conclusions
A new and robust method is described in this paper for evaluating arbitrary singular domain integrals based on RIM and by expressing the nonsingular part of the integration kernel as a series of cubic B-spline basis functions of the global distance . The distinct feature of the method is that regular, weakly, strongly, hyper-, and supersingular domain integrals are treated in a simple and unified manner. The provided numerical examples demonstrate that the presented method is accurate and robust.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors gratefully acknowledge the National Natural Science Foundation of China for financial support to this work under Grant NSFC no. 11172055.