The DRM subdomain decomposition approach to solve the two-dimensional Navier–Stokes system of equations

EABE 1173

Engineering Analysis with Boundary Elements 24 (2000) 107–119 www.elsevier.com/locate/enganabound

The DRM subdomain decomposition approach to solve the two-dimensional Navier–Stokes system of equations H. Power 1,*, R. Mingo Wessex Institute of Technology, Ashurst Lodge, Ashurst, Southampton SO40 7AA, UK

Abstract This work presents a domain decomposition boundary integral equation method for the solution of the two-dimensional Navier–Stokes system of equations. It is known that, within this topic, the standard boundary element method is at a disadvantage when compared with classical domain schemes, such as finite difference and finite elements methods. In the present approach the original domain is divided into several subdomains. In each of them the integral representation formula of a non-homogeneous Stokes flow field is applied, and on the interface of the adjacent sub-regions the full matching conditions of the problem are imposed. The domain integrals resulting from the nonhomogeneous terms of the formulation are transformed into surface integrals at the contour of each sub-region via the dual reciprocity method. Finally, some examples showing the accuracy, efficiency and flexibility of the proposed method are presented. q 2000 Elsevier Science Ltd. All rights reserved. Keywords: Domain decomposition boundary element method; Dual reciprocity method; Nonlinear viscous flow

1. Introduction The boundary element method (BEM) is now a wellestablished numerical technique for the analysis of engineering problems. The basis of the method is to transform the original partial differential equation (or system of PDEs) into an equivalent integral equation (or system of equations) by means of the corresponding Green’s theorem, and the use of a fundamental solution. In this way, the resulting integral equation is numerically solved instead of the original PDE. The use of this approach for some linear problems directly leads only to surface integrals, and then only surface elements are necessary. It is known that the BEM is at a disadvantage when dealing with nonlinear problems in comparison with classical domain schemes such as finite difference and element methods. In the boundary element formulation of a nonlinear problem, it is usual to use an integral representation formula based upon the fundamental solution of a linear partial differential equation, and express the nonlinear terms as domain integrals. In early boundary elements analyses the evaluation of domain integrals was done using cell integration, a technique which, whilst effective and general, meant * Corresponding author. Tel.: 144-1703-293223; fax: 144-1703-292853. E-mail address: [email protected] (H. Power) 1 On leave from the Instituto de Meca´nica de Fluı´dos, Universidad Central de Venezuela, Venezuela.

that the method lost its boundary only nature, introducing an additional internal discretisation. Although good results can be obtained using the cell integration technique, this approach is in general several orders of magnitude more time consuming than domain methods, making this approach impractical for large problems. Several methods have been developed to take domain integrals to the boundary in order to eliminate the need for internal cells (boundary only BEM formulations). One of the most popular to date being the dual reciprocity method (DRM), introduced by Nardini and Brebbia [1]. This method is closely related to the method of particular integrals technique (PIT) introduced by Ahmad and Banerjee [2], also used to transform domain integrals to boundary integrals. In the latter method a particular solution satisfying the non-homogeneous PDE is first found and then the remainder of the solution, satisfying the corresponding homogeneous PDE, is obtained by solving the corresponding integral equations. The boundary conditions for the homogeneous PDE must be adjusted to ensure that the total solution satisfies the boundary conditions of the original problem. The DRM also uses the concept of particular solutions, but instead of obtaining the particular solution and the homogeneous solution separately, it applies the divergence theorem to the domain integral terms and converts the domain integrals into equivalent boundary integrals. Several successful applications of the BEM to the complete Navier–Stokes system of equations based upon cell integration have been previously reported (for a good

0955-7997/00/$ - see front matter q 2000 Elsevier Science Ltd. All rights reserved. PII: S0955-799 7(99)00043-0

108

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

literature survey see Ref. [3] and Power and Partridge [4]). More recently, Zheng et al. [3], Power and Partridge [4] and Rahaim and Kassab [5] reported boundary only integral equation approaches to solve two- and threedimensional nonlinear viscous flow problems. Zheng et al. used the method of particular integrals, while Power and Partridge, and Rahaim and Kassab used the dual reciprocity approach. Zheng et al. [3] solved two-dimensional problems using the penalty function method. In this type of formulation, a penalty function is used to relate the pressure to the divergence of the velocity field, in order to reduce the Navier– Stokes equations to the Navier equations of elastostatics, with a body force given by the original convective terms. Two challenging computational problems were solved by them, namely flow in a divergent channel and plane free jet, for which complex flow patterns are present, i.e. reversed flow configuration and the presence of a free surface, respectively. The scheme used by Power and Partridge [4] was based on the Stokes’ formulation, in which, as before, the original convective terms were considered as non-homogeneous forcing terms. Their numerical scheme was applied to the case of three-dimensional steady flow inside a closed circular cylinder container, where the bottom and side walls are rotating with constant angular velocity and the top is fixed. In this case the flow configuration is also complex, due to the presence of counter rotating cells. Zheng et al. [3], and Power and Partridge’s [4] numerical solutions were limited to small Reynolds number, namely 40 and 100, respectively. In their examples, convergence becomes slower as the Reynolds number increases and sometimes the numerical solution becomes unstable. On the other hand, Rahaim and Kassab reported two-dimensional, nonisothermal solutions up to Reynolds’ number of 2300, by solving a couple system of four non-homogeneous Laplaces equations, two for the velocity components, one for the pressure and one for the temperature field. However, in contrast to the previous two works, their problem consists of slowly varying flow in a uniform channel, i.e. small convection effects, where no complex flow patterns are presented. A major problem encountered with PIT and DRM is that the resulting algebraic system consists of a series of matrix multiplications of fully populated matrices. When only a few internal points are required in the PIT or DRM approaches, the resulting computing time is in general smaller than the one corresponding to a cell integration scheme, still being costly in comparison with domain approaches. In addition, in complex problems these two approaches have been only limited to a small order of nonlinearity. When dealing with the BEM solution of large problems, it is usual to use the method of domain decomposition, in which the original domain is divided into sub-regions, and in each of them the full integral representation formula is applied. At the common interface between the adjacent sub-

regions, the corresponding full matching conditions are enforced. While the matrices, which arise in the single domain formulation, are fully populated, the sub-region formulation leads to block banded matrix systems with one block for each sub-region and overlaps between blocks when sub-regions have a common interface. When using continuous elements of high orders (more than constant), the enforcement of the matching conditions at common interfaces, or in other words the matrix assembly, leads to an over-determined system of algebraic equations (as more sub-regions are defined, the bigger is the over-determination). Several schemes are known to reduce the obtained over-determined system (for more details see Ref. [6]). The simplest approach to overcome this difficulty is to express the derivatives of the field variables at the common nodes between more than two sub-regions in terms of the variables themselves, by using the interpolation functions or a finite-difference approximation. A way of avoiding this problem is by using discontinuous elements in such common nodes, obtaining in this way a closed matrix system at the expense of having more unknown variables inside the problem’s domain. Recently Popov and Power [7] found that the DRM approximation of the internal potential of a convection diffusion problem with variable velocity could be substantially improved by using the domain decomposition scheme. Their idea of using this type of scheme to improve the accuracy of the DRM approach was suggested by Kansa and Carlson’s [8] work on radial basis function data approximations, where they observed that the best approximation is obtained when the original domain is split into matching sub-domains. In this work we proceed with Popov and Power’s idea of improving the DRM approach by using the domain decomposition scheme. Here, we implement the domain decomposition scheme to the two-dimensional version of Power and Partridge’s [4] DRM formulation for the Navier–Stokes system of equations. To avoid the problem of the over-determined system, discontinuous elements are employed at the common nodes between sub-regions, allowing us to obtain the velocity field and the traction at every point over the sub-regions’ contours and the velocity field at the internal nodes of the subregions. Also, our DRM formulation is based upon the Duchon’s thin plates splines (TPS) interpolant functions, instead of the classical 1 1 R function used by Power and Partridge. We chose these interpolant functions because they have been considered in the literature of mathematical interpolation as the optimal multivariate interpolant functions. An example is presented where it is possible to appreciate the main advantage obtained when the TPS functions are used instead of the 1 1 R function.

2. DRM sub-domain decomposition formulation Let us consider the steady state Navier–Stokes system of

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

equations for two-dimensional incompressible flow in a bounded domain V . 2ui ˆ0 2xi

m

…1†

22 ui 2p…x† 2u 2 ˆ ruj i ; 2xj 2xj 2xj 2xi

i ˆ 1; 2; x[V

…2†

with the non-slip boundary conditions ui …x† ˆ Ui …x†;

…3†

x[S

It is known that the computation of the numerical solutions of the above system of equations experiences increasing difficulty as the flow Reynolds number increases, particularly for configurations encountering flow separation. In the sub-region BEM approach, the original domain is divided into sub-domains and in each sub-domain the solution is given in terms of the corresponding integral representation formula, which has to be matched with the adjacent sub-domain according to the corresponding compatibility conditions, i.e. flow velocity and surface traction. The Navier–Stokes system of equations can be interpreted as a non-homogeneous Stokes system of equations 2ui ˆ0 2xi

m

…4†

22 ui 2p…x† 2 ˆ gi …x†; 2xi 2xj 2xj

i ˆ 1; 2; x[V

…5†

with gi …x† ˆ ruj

2ui 2xj

…6†

as the non-homogeneous term. From the integral representation formula for the Stokes system of equations (see Ref. [9]), it is found that the i component of the velocity field at a point x of the nth subregion bounded by the contour Sn that enclose the domain V n, is given by Z Z Kij …x; y†uj …y† dSy 2 Jik …x; y†tk …y† dSy cki uk …x† ˆ Sn

1

direction, Jik …x; y†

  1 …xi 2 yi †…xk 2 yk † ˆ2 ln…r†dik 1 4pm r2

Vn

Jik …x; y†gk …y† dy

…7†

for n ˆ 1; 2; …; M; where M is the total number of subregions. In the above equations we have written x instead of x and y instead of y, and for convenience this notation is used from here on, and ti(y) is the i component of the surface traction due to the flow field …~u; p†; i.e. ti …y† ˆ sij ‰u…y†; p…y†Šnj …y†: The kernel Jik …x; y† is the i component of the velocity field of the fundamental solution of Stokes system of equations (Stokeslet), located at the point y and oriented in the kth

…8†

where corresponding pressure q k is given by   1 …xk 2 yk † ; qk …x; y† ˆ 2 2p r2

…9†

where r is the distance from the point y, to any point x, i.e. r ˆ ux 2 yu; and Kij …x; y† ˆ sij ‰J~k …x; y†; qk …x; y†Šnk …y† 1 …xk 2 yk †…xi 2 yi †…xj 2 yj † nk …y† …10† pr r3 is the stress tensor corresponding to the flow J~ k ; qk: The coefficients cki have values between d ki and 0, being equal to (1/2)d ki for smooth boundaries. It is also important to point out that the above integral representation formula holds for points inside the domain V n, in which case cki ˆ dki: The final set of equations is completed by assembling the above integral equations for each sub-region, utilising the traction equilibrium and the velocity compatibility between ˆ u…2† and common interfaces of the sub-regions, i.e. u…1† i i …1† …2† sij nj ˆ 2sij nj : Until now, the above formulation is valid for both the boundary element sub-domain (BEM-SD) method and the present dual reciprocity sub-domain approach (DRM-SD). In the BEM-SD approach, standard cell integration is used to cope with the domain integrals in the integral representation formula at each sub-domain. In the present DRM-SD approach, instead of using cell integration we use the DRM approach in terms of the most efficient radial functions interpolants recently reported in the mathematical literature. Discontinuous nodes are used at the elements connecting more than two sub-regions in order to avoid the problem of an over-determined system. To express the domain integral in Eq. (7) in terms of equivalent boundary integrals, a DRM approximation is introduced. The basic idea is to expand the g~ …y† term using radial approximation functions, at each sub-region n, i.e. ˆ2

Sn

Z

109

gi …y† ˆ ruj

p X 2ui ˆ f …y; zm †am l dil ; 2yj mˆ1

l ˆ 1; 2

…11†

In this way the domain integral at each sub-domain becomes Z Vm

Jik …x; y†gk …y† dV ˆ

p X mˆ1

aml

Z Vm

Jik …x; y†f …y; zm †dkl dV …12†

k

The functions f(y,z ) are approximating functions, which depend only on the geometry of the problem, and am l are unknown coefficients to be determined. The approximation

110

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

is done at …p ˆ p1 1 p2 † nodes, with p1 boundary nodes around the contour of the sub-region and p2 nodes inside the sub-region. The most popular choice for the functions f(y,z k) in the DRM approach is f …y; zk † ˆ 1 1 R…y; zk †

…13†

where R(y,z k) is the distance between pre-specified fixed collocation point z k and a field point y, where the function is approximated, i.e. R…y; zk † ˆ uy 2 zk u: In the present case, both points are located in the sub-region of consideration. The function (13) is a member of a family of functions called radial basis functions (RBFs), which are related to the theory of mathematical interpolation. Micchelli [10] proved that for the case when the nodal points are all distinct, the matrix resulting from a radial basis function is always nonsingular. Therefore, as long as the function g~…y† is regular, the coefficients, a~ m , are well defined. In recent years, the theory of RBFs has undergone intensive research and enjoyed considerable success as a technique for interpolating multivariable data and functions. Unlike other interpolating functions, RBFs are not restricted to problems with only convex hulls or uniform data spacing. In 1982, Franke [11] published a review article evaluating virtually all the interpolation methods for scattered data sets available at that time. Among the methods tested, RBFs outperformed all the other methods regarding accuracy, stability, efficiency, memory requirement and ease implementation. In a similar study, Stead [12] examined the accuracy of partial derivative approximations over scattered date sets, and also concluded that the RBFs performed more accurately compared to other considered methods. Of the RBFs tested by Franke, Hardy’s multiquadrics (MQ) were ranked the best in accuracy, followed by Duchon’s TPS. Duchon [13] derived the TPS as an optimum solution of the interpolation problem in a certain Hilbert space via construction of a reproducing kernel. Therefore, they are the natural generalisation of cubic splines in n . 1 dimension. Even though TPS have been considered as optimal for interpolating multivariate functions, they converge only linearly (see Ref. [14]). On the contrary, the MQ functions converge exponentially and they always produce a minimal semi-norm error as proved by Madych and Nelson [15]. However, despite MQ’s excellent performance, it contains a free parameter, often referred to as the shape parameter, whose choice can greatly affect the accuracy of the approximation. The choice of the optimal shape parameter is a problem that has received the attention of many researchers. So far, this is an open problem and no mathematical theory has been developed yet for determining the optimal value. In addition, recently Kansa and Carlson [8] reported that for RBFs’ interpolation, the best result is obtained when the original domain is split into matching sub-domains. In this work we develop the DRM-SD approach for the Navier–Stokes system of equation in terms of the Stokes

fundamental solution and the use of the TPS interpolant functions. Our criterion to choose this approach was based upon the mentioned properties of the RBFs’ interpolation functions, in particular their improvement obtained by the use of the domain decomposition approach [8], and the efficiency and accuracy obtained with the DRM-SD approach to solve the convective diffusion equation, reported in the introduction of this article. Although it is possible to use different RBF interpolations in different sub-regions, according to the expected behaviour of the unknown solution, here for simplicity we only use one type of interpolation for all the sub-regions. To proceed with the DRM idea of reducing the last domain integral in Eq. (12) to equivalent boundary integrals, let us define a new auxiliary non-homogeneous Stokes field, J^ li …x; zm †; for each collocation point z k, in the following way

m

22 J^li …x; zm † 2p^ l …x; zm † 2 ˆ f …x; zm †dil 2xj 2xj 2xi

2J^ li …x; zm † ˆ0 2xi

…14†

…15†

Applying the integral representation formula to the above auxiliary non-homogeneous Stokes field, at a point x, we obtain Z Kij …x; y†J^lj …y; zm † dSy cki J^lk …x; zm † ˆ Sn

2 1

Z Sn

Z Vn

Jik …x; y†t^lk …y; zm † dSy Jik …x; y†f …y; zm †dkl dy

…16†

Substituting the last equation into Eq. (12), the domain integral can be recast in terms of a series of surface integrals, and using the resulting expression in Eq. (7), one finally arrives at a boundary only integral representation formula for the velocity field at each sub-domain Z Z Kij …x; y†uj …y† dSy 2 Jik …x; y†tk …y† dSy cki uk …x† ˆ Sn

1

Sn

p X mˆ1

1

Z Sn

aml {cki J^ lk …x; zm † 2

Z Sn

Kij …x; y†J^ lj …y; zm † dSy

Jik …x; y†t^lk …y; zm † dSy } …17†

The above formula …i ˆ 1; 2† has only boundary integrals and in their evaluation we need the additions of a series of terms, each of them representing the effect of the particular solution localised at the m pole. For the numerical solution of the problem, the surface Sn of each sub-region is discretised by means of straight-line elements. Within each element the density of the integrals in

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

the above equation is defined in terms of their nodal values and a quadratic interpolation. In this way the dimensionless form of Eq. (17) can be written in matrix form, as: …cu†i 2 Hik uk 1 Gik tk ˆ Re

p X

^ ij †Š ‰aj …Gik t^kj 2 Hik j^kj 1 …cj†

mˆ1

…18† where Hik and Gik are the standard influence matrices resulting from the integration over the boundary elements, the index i represents the collocation nodes, k the nodes at the integration elements and j the DRM collocation points in matrix notation. At this stage it is possible to invert Eq. (11) and express the value of the coefficients a^ m in terms of the nodal values of the function g~; i.e.

aml ˆ ‰Fil …ym ; yq †Š21 gi …yq †

…19†

An algorithm must be established to relate the values of g~ to the nodal values of the velocity vector u~: In the DRM approach, this is usually done by expressing the velocity field at a domain point in term of the approximation function F, i.e. ui …x† ˆ Fil …x; yn †Bnl

…20†

~m

the coefficients B are given by the unique solution of the algebraic system of equations obtained when the above equation is evaluated at each of the nodal points, x ˆ yq ; q ˆ 1; 2; …p; i.e. Bnl ˆ ‰Fil …yn ; yq †Š21 ui …yq †

…21†

It is important to point out that although the existence of ~ m is guaranteed by the radial function the coefficients B theory, their value cannot be found until the complete problem has been solved, since the velocity vector u~ , at any point in the domain V , is one of the unknowns within the problem. In this way, the derivative of the velocity field can be given by " # 2ui …x† 2 n ˆ …F …x; y †† Bnl ; …22† ui;j …x† ˆ 2xj 2xj il and hence the convective term can be written as a function of the velocity field ui as: # ) (" 2 n …F …x; y †† Bnl uj …x† …23† {ui;j }uj …x† ˆ 2xj il Substituting Eq. (21) into the above equation we have: (" # ) 2 n n q 21 q …F …x; y †† ‰Fkl …y ; y Š uk …y † uj …x† {ui;j }uj …x† ˆ 2xj il …24† At this stage it is important to point out the nonlinearity of the above relation, expressed by ui(y q)uj(x). Finally by substitution of the above equation into Eq. (19)

111

it is possible to obtain the coefficients a~ q in terms of the velocity field at the interpolation nodal points as:

aml ˆ ‰‰Fil …ym ; yt Š21 Š # (" ) 2 t n n q 21 q …F …y ; y †† ‰Fks …y ; y Š uk …y † uj …yt † …25†  2xj is To deal with the nonlinearity of the problem, a Newton– Raphson scheme is implemented and Eq. (18) can be rewritten as: ^ 2 ‰GŠ‰t^Š†a ‰HŠu 2 ‰GŠt ˆ Re…‰HŠ‰jŠ

…26†

where ‰HŠ ˆ …Hik 2 cdik †;

‰GŠ ˆ Gik ;

^ ˆ j^kj ; ‰jŠ

‰t^Š ˆ t^kj

…27†

The function to be optimised in the Newton–Raphson scheme is: ‰HŠu 2 ‰GŠt 2 Re‰CŠa ˆ F…u; q† ^ 2 ‰GŠ‰t^Š† where ‰CŠ ˆ …‰HŠ‰jŠ then  " Du # 2F ˆ ‰2FŠ 2X Dt

…28†

…29†

The coefficients of the Jacobian matrix, Jij ˆ 2Fi =2Xj ; are computed by differentiating the function F(u,t). The terms corresponding to Xi ˆ ui are: Jij ˆ Hij 2 ReDij

…30†

where Dij ˆ Cil

2 al 2uj

…31†

and the one corresponding to Xj ˆ tj are Jij ˆ 2Gij

…32†

It can be observed that if we express the system to obtain a first-order Stokes’ approximation as: Ax ˆ b;

…33†

the coefficients of the Jacobian matrix Jij, given before, and the coefficients of the above system, Aij, are almost the same except that the terms corresponding to Xj ˆ uj in Aij, which are only given by Hij instead of Eq. (30). Therefore, to compute the matrix Aij and the Jacobian matrix, it is only necessary to consider the matrix Aij and then subtract the matrix ReDij to obtain the Jacobian matrix in the iteration process. 2.1. Particular solution of the auxiliary system The DRM is essentially a generalised way of constructing

112

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

particular solutions. We first need to define the function f used in the DRM interpolation, in order to find the particular solution of the non-homogeneous Stokes system of Eqs. (14) and (15). Although a variety of functions can, in principle, be used as a basic interpolant function, it can be concluded, from the previous discussion about the radial basic function interpolation, that the best approximation is obtained with the TPS and the MQ functions. Here we use the TPS interpolant due to the mentioned difficulty with the undetermined constant in the MQ functions. In a typical problem we have p pairs of data points {…zm ; g…zm ††pmˆ1 }; which are assumed to be samples of some unknown function g(x). The Duchon TPS functions f …x† ˆ

m2

am ux 2 z u logux 2 z u 1 a 1 bx1 1 cx2 m

32 2 2 X 2 f …x† 4 5 dx1 dx2 R2 …i;j†ˆ1 2xi xj

…35†

p X

am zm1 ˆ

mˆ1

p X

am zm2 ˆ 0

…36†

mˆ1

Here, we use the above interpolant for each of the components of the vector field g~ …x†. In order to find the corresponding particular solution of the non-homogeneous Stokes system of Eqs. (14) and (15), with the TPS function as forcing term, we use the approach suggested by Power and Partridge [4] to find the corresponding three-dimensional particular solution with 1 1 R as a forcing term. Let us define the second rank tensor J^ li …x; zm † in terms of an auxiliary potential c , as follows: J^li …x; zm † ˆ

24 c 2p^l 1 ˆ0 2xk 2xk 2xi 2xl 2xi

auxu4 64 …40†

…41†

2c 2c d 2 ; 2xk 2xk il 2xi 2xl 2

23 c 2xk 2xk 2xl

…42†

Finally substituting the biharmonic potential c into Eqs. (37) and (42) and making the respective derivatives, we arrive at the following expression for the non-homogeneous auxiliary Stokes flow field:  X  p am 7 5R4 log R 2 R4 dij 3 96 mˆ1   5 a ‰3R20 dij 2 2xi xj Š 2 x^i x^j 4R2 log R 2 R2 1 3 16

b 3 ‰x …3d 2 2d1i d1j 2 d2i d2j † 1 3x22 x1 …dij 2 d1i d1j † 24 1 ij c 3 ‰x …3d 2 2d2i d2j 2 d1i d1j † 2 3x21 x2 …d1i d2j 1 d2i d1j † 1 24 2 ij  2 2 …43† 1 3x1 x2 …dij 2 d2i d2j † 2 3x2 x1 …d2i d1j 1 d1i d2j †Š 1

where x^l ˆ …xl 2 zm l †

and

R0 ˆ uxu

2

R ˆ ux 2 zm u

…37†

By differentiating the above equations, we find that the continuity Eq. (15) is satisfied by any choice of the potential c . Substituting expression (37) into Eq. (14) we obtain 24 c 24 c m dil 2 2xk 2xk 2xk 2xk 2xk 2xk 2xi 2xl

! 2

2p^l ˆ f …x†dil 2xi …38†

If we assume that c satisfies the non-homogeneous biharmonic equation:

m

1

1 J^ ji ˆ m

when the coefficients (a,b,c) satisfy the constraints

mˆ1

!

then the pressure field has to satisfy the following equation:

p^ l …x; zm † ˆ 2

2

am ˆ

ux 2 zm u6 logux 2 zm u 5ux 2 zm u6 2 cˆ am 576 3456 mˆ1 ! ! b x51 x31 x22 c x52 x21 x32 1 1 1 ; 1 24 24 2 120 2 120

…34†

interpolate the values of g(x) at the p points and minimises (in R 2) the functional

p X

p X

or

p X mˆ1

Z

whose solution is

24 c ˆ f …x† 2xk 2xk 2xk 2xk

…39†

Fig. 1. Simulated velocity field at Re ˆ 100:

…44†

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

and t^ki ˆ skij nj ˆ



 8R2 …x^i nk 1 x^j nj dik 1 x^k ni † × 2 log R

n X a 96 mˆ1    1 1 a 2 4x^i x^k x^j nj 4 log R 1 1 ‰xi nk 2 3 3 4 m

1 xj nj dik 1 xk ni Š 1

b 2 {x ‰3…n1 dik 1 nk d1i 8 1

1 ni d1k † 2 2…2n1 d1i d1k 1 n1 d2i d2k 1 n2 d1i d2k 1 n2 d1k d2i †Š 1 x22 ‰n1 dik 1 nk d1i 1 ni d1k 2 2n1 d1i d1k Š 1 2x1 x2 ‰…n2 dik 1 nk d2i 1 ni d2k † 2 2…n2 d1i d1k 1 n1 d1i d2k 1 n1 d1k d2i †Š 1

c 2 {x ‰3…n2 dik 1 nk d2i 1 ni d2k † 2 2…2n2 d2i d2k 8 2

1 n2 d1i d1k 1 n1 d2i d1k 1 n1 d2i d1k 1 n1 d1i d2k †Š 1

x21 ‰n2 dik

1 nk d2i 1 ni d2k 2 2n2 d2i d2k Š

1 2x1 x2 ‰…n1 dik 1 nk d1i 1 ni d1k † 2 2…n2 d2i d1k 1 n1 d2i d2k 1 n2 d2k d1k †Š}

…45†

Finally, it is important to observe that all the matrices, i.e. ^ t^; F, F 21, and 2F=2xl ; are geometrically defined in H, G, J; terms of the distance between the points defining the subregions, except for the last part of the TPS interpolation, which are functions of the local coordinate of the nodal points when those functions are defined in terms of a global coordinate system. If the interpolant functions are given in terms of a local coordinate system at each sub-region, it follows that all the matrices are geometrically defined, even the extended polynomials of the TPS, and therefore those matrices are identical for all the sub-regions that have the same size and shape, for which it is important to define the same local numeration of the nodes, and consequently in those sub-regions the above matrices can be evaluated in terms of only one of the sub-regions.

113

important consequence of this dependency is that the accuracy of the numerical solution, at different values of the Reynolds number, is strongly associated with the way in which the mesh discretisation is defined, requiring a higher density in the regions near the contours, where the boundary-layer is developed, as the Reynolds number increases. For high Reynolds numbers, the use of a coarse mesh often yields a convergent solution without physical meaning. By refining the computational mesh, it is possible to improve the solution of the flow phenomena, but at the expense of increasing the computational cost. In the present examples, particular attention has been given to the ability of the proposed formulation to predict complex flow behaviours at moderate Reynolds number, using a coarse uniform sub-region mesh, for which the computer code can be defined in a block structure in terms of a single sub-region, evaluating, in this way, all the other matrices in terms of this sub-region, and assembling them according to the domain decomposition used. 3.1. Laminar flow in a lid-driven square cavity Let us consider a two-dimensional square cavity whose upper boundary is moving at a constant velocity. A zero velocity condition is applied on all the walls, except at the top one, which moves with a non-dimensional horizontal velocity equal to unity. Note that the solution for Re ˆ 0 can be obtained without iteration in terms of a boundary only integral formulation. In all the cases of Re . 0; we use the solution of Re ˆ 0 as the initial guess in the iterative scheme. In Figs. 1–3, we plot the velocity field, and the vertical and horizontal centre-line velocity profiles, respectively, obtained with the present approach without the use of the domain decomposition scheme for the case of Re ˆ 100: In this test, we used 48 quadratic elements and 169 internal

3. Numerical examples To illustrate the numerical scheme described above, two benchmark problems have been chosen to be solved using the proposed DRM-SD approach at different values of the Reynolds number; namely, flow in a lid-driven square cavity and flow over a backward facing step. In the Navier–Stokes system of equations, the intensity of the nonlinearity, i.e. the convective effect, is known to be related with the magnitude of the Reynolds number. Therefore, the computational effort required for its numerical solution increases as the Reynolds number increases. An

Fig. 2. Horizontal fluid velocity distribution at the vertical centre-line of the cavity at Re ˆ 100:

114

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

this anomalous flow behaviour are given later). As can be observed, for this small value of the Reynolds’ number, our results based upon a DRM formulation are in close agreement with the benchmark values obtained by Ghia et al. [16], using a finite difference multigrid numerical scheme, even when we use only small number of internal points. For Reynolds’ numbers larger than about 100, it is not possible to achieve an accurate solution with the use of the present DRM formulation without the use of the domain decomposition approach, even if we increase the number of surface elements and internal points. For values of the Reynolds’ number higher than 200, convergence problems are encountered. As pointed out previously, at the flow regions near the corners, the flow contains singularities associated with infinite values of the velocity and surface traction. There are two different types of such anomalous flow behaviours; one associated with the flow within a corner that is confined between two intersecting stationary walls, driven by the motion of the fluid away from the walls [17], and the other due to the motion of only one of the two walls [18]. In the first case, a singular solution is encountered due to a disturbance at the corner, with the velocity diminishing away from the corner, while in the second case an infinite force is required to maintain the scraper motion of the boundary over the fixed plane wall.

Fig. 3. Vertical fluid velocity distribution at the horizontal centre-line of the cavity at Re ˆ 100:

points. Semi-discontinuous quadratic elements were used to avoid the problems at the corners of the cavity where, besides having two normal vectors at a single point, a singular behaviour of the solution is expected (more details about

1

0.5

0

-0.5

Tayl or analitical solution -1

Numerical solution fine mesh Numerical solution coarse mesh

-1.5

-6

-4

-2

0

2

4

Fig. 4. The x component of the surface traction along a corner.

6

8

10

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

115

Fig. 5. Velocity field at Re ˆ 250 obtained with the use of the function f ˆ 1 1 r; for SR ˆ 9 (where SR is the number of sub-regions).

Higdon [19] and Pozrikidis [20] found that the presence of sharp corners, i.e. Moffatt type of singularity, neither decreases the global accuracy nor discredits the overall liability of the boundary element numerical solution. In particular, they observed that there is not a problem when quadratic elements are used but it can lead to an error, which is of O(1) in the immediate vicinity of the corner, if constant elements are employed. On the other hand, Kelmanson [21,22] found that a discontinuous boundary velocity, i.e. Taylor type of singularity, might cause oscillatory solutions and slow convergence. One way to circumvent these problems is to use a high density of boundary elements in the vicinity of the singular points. Practice has shown that this is an effective remedy but not a panacea [23]. A more formal approach is to remove the singular behaviour by representing the solution in terms of a regular solution, plus the known local singular component, and to solve the integral equation for the regular component [24]. As previously commented, a simple way of avoiding this problem is by using a discontinuous node in the neighbourhood of such corners. Fig. 4 shows the obtained x component of the surface traction along the upper right corner of a long rectangular channel of depth–length ratio 1 : 3; with

Fig. 6. Velocity field at Re ˆ 250 obtained with the use of the function f ˆ r × r log…r† 1 a 1 bx 1 cy; for SR ˆ 9 (where SR is the number of subregions).

Fig. 7. Velocity field at Re ˆ 400 obtained with the use of 16 sub-regions with 12 elements at each side of each sub-region.

zero inlet and bottom velocity, constant horizontal velocity and zero vertical velocity at the top surface and a plane Couette flow at the outlet section, at zero Reynolds’ number, in such a way that a Taylor type singularity is expected at the left hand side top corner. This result was obtained by using discontinuous elements at the corners of the channel. In the figure we compare Taylor’s analytical solution with two of our numerical results, one with double the element density of the other around the upper right corner. As can be observed, as the density of the elements increases around the corner, the numerical solution tends to agree with the Taylor analytical solution. The difference between the analytical and the numerical solutions at the left hand side of the figure, which corresponds to the vertical wall at the beginning

Fig. 8. BEM internal cell distribution used by Kitagawa et al. [25] in order to obtain an accurate solution at Re ˆ 400:

116

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

Fig. 9. Horizontal fluid velocity distribution at the vertical centre line, at Re ˆ 400; for different discretisations.

of the channel, is due to the finite channel depth in our example, which is not the case in the Taylor solution. With the idea of analysing the behaviour of our DRM-SD formulation when different radial interpolation functions are employed, we carried out two numerical tests in the cavity problem, at Re ˆ 250; using a coarse mesh of nine subregions with four quadratic elements on each face and nine internal points at each sub-region. One of the tests was with the use of the function f ˆ 1 1 r and the other with augmented TPS function. In Figs. 5 and 6, we can appreciate the difference obtained with each of the interpolating functions. While the results obtained with the

Fig. 10. Vertical fluid velocity distribution at the horizontal centre line, at Re ˆ 400; obtained with the 16 sub-region mesh.

function f ˆ 1 1 r does not have any physical meaning, the solution with the TPS function produces the expected flow configuration. In order to analyse the capability of the proposed domain decomposition for predicting the flow behaviour at moderate Reynolds number, using a coarse uniform sub-region mesh, the case of Re ˆ 400 was carried out. Fig. 7 shows the obtained velocity flow field for the case of Re ˆ 400; where the original square domain was subdivided into 16 equal square sub-regions, defined each of them by 56 surface elements and 100 internal points. As can be observed, the mesh used with the present approach is substantially less refined than those used in previous cell integration approaches. As an example, in Fig. 8 we reproduced the mesh used by Kitagawa et al. [25] in the penalty function method where, 160 linear boundary elements and 157 constant internal cells were required to obtain an accurate solution at Re ˆ 400: It should be noted that in the above cell mesh, a higher level of refinement was used near the corners but this non-uniformity was not necessary in our case. Figs. 9 and 10 show a comparison between the obtained horizontal and vertical centre-line fluid velocity and those reported by Ghia et al. [16] based upon the finite-difference multi-grip solution, for the case of Re ˆ 400: In Fig. 9, we also present a convergence analysis of our scheme, when we increase the number of sub-regions and the number of internal points at each sub-region. In the figure we give the results obtained with nine equal square sub-regions with 32 surface elements for each sub-region with the use of the function f ˆ 1 1 r; and 48 surface elements for each sub-region with the use of the augmented TPS function, and the one obtained with 16 equal square sub-regions

Fig. 11. Velocity field at Re ˆ 600 obtained with a non-uniform mesh distribution.

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

117

Table 1 Comparison of reattachment lengths obtained with various methods at Re ˆ 73

Fig. 12. Vectors field at Re ˆ 50 for the channel problem with a backward facing step.

with 56 and 64 surface elements for sub-regions and the use of the augmented TPS function. As in the previous case of the DRM approach without domain decomposition, in the present DRM-SD approach, we reach a value of the Reynolds number, Re . 400; for which it was not possible to achieve an accurate solution with a uniform coarse mesh. In Fig. 11, we show the obtained velocity field with 25 non-uniform square subregions for the case of Re ˆ 600; which is in good agreement with previous cell integration formulations obtained with more refined meshes. Due to the formation of significant recirculation zones at the two bottom corners, when the Reynolds number is larger than 600, a refinement of the DRM-SD’s mesh at these corners is required in order to obtain an accurate solution beyond this value of the Reynolds number. Notice that with non-uniform meshes, we lost the advantage of calculating at once all the matrices coming from the DRM-SRD scheme. 3.2. Laminar flow over a backward facing step In this example we consider the flow in a two-dimensional channel with a downstream facing step. Due to the abrupt change in the channel width, a recirculation zone occurs even at very low Reynolds number. This flow configuration has been confirmed by experimental data and other numerical solutions (see Ref. [26] and more recently Zienkiewicz et al. [27]). The most common criterion to judge the performance of a numerical scheme for this problem is the prediction of the

Backward facing step

Reattachment lengths

Denham and Patrick [26] Taylor et al. [28] Skerget et al. [29] S–G method [27] Zienkiewicz et al. [27] Present

3.9 (experimental) 5.3 3.95 4.9 4.8 3.85

reattachment length behind the step. In this section, we compare our numerical result for the reattachment length with Denham and Patrick’s experimental results, Taylor et al.’s [28] finite element solution and Skerget et al.’s [29] boundary element solution in terms of the velocity vorticity formulation using standard cell integration, and with two recent finite element results, no upwind (S–G method) and upwind, reported by Zienkiewicz et al. [27]. The channel is defined by inlet and outlet boundaries and two non-slip walls between them. Parabolic velocity profiles are prescribed as the inlet and outlet boundary conditions, corresponding to a fully developed laminar flow in the channel, following the proposal of Skerget [29], which is the example selected in order to validate our numerical result. Therefore, special care has to be taken in the definition of the total channel length in order to guarantee the outlet conditions. As in the previous example, particular emphasis is paid to the capability of the formulation to predict complex flow behaviours with the use of a uniform coarse sub-region mesh. A channel of total length equal to eight times the inlet width with an expansion ratio of the step of 2:3 was considered. The entry is at a distance of two step heights upstream of the step, and the outlet boundary has been taken at 14 step heights downstream of the step. The computational

Fig. 13. The recirculation zone at Re ˆ 50

118

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

Fig. 14. The recirculation zone at Re ˆ 73:

domain was divided into 23 sub-regions with the following characteristics: each sub-region has 12 quadratic elements on each side, discontinuous elements at the corners of the sub-regions were used, and 50 internal points for the subregions were employed. With the above sub-region mesh we were able to obtain an accurate solution for Re ˆ 50 and Re ˆ 73: Fig. 12, for Re ˆ 50; shows the obtained velocity field. A detailed view of the recirculation zone for the example at Re ˆ 50 is presented in Fig. 13, where the obtained reattachment length can be observed. The obtained value of the reattachment length, Lr ˆ 2:7; is in good agreement with Denham and Patrick’s [26] experimental result, Le ˆ 3: Fig. 14 presents the recirculation zone and the reattachment length for the case of Re ˆ 73: As before the obtained value of the reattachment length, Lr 3:85; is in excellent agreement with Denham and Patrick’s experimental result, Le ˆ 3:9: In Table 1 it is possible to see the comparison of our results against another numerical method and the corresponding experimental value for this last case. A major advantage of a BEM formulation of the present problem, in relation with other numerical techniques, is obtained in the definition of the reattachment point, the presence of which is manifested in the change of sign of the obtained wall-shear-stress distribution. As in the case of the cavity, in order to obtain an accurate solution in the present problem for values of the Reynolds number larger than about 100, it is necessary to use a non-uniform mesh distribution.

4. Concluding remarks In this paper, we have described a numerical method, based upon the use of the Stokes’ integral representational formulae, where the original nonlinear terms in the

Navier–Stokes system of equations appear as domain integrals. The problem’s domain is then sub-divided in a multizone, and in each sub-region the integral representation is applied, adjacent sub-regions are matched according to the velocity and surface traction compatibility conditions. The domain integrals at each sub-region are treated by the DRM approximation in terms of the most efficient interpolant functions available in the mathematical literature. The proposed method is effective in solving complex flow problems, as are the ones previously presented, even with the use of a uniform coarse sub-region mesh, for which the resulting computer code can be defined in a block structure form in terms of a single sub-region. For more complicated problems or flow conditions, where an accurate solution cannot be achieved with uniform domain decomposition, the resulting non-uniform sub-region mesh is significantly less refined than those required by other BEM formulations of the problem.

Acknowledgements The authors thank the reviewers for their helpful suggestions.

References [1] Nardini D, Brebbia CA. A new approach to free vibration analysis using boundary elements. Appl Math Model 1983;7:157–62. [2] Ahmad S, Banerjee PK. A new method in vibration analysis by BEM using particular integrals. J Engng Mech Div, ASCE 1986;113. [3] Zheng R, Phan-Thien N, Coleman CJ. A boundary element approach for non-linear boundary value problems. Comput Mech 1991;8:71– 86. [4] Power H, Pratridge PW. The use of Stokes’ fundamental solution for the boundary only element formulation of the three-dimensional

H. Power, R. Mingo / Engineering Analysis with Boundary Elements 24 (2000) 107–119

[5]

[6] [7] [8] [9] [10]

[11] [12] [13]

[14] [15] [16]

[17] [18]

Navier–Stokes equations for moderate Reynolds numbers. Int J Numer Methods Engng 1994;37:1825–40. Ch. Rahain P, Kassab A. Pressure correction DRBEM solution for heat transfer and fluid flow in incompressible viscous fluids. Engng Anal 1996;18:265–72. Banerjee PK. The boundary element methods in engineering, London: McGraw-Hill, 1981. Popov V, Power H. A domain decomposition in the dual reciprocity approach. Bound Elem Commun 1996;7/1:1–5. Kansa EJ, Carlson RE. Radial basis functions: A class of grip-free, scattered data approximations. Comput Fluid Dyn J 1995;3/4:479–96. Ladyzhenskaya O. The mathematical theory of viscous incompressible flow, New York, London: Gordon and Breach, 1963. Micchelli CA. Interpolation of scattered data: distance matrices and conditionally positive definite functions. Constr Approx 1986;2:11– 22. Franke R. Scattered data interpolation: tests of some methods. Math Comp 1982;38:181–200. Stead S. Estimation of gradients from scattered data. Rocky Mount J Math 1984;14:265–79. Duchon J. In: Schempp W, Zeller K, editors. Spline minimizing rotation-invariant seminorms in Sobolev spaces, Constructive theory of functions of several variables, 571. Berlin: Springer-Verlag, 1977. pp. 85–100. Powell MJD. The uniform convergence of thin plate spline interpolation in two dimensions. Numerische Mathematik 1994;68(1):107–28. Madych WR, Nelson SA. Multivariable interpolation and conditionally positive definite functions. II Math Comput 1990;54:211–30. Ghia U, Ghia KN, Shin N. High-Re solutions for incompressible flow using the Navier–Stokes equations and a multi-grid. J Comput Phys 1982;48:387–411. Moffatt HK. Viscous and resistive eddies near a sharp corner. J Fluid Mech 1964;18:1–19. Taylor GI. On scraping viscous fluid from a plane surface. In: Schafer

[19] [20] [21] [22] [23] [24]

[25]

[26]

[27]

[28]

[29]

119

M, editor. Miszellaneen Ang. Mech., Berlin: Akademie Verlag, 1962. p. 313–5. Higdon JJL. Stokes flow in arbitrary two-dimensional domains: shear flow over ridges and cavities. J Fluid Mech 1985;159:195–226. Pozrikidis C. The flow of a liquid film along a periodic wall. J Fluid Mech 1988;188:275–300. Kelmanson MA. An integral equation method for the solution of singular slow flow problems. J Comp Phys 1983;15(1):139–58. Kelmanson MA. Boundary integral equation solution of viscous flows with free surfaces. J Engng Maths 1983;17:329–43. Pozrikidis C. Introduction to theoretical and computational fluid dynamics, Oxford: Oxford University Press, 1997. Ingham DB, Kelmanson MA. Boundary integral equation analyses of singular, potential, and biharmonic problems, Berlin: SpringerVerlag, 1984. Kitagawa K, Brebbia CA, Wrobel LC, Tanaka M. Boundary element analysis of viscous flow by penalty function formulation. Engng Anal 1986;3:194–200. Denham MK, Patrick MA. Laminar flow over a downstream-facing step in a two-dimensional flow channel. Trans Inst Chem Eng 1974;52(4):361–7. Zienkiewicz OC, Satya Sai BVK, Morgan K, Codina R. Split, characteristic based semi-implicit algorithm for laminar/turbulent incompressible flows. Int J Numer Methods Fluids 1996;23: 787–809. Taylor C, Thomas CE, Morgan K. In: Taylor C, Morgan K, editors. Analysis of turbulent flow with separation using the finite element method, recent advances in numerical methods in fluids, Computational techniques in transient and turbulent Flow, 2. Swansea: Pineridge, 1981. p. 283–325. Skerget A, Alujevic A, Brebbia CA. Analysis of laminar fluid by boundary elements, numerical methods in laminar and turbulent flow Part 1. In: Taylor, editor. Proceedings of the 4th international conference, 1985.