CC BY
34
Вычислительные технологии
Том 3, № 1, 1998
ANGULAR INTEGRAL METHOD FOR THE DIRICHLET PROBLEM
L. L. Frumin Novosibirsk State University, Russia e-mail: frumin@phys.nsu.ru
Для решения задачи Дирихле в рамках потенциальной теории предлагается модификация метода граничных элементов, основанная на замене интегрирования в уравнении потенциала двойного слоя по области границы на интегрирование по углу. Использование угловой переменной, представляющей собой естественную параметризацию границы, позволяет полностью исключить необходимость ее аппроксимации и существенно улучшает точность вычислений. Описаны результаты тестирования алгоритма на двумерных областях.
1. Introduction
Recently, the boundary elements method (BEM) has been used frequently along with the other method for solving the boundary-value problems [1, 2]. It is widely used for calculations in the theory of water waves, in problems on elastic and plastic deformation, in electrostatics and in other practical applications.
One of the problems inherent to this method is a necessity to approximate both the functions on the boundary and the boundary itself. In practice, most frequently used are the simple low-order schemes. Note, however, that a low order of boundary approximation does not justify a higher order of approximation of the function itself, as the latter does not guarantee an improvement in accuracy [3].
The present paper proposes an approach based on a substitution, in the framework of the potential method, of the integral over the domain boundary by an integral over the angle. This approach does not require any approximation of the boundary, opening new prospects for solving the problem numerically.
Let us first examine the problem setup and its solution in the framework of the BEM. As an example, we shall use the internal Dirichlet problem in a two-dimensional domain.
2. Problem setup
In the framework of the potential theory [4], finding the solution u(p),p £ Q, of the boundary-value Dirichlet problem in the domain Q £ R2, satisfying within the domain the Laplace equation and the boundary condition u = f on the boundary r, may be reduced to the problem
© L.L. Frumin, 1998.
of calculating the double layer potential m on the same boundary. For the latter, the following integral equation is valid:
+ = f(p), p,q € r, (1)
where rpq is the radius vector from the boundary point p to the point q, nq is the outward unit normal to the boundary at the point q, dlq is the length element of the boundary r. The direction of circulation is chosen so that the domain Q stays always at the left side.
Solution of the Dirichlet problem u(p) in the domain Q is determined by the integral
u(p) = 7T [M(q)(-nqrq)dlq, q € r, p € Q. (2)
2n r rpq
In the simplest case, the BEM relies upon an approximation of the domain boundary by a polyline formed by m segments boundary elements. The discrete grid functions Mj = M(pj) and fj = f (pj) are defined on these boundary elements with centres in the points pj, j = 1, ...,m. Thus, the problem for the unknown m is reduced to a system of m linear algebraic equations of the following form:
1 m
2 Mi + Y1 aij Mj = 1 = 1,-,m- (3)
j=i
The matrix aij is defined by the expression
aij = jij )hj /2nrij, i = j, (4)
where rij is a vector from pi to pj, rij = |rij|, hj is the size of the j-th boundary element, nj is the given external unit normal vector to the respective boundary element.
The Gauss theorem is used to calculate the singular diagonal elements of the aii matrix:
1 m ,
an = 2 - S aij■ (5)
j=l,i=j
The system of linear algebraic equations (3) is then solved to determine the double layer potential m
3. Angular integral method
Considering now the new algorithm, note that integration over the domain boundary in (1) may be substituted by integration over the angle 9 [4]:
dO = (flqTpq )dlq/r^. (6)
Let {pj € r}, j = 1,...,m be the set of boundary points. Let us define 9i for a fixed point pi as the angle between the abscissa of the Cartesian reference system with an origin at pi and the respective radius vector. For convenience, we shall rotate the reference system for the ordinate to coincide with the tangent to the boundary at the point pi, leaving alone for a while the eventual angular points. Thus, every point pj = pi will correspond to an angle 9ij given by
9ij = arctan | — Xl
yj- y
ANGULAR INTEGRAL METHOD FOR THE DIRICHLET PROBLEM
77
This angle defines in fact a natural boundary parametrisation, which enables to avoid its approximation by the boundary elements.
Note that, in case of a convex smooth boundary without angular points, the integration limits after rotating the reference system as explained above will be — n/2, +n/2.
Thus, the problem of finding the matrix elements aij is reduced to finding the numerical integration formula coefficients for calculation of integrals over the angle using the integrated function values in the points of an (uneven) grid 9ij, defined by the points pj, j = 1, ...,m, at a fixed pi, in the interval from —n/2 to +n/2. This problem may be solved using any suitable numerical integration formulae. In particular, rectangular and trapezoidal integration formulas were tested.
For the rectangular integration, the matrix elements aij are given by the relations:
aij = (9ij+i — 9ij )/2n,
9ij+i = arctan ( X+1—X ) , (7)
\Xj+1 — XiJ
note that 9ij+1 = +n/2 for j + 1 = i, and 9ij = —n/2 for j = i. Similarly, for the trapezoidal integration it is easy to obtain:
a
a
(Oij+1 — 0ij-1 + n)/4n for i = j, where (8)
9j+i = arctan ( X+i_Xl), j = arcta^ XOzL-XL \ , (9)
\Xj+1 xi J \Xj-1 Xi J
note that 9ij+1 = n/2 for j + 1 = i, and 9ij_1 = —n/2 for j — 1 = i.
The algorithm was tested for a known solution in a circular domain. Using a polar reference system with the angle a, whose origin coincided with the centre of the circle, the exact distribution of the double layer potential ^ was chosen in the form:
^(a) = a sin(a). (10)
Thus the right term of the equation (1) equals:
f (a) = 2(a sin(a) — 1). (11)
The Table presents a comparison of the relative calculation uncertainty e (%) for the standard BEM as described above, with that for the new algorithm as presently proposed, both applied to the test problem using various values of m.
Method Number of points
10 20 30 40 50
Form. (4), (5) (BEM) 3.339 0.787 0.346 0.194 0.124
Form. (8) (AIM) 0.658 0.164 0.073 0.041 0.026
The calculations demonstrated that the approach proposed improves significantly the accuracy of the numerical solution of the Dirichlet problem in a circular domain as compared to the traditional approach expressed by equations (4)-(5), provided the same number of points is used in the two cases. This is true even when the simplest rectangular integration is used. As is evident from the Table, the relative uncertainty for the rectangular integration is almost a factor of 5 lower as that for the calculations using formulas (4)-(5), being the same as the uncertainty for trapezoidal integration for the solution (10) in a circle.
4. Conclusion
We shall finish with a few notes regarding the non-convex domains and the three-dimensional generalisations of the algorithm.
The approximation order of the function on the boundary may be improved by, for example, using a polynomial of the order m. As regards non-convex domains, here we have to integrate a non-uniqueness function. In this case the integral should be separated into integrals over the respective uniqueness intervals. Note that here one has to integrate in an interval wider than [-n/2,n/2].
The algorithm may be also generalised for three dimensions. In this case the problem is reduced to calculations of quadrature of a two-dimensional integral over the solid angle.
Finally, the algorithm proposed offers new possibilities to deal with the angular points of the boundary, requiring a special consideration, which is out of the scope of the present paper.
The author expresses his gratitude to V. L. Frumin and L. B. Chubarov for the useful discussions and also to I. V. Khmelinski for his help in translation.
References
[1] Brebbia c.a., Telles j.c.f., wrobel l.c. Boundary Element Techniques. Springer-Verlag, Berlin and N. Y., 1984.
[2] Boundary-Integral Equation Method. Eds. t. a. cruse and f.j. rlzzo Comput. Applications in Appl. Mech., ASME, N. Y., 1975.
[3] hess j. l. Review of integral-equation techniques for solving potential-flow problems with emphasison the surface-source method. Comput. Methods in Appl. Mech. and Eng., 5, 1975, 145-196.
[4] gahoy f. d. Boundary-Value Problems (in Russian). GIFML, Moscow, 1963.
Received for publication June 17, 1997, in revised form November 30, 1997
