Aspects About the Evaluation of the Singularities When Applying the Boundary Element Method

1

Bulletin of the Transilvania University of Braşov • Vol 13(48) - 2006

Aspects about the evaluation of the singularities when applying the Boundary Element Method to solve problems of fluid flow around bodies

L. GRECU[*]

Abstract One of the most powerful numerical techniques used to solve boundary values problems for fluid mechanics is the boundary element method. When solving a problem with this method there have to be made two important steps: first, we must obtain an equivalent boundary formulation for the problem involved, in fact a boundary integral equation, and then, this boundary integral equation which is a singular one must be solved. The Boundary Element Method reduces the problem to a system of linear equations. The diagonal coefficients of the system the problem is reduced at come from singular integrals and in order to obtain a well conditioned problem their evaluation has to be made with great attention. The paper presents some methods of treating the singularities that appear when the boundary element method with different types of isoparametric boundary elements is used to solve 2 and 3-dimensionall problems of fluid flow around bodies. Based on the numerical results obtained a comparison study between different methods is made and it points out the advantages of some methods on others. For solving some of the problems presented there are used constant, linear and quadratic boundary elements.

1. INTRODUCTION

The theory of the Boundary Integral Equation for the bidimensional and three-dimensional subsonic ideal compressible fluid flow around the bodies was elaborated by Acad. Lazăr Dragoş in some papers, for example: [4], [5],[6], [7], [8].

Solutions of the boundary integral equations obtained when applying the indirect method with sources distribution, of unknown intensities f, for solving the mentioned problems, are made in [8], [11], [12], [13] using constant linear and quadratic isoparametric boundary elements and, each time, the problems are reduced to linear systems of equations. The boundary integral equations obtained are singular ones and after there discretization some of the coefficients that appear arise from singular boundary integrals. That's why one of the most difficult and important stage in solving them is the evaluation of the coefficient involved, in fact the evaluation of the singular boundary integrals that appear. Their evaluation needs a special attention and has a great practicalimportance, because they lead to the dominant and concentrated near the diagonal matrix coefficients, and they are very important for a well-conditioned behavior of the matrix.

In the following sections we'll present aspects about some techniques used to evaluate the singular integrals that appeared in the mentioned problems.

2. The bidimensional case

For the bidimensional steady subsonic ideal compressible fluid flow around a body the boundary integral equation has the form (see [8]):

(1),

where are the components of the normal unit vector outward the fluid (inward the body) in the point ,(for the subsonic flow, M= Mach number), and f is the unknown function, the intensitie of the sources, persumed to satisfay a hölder condition. The sign " ' " denotes the Cauchy principal value of the integral.

2.1 Linear isoparametric boundary elements

In order to reduce the integral equation (1) to an algebraic system, in paper ... we approximate the contour C, boundary that is considered to be smooth and closed, with a polygonal line having the segments . The extremes of the segment are noted in a local numbering.

For equation (1) becomes:

(2) .

Using a local system of coordinates which has the origin in the first node of an element we have:

, (3),

where are the form functions given by: (4)

Using isoparametric boundary elements we have, for the unknown f, the local representation:

where are the nodal values of the unknown, it means the values of f at the extremes of the boundary element , in the local numbering. The equation (2) becomes:

(5)

where

and (6)

We have an algebraic system of N equations each of them of the following form:

, where , and , (7)

The coefficients can be computed analytically if they do not have singularities. If we calculate the singular integrals occuring in (5) or (7) we can use the result of the nonsingular ones and a limit process, finding so their Cauchy principal values. So all the coefficients have expressions depending only on the coordinates of the nodes used for the boundary discretization and it can be solved with a computer code. These results with a complet set of formulas of the coefficients involved can be found in paper [11].

After solving the system we deduce the expressions of the components of the velocity field and after that the value of the local pressure coefficient.

A numerical solution obtained with a computer code made with C++, based on the method presented, is compared in the graphic bellow with the exact solution that exists for the case of an incompressible ideal fluid flow around a circular obstacle (see [7]), and with the one obtained when constant boundary elements are use to solve the singular boundary equation (1).

Fig.1 The pressure coefficient (16 nodes) Fig.2 The pressure coefficient (20 nodes)

2.2 Quadratic isoparametric boundary elements

In paper [12] we use quadratic isoparametric boundary elements to solve the singular boundary integral equation (1). In that approach the geometry and the unknown, have local a quadratic variation, on each of the boundary elements obtained after the discretization of the boundary. For obtaining the discret equation the boundary is divided into N unidimensional quadratic boundary elements, each of them with three nodes: two extreme nodes and an interior one. Considering that the discrete equation is satisfied in every node, we obtain:

. (8)

For describing the geometry and the behavior of the unknown, on a boundary element, we use a quadratic model, with the same set of basic functions, noted . Using a local system of coordinates(the intrinsic system) we have:

(9)

where ,(10)

is a matrix with a single line, ,

,

the column matrices made with the global coordinates of the nodes of the boundary element noted, and the nodal values of the unknow function. There are used two systems of notation: a global and a local one (global- is the value of f for the node number j,-and local- is the value for the node number l of the element i). Returning to the global system of notation, we obtaine the following linear algebric system:

,.(11)

For getting the matrix we need to evaluate the integrals that appeare. One of them are usual integrals, but the other are singular integrals. Doing some calculous and using some notations we deduce the following expressions for the coefficients (see [12]).

.(12)

These coefficients depend on the following integrals: ,

where , and on other coefficients depending on the nodes coordinates.

For this integrals are nonsingular integrals, and they can be evaluated with an usual numerical method or with a mathsoft, using a computer.

For the singular integrals that appear we have applied three techniques: the method of truncation the interval, the Cauchy Principal Value Integral's Method, and a regularization method which uses new coordinates and new modified shape functions that allow the new integrands to have no or only weakly singularities, and so the use of the computer for evaluating them.

2.2.1. The method of truncation the interval

Using this method, the singularity that appeares is isolated, and the integral is calculated on an interval where the integrand has no singularities (see[1]). So, it becomes an usual integral, and it can be used the computer for its evaluation.

So, the singular integrals are evaluated as follows:

a) If , has a singularity for so it will be evaluated on , where is a verry small positive number.

2) If , the singularity appeares for and so it will be evaluated on .

3) If , the singularity appeares for , and so it will be evaluated on .

After solving the system (6), so after we find the values of f for the nodes choosen for the discretization of the boundary we may also compute the velocity for these nodes. Denoting by the components of the velocity we deduce:

,, where the coefficients depend on the nodes coordinates and on the same integrals as the system's coefficients (see [10]). For eample:

, and so for the others.

If we have nonsingular integrals. If ,, the integrals are singular, but they can be evaluated using the truncation method as in case of the singular coefficients of the system's matrix.

We compare the exact solution for the incompressible ideal fluid and a circular obstacle with the numerical one obtained with the method exposed. The computer code based on this method for this particular case is made in MATHCAD, and, as we see from the graphic below, it doesnt leed to a verry good concordation, so it is not recomanded. The error between the exact solution and the numerical one is quite big and it is bigger than that obtained when we use linear and constant boundary elements for solving the problem, even if, when we use quadratic boundary elements, the boundary, and so the function, is better aproximated, and we expect a smaller error. This error is most part, as we will see, due to evaluation of the singular coefficients. We have choosen for the value 0.09, and 20 nodes for the discretization of the boundary.

Fig.3. The truncation method

2.2.2 The Cauchy Principal Value Integral's Method

Another method we have applied, based on the definition of Cauchy principal value integral( see [8], [10]), allows us to obtain a less error between the numerical and the exact solution. We start applying this second method using the following expression for the coefficients:

. (13)

For example for, with singularity for -1, after some manipulations we finally deduce that:

So, only the second integral still has a weakly singularity, and only for its evaluation we use a method that isolates the singularity. We get:

(14)

For the other coefficients we obtain, in this case, the following expressions( see[10]):

(15)

For , doing like before, we have the following expressions for:

.

. (16)

For, analogous, we have:

(17)

The coefficients of the velocity's components are evaluated in a same manner. A program in MATHCAD, based on this method, applied for the mentioned particular case, brings better results. As we can see from the graphic below (Fig.4), the error between the numerical and the exact solution is smaller then that obtained when we use the truncation method(see Fig.3) :

Fig.4. The Cauchy Principal Value Integral's Method

2.2.3. The regularization method

The last method applied is the regularization method, and was inspired by the work of M Bonnet([2]). Using the Taylor series we can replace the shape functions by some expressions with Taylor polynomials and, after we simplify some factors, we get some modified shape functions, in fact, new integrands, that have no or only weakly singularities. For their evaluations we can use the computer and a math application, for example MATHCAD. We make a short presentation of this method because it is the one that leads to the best approximation.

Let , where the value of the local coordinate, that makes . We have .

Using the Taylor serie we get:

, (18)

where , depends on but also of , and it is named the modified shape function associated with , and it is given by: . (19)

Based on this technique we obtain the expressions for the modified shape functions. Denoting by we have: , , (20)

Using the new variables and the modified shape functions, the coefficients , given by the singular integrals, may be written, each of them, like a sum of two integrals: an usual one and an integral with a weakly singularity. For evaluating both kinds of integrals we can use a computer. We obtain: , where , are the notation for the usual integral and the singular one (the complete expressions of these coefficients can be found in [10]). So, the coefficients of system (11) depends only of the coordinates of the nodes and they can be evaluated with a computer.

In order to test the method we shall consider the same problem, the incompressible ideal fluid flow around a circular obstacle. Comparisons between the numerical solution obtained with a computer code in MATHCAD, based on this method, when we use for the discretization 20 nodes, and the exact solution, are performed in Fig.5. As we can see it shows a verry good agreement, because the calculated and analytical values are very close. The absolute value of the error is performed in Fig.6.

Fig.5. The regularization method Fig.6. The absolute error

3. The three-dimensional case

A boundary singular integral equation, formulated in velocity vector terms, is obtained in [7] using an indirect integral method for the three- dimensional potential subsonic flow around bodies. In order to solve the integral equation in papers [7], [13] ther are used constant and respectively triangular liniar elements for solving the singular integral equation.

The singular integral equation for the unknown has the same form as in (1) with the only difference that here the boundary is a surface, named , and so we have a bidimensional integral to handle with. For the incompressible case we have:

(21)

When solving this integral equation using linear boundary elements it is used a change of coordinates and a parametric representation as mathematical techniques for evaluating the singular integrals. The body surface, is divided into triangles, noted ; the extremes of the panels being situated on . Considering we have, for a fixed i, the integral equation:

(22)

If is one of the triangle corner points the integral calculated on has a singularity, otherwise it doesn’t. We focus our attention on the singular integrals. Using a local system of coordinates and denoting by the corner points (nodes) of a triangle, and by the intrinsic triangular coordinates (see for example [2], [3]) , (only two independent), and using the parametric representation: , where , , and , we get

.(23)

For we consider the linear interpolation .

(24)

Considering that the triangle, noted , has a corner point in we calculate the singular integrals occuring in (22) using the following relations:

, (25)

where are the other two nodes of , and are the values of the unknown function in these nodes; . If is the area of , we obtain:

. (26)

Using the fact that the finite part of is we get the equivalent form for the right hand side term :

(27)

Finally we get: ,where

, and . (28)

Using the global notation for the corner points of the panels involved (we return to the global system of coordinates), we obtain the equation in terms of the nodal unknowns , and finally we have an algebraic system which can be writtenas follows: . (29)

In order to test the method we had considered the uniform motion in the presence of a sphere of radius , centred in the origin of the system of coordinates. In this case the integral equation (21) can be solved analytically. A solution of this equation can be found in [9], using the spherical coordinates for the nodal points situated on the sphere () and the method of succesive approximations. The expresion for the exact solution is:

(30).

Comparisons between the analytical values of the intensity on the sphere and the values calculated by means of the boundary element method are performed in Fig.7., for the control points situated on the same circle on the sphere, so for a constant .

We can observe that the calculated and analytical values of the intensity are very close.

Fig.7. The sources intensities for the control points .

REFERENCES

1. 1 ANTIA H. M. Numerical Methods for Scientists and Engineers, Birkhausen, 2002.
2. BONNE M. Boundary integral equation methods for solids and fluids, John Wiley and Sons, 1995.
3. Brebbia, C.A., Walker, S.: Boundary Element Techniques in Engineering – Butterworths, London, 1980.
4. DRAGOŞ L. Boundary Element Methods (BEM) in the theory of thin airfoils. In: Rev. Roum. Math. Pures et Appl., 34(1989), 523.
5. DRAGOŞ L., DINU A. Application of the boundary element method to the thin airfoil theory. In: AIAA Journ., 28(1990),18-22.
6. DRAGOŞ L., DINU A. Application of the Boundary Integral Equations Method, to subsonic flow past bodies and wings. In: Acta Mechanica, 86(1991), 83.
7. Dragoş, L.: Mecanica Fluidelor Vol.1 Teoria Generală Fluidul Ideal Incompresibil (Fluid Mechanics Vol.1. General Theory The Ideal Incompressible Fluid) – Editura Academiei Române, Bucureşti, 1999.
8. Dagoş, L.: Metode Matematice în Aerodinamică (Mathematical Methods in Aerodinamics) – Editura Academiei Române, Bucureşti, 2000.
9. Dragoş, L., Marcov, N.: Solutions of the Boundary Integral Equations in the Theory of Incompressible Fluids – Analele Universităţii Bucureşti (Matematică – Informatică), XLVI-1997, p.111-119.
10. GRECU L. Metoda elementelor de frontieră aplicată în mecanica fluidelor, Bucureşti: Universitatea Bucureşti, 2004. Teză de doctorat.
11. GRECU L. A solution of the boundary integral equation of the theory of infinite span airfoil in subsonic flow with linear boundary elements. In: Analele Universităţii din Busureşti Matematică, Anul LII, Nr. 2(2003), pp. 181-188.
12. Grecu L.: A Solution of the Boundary Integral Equation of the 2D Fluid Flow around Bodies using Quadratic Isoparametric Boundary Elements”, The 13th Conference on Applied and Industrial Mathematics, Piteşti, 14-16 October, 2005.
13. Grecu L.: A Solution of the Boundary Integral Equation of the Three –Dimensional Airfoil in Subsonic Flow using Linear Boundary Elements”, The 12th Conference on Applied and Industrial Mathematics, Piteşti, 15-17 October, 2004, ROMAI Journal, Vol 1, Nr. 2, 2005, pg. 101-107.

[*] University of Craiova, Faculty of Engineering and Management of Technological Systems, Drobeta Turnu Severin,