The paper presents a solution of the problem of the 2D compressible fluid flow around obstacles, based on the boundary element method in which the singular boundary integral equation is solved with cubic boundary elements. In this approach the singular boundary integral with sources distribution is considered. The numerical method is implemented into a computer code developed under Mathcad, and numerical solutions for some particular cases are obtained. Numerical solutions are compared with analytical ones for cases when the latter exist. A good agreement between them can be observed even when the number of nodes chosen for the boundary discretization is quite small, the fact that shows the great efficiency of the exposed method.
MSC: 65N38, 76M15, 76G25, 35Q35.
Keywords:boundary element method; cubic boundary element; singular boundary integral equation; compressible fluid flow
Many of the phenomena studied in the nature are described by boundary value problems, in fact, by partial differential equations or systems of partial differential equations with boundary restrictions, and only few of them can be, in general, solved analytically.
One powerful numerical technique that can be used to solve such problems is the boundary element method. This method has been widely used to solve partial differential equations or systems of partial differential equations with boundary conditions, many books being dedicated to this subject (for example, [1-5]).
The boundary element method consists in solving, by discretization, a boundary integral equation, equivalent to the mathematical model of the problem to solve, by using finite elements named boundary elements. When applying this method, two big steps have to be done: first a boundary integral formulation is obtained, and usually it represents a singular boundary integral equation or a system of singular boundary integral equations, and then this singular integral equation is solved and numerical solutions are found.
When the problem to solve is described by systems of partial differential equations with fundamental solutions, this method can be successfully applied even if the boundary conditions are not linear, as in the case of the problem considered in this paper. The method is suitable for problems with infinite domains, as those of fluid flows around obstacles, because the fundamental solutions satisfy conditions that do not involve the presence of a fictive boundary at great distance. Other numerical methods such as finite difference method, finite volume method and finite element method can be applied to solve problems of fluid flow (see, for example, ), but the obstacle needs to be placed in a grid, and so the computational effort is quite big. The boundary element method has another great advantage over these, because it reduces the dimension of the problem by one and also the necessary computational effort.
In this paper, based on this method, a numerical solution for the two-dimensional problem of the compressible fluid flow around obstacles is obtained. The problem of the fluid flow around obstacles has been the subject of many papers, among which we mention [7-10], some of them dealing only with the incompressible case. The boundary element method has also been applied to solve this problem, but in this approach, a singular boundary integral equation formulated in velocity terms is considered, and the solution is evaluated by using higher-order boundary elements. So, the primary variables of interest are directly numerically evaluated without needing a differentiation process, as in the case of other approaches in which the stream function or the potential of the fluid flow are envisaged.
This paper is focused on the second step in applying this method, namely on solving the singular boundary equation which arises when at the first step an indirect technique with sources distribution is considered.
Briefly, the problem to solve is that of a subsonic uniform compressible fluid flow perturbed by the presence of an obstacle of boundary C assumed to be smooth and closed. The goal is to find out the perturbation produced by the obstacle and the fluid action on it. The mathematical model of the problem consists in a system of partial differential equations with a nonlinear boundary condition. In dimensionless variables, it can be written as (see , chap. 4)
where and v are the components of the perturbation velocity, , denote the components of the normal unit vector outward the fluid, (for the subsonic flow), and .
Assimilating the boundary with a continuous distribution of sources, of intensity f, and assuming that f satisfies a Hölder condition on C, in , chap. 4, a singular boundary integral equation is deduced. It has the following expression:
where , denote the components of the normal unit vector outward the fluid in the point . The sign ‘′’ denotes the Cauchy principal value of the integral.
This singular boundary integral is solved in the same paper by using a collocation method. Linear boundary elements are used in paper , and quadratic ones in paper , in order to get numerical solutions. The numerical solution acquires a higher degree of accuracy when higher-order boundary elements are used to solve the boundary integral equation.
Using cubic boundary elements to solve this singular boundary integral equation, we obtain a better approximation for the solution, because cubic boundary elements offer a better approximation not only for the unknown of the problem, but also for the geometry of the boundary involved.
Using the approximation models, the problem is reduced to a linear system of equations, which has as unknowns the nodal values of the sources intensities. After solving it, the components of the perturbation velocities and the local pressure coefficient are numerically evaluated.
A difficult stage in solving singular boundary integral equations is the treatment of the singularities that arise. Different methods can be applied to overpass this difficulty (see, for example, [14,15]). Numerical evaluation of integrals of singular kernels is a very important aspect for a well-conditioned behavior of the system matrix. In paper  a regularization technique with modified shape functions has been applied. In this approach, the coefficients arising from integrals of singular kernels are evaluated with a simple method, namely the truncation method.
2 The discretized equation
In the herein boundary element approach for solving the singular boundary integral equation (2), the boundary is divided into N cubic boundary elements, noted , . Each cubic boundary element has four nodes: two extremes and two interior ones. 3N nodes are used for the boundary discretization, because each boundary element has two common nodes with the two adjacent boundary elements.
So, we obtain the discretized form of equation (2):
Considering that the discretized equation is satisfied in every node , and denoting by , we get the following equation, in fact, we get 3N equations - the linear system the problem is reduced at
We further do not use the notation referring to the Cauchy principal value of an integral, but we take into account that when the integrals in (4) have singular kernels, and they are understood in this way, and when they are the usual ones.
From the family of cubic boundary elements, we choose the isoparametric ones. The cubic isoparametric boundary element uses the same set of basic functions (see [3,4]), noted , , , , for describing the geometry and the unknown function. If the homogeneous variable ξ is used, these functions have the following expressions, similar to those given in , pp.350-351:
For such a cubic boundary element , the following approximation models stand:
where is a line matrix, , , are the column matrices made with the global coordinates of nodes, and the nodal values of the unknown function corresponding to the same boundary element. Two systems of numbering are used: a global one ( , represents the nodal value corresponding to node number j) and a local one ( , , represents the nodal value of the lth node of element number j).
We introduce in equation (4) the approximation models (6), and we obtain
where we denote by the Jacobian coordinate transformation.
We can write (7) as
Further, we get
Returning to the global system of numbering, we obtain the following linear algebraic system:
3 Coefficients evaluation
In order to completely solve the problem and to obtain the numerical solution, we need to compute the coefficients of A, so we need to evaluate the integrals that appear. Some of them are usual integrals, and we can use any mathematical software to compute them, but the others are integrals of singular kernels. For the evaluation of the singular integrals, we use in this paper the truncation method, the method that brings good results when the integrand does not oscillate near the singularity (see , chap. 6), as in the following case.
Replacing relations (5) and (6) in (9), then doing some calculus and notations, we get the expressions for the functions in (11). For example, for , we get
Doing the same for , and , the following relations hold:
The above relations allow us to evaluate the coefficients given by relation (11) with a computer code, because their expressions depend only on the nodes coordinates. In order to find the expressions of coefficients arising from integrals considered in a Cauchy principal value sense, for which we apply the truncation technique, we need to mention the connection between the two systems of numbering used in this approach: the global and the local one. This also has to be done for simplifying the implementation of the method exposed into a computer code.
Table 1 shows the relations between the two systems of numbering.
Table 1. Relation between the local and the global system of numbering
For the jth boundary element , , the nodes have, in the global system of numbering, the numbers , , 3j, , understanding that node number is node number 1, the boundary being closed. So, the coefficients given by relations (11) come from singular integrals only when i takes one of the above values ; otherwise they come from the usual ones.
Using relations (13) and (14), we have, if , the following expressions for the nonsingular coefficients given by (11):
In order to use the truncation method for evaluating the singular coefficients, we consider a very small parameter to reduce the domain of integration by eliminating the singularity. Taking into account the cases when the singularity arises, we get the following four situations.
1. When , the singularity arises when . The coefficients have in this case the expressions
2. When , the singularity arises when . The coefficients are in this case given by
3. When , the singularity arises when . The coefficients are
4. Finally, when , the singularity arises when , and so we have
In order to obtain matrix A from relation (13), we return to the global system of numbering, and we write system (14) as follows:
Finally, the system to solve is
Solving system (25), we find the nodal values , and then we can evaluate the velocity components.
Under the same assumption, that f satisfies a Hölder condition on the boundary, the components of the perturbation velocity can be obtained. On the boundary, u, v have the following expressions (see , chap. 4):
Analogously, as in case of solving the boundary integral equation, we get the nodal expressions of the velocity components. After the boundary discretization, after introducing the approximation models (6) in the above relations, and after considering , , and doing some calculus, we obtain
Using the same notations as before, we have:
and further, the following relations:
The coefficients from the above relations are evaluated in the same manner as in the case of the matrix coefficients; analogous formulas with (19)-(22) stand in the case of singular integrals from (29) and (30).
Finally, the following expressions for the velocity field components on the boundary are deduced from (27) using the connection between the two systems of numbering, and returning to the global one:
Replacing relations (33) in (32) and then in (31), we evaluate the nodal values of the components of velocity and then the local pressure coefficient with the following expression (when ):
where γ is a fluid characteristic, a constant equal to the ratio between the specific heat at constant volume and the specific heat at constant pressure. For air it has the value .
4 Numerical results and conclusions
Using a good method for evaluating the singular integrals that appear is a stage of great practical importance, because the coefficients given by these singular integrals are dominants situated near the diagonal of the matrix, and so they play an important role in establishing a good behavior of the system.
The existence of particular cases for which the considered problem has an exact solution is of major significance because it allows us to test the developed method. In this paper we compare the numerical solution with the exact one in order to make an analytical checking of the exposed method. The analytical checking always represents a great challenge and is a very difficult step to overpass to validate the method.
We consider first the incompressible case and a circular obstacle. In , chap. 5, the exact solution for this case is presented. The components of the dimensionless velocity (u, v) on the boundary and the local pressure coefficient (cp) have the following expressions:
The presented method is implemented into a computer code made with Mathcad programming tools. It is used to find out the numerical solutions for the components of velocity and also for the local pressure coefficient and the exact values of the mentioned quantities.
The input data refer to the boundary geometry, the number of elements used for the boundary discretization, the values of eps, and Mach number. The output data are: the numerical and the exact nodal values of velocity components and of the local pressure coefficient. Graphical representations of these values are also performed.
We can choose different values for the boundary elements number N and for the truncation parameter eps. For making the comparison, we consider first 30 nodes, then 45 nodes, for the boundary discretization, and .
In Figure 1 the numerical values and the exact ones for the component of the velocity along Ox axis are represented. In Figure 2 a comparison between the numerical values of the component along Oy axis and the exact ones is made. For the nodal values of the local pressure coefficient, the comparison is made in Figure 3.
Figure 1. The nodal values of the velocity component alongOxaxis - exact (vx) and numerical (Ux) solution; circular obstacle; (a) 30 nodes; (b) 45 nodes.
Figure 2. The nodal values of the velocity component alongOyaxis - exact (vy) and numerical (Uy) solution; circular obstacle; (a) 30 nodes; (b) 45 nodes.
Figure 3. The nodal values of the local pressure coefficient - exact (cp) and numerical (Cp) solution; circular obstacle; (a) 30 nodes; (b) 45 nodes.
As we can see from the three graphics, we get a good accuracy for the numerical solution.
The computer code can deal with the case of compressible fluid flows, so numerical solutions for different values of Mach number can be obtained, cases when exact solutions are not known. In Figure 4 the numerical solution for the case of a circular obstacle, when 45 nodes are used for the boundary discretization and , is performed.
Figure 4. The numerical values of the velocity components (Ux,Uy) and the local pressure coefficient (Cp); circular obstacle; ; 45 nodes.
We test the method for the case of an elliptical obstacle too. For the incompressible fluid flow and an elliptical obstacle, the problem has an exact solution too, which can be found in , chap. 5. It is obtained with the complex potential of the movement F. If the boundary is , , , at zero incidence angle, F is given by
Thus, the dimensionless velocity components are given by the following relations:
We consider in this paper an elliptical obstacle with and .
Figure 5. The nodal values of the velocity component alongOxaxis - exact (vx) and numerical (Ux) solution; elliptical obstacle; 36 nodes.
Figure 6. The nodal values of the velocity component alongOyaxis - exact (vy) and numerical (Uy) solution; elliptical obstacle; 36 nodes.
Figure 7. The nodal values of the local pressure coefficient - exact (cp) and numerical (Cp) solution; elliptical obstacle; 36 nodes.
As in the case of a circular obstacle, we can see a great accuracy between the numerical solution and the exact one, even if using only 12 cubic boundary elements for the boundary discretization.
The comparisons have shown that the numerical values of the nodal components of the velocity and of the local pressure coefficient and the exact ones are in good agreement even for a small number of nodes used for the boundary discretization. The analytical checking validates the exposed method and also the computer code.
The effectiveness of the boundary element method is clearly dependent on the implementation of efficient and accurate integration procedures to evaluate boundary integrals of the singular kernels. The herein paper points out that the accuracy of the numerical solution is good in the case of using cubic boundary elements for the boundary discretization and the truncation method for singular integrals evaluation. Moreover, the truncation method is a method very easy to apply and to implement into a computer code.
Better results can be obtained when using more nodes for the boundary discretization, or a smaller value for eps parameter, but the computational effort is bigger when more nodes are used, and so, for obtaining a good ration accuracy/computational effort, we have to choose a certain range of values for the number of nodes.
Numerical integration is always a much more stable and precise process than numerical differentiation and, if we formulate integral equivalent representations using primary functions when applying boundary element method, no differentiation of numerical quantities is needed, which is the fact that influences the accuracy of the numerical solution.
All this recommend the boundary element method, with boundary formulation which implies primary variables, as an efficient numerical technique to solve boundary value problems.
The author declares that she has no competing interests.
Dedicated to Professor Hari M Srivastava.
The author thanks the reviewers for their valuable suggestions and useful comments, and thus for their help in the improvement of the original manuscript.
Carabineanu, A: A boundary element approach to the 2D potential flow problem around airfoils with cusped trailing edge. Comput. Methods Appl. Mech. Eng.. 129, 213–219 (1996). Publisher Full Text
Dragoş, L, Dinu, A: Application of the boundary element method to the thin airfoil theory. AIAA J.. 28, 1822 (1990). Publisher Full Text