Martin G Koller

Marc Bonnet

Raul Madariaga

Modelling of dynamical crack propagation using time-domain boundary integral equations

We study dynamic antiplane cracks in the time domain by the boundary integral equation method (BIEM) based on the integral equation for di� cement discontinuity (or crack opening displacement, COD) as a function of stress on the crack. This displacement discontinuity formulation present; the advantage, with respect to methods developed by Das and others in seismology, that it has to be solved only inside the crack. This BIEM is, however, difficult to implement numerically because of the hypersingularity of the kernel of the integral equation. Hence it is rewritten into a weakly singular form using a regularization technique proposed by Bonnet. The first step, following a method due to Sladek and Sladek, consists in converting the hypersingular integral equation for the displ&cement discontinuity into an integral equation for the displacement discontinuity and its tangential derivatives (dislocation density distribution); the latter involves a C9uchy type singular kernel.

The second step is based on the observation that the hypersingularity is related to the static component of the kernel; the static singularity is then isolated and can be expressed in terms of weakly singular integrals using a result due to Bonnet.

Alth<J1ugh numerical applications diso.::ussed !n this paper are all for the antiplane problem, the technique can be applied as well to in-plane crack dynamics.

The BIEM is implemented numeri<.-:ally using continuous linear space-time base functior.s to model the COD on the crack. in the present scheme the COD gradieu� ii"tt�&.,ulation is discontinuous at the elem� .. 1t nodes while the integral equations are collocated at the element midpoints. This leads to an overdetennined discrete problem which is solved by standard least squares methods. We use the dynamic DIEM to study a set of problems that appear in earthquake source dynamics, including the spontaneous dynamic crack propagation for a very simple rupture criterion. The numerical results compare favorably with the few exact solutions that are available. Then we demonstrate that difficulties experienced with finite difference s!mulations of spontaneous crack dynamics can be removed with the use of BIEM. The results are improved by the use of singular crack tip elements.

Introduction

Many numerical studies of shear crack pr"pagation viewed as a possible model of shallow earthquakes have been published in the past twenty years. In this article, a new numerical technique is presented which is more generally applicable and of higher precision than those used up to now. The present numerical approach is based on a regularized time-domain traction boundary integral equation (RIE). Our objective is to demonstrate the potential of this tnethod by treating the antiplane case which is simplest, both fr om a geometrical and a fracture mechanical point of view. We consider an initially straight two-dimensional antiplane (mode Ill) shear crack, embedded in an infinite linearly elastic homogeneous medium, its spon taneous propagation being governed by a simple critical stress intensity criterion. The present numerical model allows out-of-plane as well as in-plane propagation. Here the onset of a potential kinking of the crack path is studied. A focus of growing interest in seismology is the dynamic interaction between faults, both on a regional (earthquakes ��jumping" from one fault segment to another, e.g. [START_REF] Bouchon | Predictability of ground displacement a.�d velocity near an earthquake fa ult : an example: the Parkfield earthquake of 1966[END_REF]) as well as on a micro scale ("en echelon" fractures). Several works have been published on static fault interaction, see e.g. [37, 29], but to our knowledge the dynamical interaction problem has been considered only in paper [START_REF] Harris | Fault steps and the .r!ynamic rupture process : 2D numerical simulations of a spontaneously propagating shear fracture[END_REF]. Further problems of interest that may be studied with our numerical scheme are the modelling of rupture along an existing kinked fault (e.g. El Asnam [START_REF] Yielding | Active tectonics in the;: Algerian Atlas Mountains -evidence from aftershocks of the 1980 El Asnam earthquake[END_REF]) or the determination of the terminal velocity of rupture propagation on a new or considerably recrystallized fault by allowing the crack to zigzag around its overall path in analogy to what is observed in experiments with tensile cracks. Two different numerical approaches have bee n used so far in the quantitative study of dynamical earthquake source mechanics: a displacement BIE fo rmulation and several finite difference (FD) techniques. Early works were those of Burridge [START_REF] Burridge | The numerical solution of certain integrai equations with non-integrable kernels arising in the theory of crack propagation and elastic wave diffraction[END_REF], Das and Aki [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF], Andrews [START_REF] Andrews | Rupture propagation with finite stress in antiplane strain[END_REF] on BIE, Andrews [START_REF] Andrews | Rupture velocity of plane strain shear cracks[END_REF], Madariaga [START_REF] Madariaga | Dynamics of an expanding circular fa ult[END_REF], Virieux and Madariaga [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] on FD. Although both approaches were eventually applied to rather sophisticated problems, including inhomogene ous material resistance or stress fields and three-dimensional elasticity problems associated with bounded plane cracks (see [START_REF] Das | Earthquake Source Mechanics[END_REF] for further references), all studies have dealt only with single plane cracks. It is indeed virtually impossible to treat geometrically non-planar crack problems with the displacement DIE. FD techniques, on the other hand, are not limited in this respect, although the possibilities of dealing with kinked cracks are restricted by the actual configuration of the computational grid. In fact, as pointed out in [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF], the accuracy of second-order FD techniques is unsatisfactory for the study of spontaneous rupture. The main numerical problems arising concern the stress resolution near the crack tip and the numerical dispersion by the computational grid. Even the method proposed by Trifu and Radulian [START_REF] Trifu | Predicted near-fi . eld ground motion for dynamic stress-drop models[END_REF], which considerably improves the stress resolution, does not seem to er .. hance significantly the efficiency of FD techniques. In contrast, our regularized time-domain traction BIE, or related approaches like that of Hirose and Achenbach [START_REF] Hi Rose | Time-domain boundary element analysis of elastic wave interaction with a crack[END_REF], are expected to overcome these shortcomings.

Time-domain boundary integral equations for crack analysis

Throughout the present paper, the discussion is restricted to 20 antiplane time-domain elasto-dynamics (shear modulus p., specific mass p). Zero body forces and homogeneous initial (t=O) conditions are assumed. Accordingly, a domain !2 c IR 2 of boundary o.Q is considered. The coordinate directions ( e1 , e 2 , e 3) are chosen so that the only nonzero displacement and stress components are u 3 (y, r), u37(y, r) (r= I, 2) respectively (in the sequel, Greek indices always take the values 1, 2). The unit normal vector n = nrer is directed towards the exterior of !2. Hooke's law reduces to OUJ O'jy = (j y 3 = J.l -(y, t), oyr G�;=O otherwise.

(1)

2 . .f. BIE for antiplane elastodynamic crack probletns Under these assumptions, the displacement u3(z:, t) at an interior point z = z r e r eD can be expr�ssed in terms of boundary displacement u3 and traction T3 = a3 r n r [22] using the fo llowing weB-known integral repre-�entation formula: u 3 (z , t) = -f f' { U3 (y, r)n,(y) IJ r;J(Z, y, t-r)-T3( y, r) U3;3(Z, y, t-r)} d r ds,. , z�(JQ [START_REF] Achenbach | Bifurcation of a running crack in antipiane strain[END_REF] an o where U3;3(z,y, tr) denotes the only nonzero component of the impulsive antiplane Green's tensor for the infinite space [22], i.e. the antiplane displacement at yeiR2 and time t created by an impulsive antiplane point load applied at ze R 2 and time r�t, I3r;J(z,y, t-r) being the (3, r)-components of the elastic stress tensor associated with U3;3(z,y, t-r) through Hooke's ia.�• !lt The analytical expressions of antiplane Green's tensors are given by eqs. [START_REF] Yielding | Active tectonics in the;: Algerian Atlas Mountains -evidence from aftershocks of the 1980 El Asnam earthquake[END_REF] and [START_REF] Ch | A novel derivation of non-hypersingular time-domain BIEs for transient elastodynamic crack analysis[END_REF] of Appendix A. In the sequel, the field components u3, T3, U 3 ;3, E3r; 3 will be simply denoted by u, T, U, E,, the unnecessary '3' indices being dropped for convenience.

The representation formula [START_REF] Achenbach | Bifurcation of a running crack in antipiane strain[END_REF] is the basis for boundary integral equation (BIE) methods for 20 antiplane problems. Among them, the well-known and widely used displacement BIE [START_REF]Boundary Element Methods in Mechanics[END_REF][START_REF] Brebbia | Boundary Element Te chniques[END_REF] is obtained through a limiting process z-+ xeo.O, where the internal source point z is moved towards a fixed point x located on the boundary in [START_REF] Achenbach | Bifurcation of a running crack in antipiane strain[END_REF] (see e.g. [7, 12� 23]) (the distinct notations x, z will be used as a reminde1 of this limit process). The integrations involving the kernels E,(z,y, t-r) are then to be understood as Cauchy principal values (CPV). The actual numerical evaluation of CPV integrals has been a source of numerical difficulties, which were overcome recently using either regularization approaches [START_REF] Bonnet | Methode des equations integrates regularisees en elastodynamique tridimensionelle[END_REF][START_REF] Kuhn | A boundary integral equation method for radiation and scattering of elastic waveo; in three-dimensions[END_REF] or a direct numerical treatment of general CPV integrals [START_REF] Guiggiani | A general algorithm for multidimensional Cauchy principal value integrals in the boundary[END_REF].

However the BIE modelling of crack problems, in the framework of elastic fracture theory, involves further specific difficulties, which are not present when considering non-cracked domains. Let from now on n be an unbounded 2D elastic domain containing a crack (upper and lower faces C + and C-respec tively): aa=c+ u c-, c+ and c-being in fact the S8.nie geometrical curve C=(AB) in the XJXrplane with opposite unit normals n .,.. , n -(Fig. 1 ). The unit normal v�tor n = n -= (n1 , n 2 ) will be used, together with the unit tangent vector T= (rh r 2 ) defined by Let us introduce the notation q,(y, r)=u(y+, r)-u(y-. r) for the crack opening displacement (COD); in the subsequent discussions, q, will be also sometimes termed as "'slip", as is usually done in geophysics. Taking into account a.a = c+ u c-, n+ = _ ,-= -n, T+ = -T-, the representation formula (2) leads to the fo llowing integral representations of the displacement and (upon application of Hooke's Law (1)) stresses for the anti plane crack problem: u(z, t)= tf 1/>(y, r)nr(y)Er(z,y, t-r) d r dsy, z�C [START_REF] Andrews | Rupture velocity of plane strain shear cracks[END_REF] 0"3p( z, t) = -p f J' q,(y , r)n7(y) �E7(z, y, t -r) dr ds y , z�C.

In view of the well-known degeneracy of the displacement BIE, which is due to the fact that the boundary data T+, Tdo not appear explicitly in (4), an alternative possibility is the statement of a traction BIE, in order to relak explicitly the traction data to the COD [8, 9, 13, 27 , 35, 38, 46, 48, 49]. The latter is obtained by a limiting process (z�C)-+ (xeC) in [START_REF] Andrews | Dynamic plane-strain shear rupture with a slip-weakening friction law calculated by a boundary integral method[END_REF]. Due to the hypersingularity of the kernel components (cjcyp).E 7 (z, y , t•-i") for:= rand z= y , the resulting hypersingular integral over Cis understood as a finite part (FP).

Because of the hypersingularity problem, which is a considerable obstacle to the nunterical imp�ementa tion of the traction BIE, researchers formerly developed and applied alternative approaches to dynamical crack analysis, leading io BIE involving integrable singularities, but having less generality regarding the geometry of the problems under consideration. Many of these works were motivated by the study of seismological problems. For example, when considering planar cracks in infinite elastic domains, the prob lem can be stated, using symmetry with respect to the crack plane P, on the half-space n+ bounded by P (this is a particular case of "substructuring"). One can then apply the displacen1ent BIE to the domain n+, using the half-space Green's tensor satisfying a free-surface (zero traction) condition on .l". This leads to a weakly singular displacement BIE, in which the unknown is the traction Ton P-C. However, tili� method allows only the study of straight cracks, fo rbidding the consideration of out-of-plane propagatio1� of kinking, originally non-straight cracks, interaction between non-collinear cracks, etc. Using this approach, Burridge [START_REF] Burridge | The numerical solution of certain integrai equations with non-integrable kernels arising in the theory of crack propagation and elastic wave diffraction[END_REF] studied dynamically loaded fixed cracks in anti plane and plane strain, while Das and Aki [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF], Andrews [START_REF] Andrews | Dynamic plane-strain shear rupture with a slip-weakening friction law calculated by a boundary integral method[END_REF] considered propagating cracks, either spontaneously or at fixed velocity, the latter introducing a slip-dependent friction law on the crack plane. The Green's tensor associated to a given domain Q is generaily now known, the half-space case being an exception, thus the kind of approach reported here is hardly generalizable. Moreover any displacement RIE using a Green's tensor associated to a domain D' greater than Q (generally D' == IR 3 ) contains CPV integrals.

Regularization of traction BIE for 2o antiplane crack problems

The application of collocation traction BIE to the numerical study of crack problems requires careful handling of the hypersingular integrals. One approach consists in a direct use of the FP traction BIE obtained :rom (5) together with a suitable treatment of the eletnentary FP integrals [START_REF] �l Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF] in which the 1!1�pping between physical element and reference element carefully preserves the sym:�netrical vanishing exclusio�1 ne:ghborhood used to define the value of FP integrals. Another approach consists in rewriting the hypersingl.dar BIE, using a regularization approach, into an equivalent form which involves only weakly singular integrals. The regularization of general three-dimensional elastodynamjc traction BIE is treated in [START_REF] Bonnet | Regularized boundary integral equations for three-dimensional bounded or unbounded elastic bodies containing curved cracks of arbitrary shape under dynamic loading[END_REF][START_REF] Bonnet | Regularization of the displacement and traction BIE for 3D elastodynamics using indirect methods[END_REF], see also [32]. Let us finaHy mention fo t completeness that regularization can also be achieved by deriving a Galerkin-type variational BIE, see e.g. [START_REF] Becache | Resolution par une methode d'equations interales d'un probleme de diffraction d'ondes elastiques transitories par une fissure[END_REF] for time-domain crack problems.

The present work is an application of the regularization of collocation traction BIE to the numerical simulation of 2D propagating cracks under antiplane shear. The regularization process [START_REF] Bonnet | Regularized boundary integral equations for three-dimensional bounded or unbounded elastic bodies containing curved cracks of arbitrary shape under dynamic loading[END_REF] leading to a weakly singular traction BIE is outlined as fo llows.

In eq. ( 6) and the sequel, the asterisk (•) denotes a CPV integral. The CPV traction BIE (6) generalizes, for the two-dimensional case, the results of Bui [START_REF] Bui | An integra! equation method for solving the problem of a plane crack of arbitrary shape[END_REF], Weaver [START_REF] Weaver | Three-dimensional crack analysis[END_REF] in elastostatics. Z!lang and Achenbach [START_REF] Ch | A new boundary integral equation fo rmulation for elastodynamic and elastostatic crack analysis[END_REF], Zhang [48] also obtained the BIE (6) for two-or three-dimensional problems using the concept of ]k conservation integral. • Then the CPV traction BIE (6) is regularized as fo llows. First eq. ( 6) is rewritten into the fo llowing equivalent equation [START_REF] Bonnet | Methode des equations integrates regularisees en elastodynamique tridimensionelle[END_REF][START_REF] Bonnet | Regularized boundary integral equations for three-dimensional bounded or unbounded elastic bodies containing curved cracks of arbitrary shape under dynamic loading[END_REF] : + I [e,p � 1/J(y, t) -e ,11 � 1/J(x , t)}�(x, y) dsy c os os _J +p I f' np(y)U(x,y, t-r)(j}(y, r) d r ds y+ e,p ! l/J(x, t) r• .!'�(x,y) d sy} (7) c 0 os 'c

± _ ± {I (J'
i!1 which .!'�(x, y) denotes the static counterpart of .!',(x, y, tr) (given by eq. ( 49) in Appendix A), i.e.

the sttess kernel associated to an antiplane static unit point fc :-ce. All integrals but the last in (7) are regular or weakly singular; this relies crucially upon the equality of the singular parts, for y = x, of the dynamic and static kernels �(x,y, t-r) and E>(x ,y) , a.nd also on the assumption that (o/os)t/J(y, t) is C 0 •a, a> 0, i.e. Holder-continuous with positive exponent, on C [32, 10]. The residual CPV integral in (7) can be rewritten [START_REF] Bonnet | Methode des equations integrates regularisees en elastodynamique tridimensionelle[END_REF][START_REF] Bonnet | Regularized boundary integral equations for three-dimensional bounded or unbounded elastic bodies containing curved cracks of arbitrary shape under dynamic loading[END_REF], by means of a.n analytic transformation using integration by parts, into a weakly singular fo rm; the result is where r = lly-xll = ( ya-xa><Ya-xa) and R(y) denotes the radius of curvature at y. The integrand of (8), including the term ( l /r)(or/on), is integrable, assuming C to be twice continuously differentiable (this assumption is necessary since (8) involves the curvature R). The line integral over C in (8) vanishes if C is a straight element. Equations ( 7) and (8) define the regularized traction BIE for two-dimensional antiplane crack problems. It is equivalent to the hypersingular traction BIE obtained from (5) but contains only integrable (in the ordinary sense) singularities, thus being suited fo r applications. This result has fu ll generality and applies to cracks of arbitrary (sufficiently regular) shape subjected to arbitrary loadings. Indeed it can be shown [START_REF] Bonnet | Regularized boundary integral equations for three-dimensional bounded or unbounded elastic bodies containing curved cracks of arbitrary shape under dynamic loading[END_REF] that equations [START_REF]Boundary Element Methods in Mechanics[END_REF] and (8) are valid, without modification, for the more general situation of piecewise regular (i.e. kinked� branched or crossed) cracks; in this case, each regular component of C has to be twice continuously differentiable. In practice, it is convenient to divide the curve C into two complementary parts � C(x) and C-L1 C(x), where LiC(x) is any subset of C containing x as an interior point, in order to perform the singularity cancellations and evaluate the residual integral (8) on a curve smaller than C (e.g. the single boundary element containing the singular point x). This leads to the foUowing developed expression of the regularized traction qT� .. obtained through insertion of ( 49), (51), (50), ( 3) and ( 8) in (7):

I * 1 { I [ 1 or 1 J } E�(x ,
r ± , x t) = 1!._ r ± (x) ' ' 21t fJ { f YP -x13 [oc/J aq, ] f YP-xp [ocf> ocf> ] x 2
-(y , t-r/c)--(y,

+� l/J(x, t){[ rp Log r]� + f [! o r _ _!_ Log r}p(y) dsy } & c r� R + f
Yf3�XfJ ol/J (y, t-rj' c) <55y

For expository convenience� the present section describes the discretization scheme for a steady crack. The numerical management of a propagating crack will be described in the next section.

1. Boundary element dis c retization scheme: description and notations

The crack curve C is approximated by J straight elements£,= [YoYd, ... , E1= [YJ-I)'J] of equal length Llx. and a constant time step Lit is used Cv0,y1 are the endpoints A, B). We denote by s the arc•length along C, so= 0, .\•1 = Lix, ... , s.� = JLix denoting the curvilinear abscissae of the subdivision pointsy0, ••• , y 1. The collo c ation points xi= 1 /2( )'1-n + y;) (i.e. the midpoints of the elements E1-... , .C'J) are introduced; their curvilinear abscissae are s;-112=s;-1 +L1x/2. The regularized traction BIE (9) is collocated at the points X; and the time steps t = t 1 = .:1 t, ... , t = t N = N Lit, L1 C( x) being chosen as the boundary element E; containing the current collocation point X;, while the COD 4J is discretized in the usual manner:

J-I ,\' (jJ(y, t) = I L v1(s)8m(t}tf>�1• j� I m� I (10)
The functions v1 and Om are chosen piecewise linear and satisfy so that the v1 and (}m can be written using a single basis fu nction S defined on ( -1, 1] :

with . (s-s• ) ( s ) v ' (s)=S Ax :J =S Ax -j, S(u)= {1 +u, -1 �u�O 1-u, O�u� 1. (t-tm) ( t ) O m(t)=S Tt =S At -n1
(ll)

(12) (13)
In order to incorporate the square-root spatial variation of the COD near the end-points of the crack, however, the functions v1(s) and v1 -•(s) are modified as fo llows:

v J-t (s)= { S(s /t1 x-(J-l)) , (J-2)Lix�s�(J-l)Ax [J -s/ Lix]11 2 , (J -1)Ax�s�JL1x. (14) (15)
The present discretization scheme yields J equations fo r J-1 scalar unknowns. This is due to the fa ct that collocation at element end points is impossible using the present choice ( 12) of shape functions v1 because.

as has been mentioned above, the CPV integral eq. ( 6) as well as its regularized counterpart (9) are valid only fo r collocation points x at which (ojos)c/J (y, t) is C 0 •a. This restriction extends to C 0 polynomial interpolation of any degree, the number of degrees of freedom and of equations matching only for piecewise constant spatial interpolation.

The ratio cLit/Ax has been chosen equal to 1/2. This value is smaller than l/�2 and therefore lies within the stability range of corresponding two-dimensional FD methods [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF]. Furthermore, this choice allows one to perform analytically many of the elementary integrations, especially fo r the singular terms or fo r straight parts of the crack. In another BElvi investigation of transient crack problem• [START_REF] Nishirnura | Boundary integral equation method for clastodynamic crack problems[END_REF], the authors chose c At/ Lix= I and reported good results.

For a given length L of the crack, the insertion of (1 0) into ( 9), with t = t n =-= nLi t and x =X;, leads to the fo llowing system of equations:

or 2 I -l [ n -1 J � T�' + . ' H '• '• •" .J.. '! + ' H �'. • m r�.�n = 0 I LJ I} f.J'J LJ I} 'f'J ' JJ j=l m =l n=l, ... ,N 21t n-1 -T" + H"•"tjl' + L H"• m ; m =0, n = 1, ... ' N J.1
m=l in matrix form. The coefficients H '// ' of ( 16) are given by

H '• ' • •m = HV'?. •m + H G'?.-m + HA �� .. m I} lj 1) I} (16) (17)
HGR'/; m = r ,s(x) f Yp � X p v' 1(s)S(n-m-rjcAt) ds,.
c-£, r HGS'/;m=O, m<n-1

HGS'l/ -1 = -r11(x) f YP � X p v'1 (s) [ S( -rjcL1t) -1] ds_.. E, r 1 r . rJ( S(n-m-r) J HA'!;m= rp(x) 2 r p (y)v1(s) 2 _ 2 2 /2 d r dsy . c , c .._ r.•c ( r r I c ) 3.2.
Comments about the computation of the elementary integrals ( 18) to [START_REF] �l Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF] Equations ( 18) to ( 25) deserve some comments and a closer examination.

I. Equations ( 23)-( 24) take into account the specific value 1/2 of the ratio cAt/ L1x used here.

(

) ( 19
) 20
(2 1)

(22) (23) (24)
(25)

2. The substitution [START_REF] Brebbia | Boundary Element Te chniques[END_REF] has been made for the time basis functions Bm in [START_REF] Dcnpsey | Dynamic crack branching in brittle solids[END_REF], (21), ( 22). ( 23), ( 24), [START_REF] �l Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF], in order to make apparent that the elementary integrals H'/;m possess the translation property

H n_ .m = nn. -m . O I) IJ (26)
which is the discrete equivalent of the convolution structure of the time-domain BIE. H'!; m will be accordingly denoted by Hijm.

3.

The process of removing the spatial singularity of the Stokes tensors appears in the elementary integrals HGS'// and HGS'!/" -1 ) eqs. (23 )--( 24): the 1 jr singularity is cancelled out by vj(s)-vj(s; -1 / 2 ) or S( -r/ ell t)i. Tnese integrais can be evaiuated anaiyticaiiy for the present discretlzation method.

First, the elements are straight segments and x is located at the center of an element� as a consequence, the bracketed tenn (i.e. the residual CPV integral (8)) vanishes. Then, our spatial interpolation choice fo r q, implies that

[£;</>(y, r)-:s </>(x,l) ]=o on£,, 2�i�J-2. (27)
Equation ( 27) does not hold for x = x. and x = XJ i.e. fo r the two end elements, because of the use of the crack-tip basis functions (14) which are not linear. Finally, the remaining integrals in (23)-( 24) are evaluated analytically in Appendix B. The analytic expressions (70) are the materialization, for the present fo rmulation, of our singularity extraction approach which led to BIEs ( 7) and ( 9); all other elementary integrals are regular. The only nonzero terms in (70) are associated to the crack-tip basis functions; this is specific of the present interpolation scheme [START_REF] Brebbia | Boundary Element Te chniques[END_REF] and does not extend to higher-degree basis functions or collocation at non-centered points.

4. The inner integrals with respect to time in eqs. ( 19) and ( 25) can be straightforwardly expressed in terms of elementary functions, irrespective of the actual shape of the crack; this results in eqs. ( 52) and (53) of Appendix A.

5. The outer spatial integrations over the elements E; in (19), (2 1) and ( 25) are to be performed. At this stage a technical difficulty arises, because the integrands, considered as fu nctions of s, have discontinuous derivatives at points such that r I c = kt1 t (k integer). Therefore the standard numerical integration methods, e.g. Gaussian, are likely to give poor r?-su\ts if applied over the whole element E;. Instead one has to partition each element E; into subeleme&lts E�: E;= U £� with E�= {yeE;, (k-l)L1t�rj c �ki1t}.

(28)

Each collocation point x is associated to a different partition [START_REF] Ida | Cohesive fo rce across the tip of a longitudinal shear crack and Griffith's specific surface energy[END_REF]. Giv�n x and E;, [START_REF] Ida | Cohesive fo rce across the tip of a longitudinal shear crack and Griffith's specific surface energy[END_REF] involves obviously a finite number of subelements E� . Then an elementary integral over E; is the sum of "subinte grals" over the £7, each "subintegral" being evaluated using standard numerical methods, or analytical methods in special cases. This technical difficulty is common to all BEM development using retarded potentials and time marching schemes. The partition (28) is clearly a purely geometrical problem (find the intersection points of a given curve C with circles centered at a given point x and �f equally spaced radii). It is reasonably tractable fo r two-dimensional problems. On the contrary, its three-dimensional equivalent (find the intersection curves of a given surface S with spheres centered at a given point x and of equally spaced radii) is extremely complicated and, to our best knowledge, a numerically efficient treatment of this particular point is still an open problem (see however [START_REF] Karabaiis | A stmplified 3-D time-domain BEM fo r dynam ic soil-structure interaction problems[END_REF], where the spatial integrations are performed by partitioning the elements into very small rectangular subelements). 6. The spatial elementary integrations in [START_REF] Dcnpsey | Dynamic crack branching in brittle solids[END_REF] to [START_REF] �l Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF], on each subelement, are done using a numerical method. Integrals containing the linear basis functions vj(s) (2 �j�J-2) are evaluated using an ordi nary Gaussian method (we took 10 points per subsegment), while a Gauss-Chebyshev method has been used fo r the crack-tip basis functions in order to integrate accurately the square-root singularity of the COD gradients.

In the particular case of straight cracks, the elementary integrals have for any value of ci1t/ Ax the following spatial translation prope,.ty: H�•-m =H7�;"= -H'j;i�i [START_REF] Kato | Mechanical interaction between neighboring active fa ults an application to the Atera fa ult. central Japan[END_REF] which stems fr om the uniformity of the space grid. Thus the values of the elementary integrals have to be computed for only one collocation point. The translation property [START_REF] Kato | Mechanical interaction between neighboring active fa ults an application to the Atera fa ult. central Japan[END_REF] does not apply to integrals involving the crack-tip basis fu nctions (i.e. fo r j= 1 or j=J -I), which must therefore be evaluated for each collocation point. All the elementary integrals for a straight crack have been evaluated analytically, these calculations being less complicated when a simple value fo r the ratio ci1 t I L1 x (like the present choice 1 /2) is used.

Time-marching scheme

The system of equations ( 17) is rewritten as

H o tjl = -2 7t rn-n I I Hn-m q,m' n = 1 ' .. . ' N. IJ m= 1
The matrix relation [START_REF] Karabaiis | A stmplified 3-D time-domain BEM fo r dynam ic soil-structure interaction problems[END_REF] exhibits the features classically encountered in time-marching schemes:

1. At each time step, only one new matrix has to be computed, namely H n-•. 2. The whole process needs the factorization of H0 oniy. Due to the discretization scheme introduced in Section 3;the system (30) has J equations for J-1 unknowns and is hence overdetermined. As a consequence, ( 30) is solved for q,n in the I� east-squares sense, as

(31)
The matrix .H 0 is fa ctorized only once, by means of the Householder algorithm [START_REF] Stoer | troduction to Numerical Analysis[END_REF] and using the subroutine SQRDC of the software library LINPACK [21]. This method uses orthogonal transformations. Hence the condition number of H0 is preserved, while it would have been squared, and hence deteriorated, using the classical normal equations of linear least-squares.

Then ( 31) is solved at each time step, using the previously fa ctorized matrix H0 and the LINPACK subroutine SQRSL [21].

Numerical simulation of the crack growth

We consider here the unilatr:ral propagation of an initially straight crack, where tip A is kept fixed. Tip B propagates either with a prescribed rupture velocity or spontaneously (i.e. according to a dynamical rupture criterion, similarly to [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF] or [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF]). The propagation is simulated by adding new elements to the crack ahead of the current location of the tip, which thus advances by discrete jumps of length Lix.

Numerical simulation of the crack propagation

Crack tip propagation is managed as follows. The crack and the COD are initially discretized as described in Section 3, using Jo elements of length Ax. The propagation is simulated by adding new elements to the crack ahead of the current location of the tip, which thus advances by discrete jumps of length Ax.

Let us consider the situation where k-1 (le� 1) new eiements have been added to the initial crack since propagation started: the current crack length is thus 1k-1Ax=(J0+k-l)Ax, and the current crack configuraiton is laibelled ( ck-1 ) . The time-marching process is continued on configuration ( ck-I) until the tnne step for which the rupture criterion (to be described in Sections 4.2 and 4.3 below) is fulfilled. This time step is labelled nb and a new straight element (the kth) is immediately added ahead of the current location of !he propagating tip. As a consequence of the linear time interpolation (13), a new square root singularity buiids up at the new crack tip while the singularity a1t the previous tip vanishes linearly in time, within the next timt: interval [n�c.1t, (n�c+ l)L1t]. Then this new configur�tion (Cd is maintained until the next time step n�c + 1 at which the rupture criterion is verified again.

In-plane spontaneous propagation

For the sake of simplidty and comparison, Irwin's fr acture criterion, generalized to dynamical problems [START_REF] Freund | Dynamic. Fracture �Uechanics[END_REF], is used. It states that the crack tip propagates in such a manner that its stress intensity factor (SIF) K(t) equals at any timet a critical value� (called fr acture toughness), which is a material constant usually assumed to be independent of rupture velocity. Close enough to the crack tip, at a distance x, in front of it, the stress 0"23 is related to K(t), under very general assumptions [START_REF] Freund | Dynamic. Fracture �Uechanics[END_REF], by 0"23(x,) == K(2nx,)-1 1 2 + 0( 1).

(32)

Accordingly, �can be approximately translated, fo r a fixed small x,, into a critical stress level a23(x,) == O"c(x,) which has been used fo r the present numerical computations.

In practice, however, the stress is usually not calculated close enough to the tip (a quantitative study of the spatial extent of the region in which the stress field is actually dominated by its singular crack-tip part has been made in [START_REF] Kobayashi | Initiation, propagation and kin king of an antiplane crack[END_REF]). Thus eq. (32) associates to the physically correct value O"c(x,) an estimated value K, of the fracture toughness which differs fr om its physical value Kc by an adjustment fa ctor o: K, = (I+ «5 )Kc( [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF], eq. ( 21)). The adjustment term o is necessary fo r comparing numerical results (obtained for the critical stress O"c, related by (32) to the estimated K,) to analytical solutions (expressed in terms of Kc). Indeed this correction was used fo r Figs. 4 and5 of [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] and Figs. 5 and6 of the present paper. This adjustment was unavoidable in [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF] or [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] since only a r .. 1ean value of the stress over the closest grid element in fr ont of the crack tip was numerically available. Using the present method, it would have been possible to compute the stresses at points arbitrarily close to the crack tip usjng the integral representation fo rmula (5) in our numerical scheme ; however, they were actually evaluated, fo r practical reasons, at x, = L1x/2, mid-points of (future) boundary elements.

Two ways of discretizing irwin's criterion have been tested. The first or.�, :•eferred to as Hsimple", strictiy corresponds to what was done in [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF] and [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] ; according to it, the next jump of the discretized crack occurs at the time nk+ ,At, nk+, being the first time step such that [START_REF] Kobayashi | Initiation, propagation and kin king of an antiplane crack[END_REF] Although reasonable at first sight-and seemingly successful in [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF] and [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] -this procedure does not give a strict discretization of Irwin's criterion as fo rmulated above, and leads to unsatisfactory results in the present study (see Fig. 6). To see that, the time-average of 0"23(x,, t) since the Jast crack jump, denoted as a23(x,' T' � Ilk), is introduced:

(34)
The "true" propa_gation process is such that d"23(x,� r, nd = O'"c(x,) for every r. However. immediately after each jump �cc�rring du;ing the stepwise numerical simulation of the propagation, the stress 0"23(xr, t) in front of the new element drops to a value much smaller than O"c(x/); thus the time-average a23(x, ,nk+IL1t,nk) (34) at r=nk+tAt is smaller than O"c(x,). In this sense, the "simple" criterion does not adequately reproduce Irwin's criterion, especially for simulating the low-velocity initiation phase of the propagation. This remark remains equauy valid even in the limit Ax � 0 if ihe ratio � t/ Lix is kept constant, as is usuaHy done.

An improved version of the criterion, referred to as '"averaged", has therefore been established. Rupture propagation is initiated after the first time step at which (33) is fulfilled. However, a further crack tip jump is only implemented when the time-averaged stress (34) reaches the critical value, i.e. at a timet such that [START_REF] Nishirnura | Boundary integral equation method for clastodynamic crack problems[END_REF] Let this occur during the nth time interval: a new crack element is introduced either afte1 the (current) nth time step (nk+J =n) if t> (n-�)Lit, or after the (previous) (n-1)th time step (nk+t =n-1) if t�(n-!)Lit; in the latter case the (current) nth time step is recomputed with the new configuration of the crack.

Study of the onset of kinking

For the static anti plane crack, the stress component a 63(p, 6) varies like cos( 6 /2) (p being kept constant, (p, 9) denoting polar coordinates with origin B, 6 = 0 being the crack direction before kinking, see Fig. 2):

it has a maximum in the crack direction 6 = 0, so that a quasistatic crack under an tip lane shear should propagate along its initial direction. During dynamical propagation, however, it is known [START_REF] Freund | Dynamic. Fracture �Uechanics[END_REF] that the maximum of U93(p, 6, t) become flatter as v(t) increases and eventually occurs for a nonzero angle 9 (see Table 1); this may cause the appearance of a kink in the propagation path. Accordingly, in order to study the onset of kinking, the averaged stress O"rn(p., 0, t, nk) ( 34) is calculated fo r p=x, around the current crack tip B(t), the rupture criterion becoming [START_REF] Kuhn | A boundary integral equation method for radiation and scattering of elastic waveo; in three-dimensions[END_REF] If an angle OK =1:0 satisfies (36), OK is considered as a potential direction of kinking. In this case, a new crack element is added, in the direction OK, to the current configuration, according to the procedure developed for in-plane propagation. Further time steps are calculated for this kinked crack, until a new angle OK satisfying criterion ( 36) is found again. If OK� OK, the onset of kinking is quaiified as stable and the numerical simulation terminates; the p ropagation velocity observed just before kinking is then tentatively interpreted as a terminal velocity. If OK� OK, the kinking is assumed to be unstable; the oblique element is then replaced by a straight one and the corresponding time steps are recalculated.

Numerical results and discussion

Our numerical computations and the subsequent discussions are made using the following non-dimen sional quantities:

u u'=-, J.l t t'= ' J.l K'=____£_ pJiiL ' ..I., = 1!. 'Y L ' c t'=t- L ( 37
)
where L is a characteristic length of the problem -either the total length of an initial finite crack or the length .dx of a boundary element for a semi-infinite crack (where no physical characteristic length exists).

For propagating cracks, I (t) and v(t) will respectively denote the arc-length B(O)B(t) and the crack tip velocity.

Unless stated otherwise, the crack faces are assumed to be stress-free, while the only nonzero component of arer is ( uli)' = 1; thus the superposition method leads to consider t3 = 1 in the BIE (9). The geophysically more realistic case of a stress drop to a constant dynamic friction level could in principle be achieved by superposing a constant stress field. However, a careful look at what happens physically if v(t) becomes negative would then be mandatory.

Das and Aki [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF] introduced a dimensionless quantity S which normalizes the difference between the mean stress within the element which is closest to the crack tip and the homogeneous fa r-field stress cr21

with respect to the stress drop, in our case equal to the far field stress itself. In the case of elastic, i.e. ideally brittle crack propagation, as treated in [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF][START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] and in the present paper, this quantity S (called Tu in [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF])

is solely an indication of how fine the numerical grid has been chosen. S becomes physically significant only if a cohesive zone is introduced at the crack tip, as in [START_REF] Andrews | Dynamic plane-strain shear rupture with a slip-weakening friction law calculated by a boundary integral method[END_REF]. In the present context, S is given, using [START_REF] Segall | Mechanics of discontinuous fa ults[END_REF] and the notations of 4, by 1 + S= .j2u�(x, = .dx/2) = 2K;.

Stabilization of the solution

The first attempts to use our numerical method produced oscillations of the numerically computed COD with respect to the space and time variables. The amplitude of these oscillations increases along with time.

Figure 3 (slip history of an instantaneously appearing stationary crack, modelled with 20 boundary ele ments) shows an example oi such oscillations, where the spatial frequency is the highest allowed by the mesh (note: in all figures, the slip on the crack-tip elements shows a seemingly linear variation, but the actual computations use indeed the square-root variation). This phenomenon gets worse as more boundary elements are used.

At first sight this phenomenon could be thought of as a manifestation of an instability of the time rrtarching scheme, i.e. with respect to the time variable. However, the fa ct that the numerical COD oscillates with r espect to the space va r iable rather suggests a spatial instability. This has been confirmed by imple menting a static version of the traction BIE using the same space discretization: the static COD comput�d with the static BEM program showed the same kind of oscillations. (

) 40
Thus, as the J collocation points (for C=[-1, !] and �x=21J) are sx=s;-112 =(2i-J-l)/J, i= l, ... , J in our approach, one has, taking a= Jk1t (k integer) in ( 40):

.. .tf [cos Jkns, S;-112] = o( 1 I Jk1t).

(

) 41
In other words, the perturbation A cos Jkrr.s added to cp '(s) produces almost zero tractions at the collocation abscissae S;-, 12, and is therefore likely to pollute the numerical solution of (9) or (39). Furthermore, in view of ( 41 ) , , a finer mesh leads to increased pollution, consistently with what was observed. Such behavior is apparent m Fig. 3 (with k = i ). Had the eiement nodes been taken as collocati0n points (Le. sx = s;= (2i-J)/J� i= 1, ... , J -1), eq. {40) would give. with a =Jkrr./2, .Jt'[sin 1Jk1ts, s;]=o(2/Jk1t) [START_REF] Tikhonov | Solutions to Ill-posed Problems[END_REF] and hence lead to the sar.1e kind of instability.

The suggests that the static traction BIE, collocated at regularly spaced points, �eads to spatia] oscillations of th(! numerical COD. and hence that the observed oscillations are of a spatial �:ature.

The initial problem of numerically solving the BIE (9) has then been reformulated as a penalized least squares problem, using the stabilizing fu nctional (in the sense of Tikhonov [START_REF] Tikhonov | Solutions to Ill-posed Problems[END_REF]) : J ,02 1 2 L f 10 1 2 2!2( 1/J , :. ) = a , L3 c !&2 t/> '(s, t n ) ds y + a 2 C c O t t/> ' (s, t .), ds y [START_REF] Trifu | Predicted near-fi . eld ground motion for dynamic stress-drop models[END_REF] in which the scaling (37) has been used. The fu nctional D(</J, t) (43) allows smoothing with respect to space (using a. :;':0) or time (using a 2 �0) in order to numerically test both possibilities. Its discretized expression is [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] and the discrete least squares problem ( 44) becomes

"t'l" umerical tests using (45) with various values of the non-dimensional coefficients a1 • a2 have been per fo rmed. The results show that the oscillations of the numerical COD disappear by taking a 1 > 0, a 2 = 0 but persist if a,= 0, a2 > 0 is used instead. This is a confirmation of the above analysis. In all the fo rthcoming numerical results, the penalty method has been usi!d with a 1 = 0.5L1x 3 1 L 3 = 0.5/(Jo)\ a2 ::.:0; a, is a regula;ization parameter in the sense of Tikhonov, smali compared to unjt y .

Comparisons with analytical results

Stationary crack

A stationary crack of length L appears suddenly at t = 0. This problem has an exact solution u(y, t), a(y, t) fo r t�t1(y), t1 (y) being the time at which the information of the finiteness of the crack reaches y [START_REF] Kobayashi | Initiation, propagation and kin king of an antiplane crack[END_REF] (tj(y) = 1 at the edges) ; and then again in the static limit [41]. Figure 4 shows the numerical results for the stress at near-tip points (p, 0) with p=0.01L and 0=0°, 15°) ... , i5°, 90°, computed using 50 boundary elements; they are within 1 o/o of the exact solution, shown by continuous lines. The slip and the f!:l �tnr nf .fl,rn<:l""•"'Ql ,.. , .:. o::o h ,.. ,.. f ro.f' *�"'"' �Tr '""'"" n"'+ .. .. . 1 .... .. . A 1-f11\ ,. .. ...,. .... . L ... . :_ ,,.l ... :.t... :-10/ "' """''' ... """" •• ••�=-� •"'-l .... .. .. '-'• ... J ...... ... . 1"".a v ..." "f..,.aavvL va Lu"' ..., ... . \'-'"'a."" va.au"' . �1 "' L'Jl a.a"' VVLCUU\,.. \.l VVli.UU l 1 /U """''-'Ul Cl'-'.)' u.-,1u0 only 20 elements ; in this case the stress at near-tip points (p = 0.01 L) is overestimated by approximately 5°/o . As a general trend, the f< �wer elements, the m.ore overestima�ed the stress, and to a lesser extent the slip.

Semi-infinite crack moving at prescribed velocity

A semi-infinite crack suddenly appears along the negative x1-axis and then propagates at a kinematically imposed constant rupture velocity v. The exact solution to this problem, fo und by Kostrov [31], has become .f£_

- Analytical 6 Numerical 0 LJ----�--�----�----�--�----� 0 3C 60 90
Angle Viewed From Crnck Tip (Degrees) Fig. 4. Shear stress u83 for an instantaneously appearing stationary crack around its tip, at a distance cf 0.01 times its total length (50 boundary elements).

a standard benchmark, which has been used in both [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF] and [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] . Our numerical results for the displace ment and stress history fo r a crack that instantaneously appears along x1 < 0 and extends at half the shear wave velocity are compared with Kostrov ' s exact solution [START_REF] Kostrov | Unsteady propagation of longitudinal shear cracks•[END_REF] and the FD results of Virieux & Madariaga ([44 ] , Fig. 2(a,b)) on Fig. 5. Kostrov's and our values are calculated at the points actually indicated, whereas the results of [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] are mean values over one grid element. Our numerical solution agrees very well with the analytical one. Since both slip and stress could not be computed simultaneously on the crack plane in L44], slip was evaluated for positions slightly off the plane. Its deviation from the exact solution in Fig. 5(a) might therefore be interpreted as solely due to this minor handicap. On the contrary, the stress diagram (Fig.

IJ

value should be J2 times larger than the actual stress values (computed at the grid midpoint) which are displayed fo r both Kostrov's and our solution. From this remark, it becomes obvious that the stress singularity is significantly better resolved by the present BIE method. It is worth noting that, owing to the self-similarity of the problem at hand, using a finer FD grid would lead to a rescaled but otherwise identical solution and could therefore not compensate fo r the discrepancy.

The accuracy of the slip values obtained in [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF]Fig. 6] does not seem to differ significantly fr om the present results. However, [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF]Fig. 7] shows a comparatively poorer stress resolution.

Spontaneously propagating crack

The only problem of this kind having a known analytical solution is the case of a suddenly appearing straight semi-infinite crack under antiplane shear (Kostrov [31]). The crack tip location /(t) has been computed using the present BIE method ; our results are displayed in Fig. 6, together with corresponding numerical results from [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF], and compared to Kostrov's solution. The latter has been calculated with slightly readjusted stress intensity factors Kc of 0.50, 1.85, 3.93, 5.90 respectively (for the reason explained in Section 4.2), chosen to rninimize the difference in slope towards the right-hand end of the parts of the curves depicted in Fig. 6. The present numerical evaluation of the early acceleration phase is excellent for K; �2 and if the averaged rupture criterion is used. Fairly good (K; = 1) to excellent (K� � 2) agreement is found when the averaged rupture criterion (see Section 4) is t1 sed. On the contrary, the simple rupture criterion makes the crack accelerate much too ouicklv_ as it is aonarent for the case K� = 2 in Fig. 6. K; = 1 should be considered as a lower bound fo r ;�� d

-;e • s�lts, whe;�as K; = 6 has been ch�sen as an �pper limit fo r the present situation of an infinite crack :

the time of fracture initiation grows with u<;) 2 , and this increases not only the number of time steps but also the number of spatial elements needed to siluulate effectively a semi-infinite crack (the K� = l and K; = 4 runs of Fig. 6 need respectively 2 and 15 minutes on a Cray 2 computer) .

In [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF], satisfactory results were obtained only for 2�K; �4; the curves fo r these limit values are included in Fig. 6. The poor stress resolution of the FI) technique leads to a delayed fr acture initiation. This delay is later compensated by a faster acceleration due to the application of the simple rupture criterion. The results for the same problem in [ 17, Fig. 13] are only evaluated fo r very short distances of rupture propaga tion. A quantitative comparison of accuracies is hardly possible using these results, particularly in view of the parameter adjustments involved. However, at least the curves fo r S=4 (K �=2.5) and S=5 (K;=3) in [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF]Fig. 13] show clearly a tendency to accelerate too quickly, probably because of the use of the simple rupture criterion.

Spontaneous growth of a fi nite crack

Following [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF], two kinds of rupture nucleation for a crack of initial length L are considered here : either the crack appears suddenly at t = 0 within a homogeneous stress field, or it is first loaded statically so that the critical stress criterion is just reached. In the latter case the crack is in metastable equilibrium, and the rupture is started by imposing a slight increase of its length, one element being added at the tip B at t = 0. This is physically more satisfactory but requires more computer time since the initial prestress field has to be calcul(l �ted up to a time at which it can be considered as static.

Figure 7 shows I ( t) fo r a finite crack unilaterally extending to twice its initial length. Our results are compared with analogous ones from [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] ; in both studies the initial crack was modelled with 40 elements or grid points, respectively. Although only two-sided symmetrical propagation was studied in [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF], a one sided extension would have given exactly the same crack tip propagation within the time interval displayed in Fig. 7. For the instantaneous crack, Kostrov ' s solution for a semi-infinite crack is included, remaining valid up to t' � 1.3 when the information of the finite length reaches the running crack tip ; after that, the numerical solutions slowly start to Jag behind, as expected. In the case of the static precrack, our solution shows a very slow acceleration phase over a first short running distance ; this is in complete contrast to Virieux and Madariaga's results. 4

The SIF of the static prestress actuaJly obtained by the two numerical methods w�re respectively used as the critical value Kc. With the present normalization, this was K;=0.485 in [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] , 0.519 here, the corre sponding analytical value being /(�= 0.5 (a small deviation in the direction of our result is expected based on the approximate calculation of I<: from the stress as explained in Section 4 and is not a direct measure of computing precision). Renormalized fo r a semi-infinite crack (multiplied by J Lf t1x), the analytical value would be K; =3.16, i.e. within the best range for the approach of [START_REF] Virieux | Dynamic fa ulting studied by a finite difference method[END_REF] , see Section 5.2.

In the case of rupture propagation fr om an instantaneously appearing crack, both methods seem to give roughly equivalent results, although the present one resolves more accurately the early but short phase of noticeable acceleration. However, the situation drastically changes when looking at the physically more realistic rupture initiation from a metastable static precrack. At the time when an additional crack element is added to the static crack� the shear stress in front of the new tip is much below the critical stress u c. Within only a few time steps, the stress will ultimately become larger than Uc , but only slightly, since the added element is small compared to the initial crack length. It therefore takes many additional time steps until the initial stress deficiency is compensated and the fr acture criterion [START_REF] Nishirnura | Boundary integral equation method for clastodynamic crack problems[END_REF] satisfied. The slip of the crack nucleating fr om an initial static situation is shown in Fig. 8(a) at eight equally spaced instants t ' = 0.5, 1, ... , 4 (the rupture initiation taking place at t ' = 0), for the crack studied in our curve s, Fig. 7. The crack is arrested once it has doubled its initial length. Note that the slip is measured from the initial state at t = 0 since only this would be detectable from purely seismological measurements.

The stress field u23 in fr ont of the actual crack tip is shown in Fig. 8(b) at the same equally spac-ed times than in Fig. 8(a). The heights of the stress peaks vary because the times at which they are depicted are accidentaUy related to the insta1:1t� f'f crack tip jumps. Nevertheless, it becomes evident that the rupture criterion [START_REF] Nishirnura | Boundary integral equation method for clastodynamic crack problems[END_REF] is satisfied by relatively low maximum strc�s values acting during a long time interval in early phase and by quickly growing stress acting for a shorter time once the crack moves fa ster.

Figure 9 shows. I (t) for a few examples of one-sided crack propagation in inhomogeneous media. In all cases, rupture nucleates from a metastable static precrack of length L embedded in a homogeneous stress field ( u2{)' = 1. The initial crack is modelled by 20 boundary elements. This turns out to be suffi cient since the numerical Kc corresponds to a value K; = 0.4 7 4; this in turn would give a value K; = 2. 12 in the semi infinite case, which lies within the range of good results.

• Curve 1 corresponds to fully homogeneous conditi�ns and therefore represents the identical physical situation as our s curve in Fig. 7. The difference is that only half as many boundary elements as in Fig. 7 are used here. Since rupture is nucleated by additioat of cne new boundary element in fr ont of the static tip, the initial perturbation is now twice as big. The consequence is that the cr��k accelerates fa ster at the very beginning. However, fo r 1 ;::. 0.2L, the c1rack tip lvcation histories in Figs. 7 and9 are virtually identical, the total running times between /= 0.2L and /= L being equal within about ± 1°/o . Thus the results are independent of the element size, except for the effect of the element size dependent rupture nucleation which manifests itself only over a short propagation distance. • Curve 2 represents a crack propagating within a homogeneous stress field and through a barrier of higher material resistance (the critical stress O'c being J2 times larger for L�l� 1.2L than in the surround ing region), the elastic constants being unchanged. A quasistatically propagating crack 1reaching this barrier would find itseif in a new metastable equilibrium. The dynamical crack, however, breaks through with a low ili•::an velocity of less than 0.2c. Once ' �, ne barrier is broken, v(t) seems to imm,�diately take up the value co1 responding to the case without barrier.

e Curve 3 shows a crack running through a material of homogeneous toughness and entering a decreasing stress environment after having reached the length 2L, in which ( <r2{)' decreases linearly fron1 1.0 to 0.0 betwt�en the positions/= L and I= 3L; v(t) decreases for /"?;:. La nd the • ;--� .. :.;: ultimately stops at /= 3.05L (its stopping near the point where the stress vanishes being due to a coincidence in the choice of the parameters).

• Curve 4 shows the combined effect of barrier and decreasing stress field. Owing to the barrier, the crack entering the dec:reasing stress field loses some of its 'dynamical impetus' and ultimately stops after a slightly shorter running distance than in case 3.

The slip �long the crack and the stress straight in front of it associated to curve 2 of Fig. 9 Figures 10 and 11 reveal the main weakness of our numerical solution of crack growth : the stepwise propagation produces high-frequency waves, which would indeed physically exist if the real propagation process were stepwise C�pseudo-physical waves " ). These waves have caused no real problem so fa r, but wiH turn out to be rather annoying when the onset of crack path kin king is studied.

4. !Vumerical study of the onset of kinking

A dynamically propagating crack under antiplane shear may experience kinking (see Section 4.2) if 693(p, 6, t) is maximized (for fixed p) at a nonzero angle OK(t). The numerical determination of OK(t) will bt: highly sensitive to the high-frequency disturbances, for they reach the stress calculation points (p, 0) at ti1nes which depend on 0.

Here the kinking of a crack nucleating from an initially static precrack modelled with 20 elements (case 1 of Fig. 9), is considered. After each time step, 6'tn(p, 0, t) is computed fo r p = Ax/2 and e = OC , 15°, ... , 75o, 90°. Once v(t) has reached a value of 2c /3, a fe w kinking attempts were observed, the n. pture criterion being satisfied at (} = 15°. However, it was only when reaching v( t) = 6c /1, after a running d\stance l= 3.5L, that the kinking turned out to be "stable" (in the sense of the procedure described in Section 4). This occurred for 8 =45 ° or 60 ° (a choice between these two angles would be based on unsignific,ant differences in the value of 60 3 (p, 8, t)). The same study for a suddenly appearing crack and using 40 • :!lements did not introduce a significant change in these results. The time elapsed between the instant of kinking and the instant at which [START_REF] Kuhn | A boundary integral equation method for radiation and scattering of elastic waveo; in three-dimensions[END_REF] is satisfied in front of the kinked element in our numerical calculations suggests that v(t) does not slow down significantly due to kinking, at least at its onset.

No analytical solution is known which exactly corresponds to the kinking problem at hand. The symmet ric bifurcation of a running crack was treated by Achenbach [START_REF] Achenbach | Bifurcation of a running crack in antipiane strain[END_REF], whose results were corrected later by Dempsey et al. [START_REF] 0-.. Rnpsey | Dynamic effects in mode Ill crack bifurcation[END_REF] . :Bifurcation was found to occur at v( t) = 0.6c with an angle of 90°; immediately after branching, v(t) was extremely low (about 0.02c), but was rapidly growing ..

Experimental results, though done in mode I (see [START_REF] Dcnpsey | Dynamic crack branching in brittle solids[END_REF] or [START_REF] 0-.. Rnpsey | Dynamic effects in mode Ill crack bifurcation[END_REF] fo r references), usually show kinking angles from 10° to 45° at v(t) �0.5c to 0.6c. The crack velocity immediately after kinking generally shows little or no change from the pre-branching value. The present numerical results are thus not farther away fr om what is expected in experimental evidence than the conclusions of the above-mentioned �nalytical works. However, the kinking velocity seems to be highly overestimated here.• Figure 12 shows how the stress a 93 in front of the "running" crack tip typically builds up, frmn time step to time step after a crack tip jump, in the case of one spatial jump per three time steps (i.e. v = 2c/3). The ratios of out-of-plane ( 8 # 0) to in-plane ( 8 = 0) stresses show oscillations of magnitude much higher than the small differences to 1 that are ex}A �ted to occur fo r them (see Table 1). Since the time-dependent stresses are properly evaluated around a fixed crack (see Fig. 4), the inaccuracies of these ratios are certainly a consequence of the high-frequency perturbations induced by the stepwise numerical crack propagation. This explains, in our opinion, the present unsatisfactory numerical value of the kinking velocity. This situation is not improved when finer meshes are used as long as �t/ �x remains constant ; moreover, using several values of At I� x would require too many calculations. This situation, however, is expected to improve considerably once a cohesive zone [START_REF] Ida | Cohesive fo rce across the tip of a longitudinal shear crack and Griffith's specific surface energy[END_REF] will be introduced . . .. . . .

• 6

First time step We are aware of the fa ct that the kinking results presented here are tentative. Other criteria, particularly energetic ones, would have to be tested too, but this was beyond the scope of the present work. Furthermore, as pointed out by Freund [START_REF] Freund | Dynamic. Fracture �Uechanics[END_REF], it might be impossible to establish a realistic branching criterion without taking into account a finite process zone.

Conclusion

We have proposed a new boundary integral equation for solving two-dimensional transient dynamical crack problems of general geometry. This method is based on a regularization of the time-domain traction BIE of the crack problem. The singularities of the traction BIE are removed by a technique proposed by Bonnet [START_REF] Bonnet | Methode des equations integrates regularisees en elastodynamique tridimensionelle[END_REF][START_REF] Bonnet | Regularized boundary integral equations for three-dimensional bounded or unbounded elastic bodies containing curved cracks of arbitrary shape under dynamic loading[END_REF] and implemented here for the particular case of antiplane loading. Special crack tip elements similar to those used by Hirose and Achenback [START_REF] Hi Rose | Time-domain boundary element analysis of elastic wave interaction with a crack[END_REF] were used in order to improve the stress field resolution near the borders of the crack. Several problems of interest in elastodynamic crack theory were solved with our method and the numerical results obtained were compared to solutions obtained by other numerical methods and with analytical solutions in the few cases where these are available.

Our goal is to develop a stable numerical method for solving dynamical crack-growth problems with a view to applications in seismology and metallurgy. Classical methods like finite differences or finite elements are difficult to apply to the study of crack growth because of their intrinsic numerical dispersion of high frequency elastic waves. Integral equations, on the other hand, give much better results as shown by [START_REF] Andrews | Dynamic plane-strain shear rupture with a slip-weakening friction law calculated by a boundary integral method[END_REF] and [START_REF] Das | Earthquake Source Mechanics[END_REF]. However, the classical fo rmulntion used by seismolog�_,ts, proposed by Das and Aki [START_REF] Das | A numerical study of two-J.imensional spontaneous rupture propagation[END_REF], can only be used to solve problems with axial symmetry, so that only flat collinear cracks may be studied. In seismology, however, the most important problem is understanding the role of crack interaction and geometry in controlling rupture growth. The method proposed here can deal with several cracks of general shape that dynamically interact and is therefore a first step in the direction of studying complex crack growth.

A difficulty that has been only partially solved here is that of the estimation of the stress field directly in front of the rupture front. In our fo rmulation the crack jumps by a grid space whenever a simple stress rupture criterion is satisfied, often attributed to Irwin in the seismological literature. In future we plan to introduce a more realistic slip-or velocity-weakening criterion ( [28] or [J 8]) in order to calcnl�te fr�«:'t \l rr.: energy and establish a criterion based on more appropriate energy dissipation considerations. The same problem arises fo r the computation of hoop stresses in frr ,t of the rupture f:-ont in studies of crack branching and kinking. In this case, we believe that results will improve significantly if we use collocation of the BIE at points closer to the rupture front in the crack-tip element.

( 3 )

 3 where e,p denotes the two-dimensional permutation symbol. The curvilinear abscissa along C is denoted by s, with s = 0 for the end point A. The faces of the crack are loaded by tractions T +, T-(the crack faces are physically considered free of tractions, rut one often puts T± = -CT'i;n, in order to consider the perturba tion induced by the presence of the crack on a reference elastic state ,rer , arer known everywhere in .C2).

Fig. 1 .

 1 Fig. 1. Two-dimensional crack problem: geometry and notations.

•

  The tir.st �tep is due to Sladek & Sladek [3R]. The inner product of the representation formula (5) of O"J r( Z , t) with n ± (x) (for fixed z� C and xe C) is tc:ken. Then the resulting equation is integrated by parts with respect toy (no singularity occurs at this stage). Finally, in order to relate the given tractions T± on C to the unknown COD l/J, the limiting process (z�C) � (xe C ) is carried out, giving as a result T±(x, t) = -pni{I* f' I,.. (x, y, t-r)e,p � l/J(y, r) d r d s y

  ,t)--pnp(x) Ir(x,y,t-r)e,p-l/J(y,T)dr-Ir(x,y) erp-tf>(y,t) ds y

  y)dsJ.=--[r, Logr]�+ -------:\ Logr n,(y)ds).

  c-.1C(.r> r os which serves as a basis for the numerical work presented in this paper.

3 .

 3 Boundary element discretization for cracks in antiplane strain

  with R"" t:: '!: m = HGR '!:m + HGS '!. • m �,lj lj lj

  new right-hand side has to be computed at each time step, involving all the past values of the COD.

Fig. 2 .

 2 Fig. 2. Polar coordinates emanating from the propagating crack tip: (a) b�• Jre kinking, (b) after kinking.

  ratio a03(0)/cr113(Q), against e and vjc •�-------

Fig. 3 .

 3 Fig. 3. Numerical solution without smoothing for an instantaneously appearing stationary crack (20 boundary elements).

S

  (b)) reveals a fundamental inaccuracy of the FD technique used in [441 . The seen1ingly excellent stress reso�ution is illusory bece!use the values shown in Fig. 5(b) ar� mt"'�n values over the FD grid inten:alclosest to the crack tip. Howe;:ver, assuming an approxin1ate stress distribution (eq, (28)) there, this rnean

Fig. 5 .

 5 Fig. 5. Disph!r•-ment and stress history for a crack that instantaneously appears al0ng x1 < 0 and extends at ,_, = c j2.

Fig. 6 .

 6 Fig.6. Comparison of exact[START_REF] Kostrov | Unsteady propagation of longitudinal shear cracks•[END_REF] and numerical ({44] and present) solutions for the crack tip location during the spontaneous rupture propagation of a semi-infinite crack.

1 : 3 Fig. 7 .

 137 Fig. 7. One-sided spontaneous rupture fr om an initial finite crack (40 boundary elements).

2 '

 2

Fig. 8 .

 8 Fig. 8. Slip and stress fields for tl = 0.5, ti = l, . .. , t8 = 4, for the case of our curve s, Fig. 7.

Fig. 9 .

 9 Fig. 9. Crack tip location for one-sided spontaneous rupture propagation, in different environments as indicated.

  Figures ll(a), (b) show respectively the slip and the stress fo r curve 3 of Fig.9{decreasing ambient stress field). The crack decelerat"'s then stops ; the stress concentration at Ba t the end of the propagation is as strong as the initial static concentration.

Fig. 10 .

 10 Fig. 10. Slip and stress histories of case 2, Fig. 9 (crack breaking through a barrier).

Fig. I 1 .

 1 Fig. I 1. Slip and stress histories oi case 3 in Fig. 9. i.e. of a crack running into a zon� of decreasing ambient stress field.

  . ---. A -. •• -•-to. •. 6••••••••-•• •••6 •••-•----••6•• • -•• ... •••• • A • • ••• : • : - • ::: _.6 • • • 6••••• .. � Second time step A .

Fig. 12 .

 12 Fig. 12. Shear stress development around the "nmmng" crack tip. from the time step to time step after a crack tip jump. l he ratio rr1nW) cr,n(O) does not c' evelop smoothly.

  t) dsy +

	.1C(.r)	c	os	os	.1Q.r>	2 r	-(y, t) --(x, t) dsy os os

i ' r f' ..

dr -2 rp(YJ '• </J(y, t-r) 2 2 2 112 dsy c j c r/c (t' -r jc)

Acknowledgments

The support of the first author by the Swiss National Science Foundation and the free use of the infrastructure at the Laboratoire de Geophysique Interne et Tectonophysique (LGIT), Grenoble, are gratefully acknowledged. All numerical results presented here were obtained using the Cray 2 supercomputer of the Centre de Calcul Vectoriel pour la Recherche, Ecole Polytechnique, Palaiseau, France.

The authors are indebted to Prof. Jean Virieux, of Universite de Nice, France, fo r helpful discussions. The stimuiating atmosphere of mutuai heip at the LGiT has been of great heip fo r the present work. They also wish to thank the referee fo r his careful review of this paper and helpful comn1ents.

Appendix A. Two-dimensional Green tensors for aiiltiplane problems

The nonzero components of the impulsive displacement and stress Green's tensor and of the static stress Green's tensor, fo r the infinite space under antiplane strain, are given by

respectively. From eqs. ( 47) and (48), the analytic expression of convol J�ion integrals which appear in the time-domain BIEs can be de�•ived as fo llows [22] :

j U(z, y, t-r).f (r) ctr = 2 1 ( 2_ 2/ 2)1.�2 dr f Er(z, y, t-r)f( r) dr 0 (50)

(5 1)

If one inserts f( t) = S ( t 1 L1 t-m) in (51) and /( t) = S( t I L1 t-m) in (50), the integrals with respect to r have the fo llowing analytic expressions:

f' "

= [[r-rfct }r jc+ A+( r)r: �' >"' -[ [ r -r /et +1 r jc+ A+( r) J:t Mr Appendix B. Evaluation of the elementary integrals HGS In order to evaluate HGS� and HGS �, two auxiliary integrals / iJ and JiJ are considered : Consideration of all the possible cases for v'1(s) in (55), (57) leads to :

One has v• • i (s)= I , v' ; -1 (s) = -l, v' i (s) =O (j ¥= i-1, i). Thus

!iJ =O

the latter being a consequence of property ( 12) of the basis fu nctions v1 .

• i= i.

(

(61 )

.j2 --= --Log( I + ,/2) (63)

L\x Jlj = O, j'�2.

(64)

A similar calculation givr�s

Finally, from eqs. ( 66) and (68) :

HGS�= 'y(X) r ,l:'x _� X y [[ vj (s) -vj (S;-1/2)] + v'j(s)[S( -r /d t) -I]] ds,,

•.l p, ,.

"'

(66) (i, j) :F(l , l ), (J, J-1 ) .

(70)