Computer Implementation of a Coupled Boundary and Finite Element Methods for the Steady Exterior Oseen Problem

Abstract

We present a numerical technique based on the coupling of boundary and finite element methods for the steady Oseen equations in an unbounded plane domain. The present paper deals with the implementation of the coupled program in the two-dimensional case. Computational results are given for a particular problem which can be seen as a good test case for the accuracy of the method.

1. Introduction

The coupling of boundary and finite element methods has recently been shown to be a very effective tool for solving a certain class of physical problems with infinite (or even large) domains, for which the traditional numerical analysis techniques are unsuitable (cf. [1–13]). In particular, Sequeira considered the theoretical aspect of the coupled boundary and finite element methods for solving the steady exterior Stokes problem in its continuous and approximate forms in [4] and provided a computer implementation of the coupled boundary and finite element methods for solving the steady exterior Stokes problem in [5]; He and Mattheij have studied the theoretical aspect of the coupled boundary and finite element methods for solving the steady exterior Oseen problem in its continuous and approximate forms in [6]; furthermore, He has studied the posedness and convergence of the coupled boundary and finite element methods for solving the nonsteady exterior Oseen coupled problem in [7].

On this subject we have studied in detail the numerical solutions of the two-dimensional exterior Oseen equations for a steady-state incompressible viscous flow. Essentially, the coupling method involves the choice of an artificial smooth boundary separating an interior inhomogeneous region from an exterior homogeneous region. An integral equation over this interface, representing the solution in the exterior region in terms of a single-layer potential, is incorporated into a variational formulation in the primitive variables velocity-pressure for the interior region. This allows a discretization along the artificial boundary together with a typical discretization by finite elements to be employed. This paper is concerned with the implementation of these coupled boundary and finite element methods for the steady Oseen problem in a completely general form and without using the standard finite element software.

One of the difficulties encountered in assessing the performance of the algorithm for the approximate solution of the exterior Oseen equations is that of finding a suitable analytical solution with which a comparison may be made. The flow, due to an infinitely long circular cylinder, rotating uniformly about its axis in an infinite mass of viscous incompressible fluid, is one such solution in the two-dimensional case.

The major aim of the present work is the development of a computational program for our coupled boundary element-finite element methods. In this sense the study of the previous example may serve as a test case of the applicability of this technique to more complicated models. A series of numerical results demonstrates the accuracy of the method.

2. The Continuous Coupled Problem

Let be a simply-connected bounded domain in assumed to have a smooth boundary and let denote the complement of . The steady-state Oseen system of equations, governing the flow of a viscous incompressible fluid in an unbounded domain, may be reduced into the form (see [6]) where is the velocity, is the pressure, represents the density of body forces with a compact support in , is the dynamic viscosity of the flow, and is a constant vector. The boundary condition satisfies , being the unit outward normal to .

One of the main difficulties in solving is to deal with an infinite domain. Now, we introduce an artificial smooth boundary separating an exterior homogeneous region from an interior inhomogeneous region which contains the support of in . So, problem is decomposed into two joined problems: and a coupling condition on the interface where Then an appropriate coupling technique between boundary element and finite element may be used to solve .

Since the numerical analysis of this coupling procedure has already been developed, for the sake of brevity, only its essential features will be presented here. For additional mathematical details, the reader is referred to work [6]. Here we use the notations given in [6].

Let us introduce the Hilbert spaces with the standard norms.

We recall the mixed variational formulation of our problem described in [6]. It consists in finding such that with the conditions: and on , where the supplementary unknown appearing as the density of the single-layer potential, which represents the solution in the exterior domain, is identified with the local stress force of the flow; here denotes the unit outward (from ) normal to .

We have where denotes the duality pairing between and or between and , and is the trace operator of the zero-order such that . Here is a 2nd order tensor which satisfies where represents the coordinate axis, and is Dirac delta function. Recalling [7, 8], we know that where is essentially the zero-order Hankel function of the first kind. Here

Remark 2.1. It is worth noting that the pressure in , arising from the solution of our problem in the bounded subdomain , is determined up to an additive constant and therefore is unique up to an additive vector proportional to the normal to . However, by appropriately assembling the two problems in and , the involved constant must be determined. To overcome this difficulty in practical terms, the strategy adopted was to impose an extra coupling condition

We must point out, in contrast to other methods, when the Oseen problem is formulated in this manner, the stress force distribution, which is normally a quantity of interest in such calculations, is determined directly, and the accurate results are shown in Section 6.

3. The Approximate Coupled Problem

For the numerical approximation of our problem, we construct and study a finite element method based on the mixed variational formulation developed in Section 3. It involves Lagrangian finite elements which are conforming both in velocity and pressure but where the incompressibility condition is poorly approximated.

For simplicity we restrict here the discussion to the case where has polygonal boundaries, but the results can be easily extended to a general curved domain, by introducing an approximate boundary . For further details we also refer to [4, 5].

From now on, will be a real positive parameter tending to 0. We introduce three finite-dimensional spaces , , and such that We define In addition, we introduce the subspace of given by

With these spaces, problem is approximated by the following. for all , with the conditions: and on , where is an approximation of on .

We now construct finite dimensional spaces , , and such that the assumptions given in [6] are satisfied. Let be a polygonal domain and be a uniformly regular family of triangulations of made of triangles with no more than one side on , whose diameters are bounded by . We suppose that the family parameterized by the mesh size is affine of class , and uniformly regular as tends to zero, in the sense of [14], namely where and are two positive constants and

We choose the following finite element spaces (refer to [4, 5]): where is the polynomial in two real variables of degree and , () are the straightline segments of the interfacial boundary .

Recalling [6], there hold the following optimal error estimates: where is a positive constant.

4. Linear System

Here some statements are borrowed from Sequeira [5] with some changes. Let us recall the regular triangulation of made up of triangles with no more than one side located on and the finite element spaces defined by . and being spaces of the continuous functions which are quadratic and affine in each , respectively, it is natural to choose, as degrees of freedom, for a function its values at the vertices and midpoints of the triangulation and for its values at the vertices of , with the constraint of the interelement continuity. On the other hand, is the space of continuous and piecewise linear functions on with . We can construct its basis as follows: let us denote by the straightline segments of the boundary (of equal length , to simplify ), and by the set of continuous functions on , such that , and ; more precisely, with a parametric representation of each (for , , ), we define Then is generated by the set of functions , continuous on and piecewise affine, satisfying Now, if denotes the total number of nodal points (vertices and midsides ) of , the solution of the approximate problem reduces to the solution of a linear system of order . It takes the following matrix form: with and where the arising matrix may be written symbolically in a partioned form, exhibiting its natural block structure with Here is a sparse submatrix of finite element type and its coefficients, dealing only with the velocity and the pressure, are evaluated by numerical quadratures in each finite element. is produced by the boundary element treatment of the infinite subdomain: it is a full, symmetric submatrix and the singularities of its integrals are removed by exact calculation. It contains one supplementary row and column due to the discretization of the coupling condition . Finally, the rectangular submatrices and of the coupling coefficients only involve the degrees of freedom on the interfacial boundary nodes connecting to the finite element mesh, and are also sparse.

To summarize, the computational structure of the coupled system is very different and this leads to some difficulties on solving it. As can be expected, the finite element system is typically large but sparse and the boundary element system is small but dense. It is therefore of interest to design solution methods that exploit these attributes to maximum advantage.

Before proceeding, let us give a short look at the numerical effort needed to derive the boundary interface nodal coefficients to be assembled. As with finite elements, a global numbering system is used for these nodes.

We start with the boundary element terms in order to obtain submatrix . They are of the following form: Taking into account the above definition of the finite element space and the construction of its basis functions -, it is easy to see that, putting we must actually calculate and use the relations To write the above coefficient , we need to add the contribution from two adjoining elements and , , and . We define which become, in terms of parameters of the straightline segments, where , and .

Computing these integrals and using the relationship it is a simple matter to derive the full matrix (refer to [5] for the more details).

Let us now examine the coupling terms where , . According to -, we can write

5. Resolution of the Matrix System

As it has been noted previously, once the elemental matrix calculations have taken into consideration all the internal and boundary interfacial nodes, ensuring compatibly between the finite and boundary element meshes, the coupled analysis is carried out as in the standard finite element process. Of course, the global assembly and solution procedure we have used do not ignore the large zero blocks that arise, in order to increase the computational efficiency of this method.

We recall that problem produces a global system of linear equations which can be represented in a condensed matrix notation as and, more explicitly, in terms of the block structure of matrix where the solution vector is such that represents the velocity and pressure at all nodes and represents the interface unknown .

is a nonsingular matrix; can be expressed as Hence, we need to compute and by solving the matrix equations Then we solve the equation by a Gaussian elimination algorithm with partial pivoting and obtain from (5.2) that

6. Numerical Tests

The performance of the numerical model described above has been tested on the traditional example of the Oseen flow past a rotating infinitely long circular cylinder of the radius . This numerical example is an extension of the performance of the numerical model describing the traditional example of the Stokes flow past a rotating infinitely long circular cylinder of the radius ; see Sequeira [5]. For a specific application, an exact solution of equations may be sought as follows: assuming, to simplify, that the stream function associated with the flow only depends on the distance to the origin (taken in the axis of the cylinder in an -cartesian plane normal to it), we write the velocity distribution is then and the boundary condition for a cylinder of the radius . And we take . Then the relation between the pressure in the fluid and the external force acting on the cylinder is, in this case, given by since . Then, we may write corresponding to The interfacial boundary introduced in the fluid region is a circle of radius . The bounded domain is then a circular annulus limited by and .

From the expressions and for and , we can derive the local stress force of the flow at the interfacial boundary

Computations were carried out with and in four basic meshes of decreasing size , , , and (see Table 1), consisting of concentric layers of finite elements, and the finest finite and boundary element mesh is shown in Table 1.

Table 1 

Representative parameters of four basic meshes.

In Tables 2 and 3, we give a summary of computational results for and , respectively.

Table 2 

Relative error of numerical solution on four basic meshes ().

Table 3 

Relative error of numerical solution on four basic meshes ().

From Tables 2 and 3, one may observe that for the finer meshes we can obtain the better accuracy for the approximate velocity , pressure , and local stress force of the fluid flow. Hence, the algorithm is actually well behaved for the velocity, the pressure, and the local stress force of the fluid flow. This is confirmed by the accuracy of the obtained results.

Acknowledgments

This work is supported by the Natural Science Foundation of China (no. 10971166) and the National Basic Research Program (no. 2005CB321703) and Sichuan Science and Technology project (no. 05GG006-006-2).