A numerical Green's function and dual reciprocity BEM method to solve elastodynamic crack problems

Engineering Analysis with Boundary Elements 29 (2005) 204–209 www.elsevier.com/locate/enganabound

A numerical Green’s function and dual reciprocity BEM method to solve elastodynamic crack problems C.A.R. Vera-Tudela, J.C.F. Telles* COPPE—Universidade Federal do Rio de Janeiro, Programa de Engenharia Civil, Caixa Postal 68506, CEP 21945-970 Rio de Janeiro, RJ, Brazil Received 25 September 2004; revised 13 January 2005; accepted 13 January 2005 Available online 10 March 2005

Abstract A computational model based on the numerical Green’s function (NGF) and the dual reciprocity boundary element method (DR-BEM) is presented for the study of elastodynamic fracture mechanics problems. The numerical Green’s function, corresponding to an embedded crack within the infinite medium, is introduced into a boundary element formulation, as the fundamental solution, to calculate the unknown external boundary displacements and tractions and in post-processing determine the crack opening displacements (COD). The domain inertial integral present in the elastodynamic equation is transformed into a boundary integral one by the use of the dual reciprocity technique. The dynamic stress intensity factors (SIF), computed through crack opening displacement values, are obtained for several numerical examples, indicating a good agreement with existing solutions. q 2005 Elsevier Ltd. All rights reserved. Keywords: Numerical Green’s function; Dual reciprocity; Dynamical crack problems

1. Introduction Accurate crack analysis is important for health monitoring and quality control in many critical applications pertaining to civil, aeronautical, nuclear and mechanical engineering areas. The appropriate numerical modelling of linear elastic fracture problems, however, strongly depends on the appropriate simulation of the singular stress field in the vicinity of the crack tip. With regards to the boundary element method (BEM), another difficulty arises: the crack surfaces geometrically coincide in the numerical model, causing a degeneration of the boundary integral equation. This problem is usually avoided by three alternative approaches: applying the sub-regions technique [1], employing the mixed or dual formulation [2,3,4], or using the associated Green’s function [2,5]. Crack modelling using the Green’s function possesses the advantage of opening for the possibility of implementing only the classical boundary integral equation to represent * Corresponding author. Tel.: C55 21 2562 7383; fax: 55 21 2562 8464/55 21 2290 6626. E-mail address: [email protected] (J.C.F. Telles). 0955-7997/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2005.01.004

and solve fracture mechanics applications. An important feature of this formulation is the inherent elimination of problem COD integrations, avoiding standard element discretization over crack surfaces. The numerical Green’s function boundary element procedure (NGF) [5,6] maintains most of the benefits of the closed form Green’s function implementation, with the advantage of being a lot more general and competitively efficient, even in comparison with simple geometry cracks where analytical Green’s functions can be found [5]. In addition to the actual crack simulation procedure, elastodynamic fracture mechanics problems may require alternative means of dealing with the domain inertial terms that appears if a non-elastodynamic fundamental solution is adopted [7]. In this respect, the dual reciprocity technique [8,9] appears as a practical alternative to implement in conjunction with the NGF approach to solve elastodynamic crack problems using an elastostatic fundamental solution. The present paper discusses the development of the numerical Green’s function idea in conjunction with the dual reciprocity technique to solve linear elastodynamic crack problems. The text is organized in such a way as to emphasize the new formulation aspects of the combination of both procedures and introduces a series of examples to illustrate the applicability of the NGF/DR-BEM.

C.A.R. Vera-Tudela, J.C.F. Telles / Engineering Analysis with Boundary Elements 29 (2005) 204–209

pK

-p K

ik

Γ

Ω

ik

f

A

= ξ

205

Fi

Ω

Ω

+

Γ

Fi

ξ

B f

B

A

Fig. 1. Unit point Pi applied within infinite space with a traction-free crack. (A) Kelvin’s problem; (B) Complementary problem.

2. The numerical Green’s function The Green’s function boundary integral equation, for crack problems, is obtained if the fundamental solution is written in terms of a superposition of the Kelvin solution plus a complementary part which provides satisfaction of the traction-free requirement over the crack surface. Fig. 1 illustrates this idea and the traction-free condition removes the need to formally discretize the crack surface in the boundary element formulation [5]. The Green’s function can be represented as k c G k c uG ij ðx; xÞZ uij ðx; xÞ C uij ðx; xÞ pij ðx; xÞ Z pij ðx; xÞ C pij ðx; xÞ

(1) uG ij ðx; xÞ

pG ij ðx; xÞ

and are the fundamental displacewhere ment and tractions, in j direction at the field point x due to a unit point load applied at the source point x in the i direction. The superscripts k and c stand for Kelvin and complementary components, respectively, of the fundamental problem. The Kelvin components ukij ðx; xÞ and pkij ðx; xÞ are known [10], the complementary displacements and tractions, ucij ðx; xÞ and pcij ðx; xÞ, are the unknowns of problem (B) (see Fig. 1). This solution can be written in terms of boundary integral equations as shown bellow [6]: ð ucij ðx; xÞ Z pkjk ðx; zÞcik ðx; zÞdGðzÞ; x ;GI (2) GI

pcij ðx; xÞ Z

ð

Fig. 1 shows that the natural conditions of problem (B) are prescribed, hence the limit of Eq. (3) as x/GI can be obtained to produce a hyper-singular boundary integral equation for the problem. This procedure has been discussed before [5,6] and yields in the limit:
Z Fp pcij ðx; zÞ


zÞcik ðx; zÞdGðzÞ; z
2GI Pkjk ðz;

(4)

GI

Ð where the symbol Fp indicates Hadamard’s finite part integral. k

Since pcij ðx; zÞZKp ij ðx; zÞ, Eq. (4) produces a hypersingular boundary integral equation for the unknown crack opening displacements: ð Fp


zÞcik ðx; zÞdGðzÞ Z Kpkij ðx; zÞ
Pkjk ðz;

(5)

GI

Eq. (5) can now be solved by a standard weighted residual method, using the point collocation technique. Therefore, using the Dirac delta function to weight the equation at a certain number of points z
m ðmZ 1; .; MÞ on GI, the following equations are found: ð Fp

Pkjk ðx; zÞcik ðx; zÞdGðzÞ; x ;GI

ð

Pkjk ðz
m ; zÞcik ðx; zÞdGðzÞ Z Kpkij ðx; z
m Þ; m Z 1; .; M

GI

(3)

(6)

where cik ðx; zÞZ ucik ðx; zs ÞK ucik ðx; zI Þ is the crack opening displacements of the Green’s function. Here, S and I stand for ‘superior’ and ‘inferior’ surfaces of the crack as indicated in Fig. 2. The sign of the integrals depend on the chosen surface of integration [5]; in this case GI has been adopted. The integral Eq. (3) is originated from the hypersingular or traction formulation. Eqs (2) and (3) produce the complementary displacements and tractions as a function of the fundamental crack openings. In order to solve the problem, it is therefore necessary to determine the crack opening displacements.

Standard Gauss quadrature, in principle, cannot calculate finite part integrals. Hence, a rectifying scheme can be

GI

ΓS Ω ΓI Fig. 2. Infinite domain with crack boundary GF Z GI g GS .

206

C.A.R. Vera-Tudela, J.C.F. Telles / Engineering Analysis with Boundary Elements 29 (2005) 204–209

adopted in the following fashion [6]: N X

ones (MZN), leading to N X

k jJ jn Pjk ðz
m ; zn Þcik ðx; zn ÞWn K Eij ðx; z
m Þ

k jJ jn Pjk ðz
m ; zn Þcik ðx; zn ÞWn K Eij ðx; z
m Þ

nZ1

nZ1

Z Kpkij ðx; z
m Þ; m Z 1; .; M

(7)

where jJ jn is the Jacobian of the transformation to the standard quadrature interval and zn the corresponding point at the Gauss station n, Wn is the associated weighting factor at this station and N is the total number of integration points. The correcting term, associated to the singular behaviour of the integrand, is Eij ðx; z
m Þ. In order to define Eij, the crack opening displacements cik can be expanded in Taylor’s series about the point z
m : vc ðx; z
m Þ 1 ½GðzÞ K Gðz
m Þ C cik ðx; zÞ Z cik ðx; z
m Þ C ik vGðzÞ 2 !

v2 cik ðx; z
m Þ ½GðzÞ K Gðz
m Þ2 C . vGðzÞ2

(8)

Then, the next step is to substitute Eq. (8) in Eq. (6). Since Pkjk Z OðrK2 Þ for 2-D applications, only the first two terms of the series produce singular behaviour; the first generates a finite part integral and the second a Cauchy principal value. Hence, Eij can be defined as

Z Kpkij ðx; z
m Þ; m Z 1; .; N

Careful examination of index expansions in Eq. (12) indicates that the system of 4N equations can be subdivides further into two systems, of dimensions 2N, sharing the same square matrix S, as is showed in the next equation: Scij ðxÞ Z pij ðxÞ

vcik ðx; z
m Þ ð2Þ ejk vGðzÞ

k uG ij ðx; xÞ Z uij ðx; xÞ C

k jJ jn pjk ðx; zn Þcik ðx; zn ÞWn

nZ1

(14) Z pkij ðx; xÞ C

N X

k jJ jn Pjk ðx; zn Þcik ðx; zn ÞWn

nZ1

3. The dual reciprocity BEM for crack problems N X

k jJ jn Pjk ðz
m ; zn ÞWn K Fp

nZ1

ð

Pkjk ðz
m ; zÞdGðzÞ

(10)

The dynamic equilibrium of a body without body force can be expressed as:

GI

Gui;jj C and eð2Þ jk Z

N X

(9)

where eð1Þ jk Z

(13)

It is important to note that matrix S is only a function of crack geometry and remains the same for any position and direction of the unit point load. Furthermore, Eq. (13) is repeatedly solved applying just back substitution in a Gauss solution routine to calculate all needed cij(x) vectors. Finally, the numerical Green’s function, for a crack simulation, is written using Eqs (2) and (3) and the calculated fundamental crack opening displacements at the discrete standard Gauss point positions:

pG ij ðx; xÞ Eij ðx; z
m Þ Z cik ðx; z
m Þeð1Þ jk C

(12)

N X

k jJ jn Pjk ðz
m ; zn Þ½Gðzn Þ K Gðz
m ÞWn K Cp

nZ1

ð !

Pkjk ðz
m ; zÞ½GðzÞ K Gðz
m ÞdGðzÞ

(11)

GI

Ð The finite part and Cauchy principal value (Cp ) indicated as the last term of Eqs (10) and (11) are to be computed either analytically or using accurate numerical schemes. In order to generate a square system of equations for cik at the Gauss point positions, the number of collocation points is taken to be the same as the number of Gauss integration

G u Z r€u i ð1 K 2nÞ j;ji

(15)

where n is Poisson’s ration, G is the shear modulus, r is the mass density and the dots denote differentiation with respect to time. The dual reciprocity procedure is introduced to transform the inertial integral involving the acceleration term u€ j , in a boundary integral [9]. This acceleration term can be approximated by a series of geometric and time functions f and a: u€ j ðx; tÞ Z

n X

m am j ðtÞf ðxÞ

(16)

mZ1

where n is the number of boundary plus internal nodes, m indicate the nodal point number, f is an approximate function and aj are the coefficients to be determined with the collocation method.

C.A.R. Vera-Tudela, J.C.F. Telles / Engineering Analysis with Boundary Elements 29 (2005) 204–209

The boundary integral equation for a crack problem can be written as: ð

ð

pG ij uj dG K

Cij uj C GE CGF

ð n X B m Zr @Cij u^ kj C mZ1

Discretizing the boundary in Ne external elements, for general loaded crack problems, Eq. (17) takes de form: Cu C

Ne X

uG ij pj dG 1

ð

^m pG ij u kj dG K

GE CGF

Ne X

hj uj K

jZ1

GE CGF

0

Zr

C m ^m uG ij p kj dGAak

GE CGF

(17)

207

n X

gj pj K ½uG pGF

jZ1

m

Cu^ C

mZ1

Ne X

hj u^ m j

K

jZ1

Ne X

! gj p^ m j

m

K ½cp^ GI am

jZ1

(21) Finally, from the application of Eq. (21) to all nodes and internal points, a system of equations is obtained in the form, Hu K Gp K g Z rðHu^ K Gp^ K lÞa

^m where u^ m kj and p kj are de displacement and traction particular E solutions, G the external boundary of the problem and GF the crack boundary. Eq. (17) is therefore rewritten as:

GE n X mZ1

ð K

pG ij uj dG K

GF

0

B m @Cij u^kj C

ð

ð

GE

^m pG ij u kj dG C

ð

uG ij pj dG K

GE

GE

^m uG ij p kj dG K

ð

Hu K Gp K g Z rðHu^ K Gp^ K lÞfK1 u€ uG ij pj dG

GF

ð

^m pG ij p kj dG

C m ^m uG ij p kj dGAak

ð18Þ

GF

^m uG ij p kj dG Z

GF

ð

^m uG ij p kj dG Z

ð

Z

ð



GS uGI p^m ij K uij kj dG Z

GI

4. Numerical examples Three numerical examples are presented, in dimensionless form, to illustrate de applicability of the formulation. P=1

 mI GS mI ^ ^ K u p p uGI dG ij kj ij kj

GI

GI CGS

Upon introduction of the boundary conditions prescribed over the external boundary, the system of equations of Eq. (23) can be solved by a suitable time marching procedure, such as the Houbolt scheme [9] adopted here. For frequency (u) domain problems, Ku2u substitutes for u€ in Eq. (23).

GF

1

The second and sixth integrals in Eq. (18) now vanish due to the traction-free condition of the Green’s function. The fourth integral is also zero for unloaded cracks. This condition, however, is found different in the last integral of Eq. (18). Here, the integral do not vanish but can be simplified as follows: ð

ð

2

Fig. 3. Central crack in a plate.

cij p^m kj dG

ð19Þ 80

GI

70

ð

G G Cij uj C pG ij uj dGK uij pj dGK uij pj dG GE

0

GE

GE

30

10

GF

GE

40

20

1

ð ð ð n X B C m m G m G m ^ ^ ^ Zr C C p dGK u dGK cij p^ m u u p @ ij kj ij kj ij kj kj dGAak mZ1

Chirino et al.(1994) and Barra (1996)

50

k1/k0

ð

NGF/DR-BEM

60

In this way, the boundary integral equation for a crack problem can be written as: ð

(23)

1

pG ij uj dG C

Cij uj C

Zr

ð

where g is the vector corresponding to the prescribed tractions over the crack surface and l stands for the last term of Eq. (21). Remembering that by definition aZ fK1 u€ , the system of equations is

2a

ð

(22)

GI

(20)

0 –10

0

0.2

0.4

0.6

0.8

1

1.2

wa/c2

Fig. 4. Normalized K1/K0 for center cracked plate under uniform harmonic traction.

C.A.R. Vera-Tudela, J.C.F. Telles / Engineering Analysis with Boundary Elements 29 (2005) 204–209

0.5h

2b

p

2a

2h

2a

1.5h

2b

208

Fig. 7. Geometry of Chen’s problem.

Fig. 5. Eccentric crack with Heaviside load.

Here, the dynamic SIF values are presented in a pnonffiffiffiffiffiffi dimensional form divided by the static SIF, K0 Z s ap, where s is the remote applied stress and 2a is the total crack length. In all examples the fundamental crack opening integration has been carried out using 20 Gauss points and four internal nodal points, evenly distributed in the vicinity of the respective cracks, have been adopted.

5. Frequency response A rectangular plate is studied and the geometry is shown in Fig. 3, the material properties are Young modulus EZ 1000, Poisson’s ratio nZ0.3 and mass density rZ1.0. The external boundary was discretized using 80 uniform quadratic elements. The plate is loaded along the opposite boundaries by uniform time harmonic tractions. In Fig. 4, the numerical results of the problem are plotted and are compared with those obtained by Chirino et al. [11] and Barra and Telles [12]. The peaks represent the symmetric mode resonance frequencies.

6. Analytical limited domain solution

and boundary conditions are shown in Fig. 5. There is a stage between the initial loading and the onset of external boundary reflexions reaching the crack, where the problem results coincide with those of the companion infinite domain problem [12]. The dimensions of the plate are bZ5, hZ10 and crack length 2aZ2.4. The material constants are: EZ7.43, nZ0.3 and rZ1. There are 24 quadratic boundary elements of equal size L and the time stepffi is DtZ0.35 (i.e., pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi CpDtZ0.40, CpZ 2Gð1K nÞ=ðrð1K 2nÞ). The results are presented in Fig. 6 and compared with those obtained in References [11] and [12].

7. Chen’s problem A finite plate with a single central crack is loaded along opposite boundaries, at time tZ0, by Heaviside tractions p as shown in Fig. 7. This problem was solved by Chen in 1975 [13] using finite differences. It is a classical problem and serves as a reference to several authors. The dimensions of the plate are bZ5 and hZ10; the crack length is 2aZ2.4. The material is linear elastic with the following properties:

In this example a crack, in a limited domain, is subjected to unit distributed Heaviside tractions. The geometry

3 Chen

2.5

1.6 1.4

NGF/DR-BEM

1.2

NGF/Freq. DomainBEM

NGF/DR-BEM

2

k1/k0

1 k1/k0

p

0.8

1.5 1

0.6 0.5

0.4 0.2

0

0

0

1

2

–0.2

0

3

1

2

3

–0.5

cpt/L Fig. 6. Normalized time response for K1/K0.

ct/L Fig. 8. Dynamic stress intensity factor for Chen’s problem.

4

C.A.R. Vera-Tudela, J.C.F. Telles / Engineering Analysis with Boundary Elements 29 (2005) 204–209

EZ7.43, nZ0.3 and rZ1.0. The boundary is discretized using 80 uniform quadratic elements and DtZ0.20. The boundary element formulation with the numerical Green’s function and dual reciprocity method (FGN/DR-BEM) is compared with those obtained by Chen in Fig. 8 where the normalized dynamic SIF values are depicted.

8. Conclusions A numerical Green’s function in conjunction with the dual reciprocity boundary element procedure has been described for 2D elastodynamic fracture mechanics problems. The accuracy of the NGF/DR-BEM can be considered satisfactory as indicated in the examples. The numerical Green’s function eliminates the problem unknowns over the crack boundaries and integration along the crack surface is only necessary when non-zero prescribed tractions exist, contributing only to the independent term of the final system of equations.

Acknowledgements This work has been carried out with the support of the CNPq (Brazilian Council for Scientific and Technological Research).

209

References [1] Blandford GE, Ingraffea AR, Ligget JA. Two-dimensional stress intensity factor computation using the boundary element method. Int J Numer Methods Eng 1981;17:387–404. [2] Cruse TA. Boundary element analysis in computational fracture mechanics. Dordrecht: Kluwer; 1998. [3] Guimara˜es S, Telles JCF. On the hyper-singular boundary element formulation for fracture mechanics applications. Eng Anal Bound Elem 1994;13:353–63. [4] Portela A, Aliabadi MH, Rooke DP. Dual boundary element analysis of cracked plates: singularity subtraction technique. Int J Fract 1992; 55:17–28. [5] Telles JCF, Castor GS, Guimara˜es S. A numerical Green’s function approach for boundary elements applied to fracture mechanics. Int J Numer Methods Eng 1995;38:3259–74. [6] Telles JCF, Guimara˜es S. Green’s function: a numerical generation for fracture mechanics problems via boundary elements. Comput Methods Appl Mech Eng 2000;188:847–58. [7] Beskos DE. Boundary element methods in dynamic analysis: Part II (1986–1996). Appl Mech Rev 1997;50:149–97. [8] Gaul M, Kogl M, Wagner M. Boundary element methods for engineers and scientists. Berlin: Springer; 2003. [9] Partridge PW, Brebbia CA, Wobel LC. The dual reciprocity boundary element method. London: Computational Mechanics Publications; 1992. [10] Brebbia CA, Telles JCF, Wrobel LC. Boundary elements techniques: theory and applications in engineering. Berlin: Springer; 1984. [11] Chirino F, Gallego R, Sa´ez A, Domı´nguez J. A comparative study of the three boundary element approaches to transient dynamic crack problems. Eng Anal Bound Elem 1994;13:11–19. [12] Barra LPS, Telles JCF. A hyper-singular numerical Green’s function generation for BEM applied to dynamic SIF problems. Eng Anal Bound Elem 1999;23:77–87. [13] Chen YM. Numerical computation of dynamic stress intensity factor by Lagrangian finite-difference method. Eng Fract Mech 1975;7: 653–60.