The boundary element method applied to incompressible viscous fluid flow

TECHNICAL PAPERS

The boundary element method applied to incompressible viscous fluid flow

R. G. R. Camacho; J. R. Barbosa

Instituto Tecnológico de Aeronáutica; 12228-900 São José dos Campos, SP. Brazil; [e-mail]; [e-mail]; [e-mail]

ABSTRACT

An Integral equation formulation for steady flow of a viscous fluid is presented based on the boundary element method. The continuity, Navier-Stokes and energy equations are used for calculation of the flow field. The governing differential equations, in terms of primitive variables, are derived using velocity-pressure-temperature. The calculation of fundamental solutions and solutions tensor is showed. Applications to simple flow cases, such as the driven cavity, step, deep cavity and channel of multiple obstacles are presented. Convergence difficulties are indicated, which have limited the applications to flows of low Reynolds numbers.

Keywords: Boundary elements, fundamental solution, integral equations

Introduction

The need of solution of the system of partial differential equations which model the flow of a fluid in channels such as pipes, blade passages, nozzles and others appeared the very first day the fluid flow was modeled. The difficulties involved in obtaining closed solutions, even for very simple flows, required the development of clever techniques, but only with the application of numerical solutions to that system of equations, some flows of practical interest were calculated.

Several computational techniques have been used: finite difference, finite element, finite volume and boundary element to name the most known. As new algorithms were discovered and faster computers were produced, each of those methods evolved in all areas in the past years. Finite difference methods have been, implemented to solve flow problems. Finite elements gained attention in the past decades; in the seventies it was still crawling. Both are bases for commercial codes for the solution of flows of almost every kind. Computer effort has been limiting the application of the numerical methods in the sense that every new discovered method of solution claims reduction in CPU time and storage requirements.

The boundary element method, nevertheless, has progressed differently depending on the areas where it has been applied. It has been developed fastest in areas related to solid mechanics and acoustics problems, (Brebbia and Walker 1980; Brebbia, Telles and Wrobel 1984; Banerjee and Butterfield 1981) and slowest in the fluid mechanics.

A didactical approach is used in this work. The method of the boundary elements is applied to fluid problems, aiming also at introducing the methodology to new users. The computational implementation is based on the Kakuda and Tosaka (1988) reports. There, the boundary element method uses a reformulation of the unsteady Navier Stokes equations in terms of velocity components only, by making use of the penalty function method, an approach successfully applied to flow analysis with finite element. The effectiveness of this method was illustrated by several numerical examples. Tosaka and Onishi (1985, 1986) proposed new integral representations for the Navier Stokes equations for both steady and unsteady flow problems. The workability and validity of the methodology developed therein were shown with several numerical results for steady problems (Tosaka, Kakuda and Onishi (1985); Tosaka and Kakuda (1986); Tosaka (1986)).

Although integral methods were available many decades ago for the application to flow problems of practical interest, a comprehensive study of the formulation and application to flow problems are still being considered more recently, as they are expected to alleviate sensibly the storage and hopefully CPU time, Despite this apparent advantage, requiring less computational effort when volume integrals are transformed into surface integrals, some disadvantages arise, such as higher mathematical complexity in order to get an usable computational formulation; the need for the calculation of singular integrals; dense matrices whose inversion is more time consuming when comparable with the banded matrices in the finite difference and finite element schemes.

In the section Application below the application of the boundary element method to the following fluid problems are shown: a) stepped channel, b) box with moving lid, c) channel flow with multiple obstacles, d) deep cavity flow e) channel flow with heat transfer.

Nomenclature

r = density

µ = dynamic viscosity

P = pressure

T = temperature

n = absolute viscosity

k = thermal conductivity

c_n = specific heat at constant volume

c_p = specific heat at constant pressure

Statement of the Problem

Let W be a domain in R² and G its closed boundary; n_i be the outer normal vector to the boundary; the fluid be a perfect gas, incompressible and viscous; (x,y) be a point of W in R². The steady state conservation equations in cartesian co-ordinates can be written as:

Conservation of Mass:

Conservation of Momentum:

x - direction

y-direction

Conservation of Energy:

Let the following change of variables take effect in the conservation equations:

where ¥ refers to the far stream condition; and Re and Pr are the Reynolds and Prandtl numbers, respectively.

Then, the conservation equations become:

Conservation of Mass

Conservation of Momentum

x-direction

y-direction

Conservation of Energy

For the sake of simplicity, the asterisk will be dropped in what follows.

The independent variables of the problem are u, n, T and P. It is possible to rewrite the conservation equations in matrix form as

where [L] is a linear partial differential operator, {U}={u n T P}^T is the vector of the unknowns and {B} the vector of nonlinear convected terms. Depending on the assumptions made, [L] and {B} can take different forms. For instance, vector {B} can be linearised and the linear terms included in [L].

Let, for the moment, all non-linear terms are included into {B}. Then

Where:

The Method

Equation (7) has no known solution. Let be an approximate solution in the sense of , that is, differs from very little but it is not equal to .

In order to derive the integral equation for steady-state problem, the method of weighted residual is applied. The weighted residual statement for Eq. (7) can be expressed as

A possible solution can be obtained provided the W_IK as an appropriate weight function. It will be shown later that W_IK is chosen as the fundamental solution tensor for the adjoint of L_IJ.

Hörmander's (1965) (Banerjee, 1994) method is used for the calculation of the weight function W_IK tensor and fundamental solution. Although it does not provide W_IK directly, it allows, as a first step, the combination of several partial differential operators L_IJinto a single differential operator, from which the tensor W_IKis calculated. The weight tensor W_IKor the fundamental solution may be determined as a solution of steady Stokes problem with heat transfer:

where d(x-y) is the Dirac delta function and is the adjoint operator of .

Hörmander's method is simultaneously applied to the continuity, Navier Stokes and energy equations for steady, incompressible flow:

Multiplication of equation (12), to the left, by results:

where, and

whose terms are ; m_ij are the minors of [L], and

Thus,

Let

then

whose solution f^* is:

where denotes the distance between x and y (Tosaka, 1989).

Therefore, the fundamental solution tensor W_JK can be determined explicitly from equations (16) and (18) in conjunction with (19) as follows:

Important to notice that, the tensor W_JKcalculation is only determined analytically, this way, it is important to be careful in the obtaining of these equations.

Discretization

Let the Green-Gauss theorem be applied to equation (10) so that the domain integral is transformed into an integral over the contour G, divided into n_e boundary elements. Then,

which tells that the system of differential equations has been transformed into a system of algebraic equations (21) that involves the values of the variables at each boundary element. If one finds the values of the variables at the elements in the boundary, the solution in the boundary is then obtained.

After the application of the Green-Gauss theorem and integrating by parts over W, one arrives at the following equation that holds for every boundary element as well:

where summation is implied by repeating indices.

It is work mentioning that right hand side of equation (22) comprises integrals over the boundary and over the domain, these due to the non-linear convective terms, B_I.

In equation (22), C_IK is the tensor coefficient dependent on the geometry of the boundary. Its value is ½, 1 or 0, provided the point y lies over a locally regular boundary, within the domain or outside the boundary, respectively. If y lies at a corner of G its value is a/2p, where a is the angle formed by the left and right tangents to G. Also,

where coma (,) stands for derivative with respect to the following index and summation is implied by repeating indices.

From equation (24) the values of S_IKare calculated:

For constant boundary element, one has

Substituting the indicated expressions into equation (22):

Equation (28) can be rewritten after substitution for the constant terms listed in equation (27), from what results a system of algebraic equations.

Numeric Implementation

Boundary: Let G be the boundary, divided into m constant elements, with the collocation points (nodes) located at mid position of each element. Application of equation (28) to m nodes gives a set of 4m equations with 4m unknowns. For the solution of this system of equations two auxiliary matrices are assembled, for each element:

from which

Integrals (31) and (32) are carried out numerically, using one-dimensional Gaussian quadrature, if X¹ Y. When X=Y, the integrands of (31) and (32) become singular, requiring the calculation in the sense of Cauchy principal value. Among several techniques available to perform these calculations, in this work the method of Telles (1987) was chosen.

Domain: To calculate the integrals over W, the domain is subdivided into M elements by an appropriate net. Triangular cells will be used in this work. Let w_j be the Gauss weight function at point j, S_We the area of element e, and J the number of Gauss integration points. Then

Gauss quadrature with seven integration points in each triangular cell of the sub domain, and the Hammer technique, as described by Partridge et al (1992), are used to determine the domain integral.

In equation (28) the values u_i and T are known; t_i and q_i are unknown gradients of velocity and temperature.

Defining

equation (29) can be rewritten as

Defining

one arrives at

Boundary conditions: For the application of the boundary conditions to equation (37), it is worth noting that elements of d and of t have some prescribed values. It is therefore convenient to rearrange d and t in such a way that the unknowns come first and then the prescribed values, that is,

Rearrangement of matrices , [G] and {D} accordingly, results in

Rearrangement of equation (39) such that only the unknowns are in the left hand side of the equation, gives

Matrix [A] is dense so that inversion is time consuming. The inversion is carried out using the Gauss elimination algorithm. Solution of equation (40) gives the values of the unknowns at the boundary.

Computer Program: A modular computer program has been developed that is able to handle geometries composed of rectangles, written in FORTRAN and run in a 2.0 GHz personal computer.

The computer program implementation was carried out with the following steps:

- Definition of the geometry by a combination of rectangles. - Boundary discretization using elements of same size. - Domain grid generation using triangular elements. - Numbering elements counterclockwise. - Imposition of the boundary conditions and initialization of domain variables (velocity, temperature and pressure) using reasonable guesses according to the problem being solved. - Assembly of matrices and for each element e. - Assembly of matrices [G] and [H] for all elements on the boundary. - Numeric evaluation of domain integrals (equation (33)). - Solution of equation (40) for the determination of the variables at the boundary. - Solution of equation (22) at internal nodes, with C_KI=1.

More details about computational implementation can be found in the Santos (1998) and Ramirez et al (2004)

Application

For the demonstration of applicability of the method, five problems were chosen:

a) the recirculating flow in a square cavity driven by a lid sliding at uniform velocity;

b) the flow facing a forward step;

c) the flow over a deep cavity and

d) the flow over a deep cavity with the upper surface at a higher temperature;

e) the flow in a channel with multiple obstacles.

Driven cavity flow: The flow in the box is depicted using streamlines, as shown in Fig. 1. The boundary conditions are the no-slip in the box boundaries, that is, zero at the non moving surfaces and the velocity of the moving slid at the upper surface. Constant temperature was set on the boundary. Grid for Fig. 1a is 40x40 and for Fig. 1b is 30x30. Criteria of convergence were based on the difference e of the previous and the actual calculated values for velocities, pressure and temperature. Convergence was achieved up to Reynolds number of 400. Recirculation is detected at the bottom-right of the cavity.

Flow in a stepped channel: Streamlines of the flow in a forward facing step is shown in Fig. 2a. Boundary conditions are: parabolic distribution of velocities at inlet, no-slip condition on the walls and constant wall temperature. The results shown are for a grid of 26x30. The predicted reattachment point is in agreement with other predicted numeric methods such as the finite volume methods (Fig 2.b), (Rocamora F, 2002)

Deep cavity flow: The streamlines for Re=10, grid 40´40 and e=0.0001, are shown in Fig. 3a, which is a good result. In Fig. 3b one can also see satisfactory results with a moderately refined mesh with a grid of 20´20 and e=0.00001, so that it is possible to obtain good results in shorter processing time. For methodology validation, the finite volume method was used 40´40 elements with refined grid in the corner regions (Figure 3c). The CPU time for the method of boundary elements is smaller for the boundary element method. In Fig. 4 are shown the velocities at the middle of the cavity, obtained from BEM and FVM. It is important to notice that the Boundary Element Method with a grid of 20´20 gives good result when compared to the Method of Finite Volumes. Since convergence is achieved only for low Reynolds number, this formulation is not appropriate for the study of flow in turbomachines, one of our goals. Therefore, research is being carried out order to linearize {B} (Eq.10) and to incorporate those non-linear terms to the linear operator [L]. This method may become very attractive since it is expected to sensible reduce computational cost.

Deep cavity flow (heated upper surface): Although the energy equation had been used in the four previous applications, they were constant temperature applications. In this example, it is shown the variation of fluid temperature in Fig. 5, that shows the temperature contours when the top surface temperature is higher than the temperature of the other surfaces. Again, good results were obtained and no additional CPU time, compared to the previous application, was required.

Channel flow with multiple obstacles:Figures 6 show the streamlines of the flow in a channel with multiple obstacles, for Reynolds numbers 20.

Conclusion

The boundary element method can be applied to the calculation of incompressible viscous flows as demonstrated above. Although coarse grids were used, the results are quite satisfactory, capturing regions of reverse flows. Comparison of CPU times for the BEM and for FVM indicates that the BEM calculations are much faster. It is therefore important to investigate other applications, such as the flow in blade passages, which is the ultimate goal for the research under way. Flow in blade passages require a method that converges for much higher Reynolds number. Such a method is being developed, using linearization of the convective terms to allow the convective information to be incorporated into the fundamental solution and therefore achieve convergence for higher Reynolds numbers.

Acknowledgements

The authors gratefully acknowledge the financial support from FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo), El Paso Reference Center for Gas Turbines and ANEEL (Brazilian federal law 9991/2000).

Paper accepted July, 2005.

Technical Editor: Aristeu da Silveira Neto.

• Banerjee, P.K., 1994, "The Boundary Element Methods in Engineering", McGraw-Hill book Company Europe.
• Banerjee, P.K., Butterfield, R. 1981 "Boundary Element Methods in Engineering Science", McGraw-Hill; New York.
• Brebbia C. A., Telles T. C. F., Wrobel, L.C., 1984, "Boundary Element Techniques. Berlin. Heidelberg, New York: Springer.
• Brebbia, C. A., Walker. S., 1980, Boundary Element Techniques in Engineering. Newnes-Butterworths.
• Hörmander, L., 1965, "Linear Partial Differential Operators", Springer-Verlag, Berlin.
• Kakuda, K., Tosaka, N., 1988, "The Generalized Boundary Element Approach For Viscous Fluid Flow Problems", Boundary Element Methods in Applied Mechanics, Editors: M. Tanaka & T. A. Cruse, Pergamon Press, pp 305-314.
• Partridge P.W, Brebbia, C.A., Wrobell L. C., 1992, "The Dual Reciprocity Boundary Element Method", Computational Mechanics Publications, Southampton.
• Ramirez R.G. C., Barbosa J. R., 2004, "The Boundary Elements Method Applied to Incompressible Viscous Fluid Flow Problems", Anais em CD, III Congresso Nacional de Engenharia Mecânica, Belem Para.
• Rocamora, Jr F. D, De Lemos M J S,. "Heat Transfer Characteristic of Laminar and Turbulent Flows Past a Back Step in a Channel with a Porous Insert", 9 th Brazilian Congress of Thermal Engineering and Sciences, Caxambu-Minas Gerais, Brasil, 2002.
• Santos, Maria de Fatima de Castro Lacaz.,1998, "O Método de Elementos de Contorno Aplicado a Problemas de Escoamento de Fluidos". Tese Doutorado (Engenharia Aeronáutica e Mecânica) Instituto Tecnológico da Aeronáutica, (Orientador) J. R. Barbosa.
• Telles, J.C.F.A., 1987, "A Self-Adaptive Co-Ordinate Trans-formation for Efficient Numerical Evaluation of General Boundary Element Integrals, Int. J. Numeric Methods in Engineering, 24: pp 956-973.
• Tosaka, N. Kakuda, K. & Onishi, K., 1985, "Boundary Element Analysis of Steady Viscous Flows Based On P-V-U Formulation", Proc. 7nth Conference of BEM in Engineering, Lake Como, 71-80(9), Italy.
• Tosaka, N., 1986, "Numerical Methods for Viscous Flow Problems Using an Integral Equation . In: Wang, S.Y., Shen, H.W., Ding, L. Z. (eds): River sedimentation, pp. 1514-1525: The University of Mississippi.
• Tosaka, N., Kakuda, K. 1986, "Numerical Solution of Steady incompressible Viscous Flow Problems by the Integral Equation Method". In: Shaw. R. P., Periaux, J., Chaudouet, A., Wu, J., Marino, C., Brebbia, C.A. (eds): Innovative Numerical Method in Engineering, pp. 211-222. Berlin, Heidelberg, New York: Springer.
• Tosaka, N., Onishi, K., 1985, "Boundary Integral Equation Formulation for Steady Navier Stokes Equations using the Stokes fundamental Solutions". Eng. Analysis 2, 128-132.
• Tosaka, N., Onishi, K., 1986, "Integral Equation Method for Thermal Fluid Problems". In: Yagawa. G., Atluri, S. N. (eds): Computational Mechanics 86. vol. 2, pp. XI-103-XI-108. Berlin, Heidelberg, New York: Springer.

Publication Dates

• Publication in this collection
  02 Jan 2006
• Date of issue
  Dec 2005

History

• Accepted
  July 2005
• Received
  July 2005