Numerical integration in the vicinity of a logarithmic singularity

Numerical integration in the vicinity of a logarithmic singularity JEAN-LOUIS MIGEOT National Fund for Scientific Research, Universitd Libre de Bruxelles, Continuum Mechanics Department - CP 194/5, 50 av. Roosevelt, 1050 Bruxelles, Belgium

This paper presents a study of the numerical integration of the (ln r) function in the vicinity of the singularity. This problem occurs for instance in the solution of field problems by the boundary element method. 1 Two methods are considered, the Gauss-Legendre method and the analytical one. As a first step, the computing time of both methods are compared and an upper limit of the number of economic Gauss points is found. We then compute the error for various observation points and plot the iso-e curves from which we derive a rule giving the number of Gauss points necessary to obtain a given precision as a function of the relative distance to the mid-point of the element. More surprising is the occurrence of lines of vanishing error. This last point is of great importance, in the boundary integral equation method for instance, as it allows important time economy in the integration.

INTRODUCTION In the BIEM (Boundary Integral Equation Method) solution o f both elastostatic and electrostatic problems, the following integrals must be computed: ~

f

f(x) In 1 ds

c

in which r is the distance between the point (X, Y) in the domain where the field is computed, which we call the observation point, and a current point (x, y ) of the boundary element C (Fig. 1) and f(X) a polynomial of degree

0 (constant element) or 1 (linear element) of the parameter X = 8/l (Fig. 1) varying on C between 0 and 1. We callf(X) the shape function of the element. Analytical expressions of these integrals can be found for constant and linear shape functions defined over geometrically linear elements. 2 Analytical integration is obviously more precise than the numerical one but it is quite costly. In the next paragraph, we compare the computation time of both methods.

COMPARISON OF THE COMPUTATION TIME For each method, i.e. the analytic method and the GaussLegendre method with 1, 2, 3, 4 and 5 points, we have computed the same integral 1000 and 4000 times and from this we have computed the time spent for one integral in each case assuming the time varies linearly with the number of integrals. The results are given in Table 1 (in this table, ,~ stands for analytical). Our conclusion is that on the CDC CYBER 170-750 we have used, the Gauss-Legendre method is economically interesting when the number of Gauss points is less than 4.

Table 1. Computation time - comparison between numerical and analytical procedures Constant element Linear element Number of Gauss Time for Time for Time for Time for Time for Time for points 1000 / 4000 f 1/ 1000 / 4000 / 1f J"

X

Figure 1. Boundary element C, current point {x, y) and observation point (x, y) Accepted February 1985. Discussion closes August 1985.

92

Engineering Analysis, 1985, Vol. 2, No. 2

1 2 3 4 5

0.264 0.396 0.511 0.617 0.684

0.768 1.68 E-4 1.206 2.70 E-4 1.595 3.60 E-4 2.010 4.60E-4 2.484 6.00 E-4

0.282 0.372 0.492 0.606 0.738

0.786 1.158 1.626 2.034 2.436

0.438

1.488

0.456

1 . 5 9 6 3.80E-4

3.50 E-4

1.68E-4 2.62 E-4 3.78 E-4 4.76E-4 5.66 E-4

0264-682X/85/020092-3 $2.00 © 1985 CML Publications

Numerical integration in the vicinity o f a logarithmic singularity: J.-L. Migeot A PRECISION TEST If A is the analytical value of an integral (with constant or linear shape function) and N is the Gauss-Legendre approxition, it is classical to say that the relative error: =

A--N

is a good measure of the approximation involved in the numerical process. However, e is not a good test in our particular case as it can be increased or decreased as much as wanted simply by changing the scale of the problem (i.e. multiplying all co-ordinates by a constant)) Therefore, we propose to use

Table 2. Number of Gauss points ensuring a given precision as a function of a relative distance between the observation point and the centre of the element Constant element

Linear element

Number R R R R of Gauss --(~,< 10-s) - - ( 7 < 10-6) --(~'< 10-s) --(~'< 10-6) points l l l l 4 3 2 1

0.95 1.4 3.5 65

1.2 2.1 6.5 210

0.90 1.5 5.8 (10 000)

1.2 2.4 12.5

A--N l as precision test as 3' is not sensitive to any scaling)

LINES OF NULL ERROR We compute 3' at several observation points near the element and then plot the iso-7 lines and in particular, those corresponding to 7 = 0. Figure 2 shows for instance 7 = 0 lines for a constant element and 4 Gauss points. Figure 3 gives the same result for a linear shape function. We observed that, for n Gauss points there exists curves (2n for constant element and 2 n - - 1 for linear element) along which 3' is equal to zero. It means that, for observation points lying on these curves, the value obtained using a n Gauss point integration scheme is equal to the analytical value although the computation time is much reduced. The preceding results are of particular interest in a BIEM program. Indeed, in this method, only the boundary values

of the unknown are obtained directly by solving the system of linear equations. To obtain displacements, potential, stresses, celerity or any other field values at a point lying inside the domain, integrals over the whole boundary (discretised in elements) are to be computed. It happens frequently that the stresses for instance are needed at a point very close to the boundary, and in this case, near a logarithmic singularity. It is here that the practical interest of the 7 = 0 curves appears. Indeed, often, one is interested to compute the field values at a given distance of the boundary. It is suggested to shift the point, parallel to the boundary element, such as to bring it on a 3' = 0 line. This method leads to an important economy of machine time without any loss of accuracy. Note that full advantage of the economy is gained only for constant elements as in BIEM model using linear elements, both integrals:

f lX lfn r d S ; c

Figure 2. Lines o f vanishing error (4 points rule, constant shape function)

1

(1-X)ln-dSr c

have to be evaluated on each element for each observation point and the lines of vanishing error are not the same for both integrals (they are symmetric with respect to the element's symmetry axis). However, even in that case, we should select our observation point on a line of vanishing error for one of the two integrals and then compute the other one analytically or with a Gauss-Legendre scheme of sufficiently high degree (see next paragraph). Tables giving the co-ordinates of several points of vanishing error are available from the author.

NUMBER OF GAUSS POINTS TO USE

Figure 3. Lines o f vanishing error (4 points rule, linear shape function)

We cannot always choose the position of the observation point, for instance when this observation point is a node of the boundary element grid. In this case, it is interesting to know how many Gauss points are necessary to ensure a given precision on the integration. Therefore, we have established the following table giving the number of Gauss points ensuring 3' < 10-s or 7 < 10-4 for a given relative distance of the observation point to the centre o f the element (see Fig. 4 and Table 2). The following formulae provide the same results but in a way more appropriated for practical applications, they give the number of Gauss points ensuring a given precision

Engineering Analysis, 1985, Vol. 2, No. 2

93

Numerical integration in the vicinity o f a logarithmic' singulariO,." J.-L. Migeot linear element int{--0.088747679207163 + 1.020406179125 z --3.724264231256 + 7.590005861177) 7 < 10 -6

constant element int{--0.040537147455563 + 0.5640260042862

R

--2.633101001036 + 7.417571945876}

linear elem en t int(--O.O 13827544482863 + 0.2876078407262 -- 1.729339911536 + 6.684946600066}

(int stands for nearest integer value).

ACKNOWLEDGMENT

Figure 4. Distance between the observation point and the centre o f the element

This research has been made possible by a fellowship from the Belgian Fonds National de la Re~zherche Scientifique (FNRS). Its support is gratefully acknowledged.

(7 < 10 -5 or "r < 10-6) as a function of the relative distance

REFERENCES

5 =R/l:

1 Brebbia, C. A., TeUes, J. C. and Wrobel, L. C. Boundary Element Techniques, Springer-Verlag, 1984 2 Migeot, J.-L. Resolution de differents problemes de champs par la mcthode des elements de frontieres. End of study thesis ULB Continumn mechanics - 1984 3 Migeot, J.-L. The effect of scaling on the boundary elemcnt solution of 2D potential problems. ULB Continuum mechanics Report (to be published)

3 ' < 10 -s

constant element int(--O.3090907265474fi 3 + 2.4928990820862 -- 6.78467219406 6 + 9.460603824453}

94

Engineering Analysis, 1985, Vol. 2, No. 2