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

Authors: Minfu Feng; Yinnian He; Ruiting Ren
Journal: Mathematical Problems in Engineering
 (2010)

---

## 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.


---

## Body

## 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Ω0 be a simply-connected bounded domain in R2 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])
(O-S)-νΔu+(w0·∇)u+∇p=finΩ,divu=0inΩ,u=u0onΓ,u(x)→0asx→∞,
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 w0=(w̅0,0) is a constant vector. The boundary condition satisfies ∫Γu0·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 from an interior inhomogeneous region Ω1 which contains the support of f in Ω. So, problem (O-S) is decomposed into two joined problems:
(O-S1)-νΔu+(w0·∇)u+∇p=finΩ1,divu=0inΩ1,u=u0onΓ,(O-S2)-νΔu+(w0·∇)u+∇p=0inΩ2,divu=0inΩ2,u(x)→0asx→∞,
and a coupling condition on the interface Γ1(O-S3)u-=u+,λ-=λ+,
where
(2.1)(u-,p-)=limx→Γ1(u(x),p(x))|Ω1,(u+,p+)=limx→Γ1(u(x),p(x))|Ω2,λ-=σ(u-,p-)·n̅|Γ1,λ+=σ(u+,p+)·n̅|Γ1,σ(u,p)=(-δijp+2νεij)2×2,εij=12(∂ui∂xj+∂uj∂xi),i,j=1,2.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(2.2)W={v∈H1(Ω1)2;v=0onΓ},M=L02(Ω1)={q∈L2(Ω1);∫Ω1q(x)dx=0},T=H0-1/2(Γ1)2,H0-1/2(Γ1)={t∈H-1/2(Γ1);∫Γ1tds=0},with the standard norms.We recall the mixed variational formulation of our problem described in [6]. It consists in finding (u,p,λ)∈H1(Ω1)2×T×M such that
(Q)a(u,v)+a1(w0,u,v)+〈γ0v,λ〉-(p,divv)=(f,v)∀v∈W,2b(λ,μ)+〈γ0u,μ〉-2〈G(γ0u),μ〉=0∀μ∈T,(q,divu)=0∀q∈M,
with the conditions: 〈n̅,λ〉=0 and u=u0 on Γ, where the supplementary unknown(2.3)λ=σ(u,p)·n̅∣Γ1appearing 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; heren̅ denotes the unit outward (from Ω2) normal to Γ1.We have(2.4)Skμ(x)=∫Γ1Uk(x-y)·μ(y)dsy∀μ∈H-1/2(Γ1)2,Gkv0(x)=∫Γ1v0(y)·σ(Uk,Pk)(x-y)·n̅(y)dsy+∫Γ1(v0(y)·Uk(x-y))(w0·n̅(y))dsy∀v0∈H1/2(Γ1)2,a(u,v)=2ν∑i,j=12∫Ω1εij(u)εij(v)dx,a1(u,v,w)=∫Ω1(u·∇)v·wdx,b(λ,μ)=〈Sλ,μ〉,〈G(γ0u),μ〉=∑k=12〈Gk(γ0u),μk〉,where〈·,·〉 denotes the duality pairing between H1/2(Γ1)2 and H-1/2(Γ1)2 or between H1/2(Γ1) and H-1/2(Γ1), and γ0 is the trace operator of the zero-order such that γ0u=u|Γ1. Here (Uk,Pk) is a 2nd order tensor which satisfies
(2.5)νΔUk(x)-(w0·∇)Uk(x)-∇Pk(x)=δ(x)ek,divUk(x)=0,k=1,2,whereek represents the coordinate axis, and δ(x) is Dirac delta function. Recalling [7, 8], we know that
(2.6)Uki=(δkiΔ-∂2∂xk∂xi)Φ(x),Pk(x)=xk2π|x|2,Φ(x)=-12πw̅0∫0x1{log(τ2+x22)1/2+exp{-w̅02ντ}H(|w̅0|ν(τ2+x22)1/2)}dτ,where(2.7)H(t)=-logt4∑m=0∞(t/4)2m(m!)2-γ+∑m=1∞(t/4)2m(m!)2[1+12+⋯+1m-γ]is essentially the zero-order Hankel function of the first kind. Here(2.8)γ=Euler′sconstant=0.5772157.Remark 2.1.
It is worth noting that the pressurep∈L2(Ω1) in (Qh), 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
(2.9)∫Γ1λ·n̅ds=0.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 Section6.


## 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 Section3. 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 Xh, Sh, and Mh such that
(3.1)Xh⊂H1(Ω1);Sh⊂H-1/2(Γ1);Mh⊂H1(Ω1).We define(3.2)Wh=Xh2∩W,Th=(Sh∩H01/2(Γ1))2⊂T.In addition, we introduce the subspaceW0h of Wh given by
(3.3)W0h={vh∈Wh;(qh,divvh)=0∀qh∈Mh}.With these spaces, problem (Q) is approximated by the following. 						(Qh)find(uh,λh,ph)∈Wh×Th×Mhsuchthata(uh,vh)+a1(w0,uh,vh)+〈γ0vh,λh〉-(ph,divvh)=(f,vh),2b(λh,μh)+〈γ0uh,μh〉-2〈G(γ0uh),μh〉=0,(qh,divuh)=0,
for all (vh,μh,qh)∈Wh×Th×Mh, with the conditions: 〈n̅,λh〉=0 and uh=u0h on Γh, where u0h is an approximation of u0 on Γ.We now construct finite dimensional spacesXh, Sh, and Mh such that the assumptions (H1)–(H5) 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 hK are bounded by h. We suppose that the family {τh} parameterized by the mesh size h is affine of class C0, and uniformly regular as h tends to zero, in the sense of [14], namely
(3.4)α̅h≤hK≤β̅ρK,∀K∈τh,whereα̅ and β̅ are two positive constants and
(3.5)ρK=sup{diameterofB;BisacirclecontainedinK}.We choose the following finite element spaces (refer to [4, 5]):
(3.6)Xh={ωh∈C0(Ω¯1);ωh∣K∈P2,∀K∈τh}⊂H1(Ω1),Sh={νh∈C0(Γ1);νh∣si∈P1,∀1≤i≤n}⊂H1/2(Γ1),Mh={qh∈C0(Ω̅1)∩M;qh∣K∈P1,∀K∈τh}⊂H1(Ω1)∩M,
where Pl is the polynomial in two real variables of degree ≤l and si, 1≤i≤n (|si|≤h) are the straightline segments of the interfacial boundary Γ1.Recalling [6], there hold the following optimal error estimates:
(3.7)∥u-uh∥1,Ω1+∥λ-λh∥-1/2,Γ1+∥p-ph∥0,Ω1≤ch,
where c is a positive constant.


## 4. 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). Xh and Mh 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∈Xh its values at the vertices and midpoints of the triangulation τh and for Mh its values at the vertices of τh, with the constraint of the interelement continuity. On the other hand, Th=Sh2∩T is the space of continuous and piecewise linear functions μh on Γ1 with ∫Γ1μhds=0. We can construct its basis as follows: let us denote by Ai-1Ai¯=Γ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(Ai)=1, π1(Aj)=0(j≠i) and ∫Γ1πids=1; more precisely, with a parametric representation of each Γi (for x∈Γi, Ai-1x¯=τAi-1Ai¯, τ∈[0,1]), we define
(4.1)πi(x(τ))={τhonΓi,1-τhonΓi+1,0elsewhere.
Then Sh is generated by the set of n-1 functions {ti}, continuous on Γ1 and piecewise affine, satisfying
(4.2)ti=πi-πi+1,1≤i≤n-1.
Now, if N=Nv+Nm denotes the total number of nodal points (vertices Nv and midsides Nm) of τh, the solution of the approximate problem (Qh) reduces to the solution of a linear system of order 3Nv+2Nm+2(n-1)+1. It takes the following matrix form:
(4.3)Mx=b,
with
(4.4)x=(uh,ph,λh)′,b=(f,0,0)′,and where the arising matrixM may be written symbolically in a partioned form, exhibiting its natural block structure
(4.5)M=[ABCB′00C′0D]orM=[ACC′D],
with
(4.6)A=[ABB′0],C=[C0],D=D.HereA 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 submatrixD. They are of the following form:
(4.7)Dik=∬Γ1Uki(x-y)λhi(x)μhk(y)dsxdsy,i,k=1,2.
Taking into account the above definition of the finite element space Th and the construction of its basis functions (4.1)-(4.2), it is easy to see that, putting i,j=1,…,n-1(4.8)D=[[dij1][dij2][dij2][dij3]],
we must actually calculate
(4.9)d̅ij1=∬Γ1U11(x-y)πi(x)πj(y)dsxdsy,d̅ij2=∬Γ1U12(x-y)πi(x)πj(y)dsxdsy,d̅ij3=∬Γ1U22(x-y)πi(x)πj(y)dsxdsy,and use the relations(4.10)dijk=d̅ijk-d̅i,j+1k-d̅i+1,jk+d̅i+1,j+1k,k=1,2,3.To write the above coefficientd̅ijk, we need to add the contribution from two adjoining elements Γ1i and Γ1i+1, Γ1j, and Γ1j+1. We define
(4.11)αij=∫Γ1i∫Γ1jK(x-y)πi(x)πj(y)dsxdsy,βij=∫Γ1i∫Γ1j+1K(x-y)πi(x)πj(y)dsxdsy,γij=∫Γ1i+1∫Γ1j+1K(x-y)πi(x)πj(y)dsxdsy,which become, in terms of parametersσ,τ∈[0,1] of the straightline segments,
(4.12)αij=∬01στK(x(σ)-y(τ))dσdτ,βij=∬01σ(1-τ)K(x(σ)-y(τ))dσdτ,γij=∬01(1-σ)(1-τ)K(x(σ)-y(τ))dσdτ,whereK(x)=U11(x),U12(x), and U22(x).Computing these integrals and using the relationship(4.13)d̅ij=αij+βij+βji+γij,it is a simple matter to derive the full matrixD (refer to [5] for the more details).Let us now examine the coupling terms(4.14)〈vh,λh〉=∫Γ1vhλhds,〈G(γ0uh),μh〉=∑i,j,k=12∬Γ1σij(Uk,Pk)(x-y)n̅j(x)uhi(x)μhk(y)dsxdsy,
where vh=(vh1,vh2), μh=(μh1,μh2). According to (4.1)-(4.2), we can write
(4.15)∫Γ1vhitjds=∫Γ1vhiπjds-∫Γ1vhiπj+1ds=∫Γ1jvhiπjds+∫Γ1j+1vhiπjds-∫Γ1j+1vhiπj+1ds-∫Γ1j+2vhiπj+1dsi=1,2,j=1,…,n-1.


## 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 (Qh) produces a global system of linear equations which can be represented in a condensed matrix notation as
(5.1)Mx=band, more explicitly, in terms of the block structure of matrixM(5.2)[ACC′D][XY]=[b0],
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
(5.3)X=A-1(b-CY),(5.4)(D-C′A-1C)Y=-C′A-1b.
Hence, we need to compute A-1b and A-1C by solving the matrix equations
(5.5)Ab̅=b,AC̅=C.
Then we solve the equation
(5.6)(D-C′C̅)Y=-C′b̅,
by a Gaussian elimination algorithm with partial pivoting and obtain from (5.3) that
(5.7)X=b̅-C̅Y.


## 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 radiusR. 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 (x1,x2)-cartesian plane normal to it), we write
(6.1)ψ=logr,withr=x12+x22;the velocity distribution is then(6.2)u=ro̅tψ=(x2r2,-x1r2)
and the boundary condition
(6.3)u=(x2R2,-x1R2)atr=R
for a cylinder of the radius R. And we take w0=(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.4)(∂u1∂x1,∂u2∂x1)+∇p=f,sinceΔu=0. Then, we may write
(6.5)p=-x2r2+{R2-r,forR<r≤R2,0,forr≥R2
corresponding to
(6.6)f={-(x1,x2)r,forR<r≤R2,0,forr≥R2.
The interfacial boundary introduced in the fluid region is a circle of radius R2. The bounded domain Ω1 is then a circular annulus limited by R and R2.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)λ=(x1x2R23+2νx2R23,x22R23-2νx1R23).Computations were carried out withR=1 and R2=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.Table 1
Representative parameters of four basic meshes.
MeshesNLNNENONVPNT3×63636841922024×848641443283425×105101002205005186×12612144312708730Here,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.In Tables2 and 3, we give a summary of computational results for ν=1 and ν=0.5, respectively.Table 2
Relative error of numerical solution on four basic meshes (ν=1).
Meshessize∥u-uh∥L2/∥u∥L2∥∇(u-uh)∥L2/∥∇u∥L2∥p-ph∥L2/∥p∥L2∥λ-λh∥L2/∥λ∥L23×6h0.05242730.06641720.14270250.30906924×83/4h0.03423850.05472780.11772960.22758215×103/5h0.02436890.04742190.09989180.18136256×121/2h0.01845960.04276310.08764850.1488314Table 3
Relative error of numerical solution on four basic meshes (ν=0.5).
Meshessize∥u-uh∥L2/∥u∥L2∥∇(u-uh)∥L2/∥∇u∥L2∥p-ph∥L2/∥p∥L2∥λ-λh∥L2/∥λ∥L23×6h0.05357820.06626080.14956500.31026004×83/4h0.03589740.05802170.12778900.23047955×103/5h0.02625480.05160520.11217450.19047526×121/2h0.02413130.04574480.08764850.1675414From Tables2 and 3, one may observe that for the finer meshes we can obtain the better accuracy for the approximate velocity uh, pressure ph, 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.



---

*Source: 845984-2010-02-01.xml*