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

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.


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 .

Mathematical Problems in Engineering
On this subject we have studied in detail the numerical solutions of the twodimensional 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.

The Continuous Coupled Problem
Let Ω 0 be a simply-connected bounded domain in R 2 assumed to have a smooth boundary Γ and let Ω denote the complement of Ω 0 ∪ Γ.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 u is the velocity, p is the pressure, f represents the density of body forces with a compact support in Ω, ν > 0 is the dynamic viscosity of the flow, and w 0 w 0 , 0 is a constant vector.The boundary condition satisfies Γ u 0 • n ds 0, n being the unit outward normal to Γ.
One of the main difficulties in solving O − S is to deal with an infinite domain.Now, we introduce an artificial smooth boundary Γ 1 separating an exterior homogeneous region Ω 2 Mathematical Problems in Engineering 3 from an interior inhomogeneous region Ω 1 which contains the support of f in Ω.So, problem O − S is decomposed into two joined problems: and a coupling condition on the interface ∂u i ∂x j ∂u j ∂x i , i,j 1, 2.

2.1
Then an appropriate coupling technique between boundary element and finite element may be used to solve O − S .
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 u, p, λ with the conditions: n, λ 0 and u u 0 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 n denotes the unit outward from Ω 2 normal to Γ 1 .
Here U k , P k is a 2nd order tensor which satisfies where e k represents the coordinate axis, and δ x is Dirac delta function.Recalling 7, 8 , we know that where is essentially the zero-order Hankel function of the first kind.Here γ Euler s constant 0.5772157.

2.8
Remark 2.1.It is worth noting that the pressure p ∈ L 2 Ω 1 in Q h , arising from the solution of our problem in the bounded subdomain Ω 1 , is determined up to an additive constant and therefore λ ∈ T is unique up to an additive vector cn proportional to the normal to Γ 1 .However, by appropriately assembling the two problems in Ω 1 and Ω 2 , 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.

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 Ω 1 has polygonal boundaries, but the results can be easily extended to a general curved domain, by introducing an approximate boundary Γ h ∪ Γ 1h .For further details we also refer to 4, 5 .
From now on, h will be a real positive parameter tending to 0. We introduce three finite-dimensional spaces X h , S h , and M h such that

Mathematical Problems in Engineering
We define In addition, we introduce the subspace W 0h of W h given by With these spaces, problem Q is approximated by the following.
, with the conditions: n, λ h 0 and u h u 0h on Γ h , where u 0h is an approximation of u 0 on Γ.
We now construct finite dimensional spaces X h , S h , and M h such that the assumptions H 1 -H 5 given in 6 are satisfied.Let Ω 1 be a polygonal domain and {τ h } be a uniformly regular family of triangulations of Ω 1 made of triangles K with no more than one side on ∂Ω 1 , whose diameters h K are bounded by h.We suppose that the family {τ h } parameterized by the mesh size h is affine of class C 0 , and uniformly regular as h tends to zero, in the sense of 14 , namely where α and β are two positive constants and ρ K sup{diameter of B; B is a circle contained in K}.

3.5
We choose the following finite element spaces refer to 4, 5 :

3.6
where P l is the polynomial in two real variables of degree ≤ l and s i , 1 ≤ i ≤ n |s i | ≤ h are the straightline segments of the interfacial boundary Γ 1 .

Mathematical Problems in Engineering 7
Recalling 6 , there hold the following optimal error estimates: where c is a positive constant.

Linear System
Here some statements are borrowed from Sequeira 5 with some changes.Let us recall the regular triangulation τ h of Ω 1 made up of triangles K with no more than one side located on ∂Ω 1 and the finite element spaces defined by 3.6 .X h and M h being spaces of the continuous functions which are quadratic and affine in each K ∈ τ h , respectively, it is natural to choose, as degrees of freedom, for a function ω h ∈ X h its values at the vertices and midpoints of the triangulation τ h and for M h its values at the vertices of τ h , with the constraint of the interelement continuity.On the other hand, T h S 2 h ∩ T is the space of continuous and piecewise linear functions μ h on Γ 1 with Γ 1 μ h ds 0. We can construct its basis as follows: let us denote by A i−1 A i Γ i 1 ≤ i ≤ n the n straightline segments of the boundary Γ 1 of equal length h, to simplify , and by {π i } 1 ≤ i ≤ n the set of continuous functions on Γ 1 , such that π i A i 1, π 1 A j 0 j / i and Γ 1 π i ds 1; more precisely, with a parametric representation of each

4.1
Then S h is generated by the set of n − 1 functions {t i }, continuous on Γ 1 and piecewise affine, satisfying and where the arising matrix M may be written symbolically in a partioned form, exhibiting its natural block structure Here A 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.D 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 2.9 .Finally, the rectangular submatrices C and C 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 D. They are of the following form: Taking into account the above definition of the finite element space T h and the construction of its basis functions 4.1 -4.2 , it is easy to see that, putting i, j 1, . . ., n − 1 we must actually calculate U 22 x − y π i x π j y ds x ds y , 4.9 and use the relations To write the above coefficient d k ij , we need to add the contribution from two adjoining elements Γ i 1 and Γ i 1 1 , Γ j 1 , and Γ j 1 1 .We define which become, in terms of parameters σ, τ ∈ 0, 1 of the straightline segments, where K x U 11 x , U 12 x , and U 22 x .Computing these integrals and using the relationship it is a simple matter to derive the full matrix D refer to 5 for the more details .

Mathematical Problems in Engineering
Let us now examine the coupling terms σ ij U k , P k x − y n j x u hi x μ hk y ds x ds y ,

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 Q h produces a global system of linear equations which can be represented in a condensed matrix notation as where the solution vector is such that X represents the velocity and pressure at all nodes and Y represents the interface unknown λ h .A is a nonsingular matrix; 5.2 can be expressed as

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 R.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 R; see Sequeira 5 .For a specific application, an exact solution of equations O − S 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 x 1 , x 2 -cartesian plane normal to it , we write ψ log r, with r x 2 1 x 2 2 ; 6.1 the velocity distribution is then and the boundary condition for a cylinder of the radius R.And we take w 0 1, 0 .Then the relation between the pressure in the fluid and the external force acting on the cylinder is, in this case, given by

6.6
The interfacial boundary introduced in the fluid region is a circle of radius R 2 .The bounded domain Ω 1 is then a circular annulus limited by R and R 2 .
From the expressions 6.2 and 6.5 for u and p, we can derive the local stress force of the flow at the interfacial boundary

6.7
Computations were carried out with R 1 and R 2 2 in four basic meshes of decreasing size h, 3/4h, 3/5h, and 1/2h see Table 1 , consisting of concentric layers of finite elements, and the finest finite and boundary element mesh is shown in Table 1.
In Tables 2 and 3, we give a summary of computational results for ν 1 and ν 0.5, respectively.
From Tables 2 and 3, one may observe that for the finer meshes we can obtain the better accuracy for the approximate velocity u h , pressure p h , and local stress force λ h 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.
explicitly, in terms of the block structure of matrix M

Hence, we need to compute A − 1 b
and A −1 C by solving the matrix equations the equation D − C C Y −C b, 5.6 by a Gaussian elimination algorithm with partial pivoting and obtain from 5.3 that X b − CY. 5.7
the total number of nodal points vertices N v and midsides N m of τ h , the solution of the approximate problem Q h reduces to the solution of a linear system of order 3N v 2N m 2 n − 1 1.It takes the following matrix form:

Table 1 :
Representative parameters of four basic meshes.NL: number of Fe-layers; N: number of points on each Fe-layer; NE: number of elements; NO: number of nodal points; NVP: number of degrees of freedom associated with the velocity and pressure; NT: total number of degrees of freedom.

Table 2 :
Relative error of numerical solution on four basic meshes ν 1 .

Table 3 :
Relative error of numerical solution on four basic meshes ν 0.5 .