Boundary element analysis of an integral equation for three-dmensional crack problems

In another paper, the authors proposed an integral equation for arbitrary shaped three-dimensional cracks. In the present paper, a discretization of this equation using a tensor formalism is formulated. This approach has the advantage of providing the displacement discontinuity vector in the local basis which varies as a function of the point of the crack surface. This also facilitates the computation of the stress intensity factors along the crack edge. Numerical examples reported for a circular crack and a semi-elliptical surface crack in a cylindrical bar show that one can obtain good resuhs, using few Gaussian poilnts and no singular elements.


INTRODUCTION
The integral equation method has been thoroughly experienced in structure analysis. Integral equations derived from the Somigliana representation were successfully applied to threedimensional problems, such as in References 1 and 2. In the latter reference, a particular symmetrical problem of cracked bodies was also investigated. Integral equations with kernels containing singular solutions, usually referred to as Kupradze elastic potentials, 3 are particularly well suited to crack analysis. The use of these potentials in linear fracture mechanics is advantageous as it permits one to obtain directly stress intensity factors from the computed vector density. Embedded cracks in an infinite medium were studied in References 4-7. In particular, Bui 4 showed that for plane cracks the mode I is entirely uncoupled from modes II and III. Later on, integral equations for three-dimensional cracks were proposed; 8 • 9 in the last reference use was made of a Kupradze double layer potential, and the unknown was the vector density directly related to the displacement discontinuity through the crack surface. In this paper, we carry out the discretization of the integral equation proposed in Reference 9, with a view to studying imbedded or surface crack problems.
For numerical purposes, the main difficulty is the reckoning of singular integrals defined in the sense of the principal value. Cruse 1 succeeded in giving the closed form of the principal value integral, using planar boundary elements. A more general way to evaluate two-dimensional singular integrals was investigated by K.azantzakis and Theocaris, 10 who also reduced these to one-dimensional finite-part integrals. 11 Among the most recent works, one can find References 12 and 13 where exact expressions for some specific two-dimensional finite!-part integrals are derived. 1 As for Kupradze potentials, Bui 4 took a small polygon instead of a circle to perform principal value integrals, it seems, however, that this procedure cannot be used with ease. Lastly, one can find in Reference 14 an improvement of the evaluatjon of integrals of order 1 /r, 1 fr 2 , 1 fr 3 , by means of a minimum number of Gaussian points. As the considered distance r, small though it may be, remains always finite, these integrals are actually not singular ones. In any case, difficulties must certainly be expected when dealing with principal value integrals defined on curved surfaces. In our analysis, we limit consideration to classical Gaussian quadrature, with no actual principal value computations. Nevertheless, two numerical examples given in the last section show that one can obtain fair results, using few Gaussian points, no singular integral computations and in particular, no singular elements around the crack edge. First, we deal with the embedded pennyshaped crack. The numerical results obtained proved to fit well to the analytical solution of Sneddon. 15 Then we treat the problem of a semi-elliptical crack in a cylindrical bar under combined loads. In the opening mode, the numerical results can be compared to the existing results 16 -18 which involved analogous crack geometries.

THEORETICAL BACKGROUND
Let us consider the problem of an elastic solid ~ containing an arbitrary shaped, imbedded or surface crack S, with the mixed boundary conditions: or: y 0 eSuD: where the superscript g stands for given, and the double sign + is related to the surface orientation, locally defined by the normal nyo at each point y 0 .
The displacement discontinuity on the crack is obtained from [u(y 0 )] =u(yri)-u(yi))='P<yo) Using (3), one can calculate the stress intensity factors which are now directly related to the density cp.

Integral equations
It was shown in Reference 9 that, under the assumptions SeC 1 ·cz, O<a~ 1, <peCl.lf(S), 0<{3~ l, the integral equation for a general crack is expressed by Vy0 ES, t(y 0 , O,.ol = 16 7<(~ _ v 2 ) pv 1 (2(<1>. •• e, F. ,)n,0 -(1-2 v)H+ .•• e, n,.o) F., +(F,v·ny 0 )<P,u A e,)+3(e,-<P,u)((F,v, ny 0 , e,)e,+(n 10 ·e,)e, A F,v) -3(e, ·<P,v)((F,u, Dy 0 , e,)e,+(ny 0 ·e,)e, 1\ F,u)) dudv with now ( Figure l): r=lly-y 0 ll,e,=(y-yo)/IIY-Yoll F: a parametrization of S, defined on a domain A in IR 2 : t: the unit vector tangent to the crack edge aS, oriented according to the orientation of S. One should notice that the symbol pv before the two-dimensional integral relates to the principal value on S (not on A), and that the expression still holds in the general case when F.u is not Figure 1. General crack configuration with respectively a particular and an ordinary point, Yo and y perpendicular to F. v· The kernel of the line integral has been expressed for brevity as the product of a linear operator~ and the density q». Such line integrals appear when surface crack problems are involved.
In the next section, we shall concern ourselves with the discretization of the two-dimensional integral. As previously mentioned, this discretization will not include the study of the principal value.

Shape functions
The co-ordinates transformation is given by -r: e --+(u, v) = ( (N (e)){ u }e, (N(e) ){ v} e) (5) where e is a point of the reference element and the superscript e relates to an ordinary element. In this paper, we are interested to the most commonly encountered co-ordinate systems: the Cartesian, the polar and the cylindrical ones. In Cartesian co-ordinates (u, v)=(y 1 , y 2 ), in polar co-ordinates (u, v) = (p, 8) (Figure 2) and in cylindrical co-ordinates, with p constant, (u, v) = (8, y 3 ) ( Figure 3).

Interpolation of~
In practice, it is advantageous to obtain the components of the density~ in the local basis which varies with the pointy of S. For this reason, we attach much importance to the choice of the coordinates system and we will consider different ways to interpolate the components of q, in the local basis. Use will be made of isoparametric elements.
In Cartesian co-ordinates, for instance, the local basis b is invariant; it coincides with the global, --fixed basis e of the space E 3 :  The density q, is considered in the fixed basis: Each of the components of ~ is then interpolated classically by In cylindrical co-ordinates, the local basis varies with the point y: The density will be considered in the variable basis: cl » = 4> PeP+ l/Joe6 + ¢3 e3 and its components are now interpolated by iE{p, e, 3} Relations (6) in turn yield where c, denotes the canonical basis of IR". The letters b, c" added in small characters recall the bases under consideration. In the case of four-node elements, the matrix [ JV] can be explicited as fol1ows: N4 I I and the density t1l in the case of cylindrical co-ordinates can be expressed by One can then express the derivatives of t1l with respect to the parameters u and v: where The matrical representation of .# 1 in the bases b and c 12 (case of four-node elements) is

Discretization of the integral equation
(8a) (8b) (9) We study next how to transform each of the terms of the right-hand side of (4). First, let us consider the terms containing «Jl. u· We have (cp,u, e, F,Jny 0 =[cp,u·(er 1\ F,v)]ny 0 Similarly As for the term «11. u 1\ e, we have (11) where an implicit sum is implied on the repeated subscript iE{l, 2, 3} or {p, (}, 3}, b 0 =(b 0 J is the basis related to y 0 . It should be noted that the components of the given stress t(y 0 , ny 0 ) in the lefthand side of(4) are obtained in the basis b 0 . Furthermore, in interacting or surface crack problems, b 0 may differ or not from b, according to the relative position of the particular point y 0 and the ordinary one y. Now in virtue of (11): In the same manner, the fourh terms of (4) can be expressed by Eventually, relations (10) are expressed in matrical form by It should be noted that, in (12d) for instance, as the given stress in the left-hand side is expressed in the basis b 0 , the vector e,. between braces { } (column-matrix, see Appendix) must be expressed in the basis b 0 , whereas the same e,. between angle brackets ( ) (row-matrix) must be expressed in the basis b.
The terms containing «P. v are discretized in the same manner; they yield relations analogous to (12) where «P.v takes the place of <l»,u, and F,u of F,v· Thus, for each particular point y 0 of the crack, we have discretized the integral equation using a tensor formalism, equations (10). This equation being discretized as a vectorial equation, and not as three separate scalar equations, this formalism allows us a compact programmation for numerical purposes. One has to solve an algebraic system in the form

I. Embedded circular crack
First, the program has been tested on the simple case of a circular crack centred at 0, with radius a, under a uniform pressure p, the analytical solution of which is well known: 15 As for any plane crack imbedded in the infinite medium, the three modes are uncoupled. The stress intensity factor in the opening mode is given by The natural choice of parameters for the penny-shaped crack is (u, v)=(p, 0). One should notice that this choice yields a co-ordinates transformation which is well suited to the crack geometry: the crack edge is especially better shaped than with the usual grid pattern, as shown in Figure 4.
Moreover, this choice of parameters facilitates the computation of the stress intensity factors which require to be taken as a limit when moving along a direction perpendicular to the crack edge. The given stress vector is expressed at each point Yo in a local basis b 0 =(ep 0 ; e 80 ; e 3 ), whereas the vector density 4> at a pointy in b=(ep; e 0 ; e 3 ); this situation is shown in Figure 5. Figure 6 represents a quarter of a mesh with 181 interior nodal points (7 nodes along the radius and 36 along the polar angle (J, each sector being equal to 1 0°). The numerical integration was performed with 2 x 2 Gaussian points.
As expected, the numerical values of the radial and tangential displacement discontinuities q, P and ¢ 9 are zero; moreover, all the nodes of the same radius have identical ¢ 3 values. It is interesting to note that, in mixed modes problems, the formulation herein proposed permits us to obtain directly the radial and tangential components of cf» (l/Jp and c/> 9 ) without additional calculations, because cf» is expressed in the local polar basis. Table I and Figure 7 give respectively the values of c/> 3 normalized by the theoretical maximum value (c/> 3 max)theor=8(1-v 2 )pa/nE, and its curve compared to Sneddon's solution. It should be noticed that one obtains in fact   E/(16n(l-v 2 ))·c/> 3 directly from the program, so that the Young's modulus E datum is not necessary. The relative error on </> 3 is maximum at the crack centre and equals 8·75 per cent; that of the computed stress intensity factor K 1 is 3·97 per cent. Figure 8 represents a finer mesh with 529 interior nodes (13 along the radius and 48 along the polar angel 8, each sector being equal to 7·5°). Use is always made of 2 x 2 Gaussian points. Table II and Figure 7 give the values of ¢ 3 normalized by the theoretical maximum ¢rvalue. The relative error on ¢ 3 at the centre is 2·97 per cent; that of the computed stress intensity factor K 1 is 0· 25 per cent.

Semi-elliptical crack in a cylindrical bar under combined loads
Let us consider a cylindrical bar with radius R, containing a semi-elliptical crack with semi-axes a and b. The crack occurs perpendicularly in the lateral surface of the cylinder, along a line which will be referred to as the surface-line of the crack (Figure 9). Three types of loads will be applied uniquely on the crack, the external surface of the bar being left free: a uniform pressure, a linear pressure along the y 1 -axis, and a shear distribution corresponding to the torsion of the bar (Figure 1 0).
According to the discussion in Reference 9, the surface crack problem can be investigated by considering a two crack problem: the first crack is the studied one, the second is a cylindershaped crack high enough to simulate a free bar. We introduce the notations:    where the kernels in (l3a) and (13b), given by (4), are not made explicit because of their length; S~~f and S~~~ denote respectively the upper and lower part of the cylindrical crack, separated by the section containing Scr· The relation (13e) holds if we choose the normal to Scr upward with respect to the crack and the normal to Scyl outward with respect to the cylinder. With another choice of surfaces orientations, however, this relation still holds but the signs before each term may change. Furthermore, one can easily verify that, in various surface crack problems, one or more relations of type (13e) may be obtained. Relations of type (l3e) allow one to verify that, though every contour integral is different from zero, there always appears a set of thesewhich cancel one another.
For numerical purposes, a crack depth-radius ratio b/R=0·4 and a semi-axes ratio b/a= 5/1 ~0·714 are considered. The cylinder lateral surface is divided into 748 elements with 731 interior nodes, while the semi-elliptical crack contains but 58 elements and 60 interior nodes ( Figure 11). So large a node number on the free surface is due to the use of four-node elements with their sides along the parallel of the cylinder or along its axis. It should be wiser to use curved cylindrical elements (i.e. elements with distinct cylindrical nodal values) so as to realize a relatively refined mesh around the surface-line together with a coarse mesh far away, and to reduce notably the node number on the free surface. We have, however, not chosen such elements because of the difficulties caused to automatic mesh generations. In fact, the sole aim of this paper is to show the feasibility and the validity of the proposed approach. Anyway, the interesting feature is the little node number on the crack itself.
Like the previous numerical example, use is made of four-node elements and the numerical quadrature is performed with 2 x 2 Gaussian points. Tables III, IV and V give the 4>rvalues, ie{l, 2, 3}, resulting from three loading types. From Tables III and IV, one can note that, under uniform or linear pressure loadings, the crack is seen to be in the opening mode. Table V shows on the contrary that, under torsion loadings, the crack deforms in a mixed mode II+ III. With the same maximum stress values (i.e. u=umax• see notation, Figure 10), ¢ 3 -values due to uniform pressures are found to be greater than those due to linear pressure.
It is important to note that one obtains Cartesian components of q, over the crack, and cylindrical components on the cylinder. Some interesting results can be deduced from the <Pivalues, ie { p, e, 3}, over the free surface. In the opening mode for instance, the point belongs to the upper crack face and located at the centre of the ellipse has, besides a positive ¢ 3 , a negative ¢P.
This means that it moves upward with respect to the crack and outward with respect to the cylinder, as to some extent predicted by physical considerations. Figures 12, 13 and 14 represent the crack in non-deformed and deformed states, where for a better representation deformations are emphasized by a multiplier equal to 5.
It is rather attractive to represent the stress intensity factor variations along the crack edge, but in surface crack problems, the stress intensity factor definition is almost difficult to establish.  Figure 11. Surface crack mesh: 58 elements, 60 interior nodes. The elements are shaped into the elliptic co-ordinates defined on the ellipse Table III. Components of E/(16n(l-v 2 ))·cPcr resulting from a uniform pressure: ¢ 1 =¢ 2 =0, ¢ 3 #:0; (v =0·3, b/R =0·4, bja = 5/7, b= 5, u= 1, the unit of a, b, R is [L], that of u is [F· L-2 ], that of the above Node 16n(l-v 2 ) ·¢ 3 Node 16n(l-v 2 ) ·¢ 3 Node 16tt(l-v 2 )·¢ 3 Node 16n(l-v 2 )·c/> 3  Table V. Components of E/(16n(l-v 2 ))·_,cr resulting from a torsion shear distribution with tmax=l:c/>3=0, l/J1, c/>2:-#:0 Generally speaking, the elastic state varies from plane strain in the interior to plane stress at the free surface. In the following, we shall limit our consideration to the stress intensity factors at the The stress intensity factors are normalized by u max fo or -rmax fo related to a slit crack with length 2b in the infinite medium, under a remote tensile or shear loading. In Reference 16, stress intensity factors are computed by means of the boundary integral equation method for semi-elliptical cracks in cylindrical bars of diameter equal to 12 mm, in tension or in bending. The crack geometry corresponding to that considered herein is bfa =2·4 mm/l36 mm=5/7, b/2R=2·4 mm/ 12 mm=0·2, for which the stress intensity factor can be obtained by interpolation: in tensile loading with a= 1000 MPa, K, = 47·10 MPajffi --.KJ!(amaxfo)=0·569 in pure bending with a max= 953·4 M Pa, · One can find in Reference 1 7 a list of stress intensity factors for various crack and solid geometries. The so-called surface crack in solid cylinder is studied under tensile loading and pure bending ( Figure 15). These actually are almond-shaped cracks which occur perpendicularly to the lateral surface of the cylinder. For lack of matter for comparisons, we consider, however, the crack of this type with b/ R =0·4; at the deepest point of the crack, the computed values are expected to be close to the above ones. The stress intensity factor is then expressed in the additive form 1 7 which yields in tensile loading: G = 0·92(2/n)· 1/cos fJ· [tan fJ/fJ] 1  Recently, in Reference 18, use is made of the finite element method to study semi-elliptical cracks in cylindrical bars under tensile loadings. When computing the stress intensity factor K 1 , the assumption that a plane strain state holds in the neighbourhood of the crack is made in order to apply the virtual crack extension method. Charts for K 1 at the centre of the crack are given for b/a=O·Ol, 0·2, 0·5, 1·0, 2·0 and for various bfR ratios. An extrapolation from these K 1 -charts then gives us approximately the interval containing the K 1 /((1J'i{h) ratio: for b/ R =0·4, bfa= 5/7, with a uniform pressure (1, The discrepancy between the above results is mainly due to the crack geometry differences; nevertheless, the values obtained are of the same order. For precise comparisons to be made, ·~ 2R. Figure 15. Surface crack in a solid cylinder (figure from Reference 17) investigations are being pursued to obtain more complete abacuses of stress intensity factors, reported for various crack configurations.

CONCLUDING REMARKS
Equation (4) has been discretized and written in a tensor form, equation (10). Whereas lengthy expressions should be expected if (4) had been dealt with as three separate scalar equations, this approach yields a concise form of the discretized equation and thus allows a compact way of programming. Moreover, a crack parametrization being chosen; this enables one to obtain directly the components ofthe displacement discontinuity vector at each point, as well as the stress intensity factors along the crack edge, in a local basis. Unlike the usual finite element method, no singular element is required near the crack edge. The method herein proposed can also easily be applied to non-symmetrical crack configuration problems for which the boundary integral equation method becomes inefficient. Moreover, since no a priori assumptions of a plain strain or plain stress are to be made, the integral equation method can be used to investigate the stress singularity in the neighbourhood of the crack.
The example on the circular crack shows that one can obtain good results, using only four-node elements and 2 x 2 Gaussian points. The problem of the semi-elliptical crack in a cylindrical bar was also investigated with combined loadings. The numerical results obtained in opening mode were found to be in the range of different existing results for analogous crack geometries. The proposed method also allows one to obtain results for more complex loadings, in particular for torsion loading to which few works in the literature are devoted. APPENDIX Notation b=(b 1 ; b 2 ; b 3 ) a local basis at a pointy in E 3 c,. the canonical basis of ~~~ f} the elastic body characterized by constants (E, v) of} the boundary of q) e=(e 1 ; e 2 ; e 3 ) the fixed basis of E 3 E 3 the three-dimensional Euclidean space F a parametrization of S, defined on a domain A in IR 2 : A3(u, v)-+ y = F(u, v)ES. In general, F,u may not be perpendicular to F.v I unit tensor ny normal to S at yeS S the crack surface as the crack edge t(y, n,) stress vector at yeS, with respect to normal.n, T Kupradze tensor y, y 0 respectively an ordinary, a particular point of S pvJ., surface integral on S, defined in the principal value sense {a} column-matrix of components of a vector a in a basis. One may encounter b{ a} which emphasizes the basis bin question (a) row-matrix, = {a} T, transpose of {a} [%] matrical representation of a linear transformation .AI from E 3 into IRn, also Cn denoted more precisely by e [AI'] =(a " b)·c (a, b, c) a®b tensor product defined by Vc E IR 3 , (a® b)c = a(b·c)