- Email: [email protected]
A recursive application of the integral equation in the boundary element method
Engineering Analysis with Boundary Elements 35 (2011) 77–84
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A recursive application of the integral equation in the boundary element method C.F. Loeffler n ´ria, ES – CEP 29075-910, Brazil Mechanical Engineering Department, Federal University of Espı´rito Santo, Av. Fernando Ferrari, 514 Vito
a r t i c l e in fo
abstract
Article history: Received 5 August 2009 Accepted 26 May 2010 Available online 19 June 2010
This paper presents a recursive application of the governing integral equation aimed at improving the accuracy of numerical results of the boundary element method (BEM). Usually, only the results at internal domain points when using BEM are found using this approach, since the nodal boundary values have already been calculated. Here, it is shown that the same idea can be used to obtain better accuracy for the boundary results as well. Instead of locating the new source points inside the domain, they are positioned on the boundary, with different coordinates to the nodal points. The procedure is certainly general, but will be presented using as an example the two dimensional Laplace equation, for the sake of simplicity to point out the main concepts and numerical aspects of the method proposed, especially due to the determination of directional derivatives of the primal variable, which is part in hypersingular BEM theory. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Boundary element method Method of weighted residuals Scalar potential problems
1. Introduction Due to the mathematical modeling advances, improved computational capability and data storage, some traditional numerical techniques have been subjected to adaptations in order to achieve better performance. After procedures such as the relocation of discretization points, adaptive mesh refinement and adaptive shape functions, in some numerical methods an iterative solution scheme has been introduced in order to improve the numerical accuracy. Following this tendency, the recursive use of the boundary integral equation is here presented as an auxiliary resource to improve boundary element method (BEM) performance. It is not a completely iterative technique, but a simple procedure based on a common scheme to calculate internal values with BEM, where the boundary integral equation is used twice.
Laplace Equation can be written [1] as Z Z uðxÞ ¼ u ðx; XÞqðXÞ dG q ðx; XÞuðXÞ dG G
In Eq. (1), q(X) is the potential normal derivative. The essential (u ¼ u) and natural (u, i ni ¼ q) boundary conditions are defined on the boundary G ¼ Gu + Gq. The external unit normal vector at point X¼X(xi) is ni. For the usual direct BEM formulation, the auxiliary functions concern the solution of the related problem governed by a Poisson Equation [2], for which a unit singular source exists at X ¼ x for an infinite domain O (X). When the source point is positioned on the boundary, the integral equation for source points located on the boundary is given by Z Z cðxÞuðxÞ ¼ u ðx; XÞqðXÞ dGCPV q ðx; XÞuðXÞ dG ð2Þ G
2. Laplace boundary integral equation Considering an auxiliary functions un(x;X) and its normal derivative qn(x;X), where x is an internal domain point, and considering the basic mathematical tools of the theory of integral equations, the boundary integral equation associated to the n
Tel.: + 55 27 33352669. E-mail address: carlosloeffl[email protected]
0955-7997/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.05.012
ð1Þ
G
G
An ingenious mathematical procedure defines the c(x) value, which is equal to 0.5 for smooth boundaries [3]. When the source point is positioned on the boundary, un(x;X) and qn(x;X) functions present singularities. In spite of being discontinuous, un(x;X) is integrable along the boundary, but the integral of qn(x;X) exists in the Cauchy principal value (CPV) sense. According to the well-known BEM procedure [4], the boundary integral is discretized, generating a system of equations given in matrix notation as HYGQ ¼ 0
ð3Þ
78
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
3. Recursive boundary element procedure for potential variable The solution of BEM system of Eq. (3) allows the potential and its normal derivative for all N boundary nodes to be found. With regard to internal unknowns, it is possible to determinate easily their values putting the source point inside the domain and using the integral equation (1) again, but including the already calculated boundary variables. In general, the accuracy of numerical values of internal variables is usually better than the boundary ones. The mathematical reason for this behavior can be understood considering weighted residual fundamentals [2], since the approximated form of the boundary integral equation (1) can be expressed in a ^ that is ˆ and qE q, weighted residual form, taken as uEu Z Z Z ðx; XÞ dG ¼ 0 ^ ˆ ii u ðx; XÞ dO ¼ ˆ ðx; XÞ dG þ u, ðuuÞq ðqqÞu Gu
O
Gq
ð4Þ In Eq. (4), the left hand side integral represents the residual minimization in the domain and the right hand side integrals perform residual minimization of the essential and natural conditions on the boundary, respectively. Therefore, the recursive use of the integral equation to determinate numerical values for internal variables can be interpreted as a new process of weighting residuals, which expects a minimization of the errors previously performed, reaching best performance in the new collocation points, commonly located inside the domain. The same idea can be applied to recalculate boundary values with different positions on the boundary element. Thus, the recursive form of Eq. (4), considering previously calculated boundary values and putting new source point xi again on the boundary is given simply by Z Z N N X X i i Qke fk u ðxi ; XÞ dG Uke fk q ðxi ; XÞ dG ð5Þ cðx Þuðx Þ ¼ Ge
e¼1
e¼1
Ge
In Eq. (5), fk (X) are the shape functions, Qke and Uke are the nodal values on the boundary. The position of new source points must not coincide with former boundary nodes, because these points are just object of residual minimization; the recursive use of the integral equation is more effective using points situated between two adjacent boundary nodes.
4. Recursive boundary element procedure for potential spatial derivative Not only can the potential values be recalculated by the recursive procedure, but spatial potential derivatives related to xj (x) directions (or normal and tangential directions) can also be obtained by the same procedure. The mathematical foundation for determination of spatial derivatives follows same steps of the hyper-singular BEM formulation [4]. Let one consider, initially, an internal source point x. Taking the Cartesian derivatives of the integral equation (4) and using the Leibniz rule, it follows that Z Z u, i ðxÞ ¼ qi ðxÞ ¼ u, i ðx; XÞqðXÞ dG q, i ðx; XÞuðXÞ dG ð6Þ G
G
Considering smooth boundaries, s(x) is equal to 0.5 and the former integral involves only the potential derivatives in direction xi due to w(x) being null in this condition; however, for nonsmooth boundaries, the potential derivatives related to two different directions appear coupled in the same equation, as indicated by the permutation tensor e3ji. Concerning the boundary integrals, none of these is convergent in general; however, these two limits always exist when considered together, in CPV sense. Provided that the boundary values were previously calculated and considering smooth boundaries, Eq. (9) can be rewritten and applied to new xi points on the boundary: ( N Z X j j j j ½Uke jk uðx Þpi ðx ; XÞ dG sðx Þqi ðx Þ ¼ CPV e¼1
N X e¼1
Ge
Z Ge
) j Qke jk qi ðx ; XÞ dG
ð10Þ
According to the weighted residuals approach, some differences appear in the comparison between the former Eq. (10) and the equivalent sentence for potentials, given by Eq. (4). In spite of both integral equations making use the previous boundary values calculated on the functional nodes, the weighting functions are different, comprising different bases to residual orthogonalisation. This is an important factor for the difference in the numerical accuracy of recursive potential and derivative potential computations, such as occurs in the performance comparison between traditional singular and hyper-singular boundary techniques. However, the most important factor is the presence of u(x) subtracting the nodal variable Uke . This term, presented by convenience in Eq. (11), is responsible for the desired solution of the hyper-singular integral in a CPV sense, but its mathematical structure cannot be deduced from the basic weighted residual principles. N Z X j j j f ðx ; XÞ ¼ ½Uke fk uðx Þpi ðx ; XÞ dG ð11Þ e¼1
Ge
In many cases, due to the slow numerical convergence rate of this term, direct hyper-singular formulation presents worse performance than classical BEM formulation. The same behavior occurs with the recursive procedure for computation of the spatial derivative of the potential; thus, an effective strategy (presented next) must be used to guarantee good accuracy for the recursive technique, when potentials are not constant or vary in a different order to the shape functions used in the discretization procedure.
G
where: u, i ðx; XÞ ¼ qi ¼ ð1=2prÞr, i
analyzed, such as is done for the singular boundary integral. However, both kernels present singularities, which request more elaborated mathematical treatment [5,6]. Thus, locating the source point on the boundary, the hyper-singular integral equation for source points is given by Z sðxÞqi ðxÞ þ wðxÞe3ji qj ðxÞ ¼ CPV pi ðx; XÞ½uðXÞuðxÞ dG G Z qi ðx; XÞqðXÞ dG ð9Þ
2
q, i ðx; XÞ ¼ pi ¼ ð1=2pr Þ½2r, i r, j nj rj,i nj
5. Numerical aspects ð7Þ ð8Þ
To deduce the hyper-singular boundary integral equation, an augmented infinitesimal domain around the source point must be generated, and the effects created by this additional boundary
The existence of integrals in CPV sense for general cases to determine potential derivatives is guaranteed only when the two integrals on right hand side of Eq. (10) are analyzed together. Although the theory of hyper-singular integrals treatment is not the purpose of this paper, the main aspects are considered next for completeness.
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
First, the approach is shown when the integral contains the spatial derivative of the potential fundamental solution. Concerning general cases, initially it is important to divide this integral in two parts: the left hand side (LS) and the right hand side (RS) of the source point, using Hadamard’s methodology to treat the singularity [7], adding and subtracting a constant flux q(x), such as shown next: (Z Z 1 qLS ðXÞqLS ðxÞ LS LS r, i 9J9 dr qðXÞqi dG ¼ 2 r p sing LS Z qRS ðXÞqRS ðxÞ RS RS r, i 9J9 dr þ r RS LS
LS þ qLS ðxÞr, LS i 9J9 ½lnðL Þlim lnðrÞ r-0
RS
)
RS þ qRS ðxÞr, RS i 9J9 ½lnðL Þlim lnðrÞ
ð12Þ
r-0
In the previous equation, L is the size of the boundary element segment, J is the Jacobian matrix and r, i is a definite trigonometric quotient along the singular element. The limits on right hand side and left hand side when r tends to zero comprise part of an expression which is annulled in composition with other singular integral, presented next. In the direct hyper-singular formulations it is common for the obedience to complete the form of Eq. (12), but with the recursive procedure some simplifications are done to make this more accessible. Using linear shape functions one can locate the recursive source points at the center of the boundary element, making the previous equation simpler, because the size of segments are equal and no discontinuity occurs at the source point. Beyond this, a single function describes the behavior of q(x) and u(x) along the element, letting Eq. (12) to be expressed directly by "Z Z 1 1 qðXÞqi dG ¼ Q p jp ðXÞ r, LS 9J9 dr 2p LS r i sing Z 1 þ Q p jp ðXÞ r, RS 9J9dr ð13Þ r i RS Using constant boundary elements, a single function describes both the potential and flux variable along the element, but the position of source point in connection with adjacent elements requires the evaluation of the two last terms present on the right hand side of Eq. (12) unless the flux is constant and the boundary elements have the same size. The position of a source point between the connection point and the center of element must be avoided, since the next hyper-singular term is sensitive for this strategy, hindering good accuracy. Now it is shown the mathematical approach for the integral that contains the spatial derivative of the fundamental flux, which is hyper-singular. According to Mansur et al. [5], the same strategy adopted for the former integral must be implemented to guarantee the existence of hyper-singular integration in the CPV sense. Thus, considering straight elements, for which the first term in bracket of Eq. (8) is null, the equation becomes: Z
(Z 1 1 ½hLS ðXÞtLS ðxÞ LS LS ni 9J9 dr uðXÞuðxÞ 2 ni 9J9 dG ¼ 2p r r sing LS Z RS ½h ðXÞtRS ðxÞ RS RS ni 9J9 dr þ r RS LS
LS þhLS ðxÞnLS i 9J9 ½lnðL Þlim lnðrÞ r-0
RS
)
RS þhRS ðxÞnRS i 9J9 ½lnðL Þlim lnðrÞ r-0
ð14Þ
79
In the last equation, t(x) is the tangential potential derivative, found by taking the following limit on each side S:
tS ðxÞ ¼ lim r-0
½uðXÞuðxÞ r
ð15Þ
For brevity, an auxiliary function h(X) was also considered, and is defined by hðXÞ ¼
uðXÞuðxÞ r
ð16Þ
Eq. (14) explicit the general sense of CPV in hyper-singular integral equation, but it is also interesting to allow numerical evaluation of hyper-singular integrals, since Kutt integration points of first order [8] can be implemented, instead of complex second order ones. The evaluation of the former hyper-singular integral can be simplified in the recursive procedure when linear elements are employed, due to the reasons already presented for the analysis of the singular integral. The same occurs with constant elements, considering also the precautions previously exposed. For these simplified conditions, the hyper-singular integral can be solved analytically with advantage, since the concept of finite part integral can be used directly to evaluate the singularity of 1/r2, that is Z Z 1 p p 1 uðXÞuðxÞ 2 ni dG ¼ U j ðXÞuðxÞ 2 ni 9J9 dr ð17Þ r r sing sing It is important to point out that excluding the singular point, where a basic procedure of finite part integration is employed, the following result is valid along all boundaries: Z 1 n dG ¼ 0 ð18Þ uðxÞ 2 i Gr
6. Numerical simulations Three examples are given to demonstrate the performance of the recursive boundary element procedure, from now on referred as RBE, in comparison with traditional or direct boundary element solution (DBE). The use of coarse boundary element meshes in all examples is deliberate, because the main advantage of RBE would be to offer better accuracy using previous numerical results achieved with a restricted number of boundary element nodes. Due to continuity problems, RBE procedure is interesting to apply with linear elements, but in the first and second examples RBE is used with constant elements, exclusively with the purpose to demonstrating its performance. 6.1. First example—one-dimensional heat transfer As an introduction to the recursive procedure, a very simple problem is presented. A square domain is insulated along the superior and inferior horizontal edges, with a constant temperature and diffusive flux on the other edges, as shown in Fig. 1. Sixteen constant elements with the same size were used to discretize the boundary. Since thermal conductivity is unitary, the potential spatial derivatives are referred to as flux from this point on. Table 1 shows the numerical temperature results for the functional boundary nodes where flux values are prescribed, using DBE. The corresponding analytical results and the percentage error in absolute values are also included. Table 2 shows the numerical and analytical results for the normal flux on the boundary nodes. Due to the symmetry of results, only two nodal points are presented.
80
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84 Table 4 RBE normal flux in geometrical nodes. X
Y
NUM.
ANAL.
ER%
0.0 0.0
0.75 0.50
0.9823 0.9833
1.0 1.0
1.77 1.67
Fig. 1. Square domain with one-dimensional diffusive heat transfer.
Table 1 DBE temperatures at the functional boundary nodes.
Fig. 2. Semi-circle domain submitted to discontinuous temperature field.
X
Y
NUM.
ANAL.
ERROR%
0.125 0.375 0.625 0.875 0.100 0.100
0.000 0.000 0.000 0.000 0.125 0.375
0.1165 0.3661 0.6143 0.8616 0.9761 0.9864
0.125 0.375 0.625 0.875 1.0 1.0
6.80 2.37 1.71 1.53 2.39 1.36
Table 2 DBE normal flux at the functional boundary nodes. X
Y
NUM.
ANAL.
ER%
0.0 0.0
0.875 0.625
1.0350 0.9658
1.00 1.00
3.50 3.42
Table 5 DBE temperatures in the nodal points located on the circular boundary.
y (o )
NUM.
ANAL.
ERROR%
10 30 50 70 85 95 110 130 150 170
0.0517 0.1644 0.2764 0.3880 0.4719 0.5281 0.6120 0.7236 0.8355 0.9483
0.0555 0.1667 0.2778 0.3889 0.4722 0.5277 0.6111 0.7222 0.8333 0.9444
6.85 1.38 0.50 0.23 0.06 0.08 0.15 0.19 0.26 0.41
6.2. Second example—heat conduction in a semi-circle with singularity
Table 3 RBE numerical temperatures in geometrical nodes. X
Y
NUM.
ANAL.
ER%
0.25 0.50 0.75 1.0 1.0
0.0 0.0 0.0 0.25 0.5
0.2415 0.4903 0.7381 0.9840 0.9870
0.25 0.50 0.75 1.0 1.0
3.40 1.94 1.59 1.60 1.30
In Table 3 the temperature results obtained with the RBE are given. For convenience, the position of these source points was chosen to coincide with geometrical nodes of the boundary elements. Similarly, Table 4 shows recursive results for normal flux. RBE results are presented by source points located along the lines where the same variables were not previously prescribed. RBE results are more accurate than traditional direct process, especially for the flux results. The average percentage error for the flux results at the source points drops from 3.46% using DBE to 1.72% using RBE. For the temperature results, the average error decreases from 2.69% using DBE to 1.96% using RBE. Note that the use of new source points coinciding with the former nodal points for RBE modifies the original prescribed (or already calculated value) of a negligible value, which is due to the redistribution effect of weighted residuals.
The performance of RBE can also be evaluated for problems with high potential gradients, due to a singularity. A semi-circular domain was taken, insulated on the curved boundary and subject to a prescribed temperature along the complete horizontal boundary. Temperature values are discontinuous at r ¼0, as shown in Fig. 2. The BEM mesh used is comprised of 22 constant elements, with different sizes. For convenience, geometric nodes of boundary elements were taken as new points to RBE, excluding the horizontal corner points and the singularity point. Table 5 presents analytical and numerical DBE temperature results for the nodes at which the normal flux is prescribed, all located on the circular boundary. Despite the small number of boundary elements used, the numerical results show a very good accuracy, with an average percentage error being around 1.01%. Table 6 shows normal flux results for points on the horizontal boundary. In comparison with the temperature results, the expressive reduction of numerical accuracy of DBE should be noted for flux determination. Table 7 contains temperature results on the circular boundary obtained with RBE. It must be noted that the points chosen for RBE procedure are the meeting of two straight boundary elements, thus it necessary to determine the correct angle between them to evaluate c(x).
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
81
Table 6 DBE normal Flux in the nodal points located on the horizontal line. (X,0)
NUM.
ANAL.
ERROR%
3.50 2.50 1.75 1.25 0.75 0.25
0.0963 0.1251 0.1819 0.2502 0.3592 0.1557
0.0909 0.1273 0.1819 0.2546 0.4244 0.1273
5.89 1.73 0.00 1.73 15.36 22.31
Table 7 RBE temperatures in source points located on the circular boundary. X
Y
NUM.
ANAL.
ER%
3.76 3.06 2.00 .695 0.00 .695 2.00 3.06 3.76
1.37 2.57 3.46 3.94 4.00 3.94 3.46 2.57 1.37
0.1081 0.2204 0.3323 0.4437 0.5002 0.5545 0.6676 0.7794 0.8922
0.1111 0.2222 0.3333 0.4444 0.5000 0.5555 0.6667 0.7778 0.8889
2.70 0.81 0.30 0.16 0.00 0.18 0.14 0.20 0.37
Table 8 RBE normal fluxes in points of horizontal line. (X,0)
NUM.
ANAL.
ERROR%
3.0 2.0 1.5 1.0
0.1056 0.1587 0.2119 0.3181
.1061 .1591 .2122 .3183
0.47 0.25 0.14 0.06
Table 9 Normal flux comparison between DRB with linear boundary elements and RBE using constant elements. (X,0)
RBE CONST.
DBE LINEAR
ANAL.
0.80 0.45 0.17 0.09
0.3977 0.7072 1.872 3.537
0.3689 0.6732 1.472 5.975
0.3979 0.7073 1.819 3.537
The average percentage error obtained is around 0.54%. As was seen in the first example, RBE improved significantly the accuracy of direct numerical results for temperature, but not so significantly for flux results, which follows. In Table 8 the numerical results for the normal flux in points located at the horizontal line calculated recursively by RBE procedure are presented. In this case, the performance of RBE procedure was very effective, more accurate than the DBE results. To check the performance for flux determination, new points were introduced closer to the singularity and the results compared with results obtained with 25 isoparametric linear boundary elements. Results for both simulations are presented in Table 9. Note that RBE numerical results were more accurate for normal flux than the linear BEM, particularly at points located near the singularity, validating the good performance of the proposed procedure in this case.
Fig. 3. Square domain subject to non-uniform temperatures and fluxes.
Table 10 DBE temperatures where normal fluxes are prescribed. X
Y
NUM.
ANAL.
ER%
1.0 1.0 1.0 0.0 0.0 0.0
0.25 0.50 0.75 0.75 0.50 0.25
0.3820 0.7403 1.0532 0.6833 0.4805 0.2479
0.3818 0.7398 1.0518 0.6816 0.4794 0.2474
0.05 0.07 0.14 0.25 0.23 0.20
6.3. Third example—square control volume with variable prescribed conditions Fig. 3 presents the geometrical and physical characteristics of the third example, in which the main feature is that boundary conditions are not constant along the boundaries. Different from the previous problems, in this case the use of RBE with constant boundary elements in which the source points coincide with geometrical nodes is not suitable, since the values of potential and normal derivatives are discontinuous and the hyper-singular integrals are not convergent in CPV sense. An alternative procedure in which the source point is located at an intermediary position between the geometric and functional nodes does not produce good results for fluxes. This is probably due to the high sensitivity of the integral term given by Eq. (11). However, linear boundary elements are the most suitable boundary elements for implementing the recursive procedure. The position of recursive source points at the element center is very suitable, particularly due to the value of s(x), which is always equal to 0.5, and non-occurring discontinuity problems. Beyond this, at the center of the boundary elements the obedience to principal Cauchy conditions are more easily satisfied, because the size of segments on either side of each source point are equal. Table 10 shows potential results for nodes where normal derivatives are prescribed. Sixteen homogeneous linear elements are used, with double nodes at each corner, adding a total of twenty functional nodes. The average percentage error shown in Table 10 is 0.16%. For a more suitable comparison with recursive results, values for the double nodes were not introduced in the composition of the previous table, since they receive influence from the prescribed potential condition.
82
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
Table 11 presents recursive results for potential at source points located between the geometric nodes. The error for RBE averages 0.14%. RBE demonstrates better performance, especially at points located along the right vertical line, following the same tendency noticed in previous problems. Table 12 presents the results for normal flux. As previously stated, values of double nodes are not included. Results for recursive normal derivatives are presented in Table 13. Despite the good performance of the direct procedure (0.65%), the accuracy of the recursive procedure is still superior (0.25%). It must be noted that RBE requires a special integration scheme of Eq. (11) for accurate numerical computation of flux results, because this integral term has a low rate of convergence. Considering linear elements, no problems exist when potentials present constant or linear variation, as occurred with the previous
Table 11 RBE temperatures in central points along the elements. X
Y
NUM.
ANAL.
ER%
1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0
0.125 0.375 0.625 0.875 0.875 0.625 0.375 0.125
0.1924 0.5654 0.9035 1.1860 0.7696 0.5865 0.3671 0.1249
0.1924 0.5652 0.9028 1.1844 0.7675 0.5851 0.3663 0.1247
0.02 0.04 0.07 0.15 0.26 0.24 0.21 0.16
Table 12 DBE absolute values of normal fluxes where potential is prescribed. X
Y
NUM.
ANAL.
ERR%
0.25 0.50 0.75 0.75 0.50 0.25
0.0 0.0 0.0 1.0 1.0 1.0
1.0261 1.1220 1.2863 0.7072 0.6076 0.5621
1.0314 1.1276 1.2946 0.6995 0.6092 0.5573
0.52 0.50 0.64 1.10 0.27 0.86
Table 13 RBE absolute values of normal fluxes. X
Y
NUM.
ANAL.
ER%
0.125 0.375 0.625 0.875 0.875 0.625 0.375 0.125
0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
1.0075 1.0711 1.2011 1.4054 0.7679 0.6464 0.5743 0.5441
1.0078 1.0711 1.2017 1.4078 0.7606 0.6493 0.5787 0.5445
0.20 0.01 0.05 0.17 0.29 0.45 0.76 0.07
two examples, since the result of the aforementioned integral depends on the average nodal boundary values. However, when higher potential gradients are imposed, a more exact evaluation of this integral term is necessary. This is easily performed with RBE resources, through the computation of potential value u(x) at the element center, computed recursively, weighed by boundary nodal values at each extremity, given that both values are available. This proposed strategy improves the accuracy of the referred integral, considering additional information about the potential behavior along the element.
7. Recursive procedure and adaptive techniques Several procedures have been used to estimate BEM error. Considering the collocation BEM, the boundary integral equation procedure is closely related to the proposed recursive strategy. However, despite the similarities, there is a conceptual difference. Using collocation BEM, the residuals disappear at the initial collocation points, requiring the use of another technique for the computation of the necessary residual of the integral equation. In general, considering conforming linear boundary elements, some authors propose to estimate residuals placed in the additional collocation points between the nodal points, then solving the integral equation again. The residuals are achieved through the resulting difference between the solution in each new collocation point and values determined by the shape function along the boundary [9]. Usually, it is assumed that the error is maximum at the center of each element and zero at both ends [10]. If constant elements are used for the analysis, the additional collocation points are taken at the other place rather than the middle of each element. Whatever it is, the recursive procedure is used just to achieve missing information related to residuals, neglecting the sense of error minimization, which is implicit when the integral equation is used recursively. Commonly, residual and error behaviors are different. In the collocation weighted residuals method, the residuals are null at the collocation points, but the numerical solution is different from the exact solution. Using the collocation BEM it is not possible to affirm that the error is minimum at the collocation points or maximum between them. More accurate results in new collocation points do not indicate that the previous error is larger, but that it can be reduced with a recursive procedure, according to the achieved results. As a matter of fact, at element and element neighborhood level residual and error behaviors are not the same. The recursive procedure results indicate that when the integral residual sentence is applied again, a reduction of the errors committed at the new collocation points can possibly be achieved. For instance, taking a particular boundary element on the vertical line in example 3, results can be found by interpolation at some points located between two nodal points, shown in Table 14. The results achieved by recursive procedure, for the same points, are also presented in the same table.
Table 14 Comparison among analytical, interpolated and recursive values calculated at particular boundary element. X
Y
ANALYTICAL VALUES
INTERPOLATED VALUES
INTERPOLATED ERROR%
RECURSIVE VALUES
RECURS. ERROR%
1.0 1.0 1.0 1.0 1.0 1.0 1.0
0.27 0.30 0.35 .375 0.40 0.45 0.48
0.41159 0.45601 0.52912 0.56519 0.60090 0.67119 0.71256
0.41066 0.45366 0.52532 0.56100 0.59698 0.66864 0.71164
0.224 0.516 0.718 0.746 0.653 0.379 0.130
0.41179 0.45620 0.52931 0.56541 0.60117 0.67159 0.71306
0.048 0.041 0.037 0.039 0.044 0.060 0.070
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
83
Table 15 Absolute values for the residuals calculated at particular boundary element. X
Y
INTERPOLATED VALUES
RECURSIVE VALUES
ABSOLUTE RESIDUALS
1.0 1.0 1.0 1.0 1.0 1.0 1.0
0.27 0.30 0.35 0.375 0.40 0.45 0.48
0.41066 0.45366 0.52532 0.56100 0.59698 0.66864 0.71164
0.41179 0.45620 0.52931 0.56541 0.60117 0.67159 0.71306
0.00113 0.00254 0.00399 0.00441 0.00419 0.00295 0.00142
Fig. 4. Percentage error curve for RBE temperatures at points located along the right vertical line in third example.
Fig. 5. Percentage error curve for RBE temperatures at points located along the left vertical line in third example.
These results show better accuracy of the recursive procedure when compared to interpolated values. According to the standard methodology used in adaptive techniques, the residuals are achieved by subtraction of interpolated and recursive values, as shown in Table 15.
Fig. 6. Percentage error curve for RBE temperatures at points located along the horizontal line in first example.
The residuals distribution indicates higher values near the element center, according to the behavior of interpolated error along the boundary (see Table 14); however, it is very different from the recursive error distribution. At this moment, it is important to point out that searching strictly to establish a good residuals measure, adaptive technique research did not detect the better performance of the recursive integral equation, producing better results at new collocation points. For convenience, continuing with example 3, a continuous curve for the recursive temperature percentage errors is plotted along the two vertical lines in which the normal flux is prescribed. Initially is presented in Fig. 4, the results concerning for the right vertical line, where the flux is variable. Placed at the center of each boundary element, percentage error temperature results for basic recursive points are highlighted, such as the percentage error temperature results for the nodal points, according the presented legend. Considering 16 linear boundary elements, there is a clear tendency for the reduction of the percentage errors between the nodal points, indicating the proposed error minimization. This conclusion is not valid at corner proximities, where the recursive procedure does not seem to be effective. The mesh refinement smoothes the percentage error curve, such as is shown for the mesh using 32 homogeneous linear boundary elements. The behavior of temperature percentage error curve is different for the left vertical line, as shown in Fig. 5. There are no significant points of minimum between two direct collocation points. The recursive procedure does not produce much better
84
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
Fig. 7. Percentage error curve for RBE temperatures at points located along the horizontal line in first example using special error measure.
results than the direct ones, which is probably due to the uniformity of the prescribed flux condition along the boundary. Like an adaptive procedure, the efficiency of the recursive procedure depends on how distant the numerical solution is in respect to the exact solution. In the first example, the exact solution for the temperature is a straight line. Thus, despite using constant elements, the boundary element model is effective. Then, observe that the temperature error obtained at the recursive points located at the geometric nodes is only slightly below the average of errors produced between two consecutive nodal points. Fig. 6 shows the behavior of recursive percentage error curve obtained along the horizontal line. As aforementioned, the recursive procedure does not mitigate the corner problems created especially by the use of constant boundary elements; however, the high errors observed on the temperature recursive curve shown in Fig. 6 on the points near the corner are produced mainly due to the low values of analytical solution placed on the denominator of the standard expression of the percentage error. Another way to measure percentage error is by using the greatest value of analytical potential in the denominator. This produces the curve shown in Fig. 7.
8. Conclusions It is a well-known concept that the recursive use of boundary integral equation to determine values at internal domain points produces more accurate numerical results than the results on boundary nodes. However, the use of this strategy to recalculate the boundary values is not given even in the specialized BEM literature. Since the mathematical structure of BEM integral equation consists of an internal product, where residuals are forced to purge out the solution subspace, the recursive use of this equation with new source points can result in error reduction—this is the
aim of the present paper. This aspect is not considered in adaptive research. Despite the similarities, the strategy used in adaptive techniques is different: other collocation points for computation of the residuals are taken only to identify potential new points where the mesh needs to be refined, the element order to be increased or then to replace the nodal points, according to the adaptive technique chosen. According to the preliminary tests implemented, the recursive procedure acts effectively without any operational restrictions to improve the direct boundary integral solution, depending on how distant is the obtained numerical solution with respect to the exact solution. Thus, it is expected that RBE can be successfully used in more difficult scalar cases and especially in elasticity problems. Of course, the RBE is not intended to replace or overcome the adaptive techniques, which comprise a more elaborate numerical procedure. RBE procedure is especially convenient with linear isoparametric boundary elements, but it can be implemented with constant or higher order elements. To solve the basic variable recursively, the procedure is very simple; however, considering spatial derivatives, hyper-singular analysis becomes necessary to guarantee the convergence of complete integral expression. Finally, concerning the computational cost the RBE procedure is not expensive. Just for comparison, when using the adaptive schemes, normally a new matrix equation system will be solved. The total processing time is thus composed of the time to solve the original system, the time to evaluate the estimation errors, and the time to evaluate the definitive solution, if the procedure has already converged. In RBE the computational cost is equivalent to the expenditure necessary to calculate the solution at internal domain points, the number of new collocation points being smaller than the number of nodal points, since some values are known and prescribed. Beyond this, input data for recursive procedure is very easy to implement: the original boundary nodes can be used to generate automatically the new collocation points.
References [1] Brebbia CA, Telles JC, Wrobel LC. Boundary element techniques. Berlin: Springer Verlag; 1984. [2] Brebbia CA. The boundary element method for engineers.. London: Pentech Press; 1978. [3] Wrobel LC, Aliabadi MH. The boundary element method. Chichester: Wiley; 2002. [4] Brebbia CA, Domı´nguez J. The boundary element method—an introductory course. WIT Press; 1998. [5] Mansur WJ, Fleury Jr P, Azevedo JPS A. vector approach to the hyper-singular BEM formulation for Laplace’s equation in 2D. International Journal of BEM Communications 1997;8:239–50. [6] Telles JCF, Prado AA. Hyper-singular formulation for 2-D potential problems. In: Aliabadi MH, Brebbia CA, editors. Advanced formulations in boundary element methods. London: Elsevier; 1993. [Chapter 9]. [7] Pessolani RV. HP adaptive hierarchical formulation for elasticity with BEM. PhD thesis in civil engineering, Rio de Janeiro: COPEE-UFRJ; 1995 [in Portuguese]. [8] Kutt HR. On the numerical evaluation of finite part integrals involving an algebraic singularity. National Research Institute for Mathematical Sciences: Special report WISK 179. Pretoria, 1975. [9] Pessolani RV. An hp-adaptive hierarchical formulation for the boundary element method applied to elasticity in two dimensions. Journal of the Brazilian Society of Mechanical Sciences 2002;24(1):23–45. [10] Kita E, Kamiya N. Error estimation and adaptive mesh refinement in boundary element method, an overview. Engineering Analysis with Boundary Elements 2001;25:479–95.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A recursive application of the integral equation in the boundary element method C.F. Loeffler n ´ria, ES – CEP 29075-910, Brazil Mechanical Engineering Department, Federal University of Espı´rito Santo, Av. Fernando Ferrari, 514 Vito
a r t i c l e in fo
abstract
Article history: Received 5 August 2009 Accepted 26 May 2010 Available online 19 June 2010
This paper presents a recursive application of the governing integral equation aimed at improving the accuracy of numerical results of the boundary element method (BEM). Usually, only the results at internal domain points when using BEM are found using this approach, since the nodal boundary values have already been calculated. Here, it is shown that the same idea can be used to obtain better accuracy for the boundary results as well. Instead of locating the new source points inside the domain, they are positioned on the boundary, with different coordinates to the nodal points. The procedure is certainly general, but will be presented using as an example the two dimensional Laplace equation, for the sake of simplicity to point out the main concepts and numerical aspects of the method proposed, especially due to the determination of directional derivatives of the primal variable, which is part in hypersingular BEM theory. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Boundary element method Method of weighted residuals Scalar potential problems
1. Introduction Due to the mathematical modeling advances, improved computational capability and data storage, some traditional numerical techniques have been subjected to adaptations in order to achieve better performance. After procedures such as the relocation of discretization points, adaptive mesh refinement and adaptive shape functions, in some numerical methods an iterative solution scheme has been introduced in order to improve the numerical accuracy. Following this tendency, the recursive use of the boundary integral equation is here presented as an auxiliary resource to improve boundary element method (BEM) performance. It is not a completely iterative technique, but a simple procedure based on a common scheme to calculate internal values with BEM, where the boundary integral equation is used twice.
Laplace Equation can be written [1] as Z Z uðxÞ ¼ u ðx; XÞqðXÞ dG q ðx; XÞuðXÞ dG G
In Eq. (1), q(X) is the potential normal derivative. The essential (u ¼ u) and natural (u, i ni ¼ q) boundary conditions are defined on the boundary G ¼ Gu + Gq. The external unit normal vector at point X¼X(xi) is ni. For the usual direct BEM formulation, the auxiliary functions concern the solution of the related problem governed by a Poisson Equation [2], for which a unit singular source exists at X ¼ x for an infinite domain O (X). When the source point is positioned on the boundary, the integral equation for source points located on the boundary is given by Z Z cðxÞuðxÞ ¼ u ðx; XÞqðXÞ dGCPV q ðx; XÞuðXÞ dG ð2Þ G
2. Laplace boundary integral equation Considering an auxiliary functions un(x;X) and its normal derivative qn(x;X), where x is an internal domain point, and considering the basic mathematical tools of the theory of integral equations, the boundary integral equation associated to the n
Tel.: + 55 27 33352669. E-mail address: carlosloeffl[email protected]
0955-7997/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.05.012
ð1Þ
G
G
An ingenious mathematical procedure defines the c(x) value, which is equal to 0.5 for smooth boundaries [3]. When the source point is positioned on the boundary, un(x;X) and qn(x;X) functions present singularities. In spite of being discontinuous, un(x;X) is integrable along the boundary, but the integral of qn(x;X) exists in the Cauchy principal value (CPV) sense. According to the well-known BEM procedure [4], the boundary integral is discretized, generating a system of equations given in matrix notation as HYGQ ¼ 0
ð3Þ
78
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
3. Recursive boundary element procedure for potential variable The solution of BEM system of Eq. (3) allows the potential and its normal derivative for all N boundary nodes to be found. With regard to internal unknowns, it is possible to determinate easily their values putting the source point inside the domain and using the integral equation (1) again, but including the already calculated boundary variables. In general, the accuracy of numerical values of internal variables is usually better than the boundary ones. The mathematical reason for this behavior can be understood considering weighted residual fundamentals [2], since the approximated form of the boundary integral equation (1) can be expressed in a ^ that is ˆ and qE q, weighted residual form, taken as uEu Z Z Z ðx; XÞ dG ¼ 0 ^ ˆ ii u ðx; XÞ dO ¼ ˆ ðx; XÞ dG þ u, ðuuÞq ðqqÞu Gu
O
Gq
ð4Þ In Eq. (4), the left hand side integral represents the residual minimization in the domain and the right hand side integrals perform residual minimization of the essential and natural conditions on the boundary, respectively. Therefore, the recursive use of the integral equation to determinate numerical values for internal variables can be interpreted as a new process of weighting residuals, which expects a minimization of the errors previously performed, reaching best performance in the new collocation points, commonly located inside the domain. The same idea can be applied to recalculate boundary values with different positions on the boundary element. Thus, the recursive form of Eq. (4), considering previously calculated boundary values and putting new source point xi again on the boundary is given simply by Z Z N N X X i i Qke fk u ðxi ; XÞ dG Uke fk q ðxi ; XÞ dG ð5Þ cðx Þuðx Þ ¼ Ge
e¼1
e¼1
Ge
In Eq. (5), fk (X) are the shape functions, Qke and Uke are the nodal values on the boundary. The position of new source points must not coincide with former boundary nodes, because these points are just object of residual minimization; the recursive use of the integral equation is more effective using points situated between two adjacent boundary nodes.
4. Recursive boundary element procedure for potential spatial derivative Not only can the potential values be recalculated by the recursive procedure, but spatial potential derivatives related to xj (x) directions (or normal and tangential directions) can also be obtained by the same procedure. The mathematical foundation for determination of spatial derivatives follows same steps of the hyper-singular BEM formulation [4]. Let one consider, initially, an internal source point x. Taking the Cartesian derivatives of the integral equation (4) and using the Leibniz rule, it follows that Z Z u, i ðxÞ ¼ qi ðxÞ ¼ u, i ðx; XÞqðXÞ dG q, i ðx; XÞuðXÞ dG ð6Þ G
G
Considering smooth boundaries, s(x) is equal to 0.5 and the former integral involves only the potential derivatives in direction xi due to w(x) being null in this condition; however, for nonsmooth boundaries, the potential derivatives related to two different directions appear coupled in the same equation, as indicated by the permutation tensor e3ji. Concerning the boundary integrals, none of these is convergent in general; however, these two limits always exist when considered together, in CPV sense. Provided that the boundary values were previously calculated and considering smooth boundaries, Eq. (9) can be rewritten and applied to new xi points on the boundary: ( N Z X j j j j ½Uke jk uðx Þpi ðx ; XÞ dG sðx Þqi ðx Þ ¼ CPV e¼1
N X e¼1
Ge
Z Ge
) j Qke jk qi ðx ; XÞ dG
ð10Þ
According to the weighted residuals approach, some differences appear in the comparison between the former Eq. (10) and the equivalent sentence for potentials, given by Eq. (4). In spite of both integral equations making use the previous boundary values calculated on the functional nodes, the weighting functions are different, comprising different bases to residual orthogonalisation. This is an important factor for the difference in the numerical accuracy of recursive potential and derivative potential computations, such as occurs in the performance comparison between traditional singular and hyper-singular boundary techniques. However, the most important factor is the presence of u(x) subtracting the nodal variable Uke . This term, presented by convenience in Eq. (11), is responsible for the desired solution of the hyper-singular integral in a CPV sense, but its mathematical structure cannot be deduced from the basic weighted residual principles. N Z X j j j f ðx ; XÞ ¼ ½Uke fk uðx Þpi ðx ; XÞ dG ð11Þ e¼1
Ge
In many cases, due to the slow numerical convergence rate of this term, direct hyper-singular formulation presents worse performance than classical BEM formulation. The same behavior occurs with the recursive procedure for computation of the spatial derivative of the potential; thus, an effective strategy (presented next) must be used to guarantee good accuracy for the recursive technique, when potentials are not constant or vary in a different order to the shape functions used in the discretization procedure.
G
where: u, i ðx; XÞ ¼ qi ¼ ð1=2prÞr, i
analyzed, such as is done for the singular boundary integral. However, both kernels present singularities, which request more elaborated mathematical treatment [5,6]. Thus, locating the source point on the boundary, the hyper-singular integral equation for source points is given by Z sðxÞqi ðxÞ þ wðxÞe3ji qj ðxÞ ¼ CPV pi ðx; XÞ½uðXÞuðxÞ dG G Z qi ðx; XÞqðXÞ dG ð9Þ
2
q, i ðx; XÞ ¼ pi ¼ ð1=2pr Þ½2r, i r, j nj rj,i nj
5. Numerical aspects ð7Þ ð8Þ
To deduce the hyper-singular boundary integral equation, an augmented infinitesimal domain around the source point must be generated, and the effects created by this additional boundary
The existence of integrals in CPV sense for general cases to determine potential derivatives is guaranteed only when the two integrals on right hand side of Eq. (10) are analyzed together. Although the theory of hyper-singular integrals treatment is not the purpose of this paper, the main aspects are considered next for completeness.
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
First, the approach is shown when the integral contains the spatial derivative of the potential fundamental solution. Concerning general cases, initially it is important to divide this integral in two parts: the left hand side (LS) and the right hand side (RS) of the source point, using Hadamard’s methodology to treat the singularity [7], adding and subtracting a constant flux q(x), such as shown next: (Z Z 1 qLS ðXÞqLS ðxÞ LS LS r, i 9J9 dr qðXÞqi dG ¼ 2 r p sing LS Z qRS ðXÞqRS ðxÞ RS RS r, i 9J9 dr þ r RS LS
LS þ qLS ðxÞr, LS i 9J9 ½lnðL Þlim lnðrÞ r-0
RS
)
RS þ qRS ðxÞr, RS i 9J9 ½lnðL Þlim lnðrÞ
ð12Þ
r-0
In the previous equation, L is the size of the boundary element segment, J is the Jacobian matrix and r, i is a definite trigonometric quotient along the singular element. The limits on right hand side and left hand side when r tends to zero comprise part of an expression which is annulled in composition with other singular integral, presented next. In the direct hyper-singular formulations it is common for the obedience to complete the form of Eq. (12), but with the recursive procedure some simplifications are done to make this more accessible. Using linear shape functions one can locate the recursive source points at the center of the boundary element, making the previous equation simpler, because the size of segments are equal and no discontinuity occurs at the source point. Beyond this, a single function describes the behavior of q(x) and u(x) along the element, letting Eq. (12) to be expressed directly by "Z Z 1 1 qðXÞqi dG ¼ Q p jp ðXÞ r, LS 9J9 dr 2p LS r i sing Z 1 þ Q p jp ðXÞ r, RS 9J9dr ð13Þ r i RS Using constant boundary elements, a single function describes both the potential and flux variable along the element, but the position of source point in connection with adjacent elements requires the evaluation of the two last terms present on the right hand side of Eq. (12) unless the flux is constant and the boundary elements have the same size. The position of a source point between the connection point and the center of element must be avoided, since the next hyper-singular term is sensitive for this strategy, hindering good accuracy. Now it is shown the mathematical approach for the integral that contains the spatial derivative of the fundamental flux, which is hyper-singular. According to Mansur et al. [5], the same strategy adopted for the former integral must be implemented to guarantee the existence of hyper-singular integration in the CPV sense. Thus, considering straight elements, for which the first term in bracket of Eq. (8) is null, the equation becomes: Z
(Z 1 1 ½hLS ðXÞtLS ðxÞ LS LS ni 9J9 dr uðXÞuðxÞ 2 ni 9J9 dG ¼ 2p r r sing LS Z RS ½h ðXÞtRS ðxÞ RS RS ni 9J9 dr þ r RS LS
LS þhLS ðxÞnLS i 9J9 ½lnðL Þlim lnðrÞ r-0
RS
)
RS þhRS ðxÞnRS i 9J9 ½lnðL Þlim lnðrÞ r-0
ð14Þ
79
In the last equation, t(x) is the tangential potential derivative, found by taking the following limit on each side S:
tS ðxÞ ¼ lim r-0
½uðXÞuðxÞ r
ð15Þ
For brevity, an auxiliary function h(X) was also considered, and is defined by hðXÞ ¼
uðXÞuðxÞ r
ð16Þ
Eq. (14) explicit the general sense of CPV in hyper-singular integral equation, but it is also interesting to allow numerical evaluation of hyper-singular integrals, since Kutt integration points of first order [8] can be implemented, instead of complex second order ones. The evaluation of the former hyper-singular integral can be simplified in the recursive procedure when linear elements are employed, due to the reasons already presented for the analysis of the singular integral. The same occurs with constant elements, considering also the precautions previously exposed. For these simplified conditions, the hyper-singular integral can be solved analytically with advantage, since the concept of finite part integral can be used directly to evaluate the singularity of 1/r2, that is Z Z 1 p p 1 uðXÞuðxÞ 2 ni dG ¼ U j ðXÞuðxÞ 2 ni 9J9 dr ð17Þ r r sing sing It is important to point out that excluding the singular point, where a basic procedure of finite part integration is employed, the following result is valid along all boundaries: Z 1 n dG ¼ 0 ð18Þ uðxÞ 2 i Gr
6. Numerical simulations Three examples are given to demonstrate the performance of the recursive boundary element procedure, from now on referred as RBE, in comparison with traditional or direct boundary element solution (DBE). The use of coarse boundary element meshes in all examples is deliberate, because the main advantage of RBE would be to offer better accuracy using previous numerical results achieved with a restricted number of boundary element nodes. Due to continuity problems, RBE procedure is interesting to apply with linear elements, but in the first and second examples RBE is used with constant elements, exclusively with the purpose to demonstrating its performance. 6.1. First example—one-dimensional heat transfer As an introduction to the recursive procedure, a very simple problem is presented. A square domain is insulated along the superior and inferior horizontal edges, with a constant temperature and diffusive flux on the other edges, as shown in Fig. 1. Sixteen constant elements with the same size were used to discretize the boundary. Since thermal conductivity is unitary, the potential spatial derivatives are referred to as flux from this point on. Table 1 shows the numerical temperature results for the functional boundary nodes where flux values are prescribed, using DBE. The corresponding analytical results and the percentage error in absolute values are also included. Table 2 shows the numerical and analytical results for the normal flux on the boundary nodes. Due to the symmetry of results, only two nodal points are presented.
80
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84 Table 4 RBE normal flux in geometrical nodes. X
Y
NUM.
ANAL.
ER%
0.0 0.0
0.75 0.50
0.9823 0.9833
1.0 1.0
1.77 1.67
Fig. 1. Square domain with one-dimensional diffusive heat transfer.
Table 1 DBE temperatures at the functional boundary nodes.
Fig. 2. Semi-circle domain submitted to discontinuous temperature field.
X
Y
NUM.
ANAL.
ERROR%
0.125 0.375 0.625 0.875 0.100 0.100
0.000 0.000 0.000 0.000 0.125 0.375
0.1165 0.3661 0.6143 0.8616 0.9761 0.9864
0.125 0.375 0.625 0.875 1.0 1.0
6.80 2.37 1.71 1.53 2.39 1.36
Table 2 DBE normal flux at the functional boundary nodes. X
Y
NUM.
ANAL.
ER%
0.0 0.0
0.875 0.625
1.0350 0.9658
1.00 1.00
3.50 3.42
Table 5 DBE temperatures in the nodal points located on the circular boundary.
y (o )
NUM.
ANAL.
ERROR%
10 30 50 70 85 95 110 130 150 170
0.0517 0.1644 0.2764 0.3880 0.4719 0.5281 0.6120 0.7236 0.8355 0.9483
0.0555 0.1667 0.2778 0.3889 0.4722 0.5277 0.6111 0.7222 0.8333 0.9444
6.85 1.38 0.50 0.23 0.06 0.08 0.15 0.19 0.26 0.41
6.2. Second example—heat conduction in a semi-circle with singularity
Table 3 RBE numerical temperatures in geometrical nodes. X
Y
NUM.
ANAL.
ER%
0.25 0.50 0.75 1.0 1.0
0.0 0.0 0.0 0.25 0.5
0.2415 0.4903 0.7381 0.9840 0.9870
0.25 0.50 0.75 1.0 1.0
3.40 1.94 1.59 1.60 1.30
In Table 3 the temperature results obtained with the RBE are given. For convenience, the position of these source points was chosen to coincide with geometrical nodes of the boundary elements. Similarly, Table 4 shows recursive results for normal flux. RBE results are presented by source points located along the lines where the same variables were not previously prescribed. RBE results are more accurate than traditional direct process, especially for the flux results. The average percentage error for the flux results at the source points drops from 3.46% using DBE to 1.72% using RBE. For the temperature results, the average error decreases from 2.69% using DBE to 1.96% using RBE. Note that the use of new source points coinciding with the former nodal points for RBE modifies the original prescribed (or already calculated value) of a negligible value, which is due to the redistribution effect of weighted residuals.
The performance of RBE can also be evaluated for problems with high potential gradients, due to a singularity. A semi-circular domain was taken, insulated on the curved boundary and subject to a prescribed temperature along the complete horizontal boundary. Temperature values are discontinuous at r ¼0, as shown in Fig. 2. The BEM mesh used is comprised of 22 constant elements, with different sizes. For convenience, geometric nodes of boundary elements were taken as new points to RBE, excluding the horizontal corner points and the singularity point. Table 5 presents analytical and numerical DBE temperature results for the nodes at which the normal flux is prescribed, all located on the circular boundary. Despite the small number of boundary elements used, the numerical results show a very good accuracy, with an average percentage error being around 1.01%. Table 6 shows normal flux results for points on the horizontal boundary. In comparison with the temperature results, the expressive reduction of numerical accuracy of DBE should be noted for flux determination. Table 7 contains temperature results on the circular boundary obtained with RBE. It must be noted that the points chosen for RBE procedure are the meeting of two straight boundary elements, thus it necessary to determine the correct angle between them to evaluate c(x).
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
81
Table 6 DBE normal Flux in the nodal points located on the horizontal line. (X,0)
NUM.
ANAL.
ERROR%
3.50 2.50 1.75 1.25 0.75 0.25
0.0963 0.1251 0.1819 0.2502 0.3592 0.1557
0.0909 0.1273 0.1819 0.2546 0.4244 0.1273
5.89 1.73 0.00 1.73 15.36 22.31
Table 7 RBE temperatures in source points located on the circular boundary. X
Y
NUM.
ANAL.
ER%
3.76 3.06 2.00 .695 0.00 .695 2.00 3.06 3.76
1.37 2.57 3.46 3.94 4.00 3.94 3.46 2.57 1.37
0.1081 0.2204 0.3323 0.4437 0.5002 0.5545 0.6676 0.7794 0.8922
0.1111 0.2222 0.3333 0.4444 0.5000 0.5555 0.6667 0.7778 0.8889
2.70 0.81 0.30 0.16 0.00 0.18 0.14 0.20 0.37
Table 8 RBE normal fluxes in points of horizontal line. (X,0)
NUM.
ANAL.
ERROR%
3.0 2.0 1.5 1.0
0.1056 0.1587 0.2119 0.3181
.1061 .1591 .2122 .3183
0.47 0.25 0.14 0.06
Table 9 Normal flux comparison between DRB with linear boundary elements and RBE using constant elements. (X,0)
RBE CONST.
DBE LINEAR
ANAL.
0.80 0.45 0.17 0.09
0.3977 0.7072 1.872 3.537
0.3689 0.6732 1.472 5.975
0.3979 0.7073 1.819 3.537
The average percentage error obtained is around 0.54%. As was seen in the first example, RBE improved significantly the accuracy of direct numerical results for temperature, but not so significantly for flux results, which follows. In Table 8 the numerical results for the normal flux in points located at the horizontal line calculated recursively by RBE procedure are presented. In this case, the performance of RBE procedure was very effective, more accurate than the DBE results. To check the performance for flux determination, new points were introduced closer to the singularity and the results compared with results obtained with 25 isoparametric linear boundary elements. Results for both simulations are presented in Table 9. Note that RBE numerical results were more accurate for normal flux than the linear BEM, particularly at points located near the singularity, validating the good performance of the proposed procedure in this case.
Fig. 3. Square domain subject to non-uniform temperatures and fluxes.
Table 10 DBE temperatures where normal fluxes are prescribed. X
Y
NUM.
ANAL.
ER%
1.0 1.0 1.0 0.0 0.0 0.0
0.25 0.50 0.75 0.75 0.50 0.25
0.3820 0.7403 1.0532 0.6833 0.4805 0.2479
0.3818 0.7398 1.0518 0.6816 0.4794 0.2474
0.05 0.07 0.14 0.25 0.23 0.20
6.3. Third example—square control volume with variable prescribed conditions Fig. 3 presents the geometrical and physical characteristics of the third example, in which the main feature is that boundary conditions are not constant along the boundaries. Different from the previous problems, in this case the use of RBE with constant boundary elements in which the source points coincide with geometrical nodes is not suitable, since the values of potential and normal derivatives are discontinuous and the hyper-singular integrals are not convergent in CPV sense. An alternative procedure in which the source point is located at an intermediary position between the geometric and functional nodes does not produce good results for fluxes. This is probably due to the high sensitivity of the integral term given by Eq. (11). However, linear boundary elements are the most suitable boundary elements for implementing the recursive procedure. The position of recursive source points at the element center is very suitable, particularly due to the value of s(x), which is always equal to 0.5, and non-occurring discontinuity problems. Beyond this, at the center of the boundary elements the obedience to principal Cauchy conditions are more easily satisfied, because the size of segments on either side of each source point are equal. Table 10 shows potential results for nodes where normal derivatives are prescribed. Sixteen homogeneous linear elements are used, with double nodes at each corner, adding a total of twenty functional nodes. The average percentage error shown in Table 10 is 0.16%. For a more suitable comparison with recursive results, values for the double nodes were not introduced in the composition of the previous table, since they receive influence from the prescribed potential condition.
82
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
Table 11 presents recursive results for potential at source points located between the geometric nodes. The error for RBE averages 0.14%. RBE demonstrates better performance, especially at points located along the right vertical line, following the same tendency noticed in previous problems. Table 12 presents the results for normal flux. As previously stated, values of double nodes are not included. Results for recursive normal derivatives are presented in Table 13. Despite the good performance of the direct procedure (0.65%), the accuracy of the recursive procedure is still superior (0.25%). It must be noted that RBE requires a special integration scheme of Eq. (11) for accurate numerical computation of flux results, because this integral term has a low rate of convergence. Considering linear elements, no problems exist when potentials present constant or linear variation, as occurred with the previous
Table 11 RBE temperatures in central points along the elements. X
Y
NUM.
ANAL.
ER%
1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0
0.125 0.375 0.625 0.875 0.875 0.625 0.375 0.125
0.1924 0.5654 0.9035 1.1860 0.7696 0.5865 0.3671 0.1249
0.1924 0.5652 0.9028 1.1844 0.7675 0.5851 0.3663 0.1247
0.02 0.04 0.07 0.15 0.26 0.24 0.21 0.16
Table 12 DBE absolute values of normal fluxes where potential is prescribed. X
Y
NUM.
ANAL.
ERR%
0.25 0.50 0.75 0.75 0.50 0.25
0.0 0.0 0.0 1.0 1.0 1.0
1.0261 1.1220 1.2863 0.7072 0.6076 0.5621
1.0314 1.1276 1.2946 0.6995 0.6092 0.5573
0.52 0.50 0.64 1.10 0.27 0.86
Table 13 RBE absolute values of normal fluxes. X
Y
NUM.
ANAL.
ER%
0.125 0.375 0.625 0.875 0.875 0.625 0.375 0.125
0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
1.0075 1.0711 1.2011 1.4054 0.7679 0.6464 0.5743 0.5441
1.0078 1.0711 1.2017 1.4078 0.7606 0.6493 0.5787 0.5445
0.20 0.01 0.05 0.17 0.29 0.45 0.76 0.07
two examples, since the result of the aforementioned integral depends on the average nodal boundary values. However, when higher potential gradients are imposed, a more exact evaluation of this integral term is necessary. This is easily performed with RBE resources, through the computation of potential value u(x) at the element center, computed recursively, weighed by boundary nodal values at each extremity, given that both values are available. This proposed strategy improves the accuracy of the referred integral, considering additional information about the potential behavior along the element.
7. Recursive procedure and adaptive techniques Several procedures have been used to estimate BEM error. Considering the collocation BEM, the boundary integral equation procedure is closely related to the proposed recursive strategy. However, despite the similarities, there is a conceptual difference. Using collocation BEM, the residuals disappear at the initial collocation points, requiring the use of another technique for the computation of the necessary residual of the integral equation. In general, considering conforming linear boundary elements, some authors propose to estimate residuals placed in the additional collocation points between the nodal points, then solving the integral equation again. The residuals are achieved through the resulting difference between the solution in each new collocation point and values determined by the shape function along the boundary [9]. Usually, it is assumed that the error is maximum at the center of each element and zero at both ends [10]. If constant elements are used for the analysis, the additional collocation points are taken at the other place rather than the middle of each element. Whatever it is, the recursive procedure is used just to achieve missing information related to residuals, neglecting the sense of error minimization, which is implicit when the integral equation is used recursively. Commonly, residual and error behaviors are different. In the collocation weighted residuals method, the residuals are null at the collocation points, but the numerical solution is different from the exact solution. Using the collocation BEM it is not possible to affirm that the error is minimum at the collocation points or maximum between them. More accurate results in new collocation points do not indicate that the previous error is larger, but that it can be reduced with a recursive procedure, according to the achieved results. As a matter of fact, at element and element neighborhood level residual and error behaviors are not the same. The recursive procedure results indicate that when the integral residual sentence is applied again, a reduction of the errors committed at the new collocation points can possibly be achieved. For instance, taking a particular boundary element on the vertical line in example 3, results can be found by interpolation at some points located between two nodal points, shown in Table 14. The results achieved by recursive procedure, for the same points, are also presented in the same table.
Table 14 Comparison among analytical, interpolated and recursive values calculated at particular boundary element. X
Y
ANALYTICAL VALUES
INTERPOLATED VALUES
INTERPOLATED ERROR%
RECURSIVE VALUES
RECURS. ERROR%
1.0 1.0 1.0 1.0 1.0 1.0 1.0
0.27 0.30 0.35 .375 0.40 0.45 0.48
0.41159 0.45601 0.52912 0.56519 0.60090 0.67119 0.71256
0.41066 0.45366 0.52532 0.56100 0.59698 0.66864 0.71164
0.224 0.516 0.718 0.746 0.653 0.379 0.130
0.41179 0.45620 0.52931 0.56541 0.60117 0.67159 0.71306
0.048 0.041 0.037 0.039 0.044 0.060 0.070
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
83
Table 15 Absolute values for the residuals calculated at particular boundary element. X
Y
INTERPOLATED VALUES
RECURSIVE VALUES
ABSOLUTE RESIDUALS
1.0 1.0 1.0 1.0 1.0 1.0 1.0
0.27 0.30 0.35 0.375 0.40 0.45 0.48
0.41066 0.45366 0.52532 0.56100 0.59698 0.66864 0.71164
0.41179 0.45620 0.52931 0.56541 0.60117 0.67159 0.71306
0.00113 0.00254 0.00399 0.00441 0.00419 0.00295 0.00142
Fig. 4. Percentage error curve for RBE temperatures at points located along the right vertical line in third example.
Fig. 5. Percentage error curve for RBE temperatures at points located along the left vertical line in third example.
These results show better accuracy of the recursive procedure when compared to interpolated values. According to the standard methodology used in adaptive techniques, the residuals are achieved by subtraction of interpolated and recursive values, as shown in Table 15.
Fig. 6. Percentage error curve for RBE temperatures at points located along the horizontal line in first example.
The residuals distribution indicates higher values near the element center, according to the behavior of interpolated error along the boundary (see Table 14); however, it is very different from the recursive error distribution. At this moment, it is important to point out that searching strictly to establish a good residuals measure, adaptive technique research did not detect the better performance of the recursive integral equation, producing better results at new collocation points. For convenience, continuing with example 3, a continuous curve for the recursive temperature percentage errors is plotted along the two vertical lines in which the normal flux is prescribed. Initially is presented in Fig. 4, the results concerning for the right vertical line, where the flux is variable. Placed at the center of each boundary element, percentage error temperature results for basic recursive points are highlighted, such as the percentage error temperature results for the nodal points, according the presented legend. Considering 16 linear boundary elements, there is a clear tendency for the reduction of the percentage errors between the nodal points, indicating the proposed error minimization. This conclusion is not valid at corner proximities, where the recursive procedure does not seem to be effective. The mesh refinement smoothes the percentage error curve, such as is shown for the mesh using 32 homogeneous linear boundary elements. The behavior of temperature percentage error curve is different for the left vertical line, as shown in Fig. 5. There are no significant points of minimum between two direct collocation points. The recursive procedure does not produce much better
84
C.F. Loeffler / Engineering Analysis with Boundary Elements 35 (2011) 77–84
Fig. 7. Percentage error curve for RBE temperatures at points located along the horizontal line in first example using special error measure.
results than the direct ones, which is probably due to the uniformity of the prescribed flux condition along the boundary. Like an adaptive procedure, the efficiency of the recursive procedure depends on how distant the numerical solution is in respect to the exact solution. In the first example, the exact solution for the temperature is a straight line. Thus, despite using constant elements, the boundary element model is effective. Then, observe that the temperature error obtained at the recursive points located at the geometric nodes is only slightly below the average of errors produced between two consecutive nodal points. Fig. 6 shows the behavior of recursive percentage error curve obtained along the horizontal line. As aforementioned, the recursive procedure does not mitigate the corner problems created especially by the use of constant boundary elements; however, the high errors observed on the temperature recursive curve shown in Fig. 6 on the points near the corner are produced mainly due to the low values of analytical solution placed on the denominator of the standard expression of the percentage error. Another way to measure percentage error is by using the greatest value of analytical potential in the denominator. This produces the curve shown in Fig. 7.
8. Conclusions It is a well-known concept that the recursive use of boundary integral equation to determine values at internal domain points produces more accurate numerical results than the results on boundary nodes. However, the use of this strategy to recalculate the boundary values is not given even in the specialized BEM literature. Since the mathematical structure of BEM integral equation consists of an internal product, where residuals are forced to purge out the solution subspace, the recursive use of this equation with new source points can result in error reduction—this is the
aim of the present paper. This aspect is not considered in adaptive research. Despite the similarities, the strategy used in adaptive techniques is different: other collocation points for computation of the residuals are taken only to identify potential new points where the mesh needs to be refined, the element order to be increased or then to replace the nodal points, according to the adaptive technique chosen. According to the preliminary tests implemented, the recursive procedure acts effectively without any operational restrictions to improve the direct boundary integral solution, depending on how distant is the obtained numerical solution with respect to the exact solution. Thus, it is expected that RBE can be successfully used in more difficult scalar cases and especially in elasticity problems. Of course, the RBE is not intended to replace or overcome the adaptive techniques, which comprise a more elaborate numerical procedure. RBE procedure is especially convenient with linear isoparametric boundary elements, but it can be implemented with constant or higher order elements. To solve the basic variable recursively, the procedure is very simple; however, considering spatial derivatives, hyper-singular analysis becomes necessary to guarantee the convergence of complete integral expression. Finally, concerning the computational cost the RBE procedure is not expensive. Just for comparison, when using the adaptive schemes, normally a new matrix equation system will be solved. The total processing time is thus composed of the time to solve the original system, the time to evaluate the estimation errors, and the time to evaluate the definitive solution, if the procedure has already converged. In RBE the computational cost is equivalent to the expenditure necessary to calculate the solution at internal domain points, the number of new collocation points being smaller than the number of nodal points, since some values are known and prescribed. Beyond this, input data for recursive procedure is very easy to implement: the original boundary nodes can be used to generate automatically the new collocation points.
References [1] Brebbia CA, Telles JC, Wrobel LC. Boundary element techniques. Berlin: Springer Verlag; 1984. [2] Brebbia CA. The boundary element method for engineers.. London: Pentech Press; 1978. [3] Wrobel LC, Aliabadi MH. The boundary element method. Chichester: Wiley; 2002. [4] Brebbia CA, Domı´nguez J. The boundary element method—an introductory course. WIT Press; 1998. [5] Mansur WJ, Fleury Jr P, Azevedo JPS A. vector approach to the hyper-singular BEM formulation for Laplace’s equation in 2D. International Journal of BEM Communications 1997;8:239–50. [6] Telles JCF, Prado AA. Hyper-singular formulation for 2-D potential problems. In: Aliabadi MH, Brebbia CA, editors. Advanced formulations in boundary element methods. London: Elsevier; 1993. [Chapter 9]. [7] Pessolani RV. HP adaptive hierarchical formulation for elasticity with BEM. PhD thesis in civil engineering, Rio de Janeiro: COPEE-UFRJ; 1995 [in Portuguese]. [8] Kutt HR. On the numerical evaluation of finite part integrals involving an algebraic singularity. National Research Institute for Mathematical Sciences: Special report WISK 179. Pretoria, 1975. [9] Pessolani RV. An hp-adaptive hierarchical formulation for the boundary element method applied to elasticity in two dimensions. Journal of the Brazilian Society of Mechanical Sciences 2002;24(1):23–45. [10] Kita E, Kamiya N. Error estimation and adaptive mesh refinement in boundary element method, an overview. Engineering Analysis with Boundary Elements 2001;25:479–95.