Elastodynamic problems by meshless local integral method: Analytical formulation

Engineering Analysis with Boundary Elements 37 (2013) 805–811

Contents lists available at SciVerse ScienceDirect

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

Elastodynamic problems by meshless local integral method: Analytical formulation P.H. Wen a, M.H. Aliabadi b,n a b

School of Engineering and Materials Science, Queen Mary, University of London, London E1 4NS, UK Department of Aeronautics, Imperial College, London SW7 2BY, UK

a r t i c l e i n f o

abstract

Article history: Received 22 August 2012 Accepted 30 January 2013 Available online 27 March 2013

In this paper, analytical forms of integrals in the meshless local integral equation method in the Laplace space are derived and implemented for elastodynamic problems. The meshless approximation based on the radial basis function (RBF) is employed for implementation of displacements. A weak form of governing equations with a unit test function is transformed into local integral equations. A completed set of the local boundary integrals are obtained in closed form. As the closed form of the local boundary integrals are obtained, there are no domain or boundary integrals to be calculated numerically. Several examples including dynamic fracture mechanics problems are presented to demonstrate the accuracy of the proposed method in comparison with analytical solutions and the boundary element method. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Meshless local integral equation method Radial basis function Dynamic fracture mechanics

1. Introduction The boundary element method (BEM) is a well-established technique for analysis of certain engineering problem such as fracture mechanics [1–4] and acoustics [5–8]. For certain nonlinear problems such as plasticity [9–11] and creep [12], BEM normally requires some internal discretization. In recent years, the computational mechanics community has turned its attention to so-called mesh reduction methods. These mesh reduction methods (commonly referred to as Meshless or Meshfree) have received much interest since Nayroles et al. [13] proposed the diffuse element method. Later, Belyschko et al. [14] and Liu et al. [15] proposed element-free Galerkin method and reproducing kernel particle methods, respectively. MLPG is reported to provide a rational basis for constructing meshless methods with a greater degree of flexibility. Local Boundary Integral Equation method (LBIE) with moving least square and polynomial radial basis function (RBF) has been developed by Sladek et al. [16]. However, Galerkin-base meshless methods, except MLGP presented by Atluri [17] still include several awkward implementation features such as numerical integrations in the local domain. Other recent developments can be found in Refs. [18–46]. Recently, Sladeks developed an analytical integration method using MLS interpolation (see [40,41]) for non-homogeneous 2D elasticity. However, the analytical solutions are complex. Wen and Aliabadi [42] presented the meshless local integral equation method for the functionally graded Reissner’s plate in analytical formulation.

More recently, Soares et al. [43,44] have presented a time domain method with MLS interpolation for elastodynamic and poroelastic problems. The application of mesh free method to linear elastostatic fracture mechanics, i.e. evaluation of stress intensity factors and analysis of crack growth, was demonstrated by Fleming et al. [45] and Rao et al. [46] using enriched basis function in the moving least square interpolation. However, their method is computationally time consuming because the coefficient matrix must be inverted at each Gauss integration point. To overcome this difficulty, Wen et al. [35] developed a mesh free Galerkin method using enriched radial basis functions. In this paper, a meshless local integral method is presented for two-dimensional dynamic problems. With the use of radial basis functions, analytical solutions for all domain integrals in the weak form are derived. The weak formulation for the governing equations with a unit test function is obtained exactly for the local domain integrals. As the closed form of the local boundary integrals are obtained, the computational time are reduced significantly. A numerical inversion technique, Durbin’s inversion method, is applied to determine each variable in the time domain. The accuracy of the proposed method is illustrated by comparing the numerical results with available analytical solution and boundary element method. 2. The meshless local integral method Consider a linear elastic body in three dimensional domain O with boundary G. The governing equations of motion can be written as

n

Corresponding author. E-mail address: [email protected] (M.H. Aliabadi).

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

sij,j þf i ¼ ru€ i

ð1Þ

806

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 37 (2013) 805–811

where sij are stresses, f i is the body force, r is the density of material and u€ i is the acceleration ði ¼ 1, 2j ¼ 1, 2 for 2DÞ and ði ¼ 1, 2, 3 j ¼ 1, 2, 3 for 3DÞ. Recalling Hooke’s law for plane stress problem, one has     E @u1 @u2 E @u2 @u1 s11 ¼ þ n s ¼ þ n , , 22 @x2 @x1 1n2 @x1 1n2 @x2   E @u1 @u2 s12 ¼ þ ð2Þ 2ð1 þ nÞ @x2 @x1

domain, the transform of a function f ðtÞ is defined as Z 1 f~ ðsÞ ¼ f ðtÞept dt

for two-dimensional problem, where E is Young’s modulus, n is Possion ratio and m ¼ E=2ð1 þ nÞ the shear modulus. The boundary conditions are given as

3. Radial basis function approximation

ui ¼ u0i

on

Gu

t i ¼ sij nj ¼ t 0i

on

Gt

ð3Þ

t 0i

u0i

and are the prescribed displacements and tractions, in which respectively on the displacement boundary GD and on the traction boundary GT , and ni is the unit normal outward to the boundary G. In the local boundary integral equation approach, the weak form of differential equation over a local integral domain Os can be written as Z ðsij,j þ f i ru€ i Þuni dO ¼ 0 ð4Þ Os

where uni is a test function. By use of divergence theorem, Eq. (4) can be rewritten in a symmetric weak form as Z Z sij nj uni dG ðsij uni,j f i uni þ ru€ i ÞdO ¼ 0 ð5Þ Gs

Os

If there is an intersection between the local boundary and the global boundary, a local symmetric weak form in linear elasticity may be written as, see [42] Z Z Z Z Z sij uni,j dO ti uni dG ti uni dG ¼ t0i uni dG þ ðf i ru€ i Þuni dO Ls

Os

GD

GT

Os

ð6Þ in which, Ls is the other part of the local boundary inside the local integral domain Os ; GD is the intersection between the local boundary Gs and the global displacement boundary; GT is a part of the traction boundary as shown in Fig. 1. The local weak forms in (5) and (6) are a starting point to derive local boundary integral equations if appropriate test functions are selected. A step function can be used as the test functions uni in each integral domain ( ji ðxÞ at x ðOs [ Gs Þ : ð7Þ uni ðxÞ ¼ 0 at x= 2Os

where p is a Laplace parameter. Eq. (6) can be rewritten as Z Z s~ ij nj uni dG ðs~ ij uni,j f~ i uni þp2 ru~ i ÞdO ¼ 0 Gs

Γ

node x

Local integral domain Ω s

Ω

Node in support domain ξk

ð9Þ

Os

A local domain Os shown in Fig. 1, which is the neighborhood of a point x ( ¼ fx1 ,x2 g) and is considered as the domain of definition of the RBF approximation for the trail function at x and also called as support domain to an arbitrary point x. To interpolate the function u over a number of randomly distributed nodes in the local support domain Os x½ ¼ fn1 , n2 ,. . ., nK g, nk ¼ ðxk1 , xk2 Þ, k ¼ 1,2,. . .,K, one has function u at the point x as ui ðxÞ ¼

K X

Rk ðx, nÞak þ

T X

P t ðxÞbt ¼ RðxÞa þ PðxÞb

ð10Þ

t¼1

k¼1

where RðxÞ ¼ fR1 ðx, nÞ,R2 ðx, nÞ,. . .,RK ðx, nÞg is a set of radial basis T functions centered around the point x, fak gKk ¼ 1 and fbt gt ¼ 1 are the unknown coefficients to be determined and fPt gTt ¼ 1 is a basis for PT1 , the set of d-variate polynomials of degree rT1. If the radial basis function is selected as multi-quadrics [48,50] qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð11Þ Rk ðx, nk Þ ¼ c2 þ ðx1 xk1 Þ2 þ ðx2 xk2 Þ2 where c is a free parameter and with the constraints and K X

1 rt r T

P t ðnk Þak ¼ 0,

ð12Þ

k¼1

Following polynomials are considered (T ¼ 6) P ¼ f1,x1 ,x2 ,x21 ,x1 x2 ,x22 g

ð13Þ

Then a set of linear equations can be written in a matrix form as R0 a þ P0 b ¼ u, where matrix 2 R ðn Þ

P0 a ¼ 0

ð14Þ 3

R2 ðn1 Þ R2 ðn2 Þ

... ...

:

:

...

: :

: :

... ...

: :

R1 ðnK Þ

R2 ðnK Þ

...

RK ðnK Þ

1

1

6 R1 ðn2 Þ

6 R0 ðnÞ ¼ 6 6 4

RK ðn1 Þ RK ðn2 Þ 7

7 7 7 5

:

ð15Þ KK

and 2P

where fi ðxÞ is arbitrary function. In the Laplace transform support domain of x

ð8Þ

0

3

1 ðn1 Þ

P 2 ðn1 Þ

...

P T ðn1 Þ

6 P1 ðn2 Þ 6 : P0 ðnÞ ¼ 6 6 : 4

P 2 ðn2 Þ

...

: :

... ...

P T ðn2 Þ 7

:

:

...

:

P 1 ðnK Þ

P2 ðnK Þ

...

P T ðnK Þ

7 7 7 5

: :

:

ð16Þ

KT

Solving equations in (14) gives b ¼ bu,

b

a ¼ au

1 ¼ ðPT0 R1 P0 T R1 0 P0 Þ 0 ,

T 1 1 a ¼ R1 P0 T R1 0 ½IP0 ðP0 R 0 P0 Þ 0 

ð17Þ Ls

Γs=ΓD+ΓT Fig. 1. Arbitrary distributed node, support domain of x, local integral domain for weak formulation.

where I denotes the diagonal unit matrix, matrices a, b, a and b are scale ðK  1Þ, ðT  1Þ, ðK  KÞ and ðT  KÞ respectively. Substituting the coefficients a and b from (17) into (10), we obtain the approximation of u in terms of the nodal values: uðxÞ ¼ UðxÞu, UðxÞ ¼ RðxÞa þ PðxÞb:

ð18Þ

also UðxÞ is called as shape function, matrix RðxÞ and PðxÞ are scale

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 37 (2013) 805–811

ð1  KÞ and ð1  6Þ matrix respectively. It is worth noting that the shape function depends uniquely on the distribution of scattered nodes within the support domain and it has the Kronecker Delta property.

807

d2 ¼ ðxlb1 xi1 Þcosbl þ ðxlb2 xi2 Þsinbl þ r 2 If the area of local integral domain is small, the domain integral can be approximated as Ii ¼ p2 ru~ i Os

4. Analytical forms of domain integrals in Laplace space Consider a unit test function, i.e. ji ðxÞ ¼ 1 and the local domain is enclosed by several straight lines as shown in Fig. 2, therefore, the local boundary integral equation becomes: Z Z Z L X s~ ij nj dG ¼ nlj s~ ij dG þ ðf~ i p2 ru~ i ÞdO ¼ 0 ð19Þ Gs

Gl

l¼1

Os

where L is number of straight line. Suppose there are nodes both in the domain and on the boundary, M ¼ M O þ MT þ M D , where M O indicates the number of nodes collocated in domain, M T and MD are numbers of nodes on the traction/displacement boundaries and consider the radial basis function interpolation in (18) and relationship between stress and strain in (2), then (4) becomes, see [42]: "   # K X L K  T  X X X E E l l l l ~ ðjÞ F n þ m F n a þ G n þ m G n ij 2il 2 2tl 2 btj u 1 þ 1n2 1il 1 1n2 1tl 1 t¼1 j¼1l¼1 i¼1 "   # K X L K  T  X X X En En l l l l ~ ðjÞ F n þ m F n a þ G n þ m G n ij 1il 2 1tl 2 btj u 2 ¼ I1 1n2 2il 1 1n2 2tl 1 t¼1 j¼1l¼1 i¼1

ð20aÞ "   # K X L K  T  X X X En En l l l l ~ ðjÞ F n þ m F n a þ G n þ m G n ij 2il 1 2tl 1 btj u 1 þ 1n2 1il 2 1n2 1tl 2 t¼1 j¼1l¼1 i¼1 "

#

  K X L K  T  X X X E E F nl þ mF 1il nl1 aij þ G nl þ mG1tl nl1 btj u~ ðjÞ 2 ¼ I2 2 2il 2 2 2tl 2 n n 1 1 t¼1 j¼1l¼1 i¼1

ð20bÞ R

for xk , k ¼ 1,2,. . .M O , where the integral functions Ii ¼ Os ðf~ i þ 2 ~ p ru i ÞdO and Z sl Z sl @Ri @P t ds, Gjtl ¼ ds ð21Þ F jil ¼ 0 @xj 0 @xj Consider form

nl1

¼

sinbl , nl2

¼ cosbl , we have solutions in closed

h i d1 F 1il ¼ ðr 2 r 1 Þcosbl  ðxla1 xi1 Þsinbl ðxla2 xi2 Þcosbl sinbl ln d2 h i d1 l l F 2il ¼ ðr 2 r 1 Þsinbl  ðxa2 xi2 Þcosbl ðxa1 xi1 Þsinbl cosbl ln d2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c2 þðxla1 xi1 Þ2 þ ðxla2 xi2 Þ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r 2 ¼ c2 þ ðxlb2 xi1 Þ2 þ ðxlb2 xi2 Þ2

r1 ¼

ð22Þ

ð23Þ

However, for a rectangular local integral domain as shown in Fig. 3, we have its analytical form as " # K K T X X X ðjÞ J i aij þ Ht btj ðf~ n þ u~ ðjÞ ð24Þ In ¼ n Þ j¼1

i¼1

t¼1

where   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 c x2 ln x1 þ c2 þðx1 xi1 Þ2 þ ðx2 xi2 Þ2 Ji ¼ 2  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ c2 x1 ln x2 þ c2 þ ðx1 xi1 Þ2 þ ðx2 xi2 Þ2 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ x1 x2 c2 þ ðx1 xi1 Þ2 þðx2 xi2 Þ2 3  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ x32 ln x1 þ c2 þ ðx1 xi1 Þ2 þðx2 xi2 Þ2 6  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ x31 ln x2 þ c2 þ ðx1 xi1 Þ2 þðx2 xi2 Þ2 6 0

c3 B þ tan1 @ 3

1 + cx1 C x þ D =2 x þ D =2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA jx1 D 1=2 jx2 D 2=2 1 1 2 2 c2 þ x22 þ x2 c2 þ ðx1 xi1 Þ2 þ ðx2 xi2 Þ2

and x þ D =2 x þ D =2

H1 ¼ x1 x2 jx11 D11=2 jx22 D22=2 , x1 þ D1 =2 x þ D =2 1 x1 x22 x D =2 jx22 D22=2 , 1 1 2 x1 þ D1 =2 x þ D =2 1 H5 ¼ x21 x22 x D =2 jx22 D22=2 , 1 1 4 H3 ¼

1 2 x1 þ D1 =2 x2 þ D2 =2 x x2 j 2 1 x1 D1 =2 x2 D2 =2 1 3 x1 þ D1 =2 x2 þ D2 =2 H4 ¼ x1 x2 x D =2 jx2 D2 =2 1 1 3 x1 þ D1 =2 x þ D =2 1 H6 ¼ x1 x32 x D =2 jx22 D22=2 1 1 3

H2 ¼

For the nodes on the traction boundary, (9) should be introduced: Z Z 0 t~ i dG t~ i dG ¼ Ii for xk k ¼ 1,2,. . .,MT ð25Þ G GT

GT

For the displacement boundary nodes, we can introduce the displacement equation directly, i.e. u~ i ðnk Þ ¼ u~ 0i , k ¼ 1,2,. . .M D . Therefore, there are total 2  ðM O þ M T þM D Þ linear algebraic equations which are used to determine the same number unknowns of displacements either in the domain or on the traction boundary. For linear test function, the analytical formulations were given by Wen and Aliabadi [42]. In the Laplace transform domain, a total number of Lþ1 samples in the transformation space sk , k ¼ 1, 2,. . ., L, are selected. Physical values are calculated for these transform parameters and the real value at time t must be obtained by an inverse transform. Here, the method given by Durbin [49] is used. Demonstration of

d1 ¼ ðxla1 xi1 Þcosbl þðxla2 xi2 Þsinbl þ r 1

L

b (xlb1 , xlb2)

Ωs

1

sl

2

l x1

x2

(xl1 , xl2)

xk

x2

xk

Δ1

βl

a

nl

(xla1 , xla2)

Fig. 2. Local integral domain with straight boundary lines.

x1 Fig. 3. Rectangular local integral domain.

Δ2

808

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 37 (2013) 805–811

the Durbin’s inverse method was made by Wen et al. [50] for elastodynamic problems with elasticity wave propagation. The calculation formula used is shown below: f ðtÞ ¼

"    #  L X 2eZt 1 2kp 2kpt  f~ ðZÞ þ i exp Re f~ Z þ 2 T T T k¼0

ð26Þ

where f ðpk Þ stands for the transformed variables in the Laplace pffiffiffiffiffiffiffi space for parameters pk ¼ Z þ 2kp i=T ð i ¼ 1Þ. The selection of parameters Z and T only slightly affects the accuracy of the numerical solution. In our computations, we have chosen Z ¼ 5=t0 and T=t 0 ¼ 20 in the following examples, where t 0 ¼ b=c1 is the unit of time, b is specified length and longitudinal pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi wave speed c1 ð ¼ Eð1nÞ=½ð1 þ nÞð12nÞrÞ. In addition, the dynamic stress intensity factor is evaluated by crack opening displacement (COD) in the Laplace transform domain. For mode I fracture, the stress intensity factor for plan stress problem is related to crack opening displacement, in the transformed domain, as following [35]: K~ I ¼

E 8ð1n2 Þ

sffiffiffiffiffiffi 2p Du~ 2 , r0

Du~ 2 ¼ u~ 2þ u~  2:

ð27Þ

where r 0 indicates the distance between the collocation point and crack tip, Du~ 2 is the Crack Opening Displacement (COD) in the Laplace transform domain.

5.1. Rectangular sheet under dynamic load A square plate of width a subjected to a Heaviside load s0 HðtÞ on the top edge with fixed edge on the bottom is considered, where HðtÞ is Heaviside function. The regular distribution of collocation point (11  11 points) is shown in Fig. 2(a and b) and the Poisson ratio is taken as zero to allow comparison with the exact solution. The sample number (L þ1) in the transformed domain is taken to be 51 in this example. Free parameters n and K are selected as 1 and 16, and the local integral size D1 ¼ D2 ¼ D as shown in Fig. 4(b). Analytical formulations in (20a and b) are used for the domain integral of the acceleration/force terms. The dynamic displacement on the top (A) and middle of plate (B) and the stresses on the bottom (C) and middle (B) of plate are plotted in the Figs. 5 and 6, respectively. Good agreement with analytical solution has been achieved particularly for the displacement. Similar results can be obtained for large number of collocation point in the support domain and high density of collocation point distribution, i.e. ð31  31Þ and ð41  41Þ. In addition, the approximation approach of domain integral in (23) is tested and the results are very close as well. Therefore, the approximation approach in (23) is used in the following example.

5.2. A single central crack in rectangular plate under tension Consider a rectangular plate of width 2b and length 2h with a centrally located crack of length 2a. It is loaded dynamically in the direction perpendicular to the crack by a uniform tension s0 HðtÞ on the top and the bottom. Due to symmetry, a quarter of plate is

5. Numerical examples 2.5 u2(A)

Exact (A)

u2(B)

1.5 1.0 0.5 c1t / b 0.0 0

2

4

6

8

10

12

B x2 x1

16

Fig. 5. Normalized displacement u2 =b on the top (A) and middle (B) of a square plate subjected to uniform tensile load on the top against the normalized time c1 t=b.

Δ

h

14

-0.5

σ 0 H(t)

A

Exact (B)

2.0

u2/b

In this section, the application of proposed meshless method for two-dimensional elasticity dynamic problems is demonstrated. Although the nodes can be arbitrary or uniformly distributed in the domain, in this paper, all collocation points are the uniformly distributed in the domain. Normalized free parameter in the radial basis function c=D ¼ 10n , where D indicates the minimum distance between the nodes in the local integral domain and n is free parameter (n is chosen as zero in general cases). The support domain is selected as a circle of radius R centered at field point x, which is determined such that the minimum number of nodes in the support domain K Z 16. The mode I stress intensity factor in the time domain is calculated using COD technique.

support domain (R)

local integral domain

C

Fig. 4. Square plate subjected to uniform tensile load on the top: (a) geometry and coordinate and (b) regular distribution of node.

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 37 (2013) 805–811

and 21  41 respectively. Figs. 10–12 show the normalized stress intensity factors various against the normalized time c1 t=b for different geometry of plate. In addition, the results given by Wen et al. [34] using the mesh free Perov–Galerkin method and the indirect boundary element method (fictitious load method [47]) are presented in the same figures for comparison. Apparently before the 1.6

2

1.2

3

1 K I / σ0√ π a

considered as shown in Fig. 7. Here Poisson’s ratio n ¼ 0:3 and Young’s modulus ¼1 unit. Firstly we observe the accuracy of stress intensity factor with the density of collocation point. Fig. 8 shows the normalized stress intensity factor obtained by the nodal values of COD near the crack tip for the free parameters c=D ¼ 1ðn ¼ 0Þ, the local integral size D1 ¼ D2 ¼ D and K ¼ 16. Obviously the nearest point (first) should not be considered as the singularity of stress at the crack tip and the points from the second to the sixth can be used to obtain high accuracy results. The analytical static solution for a square plate pffiffiffiffiffiffi b=h ¼ 1 containing a central crack, if a=b ¼ 0:5, is K I ¼ 1:325s0 pa [51] and the result by the second nodal value with 21  21 node pffiffiffiffiffiffi distribution is K I ¼ 1:331s0 pa. Therefore, the relative error can be expected to be less than 1% for elastostatic problems. Therefore, we use the second nodal value with 21  21 node distribution to evaluate stress intensity factor in the following examples. Secondly, we observe the sensitivity of the free parameter selection, i.e. factor c=D and K, if the local integral size D1 ¼ D2 ¼ D with 21  21 nodes. The contours of relative error in percentage are shown in Fig. 9. In the region of 0:5 on o 0:5 and K 410, we can obtain the numerical solution with degree of accuracy which is less than 5%. Three geometries of rectangular plate are considered in this example, i.e. b=h ¼ 0:5, b=h ¼ 1 and b=h ¼ 2 while a=b ¼ 0:5. The total numbers of nodes for each case are selected as 21  11, 21  21

809

Series1 21×21 0.8

Series2 3

2

31×31

1

Series3 41×41

0.4

Series4 r0 / b

0.0 0.00

0.05

0.10

0.15

0.20

0.25

0.30

Fig. 8. Stress intensity factor various with normalized distance r 0 =b between collocation point and crack tip.

3.0 u2(C)

2.5

u2(B)

Exact (B)

Exact (C)

%

1.0

σ 22/b

2.0 1.5

0.5

1.0 0.5 c1t / b 0

2

4

6

8

10

12

14

16

-0.5

Fig. 6. Normalized stress s22 =s0 on the middle (B) and bottom (C) of a square plate subjected to uniform tensile load on the top against the normalized time c1 t=b.

n

0.0

0.0

-0.5

σ0H (t) -1.0 10

15

20 K

25

30

Fig. 9. Contour of the relative error (%) various with free parameters: the number of point in support domain K and n, c=D ¼ 10n , in RBF.

5.0 This paper

h

4.0

MFPG [34] BEM [47]

x2 crack tip

a x1

KI/σ0√πa

3.0 2.0 1.0

c1t / b

0.0 0

b Fig. 7. Rectangular plate with a central crack of length 2a under tension s0 HðtÞ.

2

4

6

8

10

12

14

16

-1.0

pffiffiffiffiffiffi Fig. 10. Normalized stress intensity factor K I =s0 pa against the normalized time c1 t=b for b=h ¼ 0:5.

810

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 37 (2013) 805–811 3.5 This paper

3.0

MFPG [34]

2.5

BEM [47]

KI/σ0√πa

2.0 1.5 1.0 0.5

c1t / b

0.0 -0.5

0

2

4

6

8

10

12

14

16

-1.0

pffiffiffiffiffiffi Fig. 11. Normalized stress intensity factor K I =s0 pa against the normalized time c1 t=b for b=h ¼ 1:0. 3.0 2.5 This paper

K I / σ 0 √π a

2.0

MFPG [34] BEM [47]

1.5 1.0 0.5

c1t / b

0.0 0

2

4

6

8

10

12

14

16

-0.5 -1.0

pffiffiffiffiffiffi Fig. 12. Normalized stress intensity factor K I =s0 pa against the normalized time c1 t=b for b=h ¼ 0:5.

arrival time of dilatation wave traveling from the top of plate, the stress intensity factor should remain to be zero. The agreement between these solutions is considered to be good.

6. Conclusion This paper presents a meshless local integral method to twodimensional elastodynamic fracture problems by the Laplace transform technique. Considering a local integral domain with local support domain, the analytical formulations were derived in the Laplace transform domain. Durbin’s inversion method was applied to determine all variable in the time domain. The dynamic stress intensity factor of mode I was evaluated by using the COD technique. We can conclude with the following observations: (1) meshless local integral method is valid to deal with elastodynamic problem in Laplace space; (2) analytical formulations of all integrals save CPU time; (3) numerical solutions are stable and convergent for suitable selection of free parameters; (4) this method can be extended to elastoplastic, plate/shell bending and nonlinear problems directly. References [1] Aliabadi MH. The boundary element method.Applications in solids and structures, vol. 2. Chicester: Wiley; 2002. [2] Portela A, Aliabadi MH, Rooke DP. The dual boundary element method: effective implementation for crack problems. Int J Numer Methods Eng 1992;33:1269–87. [3] Portela A, Aliabadi MH, Rooke DP. Dual boundary element analysis of cracked plates: singularity subtraction technique. Int J Fract 1992;55:17–28. [4] Wen PH. Dynamic fracture mechanics: displacement discontinuity method.Topics in engineering, vol. 29. Southampton: Computational Mechanics Publications; 1996.

[5] Wrobel LC. The boundary element method.Application to thermo-fluids, vol. 1. Chicesater: Wiley; 2002. [6] Brancati A, Aliabadi MH, Bendetti I. Hierarchical adaptive cross approximation GMRES technique for solution of acoustic problems using the boundary element method. CMES Comput Model Eng Sci 2009;43:149–72. [7] Brancati A, Aliabadi MH, Mallardo MHV. A BEM sensitivity formulation for three-dimensional active noise control. Int J Numer Methods Eng 2012;90(9):1183–206. [8] Mallardo V, Aliabadi MH. An adaptive fast multipole approach to 2D wave propagation. CMES Comput Model Eng Sci 2012;87(2):77–96. [9] Telles JCF, Brebbia CA. On the application of the boundary element method to plasticity. Appl Math Model 1979;3:466–70. [10] Leitao V, Aliabadi MH, Rooke DP. The dual boundary element formulation for elastoplastic fracture mechanics. Int J Numer Methods Eng 1995;38:315–33. [11] Mallardo V. Integral equations and nonlocal damage theory: a numerical implementation using the BDEM. Int J Fract 2009;157(1–2):13–32. [12] Pindea E, Aliabadi MH. Dual boundary element analysis for time-dependent fracture problems in creeping materials. Key Eng Mater 2008;383:109–21. [13] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 1992;10:307–18. [14] Belytschko T, Lu YY, Gu L. Element-free Galerkin method. Int J Numer Methods Eng 1994;37:229–56. [15] Liu WK, Jun S, Zhang Y. Reproducing kernel particle methods. Int J Numer Methods Eng 1995;20:1081–106. [16] Sladek V, Sladek. J, Zhang Ch. Comparative study of meshless approximations in local integral equation method. CMC Comput Mater Continua 2006;4:177–88. [17] Atluri SN. The meshless method (MLPG) for domain and BIE discretizations. Forsyth, GA, USA: Tech Science Press; 2004. [18] Wen PH, Aliabadi MH. An improved meshless collocation method for elastostatic and elastodynamic problems. Commun Numer Methods Eng 2008;24(8):635–51. [19] Zhang J, Yao Zhenhan, Tanaka M. The meshless regular hybrid boundary node method for 2D linear elasticity. Eng Anal Boundary Elem 2003;27:259–68. [20] Zhang Xi, Yao Zhenhan, Zhang Zhangfei. Application of MLPG in large deformation analysis. Acta Mech Sin (English Series) 2006;22:331–40. [21] Miers LS, Telles JCF. On the NGF procedure for LBIE elastostatic fracture mechanics. Comput Model Eng Sci 2006;14:161–9. [22] Khosravifard A, Hematiyan MR, Marin L. Nonlinear transient heat conduction analysis of functionally graded materials in the presence of heat sources using an improved meshless radial point interpolation method. Appl Math Model 2011;35:4157–74. [23] Oh Hae-Soo, Davis C, Kim JG, Kwon Y-H. Reproducing polynomial particle methods for boundary integral equations. Comput Mech 2011:1–19. [24] Li Xiaolin, Li Shuling. A meshless method for nonhomogeneous polyharmonic problems using method of fundamental solution coupled with quasiinterpolation technique. Appl Math Model 2011;35:3698–709. [25] Barbieri E, Meo MA. A meshless cohesive segments method for crack initiation and propagation in composites. Appl Compos Mater 2011;18:45–63. [26] Skouras ED, Bourantas GC, Loukopoulos VC, Nikiforidis GC. Truly meshless localized type techniques for the steady-state heat conduction problems for isotropic and functionally graded materials. Eng Anal Boundary Elem 2011;35:452–64. [27] M.B., Shariati M, Eslami MR, Hassani B. Meshless analysis of cracked functionally graded materials under thermal shock. Mechanika 2010;4:20–7. [28] Ferronato M, Mazzia A, Pini GA. Finite Element enrichment technique by the Meshless Local Petrov–Galerkin method. Comput Model Eng Sci 2010;62:205–23. [29] Wen PH, Aliabadi MH. Elastic moduli of woven fabric composite by meshless local Petrov–Galerkin (MLPG) method. Comput Model Eng Sci 2010;61:133–54. [30] Wen PH, Aliabadi MH, Liu YW. Meshless method for crack analysis in functionally graded materials with enriched radial base functions. Comput Model Eng Sci 2008;30:133–47. [31] Wen PH, Aliabadi MH. An improved meshless collocation method for elastostatic and elastodynamic problems. Commun Numer Methods Eng 2008;24:635–51. [32] Sladek J, Sladek V, Wen PH, Aliabadi MH. Meshless local Petrov–Galerkin (MLPG) method for shear deformable shells analysis. Comput Model Eng Sci 2006;13:103–17. [33] Wen PH, Aliabadi MH. A variational approach for evaluation of stress intensity factors using the element free Galerkin method. Int J Solids Struct 2011;48:1171–9. [34] Wen PH, Aliabadi MH. Evaluation of mixed-mode stress intensity factors by the mesh-free Galerkin method: static and dynamic. J Strain Anal Eng Des 2009;44:273–86. [35] Wen PH, Aliabadi MH. Meshless method with enriched radial basis functions for fracture mechanics. SDHM Struct Durability Health Monit 2007;3:107–19. [36] Lancaster P, Salkauskas K. Surfaces generated by moving least square methods. Math Comp 1981;37:141–58. [37] Golgerg MA, Chen CS, Karur SR. Improved multiquadric approximation for partial differential equations. Eng Anal Boundary Elem 1996;18:9–17. [38] Li LY, Wen PH, Aliabadi MH. eshfree modeling and homogenization of 3D orthogonal woven composites. Compos Sci and Technol 2011;71(15):1777–88. [39] Wen PH, P.H, Aliabadi MH. Damage mechanics analysis of plain woven fabric composite micromechanical model for mesh-free simulations. J Compos Mater 2012;46(18):2239–53.

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 37 (2013) 805–811

[40] Sladek V, Sladek J. Local integral equations implemented by MLS— approximation and analytical integrations. Eng Anal Boundary Elem 2010;34:904–13. [41] Sladek V, Sladek J, Zhang Ch. On increasing computational efficiency of local integral equation method combined with meshless implementations. CMES Comput Model Eng Sci 2010;63:243–63. [42] Wen PH, Aliabadi MH. Analysis of functionally graded plates by meshless method: a purely analytical formulation. Eng Anal Boundary Elem 2012;36:639–50. [43] Soares D, Sladek V, Sladek J. Modified meshless local Petrov–Galerkin formulations for elastodynamics. Int J Numer Methods Eng 2012;90:1508–28. [44] Soares Jr. D, Sladek V, Sladek J, Zmindak M, Medvecky S. Porous media analysis by modified MLPG formulations. CMC Comput Mater Continua 2012;27:101–27.

811

[45] Fleming M, Chu YA, Moran B, Belytschko T, Lu YY, Gu L. Enriched elementfree Galerkin methods for crack-tip fields. Int J Numer Methods Eng 1997;40:1483–504. [46] Rao BN, Rahman SA. Coupled meshless-finite element method for fracture analysis of cracks. Int J Pressure Vessels Piping 2001;78:647–57. [47] Wen PH, Aliabadi MH, Rooke D. The influence of elastic waves on dynamic stress intensity factors (two dimensional problem). Arch Appl Mech 1996;66(5):326–35. [48] Hon YC, Mao XZ. A multiquadric interpolation method for solving initial value problems. J. Sci Comput 1997;12:51–5. [49] Durbin F. Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate’s method. Comput J. 1975;17:371–6. [50] Hardy RL. Multiquadric equations of topography and other irregular surface. J Geophys Res 1971;76:1905–15. [51] Rooke DP, Cartwright DJ. Compendium of stress intensity factors. London: Her Majesty’s Stationery Office; 1976.