Comparative study of time and frequency domain BEM approaches in frictional contact problem for antiplane crack under harmonic loading

Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Comparative study of time and frequency domain BEM approaches in frictional contact problem for antiplane crack under harmonic loading V.V. Zozulya n Centro de Investigación Científica de Yucatán A.C., Calle 43, No. 130, Colonia: Chuburná de Hidalgo, C.P. 97200, Mérida, Yucatán, México

art ic l e i nf o

a b s t r a c t

Article history: Received 16 April 2013 Accepted 5 August 2013 Available online 13 September 2013

Two different boundary element methods (BEM) for crack analysis in two dimensional (2-D) antiplane, homogeneous, isotropic and linear elastic solids by considering frictional contact of the crack edges are presented. Hypersingular boundary integral equations (BIE) in time-domain (TD) and frequency domain (FD), with corresponding elastodynamic fundamental solutions are applied for this purpose. For evaluation of the hypersingular integrals involved in BIEs a special regularization process that converts the hypersingular integrals to regular integrals is applied. Simple regular formulas for their calculation are presented. For the problems solution while considering frictional contact of the crack edges a special iterative algorithm of Udzava's type is elaborated and used. Numerical results for crack opening, frictional contact forces and dynamic stress intensity factors (SIFs) are presented and discussed for a finite III-mode crack in an infinite domain subjected to a harmonic crack-face loading and considering crack edges frictional contact interaction using the TD and FD approaches. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Fracture dynamics Crack edges contact BEM Time domain Frequency domain

1. Introduction Cracks and other structural defects are often found in materials used in engineering structures, apparatus and devises. These cracks occur in materials because of various reasons. They could appear as small flaws in the material manufacturing stage, they may arise during fabrication, or they may be the result of damage (fatigue, impact, corrosion etc.). In some situations, parts of the machines and structures contain significantly large cracks, nevertheless they work reliably. On the other hand, sometimes well-designed structures fail due to crack propagation under a significantly lower load than calculated. Therefore the development and improvement of machinery and structural design using fracture mechanics methods is very important. The inertial effects resulting from dynamic load and crack propagation need to be taken into account in many structural design situations that use methods of fracture mechanics. These two factors may occur separately or in combination. Examples are stationary cracks under dynamic loading (when the velocity of crack propagation is equal to zero) or propagating cracks under static loading. The main problem of dynamical fracture mechanics is the calculation of the SIF and J-integrals for cracked bodies. These and some other aspects of fracture dynamics methods and problems are presented in [4,7,14,38], etc.

n

Tel.: þ 52 9999428330. E-mail address: [email protected]

0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.08.006

It is important to point out that fracture dynamics problems are usually solved without taking into account the possibility of contact interaction between the opposite crack edges. Analysis of static fracture mechanics problems demonstrate that taking the crack edge contact interaction into account may significantly affect the fracture mechanics criteria. In dynamic problems, the effects of the crack edge contact interaction can significantly exceed those in the static case. Moreover, in dynamic problems it is very difficult to find classes of loads which do not cause crack edge contact interaction. The importance of taking into account the influence of crack edges contact interaction on fracture mechanics criteria has been investigated and discussed in our publications. For references see the book Guz and Zozulya [20], or review papers [21–24,58]. We would like to quote some publications where the problem of the crack edges contact interaction was formulated and methods for its solution were developed for the first time. Mathematical formulation of the elastodynamic problem for a cracked body, that takes into account the possibility of crack edge contact interaction and the formation of areas with close contact, adhesion and sliding, were reported for the first time in [44]. In [46] it was investigated specifically, yet very importantly, for the applications case of harmonic loading and representation of its solution as the Fourier series expansion was established. Algorithm for the problem solution with considerable unilateral crack edges contact interaction and friction was elaborated in [45] and adapted for the case of harmonic load in [46]. Algorithm is based on a theory of subdifferentional functionals and finding of their saddle points. Such algorithms are usually known as Uzava's type [5,37].

1500

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

Mathematical aspects of the problem and algorithm convergence were investigated in [48,50,57]. It has been shown that the algorithm may be considered a as compressive operator, acting on a special Sobolev's functional space. Briefly speaking the algorithm is the combination of two parts. The first part is the solution of an elastodynamic problem for bodies with cracks, without taking into account contact conditions. The second part is the projection of the founded solution on the set of unilateral contact restrictions and friction. The projection operators were constructed in [44,46], then used in a number of our publications, as well other authors publications. In the developed algorithm the first part is not directly defined, any numerical or even an analytical method for solution of the corresponding problem without taking into account contact conditions can be used. In our publications, the BIE method and its numerical implementation BEM were used. It was adapted for the problem solution mentioned above in [49]. Useful information on the BEM and its application to stress analysis, fracture mechanics and contact problems can be found in [3,6,10,17,26–28,33] The singularity of the integral operator kernels is one of the main problems that appears when the BIE is solved by the BEM. An approach which is based on the theory of distributions and the second Green's theorem is developed and applied for the divergent integral regularization in [47,51–53,55,56]. Simple regular formulas were elaborated and using these formulas the weakly singular, strongly and hypersingular 1-D and 2-D integrals can be considered the same way. In many important application cases, using the obtained regularized formulas divergent integrals can be calculated analytically, no numerical integration is needed. Three different BEM formulations, namely, the frequencydomain [8,9,18,40] the Laplace transform domain [13] and the time-domain BEM [12,15,19,29-31,35,36,41–43] are often applied to the elastodynamic crack analysis. Comparative study of the different BEM formulations and analysis of their accuracy and efficiency, is performed in [16,34] for elastodynamic problems without considering possibilities for contact integration. We do not know any publications, where such analysis was been done for the elastodynamic crack analysis problems by taking into account contact interaction. This paper has been written to address this shortcoming. A comparative study of the TD and the FD BEM formulations and the analysis of their accuracy and efficiency is performed in this paper for the case of the frictional contact problem for antiplane crack interacted with harmonic HS polarized waves. Numerical examples for computing the crack opening, the frictional contact forces and the dynamic SIFs are presented to compare the accuracy and the efficiency of the two different BEM formulations.

Fig. 1. Finite crack under antiplane deformation.

kinematic Cauchy relation ε3α ¼ ∂α u3

Substituting (1.2) into (1.1) the stress tensor components are defined by the displacement s3α ¼ μ∂α u3

ð1:3Þ

The stress differential equation of motion of the elastic body in this case has the form of μ∂β s3β þ b3 ¼ ρ∂2t u3

ð1:4Þ

Here and above ∂β ¼ ∂=∂xβ and ∂t ¼ ∂=∂t are derivatives with respect to the space coordinates and time, respectively, b3 is the volume force. Substituting s3α ðxα ; tÞ in Eq. (1.4) by its Hooke's law representation (1.3) we find the scalar wave equation for the displacement u3 ðxα ; tÞ to be in the form of μ∂β ∂β u3 þ b3 ¼ ρ∂2t u3

8 ðxα ; tÞ A V  ℑ

ð1:5Þ

If the problem is solved on an infinite region, then the solution for Eq. (1.5) is uniquely determined by assigning displacements and velocity vectors in the initial instant of time. Then the initial conditions are u3 ðxα ; t 0 Þ ¼ u03 ðxÞ;

∂t u3 ðxα ; t 0 Þ ¼ v03 ðxα Þ

8 xα A V

ð1:6Þ

Additional conditions at infinity have to be satisfied in this case

2. Formulation of the problem Let us consider a homogeneous linearly elastic body. It is well known (see [2,11]) that if the stress–strain state of elastic body depends on only the two coordinates xα ¼ ðx1 ; x2 Þ A V  R2 and time t A ℑ, then the main equations of elastodynamics are divided into two independent parts: plane and antiplane problem. Following [2,11] we consider the antiplane equations of elastodynamics. The elastodynamic stress–strain state in this case is defined by the following components of the stress s3α ðxα ; tÞ and strain ε3α ðxα ; tÞ tensors, which are related by Hooke's law s3α ¼ με3α

ð1:2Þ

ð1:1Þ

where μ is the shear modulus. The deformation is described by the shear component of the displacements u3 ðxα ; tÞ, which is related to the strain tensor by the

ð1:7Þ u3 ðxα ; tÞ ¼ Oðlnðr 1 ÞÞ; s3β ðxα ; tÞ ¼ Oðr 1 Þ for r-1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 Here r ¼ jxα j ¼ x1 þx2 is the distance in the 2-D Euclidian spaces. If the body occupies a finite region V, it is necessary to establish boundary conditions. We suppose that the region V is an open bounded subset of the 2-D Euclidean space ℜ2 with a C 1;1 Lyapunov's class regular boundary ∂V. The boundary contains two parts ∂V u and ∂V p such that ∂V u \ ∂V p ¼ ∅ and ∂V u [ ∂V p ¼ ∂V . On the part ∂Vu are prescribed displacements u3 ðxα ; tÞ of the body points and on the part ∂V p are prescribed tractions p3 ðxα ; tÞ respectively. Then the mixed boundary conditions are u3 ðxα ; tÞ ¼ ϕ3 ðxα ; tÞ

8 xα A ∂V u

8t Aℑ

p3 ðxα ; tÞ ¼ s3β ðxα ; tÞnβ ðxα Þ ¼ P n ½u3 ðxα ; tÞ ¼ ψ 3 ðxα ; tÞ

8 xα A ∂V p

8t Aℑ

ð1:8Þ

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

The differential operator P n : u3 -p3 is called traction operator. It transforms the displacements into the tractions. For the homogeneous isotropic elastic medium in the case of antiplane deformation has the form of P n ¼ μ∂n

ð1:9Þ

Here ∂n ¼ ni ∂i is a derivative in the direction of the vector nðxα Þ normal to the surface ∂V, ni are components of the outward unit normal vector. The body may contain an arbitrary oriented crack, which is described by its opposite surfaces Ω þ and Ω . Let us consider in more detail the unbounded homogeneous isotropic elastic body in 3-D Euclidean space with a finite crack located in the plane R2 ¼ fx : x3 ¼ 0g. We assume that displacements of the body points and their gradients are small, so the crack surface can be described by its Cartesian coordinates   Ω ¼ x : l r x1 r l; x2 ¼ 0; 1 r x3 r 1 ð1:10Þ A harmonic horizontally polarised shear SH-wave with frequency ω propagates in the plane R2 . The shear axis and axis Ox3 are coinciding as is shown in Fig. 1. The incident wave is defined by the potential function ψðxα ; tÞ ¼ ψ 0 eiðk2 n U xα ωtÞ ;

ð1:11Þ

where p ffiffiffiffiffiffiffiffi ψ 0 is the amplitude, k2 ¼ ω=c2 is the wave number, c2 ¼ μ=ρ is the velocity of the SH-wave, ω ¼ 2π=T is the frequency, T is the period of wave propagation, μ are the Lame constant, and ρ is the density of the material, n ¼ ð cos α ; sin αÞ is the unit vector, normal to the front of the incidence wave, α is the angle of the incident wave. This wave generates the stress–strain state that depends on two space coordinates xα  R2 and time t A ℑ which is called the antiplane deformation [2,11]. Wave propagation in a cracked body is a classical diffraction problem [2,11,20]. Usually this problem may be divided into two separate problems: the problem for incident waves and the problem for reflection waves. Obviously, the problem for incident wave is trivial in the case under consideration. If the wave function ψðxα ; tÞ is known, then the components of the stress tensor and displacements of the incident wave are determined in the form of the vector under action of u3 ¼ ∂2 ψ ¼ k2 n2 Refiψ 0 eiðk2 n U xα ωtÞ g; s3α ¼ μ∂α ∂2 ψ ¼ μk2 nα n2 Refψ 0 eiðk2 n U xα ωtÞ g 2

ð1:12Þ

Therefore we will pay more attention to the solution of the problem for reflected waves. On the crack's edges we have n1 ¼ 0 and x2 ¼ 0, therefore the load caused by the incident wave has the form iðk2 x1 ωtÞ

p3 ðxα ; tÞ ¼ s32 ðx1 ; tÞ ¼ p0 Refe

g;

p0 ¼

2 μk2 ψ 0

ð1:13Þ

We want to consider the possibility for the crack faces contact interaction. Therefore we suppose that an initial compressive load is applied to the crack surface in the x2 . In this case under the action of the harmonical loading frictional contact interaction of the opposite crack edges occurred. Considering the crack edges contact interaction, the load on the crack edges has the form of ( p3 ðxα ; tÞ 8 xα 2 = Ωe ps3 ðxα ; tÞ ¼ ð1:14Þ p3 ðxα ; tÞ þ q3 ðxα ; tÞ 8 xα A Ωe þ



where q3 ðxα ; tÞ is a frictional contact forces, Ωe ¼ Ω \ Ω is a region of close frictional contact, which varies with time. The force of the crack edges contact interaction q3 ðxα ; tÞ and displacement discontinuity Δu3 ðxα ; tÞ ¼ u3þ ðxα ; tÞu 3 ðx α ; tÞ should satisfy the contact constrains   q  r kτ q - ∂t Δu3 ¼ 0; n  3 q3  ¼ kτ qn - ∂t Δu3 ¼ λτ q3 8 xα A Ωe 8 t A ℑ ð1:15Þ

1501

where kτ and λτ are coefficients dependent on the contacting surfaces properties, qn ðxα ; tÞ is normal to the crack surface force of contact interaction. Here we use contact conditions (1.15) in the form of Coulomb friction, which are widely used for investigation of the elastodynamic contact problems with friction. For references see [5,20,37], also following [5,20,37,61,62], for the problem under consideration we assume that it is known beforehand. We will consider the solution of the above formulated elastodynamic contact problem for cracked body using the TD and FD BEM and analyze applicability, accuracy and efficiency, while comparing them and studying their advantages and disadvantages.

3. Integral equations and fundamental solutions in the TD The starting point in our consideration of the BIE in antiplane elastodynamics will be Betty–Rayleigh reciprocal theorem [2,11]. It presents the relation between two elastodynamic states of the elastic body. Usually one state refers to the main state and another to the secondary state. In order to distinguish between these two elastodynamic states we supplied the values which correspond to the secondary elastodynamic states with the mark “′”. For the case of zeros body forces, initial displacements and velocity the Betty– Rayleigh reciprocal theorem may be presented in the form of Z Z Z Z b3 ðxα ; τÞu′3 ðxα ; tτÞdV dτ þ p3 ðxα ; τÞu′3 ðxα ; tτÞdS dτ ℑ V ℑ ∂V [ Ω Z Z b03 ðxα ; tτÞu3 ðxα ; τÞdV dτ ¼ ℑ V Z Z p03 ðxα ; tτÞu3 ðxα ; τÞdS dτ ð2:1Þ þ ℑ ∂V [ Ω

From this theorem follows Somigliana's integral representation for the displacements. In order to obtain it we consider the main elastodynamic state, which corresponds to the problem under consideration and the secondary elastodynamic state that corresponds to an infinite region subjected to the unit impulse applied at the time t at the point xα . Such impulse can be represented by the Dirac delta functions, which depend on the space coordinates xα and уα and time in the form of b′i ðx; tτÞ ¼ δij δðxα уα ÞδðtτÞ

ð2:2Þ

We denoted the displacements and tractions that correspond to the secondary state in the infinite region as u′i ðxα ; tÞ : U ij ðxα уα ; tτÞ and p′i ðxα ; tÞ : W ij ðxα ; уα ; tτÞ

ð2:3Þ

Now the Betty–Rayleigh reciprocal theorem (2.1) is applied to two elastodynamic states: the actual state and the auxiliary state, which correspond to the action of the impulse load (2.2) in the infinite region. Taking into account properties of the delta function Somigliana's integral representation for the displacements is obtained in the form of Z Z  u3 ðyα ; tÞ ¼ p3 ðxα ; τÞU 3 ðxα yα ; tτÞ ℑ ∂V  u3 ðxα ; τÞW 3 ðxα ; yα ; tτÞ dS dτ Z Z Δu3 ðxα ; τÞW 3 ðxα ; yα ; tτÞdS dτ ð2:4Þ  ℑ Ω

where U 3 ðxα yα ; tτÞ is the fundamental solution for the wave equation (1.5), it corresponds to the displacement at point yα at time instant t, when the concentrated impulse (2.2) is applied at point xα at time instant τ, W 3 ðxα ; yα ; tτÞ is the fundamental solution for traction, it can be obtained from the fundamental solution for displacements applying traction operator (1.9). The fundamental solution for displacements can be found for example in [2,10,11]. It can be represented by the following

1502

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

equations:



U 3 ðxα yα ; t′Þ ¼

1 1 Hðt′r=c2 Þ 2πμ L

ð2:5Þ

   where Hðt′r=c2 Þ is the Heaviside step function, r xα ; yα ¼ xα yα  is the distance between the collocation point yα and the observation point xα . It is defined by the equation qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rðxα ; yα Þ ¼ ðx1 y1 Þ2 þðx2 y2 Þ2 ð2:6Þ For compactness we also introduce the following notations qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t′ ¼ tτ, L ¼ ðt′Þ2 ðr=c2 Þ2 . Applying traction operator (1.9) to the fundamental solution for displacements (2.5) we obtain fundamental solution for traction in the form of

 1 nα ðxα Þxα r Hðt′r=c2 Þ δðt′r=c2 Þ W 3 ðxα ; yα ; t′Þ ¼ ð2:7Þ  3 2πc2 c2 L r L Integral representation for traction can be obtained from Somigliano's integral representation (2.4) applying differential operator (1.9). As a result we have Z Z   p3 ðyα ; tÞ ¼ p3 ðxα ; τÞK 3 ðxα ; yα ; t′Þu3 ðxα ; τÞF 3 ðxα ; yα ; t′Þ dS dτ ℑ ∂V Z Z Δu3 ðxα ; τÞF 3 ðxα ; yα ; t′ÞdS dτ ð2:8Þ  ℑ Ω

Here the kernel K 3 ðxα ; yα ; t′Þ has the form of

 1 nβ ðyα Þxβ r Hðt′r=c2 Þ δðt′r=c2 Þ K 3 ðxα ; yα ; t′Þ ¼   2πc2 c2 L r L3

ð2:9Þ

The kernel F 3 ðxα ; yα ; t′Þ can be calculated applying the operator of the derivation in normal direction twice with respect to xα and yα to the fundamental solutions (2.5). As the result we obtain nα nβ xα xβ dU 3 ðxα yα Þ nα xα nβ xβ d U 3 ðxα yα Þ þ δαβ  2 dr r r r r dr 2 ð2:10Þ 2

F 3 ðxα ; yα ; t′Þ ¼

where dU 3 ðxα yα ; t′Þ r Hðt 0 r=c2 Þ δðt′r=c2 Þ ¼ 2  dr c2 L c2 L3 2 d U 3 ðxα yα ; t′Þ 2r 2 =c22 þ ðtτÞ2 Hðtτr=c2 Þ ¼ dr 2 c22 L5 2r δðtτr=c2 Þ δ′ðtτr=c2 Þ  3 þ c2 c22 L L3

ð2:11Þ

In the expressions for normal derivatives of the function  Eq. (2.10) r xα ; yα have the form of nα xα nα nβ xα xβ ð2:12Þ ; ∂n ð∂n rÞ ¼ δαβ  2 ∂n r ¼ r r r The term with delta function and its derivatives in (2.7), (2.9) and (2.11) possess singularity at t′ ¼ r=c2 , where the argument of the delta function cancels out. To avoid this problem following Gallego and Dominguez [15] we transform integral representations (2.4) and (2.8) using the relationship Z t Z t _ δðtτr=c2 Þf ðτÞuðτÞdτ ¼ Hðtτr=c2 Þðf_ ðτÞuðτÞ þ f ðτÞuðτÞÞdτ 0



Z Z h Zℑ Z∂Vh ℑ Ω

i _ n ðxα ; y ; t′Þ dS dτ Δu3 ðxα ; τÞW n3 ðxα ; yα ; t′Þ þ Δu_ 3 ðxα ; τÞW α 3

i _ n ðxα ; y ; t′Þ dS dτ Δu3 ðxα ; τÞW n3 ðxα ; yα ; t′Þ þΔu_ 3 ðxα ; τÞW α 3 ð2:14Þ

where 1 nα ðxα Þxα t′r=c2 Hðt′r=c2 Þ W n3 ðxα ; yα ; t′Þ ¼  2πc2 r L3 _ n ðxα ; y ; t′Þ ¼  1 nα ðxα Þxα 1 Hðt′r=c2 Þ W α 3 2πc2 L r

ð2:15Þ

In the same way the integral representation for traction (2.8) can be represented in the form of Z Z h n p3 ðxα ; τÞK n3 ðxα ; yα ; t′Þ þ p_ 3 ðxα ; τÞK_ 3 ðxα ; yα ; t′Þ p3 ðyα ; tÞ ¼ ℑ ∂V

n

u3 ðxα ; τÞF n3 ðxα ; yα ; t′Þu_ 3 ðxα ; τÞF_ 3 ðxα ; yα ; t′Þ i n u€ 3 ðxα ; τÞF€ 3 ðxα ; yα ; t′Þ dS dτ Z Z h n  Δu3 ðxα ; τÞF n3 ðxα ; yα ; t′ÞΔu_ 3 ðxα ; τÞF_ 3 ðxα ; yα ; t′Þ ℑ Ω i n ð2:16Þ Δu€ 3 ðxα ; τÞF€ 3 ðxα ; yα ; t′Þ dS dτ where μ nβ ðyα Þxβ K n3 ðxα ; yα ; t′Þ ¼  2πc2 r μ n ðy n β α Þxβ K_ 3 ðxα ; yα ; t′Þ ¼  2πc2 r

t′r=c2 L3

Hðt′r=c2 Þ

1 Hðt′r=c2 Þ L

ð2:17Þ

and

μ ðt′r=c2 Þ nα ðxα Þxα nβ ðyα Þxβ 1 F n3 ðxα ; yα ; t′Þ ¼  3 2πc2 r r rL

 3r Hðt′r=c2 Þ;  1þ c2 t′þ rÞ

μ 1 nα ðxα Þxα nβ ðyα Þxβ n 1 F_ 3 ðxα ; yα ; t′Þ ¼  2πc2 rL r r

 2r  1þ Hðt′r=c2 Þ; c2 t′þ r n μ nα ðxα Þxα nβ ðyα Þxβ 1 F€ 3 ðxα ; yα ; t′Þ ¼  2πc2 c2 L r r

ð2:18Þ

In the case of the straight crack in the unbounded plane presented in Fig. 1 all the above equations are significantly simplified. Indeed in this case n1 ¼ 0, n2 ¼ 1, x2 ¼ 0 and ∂V ¼ ∅. Therefore the integral equation that related the crack opening and the traction is presented in the form of Z Z n p3 ðyα ; tÞ ¼  ½Δu3 ðxα ; τÞF n3 ðxα ; yα ; t′ÞΔu_ 3 ðxα ; τÞF_ 3 ðxα ; yα ; t′ÞdS dτ ℑ Ω

ð2:19Þ with F n3 ðxα ; yα ; t′Þ ¼

μ ðt′r=c2 Þ ; 2πc2 r L3

n F_ 3 ðxα ; yα ; t′Þ ¼

μ 1 2πc2 r L

ð2:20Þ

These simplified equations will be used for the numerical solution comparative study of the problem that is considered here.

0

ð2:13Þ Here and throughout this paper the point on the top means derivative with respect to time. Then the integral representation for displacements (2.4) can be presented in the form of Z Z u3 ðyα ; tÞ ¼ ½p3 ðxα ; τÞU 3 ðxα yα ; t′Þ ℑ ∂V

n

_ ðxα ; y ; t′ÞdS dτ u3 ðxα ; τÞW n3 ðxα ; yα ; t′Þu_ 3 ðxα ; τÞW α 3

4. Discretization and BEM in the TD The BEM can be treated as the approximate method for the BIE solution, which includes an approximation of the functions that belong to some functional space, boundary of the body including crack surfaces and time interval by a discrete finite model. In our specific case of the finite crack in an unbounded region the boundary consists of the crack surface Ω which is divided into N boundary elements with Q nodes of interpolation per element.

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

Standard Lagrange interpolation polynomials have the form of ϕq ðξÞ ¼

ðξξ0 Þðξξ1 Þ…ðξξq1 Þðξξq þ 1 Þ…ðξξQ Þ ðξq ξ0 Þðξq ξ1 Þ…ðξq ξq1 Þðξq ξq þ 1 Þ…ðξq ξQ Þ

ð3:1Þ

For piece-wise constant approximation ( 1 8 x α A Ωn ; ϕq ðxα Þ ¼ = Ωn : 0 8 xα 2

ð3:2Þ

The observation time interval ℑ ¼ ½t 0 ; T is divided into M time steps with P nodes of interpolation per time step For piece-wise constant approximation ψ 0 ðtÞ ¼ Hðtt m þ ΔtÞHðtt m Þ

ð3:3Þ

where HðtÞ represents the Heaviside unit step function. For piece-wise linear approximation 1 ψ 1 ðtÞ ¼ ðRðtt m þ ΔtÞ2Rðtt m Þ þ Rðtt m ΔtÞÞ Δt

Q

Δu3 ðxα ; τÞ ¼ ∑

x α A Ωn ;

τ A ½t m1 ; t m ;

∑ Δu3 ðxqα ; t m Þϕq ðxÞψ_ p ðτÞ;

x α A Ωn ;

τ A ½t m1 ; t m ;

p¼1q¼1 Q

P

Δu_ 3 ðxα ; τÞ ¼ ∑

p¼1q¼1

P

p3 ðxα ; τÞ ¼ ∑

Q

∑ p3 ðxqα ; t m Þϕq ðxÞψ p ðτÞ;

p¼1q¼1

xα A Ωn ;

τ A ½t m1 ; t m  ð3:5Þ

Then TD discrete BEM for antiplane elastodynamics problems for straight crack in an infinite plane has the form of N

p3 ðyr ; t M Þ ¼  ∑

Q

M

∑

P

∑ ðF nq 3 ðy r ; x n ; t M t m ÞΔu3 ðx n ; t m Þ

∑

n¼1q¼1m¼0p¼1

nq F_ 3 ðyr ; xn ; t M t m ÞΔu_ 3 ðxn ; t m ÞÞ

ð3:6Þ

where F q;p 3 ðyr ; xn ; t M t m Þ ¼ q;p F_ 3 ðyr ; xn ; t M t m Þ ¼

Z

Z

tm

t m1

Z

Z

tm

Ωn

Ωn

t m1

 F n3 yr ; x; t M τ ϕq ðξÞdSψ p ðτÞdτ; n F_ 3 yr ; x; t M τ ϕq ðξÞdSψ_ p ðτÞdτ;

ð3:7Þ

The integration in (3.7) can be done analytically when the interpolation functions are simple enough. Therefore we consider here the piece-wise constant BE and corresponding interpolation polynomial in the form (3.2) and the two type of time interpolation: piece-wise (3.3) constant and piece-wise linear (3.4). First, we will evaluate the time integrals in (3.7), changing the order of space and time integration we have Z Z tm F q;p ðy ; x ; t t Þ ¼ F n3 ðyr ; x; t M τÞψ p ðτÞdτϕq ðξÞdS n m M r 3 q;p F_ 3 ðyr ; xn ; t M t m Þ ¼

Z

Ωn

Ωn

t m1

Z

tm

t m1

n F_ 3 ðyr ; x; t M τÞψ_ p ðτÞdτϕq ðξÞdS

ð3:8Þ

In the case of piece-wise constant time approximation interpolation polynomial and its derivative with respect to time are ψ 0 ðtÞ ¼ Hðtt m þ 1 ÞHðtt m Þ;

δðτt m þ 1 Þδðτt m ÞC þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AHðt M τr=c2 Þdτ ðt M τÞ2 ðr=c2 Þ2

ð3:10Þ

Three different cases arise while carrying out integration depending on propagation of the wave along the crack length. 1. For t M t m1 r r=c2 , the effect of perturbation has not arrived yet at the point xα ð3:11Þ

ð3:4Þ

Δu3 ðxqα ; t m Þϕq ðxÞψ p ðτÞ;

∑

Substituting them and kernels (2.20) in (3.8) we obtain integral that has to be evaluated in the form of 0 Z tm μ ðt M τÞr=c2 B M;m F 3 ðxα ; yα ; t M Þ ¼ @

3=2 2πc2 r tm 1 ðt M τÞ2 ðr=c2 Þ2 1

F M;m ðxα ; yα ; t M Þ ¼ 0 3

where RðtÞ represents the unit ramp function, defined as RðtÞ ¼ t if t 4 0 and 0 otherwise. Now the crack opening, its velocity and traction on the boundary element Ωn and the time step ½t m1 ; t m  can be presented in the form of P

1503

ψ_ 0 ðtÞ ¼ δðtt m þ 1 Þδðtt m Þ

ð3:9Þ

2. For t M t m rr=c2 o t M t m1 , the effect of the part of the perturbation has arrived at xα but the part of the perturbation produced for r such that r=c2 4 t n τ has not arrived at xα yet μ ðt M t m1 Þ ðxα ; yα ; t M Þ ¼  F M;m 3 Lm1 2πr 2 Here Lm ¼

ð3:12Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðt M t m Þ2 ðr=c2 Þ2

3. For t n t m o r=c2 , the effect of part of the perturbation has arrived at xα but the part of the perturbation produced for r such that r=c2 4 t M τ has not arrived at xα yet

 μ ðt M t m1 Þ ðt M t m Þ ðx ; y ; t Þ ¼   F M;m α M α 3 Lm1 Lm 2πr 2

ð3:13Þ

Now we proceed with space integration. Considering (3.12) and (3.13) we have to evaluate the following integral: Z Δn Z Δn μ 1 c2 ðt M t m1 Þ F M:m ðx1 ; yr1 ; t M Þdx1 ¼  dx1 ð3:14Þ 3 2π Δn ðx1 yr1 Þ2 Lm Δn In the case if r an integral is regular and can be evaluated numerically or analytically without any problems. Analytical evaluation is given as

 þ μ 1 Lm1 L n;r n;r n;r m1 F M;m ðr ; t Þ ¼ ðr Δ Þ ðr þ Δ Þ n n M 3 π ðr n;r Þ2 þ Δ2n t M t m1 t M t m1 ð3:15Þ Here

  r n;r ¼ xn1 yr1 ; and 7 ¼ Lm

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðt M t m1 Þ2 ððr n;r 7 Δn Þ=c2 Þ2

In the case where r ¼ n x1 A ½yn1 Δn ; yn1 þΔn  integral (3.14) is hypersingular, special treatment is needed for its evaluation. Here we use the method of the divergent integrals regularization developed in our publications [47,51,53,55]. For that purpose we present integral (3.14) in the form of Z Δn F:P: F M:m ðx1 ; yr1 ; t M Þdx1 3 Δn

μ ¼  F:P: 2π

Z

Δn

Δn

ϕðx1 Þ ðx1 yr1 Þ2

dx1 ;

x1 A ½Δn ; Δn 

ð3:16Þ

1504

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

μ ¼ 2πc2 Δt

where c2 ðt M t m1 Þ ϕðx1 Þ ¼ Lm Now the integration by parts in the sense of generalized functions gives a regularized formula for the hypersingular integral (3.16) calculation in the form of sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

2 2μ 1 Δn M;m n;n F 3 ðr ; t M Þ ¼  1 ð3:17Þ π Δn c2 ðt M t m1 Þ In the case of piece-wise linear time approximation interpolation polynomial and its derivative with respect to time are 1 ðRðtt m þ ΔtÞ2Rðtt m Þ þ Rðtt m ΔtÞÞ; Δt 1 ψ_ 1 ðtÞ ¼ ðHðtt m þ ΔtÞ2Hðtt m Þ þ Hðtt m ΔtÞÞ Δt ψ 1 ðtÞ ¼

ð3:18Þ

Substituting them and kernels (2.20) in (3.8) we obtain an integral that has to be evaluated in the form of Z tm μ ðτt m1 Þ F M;m ðx ; y ; t Þ ¼ α M α 3 2πc2 r tm1 Δt 0 1 B @

ðt M τÞr=c2

δðτt m þ 1 Þδðτt m ÞC

3=2 þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A ðt M τÞ2 ðr=c2 Þ2 ðt M τÞ ðr=c2 Þ 2

0 B @

μ 2πc2 r

Z

tm þ 1 tm

ðt m þ 1 τÞ Δt

1

ð3:19Þ

In this case four different cases arise in carrying out integration depending on the propagation of the wave along the crack length: 1. For t n t m1 r r=c2 , Δtðnm þ 1Þ r r=c2 , ðt n t m1 ¼ Δtðnm þ 1ÞÞ the effect of the perturbation has not yet arrived at point xα ð3:20Þ

2. For t n t m rr=c2 o t n t m1 , the effect of the part of the perturbation has arrived at xα but the part of the perturbation produced for r such as r=c2 4t n τ has not arrived at xα yet μ Lm1 ðxα ; yα ; t M Þ ¼  ð3:21Þ F M;m 3 2πr 2 Δt 3. For t n t m þ 1 o r=c2 rt n t m , the effect of part of the perturbation has arrived at xα but the part of the perturbation produced for r such as r=c2 4 t n τ has not arrived at xα yet μ ðLm1 2Lm Þ ðxα ; yα ; t n Þ ¼  ð3:22Þ F M;m 3 2πr 2 Δt 4. For r=c2 r t n t m þ 1 , the effect of part of the perturbation has arrived at xα but the part of the perturbation produced for r such as r=c2 4 t n τ has not arrived at xα yet μ ðLm1 2Lm þ Lm þ 1 Þ ðxα ; yα ; t n Þ ¼  F M;m 3 2πr 2 Δt

ð3:23Þ

Now we proceed with space integration. When considering (3.21)–(3.23) we have to evaluate the following integral: Z Δn F M;m ðx1 ; yr1 ; t M Þdx1 3 Δn

r 2 Δn ðx1 y1 Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c22 ðt M t m1 Þ2 ðx1 yr1 Þ2 dx1

ð3:24Þ

In the case that r a n the integral is regular and can be evaluated numerically or analytically without any problem. The analytical evaluation is given as



n;r  μ Lm1 r Δn n;r ðr ; t Þ þ arctan F M:m M  3 2πΔt r n;r Δn c2 Lm1

n;r  þ Lm1 r þ Δn arctan ð3:25Þ  n;r þ r þ Δn c2 Lm1 In the case that r ¼ n x1 A ½yn1 Δn ; yn1 þΔn integral (3.24) is hypersingular, special treatment is needed for its evaluation. Using the same approach as in the previous case we obtain a regularized formula for calculation of the hypersingular integral (3.24) in the form of sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi μ ðt M t m1 Þ Δn n;n 1 ð3:26Þ F M:m ðr ; t Þ ¼  M 3 2πΔn Δt c2 ðt M t m1 Þ The discrete BEM equation (3.6) by considering numerical calculation of the Riemann convolution integrals in (3.7) can be presented in the following matrix form:[10,19,42] ð3:27Þ

m¼0

δðτt m þ ΔtÞδðτt m ÞC

3=2 þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A 2 2 ðt M τÞ2 ðr=c2 Þ2 ðt M τÞ ðr=c2 Þ

F n3 ðxα yα ; t n Þ ¼ 0

1

M

ðt M τÞr=c2

Hðt M τr=c2 Þdτ

Δn

pM ¼  ∑ AM;m U Δum

2

Hðt M τr=c2 Þdτ þ

Z

where Mis a number of time steps and N is a number of BE.    p3 ðy1 ; t M Þ      ⋮ pM ¼  ;    p3 ðyN ; t M Þ     F M;m y ; x ; t t ⋯ F M;m ðy ; x ; t t Þ   3 m m  N M 1 M 1 1 3   ; ⋮ ⋮ AM;m ¼     F M;m y ; x ; t t M;m ⋯ F 3 ðyN ; xN ; t M t m Þ   3 m 1 M N    Δu3 ðy1 ; t m Þ      ⋮ Δum ¼  ð3:28Þ     Δu3 ðyN ; t m Þ  Now, taking into account the initial conditions from Eq. (3.27) the explicit time-stepping algorithm follows

 M1 ΔuM ¼ ðA0 Þ1 pM  ∑ AMm þ 1 U Δum ; ð3:29Þ m¼1

0 1

here ðA Þ is the inverse of the system matrix A0 at the time-step M ¼ 0. The main specific feature of the problem considered here is the presence of the frictional contact conditions (1.15) which make the problem strongly nonlinear. Special algorithms have to be used for such problems solution. Here we use Uzawa's type algorithm for the first time proposed in [45] and then further developed in [59,60] for the solution of elastodynamic contact problems for bodies with cracks. Each time interval algorithm consists of the following steps:

 the initial distribution of the contact forces q03 ðx1 ; tÞ ; 8 x A Ω ; 8 t A ℑ is assigned;

 the problem without constrains is solved and the unknown quantities on the contact surfaces Δu3 ðx1 ; tÞ are defined;

 the normal and tangential components of the vector of the contact forces are corrected to satisfy the unilateral restrictions qk3 þ 1 ðx1 ; tÞ ¼ P τ ½qk3 ðx1 ; tÞρτ ∂t Δuk3 þ 1 ðx1 ; tÞ; where   P τ q3 ¼

(

q3 ðx3 ; tÞ kτ qn ðx3 ; tÞq3 =jq3 j

if jq3 j r kτ qn ðx3 ; tÞ   if q  4 kτ q ðx3 ; tÞ 3

n

ð3:30Þ

ð3:31Þ

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513



  is the operator of the orthogonal projection onto set q3  r kτ qn ðx3 ; tÞ, coefficient ρτ has been chosen based on the conditions that give the best convergence of the algorithm; proceed to the next step of the iteration.

Algorithms of such kind are widely used to solve different constrained optimization problems; for example, elastostatic contact problems (see [5,37]). In the elastostatic case the algorithm has well established theoretical basis with proof of the existence, uniqueness and convergence of the solution (see [54]). In the case of dynamical problems no well-established results of the existence, uniqueness and convergence are available. Variational formulation by considering non-smooth functional relations and investigation of such algorithms convergence in the case of elastodynamic contact problems have been presented in [57]. Particularly we find that coefficient ρτ has a significant influence on algorithm convergence. In the case that ρτ is too small algorithm can converge too slow, in the case that it is too big algorithm can be divergent.

5. Integral equations and fundamental solutions in the FD Let us consider the problem of a harmonic horizontally polarised shear SH-wave with frequency ω propagating in the plane R2 in the FD. In the problem related to reflected waves load is applied to the crack edges as it is shown in Fig. 1. Because of contact constrains (1.15) the problem under consideration is nonlinear. Therefore, as it was shown in [50] and our other publications, the problem for reflected waves presents the periodic steady state process, but not the harmonic process. As a result, components of the stress–strain state, caused by the reflected waves cannot be represented as functions of coordinates xα , multiplied by factor eiωt , as it usually done in elastodynamics in the case of the harmonic loading [2,11]. That is why we have to expand components of the displacement vector and stress tensor into Fourier series with the parameter of loading ω 1  u3 ðxα ; tÞ ¼ Re ∑ uk3 ðxα Þeiωk t ; 1

1  s3β ðxα ; tÞ ¼ Re ∑ sk3β ðxα Þeiωk t

ð4:1Þ

1

uk3 ðxα Þ ¼

ωk 2π

sk3β ðxα Þ ¼

T

u3 ðxα ; tÞeiωk t dt;

0

ωk 2π

Z

T 0

s3β ðxα ; tÞeiωk t dt

ð4:2Þ

In the same way, traction p3 ðx1 ; tÞ on the crack edges and their opening Δu3 ðx1 ; tÞ may be expanded into Fourier series 1  p3 ðx1 ; tÞ ¼ Re ∑ pk3 ðx1 Þeiωk t ; 1

1  Δu3 ðx1 ; tÞ ¼ Re ∑ Δuk3 ðx1 Þeiωk t

ð4:3Þ

1

pk3 ðx1 Þ ¼

ωk 2π

Δuk3 ðx1 Þ ¼

ωk 2π

T

0

p3 ðx1 ; tÞeiωk t dt

Z

T 0

Δu3 ðx1 ; tÞeiωk t dt

∂β ∂β uk3 þ

ω2k k ρ k u þ b ¼0 c22 3 μ 3

for k ¼ 0; 7 1; 7 2; :::; 1

8 xα A Ω

ð4:5Þ

In [20,21,23, 24] it was shown that the Fourier series expansions of the displacement discontinuity Δuk3 ðxα Þ and traction pk3 ðxα Þ are related by the boundary integral equations (BIE) in the form of Z pk3 ðxα Þ ¼ F:P: F 3 ðxα yα ; ωk ÞΔuk3 ðyα ÞdΩ Ω

for k ¼ 0; 7 1; 7 2; :::; 7 1

8 xα A Ω

ð4:6Þ

The kernels F 3 ðxα yα ; ωk Þ may be obtained from the fundamental solutions U 3 ðxα yα ; ωk Þ for the steady-state wave equations. The fundamental solutions for the steady-state wave equation is well known and may be found anywhere, for example in [2,11,10]. It is in the form of U 3 ðxα yα ; ωk Þ ¼

i ð1Þ k H ðl Þ; 4μ 0 2

k

l2 ¼ rωk =c2

ð4:7Þ

where H 10 ðzÞ is the Bessel function of the third kind and of zero order (Hankel function) [1]. Applying the operator of derivation in normal direction (1.9) twice with respect to xα and yα respectively to the fundamental solutions U 3 ðxα yα ; ωk Þ we obtain F 3 ðxα ; yα ; ωk Þ ¼ μ2 ∂n ∂n U 3 ðxα yα ; ωk Þ

ð4:8Þ

where ∂n ¼ nα ∂α is the normal derivative. The normal derivatives of the fundamental U 3 ðxα yα ; ωk Þ are calculated by the equations

solutions

nα xα dU 3 ðr; ωk Þ ; dr r

nα nβ xα xβ dU 3 ðxα yα Þ ∂n ∂n U 3 ðxα yα Þ ¼ δαβ  2 dr r r

∂n U 3 ðxα yα ; ωk Þ ¼

2

þ

nα xα nβ xβ d U 3 ðxα yα Þ r r dr 2

ð4:9Þ

Taking into account that in the case of the plane finite crack shown in Fig. 1 vectors of outward normal are n1 ¼ 0, n2 ¼ 1 and x2 ¼ 0, the normal derivative for the fundamental solution U 3 ðxα yα ; ωk Þ have the form of ∂n ∂n U 3 ðxα yα ; ωk Þ ¼

dU 3 ðr; ωk Þ r dr

ð4:10Þ

we obtain the kernel F 3 ðxα yα ; ωk Þ. For our specific case it has the form of iω k F 3 ðxα yα ; ωk Þ ¼ μ k H ð1Þ ðl Þ rc2 1 1

ð4:12Þ

The kernel is a complex value function. It may be represented in the form of Im F 3 ðxα yα ; ωk Þ ¼ F Re 3 ðx α y α ; ωk Þ þ iF 3 ðx α yα ; ωk Þ

ð4:13Þ

In order to separate real and imaginary parts of the fundamental solutions we substitute the Hankel functions by the Bessel functions of the first and the second kind using the relation

where Z

from of

Using the formula for the derivative of the Bessel function H 10 ðzÞ) [1]



 d ð1Þ ωr ω ωr H0 ¼  H ð1Þ ð4:11Þ dr c2 c2 1 c2

where Z

1505

ð4:4Þ

Inserting the Fourier series expansions defined by (4.1) into governing equation (1.5), instead of the one wave equation we obtain a countable set of steady-state wave equations in the

H ð1Þ ν ¼ J v ðzÞ þ iY v ðzÞ

ð4:14Þ

Now the real and imaginary parts of the fundamental solution have the form of ωk Y 1 ðl2 Þ; F Re 3 ðx α yα ; ωk Þ ¼ μ rc2

F Im 3 ðx α y α ; ωk Þ ¼ μ

ωk J ðl2 Þ rc2 1

ð4:15Þ

1506

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

Let us consider in detail the structure of the kernel F 3 ðxα yα ; ωk Þ and analyze its behavior when xα -yα . The Bessel functions of the first and second kind and the first order may be represented by the series expansion [1]   z 1 z 2 1 z 4 J 1 ðzÞ ¼ 1 þ ⋯ 2 1!2! 2 2!3! 2  2 z 1 z 2 1ð1 þ 1=2Þ Y 1 ðzÞ ¼ ðlnðz=2Þ þ γÞJ 1 ðzÞ þ π 2π 1!2! 2  1 z 4 2 ð4:16Þ ⋯  þ ð2 þ 1 þ1=3 2!3! 2 πz From these representations it follows that J 1 ðzÞ-z and Y 1 ðzÞ-1=z for z-0 and therefore for 1 F Re 3 ðx α yα ; ωk Þ- 2 ; r

F Im 3 ðx α y α ; ωk Þ-0

for xα -yα

ð4:17Þ

The above analysis follows that the kernel (4.15) of the integral equation (4.6) is hypersingular and special treatment is needed for its calculation. We will consider it in the sense of the Hadamard's finite part using the method of divergent integrals regularization based on the theory of distributions developed in our publications [47,51,53,55] for 2-D and 3-D elastodynamic problems. From the above analysis it follows that the BIEs (4.6) establishes a relation between complex functions pk3 ðxα Þ and Δuk3 ðxα Þ. Separating real and imaginary parts

will be used. For that purpose we present the kernel F Re 3 ðx α yα ; ωk Þ in the form of Z Z St n r r r r ½F Re F St F Re k ðx y ; Δn Þ ¼ 3 ðxy ; kÞF 3 ðxy Þdx þ F:P: 3 ðxy Þdx 2Δn

ð4:22Þ F St 3 ðx α yα Þ ¼

4μ Δn n r F St 3 ðx y ; Δn Þ ¼  π ðxn yr Þ2 Δ2n

Ω

The BIEs (4.19) are solved numerically, transforming them into a linear finite-dimensional system of the BEM equations using the collocation method. For zero order interpolation polynomial (3.2) with piecewise constant approximation of the boundary elements, the system of boundary element equations has the form of

Δn

Y n1 ðm; n; cα Þ ¼

N

n¼1

Im F Re 3 ðx n yr ; ωk ÞΔu3 ðxn ; kÞ

Their coefficients are defined by the equations Z F Re F Re k ðx n yr ; Δn Þ ¼ 3 ðxy r ; ωk ÞdΩ; Z2Δn F Im F Im k ðx n y r ; Δn Þ ¼ 3 ðxy r ; ωk ÞdΩ

ð4:20Þ

ð4:21Þ

2Δn

Let us consider in detail how to calculate these coefficients. If points xα and yα belong to different boundary elements, the integrals in (4.21) are non-singular and their calculation gives no difficulties. These can be calculated analytically or numerically using the Gauss quadrature formulas. If points xα and yα , belong to the same boundary elements some integrals in (4.21) are hypersingular. For calculation of these hypersingular integrals above mentioned regularization method

Z

Δn Δn

2

1

2Δn þ ∑ ð1Þ k¼0

mn k γ k ðΔn ; cα Þ

#

ðk!Þðk þ 1Þ!

n

ðω=cα rÞY 1 ðlα Þdβ ¼ γJ n1 ðm; n; cα Þðω=2cα Þ2 1

 2Δn þ ∑ ð1Þk k¼0

! γ mn k ðΔn ; cα Þ ðψðk þ 1Þ þ ψðk þ 2ÞÞ k!ðk þ 1Þ!

ð1Þk k ¼ 0 k!ðk þ 1Þ! 1

þ 2ðω=2cα Þ2 ∑

mn ð4:24Þ ðηmn k ðΔn ; cα Þγ k ðΔn ; cα Þ=ð2k þ 1ÞÞ       where r n ¼ ym xn þ β, r þ Δ ¼ ym xn þ Δn , r Δ ¼ ym xn Δn  and

Z γ mn k ðΔn ; cα Þ ¼ ηmn k ðΔn ; cα Þ ¼

N

Im Re pIm 3 ðyr ; kÞ ¼ ∑ ½F 3 ðx n yr ; ωk ÞΔu3 ðx n ; kÞ

"

¼ ðω=2cα Þ

Re Re pRe 3 ðy r ; ωk Þ ¼ ∑ ½F 3 ðx n y r ; ωk ÞΔu3 ðx n ; kÞ n¼1 Im F Im k ðx n y r ; ωk ÞΔu3 ðx n ; kÞ;

ð4:23Þ

The residuary integrals in (4.21) are regular and can be calculated using standard approaches. Since we consider the rectilinear crack, all those regular integrals can be calculated analytically using the methodology presented for the first time in [48], and then have been used in numerous of our publications (see for example [20, 22–25]) Z Δn n J n1 ðm; n; cα Þ ¼ ðω=cα rÞJ 1 ðlα Þdβ

ð4:18Þ

and using representation (4.13) for the kernels F 3 ðxα yα ; ωk Þ we obtain the system of integral equations for each k in the form of Z Re pRe F Re 3 ðy α ; kÞ ¼ F:P: 3 ðx α y α ; ωk ÞΔu3 ðx α ; kÞdΩ Z Ω Im  F Im 3 ðx α yα ; ωk ÞΔu3 ðxα ; kÞdΩ; Ω Z Re F Im pIm 3 ðyα ; kÞ ¼ 3 ðx α y α ; ωk ÞΔu3 ðx α ; kÞdΩ Ω Z Im ð4:19Þ þF:P: F Re 3 ðx α yα ; ωk ÞΔu3 ðx α ; kÞdΩ

2

Here 2μ=πr is the fundamental solution for the elastostatic problem, extracted from F Re 3 ðx α yα ; kÞ by considering the series expansion (4.16) and the asymptotic analysis for the case of ωk -0. The elastostatic fundamental solution is hypersingular. By applying the regularization technique mentioned above and calculating the finite part integral according to Hadamard we obtain

Im pk3 ðxα Þ ¼ pRe 3 ðx α ; kÞ þ ip3 ðx α ; kÞ; Im Δuk3 ðxα Þ ¼ ΔuRe 3 ðx α ; kÞ þ iΔu3 ðx α ; kÞ

2Δn

Δn

Δn

ðr n ω=2cα Þ2k dβ ¼

ðω=2cα Þ2k 2k þ 1 þ1 ðr þ Δ þ r 2k Δ Þ; 2k þ 1

ðω=2cα Þ2k 2k þ 1 þ1 ðr þ Δ lnðr þ Δ ω=2cα Þ þ r 2k lnðr Δ ω=2cα ÞÞ Δ 2k þ 1 ð4:25Þ

The use of these equations significantly reduces the time spent on the kernels (4.21) calculation. In contrast with other approaches (see for example [32]) the methods developed here demonstrate good accuracy regardless of the frequency ωk being low or high. The system BEM equation (4.20) does not lose stability with increasing ωk . For convenience and compactness the BEM equation (4.20) can be presented in the form of pk ¼ Ak UΔuk

for k ¼ 0; 7 1; 7 2; :::; 1

ð4:26Þ

where  Re   p ðy1 ; ωk Þ   3    ⋮      pRe ðy ; ωk Þ  N   3 k p ¼  Im ;  p3 ðy1 ; ωk Þ      ⋮    Im   p3 ðyN ; ωk Þ 

   ΔuRe ðy1 ; ωk Þ  3     ⋮      ΔuRe ðy ; ωk Þ  N   3 k Δu ¼   ðy1 ; ωk Þ   ΔuIm 3     ⋮     Im  Δu3 ðyN ; ωk Þ 

ð4:27Þ

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

Crack opening Δu3

1.4

1507

Crack opening Δu3

1.4 1. 100 BE 2. 10 BE

1. 2. 3. 4. 5.

4

1.2

1.2 2

5

1

t=0.15 t=0.30 t=0.45 t=1.0 t=4.0

1 1

0.8

0.8

0.6

0.6

3

2 0.4

0.4

0.2

0.2

1

0

0

1

2

3

4

0 −1

−0.5

Time t

0

0.5

1

Normalized crack length Fig. 2. Crack opening in TD for impulse load.

1. 10 BE 2. 50 BE 1. 100 BE

1.2

1. 10 BE 2. 50 BE 3. 100 BE

1.2

1

1 3

0.8

0.8

2

3

0.6

0.6

1 2 0.4

0.4 1

0.2

0.2

0 −1

−0.5

0

0.5

1

0 −1

−0.5

0

0.5

1

Normalized crack length

Normalized crack length

Fig. 3. Crack opening in the TD for different number of BE.

 Re  F ðx1 y1 ; ωk Þ  3  ⋮   Re  F ðx1 y ; ωk Þ N 3  Ak ¼  Im  F 3 ðx1 y1 ; ωk Þ   ⋮   Im  F 3 ðx1 yN ; ωk Þ

⋯ ⋯ ⋯ ⋯

F Re 3 ðxN y 1 ; ωk Þ ⋮ F Re 3 ðx N yN ; ωk Þ

F Im 3 ðx 1 y 1 ; ωk Þ ⋮ F Im 3 ðx 1 yN ; ωk Þ

F Im 3 ðx N y1 ; ωk Þ ⋮ F Im 3 ðx N y N ; ωk Þ

F Re 3 ðx 1 y1 ; ωk Þ ⋮ F Re 3 ðx1 y N ; ωk Þ

⋯ ⋯ ⋯ ⋯

  F Im 3 ðx N y1 ; ωk Þ   ⋮    F Im 3 ðx N y N ; ωk Þ   Re F 3 ðxN y1 ; ωk Þ   ⋮   Re F 3 ðxN yN ; ωk Þ 

ð4:28Þ

The algorithm for the solution of the problem in the FD includes the following steps:

 the initial distribution of the contact forces q03 ðx1 ; tÞ in the TD is assigned;

 the Fourier coefficients q3 ðx1 ; ωk Þ in the FD are calculated using formula (4.4);

 the FD BEM equations (4.20) and (4.26) are solved and Δu3 ðx1 ; ωk Þis calculated in the FD;

The problem under consideration has been transformed into the BIE (4.19), which has to be solved taking the unilateral restrictions into account (1.15). As it was mentioned before the main specific feature of the problem considered here is the presence of the frictional contact conditions (1.15) which make the problem strongly nonlinear. Here, we use the same Uzawa's type algorithm developed in [45] and slightly modified for the case of the FD solution in [46], more information about algorithms of such type and their convergence can be found in [57].

 the TD value Δun3 ðx; tÞ is calculated using formula (4.3);  the tangential components of the vector of contact forces are corrected to satisfy the unilateral restrictions qk3 þ 1 ðx1 ; tÞ ¼ P τ ½qk3 ðx1 ; tÞρτ ∂t Δuk3 þ 1 ðx1 ; tÞ where P τ ½q3  ¼

(

q3 ðx3 ; tÞ

if q3 jr kτ qn ðx3 ; tÞ

kτ qn ðx3 ; tÞq3 =jq3 j

if jq3 j 4kτ qn ðx3 ; tÞ

ð4:29Þ

ð4:30Þ

1508



V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

  is the operator of the orthogonal projection onto set q3  r kτ qn ðx3 ; tÞ, coefficient ρτ has been chosen based on the conditions that give the best convergence of the algorithm; and proceed to the next step of the iteration.

6. Verification and comparative study of the TD and FD approaches The FD approach presented here has been verified and compared with known analytical and numerical results reported in our previous publications that have been quoted in this paper many times. Therefore we will verify and compare the TD approach. First we consider a finite crack subjected by a suddenly applied unit impulse that causes an antiplane deformation as it is shown in Fig. 1. This problem has an analytical solution reported in [39]

and used in numerous publications for the verification numerical BE solution [15,43]. In Figs. 2 and 3 crack openings calculated for dimensionless quantities μ ¼ 1, c2 ¼ 1 and v ¼ 0:25 are presented for the problem considered in [39]. Fig. 2 illustrates the dependence of the crack opening on time. It is observed that the presented solution is in a good agreement with the analytical results of [39] and BE solutions. Indeed maximum of the dynamic crack opening exceed its static value by 30% and for a long time the dynamical solution is equal to the static one that can be calculated using formula qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4p l u3 ðx1 Þ ¼ 3 1ðx1 =lÞ2 ð5:1Þ μ Fig. 3 illustrates the dependence of the dynamic crack opening on a mesh refinement in the spatial discretization; 10, 50 and 100 BE have been used. The diagram on the left corresponds to a long time and a static solution and the diagram on the right

Fig. 4. Frictional contact force and crack opening for k2 ¼ 0:25 and kτ ¼ 0:2.

Contact force q3

0.5

Crack opening Δu3μ/2

2

1. linear 2. constant 3. frequency domain

1

1.5

2 0.25

1 3 0.5

0

0 −0.5

−0.25

−1 −1.5

−0.5

0

T

2T Time

3T

−2

0

T

2T Time

Fig. 5. Frictional contact force and crack opening for k2 ¼ 0:25 and kτ ¼ 0:4.

3T

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

corresponds to time t ¼ 1, when the maximal value of the dynamic crack opening occurs. We have to mention that the calculations presented in Figs. 2 and 3 correspond to the piece-wise linear time approximation. The presented results show that the solution is not mesh-sensitive, and that no essential quantities or qualitative differences were observed for different mesh and piece-wise constant and piece-wise linear time approximations. Thus, the presented results show very good agreement of the solution obtained by the proposed TD BEM with the known analytical and numerical solutions. Therefore we will solve the above formulated elastodynamic problem for the finite crack subjected to harmonic loading (1.13) using the TD BEM byh considering frictional contact conditions (1.15) and compare the result with FD solution. The case of normal incidence of the harmonic shear SH-wave which produces an antiplane deformation is studied here. For simplicity we suppose that the force normal to the crack edges is constant along the crack length for any time and equal to one, i.e.

Contact force q3

1509

qn ¼ 1. We also assume that the incident shear SH-wave has the unit amplitude and the crack is located in the material with the following mechanical characteristics: elastic modulus E ¼ 200 GPa, Poisson's ratio ν ¼ 0:25, specific density ρ ¼ 7800 kg=m3 . The crack opening and frictional contact force in the middle point of the crack (x1 ¼ 0) during three periods 3T of the wave action for different wave numbers and coefficients of friction are presented in Figs. 4–7. Calculations have been done for piece-wise constant and piece-wise linear TD and FD approaches. The diagram on the left shows that the frictional contact forces calculated by the three methods mentioned above coincide for different wave numbers and coefficients of friction, whereas crack opening slightly differ. Figs. 8–11 illustrate the frictional contact force and the crack opening distribution along the crack length during the time period T of the wave load action. It is interesting to notice that in the FD approach iterations are carried out for the whole period of time, whereas in the TD approach for each time interval. In [57] it was shown that usually

Crack opening Δu3μ/2

1.5

1. linear 2. constant 3. frequency domain

1

0.2

1

0.1

0.5 3 2

0

0

−0.1

−0.5

−0.2

−1

0

T

2T

3T

−1.5

0

T

2T

Time

3T

Time

Fig. 6. Frictional contact force and crack opening for k2 ¼ 0:75 and kτ ¼ 0:2.

Contact force q3

0.5

Crack opening Δu3μ/2

2 3

1. linear 2. constant 3. frequency domain

1

1.5 0.25

1 0.5

0

0 2 −0.5

−0.25

−1 −1.5

−0.5

0

T

2T Time

3T

−2

0

T

2T Time

Fig. 7. Frictional contact force and crack opening for k2 ¼ 0:75 and kτ ¼ 0:4.

3T

1510

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

Fig. 8. Frictional contact force for k2 ¼ 0:25 and kτ ¼ 0:2.

Fig. 9. Crack opening for k2 ¼ 0:25 and kτ ¼ 0:2.

3–6 iterations are enough for the FD algorithm convergence. Fig. 12 illustrates the convergence of the TD algorithm at each time step for the piece-wise constant and piece-wise linear approximations respectively. The SIF is the main characteristic of fracture mechanics [4,7]. In our previous publications it was shown that contact interaction affects fraction mechanics criterions. Therefore we will study the influence of the frictional contact interaction of the crack edges on the SIF here. For the tearing (or antiplane) mode III the displacement at the crack tip has the form of rffiffiffiffiffi K III 2ε u3 ¼ sin ðθ=2Þ ð5:2Þ π μ

Using limit transition for θ ¼ π we obtain pffiffiffi μ π K III ¼ limpffiffiffiffiffiu3 ðlε; tÞ ε-0 2ε

ð5:3Þ

The stress intensity factor against wave number for different coefficients of friction is presented in Fig. 13. The curves on the graphs correspond to: 1 – solution without contact and 2 and 3 – with contact for kτ ¼ 0:2 and kτ ¼ 0:4 respectively. We have to mention that the solution obtained without crack edges contact interaction (curve 1) coincides with the one presented in [32]. Our calculations show that for the same mesh and number of time steps the average time of calculation for the piece-wise constant and piece-wise linear time approximation in the TD

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

1511

Fig. 10. Frictional contact force for k2 ¼ 0:25 and kτ ¼ 0:4.

Fig. 11. Crack opening for k2 ¼ 0:25 and kτ ¼ 0:4.

1.5

1

0.5

0

Fig. 12. Number of iterations at each time step.

0

0.25

0.5

0.75

1

1.25

1.5

Fig. 13. Stress intensity factor against wave number.

1.75

2

k2l

1512

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

formulation is the same, while for the FD formulation the average time of calculation is 1.5–2.0 times less. Also the convergence of the algorithm and the accuracy of calculations in the FD formulation are much better. For example convergence of the algorithm in the TD formulation is more sensitive to the mesh refinement and coefficient ρτ in (3.31). But on the other hand the TD formulation is more general and the problem can be solved not only for harmonic or periodic load, but also for a much wider range of loads.

7. Conclusion This paper presents the comparative study of two approaches based on the TD and FD BEM analysis of antiplane crack problem in 2-D, homogeneous, isotropic and linear elastic solids by considering the frictional contact of the crack edges. Hypersingular integrals that appear in the TD and FD BIE formulations, which correspond to the elastodynamic fundamental solutions, are evaluated using special regularization algorithm that are based on the theory of generalized functions. This algorithm converts the hypersingular integrals to regular integrals, which can be computed analytically, hence no numerical calculations are needed. In order to fulfill frictional contact conditions a special iterative algorithm of Udzava's type is used. The algorithm consists of two parts. The first part of the algorithm is a solution of an elastodynamic problem for bodies with cracks, without taking into account unilateral restrictions. The second part of the algorithm is a projection of the founded solution on the set of unilateral restrictions. For the first part of the algorithm the DT and FD BEM are used. A projection operator has been constructed and further modified for better convergence in our previous publications. An analysis of the algorithm convergence for each time interval for the TD solution was done. Numerical results for the crack opening, frictional contact forces and dynamic SIF are presented and discussed for a finite III-mode crack in an infinite domain subjected to a harmonic crack-face loading while considering crack edges frictional contact interaction. It was shown that the FD solution takes 1.5–2.0 times less time for its calculation in comparison with the TD solution. But the TD formulation is more general and the problem can be solved not only for harmonic or periodic load, but also for a much wider range of loads.

Acknowledgments This work was partly supported by the Committee of Science and Technology of Mexico (CONACYT) by the Research grants (Project no. 000000000101415), which is gratefully acknowledged. References [1] Abramowitz M, Stegun IA. Handbook of mathematical functions, with formulas, graphs and mathematical tables. Applied Mathematics Series 1964 pp. 1046. [2] Achenbach JD. Wave propagation in elastic solids. Amsterdam: North-Holland Publishing Co.; 1973. [3] Aliabadi MH. Boundary element formulations in fracture mechanics. Applied Mechanics Reviews 1997;50(2):83–96. [4] Anderson TL. Fracture mechanics: fundamentals and applications. CRC Press LLC; 1995. [5] Antes H, Panagiotopoulos PD. The boundary integral approach to static and dynamic contact problems. Basel, Boston, Berlin: Birkhauser; 1992. [6] Balas J, Sladek J, Sladek V. Stress analysis by boundary element methods. Amsterdam: Elsevier; 1989. [7] Cherepanov GP. Mechanics of brittle fracture. New York: McGraw Hill; 1979. [8] Denda M, Wang C-Y, Yong YK. 2-D time-harmonic BEM for solids of general anisotropy with application to eigenvalue problems. Journal of Sound and Vibration 2003;261:247–76.

[9] Dineva P, Rangelov T, Gross D. BIEM for 2D steady-state problems in cracked anisotropic materials. Engineering Analysis with Boundary Elements 2005;29: 689–98. [10] Dominguez J. Boundary elements in dynamics. Southampton: Computational Mechanics Publications; 1993. [11] Eringen AC, Suhubi ES. Elastodynamics. Linear theory, vol. 2. New York: Academic Press; 1975. [12] Fedelinski P, Aliabadi MH, Rooke DP. A single-region time domain BEM for dynamic crack problems. International Journal of Solids and Structures 1995;32:3555–71. [13] Fedelinski P, Aliabadi MH, Rooke DP. The Laplace transform DBEM for mixedmode dynamic crack analysis. Computer and Structures 1996;59(6):1021–31. [14] Freund LВ. Dynamic fracture mechanics. New York, NY, USA: Cambridge University Press; 1998. [15] Gallego R, Dominguez J. Hypersingular BEM for transient elastodynamics. International Journal for Numerical Methods in Engineering 1996;39: 1681–705. [16] Garcia-Sanchez F, Zhang Ch. A comparative study of three BEM for transient dynamic crack analysis of 2-D anisotropic solids. Computational Mechanics 2007;40:753–69. [17] Garcia-Sanchez F, Saez A, Dominguez J. Traction boundary elements for cracks in anisotropic solids. Engineering Analysis with Boundary Elements 2004;28: 667–76. [18] Garcia-Sanchez F, Saez A, Dominguez J. Two-dimensional time-harmonic BEM for cracked anisotropic solids. Engineering Analysis with Boundary Elements 2006;30:88–99. [19] Garcıa-Sanchez F, Zhang Ch, Sladek J, Sladek V. A 2D time-domain BEM for dynamic crack problems in anisotropic solids. In: Manolis GD, Polyzos D, editors. Recent advances in boundary element methods. Springer; 2009 http: //dx.doi.org/10.1007/978-1-4020-9710-2, ISBN 978-1-4020-9709-6 e-ISBN 978-1-4020-9710-2. [20] Guz AN, Zozulya VV. Brittle fracture of materials under dynamic loading. Kiev: Naukova Dumka; 1993 ([in Russian]). [21] Guz AN, Zozulya VV. Problems of dynamic fracture mechanics with allowance for contact of the crack edges. International Applied Mechanics 1995;31 (1):1–31. [22] Guz AN, Zozulya VV. About taking into account of cracks edges contact under dynamic loads. Materials Science 1996;32(1):29–44. [23] Guz AN, Zozulya VV. Fracture dynamics with allowance for a crack edges contact interaction. International Journal of Nonlinear Sciences and Numerical Simulation 2001;2:173–233. [24] Guz AN, Zozulya VV. Elastodynamic unilateral contact problems with friction for bodies with cracks. International Applied Mechanics 2002;38 (8):895–932. [25] Guz AN, Zozulya VV. Investigation of the effect of frictional contact in III-mode crack under action of the SH-wave harmonic load. CMES: Computer Modeling in Engineering & Science 2007;22(2):119–28. [26] Guz AN, Zozulya VV. Contact problem for the flat crack under two normally incident shear H-waves with wave mode-shifting. Theoretical and Applied Fracture Mechanics 2010;54(3):189–95. [27] Guz AN, Zozulya VV. Contact problem for the flat crack under two normally incident tension-compression waves with wave mode-shifting. Engineering Analysis with Boundary Elements 2011;35(1):34–41. [28] Guz AN, Menshykov OV, Zozulya VV, Guz IA. Contact problem for the flat elliptical crack under normally incident shear wave. CMES: Computer Modeling in Engineering & Science 2007;17(3):205–14. [29] Israil ASM, Banerjee PK. Advanced time-domain formulation of BEM for twodimensional transient elastodynamics. International Journal for Numerical Methods in Engineering 1990;29:1421–40. [30] Israil ASM, Banerjee PK. Two-dimensional transient wave-propagation problems by time domain BEM. International Journal of Solids and Structures 1990;26(8):851–64. [31] Israil ASM, Banerjee PK. Interior stress calculations in 2-D time-domain transient BEM analysis. International Journal of Solids and Structures 1990;27(7): 915–27. [32] Loeber JF, Sih GC. Diffraction of antiplane shear waves by a finite crack. Journal of the Acoustical Society of America 1968;44:90–8. [33] Man KW. Contact mechanics using boundary elements. Southampton: Computational Mechanics Publications; 1994. [34] Manolis GA. Comparative study on three boundary element method approaches to problems in elastodynamics. International Journal for Numerical Methods in Engineering 1983;19:73–91. [35] Mettu SR, Nicholson JW. Computation of dynamic stress intensity factors by the time domain boundary integral equation method. II. Examples. Engineering Fracture Mechanics 1988;31(5):769–82. [36] Nicholson JW, Mettu SR. Computation of dynamic stress intensity factors by the time domain boundary integral equation method. I. Analysis. Engineering Fracture Mechanics 1988;31(5):759–67. [37] Panagiotopoulos PD. Inequality problems in mechanics and applications. Convex and non convex energy functions. Stuttgart: Birkhauser; 1985. [38] Sih GC, editor. Elastodynamic crack problems. Noordhoff: Leyden; 1977. [39] Thau SA, Lu T-H. Diffraction of transient horizontal shear waves by a finite crack and a finite rigid ribbon. International Journal Engineering Science 1970;8:857–74. [40] Zhand Ch, Gross D. On wave propagation in elastic solids with cracks. Southampton: Computational Mechanics Publications; 1997.

V.V. Zozulya / Engineering Analysis with Boundary Elements 37 (2013) 1499–1513

[41] Zhang Ch. Transient elastodynamic antiplane crack analysis of anisotropic solids. International Journal of Solids and Structures 2000;37:6107–30. [42] Zhang Ch. A 2-D time-domain BIEM for dynamic analysis of cracked orthotropic solids. CMES: Computer Modeling in Engineering & Sciences 2002;3: 381–98. [43] Zhang Ch. A 2D hypersingular time-domain traction BEM for transient elastodynamic crack analysis. Wave Motion 2002;35:17–40. [44] Zozulya VV. Dynamic problems of crack theory with regions of contact, adhesion and slip. Dokladi Akademii Nauk Ukrainskoy SSR, Seria A, Fisika, Mathamika i Techicheckie Nauki 1990;1:47–50 ([in Russian]). [45] Zozulya VV. Solvability of dynamic problems of crack theory in regions of contact, adhesion and slip. Dokladi Akademii Nauk Ukrainskoy SSR, Seria A, Fisika, Mathamika i Techicheckie Nauki 1990;3:53–5 ([in Russian]). [46] Zozulya VV. Action of a harmonic load on a crack in an infinite body with allowance for interaction of its edges. Dokladi Akademii Nauk Ukrainskoy SSR, Seria A, Fisika, Mathamika i Techicheckie Nauki 1990;4:46–9 ([in Russian]). [47] Zozulya VV. Hadamard integrals in dynamic problems of crack theory. Dokladi Akademii Nauk Ukrainskoy SSR, Seria A, Fisika, Mathamika i Techicheckie Nauki 1991;2:19–22 ([in Russian]). [48] Zozulya VV. Contact interaction between the edges of a crack and an infinite plane under a harmonic load. International Applied Mechanics 1992;28(1): 61–5. [49] Zozulya VV. Boundary functionals method in contact problems of the dynamics of bodies with cracks. Doklady Akademii Nauk Ukraine 1992;2:38–44 ([in Russian]). [50] Zozulya VV. Solution of a problems of the dynamics of bodies with cracks by the method of boundary integral equations. Doklady Akademii Nauk Ukraine 1992;3:38–43 ([in Russian]). [51] Zozulya VV. Regularization of the divergent integrals. I. General consideration. Electronic Journal of Boundary Elements 2006;4:49–57.

1513

[52] Zozulya VV. Regularization of the divergent integrals. II. Application in fracture mechanics. Electronic Journal of Boundary Elements 2006;4(2):58–66. [53] Zozulya VV. The regularization of the divergent integrals in 2-D elastostatics. Electronic Journal of Boundary Elements 2009;7(2):50–88. [54] Zozulya VV. Variational formulation and nonsmooth optimization algorithms in elastostatic contact problems for cracked body. CMES: Computer Modeling in Engineering & Sciences 2009;42(3):187–215. [55] Zozulya VV. Regularization of hypersingular integrals in 3-D fracture mechanics: triangular BE, and piecewise-constant and piecewise-linear approximations. Engineering Analysis with Boundary Elements 2010;34(2):105–13. [56] Zozulya VV. Divergent integrals in elastostatics: regularization in 3-D case. CMES: Computer Modeling in Engineering & Sciences 2010;70(3):253–350. [57] Zozulya VV. Variational formulation and nonsmooth optimization algorithms in elastodynamic contact problems for cracked body. Computer Methods in Applied Mechanics and Engineering 2011;200(5–8):525–39. [58] Zozulya VV, Gonzalez-Chi PI. Dynamic fracture mechanics with crack edges contact interaction. Engineering Analysis with Boundary Elements 2000;24 (9):643–59. [59] Zozulya VV, Menshykov OV. Use of the constrained optimization algorithms in some problems of fracture mechanics. Optimization and Engineering, 4; 365–84. [60] Zozulya VV, Men'shikov VA. Solution of tree-dimensional problems of the dynamic theory of elasticity for bodies with cracks using hypersingular integrals. International Applied Mechanics 2000;36(1):74–81. [61] Zozulya VV. Stress intensity factor in contact problem for the flat crack under antiplane shear wave. International Applied Mechanics 2007;43(9):1043–7. [62] Zozulya VV. Frequency and time domain BEM in fracture dynamic contact problems. Comparative study. In: Eberhardsteiner J. et. al., editors. European Congress on Computational Methods in Applied Sciences and Engineering, (ECCOMAS 2012), Vienna, Austria; 2012; 20 p.