Exponential basis functions in solution of problems with fully incompressible materials: A mesh-free method

Journal of Computational Physics 231 (2012) 7255–7273

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Exponential basis functions in solution of problems with fully incompressible materials: A mesh-free method S.M. Zandi a, B. Boroomand a,⇑, S. Soghrati b a b

Department of Civil Engineering, Isfahan University of Technology, Isfahan 84156-83111, Iran Department of Civil and Environmental Engineering, University of Illinois at Urbana-Champaign, 205 North Mathews Avenue, Urbana, IL 61801, USA

a r t i c l e

i n f o

Article history: Received 14 November 2011 Received in revised form 30 May 2012 Accepted 26 June 2012 Available online 11 July 2012 Keywords: Incompressible material Elasticity problem Steady state fluid Meshless method Exponential basis function Fundamental solution Discrete transformation

a b s t r a c t In this paper exponential basis functions (EBFs) satisfying the governing equations of elastic problems with incompressible materials are introduced. Due to similarity between elasticity problems and steady state fluid problems the bases found for the former problems are used for latter problems. We discuss on using single field form known as displacement/velocity based formulation and also on using a two-field form known as u-p formulation. In the first formulation we find the pressure bases through performing a limit analysis using a fictitious bulk modulus while in the second formulation the bases are found directly by considering the pressure as a separate variable. In the second formulation we directly apply the condition of incompressibility. It is shown that both formulations yield identical bases meaning that the first one may be used in a standard approach. However, it is also shown that when the incompressibility condition is applied by a Laplacian of pressure in the second formulation, some additional spurious EBFs may be obtained. Having defined appropriate bases, we follow the solution strategy recently introduced by the authors for other engineering problems. Some well-known benchmark problems are solved to show the capabilities of the method. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Although no truly incompressible material exists, there are some engineering problems in which the deviatoric part of strains (or strain rate) governs the material behaviour. In transient problems this is usually due to the fact that volumetric deformations, carrying small portion of energy, propagate relatively fast and either leave the domain or convert to other forms of energy with less effect on further deformations. Therefore many low-rate fluid flows or the so called ‘‘steady state’’ ones are sometime mathematically modelled by materials with incompressible behaviour. Studies on steady-state fluid flows with low-Reynolds number may be performed by analyzing an equivalent solid problem while replacing displacements with velocity components and using viscosity in place of shear modulus. In such problems, although the volumetric deformation is ignored, the mean stress (or the pressure) exists. This causes a sort of indeterminacy with respect to the mean stress (or the pressure). This is in contrast with the behaviour of compressible materials in which the pressure is related to the volumetric part of the strain via the bulk modulus [1]. In the literature of numerical simulations two main approaches are traceable for applying the constraint; the use of standard displacement/velocity formulation with ‘‘nearly incompressible’’ materials and the use of mixed formulation with ‘‘fully incompressible’’ material behaviour. Nevertheless, in both approaches the imposition of the constraint not only produces

⇑ Corresponding author. E-mail address: [email protected] (B. Boroomand). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.06.036

7256

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

well-known numerical instabilities but also leads to poor pressure values. Nowadays, it is well understood that the difficulty arises from the approximation used. This is the main issue that we focus on here in this paper. Focusing on problems with constant material properties, e.g. Newtonian fluid flows, the approximation we use in this paper is based on the construction of a set of residual-free exponential basis functions (EBFs). Here the term ‘‘residual-free’’ refers to the fact that the bases satisfy the governing equations, e.g. Navier–Stokes equations, and at the same time the incompressibility condition. With the use of such basis functions, the proposed method falls in the category of Trefftz methods (see [2] for a good review). On the same line are the methods using fundamental solutions, like boundary element method (BEM) and the method of fundamental solutions (MFS). In MFS the boundaries are discretized into points and the boundary conditions are satisfied through a collocation approach while the sources of the fundamental solutions are placed outside of the solution domain [3]. In the context of using MFS in incompressible Stokes flows, one may refer to studies in [4– 6]. Latest practical application of the method may be seen in [7]. In this paper, in order to find appropriate EBFs, we employ two approaches mentioned earlier; i.e. a single displacement/ velocity field formulation and a two-field u-p formulation. In the first approach we need explicit relation between the mean stress (or the pressure) and the volumetric strain. This however is not known for incompressible materials. Therefore we start from a fictitious compressible material and, through a limit analysis, show how the pressure bases can be found from displacement ones. To show the validity of the pressure bases found, we repeat the procedure with the second approach. This helps us to demonstrate that, with the aid of appropriate EBFs, the standard equations can be solved even when the bulk modulus tends to infinity (or when Poisson’s ratio tends to 0.5). This however is in contrast with other numerical methods based on displacement formulations (that fail when the material becomes fully incompressible). The immediate conclusion is that, in our approach it is not necessary to use methods with mixed formulations, iterative solutions and so on. Worthwhile to mention that, as we shall show later, for a special form of the second formulation some spurious EBFs may also be obtained. Such an observation reveals the significance of using single field formulation. Having found appropriate EBFs for fully incompressible materials, we solve benchmark problems with general boundary conditions in a meshless style. To this end, we write the approximate solution as a series of EBFs with constant coefficients. The constant coefficients are evaluated through a point collocation method on the domain boundaries via a complex discrete transformation technique. The transformation, proposed in [8–11], allows us to use a larger number of basis functions compared with the number of boundary information. The method has been applied to static and time harmonic elastic problems and described in more details in [11]. For further details the readers may refer to [12–14] for the application of the method to elasticity problems such as plates. Also, in [15] the authors have recently used the method for analyzing incompressible inviscid flow with moving boundaries using Lagrangian formulation. The layout of the paper is as follows. In Section 2 the model problem is described and two formulations are used to obtain the governing equations in three forms. In Section 3 the associated EBFs are given for the cases presented in the preceding section. In Section 4 the procedure of the solution method is described and finally in Section 5 several benchmark problems are solved to show the capability of the proposed procedure in the solution of problems with different domain shapes.

2. Governing equations We start from an elasto-static problem with compressible material and try to find the solution for a problem with fully incompressible material. As mentioned earlier, low-Reynolds steady state stokes flow problems may be considered by small changes/replacements afterward. The governing equations of an elasto-static problem in the absence of body forces may be written in a vector notation as

ST DSu ¼ 0

ð1Þ

which is to be solved with the following boundary conditions:

u ¼ uB

on Cu

ð2Þ

and

~ DSu ¼ tB n

on Ct

ð3Þ

In the above relations, u is the vector of displacements (or velocities for steady stokes flow), S is the well-known operator for defining strains (or strain rate) as e ¼ Su, D is the matrix of material constants (here in 2D elasticity problems), uB and tB are ~ contains the components the prescribed displacements and tractions on the boundaries Cu and Ct , respectively. The matrix n of the unit vector normal to the boundary for defining the tractions. Here we consider 2D isotropic problems whose matrix of material constants D can be presented as

2

D1 6 D ¼ 4 D2 0

D2 D1 0

3 0 7 0 5 D3

ð4Þ

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

7257

where in the case of plane strain isotropic problems we have

D1 ¼

Eð1  mÞ ; ð1 þ mÞð1  2mÞ

D2 ¼

Em ; ð1 þ mÞð1  2mÞ

D3 ¼

E 2ð1 þ mÞ

ð5Þ

in which E and m denote the Young’s modulus and the Poisson’s ratio of the material. It is clear that in the incompressibility limit m ! 0:5, the matrix of material constants, D, becomes indeterminate and the pressure cannot be evaluated from the strain field. Another interpretation is that when m ! 0:5 the bulk tends to infinity;

K¼

E ; 3ð1  2mÞ

lim K ¼1

ð6Þ

m!0:5

Since the mean stress (or the pressure) is to remain finite, an immediate consequence of the above relations is that the volumetric strain should vanish. This can clearly be seen from the following expression which relates the pressure part p to the volumetric strain ev for a standard compressible material,

p ¼ K ev ;

lim ev ¼0

ð7Þ

K!1

where

ev ¼ r  u

ð8Þ

In which r is the well-known gradient operator. Note that in (7) we have considered the pressure as the mean stress or vice versa (disregarding the sign) for convenience. A new interpretation of (7) may be given as follows

limðK ev Þ ¼ p

ð9Þ

m!0:5

Although this is true for exact solutions, in a numerical practice it is not achievable. However we shall use it to find EBFs associated with the pressure part. Eqs. (7) and (9) show that, at incompressibility limit, the direct relation between ev and p is lost which leads to indeterminacy of the problem with respect to the pressure. Therefore it is convenient to separate the pressure from the total stress field and treat it as an independent variable. This leads to methods based on mixed formulation and for obtaining solution one may use iterative procedures [1,16]. In the mixed formulations, by splitting the stress field into deviatoric and volumetric part, the pressure is regarded as an additional unknown and one may rewrite (1) as

ST ðDd Su þ pÞ ¼ 0

ð10Þ

where Dd is the deviatoric part of the matrix of material constants that relates the deviatoric stresses to the total strains [1]. In the above equation, for 2D problems, p is

pT ¼ ½ 1 1 0 p

ð11Þ

In the latter formulation, referred to as u-p formulation, an additional equation, known as incompressibility constraint, should be introduced:

ru¼0

ð12Þ T

By operating divergence operator on the vector Eq. (10), note that in (10) S P  rp, and considering (12) we obtain

r2 p ¼ 0

ð13Þ

2

where r is the Laplacian operator. Considering (12) or (13) results in two forms of mixed formulations which are recommended in many studies. However, as will be seen later, considering (13) leads to some invalid exponential basis which should be removed before using the series. In this paper we shall frequently refer to the following three formulae. Case 1.

(

ST DSu ¼ 0

m ! 0:5

ð14Þ

Case 2.

(

ST Dd Su þ ST p ¼ 0

ru¼0

ð15Þ

7258

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

Case 3.

(

ST Dd Su þ ST p ¼ 0

ð16Þ

r2 p ¼ 0

The boundary conditions in Case 1 are as those given in (2) and (3) but for Case 2 and 3 one may need to specify boundary values for the pressure as well. As is seen, Case 1 represents single field formulation for which many numerical methods fail in yielding solution. The reader may note that in fluid problems, with low Reynolds number, with viscosity l in hand one may choose a fictitious Poisson’s ratio to define a fictitious Young’s modulus as E ¼ 2ð1 þ mÞl to use (4) and (5) in (14)–(16). In that case the vector u in (1) represents the fluid velocity field. 3. Exponential basis functions In this section we shall find some EBFs that exactly satisfy the governing equations. The main idea, here, is to find some basis functions so that when m ! 0:5 the pressure part in (9) gives finite values. However to demonstrate that no other basis exists, we shall also use u-p formulations to find the EBFs. 3.1. EBFs for Case 1 (single field formulation) The displacement/velocity EBFs for Case 1 may be found by assuming (see also [11])

u ¼ heaxþby

ð17Þ

where x and y are the coordinates of a generic point in the problem domain, a and b are two complex values. Also T h ¼ h h1 h2 i is an appropriate vector to be determined. By inserting (17) in (14), the following relation is resulted

"

D1 a2 þ D3 b2

ðD2 þ D3 Þab

ðD2 þ D3 Þab

2

#

D3 a2 þ D1 b

h1 h2

 ¼0

ð18Þ

The non-trivial solution is obtained when the determinant of the coefficient matrix vanishes, which gives

ða2 þ b2 Þ2 ¼ 0

ð19Þ

From the above characteristic equation one may find a in terms of b or vice versa. When a is calculated in terms of b, we obtain

a ¼ ib ðFolded rootsÞ; h ¼ h i 1 iT

ð20Þ

pffiffiffiffiffiffiffi where i ¼ 1. When b is calculated in terms of a T

b ¼ ia ðFolded rootsÞ; h ¼ h i 1 i

ð21Þ

Since the characteristic Eq. (19) has folded roots, the missing EBFs are found by considering another set of solutions as follows

u ¼ ðax þ by þ cÞeaxþby

ð22Þ

where a; b and c are new arrays to be found. Substitution of (22) in (14) results in the following system of equations

Aax þ Aby þ ðAc þ Ca þ BbÞ ¼ 0

ð23Þ

where A is the same coefficient matrix as in (18) and

B¼



2D3 b

ðD2 þ D3 Þa

ðD2 þ D3 Þa

2D1 b

 ;

C¼



2D1 a

ðD2 þ D3 Þb

ðD2 þ D3 Þb

2D3 a

 ð24Þ

Since the determinant of A vanishes, one may consider

a ¼ ib ! a ¼ ah i 1 iT ; b ¼ bh i 1 iT

ð25-aÞ

and T

b ¼ ia ! a ¼ ah i 1 i ;

T

ð25-bÞ

b ¼ bh i 1 i

where a and b are two arbitrary coefficients. By substituting (25) in (23), considering definitions (5), the third array is evaluated

(

a ¼ ib;

For a – 0; b ¼ 0 ! c ¼ ahð4m  3Þ=b 0iT T

For a ¼ 0; b – 0 ! c ¼ bhið3  4mÞ=b 0i

ð26-aÞ

7259

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

and

( b ¼ ia;

For a – 0; b ¼ 0 ! c ¼ ahið3  4mÞ=a 0i

T

ð26-bÞ

For a ¼ 0; b – 0 ! c ¼ bhð4m  3Þ=a 0iT

The final solution thus may be written as (for instance a ¼ ibÞ

u ¼

( )! ( )! "         ð4m3Þ ð4m3Þ i i ib xþb y i i b b eibi xþbi y þ c2i e i i þ c3i xþ xþ c1i eibi xþbi y þ c4i eibi xþbi y 1 1 1 1 0 0 ( )! ( )! #     mÞ mÞ i i i ð34 i ð34 b b yþ yþ eibi xþbi y þ c6i eibi xþbi y þc5i 1 1 0 0

P

ð27Þ

Remark 1. The bases shown in (27) may be combined to obtain new forms of bases which of course are not independent of the original set. For instance one may write

( )!     ð4m3Þ i i ibi xþbi y b eibi xþbi y ¼ k e þ xþ 1 1 0

( )!   0 i eibi xþbi y ; xþ i ð4mb3Þ 1

k¼i

ð4m  3Þ b

ð28Þ

This is equivalent to choosing an appropriate ratio between a and b in (26). This shows that the following expression may also be used in place of (27)

" ( ) ( )! ( )! ( ) ( ) ( ) 0 0 i i i i 1 ibi xþbi y 2 ibi xþbi y 3 ibi xþbi y 4 e eibi xþbi y e e xþ xþ u¼ ci þ ci þ ci þ ci i ð4mb3Þ i ð4mb3Þ 1 1 1 1 ð29Þ ( )! ( )! # ( ) ( ) 0 0 i i eibi xþbi y þ c6i eibi xþbi y yþ yþ þc5i þ ð4mb3Þ  ð4mb3Þ 1 1 X

Other forms are still possible but the important effect in all possible forms is that the six bases should be linearly independent. The proof of linear independency may be performed through conventional mathematical methods, i.e. using the Wronskian determinant [17]. However, unlike the conventional approaches, here we are dealing with some vector functions. Fortunately in Case 1, we have six bases each having two elements. This helps us to use the following expressions to perform the study

ui ¼ 0;

@ ui ¼ 0; @x

@ ui ¼ 0 @y

ð30Þ

with ui containing a set of six bases, in (27) or (29) or any other set, associate with a given bi . Arranging (30) in a matrix form, one may write

Hi ci ¼ 0;

 ci ¼ ci1

ci2

ci3

ci4

ci5

ci6

T

ð31Þ

Here Hi is 6  6 matrix obtained from (30). For the case of detðHi Þ – 0 the independency of the vector functions is proved. Such a procedure is equivalent to defining a rectangular matrix H1i so that ui ¼ H1i ci and searching for ci – 0 when ui ¼ 0. This may be performed by arranging the matrix in a square matrix, with other elements being zero, and searching for the null space of the new square matrix. Upon obtaining null spaces with non-variable elements, the dependency of the functions is proved. This latter form may be used in a similar study in Cases 2 and 3. h Remark 2. In this study we have used a combination of EBFs to obtain four independent ones. To this end, we have considered a ¼ b ¼ 1 in (25-a) and obtained the following expression for u

8 91 2 ( ) 0( ) ( ) ( ) < iðð3þ3iÞþbi ð4þ4iÞmÞ = þi i þi þi b i 4c1 Aeibi xþbi y u¼ eibi xþbi y þ c2i eibi xþbi y þ c3i @ xþ yþ i : ; 1 1 1 1 1 91 8 0( ) 3 ( ) < iðð33iÞþbi ð44iÞmÞ = i i bi 4@ ib xþb y Ae i i 5 xþ yþ þci ; : 1 1 1 X

ð32Þ

Our experience shows that the above expression (considering similar relation for b ¼ iaÞ is sufficient for solving many problems in solid and fluid mechanics. h It is interesting to note that the first two EBFs in (27) (or (29)) are divergence-free and this means that they satisfy incompressibility condition. However, the rest of EBFs are not divergence-free which means that these two sets are responsible for volumetric strain as well as deviatoric strains. Another significant property in the above EBFs is that they are valid for all values of Poisson’s ratio including m ¼ 0:5.

7260

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

3.1.1. Pressure EBFs for single field formulation (Case 1) With the EBFs of displacements/velocities (27) in hand, one can evaluate the pressure basis functions straightforwardly by considering (9) as

E

p ¼ lim

m!0:5 3ð1  2mÞ

ru

ð33-aÞ

Note that in (33-a) both numerator and denominator tend to zero as m ! 0:5. The limit may be evaluated conveniently by the L’Hospital rule for each displacement EBF as

p ¼ lim

E r  ð@u=@ mÞ E @u ¼  lim r  @ð1  2mÞ=@ m 6 m!0:5 @m

ð33-bÞ

m!0:5 3

The displacements/velocities EBFs and the associated pressure basis functions for Case 1 are summarized in Table 1. The table inherently shows that the displacement/velocity field given in (27) produces the following pressure field (a ¼ ibÞ

p ¼ 2D3

X

c3i ieibi xþbi y þ c4i ieibi xþbi y  c5i eibxþby  c6i eibxþby



ð34Þ

Remark 3. Upon using EBFs given in (32), the following pressure bases are obtained

p ¼ 2D3

X

c3i ð1  iÞeibi xþbi y þ c4i ð1 þ iÞeibxþby



ð35Þ

We shall use these bases in the solution of numerical examples. h 3.2. EBFs for two-field u-p formulation (Case 2) The EBFs for Case 2 may be found by assuming

  u p

 axþby ; ¼ he

 ¼ hh h 1

h2

h3 i

T

ð36Þ

By inserting (36) in (15) and arranging the system of equations in a matrix form, the following equation is resulted

2

38

9

D3 ð4a2 þ 3b2 Þ=3 D3 ab=3

a > < h1 > =

a

0

6 4 D3 ab=3

7 D3 ð3a2 þ 4b2 Þ=3 b 5

b

> :

h2 h3

> ;

¼0

ð37Þ

By setting the determinant of the coefficient matrix to zero we have

ða2 þ b2 Þ2 ¼ 0

ð38Þ

When a is calculated in terms of b, we obtain

a ¼ ib ðFolded rootsÞ; h ¼ h i 1 0 iT

ð39Þ

and when b is calculated in terms of a

 ¼ h i 1 0 iT b ¼ ia ðFolded rootsÞ; h

ð40Þ

Table 1 Exponential basis functions for standard formulation (Case 1). Displacement/velocity EBFs for all Poisson’s ratios

Associated pressure EBFs for

a ¼ ib 

 i ibxþby e 1 8 91 0   < ð3  4mÞ = @ i x þ  Aeibxþby b 1 : ; 0 8 91 0   < ð3  4mÞ = @ i y þ i Aeibxþby b 1 : ; 0

b ¼ ia   i axiay e 1 ( )!   ð3  4mÞ i eaxiay x þ i a 1 0 ( )!   ð3  4mÞ i eaxiay yþ  a 1 0

0 2iD3 eibxþby

2D3 eibxþby

0 2iD3 eaxiay 2D3 eaxiay

m ¼ 0:5

7261

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

Since the characteristic Eq. (38) has folded roots, the missing EBFs are found by considering another set of solutions for the problem, that is

  u p

¼ ðax þ by þ cÞeaxþby

ð41Þ

Here again a; b and c are new arrays to be found. Substituting (41) in (15) results in the following system of equations

Aax þ Aby þ ðAc þ Ca þ BbÞ ¼ 0

ð42Þ

and again A is the same coefficient matrix as in (37) and

2

3 2D3 b D3 a=3 0 6 7 B ¼ 4 D3 a=3 8D3 b=3 1 5; 0 1 0

2

3 8D3 a=3 D3 b=3 1 6 7 C ¼ 4 D3 b=3 2D3 a 0 5 1 0 0

ð43Þ

Since the determinant of A vanishes, according to (39) and (40) one may consider

a ¼ ib ! a ¼ ah i 1 0 iT ; b ¼ bh i 1 0 iT

ð44-aÞ

and T

b ¼ ia ! a ¼ ah i 1 0 i ; b ¼ bh i 1 0 i

T

ð44-bÞ

Here again a and b are two arbitrary coefficients (here we consider a ¼ b ¼ 1Þ. By inserting (44) in (42), considering definitions (5), the third array is evaluated

a ¼ ib ! c ¼ c and

b ¼ ia ! c ¼ c

D

D

iðð1iÞþbÞ b

iðð1iÞþaÞ

a

1 2D3 ð1  iÞ

ET

1 2D3 ð1  iÞ

ð45-aÞ

ET

ð45-bÞ

where c in the above relations is also an arbitrary coefficient. The final solution thus may be written as (for a ¼ ibÞ

8 9 8 9 8 iðð1þiÞþb Þ 91 08 9 2 8 9 i > >i > > > > > = < i > = = < = bi P 6 1 < = ib xþb y B C ibi xþbi y eibi xþbi y þ c3i @ 1 x þ 1 y þ 1 u ¼ 4ci 1 e i i þ c2i 1 Ae > > > > > > > > > > : ; : ; : ; : ; : ; 0 0 0 0 2D3 ð1  iÞ 8 iðð1iÞþb Þ 91 8 9 3 08 9 i > > > > < i > < = = < i > = bi B C ibi xþbi y 7 þc4i @ 1 xþ 1 yþ 1 5 Ae > > > > : > : ; ; : > ; 0 0 2D3 ð1 þ iÞ

ð46Þ

As is seen, the pressure components are identical to those obtained earlier in the single field formulation while considering m ¼ 0:5 (see (35)). This shows that there is no missing EBF for the pressure in Case 1. 3.3. EBFs for two-field u-p formulation (Case 3) The EBFs for Case 3 may be found again by assuming (36) and inserting in (16) to obtain the following equation

2 6 4

D3 ð4a2 þ 3b2 Þ=3 D3 ab=3 0

38 9 > < h1 > = 7 D3 ð3a2 þ 4b2 Þ=3 b 5 h2 ¼ 0 > : > ; h3 a2 þ b2 0 D3 ab=3

a

ð47Þ

By setting the determinant of the coefficient matrix to zero we have

ða2 þ b2 Þ3 ¼ 0

ð48Þ

When a is calculated in terms of b, we obtain

D

a ¼ ib ðThree repeating rootsÞ; h ¼ h i 1 0 iT ; h ¼ i 1

2D3 b 3

ET

ð49Þ

and when b is calculated in terms of a

 ¼ h i 1 0 iT ; b ¼ ia ðThree repeating rootsÞ; h

¼ h

D

1 i

2D3 a 3

ET

ð50Þ

It can easily be shown that for repeating roots one may evaluate characteristic vectors as given in (44) and (45), not shown again for the sake of brevity. The final solution thus may be written as (for a ¼ ibÞ

7262

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

8 9 2 8 9 8 9 8 9 >i > > < i = < i = < i > = P 6 1 < = ib xþb y u ¼ 4ci 1 e i i þ c2i 1 eibi xþbi y þ c4i 1 eibi xþbi y þ eibi xþbi y þ c3i 1 : 2D3 bi ; : 2D3 bi ; > > : > : > ; ; 0 0 3 3 8 iðð1þiÞþb Þ 8 iðð1iÞþb Þ 8 9 8 9 91 91 3 08 9 08 9 i i i i i i > > > > > > > > < > < > < < = < > = < > = = = = bi bi 7 B C B C þc5i @ 1 x þ 1 y þ 1 xþ 1 yþ 1 Aeibi xþbi y þ c6i @ 1 Aeibi xþbi y 5 > > > > > > > > > > > > : ; : ; : : : ; : ; ; ; 0 0 0 0 2D3 ð1  iÞ 2D3 ð1 þ iÞ

ð51Þ

However, it can be seen that the third and fourth EBFs in (51) do not satisfy the incompressibility condition, for instance the third EBF gives

r  u ¼ 2c3i bi eibi xþbi y

ð52Þ

which does not vanish except for bi ¼ 0. Such spurious modes should be eliminated from the solution. This shows that for Case 3, there is a possibility of obtaining wrong solution in a numerical method. It may be noted that the Laplace equation for the pressure in Case 3 appears because of assuming full incompressibility. However for time marching algorithms using well-known numerical methods such as finite difference, the relation appears as a Poisson’s equation, i.e. with a right hand side depending on the residual of the incompressibility conditions. In that case the solution may be split into two parts, i.e. homogenous and particular parts, and thus a similar conclusion made here is valid for the homogeneous part. The discussion given here shows the necessity of employing appropriate treatment when using such numerical methods (see [18,19] for instance). 4. The solution procedure Having defined the EBFs, we now use them in a boundary point method. The solution series of Case 1 is written in a more general form as

u¼

N X cj hj eaj xþbj y

ð53Þ

j¼1

Here the counter j is used to denote that all EBFs obtained for a ¼ ib and b ¼ ia are included in the series. Also N is the total number of selected EBFs to achieve adequate accuracy. 4.1. Imposition of boundary conditions In this section, the way that the boundary conditions are imposed through a collocation approach is described [11]. To this end, by considering u in its general form (53) and applying a collocation approach at M boundary points one may write the following relation

B ¼ U

N X C j Vj

ð54Þ

j¼1

 B is a vector containing all point-wise boundary conditions, Vj is a normalized vector which contains In the above relation, U  B , and C j is proportional to the contribution of each exponential basis to the boundaries arranged in the same manner as U coefficient cj in (53) considering the normalization factor used in Vj . For instance supposing that m boundary points, out of M, are allocated to Dirichlet type of conditions and n ¼ M  m points are allocated to Neumann conditions (note that  B can be arranged in the following form: at the corners we choose two close points to define distinct normal vectors), then U

B ¼ U



ðuTB Þ1

ðuTB Þ2

. . . ðuTB Þm

j ðtT Þ1

ðtT Þ2

. . . ðtT Þn

T

ð55Þ

In a similar manner Vj can be arranged as

Vj ¼

1n T ðuj Þ1 sj

ðuTj Þ2

. . . ðuTj Þm

j ðtTj Þ1

ðtTj Þ2

. . . ðtTj Þn

oT

ð56Þ

where sj is a scaling factor for normalization and



ðuj Þk ¼ hj eaj xþbj y x¼x

k ;y¼yk

8ðxk ; yk Þ 2 Cu ; k ¼ 1; . . . ; m

ð57Þ

and



~ k ðDd Shj eaj xþbj y þ Pj Þ ðtj Þk ¼ n x¼x

k ;y¼yk

8ðxk ; yk Þ 2 Ct ; k ¼ 1; . . . ; n

ð58Þ

where PTj ¼ ½1 1 0Pj and Pj denotes the pressure EBFs corresponding to jth displacement basis function (see Table 1 or rela~ k is the outward normal at kth boundary point. The scaling factor sj is used as follows tion (35) ). In (58) n

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

sj ¼ maxðjV lj jÞ;

l ¼ 1; . . . ; M

l

7263

ð59Þ

with V lj being the lth element of Vj (in the above relations jj denotes the Hermitian length). The readers may refer to [11] for normalization with respect to material constants appearing in part of Vj as ðtj Þk . Note that with defining scaling factor sj , the coefficients cj in (53) are now related to the coefficients C j in (54) as

cj ¼

1 Cj sj

ð60Þ

Now for evaluating the coefficients C j we assume

B C j ¼ VTj RU

ð61Þ

where R plays the role of a projection matrix assumed to be suitable for all j ¼ 1; . . . ; N. By inserting (61) in (54) we find

B ¼ U

N X

 B ¼ GRU  B; Vj VTj RU

G¼

j¼1

N X Vj VTj

ð62Þ

j¼1

In the above relations, G is a symmetric M  M matrix. Since the rank of G might be less than M, R is evaluated as

R ¼ Gþ

ð63Þ

þ

where G is pseudo inverse of G. Now we can evaluate the displacements (or velocities) by the following equation

( u¼R

!

N X 1 j¼1

)

B eaj xþbj y hj VTj RU

sj

ð64Þ

In the above equation, RðÞ denotes the real part of the quantity. For more studies on the applications of the transformation used above the reader may refer to [8–11]. Note that the pressure does not appear as a separate variable during the solution process. Evaluation of pressure and stresses is a post processing task that will be discussed in Section 4.3. 4.2. Selection of a and b The parameters a and b play an important role in the variation of EBFs throughout the solution domain and thus have a prominent effect on the projection matrix R. Since these parameters can take on complex values, we define a grid of points on Gaussian plane. The solution accuracy may differ by changing the grid. A detailed discussion on the selection of a and b is out of the scope of this paper and the reader may refer to [11] for further information. In this reference, the authors have suggested two simple strategies in this regard; one with mathematical basis and another heuristically based on numerical experiments. As shown in [11], the obtained accuracy with either of these approaches is satisfactory in all the cases studied. In this paper we have used the heuristic strategy (see Appendix A). 4.3. Pressure and stress evaluation As mentioned before, pressure and stress evaluation is a post processing task in the present method. After finding the projection matrix R and subsequently the coefficients C j in (61), by the use of the pressure bases introduced in Section 3.1, one may write

p¼

N X cj P j

ð65Þ

j¼1

where Pj denotes the pressure EBFs corresponding to jth displacement basis function (see Table 1 or relation (35)). Eq. (65) can also be written as

( p¼R

N X 1 j¼1

sj

! Pj VTj

) B RU

ð66Þ

For evaluating the stress vector r, we can write

r ¼ Dd Su þ p

ð67Þ

In the above relation, p is a vector defined in (9) considering (66) for p and u is directly obtained from (64). 5. Numerical results In this section we present the results of the method applied to some benchmark problems. These numerical examples include problems in solid mechanics and also some well-known fluid flow benchmarks. Numerical results are compared

7264

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

with those obtained from analytical methods and other numerical approaches to show the capabilities of the proposed procedure.

5.1. Thick-walled long tube subjected to internal pressure As the first numerical experiment, we consider a long tube with internal pressure that may be analyzed as a plane strain case. The inner radius is r i ¼ 0:5, the outer one is ro ¼ 1:0 and tube is subjected to an internal pressure pi ¼ 1:0. Material of the tube is incompressible linear elastic with shear modulus G ¼ 1:0. The analytical solution for the radial displacement field, radial and hoop stress is [20]

ur ¼

1 pi r 2i r 2o ; 2G r2o  r2i r

rrr ¼

  pi r2i r2o ; 1  r2 r 2o  r 2i

rhh ¼

  pi r 2i r2o 1 þ r2 r 2o  r 2i

ð68Þ

For solving this pressure dominant problem we use 146 points on the halved-domain boundaries (Fig. 1). The number of basis functions is 448. Fig. 2 shows the analytical and numerical result of the radial displacement. The radial and hoop stresses are compared with the analytical solution in Fig. 3. As is seen, the numerical results are in excellent agreement with the analytical solution.

5.2. Cook’s membrane problem The Cook’s membrane problem is a plane strain tapered panel clamped at one end and subjected to an in-plane shear load on the free end (Fig. 4). The material is fully incompressible (m ¼ 0:5Þ and the Young’s modulus is E ¼ 250. We consider a constant shear stress s ¼ 100=16 producing total shear force of V =100 at the free end. For numerical simulation, 181 boundary points and 448 EBFs are used. The vertical displacement at point A is obtained as 7.59023 which is in good agreement with other studies [21]. Fig. 5 shows the pressure variation along line BC that is consistent with the results obtained in [21].

Fig. 1. Half domain and the boundary points used.

Fig. 2. Radial displacement versus radius of the tube.

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

7265

Fig. 3. Stresses in tube; (a) radial stress, (b) hoop stress.

Fig. 4. Cook’s problem definition (left) and boundary points (right).

5.3. Lid driven square cavity We choose the well-known lid driven cavity problem as the next numerical experiment [1,4,22]. Here a square domain is considered with unit side lengths (Fig. 6). The fluid/material is assumed to be fully incompressible with unit viscosity (unit elastic shear modulus). All boundaries are restrained in x and y directions, but at the top boundary a unit tangential velocity (displacement) is applied. In this problem we use 800 basis functions and 200 boundary points. The results are presented for the horizontal velocity along vertical centre line AA and for vertical velocity and pressure along horizontal centre line BB (Figs. 7 and 8). We have included the results obtained in reference [22] for comparison. It can be seen that not only the velocities obtained are of high accuracy but also the pressure field is free from oscillation.

7266

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

Fig. 5. Cook’s problem: pressure along the line BC.

Fig. 6. Lid driven cavity problem.

5.4. Rectangular cavity with wave-shaped bottom The benchmark is of the Stokes flow in a rectangular cavity with an irregular boundary [5,22]. Fig. 9 depicts the schematic diagram as well as the boundary points of a wave-shaped bottom cavity with top lid moving with a unit velocity in the xdirection. Similar to the previous problem, the fluid/material is fully incompressible with unit viscosity (unit elastic shear modulus). Also, all boundaries are restrained in x and y directions, but the top boundary experiences a unit tangential velocity (displacement). In this problem we use 560 basis functions and 148 boundary points. Fig. 10 shows the comparison of the results for u1 –y plot obtained in [5] and the presented meshless method in this paper. Fig. 11 depicts the velocity vectors and the pressure contours. The conclusions made in the previous examples seem to be valid in this problem.

5.5. Circular cavity We consider a problem of fluid flow in a circular cavity (R = 1) for which analytical solution has been given in [23]. The upper half surface of the circular cavity moves with a constant unit circumferential velocity in the counterclockwise direction and no-slip boundary conditions are imposed on the lower half surface of the circular cavity. The cavity is filled with

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

7267

Fig. 7. Comparison of velocities (displacements) along centre lines for lid driven square cavity problem; (a) horizontal velocity along AA, (b) vertical velocity along BB.

Fig. 8. Pressure along the horizontal centre line BB for square cavity.

7268

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

Fig. 9. Rectangular cavity with wave-shaped bottom; (a) problem geometry, (b) boundary points M ¼ 148.

Fig. 10. Comparison of horizontal velocity profiles along x = 2.5 for a wave-shaped bottom cavity.

fully incompressible fluid (m ¼ 0:5Þ with unit viscosity. Fig. 12 depicts the geometric configuration and the boundary conditions. The numerical results are compared with the analytical solutions reported by Hwu et al. [23]. For solving this problem we use 640 basis functions and 60 boundary points. Fig. 13-a shows the horizontal velocity profiles along vertical centre line obtained by the presented method in this paper and the analytical solution in [23]. Fig. 13-b depicts pressure distribution along horizontal centre line of the circular cavity (not reported in the literature). The velocity vectors are shown in Fig. 14. As expected, excellent agreement is seen between the two sets of results and also there is no non-physical oscillation in the pressure values. 5.6. Circular cavity with eccentric rotating cylinder The presented model is employed to solve a 2D Stokes flow confined in the two eccentric cylinders [5]. Fig. 15 illustrates the problem including the geometric configuration and the mixed boundary conditions (mixed Neumann and Dirichlet ones) associated with a half of the domain.

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

7269

Fig. 11. Velocity vectors and pressure contours for a wave-shaped bottom cavity; (a) velocity vectors, (b) pressure contours.

Fig. 12. Incompressible fluid flow in a circular cavity (R = 1).

The outer cylinder with radius R1 ¼ 1:0 is centred at the origin, and the inner cylinder with R2 ¼ 0:5 is centred at ð0:25; 0Þ. The outer cylinder rotates with a constant unit circumferential velocity in the clockwise direction and the inner cylinder is kept stationary. No-slip boundary conditions are specified on the fluid–solid interfaces. The velocity vector plots obtained by 128 boundary points and 560 basis functions are illustrated in Fig. 16. For the verification purpose, the computed velocities at different points along y ¼ 0 are shown in Fig. 17. Table 2 shows the net flow through the horizontal cross section, y ¼ 0, calculated from the results of the analytical solution, proposed method, MFS and FEM [5]. The numerical results clearly reveal the capability of the present method in simulation of Stokes flow.

6. Conclusions In this paper, exponential basis functions (EBFs) have been used for the modeling of problems with fully incompressible elastic materials. The EBFs presented here eliminate the indeterminacy effect seen in the standard displacement elasticity formulation when Poisson’s ratio becomes 0.5. Due to the resemblance between the governing equations of

7270

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

Fig. 13. Comparison of the results in a circular cavity; (a) horizontal velocity along vertical centre line, (b) pressure along horizontal centre line.

Fig. 14. Velocity vectors in a circular cavity.

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

Fig. 15. Circular cavity with eccentric rotating cylinder; (a) problem definition, (b) half domain boundary conditions.

Fig. 16. Velocity vectors in an eccentric rotating cylinder.

7271

7272

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

Fig. 17. Comparison of u2 -velocity profiles along y = 0 for the eccentric cylinders.

Table 2 Comparison of the net flow in example 5.6.

R1

0:25 u2 ðx; 0Þdx R 1 Net flow 0:75 u2 ðx; 0Þdx

Net flow

Analytical [24]

MFS 270 nodes [5]

FEM 4000 nodes [5]

Present method

0.17671

0.17688

0.17659

-0.17670

0.17671

0.17706

0.17637

0.17669

low Reynolds fluid flows and elasto-static incompressible problems, one may use the formulation to solve steady state incompressible fluid flow problems. The EBFs are constructed to solve standard elasticity/fluid equations without the necessity of using mixed methods, iterative solutions and so on. In this numerical scheme, pressure can be directly evaluated as a secondary variable after solving the problem by considering the displacements (or velocities) as primary variables. The boundary conditions are imposed through a collocation approach and thus the method can be categorized in meshless types. The proposed numerical scheme is applied to six examples: (1) Thick-walled long tube subjected to internal pressure, (2) Cook’s membrane, (3) a square cavity, (4) a wave-shaped bottom cavity, (5) a circular cavity, and (6) a circular cavity with eccentric rotating cylinder. Compared to the analytical and other available numerical solutions, our numerical experiments demonstrate that the present method is capable of producing very accurate results by using a few collocation points. The results promisingly show the capabilities of the proposed meshless method for being used as a unified formulation for solid and fluid phases in problem with fluid–structure interaction. It is however clear that the method presented in this paper, in its current form, is suitable for engineering problems with homogeneous materials. The reason lies in fact that we assume the governing equations are of constant coefficients type. Our latest studies, to be presented later, pertain to the development of a similar method with alternative bases functions for problems with inhomogeneous materials. Noteworthy also is that the method presented in this paper can be straightforwardly extended to three dimensional problems. Considering time as a new axis, the logic for extending the method to three dimensional problems may be reused to solve time-dependent problems while considering the partial differential equations in space and time. These issues are the subjects of current studies by the authors. Another relevant area of research pertains to the application of the method to problems in exterior domains for which the bases must be chosen so that the radiation conditions are satisfied (see [10] for similar discussions and studies). Appendix A The heuristic strategy for choosing a and b proposed in [11] is as follows

!   k c m a ¼    þ 2i ; L N

  ¼ 1; . . . ; M; m

 ¼ 1; . . . ; N  k

for b ¼ ia. In a similar manner, when a ¼ ib we select

ðA-1Þ

S.M. Zandi et al. / Journal of Computational Physics 231 (2012) 7255–7273

!   k c m b ¼    þ 2i ; L N

  ¼ 1; . . . ; M; m

 ¼ 1; . . . ; N  k

7273

ðA-2Þ

 NÞ  2 N2 ; c  2 R and L is a characteristic length. The following bounds are found to be In the above relations, we choose ðM; appropriate for many cases

 6 7:2 ðsay c  ¼ 2pÞ; 5:6 6 c

L ¼ 1:6 maxðLx ; Ly Þ;

 min ¼ 4; M

 min ¼ 2; N

 max ¼ 8 N

ðA-3Þ

where Lx and Ly are the dimensions of the rectangle that circumscribes the domain. According to this pattern, the number of   N.  For modelling of the problems presented in this paper we have used c  ¼ 4  5 and  ¼ 2p; N constructed EBFs is 16M  ¼ 7  10. M References [1] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, fifth ed., vol. 3, Butterworth, Heineman, 2000. [2] Z.C. Li, T.T. Lu, H.T. Huang, A.H.D. Cheng, Trefftz, collocation, and other boundary methods – a comparison, Numerical Methods for Partial Differential Equations 23 (2007) 93–144. [3] G. Fairweather, A. Karageorghis, The method of fundamental solutions for elliptic boundary value problems, Advances in Computational Mathematics 9 (1998) 69–95. [4] D.L. Young, S.J. Jane, C.M. Fan, K. Murugesan, C.C. Tsai, The method of fundamental solutions for 2D and 3D Stokes problems, Journal of Computational Physics 211 (2006) 1–8. [5] D.L. Young, C.L. Chiu, C.M. Fan, C.C. Tsai, Y.C. Lin, Method of fundamental solutions for multidimensional Stokes equations by the dual-potential formulation, European Journal of Mechanics B/Fluids 25 (2006) 877–893. [6] D.L. Young, Y.C. Lin, C.M. Fan, The method of fundamental solutions for solving incompressible Navier–Stokes problems, Engineering Analysis with Boundary Elements 33 (2009) 1031–1044. [7] F. Boselli, D. Obrist, L. Kleiser, Numerical simulation of the flow in semicircular canals with the method of fundamental solutions, PAMM, Proceedings in Applied Mathematics and Mechanics 9 (2009) 485–486. [8] B. Boroomand, F. Mossaiby, Generalization of robustness test procedure for error estimators part I: formulation for patches near kinked boundaries, International Journal for Numerical Methods in Engineering 64 (2005) 427–460. [9] B. Boroomand, F. Mossaiby, Generalization of robustness test procedure for error estimators part II: test results for error estimators using SPR and REP, International Journal for Numerical Methods in Engineering 64 (2005) 461–502. [10] B. Boroomand, F. Mossaiby, Dynamic solution of unbounded domains using finite element method: discrete Green’s functions in frequency domain, International Journal for Numerical Methods in Engineering 67 (2006) 1491–1530. [11] B. Boroomand, S. Soghrati, B. Movahedian, Exponential basis functions in solution of static and time harmonic elastic problems in a meshless style, International Journal for Numerical Methods in Engineering 81 (2010) 971–1018. [12] B. Shamsaei, B. Boroomand, Exponential basis functions in solution of laminated structures, Composite Structures 93 (2011) 2010–2019. [13] M. Shahbazi, B. Boroomand, S. Soghrati, A mesh-free method using exponential basis functions for laminates modeled by CLPT, FSDT and TSDT; part I: formulation, Composite Structures 93 (2011) 3112–3119. [14] M. Shahbazi, B. Boroomand, S. Soghrati, A mesh-free method using exponential basis functions for laminates modeled by CLPT, FSDT and TSDT; part II: implementation and results, Composite Structures 94 (2011) 84–91. [15] S.M. Zandi, B. Boroomand, S. Soghrati, Exponential basis functions in solution of incompressible fluid problems with moving free surfaces, Journal of Computational Physics 231 (2012) 505–527. [16] I. Bijelonja, I. Demirdzˇic´, S. Muzaferija, A finite volume method for incompressible linear elasticity, Computer Methods in Applied Mechanics and Engineering 195 (2006) 6378–6390. [17] F.B. Hildebrand, Advanced Calculus for Applications, second ed., Prentice-Hall Inc., New Jersey, 1976. [18] J.H. Ferziger, M. Peric´, Computational Methods for Fluid Dynamics, Springer Verlag Berlin Heidelberg, New York, 2002. [19] Y. Morinishi, T.S. Lund, O.V. Vasilyev, P. Moin, Fully conservative higher order finite difference schemes for incompressible flow, Journal of Computational Physics 143 (1998) 90–124. [20] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, McGraw-Hill, New York, 1970. [21] K.B. Nakshatrala, A. Masud, K.D. Hjelmstad, On finite element formulations for nearly incompressible linear elasticity, Computational Mechanics 41 (2008) 547–561. [22] M. Cheng, K.C. Hung, Vortex structure of steady flow in a rectangular cavity, Computers & Fluids 35 (2006) 1046–1062. [23] T.Y. Hwu, D.L. Young, Y.Y. Chen, Chaotic advections for Stokes flow in a circular cavity, Journal of Engineering Mechanics 123 (1997) 774–782. [24] B.Y. Ballal, R.S. Rivlin, Flow of a Newtonian fluid between eccentric rotating cylinders: inertial effects, Archive for Rational Mechanics and Analysis 62 (1976) 237–294.