Development of a meshless hybrid boundary node method for Stokes flows

Engineering Analysis with Boundary Elements 37 (2013) 899–908

Contents lists available at SciVerse ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Development of a meshless hybrid boundary node method for Stokes flows Fei Tan n, Youliang Zhang, Yinping Li State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China

art ic l e i nf o

a b s t r a c t

Article history: Received 5 November 2012 Accepted 18 March 2013 Available online 26 April 2013

The meshless hybrid boundary node method (HBNM) is a promising method for solving boundary value problems, and is further developed and numerically implemented for incompressible 2D and 3D Stokes flows in this paper. In this approach, a new modified variational formulation using a hybrid functional is presented. The formulation is expressed in terms of domain and boundary variables. The moving leastsquares (MLS) method is employed to approximate the boundary variables whereas the domain variables are interpolated by the fundamental solutions of Stokes equation, i.e. Stokeslets. The present method only requires scatter nodes on the surface, and is a truly boundary type meshless method as it does not require the ‘boundary element mesh’, either for the purpose of interpolation of the variables or the integration of ‘energy’. Moreover, since the primitive variables, i.e., velocity vector and pressure, are employed in this approach, the problem of finding the velocity is separated from that of finding pressure. Numerical examples are given to illustrate the implementation and performance of the present method. It is shown that the high convergence rates and accuracy can be achieved with a small number of nodes. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Stokes equations Modified variational principle Moving least-squares approximation Stokeslets Hybrid boundary node method Meshless

1. Introduction Development of efficient computational approaches for the numerical simulation of the steady and the unsteady viscous fluid flow has become a very important need in various fields of practical engineering. The numerical solution of the Stokes equations for incompressible viscous fluids has been usually addressed with classical numerical methods such as finite element method (FEM) [1,2], finite difference method (FDM) [3], boundary element method (BEM) [4,5], et al. However, their reliance on a mesh leads to complications for certain classes of problems. When the geometry of the problem becomes complicated, mesh generation, modification and re-meshing are often arduous and timeconsuming (especially in 3D). These difficulties can be overcome by the so-called meshless methods, which have attracted more and more attention during the past decades. Compared with mesh-based methods, the meshless methods avoid mesh generation and the approximation solution is constructed entirely based on a cluster of scatter nodes, and nodes can be easily added and removed without a burdensome remeshing of elements. So far a number of meshless methods have been developed by different authors, and some of them are used to solve the viscous flows, e.g. the smooth particle hydrodynamics (SPH) [6], the element free Galerkin method (EFGM) [7], the finite

n

Corresponding author. Tel./fax: +86 27 8719 8783. E-mail address: [email protected] (F. Tan).

0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.03.012

point method (FPM) [8,9], the meshless local Petrov–Galerkin method (MLPG) [10,11], the local boundary integral equation method (LBIE) [12], the differential quadrature [13,14], the radial point interpolation method (RPIM) [15–17], the method of fundamental solutions (MFS) [18–20], the Galerkin boundary node method (GBNM) [21–24] and so on. Some of these methods, in reality, are not real meshless methods, since they need background meshes for the numerical integration. Zhang et al. [25] proposed a boundary type truly meshless method—hybrid boundary node method (HBNM). It combines the moving least squares (MLS) interpolation scheme with the hybrid displacement variational formulation, and the integration is limited to a fixed local region on the boundary. Only scatter nodes on boundary of the domain are required. The HBNM has been employed to solve a wide range of problems, such as the potential problems [25], the elastostatics problems [26–28] and so on. Miao et al. [29] and Yan et al. [30] introduced the dual reciprocity method (DRM) into HBNM, and extended this method to solve the elastodynamics problems. Recently, Wang et al. [31,32] implemented the fast multipole method (FMM) to accelerate the HBNM for 3D elasticity problems. In this paper, the solution of incompressible Stokes flows will be addressed. In general, incompressible Stokes flows can be solved by using either primitive (i.e. velocity and pressure) or derived (such as vorticity and stream function) variables as primary unknowns. For 2D flows, one can transform the Stokes equations into the biharmonic equations using the steam function formulation [23], or recast the Stokes equations in velocity–vorticity form to the

900

F. Tan et al. / Engineering Analysis with Boundary Elements 37 (2013) 899–908

Laplace equation and Poisson equation [33]. Both these methods will have no immediate solutions of the pressure, and are not generally applicable to three-dimensional problems. Consequently, the 3D Stokes problems are most often solved by the methods based on primitive variables. Even for 2D problems, the use of primitive variables is quite common, especially when robust codes are needed for various applications. Therefore, methods based on primitive variables are considered here. The present paper develops and employs the HBNM for 2D and 3D Stokes flows in the velocity–pressure formulation. The discrete equations are obtained from a variational principle using a modified hybrid functional based on four independent variables, i.e. velocities and tractions on the boundary and velocities and pressure inside the domain. In this approach, the boundary variables are approximated by MLS, whereas the domain variables are interpolated by the fundamental solutions of the Stokes equations, i.e. Stokeslets. And the problem of finding the velocity is separated from that of finding the pressure. The accuracy and efficiency of the present numerical scheme are demonstrated by some numerical examples. The plan of this paper is as follows: following the first introductory section, the governing equations of incompressible viscous fluid flow and the fundamental solutions of Stokes equations are described in Section 2. The formulations of HBNM for 2D and 3D incompressible Stokes flows are developed in Section 3. Subsequently some numerical examples are given to evaluate the effectiveness and validity of the method are presented in Section 4. Finally, the paper will end with conclusions in Section 5.

2. Governing equations and Stokeslets Let us consider the steady-state problems of a Newtonian, viscous and incompressible fluid flow. The region of the flow contained in Euclidean space R2 is denoted by Ω with the boundary Γ. The steady-state flow problem is governed by the following fundamental equations and conditions in the dimensionless form: Equations of linear momentum sij;j ¼ 0

in Ω

ð1Þ

Incompressibility condition ui;i ¼ 0

in Ω

ð2Þ

Boundary conditions ui ¼ ui

on Γ u

t i ≡sij nj ¼ t i

on Γ t

ð3Þ ð4Þ

where ui is the velocity vector component, ni is the unit outward normal vector component. ui and t i are the prescribed velocities and tractions on the essential boundary Γ u and on the traction boundary Γ t , respectively. sij is the Cauchy stress tensor defined as sij ¼ −pδij þ 2μεij

ð5Þ

r ;k 2πr

ð8Þ

r ;i r ;j r ;k πr

ð9Þ

pIk ¼ − sIijk ¼

t Iik ≡sIijk nj ¼

r ;i r ;k r ;j nj πr

ð10Þ

where r ¼ jx−yI j is the distance between field point x and source point yI . Similarly, the Stokeslets for 3D flow can be obtained as [18]:  1  δ þ r ;i r ;k 8πμr ik

ð11Þ

r ;k 4πr 2

ð12Þ

3r ;i r ;j r ;k 4πr 2

ð13Þ

uIik ¼ − pIk ¼ − sIijk ¼

t Iik ≡sIijk nj ¼

3r ;i r ;k r ;j nj 4πr 2

ð14Þ

3. Development of HBNM for stokes equations 3.1. Variational formulation In this paper, the HBNM is based on a modified variational principle. For the Stokes flows problems, the functions to be independent are: velocity field in the domain, ui ; pressure in the domain p; boundary velocity field, u~ i ; boundary tractions t~ i . The modified variational functional is defined as  Z
1 μðui;j þ uj;i Þui;j d Ω Π HB ¼ Ω 2 Z Z Z ð15Þ − t i u~ i d Γ− pui;i d Ω− t~ i ðui −u~ i Þ d Γ Γt

Ω

Γ

In the above equation, the boundary velocity u~ i satisfies the essential boundary conditions, i.e. u~ i ¼ ui on Γ u . Performing variations of Eq. (15), integrating by parts and taking the constitutive equation and essential boundary conditions on Γ u into account, we have Z Z Z δΠ HB ¼ − sij;j δ ui d Ω− ui;i δp d Ω þ ðt i −t~ i Þδui d Γ Ω Γ Z ZΩ ~ ~ ð16Þ − ðui −u~ i Þδ t i d Γ− ðt i −t i Þδ u~ i d Γ Γ

Γt

The vanishing of δΠ HB for arbitrary variations δ ui and δp in Ω, δ u~ i and δ t~ on Γ, with δ u~ i ¼ 0 on Γ u , gives the following Euler equations: sij;j ¼ 0

in Ω

ð17Þ

ui;i ¼ 0

in Ω

ð18Þ

t i −t~ i ¼ 0

on Γ

ð19Þ

where p is the pressure, μ is the constant positive viscosity coefficient, δij is the Kronecker delta and

ui −u~ i ¼ 0

on Γ

ð20Þ

1 εij ¼ ðui;j þ uj;i Þ 2

t i −t~ i ¼ 0

on Γ t

ð21Þ

ð6Þ

is the deformation tensor. The fundamental solutions for the Stokes operators are called the Stokeslets, which represent the flow fields due to concentrated point forces. The Stokeslets for 2D flow can be obtained as [18]:
 1 1 δik ln þ r ;i r ;k uIik ¼ − ð7Þ 4πμ r

Consequently, the solution of the problem is now given in terms of the functions ui , p, u~ i and t~ i , which makes δΠ HB stationary. Let δΠ HB ¼ 0, the following equivalent integration equations can be obtained as Z Z − sij;j δ ui d Ω þ ðt i −t~ i Þδui d Γ ¼ 0 ð22Þ Ω

Γ

F. Tan et al. / Engineering Analysis with Boundary Elements 37 (2013) 899–908

Z Γ

ðui −u~ i Þδ t~ i d Γ ¼ 0

Z Γt

ð23Þ

ðt i −t~ i Þδ u~ i d Γ ¼ 0

ð24Þ

If the traction boundary condition, t~ i ¼ t i on Γ t is imposed, Eq. (24) will be automatically satisfied. So it can be ignored in the following analysis. It is clear that the modified variational formulation can be satisfied in any arbitrary sub-domain. The same as HBNM for potential problems [25], the sub-domain Ωs is chosen as the intersection of the domain Ω and a circle (or a sphere) centered at a boundary node sJ with the radius is r J . The test function vJ ðQ Þ is used to replace the variational part. In this case, the integral equations become Z Z − sij;j vJ ðQ Þ d Ω þ ðt i −t~ i ÞvJ ðQ Þ d Γ ¼ 0 ð25Þ Ωs

Γ s þLs

Z Γ s þLs

ðui −u~ i ÞvJ ðQ Þ d Γ ¼ 0

ð26Þ

It can be seen that the variables u~ i and t~ i on Ls in Eqs. (25) and (26) are not defined. If vJ ðQ Þ is selected in such a way that all integrals over Ls vanishes, the problem can be solved conveniently. This can be easily accomplished by using the following Gauss function as the test function 8 2 2 J =cJ Þ  < exp½−ðdJ =cJ Þ −exp½−ðr 0 ≤dJ ≤r J 1−exp½−ðr J =cJ Þ2  vJ ðQ Þ ¼ ð27Þ :0 dJ ≥r J where dJ is the distance between the integral point Q in the domain and the nodal point sJ , and cJ is a constant controlling the test function shape. From Eq. (27) it can be seen that vJ ðQ Þ vanishes on the boundary LS . Therefore, Eqs. (25) and (26) can be rewritten as follows: Z Z − sij;j vJ ðQ Þ d Ω þ ðt i −t~ i ÞvJ ðQ Þ d Γ ¼ 0 ð28Þ Ωs

Z Γs

they can be expressed as ( ) " I #( I ) N u11 uI12 x1 u1 u¼ ¼ ∑ I uI22 xI2 u2 I ¼ 1 u21 ( t¼

)

t1 t2

The above equations are the boundary local integral equations for Stokes flows problems.

N

I u~ i ¼ ∑ ΦI u^ i I¼1 N

I¼1

3

31

)

xI1

ð33Þ

xI2

32

33

38 I 9 t I13 < > x1 > = I t I23 7 5 x2 > : xI > ; tI

t I12 t I22 t I32

33

ð34Þ

3

ð35Þ

3

3.3. HBNM formulations for stokes flows As u is expressed by Eqs. (32) and (34), the term sij;j in the left hand side of Eq. (28) vanished. For 2D flows, by substituting Eqs. (30)–(33) into Eqs. (28) and (29), and omitting the vanished terms, we have #( I ) Z " I N x1 t 11 t I12 ∑ vJ ðQ Þ dΓ I xI2 t I22 I ¼ 1 Γ s t 21 #8 I 9 Z " N ΦI 0 < t^ 1 = ¼ ∑ ð36Þ v ðQ Þ dΓ 0 ΦI : t^ I ; J I ¼ 1 Γs 2 N

∑

Z " uI 11

uI12

uI21

uI22

Γs

I¼1

N

Z "

I¼1

Γs

ΦI 0

#(

xI1

)

vJ ðQ Þ dΓ xI2 #8 I 9 0 < u^ 1 = v ðQ Þ dΓ ΦI : u^ I2 ; J

ð37Þ

And for 3D flows, by substituting Eqs. (30), (31), (34) and (35) into Eqs. (28) and (29), and omitting the vanished terms, we have 2 I 38 I 9 t 11 t I12 t I13 > Z < x1 > = N 6 tI I I 7 I ∑ 4 21 t 22 t 23 5 x2 vJ ðQ Þ dΓ > : xI > ; I ¼ 1 Γs tI tI tI N

Z

¼ ∑

I¼1

32

2

ΦI 6 0 4 Γs 0

33

0 ΦI 0

3

8 9 3> t^I > > 0 > > = < 1> I ^ 0 7 5 t 2 vJ ðQ Þ dΓ > > > > ΦI > ; : t^I >

ð38Þ

3

Z

N

∑

2

uI11 6 uI 4 21 Γs uI31 N

where N stands for the number of nodes located on the boundary; I I u^ i and t^i are the fictitious nodal values; ΦI is the shape function of the MLS, which can be seen in Ref. [25,28]. The domain variables ui and t i are interpolated by the fundamental solutions of Stokes equations, i.e. Stokeslets. For 2D flows,

t I22

#(

where xI1 , xI2 and xI3 are interpolation coefficients, uIij and t Iij are the fundamental solutions with source point at a point P I , the specific expressions are Eqs. (7), (10), (11) and (14).

I¼1

ð31Þ

t I21

8 9 2 I t 11 > < t1 > = N 6 tI t ¼ ∑ 4 21 t¼ 2 > :t > ; I ¼ 1 tI

¼ ∑ I t~ i ¼ ∑ ΦI t^ i

t I12

31

31

ð30Þ

t I11

3

3.2. Variables interpolation As a truly boundary-type meshless method, the HBNM uses the MLS to approximate the boundary variables, and applies the fundamental solution interpolation to obtain the solutions in the domain. By using the MLS principle, the boundary variables u~ i and t~ i can be written as

"

I¼1

¼ ∑ ð29Þ

N

¼ ∑

ð32Þ

And for 3D flows, they can be expressed as 8 9 2 I 38 I 9 u11 uI12 uI13 > > < u1 > = < x1 > = N 6 7 u ¼ u2 ¼ ∑ 4 uI21 uI22 uI23 5 xI2 > > > : u ; I ¼ 1 uI : xI > ; uI uI

Γs

ðui −u~ i ÞvJ ðQ Þ d Γ ¼ 0

901

I¼1

Z

uI12 uI22 uI32 2

ΦI 6 0 4 Γs 0

38 I 9 uI13 < > x1 > = I 7 u23 5 xI2 vJ ðQ Þ dΓ > : xI > ; uI 33

0 ΦI 0

3

38 ^ I 9 u > 0 > > = < 1> 7 0 5 u^ I2 vJ ðQ Þ dΓ > > ; : u^ I > ΦI >

ð39Þ

3

Applying the above equations for all boundary nodes, the final discrete system of equations is obtained Tx ¼ Ht^

ð40Þ

^ Ux ¼ Hu

ð41Þ

902

F. Tan et al. / Engineering Analysis with Boundary Elements 37 (2013) 899–908

where H IJ ¼

Z "Φ J

0

0

ΦJ

t J11

t J12

t J21

t J22

Γ Is

2

Z T IJ ¼

4

Γ Is

2

Z

2

# vI ðQ Þ dΓ

3 5 vI ðQ Þ dΓ

uJ11

uJ12

uJ21

uJ22

uJ11 6 J 6 u21 U IJ ¼ 4 Γ Is uJ31

uJ12

U IJ ¼

4 Γ Is

Z

2

for 2D;

ΦJ Z 6 0 H IJ ¼ 4 Γ Is 0

0

2

t J11 6 J 6t T IJ ¼ 4 21 Γ Is t J31 Z

for 2D;

0 ΦJ

3 0 0 7 5vI ðQ Þ dΓ ΦJ

t J12

t J13

t J22 t J32

t J23 t J33

for 3D

3 7 7vI ðQ Þ dΓ 5

for 3D

4. Numerical experiments

3 5 vI ðQ Þ dΓ uJ13

uJ32

In this section, some numerical examples are presented to demonstrate the accuracy and efficiency of the proposed method. For the purpose of error estimation and convergence studies, a Euclidean norm error, normalized by jumax j is defined as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  1 N 2 ð48Þ e¼ ∑ ðuðeÞ −uðnÞ  i Þ jumax Ni¼1 i

for 2D;

3

7 uJ23 7 5vI ðQ Þ dΓ

uJ22

for 3D

uJ33

ð42Þ

where jumax j is the maximum value of uover N sample points, the superscripts ðeÞ and ðnÞ refer to the exact and numerical solutions, respectively. For simplicity, we always assume that the viscosity of the fluid is taken as one, i.e. μ ¼ 1. In all computations, unless otherwise specified, the radius of the sub-domain r J in Eq. (27) is chosen as 0:8h, and the parameter cJ is taken to be such that r J =cJ ¼ 1:2. In order to deal with the discontinuities at the corners, the nodes are not arranged at these places and the support domain for interpolation is truncated.

ð43Þ

4.1. A square cavity

N T x ¼ ½x11 ; x12 ; …; xN 1 ; x2 

for 2D;

N N T x ¼ ½x11 ; x12 ; x13 ; …; xN 1 ; x2 ; x3 

for 3D

1 1 N N t^ ¼ ½t^ 1 ; t^2 ; …; t^1 ; t^ 2 T

for 2D;

1 1 1 N N N t^ ¼ ½t^1 ; t^ 2 ; t^ 3 ; …; t^ 1 ; t^2 ; t^3 T

for 3D

^ ¼ ½u^ 11 ; u^ 12 ; …; u^ N ^N T u 1 ; u2 

for 2D;

^ ¼ ½u^ 11 ; u^ 12 ; u^ 13 ; …; u^ N ^N ^N T u 1 ; u2 ; u3 

for 3D

Γ Is

with denotes the local boundary around node sI . Rearranging Eq. (41), we have ^ x ¼ U−1 Hu

It can be seen from above derivations that the present method is a truly meshless method, as absolutely no boundary elements are needed, either for the interpolation purpose or for the integration purpose. Besides, no further integration is needed in the ‘post-processing’ step.

Substitute Eq. (42) into Eq. (40), we have ^ t^ ¼ 0 TU Hu−H −1

The equations can be solved in the same manner as the conventional BEM, except that transformations between u^ i and u~ i , t^ i and t~ i must be performed, due to that the shape functions of MLS lack the delta functional property. For any well-posed problem, either u~ i or t~ i is known on the boundary, so the corresponding fictitious nodal values u^ i and t^i can be obtained as follows: N

N

J¼1

J¼1

I u^ i ¼ ∑ RIJ u~ Ji ¼ ∑ RIJ uJi

N

N

J¼1

J¼1

ð44Þ

I J J t^i ¼ ∑ RIJ t~ i ¼ ∑ RIJ t i

ð45Þ J

where RIJ ¼ ½ΦJ ðsI Þ−1 , and uJi and t i are the prescribed solutions of boundary node J, respectively. In order to avoid the ‘boundary layer effect’ and singular integration, the source points are arranged outside boundary. And P I is determined by P I ¼ SI þ h ξ nðSI Þ

First we will consider a simple case, on a square domain Ω ¼ ½−1; 1  ½−1; 1 with an analytical solution u1 ¼ 4y3 −x2 ; u2 ¼ 4x3 þ 2xy−1; p ¼ 24xy−2x where the essential boundary conditions are imposed on all boundaries. The tractions along the boundary in a counterclockwise path from the starting point (−1, 1), in the case that five uniformly spaced nodes are used on each edge of the boundary, together with the exact solutions are shown in Figs. 1 and 2. And Figs. 3 and 4 display a comparison between the numerical results with the exact solutions for velocity, pressure and stress on the diagonal (from (−1, −1) to (1, 1), 19 uniformly spaced sample points). These figures show that the numerical solutions are in excellent agreement with the exact solutions with just a few nodes.

ð46Þ

whereh as the nodal spacing; nðSI Þ is the outward normal direction to the boundary at node SI ; ξ is the scale factor, and is chosen as 6.0 according to Ref. [26]. Apply the boundary conditions into Eq. (43), rearrange them, and obtain final system equations, then the unknown variables on the boundary can be solved. After that, the unknown vector x can be obtained by Eq. (42). The various variables inside domain can be evaluated by Eqs. (32)-(35) without further integrations. Similarly, the pressure can be calculated by N

p¼ ∑

h

I¼1 N

¼ ∑

I¼1

h

pI1

pI2

pI1

pI2

( ) i xI 1

xI2 pI3

in 2D or p

8 I 9 < x1 > = i> xI2 > : xI > ; 3

in 3D

ð47Þ Fig. 1. Traction t1 distribution around the square boundary.

F. Tan et al. / Engineering Analysis with Boundary Elements 37 (2013) 899–908

903

convergence of velocity u1 , pressure p and stress s1 on the diagonal are shown in Fig. 5. It can be observed that the convergence rates of the present method are very high, the convergence order of velocity is about 3.26, and is higher than that of the pressure and stress, 2.69. The values of the numerical solutions of velocities and pressure at point (−0.6, −0.2) are displayed in Table 1. For comparison, results from analytical method and the meshless Galerkin boundary node method (GBNM) [23] are also listed. From this table, we can find that the numerical solutions converge to their corresponding analytical solution as the number of boundary nodes increase, and the accuracy of the proposed method is higher than that of the GBNM. 4.2. A circular cavity

Fig. 2. Traction t2 distribution around the square boundary.

This problem is the well-known lid-driven Stokes flow in a circular cavity, which has also been described in papers [22,23]. The flow domain considered here is the unit circle as shown in Fig. 6. The upper half of the boundary is assumed to move with a unit tangential velocity in the counterclockwise direction, while the lower half is set with no-slip and impervious Dirichlet conditions. The analytical solution for the velocity is [22,23]
  y π 2y ð1−x2 −y2 Þð1−x2 þ y2 Þ u1 ¼ − þ tan−1 þ 2 2 π 2 1−x −y π½ð1−x2 −y2 Þ2 þ 4y2 
  x π 2y 2xyð1−x2 −y2 Þ þ tan−1 − u2 ¼ 2 2 π 2 1−x −y π½ð1−x2 −y2 Þ2 þ 4y2  Fig. 7 plots the results of velocities on the vertical centerlines of the cavity. In this analysis, 32 regularly and randomly distributed boundary nodes are used. It can be seen that the expressions of the analytical solutions are very complex, but the numerical results

Fig. 3. Velocity distribution along the diagonal (from (−1, −1) to (1, 1)).

Fig. 5. Relative errors and convergence rates for the Stokes flows on a square.

Table 1 The values of the numerical solutions of velocities and pressure at point (−0.6, −0.2). Exact solutions

Numerical solutions N¼8 Fig. 4. Distribution of pressure and stress along the diagonal (from (−1, −1) to (1, 1)).

To study the convergence of the proposed method, five different regular nodes arrangements of 5, 10, 20, 40 and 80 nodes on each edge have been used. The results of relative errors and

N ¼ 16

N ¼ 32

u1 GBNM [23] 0.0029 −0.3369 −0.3857 Present method −1.1600 −0.3959 −0.3926 u2 GBNM [23] −1.7255 −1.5929 −1.6156 Present method −3.1600 −1.6382 −1.6255 p GBNM [23] 5.2635 4.0600 4.0881 Present method 1.1997 4.0656 4.0778

N ¼64 −0.3890 −0.3920 −1.6199 −1.6241 4.0839 4.0799

−0.3920 −1.6240 4.0800

904

F. Tan et al. / Engineering Analysis with Boundary Elements 37 (2013) 899–908

y

u

1 ur

0

r

ur

u

x

0

Fig. 6. Circular cavity flow problem with boundary conditions.

Fig. 8. Relative errors and convergence rates for the lid-driven circular cavity flow.

1

y

0.5

0

-0.5 Fig. 7. Velocity u for the lid-driven circular cavity flow along the y-axis.

still capture their behavior very well. The accuracy obtained with the regularly distributed boundary nodes is slightly higher than with the randomly distributed boundary nodes. When four different regular nodes arrangements of 8, 16, 32 and 64 nodes have been used, the relative errors and convergence of velocity u1 and stress sxy are shown in Fig. 8. For comparison, results from the FEM are also depicted in this figure. During the finite element analysis, the domain is discretized using 20,000 elements. It can be seen that the convergence of the present method is obviously higher than that of the FEM. Moreover, the velocity vector is observed in Fig. 9, and the velocity component u1 and u2 contours are depicted in Figs. 10 and 11, respectively. The present computations give very good results as compared with the other solutions [15,22], and show the correct fluid motion, which demonstrates the convergence and effectiveness of the proposed approach. 4.3. Flow between two rotating cylinders The third example considered here is a Stokes flow confined in two rotating coaxial cylinder. The configuration and boundary conditions of the problems are depicted in Fig. 12. The radii of the inner and outer cylinders are R1 ¼ 0:5 and R2 ¼ 1:5, respectively. No-slip boundary conditions are specified on the fluid–solid interfaces, the inner cylinder rotates clockwise with U 1 ¼ −1, and the outer cylinder rotates counterclockwise with U 2 ¼ 2. The analytical solutions for this problem are [22] uθ ¼ c1 =r þ c2 r; ur ¼ 0; p ¼ p0 ; r∈½R1 ; R2 

-1 -1

-0.5

0

x

0.5

1

Fig. 9. Velocity vectors for the lid-driven circular cavity flow.

where c1 ¼ ½R1 R2 ðU 1 R2 −U 2 R1 Þ=ðR22 −R21 Þ, c2 ¼ ðU 2 R2 −U 1 R1 Þ=ðR22 − R21 Þ and p0 is a constant. In this study, 32 uniformly spaced nodes are used on each cylinder surface. Figs. 13 and 14 show a comparison between the numerical solutions and the exact solutions for the velocities, pressure and stresses along the x-axis. Again, it can be clearly seen that the HBNM solutions are in excellent agreement with the analytical solutions. It is worth noting that the pressure is only determined up to a constant, and in this case p0 ¼ 0, which is the same to the result of [22]. Relative errors and convergence rates of velocity u1 and stress sxy are plotted in Fig. 15. It is also shown that the convergence rates of the present method are very high, and velocity converges a little faster than the pressure and stress. This may be due to that the singularities of stresses in Eq. (9) are stronger than the velocities in Eq. (7). The same phenomenon can be found in the GBNM Ref. [21], Table 1. The velocity vector is shown in Fig. 16, and the velocity component u1 and u2 contours are depicted in Figs. 17 and 18, respectively. It can be observed that the flow patterns are

F. Tan et al. / Engineering Analysis with Boundary Elements 37 (2013) 899–908

1 -0.8

-0.6 -0. 5 -0.4 -0.3

0.5

y

-0.6 -0.5 -0.4 -0.3

-0.8 -0.7

-0.1

0

-0. 2

0.1

0.1

0.2

0.2

0

0.3

0.3

0

0.1

0.3

-0.2

-0.1

2 -0. 0

0

-0.9

-0.7

2

0.1

0

0. 1

0. 3

0

0

0.2 0

-1 -1

0. 2

0

0.3

0.

0

0

-0.5

905

0.1

-0.5

0

0.5

x

1

Fig. 13. Velocities distribution along the x-axis for Stokes flow between two rotating coaxial cylinders.

Fig. 10. Velocity u1 contour for the lid-driven circular cavity flow.

-0.2

- 0. 4

1

0.2

0.5 6

0.8

0. 4

-0.

0

0

-0.4

-0.2

y

0.6

02

0

0. 2

0

0

-0.5

0

-1 -1

-0.5

0

0

x

0.5

1

Fig. 14. Distribution of pressure and stresses along the x-axis for Stokes flow between two rotating coaxial cylinders.

Fig. 11. Velocity u2 contour for the lid-driven circular cavity flow.

Fig. 12. Schematic diagram for Stokes flow between two rotating coaxial cylinders.

Fig. 15. Relative errors and convergence rates for Stokes flow between two rotating coaxial cylinders.

906

F. Tan et al. / Engineering Analysis with Boundary Elements 37 (2013) 899–908

1

0.5

0.5

0

0

-0.2

1

0 0. 6

0. 8

1.5 0. 4

1.5

1

0

0.2

0.6 0.8

-0.6 -0.2 -0.4 0 0. 4 0 . 2 0. 6

-1.5 -1.5

-1

-0.5

0

0.5

x

1

1.5

Fig. 16. Velocity vectors for Stokes flow between two rotating coaxial cylinders.

-1.5 -1.5

-1

-0.5

0

x

0.2

0

2

-0.8

. -0

0

0.4

y

-0.2

-1

y

0.4

.4 -0-0.6

-0.8

-1

2 0.

0.8 4

-1

-0 .8

0.2

- 0.

-1

-0.6

-0.5

-1

-0.5

0

0.5

1

1.5

Fig. 18. Velocity u2 contour for Stokes flow between two rotating coaxial cylinders.

1.5 -1

-0.8

-0.6

-0 . 8 0

0. 6

1

0.5 0

0

8

-0.6 -0.4

0.2 0

0.6

0.4

0.8

-0. 2

-0.

-1

2 -0.

0. 4

-1

-1.5 -1.5

4

0

0.2

-0.5

-0. .2

Fig. 19. Schematic diagram for Stokes flow in a rectangular cavity with waveshaped bottom.

0

y

-0

0.8

-0.2

0

-1

-0.6

-0.4 -0.2 0.2 0.4

1

0.6

0.2

0. 4

0.8

1 1.2

-1

-0.5

0

x

0.5

1

1.5

Fig. 17. Velocity u1 contour for Stokes flow between two rotating coaxial cylinders.

symmetric with respect to the x-axis and y-axis due to the reversibility of the Stokes flow. 4.4. Rectangular cavity with wave-shaped bottom The capability of the proposed numerical scheme to handle irregular computational domains is demonstrated by solving a liddriven rectangular cavity with wave-shaped bottom. The parametric representation of the wave-shaped bottom is given as follows 1 0   1−x−2:5 A; 0 ≤x ≤ 5 y ¼ −0:3sinðπxÞ@ 2:5 The configuration and boundary conditions of the problem are shown in Fig. 19. In this simulation, 80 nodes on the boundary are employed. The distribution of velocity component u1 on the vertical centerline x¼2.5 is depicted in Fig. 20. The results are compared with a GBNM solution obtained using 140 nodes [22], and a MFS solution

Fig. 20. Comparison of velocity u1 profile along line x¼ 2.5 for the rectangular cavity with wave-shaped bottom.

obtained using 420 nodes [19]. Although the analytical solutions for such a complex geometry are not available, the results predicted by the HBNM are in close agreement with the results of other algorithms. This demonstrates that the HBNM is found to be computationally efficient to handle the irregular shaped domains. Figs. 21–23 depict the velocity vector, velocity component u1 and u2 contours, respectively. The numerical results indicate that the flow symmetry pattern still remains as expected for the Stokes flow.

F. Tan et al. / Engineering Analysis with Boundary Elements 37 (2013) 899–908

907

2

y

1.5 1 0.5 0 0

0.5

1

1.5

2

2.5

3

x

3.5

4

4.5

5

Fig. 21. Velocity vectors for the rectangular cavity with wave-shaped bottom.

0. 1

y

1.5 0

0.7 0.1 0.30 .2 0

-0.2

-0.1

0

0.5

1

1.5

2

-0.2

0. 1

0.6

-0.1

0

2.5

x

Fig. 25. Comparison of velocity u3 profile along x at y¼ 0.5 and z¼ 0.5 for a cubic cavity.

0

3

1

-0.1

0

0

0. 2 0

0.4

-0.3

-

.1 -0

0.5

0.7 0.30.2 0.1 0

r

0.3

0.8 0.5

0.9

0.6

-0.1 -0.2

1

0 .5

0.8 0.5

0.4

01

2 0.1

0

3.5

4

4.5

5 0.8

0. 0 5

0.1 5 01 0

-0 .

0.4

05

3

0

0 -0.05

0.

0

0.5

1

1.5

2

2.5

x

3

3.5

4

0.6

z

0

0

5

-0.2 . --0 0.3 2 -0 -0.1 - 0. 0

.2 -0 0.1

0

-0.2

0

y

0. 05 0 .0

5 0.0 1 . 0

0.05

05

.1 5 -0 -0.1

1 0.5

-0.

-0.3 -0 .2-05 0.05 0 3

. 15

0.1

1.5

0.05

0.10

1 0.

00.25

30 .35 .15

-0.

0.

0.2

2

-0.0

Fig. 22. Velocity u1 contour for the rectangular cavity with wave-shaped bottom.

4.5

0.2

5

Fig. 23. Velocity u2 contour for the rectangular cavity with wave-shaped bottom.

0

0

0.2

0.4

x

0.6

0.8

1

Fig. 26. Velocity vectors for the cubic cavity flow at the plane of y¼ 0.5.

cube is assumed to move with a unit velocity in the horizontal x1 direction. Others are set with no-slip and impervious Dirichlet conditions. Twenty-five nodes, uniformly spaced, are used on each surface of the cube. Figs. 24 and 25 show the comparisons of the velocity component u1 profiles along the vertical center lines and u3 profiles along the horizontal center lines between the present method and the MFS using the same boundary node arrangement. It can be found that close agreement between the present results and the results of the MFS is achieved. The distributions of velocity vector u1 −u3 on x–z plane at y¼ 0.5 are illustrated in Fig. 26. This figure exhibits the symmetric characteristics of the flow variables as expected for the Stokes flow and thus shows the correct fluid motion. Fig. 24. Comparison of velocity u1 profile along z at x¼ 0.5 and y¼ 0.5 for a cubic cavity.

4.5. A lid-driven cubic cavity flow For the 3D application, the well-known lid-driven Stokes flow in a ð0; 0; 0Þ  ð1; 1; 1Þ cubic cavity is considered. The top lid of the

5. Conclusions In this work, the hybrid boundary node method has been successfully extended to 2D incompressible Stokes flows. In comparing with analytical solutions and other numerical solutions, the results show that the present scheme is efficient and accurate for using a few boundary nodes.

908

F. Tan et al. / Engineering Analysis with Boundary Elements 37 (2013) 899–908

The present method is based on a modified variational functional. The boundary variables are approximated by the MLS scheme, and the domain variables are interpolated by the fundamental solutions of Stokes equations, i.e. Stokeslets. Compared with the Galerkin boundary node method (GBNM), the present method is a truly meshless method, only randomly distributed nodes are required be constructed on the boundary, and absolutely no cells are needed either for interpolation purpose or for integral purpose. Moreover, the primitive variables (i.e. velocity and pressure), not the derived variables (e.g. vorticity and stream function), are employed in this approach, so the problem of finding velocity can be separated from that of finding pressure, and the present method can also be extended to 3D problems. The demonstrated good performance of the HBNM shows that it is a powerful tool for the numerical solution of incompressible viscous fluid flows. Although 2D and 3D steady state Stokes flow problems are primarily solved, this method can be extended to nonlinear or unsteady state problems by first reducing them to linear and stationary ones by perturbation and time-stepping procedures, and then solving with the help of the dual reciprocity method (DRM).

Acknowledgments This work is supported by the “Hundred Talents Program” of the Chinese Academy of Sciences and the National Natural Science Foundation of China (No. 11272330). The financial supports are gratefully acknowledged. The authors also would like to thank anonymous reviewers for their useful comments and language editing which have greatly improved the manuscript. References [1] Girault V, Raviart PA. Finite Element Approximation of the Navier-Stokes Equations. Berlin: Springer; 1979. [2] Franca LP, Nesliturk A. On a two-level finite element method for the incompressible Navier-Stokes equations. Int J Numer Meth Eng 2001;52: 433–53. [3] Chorin AJ. A numerical method for solving incompressible viscous flow problems. J Comput Phys 1997;135:118–25. [4] Power H, Wrobel LC. Boundary Integral Methods in Fluid Mechanics. Southampton: Computational Mechanics Publications; 1995. [5] Power H, Mingo R. The DRM subdomain decomposition approach to solve the two-dimensional Navier-Stokes system of equations. Eng Anal Bound Elem 2000;24:107–19. [6] Fang JN, Parriaux A, Rentschler M, et al. Improved SPH methods for simulating free surface flows of viscous fluids. Appl Numer Math 2009;59:251–71. [7] Zhang L, Ouyang J, Zhang XH. On a two-level element-free Galerkin method for incompressible fluid flow. Appl Numer Math 2009;59:1894–904.

[8] Onate E, Sacco C, Idelsohn S. A finite point method for incompressible flow problems. Comput Vis Sci 2000;3:67–75. [9] Fang JN, Parriaux A. A regularized Lagrangian finite point method for the simulation of incompressible viscous flows. J Comput Phys 2008;227: 8894–908. [10] Lin H, Atluri SN. The meshless local Petrov-Galerkin (MLPG) method for solving incompressible Navier-Stokes equations. Comput Model Eng Sci 2001;2(2):117–42. [11] Wu XH, Tao WQ, Shen SP, et al. A stabilized MLPG method for steady state incompressible fluid flow simulation. J Comput Phys 2010;229:8564–77. [12] Sellountos EJ, Sequeira A. An advanced meshless LBIE/RBF method for solving two-dimensional incompressible fluid flow. Comput Mech 2008;41:617–31. [13] Ding H, Shu C, Yeo KS, et al. Numerical computation of three-dimensional incompressible viscous flows in the primitive variable form by local multiquadric differential quadrature method. Comput Meth Appl Mech Eng 2006;195:516–33. [14] Shan YY, Shu C, Lu ZL. Application of local MQ-DQ method to solve 3D incompressible viscous flows with curved boundary. Comput Model Eng Sci 2008;25(2):99–113. [15] Young DL, Jane SC, Lin CY, et al. Solutions of 2D and 3D Stokes laws using multiquadrics method. Eng Anal Bound Elem 2004;28:1233–43. [16] Chinchapatnam PP, Djidjeli K, Nair PB. Radial basis function meshless method for the steady incompressible Navier-Stokes equations. Int J Comput Math 2007;84:1509–26. [17] Demirkaya G, Soh CW, Ilegbusi OJ. Direct solution of Navier-Stokes equations by radial basis functions. Appl Math Model 2008;32:1848–58. [18] Alves CJS, Silvestre AL. Density results using Stokeslets and a method of fundamental solutions for the Stokes equations. Eng Anal Bound Elem 2004;28:1245–52. [19] Young DL, Jane SJ, Fan CM, et al. The method of fundamental solutions for 2D and 3D Stokes problems. J Comput Phys 2006;211:1–8. [20] Young DL, Lin CY, Fan CM, et al. The method of fundamental solutions for solving incompressible Navier-Stokes problems. Eng Anal Bound Elem 2009;33:1031–44. [21] Li XL, Zhu JL. A meshless Galerkin method for Stokes problems using boundary integral equations. Comput Meth Appl Mech Eng 2009;198:2874–85. [22] Li XL. Meshless analysis of two-dimensional Stokes flows with the Galerkin boundary node method. Eng Anal Bound Elem 2010;34:79–91. [23] Li XL. Development of a meshless Galerkin boundary node method for viscous fluid flows. Math Comput Simulat 2011;82:258–80. [24] Li XL. The meshless Galerkin boundary node method for Stokes problems in three dimensions. Int J Numer Meth Eng 2011;88:442–72. [25] Zhang JM, Yao ZH, Li H. A hybrid boundary node method. Int J Numer Meth Eng 2002;53:751–63. [26] Zhang JM, Yao ZH. The regular hybrid boundary node method for 2D linear elasticity. Eng Anal Bound Elem 2004;27:259–68. [27] Miao Y, Wang YH. Development of hybrid boundary node method in two dimensional elasticity. Eng Anal Bound Elem 2005;29:703–12. [28] Miao Y, Wang YH. Meshless analysis for three-dimensional elasticity with singular hybrid boundary node method. Appl Math Mech 2006;27(5):673–81. [29] Miao Y, Wang Q, Liao BH, et al. A dual hybrid boundary node method for 2D elastodynamics problems. Comput Model Eng Sci 2009;53(1):1–22. [30] Yan F, Wang YH, Miao Y, et al. Dual reciprocity hybrid boundary node method for free vibration analysis. J Sound Vib 2009;32:1036–57. [31] Wang Q, Miao Y, Zheng JJ. The hybrid boundary node method accelerated by fast multipole expansion technique for 3D elasticity. Comput Model Eng Sci 2010;70(2):123–51. [32] Wang Q, Miao Y, Zhu HP, et al. An O(N) fast multipole hybrid boundary node method for 3D elasticity. Comput Mater Continua 2012;28(1):1–26. [33] Tsai CC, Young DL. Cheng AHD. Meshless BEM for three-dimensional Stokes flows. Comput Model Eng Sci 2002;3(1):117–28.