Dynamic fracture mechanics with contact interaction at the crack edges

Engineering Analysis with Boundary Elements 24 (2000) 643±659

www.elsevier.com/locate/enganabound

Dynamic fracture mechanics with contact interaction at the crack edges V.V. Zozulya*, P.I. Gonzalez-Chi Centro de InvestigacioÂn CientI®ca de YucataÂn, A.C. Calle 43, No 130, Colonia, Chuburna de Hidalgo, MeÂrida, 97200 YucataÂn, Mexico Received 11 January 2000; received in revised form 13 July 2000; accepted 14 July 2000

Abstract This paper deals with elastodynamic contact problems with unilateral restrictions and friction for bodies with cracks. General dynamic loading and a relevant case of harmonic loading are considered. The mathematical aspects of this problem are brie¯y discussed. The variational formulation of the problem, variational boundary inequalities and boundary functionals are used to solve it. The problems of a tension-compression plane harmonic wave interaction with one and two co-linear cracks of ®nite length with the allowance of unilateral contact interaction at the crack edges are solved. The in¯uence of the contact interaction between the crack edges on a stress intensity factor is investigated. q 2000 Elsevier Science Ltd. All rights reserved. Keywords: Contact with friction; Crack; Boundary integral equations; Stress intensity factor

1. Introduction Almost every material used in engineering contains cracks. These cracks can arise during the material preparation, parts machining, structure manufacture or maintenance. Under load action, the existing crack surfaces in the material form interfaces where unilateral contact forces and frictional contact forces interact at the normal and tangential direction, respectively. The unilateral contact and the friction are non-linear phenomena because the boundaries between contact and non-contact regions and also between adhesion and sliding zones are a priori unknown. The conditions for unilateral contact and friction are mathematically expressed in the form of inequalities. The unilateral contact problem for an elastic body with rigid support was formulated by Signorini in 1933 [16]. The ®rst strong mathematical statement of this problem was done for elastostatics using variational inequalities by Fichera [20]. Since then, many other problems with conditions in the form of inequalities in mechanics and physics have been investigated from both mathematical and applied point of view. For the mathematical aspects of the problem, we refer to Refs. [16,32,38]. The algorithms and numerical solutions were studied in Refs. [1,21,33,38]. For publications and surveys in this topic see Refs. [28,42]. The problems with inequalities are usually solved numerically using ®nite [1,21,33] or boundary [1,2,27] * Corresponding author. E-mail address: [email protected] (V.V. Zozulya).

element methods and some kind of iteration process. If boundary conditions take the form of inequalities, as it happens in two-dimensional and three-dimensional contact problems, the boundary element method has the advantage of numerical implementation. Which means that ®nding the unknown boundary conditions solves the problem. Part of the boundary conditions has the form of inequalities; because of this the solution is reduced to the application of the boundary element method and the iteration process mentioned above. The boundary element method is based on the transformation of the boundary value and the initial boundary value problems into boundary integral equations. Previously, this method was used to investigate theoretically the existence and uniqueness of a solution [30,36] but, in recent years with the aid of computers, it has become one of the most effective and powerful numerical methods for the solution of scienti®c and engineering problems [4,5,9]. For the application of the boundary integral equations method to elastodynamics and fracture mechanics refer to Refs. [4,5,7,8,15,27]. One of the main problems arising when the boundary integral equations are used, is the singularity of the integral operator kernels. These kernels may be classi®ed as weakly singular, singular and hypersingular. The boundary integral equations with weakly singular and singular kernels were studied in Refs. [4,36]. The hypersingular integrals were introduced by Hadamard [31] and now are widely used in various problems. The application of the hypersingular integrals in fracture mechanics was studied in Refs. [27,48,57±59].

0955-7997/00/$ - see front matter q 2000 Elsevier Science Ltd. All rights reserved. PII: S 0955-799 7(00)00029-1

644

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

Nomenclature V Volume of an elastic body I Time interval 2V Body boundary 2Vu and 2Vp Boundary parts, where the displacements and surface traction are assigned u Displacements p Surface traction b Body forces N Number of cracks V n1 and V n2 Opposite edges of the crack e ij and s ij Strain and stress tensors Aij Elasticity differential operator for displacements 2 i ˆ 2=2xi Derivative with respect to coordinates 2t ˆ 2=2t Derivative with respect to time Elastic modules tensor cijkl r Material's density l and m Lame constants d ij Kroneker's symbol u0i and v0i Functions that describe the initial conditions c i and w i Functions that describe the boundary conditions on the parts 2Vp and 2Vu Du Displacements of the crack surfaces n 1 and n 2 Unit vectors normal to the crack surfaces q Forces of the crack edges contact interaction qn and qt Normal and tangential components of the contact force vector Dun and Dut Normal and tangential components of the displacement discontinuity vector h0 Initial crack opening kt and l t Coef®cients dependent on the contact surfaces properties Ve Close contact area L Laplace transform Inverse Laplace transform L 21 k Parameter of the Laplace transform F Finite Fourier transform Finite inverse Fourier transform F 21 Uij, Wij, Kij and Fij Elastodynamic Green's functions in the Laplace transform space U ik ; Wik ; Kik and Fik Elastodynamic potentials in the Laplace transform space c1 and c2 Velocity of dilatational and distortional waves Kn …li † McDonald's functions Hn …li † Hankel's functions Differential operator that transforms a displacement into surface traction Pik Differential operators D^ ik Integro-differential operators T^ ik Hl;m and Hrg Banach functional spaces g bi Maximal monotonic operator Superpotential ji bci Maximal monotonic operator conjugated to b i jci Superpotential conjugated to ji 2 Subdifferential of convex analysis F i …ui † Non-smooth functional Ku …ui † and Kq …qi † Sets of one-sided restrictions with friction C…u i ; qi † Boundary functional Pn[qn] and Pt [qt ] Operators for orthogonal projection on the sets of one-sided restrictions with friction r n and r t Coef®cients used for the best convergence of an algorithm P t, P k, Puq and Ppq Operators of projection into ®nite-dimensional functional spaces

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

645

Xgq and Ygq Finite-dimensional functional spaces Nu, Np and Ne Number of the boundary elements on 2Vu, 2Vp and V , respectively l1 and l2 Crack lengths h1 and h2 Distances between the crack centers v and T Vibration frequency and period c0 Amplitude of the tension-compression wave a angle of the incident wave E and y Young's modulus and Paisson's ratio Wave number k 1p KI and KII Stress intensity factors Mathematical aspects of the problem under consideration and the algorithms for its solution have been investigated in Refs. [10,16,17,21,32,37,38]. We have developed a different algorithm which consists in ®nding a saddle point of a subdifferentional boundary functional [46,53]. It was shown that the algorithm may be considered as a compressive operator acting on special Sobolev's functional spaces [32,40,41]. The convergence of this algorithm was proved in Refs. [46,56]. This algorithm solves an elastodynamic problem for bodies with cracks without taking into account unilateral restrictions, followed by the projection on sets of unilateral restrictions and friction. The elastodynamic problem for bodies with cracks with no unilateral restrictions is solved using the boundary integral equations in the Laplace transform methods. Then, the projection operators on a set of unilateral restrictions are constructed. The character of the kernel singularities for boundary integral operators is studied. Two regularization methods for integrals with ªstrongº singularities are considered. The ®rst method is based on the reduction of the boundary integral equations with ªstrongº singularities to integrodifferential equations. The second one is based on ®nite part integrals in the sense of Hadamard. The solution analysis for static fracture mechanics demonstrates that taking into account the contact interaction of the crack edges may signi®cantly result in fracture mechanics criterions. For fracture dynamics problems the effects occurring may considerably exceed those in statics. Moreover, in fracture dynamics, ®nding the loads that will not cause the crack edges contact interaction is a more complicated task. This problem is also essential for harmonic loading, when the steady-state regime with the harmonic time dependence of the stress±strain state is considered. For instance, the interaction of a plane harmonic tension-compression wave with a ®nite length crack is solved using such assumptions [18,22,39]. In these and other studies, it was pointed out that such an approach is incorrect; under a compressive wave action, contact interaction of the crack edges always occurs. Despite the fact that practically every researcher in the ®eld of fracture mechanics agrees on the necessity to take into account contact interactions at the crack edges, there were no studies about such problems [11]. They were investigated for the ®rst time in Refs. [24±29,45±47,49±55].

The case of harmonic load was considered in Ref. [47], where it was also shown that if contact interaction at the crack edges is accounted, the harmonic load results in a steady-state periodic process, not harmonic. The investigation of the contact interaction at crack edges with a ®nite length under harmonic loading in a plane was carried out in Refs. [26,49,51]. The in¯uence of contact interactions at the crack edges on the stress intensity factor for one crack was studied in Refs. [24,52] and for two collinear cracks in Refs. [25,50,54]. The technique evolved in these papers may be applied to elastodynamic unilateral contact problems [23]. The results obtained in previous articles are comprehensively and sequentially elucidated and some new results are also presented here. In particular, the mathematical statement and the solution methods are presented. The mechanical effects caused by contact interaction and their in¯uence on fracture mechanics criterions are investigated.

2. Statement of the problem Let us assume an elastic body in a three-dimensional Euclidean space R 3 that occupies a volume V. The boundary of the body 2V is piecewise-smooth and consists of two parts 2Vp and 2Vu, where the surface load vector p…x; t† and the displacements vector u…x; t† are assigned, respectively (Fig. 1). There are N arbitrarily oriented cracks in the body, which are described by their surfaces V n1 <

Fig. 1. Elastic body with N arbitrary oriented cracks under dynamic loading.

646

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

V n2 ; where V n1 and V n2 are the opposite edges. The body may be affected by body forces b…x; t†: We assume, that the displacements of the body points and their gradients are small, so its stress±strain state is described by linear elastodynamic equations [18,22,27]. We will relate every body point with a space point and the elastic body with volume V , R3 which will be occupied by body at moment t. Displacements and deformations of the body are described by a displacement vector ui …x; t† and a strain tensor eij …x; t†: They are connected by the Cauchy relations eij ˆ

1 2

ui …x; t0 † ˆ u0i …x†; 2t ui …x; t0 † ˆ v0i …x†; ;x [ V pi …x; t† ˆ s ij …x; t†nj …x† ˆ c i …x; t†; ;x [ 2Vp ; ;t [ I ui …x; t† ˆ wi …x; t†; ;x [ 2Vu ; ;t [ I

From this equation it follows that the strain tensor is symmetric and therefore eij ˆ eji . The components of the strain tensor must satisfy the Saint-Venant's relations: 2k 2l eij 2 2i 2l ekj ˆ 2k 2j eil 2 2i 2j ekl leading to 81 equations, but only six are independent due to the Euclidean space and strain tensor symmetries. The body stress state is described by a stress tensor s ij …x; t†. From the theorem of impulse moment balance it follows that the stress tensor is symmetric, and from the theorem of impulse balance it follows the equations of motion: 2j s ij 1 bi ˆ r22t ui ; ;x [ V; ;t [ I ˆ ‰t0 ; t1 Š If we consider an elastic body, the stress±strain relations have Hook's law form

s ij ˆ cijkl ekl In these equations, we introduce the following notations: 2i ˆ 2=2xi and 2t ˆ 2=2t; which are derivatives with respect to the space coordinates and time, respectively; r is the material density, cijkl are the components of the elastic module tensor, which for isotropic bodies has the form cijkl ˆ ldij dkl 1 m…dik djl 1 dil djk †

l and m are the Lame constants, dij is the Kroneker's symbol. The summation convention applies to repeated indexes. The above linear elastodynamic equations may be presented in the form …2:1†

The operator Aij for an anisotropic body has the form Aij ˆ cijkl 2k 2l and for an isotropic one

(2.2)

Let us formulate the conditions that must be satis®ed on the crack surfaces. We designate

V1 ˆ

N [ nˆ1

…2i uj 1 2j ui †

Aij uj 1 bi ˆ r22t ui ; ;x [ V; ;t [ I

conditions. We present these conditions in the form

V n1 ;

V2 ˆ

N [ nˆ1

V n2

The mutual displacements of the crack surfaces are characterized by the displacement discontinuity vector [4,27] Du…x; t† ˆ u1 …x; t† 2 u2 …x; t† ;x [ V 1 < V 2 ; ;t [ I The contact forces of the interaction at the crack edges are connected with the components of the stress tensor by the relations qi …x; t† ˆ 2s ij …x; t†nj …x†; ;x [ V 1 < V 2 > V e ; ;t [ I This equality should be satis®ed for each of the crack surfaces 1 1 q1 i …x; t† ˆ 2s ij …x; t†nj …x†; ;x [ V > V e ; ;t [ I 2 1 q2 i …x; t† ˆ 2s ij …x; t†nj …x†; ;x [ V > V e ; ;t [ I S S S where V e ˆ n ˆ 1N V ne ˆ Nnˆ1 V n1 > Nnˆ1 V n2 ˆ V 1 > 2 V are the areas of complete contact at the crack edges, 2 1 n1 j …x† and nj …x† are the normal vectors to the surfaces V 2 and V . By virtue of the assumption regarding to the surface properties, the equality n j1 …x† ˆ 2n2 j …x† takes place. Therefore 2 qi …x; t† ˆ q1 i …x; t† ˆ 2qi …x; t†;

;x [ V 1 < V 2 > V e ; ;t [ I The contact between the crack surfaces is supposed to be unilateral. So, the normal component of the contact force vector cannot be tensile. The limitations on the tangential components of the contact force vector and the displacement discontinuity vector depend on the method used to take into account the friction at the contact area. It is usually assumed that friction occurs in accordance with Coulomb's law [16,27,28,33,34,38,42]. Taking into account the above explanations, the onesided restrictions with friction in the form of inequalities on the crack surfaces take the form [16,27,38] Dun $ h0 ; qn $ 0; …Dun 2 h0 †qn ˆ 0;

Aij ˆ mdij 2k 2k 1 …l 1 m†2i 2j :

;x [ V 1 < V 2 ; ;t [ I

For the correct formulation of the elastodynamic problems it is necessary to assign the initial and boundary

uqt u # kt qn ! 2t Dut ˆ 0; uqt u ˆ kt qn ! 2t Dut ˆ 2lt qt …2:3†

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

where qn ; qt ; Dun ; Dut are the normal and tangential components of the contact forces and the displacement discontinuity vectors, respectively, h0 is the initial opening of cracks and kt and l t are coef®cients which dependent on the properties of the contacting surfaces. Therefore, the elastodynamic contact problem for the body with cracks is reduced to initial-boundary problems (2.1) and (2.2), with restrictions in the form of inequality (2.3). The Laplace transform with respect to time is widely used for the solution of elastodynamic problems [4,15,18,22]. We have demonstrated that this method is quite effective in the solution of the elastodynamic unilateral contact problems for bodies with cracks [26±29,53,55]. We apply the Laplace transform Z1 f …x; k† ˆ L{f …x; t†} ˆ f …x; t†e2kt dt; ;k : Re…k† . Re…k0 † 0

…2:4† to Eq. (2.1) and to the initial and boundary conditions (2.2). Here k0 is chosen from the convergence condition of integral (2.4). Taking into account, that L{22t ui …x; t†} ˆ k2 ui …x; k† 2 ku0i …x† 2 v0i …x†; instead of having the initial-boundary problems (2.1) and (2.2) we obtain the boundary-value problem with parameter k in the form Aij uj …x; k† 1 Fi …x; k† 2 rk2 ui …x; k† ˆ 0;

…2:5†

pi …x; k† ˆ s ij …x; k†ni …x† ˆ c i …x; k†; ;x [ 2Vp ; ;k : Re…k† . Re…k0 † ui …x; k† ˆ wi …x; k†; ;x [ 2Vu ; ;k : Re…k† . Re…k0 † 2

for bodies with cracks. The ®rst one is based on the fact that the problem is solved at a real-time space [45,46]. The second approach uses the Laplace transformation with respect to time to solve the problem [25,26,53,55]. The problem with the unilateral restrictions on the body boundary using Laplace transformation is considered in Ref. [23]. 3. Boundary integral equations for elastodynamics The application of the boundary integral equation method in elastodynamics and fracture mechanics was studied in Refs. [4,5,7±9,13,15,27,36], where it was shown that in elastodynamics there are several formulations of the boundary integral equation method. To solve the elastodynamic contact problems with unilateral restrictions for bodies with cracks we use direct formulation of the boundary integral equation method in the Laplace transform space. In order to obtain integral representations of the displacement vector components in the Laplace transform space, it is suf®cient to apply the Laplace integral transform (2.4) to the Somigliana elastodynamic theorem [4,14,18,27], obtaining Z ui …x; k† ˆ pj …y; k†Uji …y 2 x; k†dS 2V

2 1

;x [ V; ;k : Re…k† . Re…k0 †

ku0i …x†

v0i …x†.

2 The where F…x; k† ˆ b…x; k† 1 k ui …x; k† 2 unilateral restriction (2.3) due to their non-linearity cannot be transformed. Using this approach, the initial-boundary problems (2.1) and (2.2) with the unilateral restriction (2.3) is replaced by an in®nite set of boundary-value problem (Eq. (2.5)) with the parameter k [ {k [ C : Re…k† . Re…k0 †} and the unilateral restriction (Eq. (2.3)). Therefore, the solution of problem (2.5) should be such that, after the application of the inverse Laplace transform 1 Zc 1 i;1 f …x; k†e2kt dk; …2:6† f …x; t† ˆ L21 {f …x; k†} 2pi c 2 i;1 to Dui …x; k† and qi …x; k† their physical analogies Dui …x; t† and qi …x; t†should satisfy not only the differential equation (2.1) with initial and boundary conditions (2.2) but also, the unilateral restriction (2.3). There are at least two approaches for the solution of elastodynamic contact problems with unilateral restrictions

647

Z 2V

Z V

uj …y; k†Wji …y; x; k†dS

fj …y; k†Uji …y 2 x; k†dV

…3:1†

The vector of the internal forces may be represented in the form Z pj …y; k†Kji …y 2 x; k†dS pi …x; k† ˆ 2V

2 1

Z 2V

Z V

uj …y; k†Fji …y; x; k†dS

fj …y; k†Kji …y 2 x; k†dV

…3:2†

The kernels of the integral operators (3.1) and (3.2) are the elastodynamic Green's functions in the Laplace transform space [4,14,18,27]. We will introduce the following dynamic potentials of the elastodynamic theory in the Laplace transform space in order to compact the boundary integral equations and to facilitate their analysis Z fj …y; k†Uji …y 2 x; k†dV Uik …f; x; V† ˆ V

Uik …p; x; 2V† ˆ Wik …u; x; 2V† ˆ

Z 2V

Z 2V

pj …y; k†Uji …y 2 x; k†dS uj …y; k†Wji …y; x; k†dS

648

Kik …f; x; V† ˆ

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

Z

Kik …p; x; 2V† ˆ Fik …u; x; 2V† ˆ

form

fj …y; k†Kji …y; x; k†dV

V

Z 2V

Z 2V

Wik …y; x; k† ˆ rb…c21 2 2c22 †nk 2r Uir 1 c22 …2j Uik 1 2k Uij †nj c

pj …y; k†Kji …y; x; k†dS uj …y; k†Fji …y; x; k†dS

Kik …y; x; k† ˆ r‰…c21 2 2c22 †ni 2r Ukr 1 c22 …2j Uki 1 2i Ukj †nj Š (3.3)

Fik …y; x; k† ˆ r‰…c21 2 2c22 †ni 2r Ukr 1 c22 …2j Uki 1 2i Ukj †nj Š …3:6†

Then formulas (3.1) and (3.2) will have the form ui …x; k† ˆ Uik …p; x; 2V† 2 Wik …u; x; 2V† 1 Uik …f; x; V† pi …x; k† ˆ Kik …p; x; 2V† 2 Fik …u; x; 2V† 1 Kik …f; x; V†

…3:4†

These are the initial formulas used to construct different kinds of the boundary integral equation for elastodynamics and, in particular, for bodies with cracks and slits [4,13,27,55]. In order to apply the boundary integral equation method for the solution of the problems under consideration, it is necessary to know the fundamental solutions of the elastodynamic theory and their main properties. The fundamental Green's tensor for displacements represents the solution of the elastodynamic differential equations in the Laplace transformed space [4,14,18,27] ‰Aij …2x † 2 k2 dij ŠUjk …y 2 x; k† ˆ dij d…y 2 x† where Aij …2x † ˆ …c21 2 c22 †2i 2j 1 c22 dij 2k 2k It has the form Uij ˆ …aprc22 †21 …cd ij 2 x2i r2j r†; r 2 ˆ …yi 2 xi †…yi 2 xi † …3:5† For a three-dimensional case a ˆ 4 and 21 2l2 21 2l1 2 x ˆ ‰…3l22 2 c22 …3l22 =c1 Š=r 2 1 3l2 1 1†e 1 1 3l1 1 1†e

21 2l2 21 2l1 2 c ˆ ‰e 2l 2 1 …l22 2 c22 …l22 =c1 Š=r 2 1 l2 †e 1 1 l1 †e

for a two-dimensional case a ˆ 2 and

x ˆ K2 …l2 † 2 c22 =c21 K2 …l1 † 2 2 c ˆ K 0 …l2 † 1 l21 2 ‰K1 …l2 † 2 c2 =c1 K1 …l1 †Š

Here, l1 ˆ kr=c1 and l2 ˆ kr=c2 . c1 and c2 are velocities of dilatational and distortional waves, respectively. In the case of steady-state oscillations, it is necessary to replace in these equations k by iv . Also, when n ˆ 2 the McDonald's functions Kn …li † must be replaced by Hankel's functions Hn …li † [43]. The remaining kernels in Eq. (3.4) are obtained by applying the differential operator Pik ˆ cijlk nj 2l ˆ lni 2k 1 m…dik 2n 1 nk 2i † to Uji …y 2 x; k†. The corresponding expressions have the

As it is known [4,36] when y ! x in the three-dimensional case Uji …y; x; k† ! r 21 ; Wji …y; x; k† ! r22 ; Kji …y; x; k† ! r 22 ; Fji …y; x; k† ! r 23 and in the two-dimensional case Uji …y; x; k† ! `n…r†; Wji …y; x; k† ! r 21 ; Kji …y; x; k† ! r 21 ; Fji …y; x; k† ! r 22 The analysis of these formulas shows that the boundary potentials Uik …p; x; 2V† contain kernels with weak singularities. This means, that they are continuous in R 3 and therefore, they can be continuously extended on the boundary 2V. The potentials Wik …p; x; 2V† and Kik …p; x; 2V† contain kernels with singularities. They become discontinuous as they cross the boundary. The potentials Fik …p; x; 2V† contain kernels with strong singularities. Notice that these potentials continuously intersect the boundary. The boundary properties of potentials (3.5) and (3.6) have been studied in Refs. [4,27,30,36]. Therefore, we will only present the ®nal results here. They are expressed by the equalities Uik …p; x; 2V†^ ˆ Uik …p; x; 2V†0 ; Wik …p; x; 2V†^ ˆ 7 12 ui …x; k† 1 Wik …u; x; 2V†0 Kik …p; x; 2V†^

1 2

ˆ ^ pi …x; k† 1

Kik …p; x; 2V†0 ;

…3:7†

Fik …u; x; 2V†^ ˆ Fik …u; x; 2V†0 The symbols ª ^ º and ª 7 º denote, that we have two equalities: one with the top sign and another with the bottom sign. The superscript ª0º denotes that the direct value of the corresponding potentials should be taken on the surface 2V. The properties of the boundary potentials on the crack surface are determined by the fact that the distance between the surfaces V 1 and V 2 is small compared to the linear dimensions of the crack. According to Refs. [27,31,54] we can obtain the following boundary potential properties on

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

the crack surface

649

T^ ik ˆ l2 ni …y†2l Ulr D^ rk 1 lm‰nj …y†2l Uli D^ jk 1 nr …y†2l Ulr D^ ik

U ik …p; x; V 1 < V 2 † ˆ 0;

1 ni …y†2l Ukr D^ rl 1 ni …y†2k Ulr D^ rl Š 1 m‰nj …y†2l Uki D^ ji

Wik …u; x; V 1 < V 2 † ˆ Wik …Du; x; V†

1 nj …y†2k Uli D^ ji 1 nr …y†2l Ukr D^ il 1 nr …y†2k Ulr D^ il ŠTik

K ik …p; x; V 1

2

…3:8†

Fik …u; x; V 1 < V 2 † ˆ Fik …Du; x; V† Using formula (3.4) we can obtain the integral representations of the surface forces and displacements vectors on the surface of the body 2V. Also, we can obtain the properties of the surface potentials on the boundary (Eq. (3.7)) and on the crack surface (Eq. (3.8)). On the smooth parts of the boundary, the integral representations (Eq. (3.4)) have the following form 1 2

u i …x; k† ˆ Uik …p; x; 2V† 2 Wik …u; x; 2V† 2 Wik …Du; x; V† 1 U ik …f; x; V†; ;x [ 2V

1 2

(3.9)

4. Regularizations of the hypersingular integrals The boundary potentials F ik …p; x; 2V† are integral operators, their kernels contain a strong singularity. The integrals with such kernels do not exist as Riemann's integrals or as Lebesgue's integrals, or as principal value according to Cauchy. There are at least two approaches that can be applied to study this problem. In the ®rst one, the integral operators with hypersingular kernels are replaced by integrodifferetial operators with the kernels that are the superposition of differential operators and singular kernels. Following Refs. [4,27], instead of the boundary potentials Fik …p; x; 2V†, we will get two potentials with the following form Tik …p; x; 2V† ˆ rk2 T^ ki …p; x; 2V†

ˆ rk

2

Z 2V

Z 2V

uj …y; k†Tji …y; x; k†dS …4:1† uj …y; k†T^ ji …y; x; k†dS

…4:2†

T ik …u; x; V 1 < V 2 † ˆ Tik …Du; x; V†

pi …x; k† ˆ Kik …p; x; 2V† 2 Fik …u; x; 2V† 2 Fik …Du; x; V† ;x [ 2V

T^ ki …u; x; 2V†^ ˆ T^ ki …u; x; 2V†0 ;

T^ ki …p; x; V 1 < V 2 † ˆ T^ ki …Du; x; V†;

1 K ik …f; x; V†; ;x [ 2V

K ik …f; x; V†;

Here D^ ik ˆ ni …x†2k 2 nk …x†2i are differential operators. The kernels Tji …y; x; k† are expressed through Uji …y; x; k† and, therefore, the potentials Tik …y; x; 2V† contain a weak singularity. The kernels T^ ji …y; x; k† are expressed through the ®rst derivative of Uji …y; x; k†. The resulting potentials T^ ki …y; x; k† are the singular integrodifferetial operators. The boundary properties and the properties on the crack surfaces of potentials (4.1) are de®ned by analogy with Eqs. (3.7) and (3.8), and they are expressed by the following formulas

Tik …u; x; 2V†^ ˆ Tik …u; x; 2V†0

pi …x; k† ˆ Kik …p; x; 2V† 2 Fik …u; x; 2V† 2 Fik …Du; x; V†

1

ˆ lni …y†nr …x†Ukr 1 mnj …y†‰nj …x†Uki 1 ni …x†Ukj Š

< V † ˆ 0;

Here we use the notation from Eqs. (3.7) and (3.8). Considering Eqs. (3.7), (3.8) and (4.2), the integral representations of the surface traction and the displacement vectors at the smooth parts of the boundary 2V and crack surfaces V have the form 1 2

u i …x; k† ˆ Uik …p; x; 2V† 2 Wik …u; x; 2V† 2 Wik …Du; x; V† 1 U ik …f; x; V†

1 2

pi …x; k† ˆ Kik …p; x; 2V† 2 T^ ki …u; x; 2V† 2 Tik …u; x; 2V† 2 T^ ki …Du; x; V† 2 T ik …Du; x; V† 1 K ik …f; x; V†

pi …x; k† ˆ Kik …p; x; 2V† 2 T^ ki …u; x; 2V† 2 Tik …u; x; 2V† 2 T^ ki …Du; x; V† 2 T ik …Du; x; V† 1 K ik …f; x; V† …4:3† The second approach to study the potentials Fik …p; x; 2V†; which involve the Hadamard's ®nite-part integrals [31,44]. We have shown in Refs. [27,48,57±59] how the Hadamard's theory is applied to the elastodynamic problems for the bodies with cracks. In order to regularize the boundary potentials Fik …u; x; 2V† we transform

650

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

if we consider the corresponding integrals in the sense of the Hadamard's ®nite part and the Cauchy's main value. We will consider a variant of the boundary integral equations method which is based on formula (3.9). Taking into account the boundary and the crack edges conditions we obtain [27,28]

them to the form [27,48] Fik …u; x; 2V† ˆ Fik …u; x; 2V=2Ve † 1 1

Z 2Ve

Z 2Ve

uj …y; k†‰Fji …y; x; k† 2 Fji …y; x†ŠdS ‰uj …y; k† 2 uj …x; k† 2 2t uj …x; k†…y 2 x†ŠFji …y; x†dS

1 uj …x; k†

Z 2Ve

1 2t uj …x; k†

ˆ

Fji …y; x†dS

Z 2Ve

1 2

…y 2 x†Fji …y; x†dS

(4.4)

Only the last two terms contain singularities. Let us introduce the Cartesian coordinate system on 2Ve such, that the origin is located at the y point and the x3 axis coincides with an external normal to 2Ve ; while the other two axes lie at the tangential plane P e . We will transform the integrals with singularities to the form Z

Fji …y; x†dS ˆ

2Ve

1 Z

Z Pe

Z 2Ve

{Fji …y; x† 2 Fji ‰p…y†; xŠ}dS

ˆ

Z 2Ve

1

Fji ‰p…y†; xŠdS

Pe

‰p…y† 2 xŠFji ‰p…y†; xŠdS

where p : 2Ve ! P e is the operator for orthogonal projection from 2Ve to P e : If the 2Ve is smooth then the ®rst two integrals on the right side of these equalities are regular ones. The potentials Z Pe

‰p…y† 2 xŠFji ‰p…y†; xŠdS and

c i …x; k† 2 f…x; k†; ;x [ 2Vp

pi …x; k† 2 Kik …p; x; 2Vu † 1 Fik …u; x; 2Vp † 1 Fik …Du; x; V† ˆ f…x; k†; ;x [ 2V u

pi …x; k† 2 Kik …p; x; 2Vu † 1 Fik …u; x; 2Vp † 1 Fik …Du; x; V† ˆ f…x; k†; ;x [ V

(4.6)

where f…x; k† ˆ K ik …f; x; V† 1 Kik …c; x; 2V p † 2Fik …w; x; 2Vu †. These equations are valid ;k : Re…k† . Re…k0 †. In the case of an unbounded body, system (4.6) is transformed to the form …4:7†

5. Mathematical analysis and solution algorithm

{…y 2 x†Fji …y; x† 2 ‰p…y† 2 xŠFji ‰p…y†; xŠ}dS

Z

1 2

pi …x; k† ˆ 2Fik …Du; x; V†; ;x [ V

…y 2 x†Fji …y; x†dS

2Ve

Kik …p; x; 2Vu † 2 Fik …u; x; 2Vp † 2 Fik …Du; x; V†

Z Pe

Fji ‰p…y†; xŠdS

…xi 2 yi †…xj 2 yj †=r 5

which must be calculated over the plane domains. Such integrals were calculated for circular, triangular and rectangular domains in Refs. [27,48,57±59]. For the two-dimensional case, the problem is reduced to the calculation of the integrals r 21 and r 22 over a linear segment. It is evident that Fik …u; x; 2V† ˆ T^ ki …u; x; 2V† 1 Tik …u; x; 2V†;

r L : Hl;m g …V £ I† ! Hg …V; k†;

L21 : Hrg …V; k† ! Hl;m g …V £ I†

contain kernels that are fundamental solutions of the static theory of elasticity. As the result of this analysis, they are transformed into the integrals r 22 ; …xi 2 yi †…xj 2 yj †=r 4 ; r 23 and

The mathematical investigations of various physical problems with one-sided restrictions were described in Refs. [2,16,17,19,20,27,32,37,38,53,54,56]. It was shown in those publications that such problems do not have a classical solution, requiring special functional spaces for their mathematical investigations [3,12,16,37,38,40,41]. We will use the following functional spaces H gl;m …V £ I† and Hrg …V; k† since in the case of an arbitrary dynamic loading the approach with Laplace transformation is applied. The direct and inverse Laplace transforms act on the pair of these functional spaces in the following way

…4:5†

…5:1†

In the case of harmonic loading, the calculation of the Fourier coef®cients corresponds to the direct Laplace transformation and the summation of the Fourier series corresponds to the inverse Laplace transformation v ZT u …x; t†eivk t dt F :! uki …x† ˆ 2p 0 i …5:2† ( ) 1 X 21 k i vk t ui …x†e F :! ui …x; t† ˆ Re 21

In this case, it should be assumed that g ˆ 2iv in Eq. (5.1). We will consider the variational (ªweakº) formulation of the problem. Therefore, the initial data and the solutions

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

the functional [27,54]

belong to the functional spaces

wi [

H1=2;1=2 …V g

£ I†; c i [

H21=2;21=2 …V g

£ I†;

infsup‰C…u i ; pi †Š ˆ inf{c…u i † 1 sup‰kpn ; F n …Dun †l

bi [ H1;1 g …V £ I†

1 kpt ; F t …Dut †lŠ}

0;0 ui [ H1;1 g …V £ I†; s ij ; eij [ Hg …V £ I†

pi [ bi …ui † or pi [ 2ji …ui †:

…5:3†

These relationships may be transformed to boundary conditions for the stress tensor components ui [ bci …pi † or ui [ 2jci …pi †;

…5:4†

2jci …pi †

ˆ is a maximal monotonic operator, where conjugated to bi …pi †. Now, we will consider the variational formulation of the elastodynamic contact problems with unilateral restrictions for the bodies with cracks. For this purpose, we will use the subdifferential boundary functionals (4.3) and (4.4) and the Hamilton±Ostrogradski principle. Following Refs. [27,28,53] we will obtain the variational inequality in the form Z Z …pi 2 c i †…vi 2 ui †dS dt I

2 1

2Vp

Z Z I

Z I

2Vu

…ui 2 wi †Pik …vk 2 uk †dS dt

F i …Dvi † 2 F i …Dui †dt $ 0

…5:5†

8Z > < k t upn uuut udV; if u i $ h0 ; upt u ˆ kt upn V F i …ui † ˆ > : 1; otherwise where ui and vi [ K…ui † …2V £ I†; upt u , kt upn u ) 2t Dut ˆ 0; Ku …ui † ˆ {ui [ H1=2;0 g upt u ˆ kt upn u ) 2t Dut ˆ lt pt }

…5:7†

ui [ Ku …ui † pi [ Kp …pi †

These functional spaces and their properties were studied in Refs. [3,12,16,40,41]. The theory of variational inequalities and the convex analysis are used for the mathematical investigation of problems with boundary conditions in the form of inequalities. Information about this subject may be found in Refs. [10,17,38]. We will consider a maximal monotonic operator bi : ui ! pi and a superpotential ji such that bi ˆ 2ji [38]. Here 2 designates the subdifferential of ji. This functional represents the energy of local connection. The boundary condition in the displacements may be written in the form

bci …pi †

651

(5.6)

For the solution of the variational inequality (5.5) we will use methods from the duality theory [10,17]. The variational inequality (5.5) will be transformed to the problem of ®nding infsup‰C…u i ; pi †Š; i.e. a saddle point of

…2V £ I†; pn $ 0; Kp …pi † ˆ {pi : pi [ H21=2;0 g upt u # kt pn ; ;x [ V; ;t [ I}

…5:8†

where ( supkp n ; F n …Dun †l ˆ

0;

if Dun 2 h0 $ 0;

1;

if Dun 2 h0 , 0

pn [ Kp …pi †

8Z > < k t pn uDut udV; if up t u # kt pn ; V supkpt ; F t …Dut †l ˆ > : 1; if upt u $ kt pn pn [ Kp …pi †

(5.9)

The problem of ®nding the saddle point for the subdifferentional functional (5.7) on sets (5.6) and (5.8) is reduced to the successive determination of the inf‰C…u i ; pi †Š for preassigned values of the subdifferentional functional (5.9). At the next step of the iteration, using the known values for pi and Dui the new values of these functionals will be calculated more precisely. We will use the Udzavy-type algorithm [10,17,21,32,37] developed in Refs. [45,46,53]. The algorithm is summarized as follows. 1. The initial distribution of the contact forces on the crack edges p0i …x; t†; ;x [ V; ;t [ I is assigned. 2. The Laplace transformation of the contact forces vectors are calculated using the equation p i0 …x; k† ˆ L{p0i …x; t†}; ;x [ V; ;k : Re…k† . Re…k 0 †: 3. The system of the boundary integral equations (4.6) is solved and the unknown functions on the boundary pi …x; k†; ui …x; k† and Dui …x; k†; ;k : Re…k† . Re…k0 † are de®ned. 4. The components of the displacement discontinuity vector on V are calculated from Du i …x; k† using inverse Laplace transformation Du1i …x; t† ˆ L21 {Du1i …x; k†}; ;x [ V; ;k : Re…k† . Re…k 0 †: 5. The normal and tangential components of the contact

652

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

forces are corrected in such way to satisfy restrictions (2.3) p1n …x; t† ˆ Pn {p0n …x; t† 2 rn bDu1n …x; t† 2 h0 …x†c} p1t …x; t†

ˆ

Pt {p0t …x; t†

2

rt 2t Du1t …x; t†}

…5:10†

where ( Pn …pn † ˆ

ˆ

0;

if pn # 0

pn ;

if pn . 0

8 > < pt ; > : k t pn

pt ; upt u

Pt …pt †

if upt u # kt pn if upt u $ kt pn

are operators for orthogonal projection into sets of the restrictions pn $ 0 and upt u # kt pn . The coef®cients rn and rt are chosen based on the necessary conditions for the best algorithm convergence. 6. Proceed to the second step.

The operators of the discrete Laplace transforms (direct and inverse) may be written in the form r Lm ˆ Pk +L+Pt : Hr;s g …2V £ Im † ! Hg …V; kl †

…Lm †21 ˆ Pt +L21 +Pk : Hrg …V; kl † ! Hr;s g …2V £ Im †

The numerical calculation of the direct and inverse Laplace transforms was studied in Refs. [4,6,27,35]. For harmonic loading the direct and inverse Laplace transformations correspond to the calculation of the Fourier coef®cients and to the summation of the Fourier series according to formula (5.2). This problem is reduced to the numerical integration and to the calculation of the ®nite sum of the Fourier series. Now we will consider the approximation of the functional 21=2 …2V; k†. We divide the boundspaces H1=2 g …2V; k†and Hg ary of the body into N boundary elements 2V ˆ

N [ nˆ1

The mathematical investigation of this algorithm convergence was studied in Refs. [27,46,54,57]. Other algorithms and their mathematical investigations were studied in Refs. [1,2,10,17,21,32±34,37,42]. 6. The discrete boundary integral equations The analytical solution of a dynamic contact problem with unilateral restrictions for bodies with cracks is associated with severe dif®culties. Therefore, we will solve it using numerical methods. For this purpose, it is necessary to construct the discrete subspaces of the Sobolev's spaces r Hr;s g …2V £ I† and Hg …V; k†; the operators of a projection which act on these functional spaces. It is also necessary to construct the discrete operators for the direct and inverse Laplace transforms and the integral operators in the discrete form. These discrete functional spaces and operators are used in the numerical solution of the problem. First, we will consider the problem of time and space of the Laplace transform discretization. The direct and inverse Laplace transforms act on the Sobolev's spaces Hr;s g …2V £ I† and Hrg …V; k† according to formula (5.1). We will construct a discrete Sobolev space Hr;s g …2V £ Im † in which a time variable t takes discrete values from a set Im ˆ {t0 ; t1 ; ¼; tm }: The discrete Sobolev's space Hrg …V; kl † is constructed in the same way, in this case parameter k of the Laplace transform imparts L discrete values which must satisfy the condition Re…k1 † . Re…k0 †: The projection operators on these spaces are de®ned in the following way r;s Pt : Hr;s g …2V £ R† ! Hg …2V £ Rm †;

Pk : Hrg …V; k† ! Hrg …V; kl †

…6:1†

2Vn ; 2Vn

\

2Vk ˆ B; if n ± k

We consider Q nodes on each of the boundary elements and construct the local projection operators of the spaces 21=2 …2Vn ; k† on the ®nite-dimensional H1=2 g …2Vn ; k† and Hg g subspaces Xq …2Vn ; k† and Ygq …2Vn ; k† g Puq : H1=2 g …2Vn ; k† ! Xq …2Vn ; k†; ;x [ 2Vn

Ppq : H21=2 …2Vn ; k† ! Ygq …2Vn ; k†; ;x [ 2Vn g We de®ne the global projection operators as the sum of the local operators Punq ˆ

N X nˆ1

Puq ; Ppnq ˆ

N X nˆ1

Ppq

then Punq Ppnq

:

:

H1=2 g …2V; k†

!

H21=2 …2V; k† g

Xgq

!

Ygq

N [ nˆ1

! 2Vn ; k ; ;x [ 2V

N [ nˆ1

! 2Vn ; k ; ;x [ 2V

The local projection operators Pun and Ppn map the displacement and surface traction vectors de®ned on the boundary elements 2Vn into sets of their value at the nodes of interpolation Puq ‰ui …x; k†Š ˆ {uni …xq ; k†; q ˆ 1; ¼; Q}; ;x [ 2Vn Ppq ‰pi …x; k†Š ˆ {pni …xq ; k†; q ˆ 1; ¼; Q}; ;x [ 2Vn

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

in the same way for the global operators

Kjin …xr ; yq ; k† ˆ

Punq ‰ui …x; k†Š ˆ {uni …xq ; k†; q ˆ 1; ¼; Q; n ˆ 1; ¼; N};

Fjin …xr ; yq ; k† ˆ

;x [ 2V Ppnq ‰pi …x; k†Š ˆ {pni …xq ; k†; q ˆ 1; ¼; Q; n ˆ 1; ¼; N}; ;x [ 2V We will construct the projection operators …Pun †21 and For this purpose we will introduce the systems of basic functions wnq …x† and c nq …x† in the subspaces Xgq …2Vn ; k† and Ygq …2Vn ; k†: Then the displacement, displacement discontinuity and surface traction vectors are presented in the form

…Ppn †21 :

ui …x; k† <

N X

nˆ1 qˆ1

Dui …x; k† <

pi …x; k† <

Q X

N X

uni …xq ; k†wnq …x†;

Q X

nˆ1 qˆ1

Q N X X nˆ1 qˆ1

wnq …x† [

Nu X Q X nˆ1 qˆ1

2

c nq …x† [ Ygq …2Vn ; k†

1 2

1 2

um i …xr ; k† ˆ

N X

Q X

nˆ1 qˆ1

1

pm i …xr ; k† ˆ

Q N X X nˆ1 qˆ1

‰Kjin …xr ; yq ; k†pnj …yq ; k†

2 Fjin …xr ; yq ; k†unj …yq ; k† 2 Fjin …xr ; yq ; k†Dunj …yq ; k†Š 1 Fik …f; x; Vn †

Np Q X X

1

Ne X Q X nˆ1 qˆ1

Ujin …xr ; yq ; k† ˆ Wjin …xr ; yq ; k† ˆ

2Vn

Z 2Vn

Uji …xr ; yq ; k†c nq …y†dS; Wji …xr ; yq ; k†wnq …y†dS;

Np Q X X nˆ1 qˆ1

Fjin …xr ; yq ; k†unj …yq ; k†

Fjin …xr ; yq ; k†Dunj …xq ; k†

Nu X Q X nˆ1 qˆ1

Kjin …xr ; yq ; k†pnj …yq ; k†

Fjin …xr ; yq ; k†unj …yq ; k† Fjin …xr ; yq ; k†Dunj …xq ; k†

Nu X Q X nˆ1 qˆ1

Np Q X X

1

Ne X Q X nˆ1 qˆ1

Kjin …xr ; yq ; k†pnj …yq ; k†

Fjin …xr ; yq ; k†unj …yq ; k† Fjin …xr ; yq ; k†Dunj …xq ; k†

ˆ fni …xr ; k†; ;xr [ V

(6.4)

where

f in …xr ; k† ˆ Kik …f; xr ; k† 1

where Z

Fji …xr ; yq ; k†wnq …y†dS;

ˆ fni …xr ; k†; ;xr [ 2Vu

nˆ1 qˆ1

…6:3†

2Vn

c in …xr ; k† 2 fni …xr ; k†; ;xr [ 2Vp

nˆ1 qˆ1

1

2 Wjin …xr ; yq ; k†unj …yq ; k†

1 2

1 2

pni …xr ; k† 2

‰Ujin …xr ; yq ; k†pnj …yq ; k†

2 Wjin …xr ; yq ; k†Dunj …yq ; k†Š 1 Ujk …f; x; Vn †

Ne X Q X

pni …xr ; k† 2

…6:2† Here and below, if the node belongs to several boundary elements, it will be accounted once in these sums. Inserting expression (6.2) into (3.9) we will obtain the ®nite-dimensional representations for the displacement and surface force vectors

Z

Kji …xr ; yq ; k†c nq …y†dS;

Kjin …xr ; yq ; k†pnj …yq ; k† 2

nˆ1 qˆ1

ˆ

pni …xq ; k†c nq …x†;

2Vn

The system of linear algebraic equations for the boundary element method corresponds to the elastodynamic boundary integral equations for the bodies with cracks in the Laplace transform space. Now considering Eq. (4.6), the boundary integral equations have the form

Xgq …2Vn ; k†

Duni …xq k†wnq …x†; wnq …x† [ Xgq …V n ; k†

Z

653

2

Nu X Q X nˆ1 qˆ1

Np Q X X nˆ1 qˆ1

Kjin …xr ; yq ; k†c jn …yq ; k†

Fjin …xr ; yq ; k†wnj …yq ; k†

The boundary integral equations in the case of an unbounded body (4.7) are transformed into a system of

654

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

Here D ˆ 21 21 1 22 22 is the two-dimensional Laplace operator. The symbol ª p º means that the stresses and displacements are caused by the incident wave. We will consider the problem for re¯ected waves. The load on the crack edges caused by the incident wave has the form ( pn …x; t† ˆ

p0n

pt …x; t† ˆ

p0t

e (

Fig. 2. Plane with two ®nite length collinear cracks interacting with a harmonic tension-compression wave.

ˆ

Ne X Q X nˆ1 qˆ1

Fjin …xr ; yq ; k†Dunj …xq ; k†

ei…k1 x1 e

if ;x [ V 1

i{k1 ‰…x1 1h1 †cos a1h2 sin aŠ2vt} cos a2vt†

;

if ;x [ V 2 if ;x [ V 1

;

i{k1 ‰…x1 1h1 †cos a1h2 sin aŠ2vt}

; if ;x [ V 2

pn ˆ 2mk22 ‰1 2 2…c2 =c1 †2 cos2 aŠf0 ;

linear algebraic equations pni …xr ; k†

ei…k1 x1 cos a2vt† ;

…6:5†

pt ˆ 2mk22 …c2 =c1 †2 sin2 af0 ; p k2 ˆ v=c2 ; c2 ˆ m=r

In Eqs. (6.4) and (6.5), Nu, Np and Ne represent the number of the boundary elements on 2Vu, 2Vp and V , respectively. The equations for the boundary element method, their solution and application to science and engineering were studied in detail in Refs. [4,5,8,9,13,15].

The contact forces vector q…qn ; qt † and displacement discontinuity vector Du…Dun ; Dut † on the cracks edges must satisfy the unilateral contact conditions with friction:

7. The harmonic loading in a plane with cracks

;x [ V; ;t [ ‰0; TŠ

Let us consider two collinear cracks of length l1 and l2 in a plane R 2 ˆ {x : x3 ˆ 0} (Fig. 2). Their position is determined by the coordinates

uq t u # kt qn ! 2t Dut ˆ 0; uqt u # kt qn ! 2t Dut ˆ 2lt qt

V 1 ˆ {x : x2 ˆ 0 ; 2l1 # x1 # l1 } V 2 ˆ {x : x2 ˆ h2 ; h1 2 l1 # x1 # l1 1 h1 } where h1 and h2 are the distances between the crack centers. A harmonic tension-compression wave propagates in the plane. The incident wave is described by a potential function

c…x; t† ˆ c 0 ei…ki n´x2vt† p k1 ˆ v=c1 ; c1 ˆ …l 1 2m†=r; n ˆ …cos a; sin a† where l and m , are the Lame constants and r is the material density, v ˆ 2pT 21 is the vibration frequency, T is the vibration period, c 0 is the amplitude of the tensioncompression wave, c1 is the wave velocity and a is the angle of the incident wave. The potential function depends of two space coordinates x…x 1 ; x2 †: In this case, we deal with a two-dimensional problem. If the wave function c…x; t† is known, the complex amplitudes of the stress tensor and displacement vector are de®ned by the equations p p s 11 ˆ lDc p 1 2m21 21 c p ; s 22 ˆ lDc p 1 2m22 22 c p ; p s 12 ˆ 2m21 22 c p ; up1 ˆ 21 c p ; up12 ˆ 22 c p

Dun $ 2h0 ; qn $ 0; …Dun 1 h0 †qn ˆ 0;

…7:1† Here, we used the notation from Eq. (3.3). Due to the contact interaction, the load vector on the crack edges has the form p ˆ …p1 ; p2 †; p2 ˆ p0n 1 qn ; p1 ˆ p0t 1 qt ; ;x [ V e and q ˆ 0; ;x Ó V e where V e is a region of complete contact. The presence of the unilateral restriction (Eq. (7.1)) makes the problem non-linear. Therefore, the problem for re¯ected waves is described by periodic but not by harmonic functions. The components of the stress±strain state due to the re¯ected waves cannot be presented as functions of the coordinates multiplied by e 2ivt as it is usually done for the solution of elastodynamic problems with harmonic loading [18,22]. For the problem under consideration we will present the components of the displacement vector and stress tensor as Fourier series of the load parameter v ( ua …x; t† ˆ Re

1 X 21

(

s ab …x; t† ˆ Re

) uka …x†eivk t

1 X 21

; )

k s ab …x†eivk t

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

where vk ˆ vk; and uka …x† ˆ

integral equation (7.3). As it was mentioned earlier, these equations establish the relationship between the complex functions pka …x† and Duka …x†: The kernels of these integral equations have ªstrongº singularities. In addition, the integration area consists of two parts V ˆ V 1 < V 2 . Therefore the boundary integral equation (7.3) has the form

v ZT u …x; t†eivk t dt; 2p 0 a

k s ab …x† ˆ

v ZT s …x; t†eivk t dt 2p 0 ab

k …x† are expressed The Fourier coef®cients uka …x† and s ab through the Fourier coef®cients of displacement discontinuity Duka …x† in the form

uka …x† ˆ 2

Z V

k …x† ˆ 2 s ab

…7:2† V

21

(

1 X 21

) Duk …x†eivk t

v ZT p…x; t†eivk t dt; p …x† ˆ 2p 0 v ZT Du…x; t†eivk t dt 2p 0

The Fourier coef®cients for the coef®cients of the surface force vector pka …x† and the displacement discontinuity vector Duka …x† are connected by the following system of integral equations Z V

F ba …y; x; vk †Dukb …y†dV;

V1

Z V2

Z V2

p 2k …x† ˆ 2FP 2

Z

F11 …y; x; vk †Duk1 …y†dV F21 …y; x; vk †Duk2 …y†dV

Z V1

Z V2

Z V2

F11 …y; x; vk †Duk1 …y†dV

F22 …y; x; vk †Duk2 …y†dV

F22 …y; x; vk †Duk2 …y†dV F12 …y; x; vk †Duk1 …y†dV;

;x [ V 1 ; k ˆ 0; ^1; ^2; ¼; ^1:

k

pka …x† ˆ 2

2

2

where:

Duk …x† ˆ

2

F abg …y1 ; x; vk †Dukg …y1 †dV

In the same way, the surface load on the crack edges and its opening may be presented by their Fourier series (1 ) X k ivk t p …x†e ; p…x; t† ˆ Re

Du…x; t† ˆ Re

pk1 …x† ˆ 2FP

W ga …y1 ; x; vk †Dukg …y1 †dV

Z

655

…7:3†

k ˆ 0; ^1; ^2; ¼; ^1; The kernels W ba …y; x; vk †; Fba …y; x; vk † and F abg …y; x; vk † in Eqs. (7.2) and (7.3) are the elastodynamic fundamental solutions in the space of the Fourier series. The kernels Wba …y; x; vk † and Fba …y; x; vk † may be calculated using formula (3.6) and the Green's fundamental displacement tensor Uba …y; x; vk † from Eq. (3.5) for n ˆ 2 and k ˆ iv: To the calculate F abg …y; x; vk † it is necessary to apply the operator X ˆ ldab 2n 1 m…dan 2b 1 dbn 2a † abn

to Uyg …y; x; vk †: Let us consider in detail the structure of the boundary

For ;x [ V 2 the integral equations have the same structure with the substitution of V 1 by V 2 and vice versa. Notice, that the kernels F 12 …y; x; vk † and F21 …y; x; vk †; ;x; y [ V 1 and ;x; y [ V 2 are equal to zero. That is why, in the case of one rectilinear crack, the perpendicular to the crack edges load does not cause the tangential displacement discontinuity of the crack edges. Vice versa, the tangential to the crack edges load does not cause the perpendicular displacement discontinuity. For one curvilinear or for two and more collinear cracks these properties do not take place. The kernels of the integral equations Fab …y; x; vk † when ;x [ V 1 and ;y [ V 2 and vice versa, are smooth functions. Therefore, to calculate the coef®cients of the boundary element equations it can be done using the quadrature formula. When ;x; y [ V 1 or ;x; y [ V 2 the kernels are hypersingular and they must be considered in the sense of the Hadamard's ®nite part. The formulas for the kernels F11 …y; x; vk † and F22 …y; x; vk † after the separation of the real and imaginary parts, have the form Re …y; x; vk † F11

G ˆ p

"

y 2…1 2 y†

! vk Y01 …l1k † c2 !

! vk vk 1 k Y …l † 2 2 Y 1 …l k † rc1 1 1 rc2 1 2 # 6xRe …vk ; r; c1 ; c2 † 2 r2 1 2 2y 1 12y

656 Im F11 …y; x; vk †

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

G ˆ p

"

n 2…1 2 y†

! vk 1 k J …l † c2 0 1 !

! vk 1 k vk 1 k J …l † 2 2 J …l † rc1 1 1 rc2 1 2 # 6xIm …vk ; r; c1 ; c2 † 2 r2 " ! G y2 vk Re Y01 …l1k † F22 …y; x; vk † ˆ p 2…1 2 y 2 † c2 ! ! 2y vk vk 1 k Y …l † 1 2 Y 1 …l k † 1 1 2 y rc1 1 1 rc2 1 2 # 6xRe …vk ; r; c1 ; c2 † 1 r2 " ! G y2 vk 1 k Im J …l † F22 …y; x; vk † ˆ p 2…1 2 y 2 † c2 0 1 ! ! 2y vk 1 k vk 1 k J …l † 1 2 J …l † 1 1 2 y rc1 1 1 rc2 1 2 # 6xIm …vk ; r; c1 ; c2 † 1 r2 1 2 2y 1 12y

xRe ˆ 2Y21 …l2k † 1 c21 Y21 …l1k †=c22 ; xIm ˆ 2J21 …l2k † 1 c21 J21 …l1k †=c22 The analytical representations for these and other fundamental solutions in elastodynamics were presented in general form in Ref. [27]. The approach presented here was applied to some problems with harmonic loads on the crack edges, taking into account their contact interaction [24,25,27,28,47, 49,50,51,54]. We will consider some examples here. To simplify the computation and analysis of the results, we introduce the dimensionless coordinates x ˆ x1 =l and y ˆ x2 =l and the frequency wave number k1p ˆ vl=c1 : For the examples presented below the mechanical properties of the material are: elastic modulus E ˆ 200 GPa; Paisson ratio y ˆ 0:3 and speci®c density r ˆ 7800 kg=m 3 : The calculations were performed with k1p ˆ 0:45 and N ˆ 15 (number of boundary elements along the crack), Nt ˆ 36 (number of time intervals into which the loading period is broken up), Nf ˆ 10 (number of the Fourier terms). The dependency of the calculation accuracy of N, Nt and Nf was studied in [27]. First we will consider one crack and a harmonic tensioncompression wave, which propagates in the plane perpendicular to V as it is shown in Fig. 3. Following Refs. [27,49±52,54] we will consider the problem of re¯ected waves. Only this problem in¯uences the fracture mechanics criterions. Let a harmonic load

Fig. 3. Plane with a ®nite length crack interacting with a harmonic tensioncompression wave.

p 0 …x1 ; 0; t† ˆ cos…vt† be applied on the crack edges. The distribution of the contact forces and displacement discontinuity are presented in Fig. 4. These graphics were constructed by a computer using the results calculated with the above initial data. Many authors indicate [18,27,39] that the contact interaction of the crack edges affects the stress intensity factor (SIF). We will investigate this in¯uence for some special problems. Other examples may be found in Refs. [24± 29,49±52,54]. The stress±strain state in the neighborhood of the a crack tip is de®ned by the asymptotic formulas [11]   KI q q 3q 1 2 sin sin s 11 ˆ p cos 2 2 2 2pr   KII q q 3q 2 1 cos cos 1 p sin 2 2 2 2pr   KI q q 3q 1 1 sin sin cos s 22 ˆ p 2 2 2 2pr KII q q 3q cos sin cos 1 p 2 2 2 2pr KI q q 3q sin cos cos s 12 ˆ p 2 2 2 2pr   KII q q 3q 1 2 sin sin cos 1 p 2 2 2 2pr r   KI r q 2q 1 2 2y 1 sin cos u1 ˆ m 2p 2 2 r   K r q q 2 2 2y 1 cos 2 sin 1 II m 2p 2 2 r   K r q q 1 2 2y 1 cos 2 sin u2 ˆ I m 2p 2 2 r   K r q q 2y 2 1 1 sin 2 cos 1 II m 2p 2 2 were r and q are polar coordinates, connected with the crack tip. Usually these formulas are used to de®ne SIF. In [4,13] it is shown that more accurate results may be obtained using the equations for the displacements corresponding to the

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

657

p Fig. 6. KIm =p pl pro®les for h1 ˆ 0; k ˆ 0:45 : taking into account of a crack edges contact (solid) and without taking into account crack edges contact (dotted).

SIF exceeds the corresponding static values by 20 and 15% (with and without initial crack opening, respectively). p Some results of the dimensionless value KIm =p pl for two cracks are shown in Figs. 6 and 7. In this case, the maximum value of SIF with contact interaction and without the contact interaction of the crack edges may differ by 30% and more. Fig. 4. Contact forces and crack opening in a space and time for k ˆ 0:45:

same type of the fracture mode. Using these formulas we obtain p m 2p p Du 2 …l 2 r; t†; KI …t† ˆ lim r!0 4…1 2 y† r p m 2p p Du 1 …l 2 r; t† KII …t† ˆ lim r!0 4…1 2 y† r The problems related with the numerical calculation of the SIF were discussed in Refs. [4,13,27]. The results of the p dimensionless KIm =p pl calculation for the problem, which is shown in Fig. 3, are presented in Fig. 5. Here, KIm is a maximum value of KI …t†: Curve 1 corresponds to the problem without contact interaction. The curves 2 and 3 show the in¯uence of the contact interaction without the initial opening of the crack …h0 …x1 † ˆ 0† and with the initial opening …h0 …x1 † ˆ Du0 …1 2 ux1 u†; respectively, where Du0 is a maximum opening crack under the static load. For the problem with no contact interaction of the crack edges [18,39] these graphs show that the maximum value of the dynamic SIF exceeds the corresponding static value by 30%. For the problem with contact interaction the dynamic

p Fig. 5. KIm =p pl pro®les as functions of wave number: without taking into account crack edges contact (1). Taking into account crack edges contact: without initial opening (2) and with initial opening (3).

8. Conclusions The results discussed in this article lead to the following conclusions and recommendations for future progress: 1. The correct mathematical formulation of the elastodynamic problem for bodies with cracks must consider the possibility of contact interaction at the crack edges with the formation of close contact, adhesion and sliding areas. 2. An elastodynamic problem for arbitrary loading can be presented in the form of the boundary-value problem in a space of Laplace transforms and unilateral restrictions. 3. An elastodynamic problem for harmonic loading can be presented in the form of the boundary-value problems for Fourier coef®cients of the displacement vectors and unilateral restrictions. 4. The algorithm for the solution of the elastodynamic contact problems with unilateral restrictions consists of de®ning a saddle point for the boundary functional.

p Fig. 7. KIm =p pl pro®les for h2 ˆ 0; k ˆ 0:45 : taking into account of a crack edges contact (solid) and without taking into account crack edges contact (dotted).

658

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659

5. The singular and hypersingular integral equations may be effectively applied for the problem solution. The regularization methods for the divergent integrals used in elastodynamics have been developed. 6. It has been found, that contact interaction at the crack edges in¯uences the fracture mechanic criterions. The SIF may differ by 30% or more depending on whether we consider the contact interaction or not. The results presented here and in previous publications show how the contact interaction of the cracks edges in¯uences the fracture mechanics criterions. This must be taken into account when the strengths of structures are calculated using the fracture mechanics methods. References [1] Alibadi M, Brebbia CA. Contact mechanics, computational techniques. Southampton: Computational Mechanical Publications, 1993. [2] Antes H, Panagiotopoulos PD. The boundary integral approach to static and dynamic contact problems. Basel: Birkhauser, 1992. [3] Agranovitch MS, Vishik MI. Elliptic problems with a parameter and parabolic problems of the general form. Usp Mat Nauk 1964;19(3):53±161 (in Russian). [4] Balas J, Sladek J, Sladek V. Stress analysis by boundary element methods. Amsterdam: Elsevier, 1989. [5] Banerjee PK. Boundary element method in engineering science. New York: McGraw Hill, 1994. [6] Bellman RE, Kalaba RE, Lockett J. Numerical inversion of the Laplace transform. New York: Elsevier, 1966. [7] Beskos DE. Boundary element methods in dynamic analysis. Appl Mech Rev 1987;40:1±23. [8] Brebbia CA, Domingnez J. Boundary elements. An introductory course. New York: McGraw Hill, 1989. [9] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques. Theory and applications in engineering. Berlin: Springer, 1984. [10] Cea J. Optimization. Teorie et algorithmes. Paris: Dunod, 1971 (in French). [11] Cherepanov GP. Mechanics of brittle fracture. New York: McGraw Hill, 1979. [12] Chudirovich IYu. Method of boundary equations in dynamic problems of the elasticity theory, Kharkov, Deponted in VINITI 26.06.90. No 3649-B90, 1990, 121pp. (in Russian). [13] Cruse TA. Boundary element analysis in computational fracture mechanics. Dordrecht: Kluwer, 1988. [14] Cruse TA, Rizzo FJ. A direct formulation and numerical solution of the general transient elastodynamic problem. 1. J Math Anal Appl 1968;22(1):244±59. [15] Dominguez J. Boundary elements in dynamics. Southampton: Computational Mechanical Publications, 1993. [16] Duvaut G, Lions J-L. Les inequations en mecanigue et en physique. Paris: Dunod, 1972. [17] Ekeland I, Temam R. Convex analysis and variational problems. Amsterdam: North-Holland, 1975. [18] Eringen AC, Suhubi ES, Elastodynamics. Linear theory. New York: Academic Press, 1975. p. 343±1003. [19] Fichera G. Existence theorems in elasticity. In: Flugge S, editor. Encyclopedia of physics, vol. V1u/2. Berlin: Springer, 1972. [20] Fichera G. Boundary value problems of elasticity with unilateral constraints. In: Flugge S, editor. Encyclopedia of physics, vol. V1u/ 2. Berlin: Springer, 1972. [21] Glowinski R, Lions J-L, Tremolieres R. Numerical analysis of variational inequalities. Amsterdam: North-Holland, 1981.

[22] Guz' AN, Kubenko VD, Cherevko MA. Diffraction of elastic waves. Kiev: Naukova Dumka, 1978 (in Russian). [23] Guz' AN, Zozulya VV. Dynamic problems of the elasticity theory with restrictions in the form of ineqealities. Dokl Akad Nauk Ukr SSR, Ser A Phys Math Tech Sci 1991;5:47±50 (in Russian). [24] Guz' AN, Zozulya VV. Dynamic problem for a plane with a slit. Taking into account of interaction of edges. Dokl Acad Nauk USSR 1991;318(2):304±7 (in Russian). [25] Guz' AN, Zozulya VV. Dynamic contact problem for a plane with two cracks. Dokl Acad Nauk USSR 1991;321(2):278±80 (in Russian). [26] Guz' AN, Zozulya VV. Contact interaction between crack edges under dynamic load. Int Appl Mech 1992;28(7):407±14. [27] Guz' AN, Zozulya VV. Brittle fracture of constructive materials under dynamic loading. Kiev: Naukova Dumka, 1993 (in Russian). [28] Guz' AN, Zozulya VV. Dynamic problems of fracture mechanic with account of the contact interaction of the crack edges. Int Appl Mech 1995;31(1):1±31. [29] Guz' AN, Zozulya VV. On taking into account of cracks edges contact under a dynamic load. Mater Sci 1996;32(1):38±52. [30] Gunter NM. Theory of potential and its application to main problems of mathematical physics. Moscow: GITTL, 1953 (in Russian). [31] Hadamard J. Le probleme de cauchy et les eguations aux devivees partielles lineaires hyperboliques. Paris: Herman, 1932. [32] Hlavacek I, Haslinger J, Necas J, Lovisec J. Solution of variational inequalities in mechanics. New York: Springer, 1988. [33] Kikuchi N, Oden JT. Contact problems in elasticity. Philadelphia, PA: SIAM, Publications, 1987. [34] Kravchuk AS. Theory of contact problems with allowance for friction on the surface of contact. J Appl Mat Mech 1980;44(1):122±9. [35] Krylov VI, Skoblya NS. Methods for an approximate Fourier transform and inversion of Laplace transform. Moscow: Nauka, 1974 (in Russian). [36] Kupradze VD, Gegelia TG, Bashelishvili MO, Burkhuladze TV. Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity. Amsterdam: North-Holland, 1979. [37] Necas J, Hlavacek I. Mathematical theory of elastic and elasticoplastic bodies: an introduction. Amsterdam: Elsevier, 1981. [38] Panagiotopoulos PD. Inequality problems in mechanics and applications. Convex and non convex energy functions. Stuttgart: Birkhauser, 1985. [39] Parton VZ, Boriskovsky BG. Dynamic mechanics of fracture. Moscow: Nauka, 1985 (in Russian). [40] Slobodetskii LN. Generalized Sobolev's spaces and their application to boundary problems for differential equations in partial derivatives. Am Math Soc Trans 1958;57(2):207±75. [41] Sobolev SL. Applications of functional analysis in mathematical physics. Trans Am Math Soc, Math Moh 1963:7. [42] Telega YuI. Variational methods in mechanical contact problems. Usp Mech 1987;10(2):3±95 (in Russian). [43] Watson GN. A treatise on the theory of Bessel functions. London: Cambridge University Press, 1966. [44] Wendland WL, Stephan EP. A hypersingular boundary integral method for two-dimensional screen and crack problems. Arch Rat Mech Anal 1990;112(4):363±90. [45] Zozulya VV. Dynamic problems of the crack theory with regions of contact, engagement and sliding. Dokl Acad Nauk Ukr SS Ser A Phys Math Tech Sci 1990;1:47±50 (in Russian). [46] Zozulya VV. Solvability of dynamic problems of the theory of cracks with regions of a contact, engagement and sliding. Dokl Acad Nauk Ukr SSR Ser A Phys Math Tech Sci 1990;3:53±5 (in Russian). [47] Zozulya VV. Action of harmonic load on a crack in an in®nite body with allowance for interaction of its edges. Dokl Acad Nauk Ukr SSR Ser A Phys Math Tech Sci 1990;4:46±9 (in Russian). [48] Zozulya VV. Integrals of Hadamard type in dynamic problem of the crack theory. Dokl Acad Nauk Ukr SSR Ser A Phys Math Tech Sci 1991;2:19±22 (in Russian).

V.V. Zozulya, P.I. Gonzalez-Chi / Engineering Analysis with Boundary Elements 24 (2000) 643±659 [49] Zozulya VV. Investigation of the contact of edges of cracks interacting with a plane, longitudinal, harmonic wave. Sov Appl Mech 1991;27(12):61±6. [50] Zozulya VV. Dynamic problem for a plane with two cracks. Taking into account the contact of edges. Dokl Akad Nauk Ukr SSR Ser A Phys Math Tech Sci 1991;8:75±80 (in Russian). [51] Zozulya VV. Contact interaction between the edges of a crack and an in®nite plane under a harmonic loading. Int Appl Mech 1992;28(1):61±5. [52] Zozulya VV. Investigation of the effect of crack edge contact for loading by a harmonic wave. Int Appl Mech 1992;28(2):95±100. [53] Zozulya VV. Method of boundary functionals in contact problems of elastodynamics for bodies with cracks. Dokl Akad Nauk Ukr 1992;2:38±44 (in Russian). [54] Zozulya VV. Harmonic loading of the edges of two collinear cracks in a plane. Int Appl Mech 1992;28(3):43±7.

659

[55] Zozulya VV. Solving problems of dynamics of bodies with cracks using the method of boundary integral equations. Dokl Akad Nauk Ukr 1992;3:38±43 (in Russian). [56] Zozulya VV. Mathematical investigation of the elastodynamic contact problems with friction for bodies with cracks, nonlinear analysis. Theory, methods and applications, 2000 (in preparation). [57] Zozulya VV, Lukin AN. Solution of three-dimensional problems of fracture mechanics by the method of integral boundary equations. Int Appl Mech 1998;34(6):544±51. [58] Zozulya VV, Gonzalez-Chi PI. Weakly singular, singular and hypersingular integrals in elasticity and fracture mechanics. J Chin Inst Eng 1999;22(6):763±75. [59] Zozulya VV, Menshicov WA. Hypersingular integrals in the tree dimensional elastodynamic problems for bodies with cracks. Int Appl Mech 2000;36(1):88±94.