- Email: [email protected]
Dynamic analysis of cracks using boundary element method
Engineering Fracture Mechanics Vol. 34, No. 5/a, pp. 1051-1061,1989
Printed in Great Britain.
0013.7944189$3.00+ 0.00 0 1989Pergamon Press plc.
DYNAMIC ANALYSIS OF CRACKS USING BOUNDARY ELEMENT METHOD FRANCISCO CHIRINO and JOSE DOMINGUEZ Escuela Superior de Ingenieros Industriales, Universidad de Sevilla, Av. Reina Mercedes s/n, 41012~Sevilla, Spain Abstract-This paper presents a procedure for dynamic stress intensity factor computations using traction singular quarter-point boundary elements. Cracks in a complete space, a half-space and a finite body loaded by steady state waves are studied. Curves for elastodynamic stress intensity factors vs frequency are presented. Transient stress intensity factors are computed by means of Fourier transform. The results are compared with other authors and shown to be accurate in all cases. The dynamic stress intensity factors are computed in a very direct and easy way to implement, This versatile procedure allows for the study of problems with complex geometry that include one or several cracks.
INTRODUCTION ~NSI~E~LE attention has been given in recent years to dynamic fracture mechanics problems. The dynamic stress intensity factor plays a very important role in these problems. In his early work, De Hoop[l] presented the dynamic analysis of a semi-infinite crack in an infinite solid subjected to a step pressure. Sih and Loeber[2], and Ma1[3] computed dynamic stress intensity factors of a finite crack in an infinite plane under plane harmonic waves. Ma1[4,5] studied the difFraction of axisymmetric harmonic waves by a circular crack in an infinite medium. Transient dynamic stress intensity factors were obtained by Thau and Lu[G] for a finite crack in an infinite plane and at times prior to the arrival at the crack tip of waves rediffracted by the other tip. Sih, Embley and Ravera[7] studied the same problem as Tau and Lu applied to longer periods of time. Transient dynamic stress intensity factors for the half-plane crack in an unbounded solid were computed by Freund[8] by superposition over a fundamental solution obtained from the theory of dislocations. Dynamic stress intensity factors for cracks in a half-plane subjected to time harmonic excitation were obtained by Achenbach et af.[9], Stone et ul.[lO], Keer et a/.[1 I], Lin et ul.[12], Shah et a1.[13], and Tittmann er a1.[14], and for 3D problems by Angel and Achenbach[l5]. Van der Hijden and Neerhof@6] studied the same problem of Sih and Loeber[2] using a different computation procedure that permits an increase in the range of frequencies of the analysis. The aforementioned procedures for dynamic stress intensity factor computations are restricted to infinite or semi-infinite domains and include a numerical evaluation of integrals or a solution of integral equations. When the domain is finite the dynamic crack problem has been tackled by discreti~tion of the body using a numerical method. Chen[ 171studied the case of a centre cracked plate subjected to a step function load using finite differences. Finite elements have been used by several authors; for instance, Nishioka and Atluri[ 181,Atluri and Nishioka[ 191,Aoki et a1.[20], and Murti and Valliapan[21]. In the present paper the Boundary Element Method is used for Dynamic stress intensity factor computations using a singular quarter-point quadratic element. The method is applied to cracks in the complete plane, the half-plane and finite two-dimensional domains under harmonic load. Also, a Fourier transform algorithm is used to determine transient dynamic stress intensity factors. The Boundary Element Method has been applied to static stress intensity factor compu~tions in various ways[22-241, in particular, a singular quarter-point boundary element was used by Blandford et a1.[25], and by Martinez and Dominguez[26], This procedure is now applied to dynamic loading. The domain is always divided into subdomains by means of cuts along the crack in such a way that the two crack faces belong to the boundaries of different regions. All boundaries are discretized into elements; boundary conditions applied to external boundaries, and compatibility and equilibrium conditions to internal ones. The multidomain discretization avoids the numerical problems derived from having two displacement variables for each point along the crack. 1051
FRANCISCO CHIRINO and JOSE DOMINGUEZ
1052
Details of the fo~ulation and numerical implemen~tion of the Boundary Element Method may be seen, for instance, in the book of Brebbia and Dominguez[27& SINGUT.,AR ELEMENT
FOR STRESS
INTENSITY
FACTOR
COMPUTATIONS
The Boundary Element Method formulation for elastodynamic harmonic problems leads to the same integral equation of the static formulation, the only difference being that the variables are frequency dependent. The fundamental solution, corresponding to a harmonic point load in the unbounded plane, is a function of modified Bessel functions of the second kind (see, e.g. Dominguez and Alarcon~28}). The singular quarter-point boundary element may be formulated in a similar way to the static case. The geometry, displa~ment and traction aiong a quadratic element for two-dimensional domains may be represented as f;=CdX+6?ff+#Y
(1)
where 4 ‘, $2 and 4’ are well known quadratic shape functions written in terms of the natural coordinate 4 (Fig. l);A represents a Cartesian component of the displacement. traction or geometry along the element, and ff is the value of that variable at node j. When the quadratic element has a straight-line geometry and the mid-node is placed at a quarter of the length of the element, a simple relation between 5 and the variable Palong the element (Fig. 1) exists. The equations (1) may be written in terms of J as f;=
at+a:
J P
5+“;:
-
(2)
where
(3) This kind of element is usually called a q~rter-point element and is able to represent the ,/r behaviour of the displacement near the tip when one makes P coincide with r along one of the crack faces (Fig. 2). The singularity of the stress, and consequently of the traction, near the tip of the crack may be included in the representation of the traction by using modified quarter-point elements with singular shape functions.
where
ij are now the nodal vahtes of ‘3 t, =
3. ti,
ti divided by the nodal values of p, i.e.: if =
t’/2;
it =I$
tj
Jq
and (4) written in terms of i is now ti=a;llTi;-+a:+(i:Ji-/i
Fig. 1. Quadratic elements.
Fig. 2. Coordinates near the tip.
Dynamic analysis of cracks
1053
where a!=i,!;
$=
-i:+4i:-3i)
and
$=2i!-4if.
The singular quarter-point boundary element includes a representation of the traction by means of (6) and a representation of the displacement by means of (2). Thus, both displacements and tractions may be represented, within the vicinity of the crack tip by the first three terms of their series expansion. The first and second mode stress intensity factors can be defined by the following limits (Fig. 2). K, = lim (2n~,)‘/~a 22 x, -0
K,, = limo(27tx,)“2a,2.
(7)
If boundary discretization is carried out in such a way that the first element to cut the domain from the crack tip follows 8 = 0” and this element is a singular quarter-point boundary element, then for this element, P = xi, C,= CJ,2, f2 = oz and the nodal values for the tractions at the tip node K are: if = lim t:(p/Z)“2 = lim o,~(x,//)‘/~ i-0
if: = fi
x, -0
tf(P/l)‘12 = liyo 022(x, /l)‘12.
(8)
Thus, the stress intensity factors coincide with the traction nodal values except for a constant K, = i$(2d)i’2 K,, = i$(2d)i’2
(9)
and may be computed directly with the boundary element code. Since the above procedure makes use of the correct representation of the tractions and only depends on the nodal value at the tip, little dependence upon local boundary discretization may be expected. Blandford et a/.[251 used the singular quarter-point boundary element to compute stress intensity factors in elastostatics by means of a correlation formula of the displacements along the crack surfaces that has been previously used in finite elements[29]. They studied the effects of the size of the elements near the tip and of the use of transition elements. Martinez and Dominguez[26] showed how the use of the traction nodal values of the singular element at the crack tip (9) is substantially less sensitive to the discretization than any of the displacement correlation procedures. When working in steady state elastodynamics, it is even more important to have a procedure with very little dependence on the size of the elements near the tip. This is because it would be difficult to establish rules of discretization that guarantee accurate values of the stress intensity factors if this dependence is important and changes with frequency as can be expected. In the following, the singular quarter-point boundary element in combination with eq. (9) is used for dynamic stress intensity factor computations. DIFFRACTION OF WAVES BY A CRACK IN THE COMPLETE PLANE The problem of harmonic elastic waves that propagate in the complete x,x,-plane and are diffracted at a crack whose length is 2~2,is analysed. The direction of propagation of the waves form an angle y with the crack (Fig. 3). The stress and displacement fields due to P or SV plane waves in a two-dimensional domain can be represented in terms of two scalar potentials: --iF(x,cosy
4 =fjoexp 1 *=o EFM 3+5/6-D
P
+x,siny
+C,t) 1
(10)
FRANCISCO CHIRINO and JOSE DOMINGUEZ
1054
Incident
Waves
Fig. 3. Waves impinging on a crack.
for P-waves, and
-iF(x,cosy +x,siny +C,t) (11) s 1 I for SV-waves; where y is the angle that the wave makes with the x,-axis; C, and C, are the P and S-waves velocity, respectively; o is the frequency and &, Il/,,are amplitudes. Displacements and stresses are easily obtained by differentiation of these potentials. The analysis of wave diffraction by a crack in an infinite medium is carried out by the superposition of two problems. One, the incident field in the domain without crack; and the other, the cracked domain loaded on the crack faces by tractions equal and opposite to those that appear in the untracked plane along the lines where the crack faces are in the cracked plane (scattered field): $ =&exp
uf = uy.+ ui;
al,=a?+a,i
(12)
where ~7 and gi are the displacement and stress components for the incident field, while U, and rrgcorrespond to the scattered field. Only the second problem has to be solved because the incident field in the complete plane is known. Since there are not infinite values of the stress in the untracked plane, the stress intensity factors of the original diffraction problem are the same as in the second problem (scattered field). The Boundary Element Method formulation for elastodynamics is based on the integral representation of the boundary points displacement which, once the boundary conditions are applied, becomes an integral equation that is solved numerically by discretization of the boundary into elements (see, e.g. Dominguez and Alarcon[28]). When the problem at hand satisfies the radiation conditions the displacement representation is written for the external region in terms of integrals that only extend over the internal boundaries (see, e.g. Eringen and Suhubi[30]). This fact is the basis of an important difference between domain methods (Finite Elements, Finite Differences) and the Boundary Element Method where unbounded domains are represented by elements only on the internal boundaries and no problem exists because of spurious reflections on external boundaries. To avoid the problem of having two different values of the displacement for points on the faces of the crack, a boundary that divides the domain into two parts along the x,-axis is introduced. This boundary extends to infinity. However, the discretization is truncated at a distance to the crack tip equal to fifteen times the half-length of the crack. This can be done because elements that are far from the crack have very little effect on the solution near the tip since both the scattered field and the fundamental point load solutions satisfy the regularity conditions. Figure 4 shows the
a
1
15a 1
Fig. 4. Boundary discretization.
Dynamic analysis of cracks
1055
number of elements used in one half of the boundary. The problem is always decomposed into its symmetric and skewsymmetric parts. The two elements to which the node at the crack tip belongs are singular quarter-point elements, whereas all the others are regular quadratic elements. Stress intensity factors for P and SV waves with several angles of incidence have been computed. The values of Ki and K,, depend on the incident angle y and the Poisson’s ratio of the elastic medium. Figure 5 shows the moduli of the stress intensity factors for P-waves normalized with respect to their corresponding static values. Results are plotted vs the dimensionless frequency oa/Cp for a Poisson’s ratio 0.25. Plane strain is assumed. It can be seen in the figure how K, is maximum for y = 90” and 4, for y = 45”. In all cases the values tend to rise first, reach a maximum, and then fall. The dynamic stress intensity factor can be as much as 1.3 times the static one. The values computed using the BEM with singular quarter-point elements are compared with those obtained by Chen and Sih[31] using a method derived by Copson[32] to solve a system of integral equations. As can be seen in the figure, the agreement is very good. The values of $ and K,, for W-waves are shown in Fig. 6. As could be expected the maximum value of K, is not for y = 45” and the maximum of K,, for y = 90”. The agreement with Chen and Sih[31] is again very good.
I I
*
0.1
t
0.2
0.3
0.4
0.5
0.6
Chen
v = 0.25
watcp
0.7
& Sih B.E.M.
.
0.27
+
I
0.1
0.2
0.3
0.4
0.5
e
I
0.6
0.7
u PICp
Fig. 5. Ratio of dynamic and static stress intensity factors for incident P-wave.
FRANCISCO CHIRINO and JOSE DGMINGUEZ
1056
Chen
8, Sih B.E.M.
A
I
0.1
1
a
1
8
1
,
I
0.2
0.3
0.4
0.5
0.6
0.7 oalcp
0.6
07 wa/cp
1. 2
1.0
0.8
Chen
& Sih
-
Y
y=60°
0.6
0.1
0.2
0.3
0.4
0.5
Fig. 6. Ratio of dynamic and static stress intensity factors for incident W-wave.
CRACK
NEAR A FREE SURFACE
In recent years, Achenbach and co-workers[9,33,11,12] have studied the scattering of harmonic waves by edge or sub-surface cracks. Their formulation leads to a system of integral equations which is solved numerically. The fundamental solution of these integral equations corresponds to dislocations in the half-plane and is calculated numerically as a combination of Hankel functions and contour integrals. In the following, stress intensity factors are computed for a horizontal subsurface crack subjected to harmonic excitation. The problem has been taken from Keer et aZ.[l l] for comparison. The regular frequency domain dynamic formulation of the BEM with the complete space fundamental solution is used. Thus, an integral equation has to be solved numerically, but, the kernels of the integral are much simpler and easier to compute than those of the dislocations in the half-plane. On the other hand, the free surface boundary conditions are not satisfied by the fundamental solution and have to be enforced by discretization of the free surface into boundary elements for which zero tractions are prescribed. However, the elements far from the crack have very little effect on the stress field near the crack and accurate representation of the effect of the free surface can be obtained by discretizing a limited zone close to the crack. Due to the nature of the Boundary Element Method, the fact of extending the discretization only to a finite part of the free surface does not introduce undesired artificial reflections. Figure 7 shows the boundary element discretization used for the problem at hand. The free surface discretization and the internal boundary extend to a distance from the crack of 15~ along
1057
Dynamic analysis of cracks
Fig. 7. Boundary element discretization for subsurface crack.
the x, axis. All the elements are quadratic and those at the crack tips are quarter-point singular. The incident field is a uniform tension applied on x2 = 0 which results in uniform tension on the crack faces. The Poisson’s ratio is 0.3. Figure 8 shows the values of K, (magnitude of the ratio of the dynamic and static mode I stress intensity factor) vs the parameters d/a and ad/C,, with C, as the Rayleigh wave velocity. The results are in good agreement with those taken from the figures of Keer et aZ.[ll]. The same comparison is done in Fig. 9 for the mode II stress intensity factor. CENTRAL
CRACK IN A RECTANGULAR
PLATE
The following example has been taken from Chen and Wilkins[34] and shows how the singular quarter-point boundary element can be applied to bounded domain problems. Figure 10 presents the geometry and boundary conditions of the cracked plate. The mechanical properties are: Shear modulus p = 76.923 GPa, Poisson’s ratio v = 0.3; and density p = 5000 kg/m’. In addition to the purely elastic material, internal damping is considered by means of a complex shear modulus c~f= ~(1 + 2pi). Values of j? = 0.01, fi = 0.025 and /I = 0.05 are assumed. The load is a uniform tension applied on two opposite sides. Figure 11 shows one quarter of the plate with the boundary divided into three different zones: zone A (a + L, = 2a), zone B (L,) and zone C (L, + L., + L5). The boundary discretization was made according to the following rules: The length of the elements should be I< J./l0 in zone A, 1 < I.15 in zone B, and I < 113 in zone C, with rZbeing the length of the S-waves in the plate material. The magnitude of the ratio of the dynamic and static mode I stress intensity factor for three different values of the internal damping, is shown vs frequency in Fig. 12. Several peaks of resonance can be seen in the figure. The peaks are shifted and damped with increasing values of /J. The first peak is 20% below the first natural frequency of the untracked plate. Since the only 13 d/a = 04
12 II
B.E.M. 0.E.M. KEER
/VP
-
et al.
KEER
et aI.-----
-----
d/a=2
Fig. 8. Ratio of dynamic and static mode I stress intensity factor for subsurface crack and vertically incident P-wave.
Fig. 9. Ratio of dynamic and static mode II stress intensity factor for subsurface crack and vertically incident W-wave.
1058
FRANCISCO CHIRINO and JOSE DOMINGUEZ
v
=
P
= 76.923
V.JV
GPa
P= 5 Ton.lm3 P=
0.4
l-4
1 1
GPa
Fig. 11. Model for centre cracked plate.
Fig. 10. Centre cracked plate.
results known for the cracked plate problem correspond to transient loads, comparison will be done in the next section cone the above results are transformed into the time domain. TRANSIENT
STRESS
INTENSITY
FACTORS
Going back to the problem of diffraction of waves by a crack in an infinite medium (Fig. 3), a normal incident P-wave with a step function time variation is considered. Results for transient stress intensity factors may be obtained by using the frequency analysis of the problem carried out previously and the FFT algorithm. Thau and Lu[6] computed values of the stress intensity factors for this problem using a generalized Weiner-Hopf technique, but they were only able to study the problem until the time when a diffracted P-wave on reaching the opposite edge is rediffracted and comes back to the original edge. A comparison of the boundary element results with Thau and Lu[6] may be seen in Fig. 13 where Ki is shown vs dimensionless time. The agreement is again good, and K, reaches a value that is 30% higher than the static. Sih et al.[7j studied the same problem
DAMPING 0%
-
ua/cs .b
0.6
0.6
1
1.2
Fig. 12. Ratio of dynamic and static mode I strw intensity factor for centre cracked plate under uniform harmonic traction.
Dynamic analysis of cracks
1059
y” 0.6 .. 0.E.M.
-
THAU- LU -----
Fig. 13. Ratio of transient and static mode I stress intensity factor for crack in the complete plane under a step pressure.
using integral transforms coupled with the technique of Cagniard. They could obtain values of K, for longer periods of time, but because of numerical errors they could not represent the discontinuity in the slope, shown by the BEM and by Thau and Lu, which is a consequence of the arrival of the Rayleigh waves diffracted by the other tip. A comparison of the boundary element results with Sih et a1.[7j is shown in Fig. 14. The last problem to be analysed corresponds to the transient response of the centrally cracked plate that has already been studied in the frequency domain (Fig. 10). A uniform traction P on two opposite sides that varies with time as a Heaviside step function is considered. Again the FFT algorith is used to compute the transient response. The magnitude of the mode I stress intensity factor normalized by P(t)r au is represented for the first 12 ps in Fig. 15. Values for damping factors of l%, 2.5% and 5% are plotted. In this figure, results are compared with those obtained by Chen[l7] using finite differences, and with Baker’s[35] of the semi-infinite crack under pressure in an infinite plate. The latter are valid for times prior to the arrival of the first wave diffracted by the other tip (R,). Both consider purely elastic material. It can be said that the results am in good agreement with Chen’s and Baker’s and that, as could be expected, the existence of internal damping reduces the maximum value reached by the stress intensity factor. This reduction becomes more important as the internal damping increases. Results for the purely elastic material are not shown because they present spurious oscillations that can be explained by numerical errors in the Fourier transformation process. The use of filters could solve the problem, but, since the internal damping makes the spurious oscillations disappear, it is easier to introduce a small amount of damping in the material.
4
Fig. 14. Ratio of transient and static mode I stress intensity factor for crack in the complete plane under a step pressure.
FRANCISCO CHIRINO and JOSE ~MINGUEZ
1060
1.6 Cl.2
0.8 0.1 0
-0.5 2
L
6
a
Fig. 15. Ratio of transient mode I stress intensity factor and P(t)& uniform traction.
10
12 t)ts
for centre cracked plate under step
CONCLUSIONS A procedure for dynamic stress intensity factor computations using the Boundary Element Method in the frequency domain has been proposed. The stress intensity factors are computed from the nodal values of a traction singular quarter-point element located at the crack tip. The length of the singular and adjacent elements should be smaller than L/10, with 3, as the length of the S-waves, and also smaller than 0.6~2,with LZbeing the half-length of the crack. Elements far from the crack should not be longer than L/3. It has been shown that this direct and very simple procedure, previously used by one of the authors in eiastostatics, produces accurate results for cracks in the complete space, in the half-space, or in finite bodies. The method is very versatile, it may be easily included in any quadratic boundary element program, and can be used for many different geometries, including problems with several cracks. The study has been done in the frequency domain and it has been shown how transient stress intensity factors can be computed using the Fourier Transform. A study of the use of the singular quarter-point boundary element with the time domain formulation of the BEM for transient stress intensity factor computations, will be presented in a forthcoming paper. ,4cknowZegemenrs-The authors would like to express their gratitude to the Spanish Comision Interministerial de Ciencia y Tecnologia for supporting this work under the research grant No. PBS6-0139.
REFERENCES [I] De Hoop, Representation theorems for the displacement in an elastic solid and their applications to elastodynamics diffraction theory. Thesis, Technische Hogeschool, Delft (1959). [2] G. C. Sih and J. F. Loeber, Wave propagations in an elastic solids with a line of discontinuity or finite crack. Q. uppl. Murk 27, 193-213 (1969). [3] A. K. Mal, Interaction of elastic waves with a Griffith crack. Znt. J. Engng Sci. 8, 763-776 (1970). [4] A. K. Mal, Diffraction of elastic waves by a penny-shaped crack. Q. appl. Math. 26, 231-238 (1968). [5] A. K. Mal, Interaction of elastic waves with a penny-shaped crack. Znt. J. Z%rgngSci. g, 381-388 (1970). [6] S. A. Thau and T. H. Lu, Transient stress intensity factors for a finite crack in an elastic solid caused by a dilatational wave. Znt. J. Solids Struct. 7, 731-750 (1971). [7] G. C. Sih, G. T. Embley and R. S. Ravera, Impact response of a finite crack in plane extension. Znt. J. Soliak Struct. 8, 977-993
(1972).
[S) L. B. Freund, The stress intensity factor due to normal impact loading of the faces of a crack. ht. J. Engng Sci. 12, 179-189 (1974). 191 J. D. Achenbach, L. M. Keer and D. A. Mendel~hn, Elastodynamic analysis of an edge crack. J. uppl. Mech. 47, 551-556 (1980). (IO] S. F. Stone, M. L. Ghosh and A. K. Mal, Diffraction of antiplane shear waves by an edge crack. J. uppl. Mech. 47, 359-362 (1980). [I l] L, M. Keer, W. Lin and J. D. Achenbach, Resonance effects for a crack near a free surface. J. uppl. Mech. 51,65-70 (1984).
Dynamic analysis of cracks
1061
[12] W. Lin, L. M. Keer and J. D. Achenbach, Dynamic stress intensity factors for an inclined subsurface crack. J. uppl. Mech., 51, 773-779 (1984). [13] A. H. Shah, K. C. Wong and S. K. Datta, Surface displacements due to elastic wave scattering by buried planar and non-planar cracks. Waue Motion 7, 319-333 (1980). [14] B. R. Tittmann, L. A. Ahlberg and A. K. Mal, Rayleigh wave diffraction from surface-breaking discontinuities. Appl. Phys. Lett. 20, 1333-1335 (1986). [15] Y. C. Angel and J. D. Achenbach, Stress intensity factor for 3-D dynamic loading of a cracked halfspace. J. Efasriciry 15, 89-102 (1985). [16] J. H. M. T. Van Der Hijden and F. L. Neerhoff, Scattering of elastic waves by a plane crack of finite width. J. appf. Mech. 51, 646-651 (1984). [17] Y. M. Chen, Numerical computation of dynamic stress intensity factor by Lagrangian finite-difference method. Engng Fracture Me&. 7, 653-660 (1975). (181 T. Nishioka and S. N. Atluri, A numerical study of the use of path independent integrals in elastodynamic crack propagation. Engng Fracture Me&. 18, 23-33 (1983). [19] S. N. Atluri and T. Nishioka, Numerical studies in dynamic fracture mechanics. Znr. J. Fracture 27, 245-261 (1985). [20] S. Aoki, K. Kishimoto, H. Kondo and M. Sakata, Elastodynamic analysis of crack by finite element method using singular element. Inr. J. Fracrure 14, 59-68 (1978). [21] V. Murti and S. Valliappan, The use of quarter point element in dynamic crack analysis. Engng Fracture Mech. 23, 585-614 (1986). [22] M. D. Snyder and T. A. Cruse, Boundary integral equation analysis of anisotropic cracked plates. Inr. J. Frucrure 11,
315-328 (1975). [23] M. Stem, E. B. Becker and R. S. Dunham, A contour integral computation of mixed-mode stress intensity factors. Inr. J. Fracture 12, 359-368 (1976). [24] T. A. Cruse, Two-dimensional BIE fracture mechanics analysis. Appl. Math. Modelring 2, 287-293 (1978). [25] G. E. Blandford, A. R. Ingraffea and J. A. Liggett, Two-dimensional stress intensity factor computations using the boundary element method. Inr. J. Numer. Merh. Engng 17, 387-404 (1981). [26] J. Martinez and J. Dominguez, On the use of quarter-point boundary elements for stress intensity factor computations. Inr. J. Numer. Merh. Engng 20, 1941-1950 (1984). [271 C. A. Brebbia and J. Dominguez, Boundary Elements. An Inrroductory Course. CMP Publications McGraw-Hill, New
York (1989). [28] J. Dominguez and E. Alarcon, Elastodynamics, in Progress in Boundary Efernenr Merhocis (Edited by C. A. Brebbia). pp. 213-256. Pentech Press, Plymouth (1981). [29] C. F. Shih, H. G. de Lorenzi and M. D. German, Crack extension modeling with singular quadratic isoparametric elements. Znr. J. Fracture 12, 647651 (1976). [30] A. C. Eringen and E. S. Suhubi, Eiastodynamics. Academic Press, New York (1975). [31] E. P. Chen and G. C. Sih, Scattering waves about stationary and moving cracks, in Mechanics of Fracture: Efustodynumic Crack Problems (Edited by G. C. Sih), pp. 119-212. Noordhoff, Leyden (1977). [32] E. T. Copson, On certain dual integral equations. Proc. Glasgow Math. Assoc., Vol. 5, pp. 19-24 (1961). [33] J. D. Achenbach and R. J. Brind, Scattering of surface waves by a sub-surface crack. J. Sound Vibration 76, 43-56 (1981). [34] Y. M. Chen and M. L. Wilkins, Numerical analysis of dynamic crack problems, in Mechanics Elasrodynamic Crack Problems (Edited by G. C. Sih), pp. 295-345. Noordhoff, Leyden (1977). [35] B. R. Baker, Dynamic stresses created by a moving crack. J. appl. Mech. 29, 449-545 (1962). (Received 7 February 1989)
of Fracture:
Printed in Great Britain.
0013.7944189$3.00+ 0.00 0 1989Pergamon Press plc.
DYNAMIC ANALYSIS OF CRACKS USING BOUNDARY ELEMENT METHOD FRANCISCO CHIRINO and JOSE DOMINGUEZ Escuela Superior de Ingenieros Industriales, Universidad de Sevilla, Av. Reina Mercedes s/n, 41012~Sevilla, Spain Abstract-This paper presents a procedure for dynamic stress intensity factor computations using traction singular quarter-point boundary elements. Cracks in a complete space, a half-space and a finite body loaded by steady state waves are studied. Curves for elastodynamic stress intensity factors vs frequency are presented. Transient stress intensity factors are computed by means of Fourier transform. The results are compared with other authors and shown to be accurate in all cases. The dynamic stress intensity factors are computed in a very direct and easy way to implement, This versatile procedure allows for the study of problems with complex geometry that include one or several cracks.
INTRODUCTION ~NSI~E~LE attention has been given in recent years to dynamic fracture mechanics problems. The dynamic stress intensity factor plays a very important role in these problems. In his early work, De Hoop[l] presented the dynamic analysis of a semi-infinite crack in an infinite solid subjected to a step pressure. Sih and Loeber[2], and Ma1[3] computed dynamic stress intensity factors of a finite crack in an infinite plane under plane harmonic waves. Ma1[4,5] studied the difFraction of axisymmetric harmonic waves by a circular crack in an infinite medium. Transient dynamic stress intensity factors were obtained by Thau and Lu[G] for a finite crack in an infinite plane and at times prior to the arrival at the crack tip of waves rediffracted by the other tip. Sih, Embley and Ravera[7] studied the same problem as Tau and Lu applied to longer periods of time. Transient dynamic stress intensity factors for the half-plane crack in an unbounded solid were computed by Freund[8] by superposition over a fundamental solution obtained from the theory of dislocations. Dynamic stress intensity factors for cracks in a half-plane subjected to time harmonic excitation were obtained by Achenbach et af.[9], Stone et ul.[lO], Keer et a/.[1 I], Lin et ul.[12], Shah et a1.[13], and Tittmann er a1.[14], and for 3D problems by Angel and Achenbach[l5]. Van der Hijden and Neerhof@6] studied the same problem of Sih and Loeber[2] using a different computation procedure that permits an increase in the range of frequencies of the analysis. The aforementioned procedures for dynamic stress intensity factor computations are restricted to infinite or semi-infinite domains and include a numerical evaluation of integrals or a solution of integral equations. When the domain is finite the dynamic crack problem has been tackled by discreti~tion of the body using a numerical method. Chen[ 171studied the case of a centre cracked plate subjected to a step function load using finite differences. Finite elements have been used by several authors; for instance, Nishioka and Atluri[ 181,Atluri and Nishioka[ 191,Aoki et a1.[20], and Murti and Valliapan[21]. In the present paper the Boundary Element Method is used for Dynamic stress intensity factor computations using a singular quarter-point quadratic element. The method is applied to cracks in the complete plane, the half-plane and finite two-dimensional domains under harmonic load. Also, a Fourier transform algorithm is used to determine transient dynamic stress intensity factors. The Boundary Element Method has been applied to static stress intensity factor compu~tions in various ways[22-241, in particular, a singular quarter-point boundary element was used by Blandford et a1.[25], and by Martinez and Dominguez[26], This procedure is now applied to dynamic loading. The domain is always divided into subdomains by means of cuts along the crack in such a way that the two crack faces belong to the boundaries of different regions. All boundaries are discretized into elements; boundary conditions applied to external boundaries, and compatibility and equilibrium conditions to internal ones. The multidomain discretization avoids the numerical problems derived from having two displacement variables for each point along the crack. 1051
FRANCISCO CHIRINO and JOSE DOMINGUEZ
1052
Details of the fo~ulation and numerical implemen~tion of the Boundary Element Method may be seen, for instance, in the book of Brebbia and Dominguez[27& SINGUT.,AR ELEMENT
FOR STRESS
INTENSITY
FACTOR
COMPUTATIONS
The Boundary Element Method formulation for elastodynamic harmonic problems leads to the same integral equation of the static formulation, the only difference being that the variables are frequency dependent. The fundamental solution, corresponding to a harmonic point load in the unbounded plane, is a function of modified Bessel functions of the second kind (see, e.g. Dominguez and Alarcon~28}). The singular quarter-point boundary element may be formulated in a similar way to the static case. The geometry, displa~ment and traction aiong a quadratic element for two-dimensional domains may be represented as f;=CdX+6?ff+#Y
(1)
where 4 ‘, $2 and 4’ are well known quadratic shape functions written in terms of the natural coordinate 4 (Fig. l);A represents a Cartesian component of the displacement. traction or geometry along the element, and ff is the value of that variable at node j. When the quadratic element has a straight-line geometry and the mid-node is placed at a quarter of the length of the element, a simple relation between 5 and the variable Palong the element (Fig. 1) exists. The equations (1) may be written in terms of J as f;=
at+a:
J P
5+“;:
-
(2)
where
(3) This kind of element is usually called a q~rter-point element and is able to represent the ,/r behaviour of the displacement near the tip when one makes P coincide with r along one of the crack faces (Fig. 2). The singularity of the stress, and consequently of the traction, near the tip of the crack may be included in the representation of the traction by using modified quarter-point elements with singular shape functions.
where
ij are now the nodal vahtes of ‘3 t, =
3. ti,
ti divided by the nodal values of p, i.e.: if =
t’/2;
it =I$
tj
Jq
and (4) written in terms of i is now ti=a;llTi;-+a:+(i:Ji-/i
Fig. 1. Quadratic elements.
Fig. 2. Coordinates near the tip.
Dynamic analysis of cracks
1053
where a!=i,!;
$=
-i:+4i:-3i)
and
$=2i!-4if.
The singular quarter-point boundary element includes a representation of the traction by means of (6) and a representation of the displacement by means of (2). Thus, both displacements and tractions may be represented, within the vicinity of the crack tip by the first three terms of their series expansion. The first and second mode stress intensity factors can be defined by the following limits (Fig. 2). K, = lim (2n~,)‘/~a 22 x, -0
K,, = limo(27tx,)“2a,2.
(7)
If boundary discretization is carried out in such a way that the first element to cut the domain from the crack tip follows 8 = 0” and this element is a singular quarter-point boundary element, then for this element, P = xi, C,= CJ,2, f2 = oz and the nodal values for the tractions at the tip node K are: if = lim t:(p/Z)“2 = lim o,~(x,//)‘/~ i-0
if: = fi
x, -0
tf(P/l)‘12 = liyo 022(x, /l)‘12.
(8)
Thus, the stress intensity factors coincide with the traction nodal values except for a constant K, = i$(2d)i’2 K,, = i$(2d)i’2
(9)
and may be computed directly with the boundary element code. Since the above procedure makes use of the correct representation of the tractions and only depends on the nodal value at the tip, little dependence upon local boundary discretization may be expected. Blandford et a/.[251 used the singular quarter-point boundary element to compute stress intensity factors in elastostatics by means of a correlation formula of the displacements along the crack surfaces that has been previously used in finite elements[29]. They studied the effects of the size of the elements near the tip and of the use of transition elements. Martinez and Dominguez[26] showed how the use of the traction nodal values of the singular element at the crack tip (9) is substantially less sensitive to the discretization than any of the displacement correlation procedures. When working in steady state elastodynamics, it is even more important to have a procedure with very little dependence on the size of the elements near the tip. This is because it would be difficult to establish rules of discretization that guarantee accurate values of the stress intensity factors if this dependence is important and changes with frequency as can be expected. In the following, the singular quarter-point boundary element in combination with eq. (9) is used for dynamic stress intensity factor computations. DIFFRACTION OF WAVES BY A CRACK IN THE COMPLETE PLANE The problem of harmonic elastic waves that propagate in the complete x,x,-plane and are diffracted at a crack whose length is 2~2,is analysed. The direction of propagation of the waves form an angle y with the crack (Fig. 3). The stress and displacement fields due to P or SV plane waves in a two-dimensional domain can be represented in terms of two scalar potentials: --iF(x,cosy
4 =fjoexp 1 *=o EFM 3+5/6-D
P
+x,siny
+C,t) 1
(10)
FRANCISCO CHIRINO and JOSE DOMINGUEZ
1054
Incident
Waves
Fig. 3. Waves impinging on a crack.
for P-waves, and
-iF(x,cosy +x,siny +C,t) (11) s 1 I for SV-waves; where y is the angle that the wave makes with the x,-axis; C, and C, are the P and S-waves velocity, respectively; o is the frequency and &, Il/,,are amplitudes. Displacements and stresses are easily obtained by differentiation of these potentials. The analysis of wave diffraction by a crack in an infinite medium is carried out by the superposition of two problems. One, the incident field in the domain without crack; and the other, the cracked domain loaded on the crack faces by tractions equal and opposite to those that appear in the untracked plane along the lines where the crack faces are in the cracked plane (scattered field): $ =&exp
uf = uy.+ ui;
al,=a?+a,i
(12)
where ~7 and gi are the displacement and stress components for the incident field, while U, and rrgcorrespond to the scattered field. Only the second problem has to be solved because the incident field in the complete plane is known. Since there are not infinite values of the stress in the untracked plane, the stress intensity factors of the original diffraction problem are the same as in the second problem (scattered field). The Boundary Element Method formulation for elastodynamics is based on the integral representation of the boundary points displacement which, once the boundary conditions are applied, becomes an integral equation that is solved numerically by discretization of the boundary into elements (see, e.g. Dominguez and Alarcon[28]). When the problem at hand satisfies the radiation conditions the displacement representation is written for the external region in terms of integrals that only extend over the internal boundaries (see, e.g. Eringen and Suhubi[30]). This fact is the basis of an important difference between domain methods (Finite Elements, Finite Differences) and the Boundary Element Method where unbounded domains are represented by elements only on the internal boundaries and no problem exists because of spurious reflections on external boundaries. To avoid the problem of having two different values of the displacement for points on the faces of the crack, a boundary that divides the domain into two parts along the x,-axis is introduced. This boundary extends to infinity. However, the discretization is truncated at a distance to the crack tip equal to fifteen times the half-length of the crack. This can be done because elements that are far from the crack have very little effect on the solution near the tip since both the scattered field and the fundamental point load solutions satisfy the regularity conditions. Figure 4 shows the
a
1
15a 1
Fig. 4. Boundary discretization.
Dynamic analysis of cracks
1055
number of elements used in one half of the boundary. The problem is always decomposed into its symmetric and skewsymmetric parts. The two elements to which the node at the crack tip belongs are singular quarter-point elements, whereas all the others are regular quadratic elements. Stress intensity factors for P and SV waves with several angles of incidence have been computed. The values of Ki and K,, depend on the incident angle y and the Poisson’s ratio of the elastic medium. Figure 5 shows the moduli of the stress intensity factors for P-waves normalized with respect to their corresponding static values. Results are plotted vs the dimensionless frequency oa/Cp for a Poisson’s ratio 0.25. Plane strain is assumed. It can be seen in the figure how K, is maximum for y = 90” and 4, for y = 45”. In all cases the values tend to rise first, reach a maximum, and then fall. The dynamic stress intensity factor can be as much as 1.3 times the static one. The values computed using the BEM with singular quarter-point elements are compared with those obtained by Chen and Sih[31] using a method derived by Copson[32] to solve a system of integral equations. As can be seen in the figure, the agreement is very good. The values of $ and K,, for W-waves are shown in Fig. 6. As could be expected the maximum value of K, is not for y = 45” and the maximum of K,, for y = 90”. The agreement with Chen and Sih[31] is again very good.
I I
*
0.1
t
0.2
0.3
0.4
0.5
0.6
Chen
v = 0.25
watcp
0.7
& Sih B.E.M.
.
0.27
+
I
0.1
0.2
0.3
0.4
0.5
e
I
0.6
0.7
u PICp
Fig. 5. Ratio of dynamic and static stress intensity factors for incident P-wave.
FRANCISCO CHIRINO and JOSE DGMINGUEZ
1056
Chen
8, Sih B.E.M.
A
I
0.1
1
a
1
8
1
,
I
0.2
0.3
0.4
0.5
0.6
0.7 oalcp
0.6
07 wa/cp
1. 2
1.0
0.8
Chen
& Sih
-
Y
y=60°
0.6
0.1
0.2
0.3
0.4
0.5
Fig. 6. Ratio of dynamic and static stress intensity factors for incident W-wave.
CRACK
NEAR A FREE SURFACE
In recent years, Achenbach and co-workers[9,33,11,12] have studied the scattering of harmonic waves by edge or sub-surface cracks. Their formulation leads to a system of integral equations which is solved numerically. The fundamental solution of these integral equations corresponds to dislocations in the half-plane and is calculated numerically as a combination of Hankel functions and contour integrals. In the following, stress intensity factors are computed for a horizontal subsurface crack subjected to harmonic excitation. The problem has been taken from Keer et aZ.[l l] for comparison. The regular frequency domain dynamic formulation of the BEM with the complete space fundamental solution is used. Thus, an integral equation has to be solved numerically, but, the kernels of the integral are much simpler and easier to compute than those of the dislocations in the half-plane. On the other hand, the free surface boundary conditions are not satisfied by the fundamental solution and have to be enforced by discretization of the free surface into boundary elements for which zero tractions are prescribed. However, the elements far from the crack have very little effect on the stress field near the crack and accurate representation of the effect of the free surface can be obtained by discretizing a limited zone close to the crack. Due to the nature of the Boundary Element Method, the fact of extending the discretization only to a finite part of the free surface does not introduce undesired artificial reflections. Figure 7 shows the boundary element discretization used for the problem at hand. The free surface discretization and the internal boundary extend to a distance from the crack of 15~ along
1057
Dynamic analysis of cracks
Fig. 7. Boundary element discretization for subsurface crack.
the x, axis. All the elements are quadratic and those at the crack tips are quarter-point singular. The incident field is a uniform tension applied on x2 = 0 which results in uniform tension on the crack faces. The Poisson’s ratio is 0.3. Figure 8 shows the values of K, (magnitude of the ratio of the dynamic and static mode I stress intensity factor) vs the parameters d/a and ad/C,, with C, as the Rayleigh wave velocity. The results are in good agreement with those taken from the figures of Keer et aZ.[ll]. The same comparison is done in Fig. 9 for the mode II stress intensity factor. CENTRAL
CRACK IN A RECTANGULAR
PLATE
The following example has been taken from Chen and Wilkins[34] and shows how the singular quarter-point boundary element can be applied to bounded domain problems. Figure 10 presents the geometry and boundary conditions of the cracked plate. The mechanical properties are: Shear modulus p = 76.923 GPa, Poisson’s ratio v = 0.3; and density p = 5000 kg/m’. In addition to the purely elastic material, internal damping is considered by means of a complex shear modulus c~f= ~(1 + 2pi). Values of j? = 0.01, fi = 0.025 and /I = 0.05 are assumed. The load is a uniform tension applied on two opposite sides. Figure 11 shows one quarter of the plate with the boundary divided into three different zones: zone A (a + L, = 2a), zone B (L,) and zone C (L, + L., + L5). The boundary discretization was made according to the following rules: The length of the elements should be I< J./l0 in zone A, 1 < I.15 in zone B, and I < 113 in zone C, with rZbeing the length of the S-waves in the plate material. The magnitude of the ratio of the dynamic and static mode I stress intensity factor for three different values of the internal damping, is shown vs frequency in Fig. 12. Several peaks of resonance can be seen in the figure. The peaks are shifted and damped with increasing values of /J. The first peak is 20% below the first natural frequency of the untracked plate. Since the only 13 d/a = 04
12 II
B.E.M. 0.E.M. KEER
/VP
-
et al.
KEER
et aI.-----
-----
d/a=2
Fig. 8. Ratio of dynamic and static mode I stress intensity factor for subsurface crack and vertically incident P-wave.
Fig. 9. Ratio of dynamic and static mode II stress intensity factor for subsurface crack and vertically incident W-wave.
1058
FRANCISCO CHIRINO and JOSE DOMINGUEZ
v
=
P
= 76.923
V.JV
GPa
P= 5 Ton.lm3 P=
0.4
l-4
1 1
GPa
Fig. 11. Model for centre cracked plate.
Fig. 10. Centre cracked plate.
results known for the cracked plate problem correspond to transient loads, comparison will be done in the next section cone the above results are transformed into the time domain. TRANSIENT
STRESS
INTENSITY
FACTORS
Going back to the problem of diffraction of waves by a crack in an infinite medium (Fig. 3), a normal incident P-wave with a step function time variation is considered. Results for transient stress intensity factors may be obtained by using the frequency analysis of the problem carried out previously and the FFT algorithm. Thau and Lu[6] computed values of the stress intensity factors for this problem using a generalized Weiner-Hopf technique, but they were only able to study the problem until the time when a diffracted P-wave on reaching the opposite edge is rediffracted and comes back to the original edge. A comparison of the boundary element results with Thau and Lu[6] may be seen in Fig. 13 where Ki is shown vs dimensionless time. The agreement is again good, and K, reaches a value that is 30% higher than the static. Sih et al.[7j studied the same problem
DAMPING 0%
-
ua/cs .b
0.6
0.6
1
1.2
Fig. 12. Ratio of dynamic and static mode I strw intensity factor for centre cracked plate under uniform harmonic traction.
Dynamic analysis of cracks
1059
y” 0.6 .. 0.E.M.
-
THAU- LU -----
Fig. 13. Ratio of transient and static mode I stress intensity factor for crack in the complete plane under a step pressure.
using integral transforms coupled with the technique of Cagniard. They could obtain values of K, for longer periods of time, but because of numerical errors they could not represent the discontinuity in the slope, shown by the BEM and by Thau and Lu, which is a consequence of the arrival of the Rayleigh waves diffracted by the other tip. A comparison of the boundary element results with Sih et a1.[7j is shown in Fig. 14. The last problem to be analysed corresponds to the transient response of the centrally cracked plate that has already been studied in the frequency domain (Fig. 10). A uniform traction P on two opposite sides that varies with time as a Heaviside step function is considered. Again the FFT algorith is used to compute the transient response. The magnitude of the mode I stress intensity factor normalized by P(t)r au is represented for the first 12 ps in Fig. 15. Values for damping factors of l%, 2.5% and 5% are plotted. In this figure, results are compared with those obtained by Chen[l7] using finite differences, and with Baker’s[35] of the semi-infinite crack under pressure in an infinite plate. The latter are valid for times prior to the arrival of the first wave diffracted by the other tip (R,). Both consider purely elastic material. It can be said that the results am in good agreement with Chen’s and Baker’s and that, as could be expected, the existence of internal damping reduces the maximum value reached by the stress intensity factor. This reduction becomes more important as the internal damping increases. Results for the purely elastic material are not shown because they present spurious oscillations that can be explained by numerical errors in the Fourier transformation process. The use of filters could solve the problem, but, since the internal damping makes the spurious oscillations disappear, it is easier to introduce a small amount of damping in the material.
4
Fig. 14. Ratio of transient and static mode I stress intensity factor for crack in the complete plane under a step pressure.
FRANCISCO CHIRINO and JOSE ~MINGUEZ
1060
1.6 Cl.2
0.8 0.1 0
-0.5 2
L
6
a
Fig. 15. Ratio of transient mode I stress intensity factor and P(t)& uniform traction.
10
12 t)ts
for centre cracked plate under step
CONCLUSIONS A procedure for dynamic stress intensity factor computations using the Boundary Element Method in the frequency domain has been proposed. The stress intensity factors are computed from the nodal values of a traction singular quarter-point element located at the crack tip. The length of the singular and adjacent elements should be smaller than L/10, with 3, as the length of the S-waves, and also smaller than 0.6~2,with LZbeing the half-length of the crack. Elements far from the crack should not be longer than L/3. It has been shown that this direct and very simple procedure, previously used by one of the authors in eiastostatics, produces accurate results for cracks in the complete space, in the half-space, or in finite bodies. The method is very versatile, it may be easily included in any quadratic boundary element program, and can be used for many different geometries, including problems with several cracks. The study has been done in the frequency domain and it has been shown how transient stress intensity factors can be computed using the Fourier Transform. A study of the use of the singular quarter-point boundary element with the time domain formulation of the BEM for transient stress intensity factor computations, will be presented in a forthcoming paper. ,4cknowZegemenrs-The authors would like to express their gratitude to the Spanish Comision Interministerial de Ciencia y Tecnologia for supporting this work under the research grant No. PBS6-0139.
REFERENCES [I] De Hoop, Representation theorems for the displacement in an elastic solid and their applications to elastodynamics diffraction theory. Thesis, Technische Hogeschool, Delft (1959). [2] G. C. Sih and J. F. Loeber, Wave propagations in an elastic solids with a line of discontinuity or finite crack. Q. uppl. Murk 27, 193-213 (1969). [3] A. K. Mal, Interaction of elastic waves with a Griffith crack. Znt. J. Engng Sci. 8, 763-776 (1970). [4] A. K. Mal, Diffraction of elastic waves by a penny-shaped crack. Q. appl. Math. 26, 231-238 (1968). [5] A. K. Mal, Interaction of elastic waves with a penny-shaped crack. Znt. J. Z%rgngSci. g, 381-388 (1970). [6] S. A. Thau and T. H. Lu, Transient stress intensity factors for a finite crack in an elastic solid caused by a dilatational wave. Znt. J. Solids Struct. 7, 731-750 (1971). [7] G. C. Sih, G. T. Embley and R. S. Ravera, Impact response of a finite crack in plane extension. Znt. J. Soliak Struct. 8, 977-993
(1972).
[S) L. B. Freund, The stress intensity factor due to normal impact loading of the faces of a crack. ht. J. Engng Sci. 12, 179-189 (1974). 191 J. D. Achenbach, L. M. Keer and D. A. Mendel~hn, Elastodynamic analysis of an edge crack. J. uppl. Mech. 47, 551-556 (1980). (IO] S. F. Stone, M. L. Ghosh and A. K. Mal, Diffraction of antiplane shear waves by an edge crack. J. uppl. Mech. 47, 359-362 (1980). [I l] L, M. Keer, W. Lin and J. D. Achenbach, Resonance effects for a crack near a free surface. J. uppl. Mech. 51,65-70 (1984).
Dynamic analysis of cracks
1061
[12] W. Lin, L. M. Keer and J. D. Achenbach, Dynamic stress intensity factors for an inclined subsurface crack. J. uppl. Mech., 51, 773-779 (1984). [13] A. H. Shah, K. C. Wong and S. K. Datta, Surface displacements due to elastic wave scattering by buried planar and non-planar cracks. Waue Motion 7, 319-333 (1980). [14] B. R. Tittmann, L. A. Ahlberg and A. K. Mal, Rayleigh wave diffraction from surface-breaking discontinuities. Appl. Phys. Lett. 20, 1333-1335 (1986). [15] Y. C. Angel and J. D. Achenbach, Stress intensity factor for 3-D dynamic loading of a cracked halfspace. J. Efasriciry 15, 89-102 (1985). [16] J. H. M. T. Van Der Hijden and F. L. Neerhoff, Scattering of elastic waves by a plane crack of finite width. J. appf. Mech. 51, 646-651 (1984). [17] Y. M. Chen, Numerical computation of dynamic stress intensity factor by Lagrangian finite-difference method. Engng Fracture Me&. 7, 653-660 (1975). (181 T. Nishioka and S. N. Atluri, A numerical study of the use of path independent integrals in elastodynamic crack propagation. Engng Fracture Me&. 18, 23-33 (1983). [19] S. N. Atluri and T. Nishioka, Numerical studies in dynamic fracture mechanics. Znr. J. Fracture 27, 245-261 (1985). [20] S. Aoki, K. Kishimoto, H. Kondo and M. Sakata, Elastodynamic analysis of crack by finite element method using singular element. Inr. J. Fracrure 14, 59-68 (1978). [21] V. Murti and S. Valliappan, The use of quarter point element in dynamic crack analysis. Engng Fracture Mech. 23, 585-614 (1986). [22] M. D. Snyder and T. A. Cruse, Boundary integral equation analysis of anisotropic cracked plates. Inr. J. Frucrure 11,
315-328 (1975). [23] M. Stem, E. B. Becker and R. S. Dunham, A contour integral computation of mixed-mode stress intensity factors. Inr. J. Fracture 12, 359-368 (1976). [24] T. A. Cruse, Two-dimensional BIE fracture mechanics analysis. Appl. Math. Modelring 2, 287-293 (1978). [25] G. E. Blandford, A. R. Ingraffea and J. A. Liggett, Two-dimensional stress intensity factor computations using the boundary element method. Inr. J. Numer. Merh. Engng 17, 387-404 (1981). [26] J. Martinez and J. Dominguez, On the use of quarter-point boundary elements for stress intensity factor computations. Inr. J. Numer. Merh. Engng 20, 1941-1950 (1984). [271 C. A. Brebbia and J. Dominguez, Boundary Elements. An Inrroductory Course. CMP Publications McGraw-Hill, New
York (1989). [28] J. Dominguez and E. Alarcon, Elastodynamics, in Progress in Boundary Efernenr Merhocis (Edited by C. A. Brebbia). pp. 213-256. Pentech Press, Plymouth (1981). [29] C. F. Shih, H. G. de Lorenzi and M. D. German, Crack extension modeling with singular quadratic isoparametric elements. Znr. J. Fracture 12, 647651 (1976). [30] A. C. Eringen and E. S. Suhubi, Eiastodynamics. Academic Press, New York (1975). [31] E. P. Chen and G. C. Sih, Scattering waves about stationary and moving cracks, in Mechanics of Fracture: Efustodynumic Crack Problems (Edited by G. C. Sih), pp. 119-212. Noordhoff, Leyden (1977). [32] E. T. Copson, On certain dual integral equations. Proc. Glasgow Math. Assoc., Vol. 5, pp. 19-24 (1961). [33] J. D. Achenbach and R. J. Brind, Scattering of surface waves by a sub-surface crack. J. Sound Vibration 76, 43-56 (1981). [34] Y. M. Chen and M. L. Wilkins, Numerical analysis of dynamic crack problems, in Mechanics Elasrodynamic Crack Problems (Edited by G. C. Sih), pp. 295-345. Noordhoff, Leyden (1977). [35] B. R. Baker, Dynamic stresses created by a moving crack. J. appl. Mech. 29, 449-545 (1962). (Received 7 February 1989)
of Fracture: