Int J Fract (2013) 183:169–186

DOI 10.1007/s10704-013-9885-5

ORIGINAL PAPER

An efficient numerical method for the solution of the second boundary value problem of elasticity for 3D-bodies with cracks

S. Kanaun · A. Markov · S. Babaii

Received: 26 September 2012 / Accepted: 18 September 2013 / Published online: 15 October 2013 © Springer Science+Business Media Dordrecht 2013

Abstract The second boundary value problem of elasticity for 3D-bodies containing cracks is considered. Presentation of the solution in the form of the double layer potential reduces the problem to a system of 2D-integral equations which kernels are similar for the body boundary and crack surfaces. For discretization of these equations, Caussian approximation functions centered at a set of nodes homogeneously distributed on the body and crack surfaces are used. For such functions, calculation of the elements of the matrix of the discretized problem is reduced to five standard 1Dintegrals that can be tabulated. For planar cracks, these integrals are calculated in closed analytical forms. The method is mesh free, and for its performing, only node coordinates and surface orientations at the nodes should be defined. Calculation of stress intensity factors at the crack edges in the framework of the method is discussed. Examples of an elliptical crack, a lens-shaped crack, and a spherical body subjected to concentrated and distributed surface forces are considered. Numerical results are compared with the solutions of other authors presented in the literature. Convergence of the method with respect to the node grid steps is analyzed.

S. Kanaun (B)

The Mexican Oil Institute, Lazario Cardenes 25,

D.F., Mexico 52399, Mexico e-mail: kanaoun@itesm.mx

A. Markov · S. Babaii

Technological Institute of Higher Education of Monterrey, State of Mexico Campus, Atizapan, Edo de México 52926, Mexico

An efficient algorithm of the node grid generation is proposed.

Keywords Elasticity · Second boundary value problem · Crack problem · Gaussian approximating functions · Stress intensity factor 1 Introduction

The 3D-problem of elasticity for homogeneous bodies containing cracks has many important applications in geophysics, fracture and damage mechanics, micromechanics of cracked media, etc. By numerical solution of this problem, the boundary element method (BEM) has some advantages over the finite element method (FEM). In the BEM, the problem is reduced to 2Dintegral equations on the body and crack surfaces, and the principal unknowns are vectors that have no singularities on the crack edges. Whereas the FEM algorithms operate with 3D-displacement field and its derivatives that have singularities in the vicinity of the crack tips. Such singularities are obstacles for constructing accurate numerical solutions. Usually for numerical solution of boundary integral equations, the integration area (the body boundary and crack surfaces) is divided into boundary elements (see, e.g., Cruse 1988). Inside each element, the solution is approximated by standard functions (polynomial splines or non-uniform rational

B-splines: see Simpson et al. 2013) with unknown coefficients. The discretized problem is a linear algebraic 123 170 S. Kanaun et al. system for the coefficient of the approximation. The components of the matrix of the discretized problem are integrals over the boundary elements. For polynomial approximating functions and curvilinear boundary elements, these integrals can be calculated only numerically. As the result, a great portion of the computer time is spent on the construction of the matrix of the discretized problem. Division of the body boundary and crack surfaces into the boundary elements is an independent and laborious problem.

By analysis of fracture processes in bodies with cracks, there appears a problem of accurate calculation of stress intensity factors (SIF) on the crack edges. Several strategies were proposed to modify standard BEM and FEM in order to increase accuracy of SIF calculations (see, e.g., Farhat et al. 2001, 2003; Rabczuk et al. 2004; Duflot 2006; Simpson and Trevelyan 2011;

Barbieri et al. 2012). In all these strategies, special finite elements are introduced near crack edges for correct description of the fields in this region. Introducing such elements complicates the computational algorithms, and in the most works, the modified methods were carried out for the solution of 2D-problems only (1D-cracks).

In the present work, we develop a numerical method that essentially reduces the difficulties in the numerical solution of the integral equations of the problem under consideration. First, the problem is formulated in terms of 2D-integral equations which kernels are similar for the body boundary and crack surfaces. Second, for the numerical solution of these equations, a class of

Gaussian approximating functions is used. The theory of approximation by Gaussian and other similar functions was developed by Maz’ya and Schmidt (2007).

Employing these functions for the solution of boundary integral equations has many important advantages. The method based on these functions is mesh free, and only node coordinates and surface orientations at the nodes are necessary geometrical information for its performing. Gaussian functions essentially simplify construction of the matrix of the discretized problem. The elements of this matrix are combinations of five standard absolutely converging 1D-integrals that can be tabulated for small values of the arguments and have simple analytical asymptotics for large arguments. Thus, these elements are calculated fast.

For planar cracks of arbitrary shapes, the discretized problem is essentially simplified because the mentioned integrals are calculated in closed analytical forms. In addition, for regular node grids, the matrix of the discretizad problem has the Teoplitz structure, and only a row and a column of such a matrix should be computed. Because for Teoplitz matrices, the Fast