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

V.V. Zozulya

► To cite this version:

V.V. Zozulya. Comparative study of time and frequency domain BEM approaches in frictional contact problem for antiplane crack under harmonic loading. Engineering Analysis with Boundary Elements, 2013, 37 (11), pp.1499 - 1513. 10.1016/j.enganabound.2013.08.006. hal-01382199

HAL Id: hal-01382199

https://hal.science/hal-01382199

Submitted on 16 Oct 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comparative study of time and frequency domain BEM approaches in frictional contact problem for antiplane crack under harmonic loading

V.V. Zozulya*

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

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.

## 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.

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 subdifferential functionals and finding of their saddle points. Such algorithms are usually known as Uzava's type [5,37].
Mathematical aspects of the problem and algorithm convergence were investigated in [48,50,57]. It has been shown that the algorithm may be considered as a compressive operator, acting on a special Sobolev's functional space. Briefly speaking the algo rithm 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 con structed 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 con ditions 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 diver gent 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 frequency domain [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.

## 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_α = (x_1, x_2) \in V \subset R^2$ and time $t \in \mathfrak{F}$, 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 $\sigma_{3\alpha}(x_\alpha, t)$ and strain $\varepsilon_{3\alpha}(x_\alpha, t)$ tensors, which are related by Hooke's law

$$ \sigma_{3\alpha} = \mu\varepsilon_{3\alpha} \qquad (1.1) $$

where $\mu$ is the shear modulus.

The deformation is described by the shear component of the displacements $u_3(x_\alpha, t)$, which is related to the strain tensor by the

<figure><img src="image_2.png" /><figcaption>Fig. 1. Finite crack under antiplane deformation.</figcaption></figure>

kinematic Cauchy relation

$$ \varepsilon_{3\alpha} = \partial_{\alpha} u_3 \qquad (1.2) $$

Substituting (1.2) into (1.1) the stress tensor components are defined by the displacement

$$ \sigma_{3\alpha} = \mu \partial_{\alpha} u_3 \qquad (1.3) $$

The stress differential equation of motion of the elastic body in this case has the form of

$$ \mu \partial_{\beta} \sigma_{3\beta} + b_3 = \rho \partial_t^2 u_3 \quad \forall (x_\alpha, t) \in V \times \mathfrak{F} \qquad (1.4) $$

Here and above $\partial_\beta = \partial/\partial x_\beta$ and $\partial_t = \partial/\partial t$ are derivatives with respect to the space coordinates and time, respectively, $b_3$ is the volume force.

Substituting $\sigma_{3\alpha}(x_\alpha, t)$ in Eq. (1.4) by its Hooke's law representa tion (1.3) we find the scalar wave equation for the displacement $u_3(x_\alpha, t)$ to be in the form of

$$ \mu \partial_{\beta} \partial_{\beta} u_3 + b_3 = \rho \partial_t^2 u_3 \quad \forall (x_\alpha, t) \in V \times \mathfrak{F} \qquad (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

$$ u_3(x_\alpha, t_0) = u_3^0(x), \quad \partial_t u_3(x_\alpha, t_0) = v_3^0(x_\alpha) \quad \forall x_\alpha \in V \qquad (1.6) $$

Additional conditions at infinity have to be satisfied in this case

$$ u_3(x_\alpha, t) = O(\ln(r^{-1})), \quad \sigma_{3\beta}(x_\alpha, t) = O(r^{-1}) \quad \text{for } r \to \infty \qquad (1.7) $$

Here $r = |x_\alpha| = \sqrt{x_1^2 + x_2^2}$ 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 $\mathbb{R}^2$ with a $C^{1,1}$ Lyapunov's class regular boundary $\partial V$. The boundary contains two parts $\partial V_u$ and $\partial V_p$ such that $\partial V_u \cap \partial V_p = \emptyset$ and $\partial V_u \cup \partial V_p = \partial V$. On the part $\partial V_u$ are prescribed displacements $u_3(x_\alpha, t)$ of the body points and on the part $\partial V_p$ are prescribed tractions $p_3(x_\alpha, t)$ respec tively. Then the mixed boundary conditions are

$$ \begin{aligned} u_3(x_\alpha, t) &= \phi_3(x_\alpha, t) \quad \forall x_\alpha \in \partial V_u \quad \forall t \in \mathfrak{F} \\ p_3(x_\alpha, t) &= \sigma_{3\beta}(x_\alpha, t) n_\beta(x_\alpha) = P_n[u_3(x_\alpha, t)] \\ &= \psi_3(x_\alpha, t) \quad \forall x_\alpha \in \partial V_p \quad \forall t \in \mathfrak{F} \end{aligned} \qquad (1.8) $$
equations:

$$U_3(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{1}{2\pi\mu L} H(t' \ r/c_2) \qquad (2.5)$$

where $H(t' \ r/c_2)$ is the Heaviside step function, $r(\mathbf{x}_a, \mathbf{y}_a) = |\mathbf{x}_a \ \mathbf{y}_a|$ is the distance between the collocation point $\mathbf{y}_a$ and the observation point $\mathbf{x}_a$. It is defined by the equation

$$r(\mathbf{x}_a, \mathbf{y}_a) = \sqrt{(x_1 \ y_1)^2 + (x_2 \ y_2)^2} \qquad (2.6)$$

For compactness we also introduce the following notations
$t' = t \ \tau, L = \sqrt{(t')^2 \ (r/c_2)^2}$.

Applying traction operator (1.9) to the fundamental solution for displacements (2.5) we obtain fundamental solution for traction in the form of

$$W_3(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{1}{2\pi c_2} \frac{n_\alpha(\mathbf{x}_a)\mathbf{x}_a}{r} \left( \frac{r H(t' \ r/c_2)}{c_2 L^3} \frac{\delta(t' \ r/c_2)}{L} \right) \qquad (2.7)$$

Integral representation for traction can be obtained from Somigliano's integral representation (2.4) applying differential operator (1.9). As a result we have

$$p_3(\mathbf{y}_a, t) = \int_{\mathfrak{S}, \partial V} [p_3(\mathbf{x}_a, \tau) K_3(\mathbf{x}_a, \mathbf{y}_a, t')] u_3(\mathbf{x}_a, \tau) F_3(\mathbf{x}_a, \mathbf{y}_a, t')] dS \ d\tau \\ \int_{\mathfrak{S}, \Omega} \Delta u_3(\mathbf{x}_a, \tau) F_3(\mathbf{x}_a, \mathbf{y}_a, t')] dS \ d\tau \qquad (2.8)$$

Here the kernel $K_3(\mathbf{x}_a, \mathbf{y}_a, t')$ has the form of

$$K_3(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{1}{2\pi c_2} \frac{n_\beta(\mathbf{y}_a)\mathbf{x}_\beta}{r} \left( \frac{r H(t' \ r/c_2)}{c_2 L^3} \frac{\delta(t' \ r/c_2)}{L} \right) \qquad (2.9)$$

The kernel $F_3(\mathbf{x}_a, \mathbf{y}_a, t')$ can be calculated applying the operator of the derivation in normal direction twice with respect to $\mathbf{x}_a$ and $\mathbf{y}_a$ to the fundamental solutions (2.5). As the result we obtain

$$F_3(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{n_\alpha n_\beta}{r} \left( \delta_{a\beta} \frac{\mathbf{x}_a \mathbf{x}_\beta}{r^2} \right) \frac{dU_3(\mathbf{x}_a, \mathbf{y}_a)}{dr} + \frac{n_\alpha \mathbf{x}_a}{r} \frac{n_\beta \mathbf{x}_\beta}{r} \frac{d^2 U_3(\mathbf{x}_a, \mathbf{y}_a)}{dr^2} \qquad (2.10)$$

where

$$\frac{dU_3(\mathbf{x}_a, \mathbf{y}_a, t')}{dr} = \frac{r H(t' \ r/c_2)}{c_2^2 L^3} \frac{\delta(t' \ r/c_2)}{c_2 L} \\ \frac{d^2 U_3(\mathbf{x}_a, \mathbf{y}_a, t')}{dr^2} = \frac{2r^2/c_2^2 + (t \ \tau)^2}{c_2^2 L^5} \frac{H(t \ \tau \ r/c_2)}{L} \\ \frac{2r \delta(t \ \tau \ r/c_2)}{c_2^3 L^3} + \frac{\delta'(t \ \tau \ r/c_2)}{c_2^2 L} \qquad (2.11)$$

In Eq. (2.10) the expressions for normal derivatives of the function $r(\mathbf{x}_a, \mathbf{y}_a)$ have the form of

$$\partial_n r = \frac{n_\alpha \mathbf{x}_\alpha}{r}, \quad \partial_n(\partial_n r) = \frac{n_\alpha n_\beta}{r} \left( \delta_{a\beta} \frac{\mathbf{x}_a \mathbf{x}_\beta}{r^2} \right) \qquad (2.12)$$

The term with delta function and its derivatives in (2.7), (2.9) and (2.11) possess singularity at $t' = r/c_2$, 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

$$\int_{0}^{t} \delta(t \ \tau \ r/c_2) f(\tau) u(\tau) d\tau = \int_{0}^{t} H(t \ \tau \ r/c_2) (f(\tau) u(\tau) + f(\tau) \dot{u}(\tau)) d\tau \qquad (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

$$u_3(\mathbf{y}_a, t) = \int_{\mathfrak{S}, \partial V} [p_3(\mathbf{x}_a, \tau) U_3(\mathbf{x}_a, \mathbf{y}_a, t') \\ u_3(\mathbf{x}_a, \tau) W_3^*(\mathbf{x}_a, \mathbf{y}_a, t')] \dot{u}_3(\mathbf{x}_a, \tau) \dot{W}_3^*(\mathbf{x}_a, \mathbf{y}_a, t')] dS \ d\tau$$

$$\int_{\mathfrak{S}, \partial V} [\Delta u_3(\mathbf{x}_a, \tau) W_3^*(\mathbf{x}_a, \mathbf{y}_a, t') + \Delta \dot{u}_3(\mathbf{x}_a, \tau) \dot{W}_3^*(\mathbf{x}_a, \mathbf{y}_a, t')] dS \ d\tau \\ \int_{\mathfrak{S}, \Omega} [\Delta u_3(\mathbf{x}_a, \tau) W_3^*(\mathbf{x}_a, \mathbf{y}_a, t') + \Delta \dot{u}_3(\mathbf{x}_a, \tau) \dot{W}_3^*(\mathbf{x}_a, \mathbf{y}_a, t')] dS \ d\tau \qquad (2.14)$$

where

$$W_3^*(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{1}{2\pi c_2} \frac{n_\alpha(\mathbf{x}_a)\mathbf{x}_a}{r} \frac{t'}{L^3} \frac{r}{L} H(t' \ r/c_2) \\ \dot{W}_3^*(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{1}{2\pi c_2} \frac{n_\alpha(\mathbf{x}_a)\mathbf{x}_a}{r} \frac{1}{L} H(t' \ r/c_2) \qquad (2.15)$$

In the same way the integral representation for traction (2.8) can be represented in the form of

$$p_3(\mathbf{y}_a, t) = \int_{\mathfrak{S}, \partial V} [p_3(\mathbf{x}_a, \tau) K_3^*(\mathbf{x}_a, \mathbf{y}_a, t') + \dot{p}_3(\mathbf{x}_a, \tau) \dot{K}_3^*(\mathbf{x}_a, \mathbf{y}_a, t')] \\ u_3(\mathbf{x}_a, \tau) F_3^*(\mathbf{x}_a, \mathbf{y}_a, t') \dot{u}_3(\mathbf{x}_a, \tau) \dot{F}_3^*(\mathbf{x}_a, \mathbf{y}_a, t') \\ u_3(\mathbf{x}_a, \tau) F_3^*(\mathbf{x}_a, \mathbf{y}_a, t')] dS \ d\tau \\ \int_{\mathfrak{S}, \Omega} [\Delta u_3(\mathbf{x}_a, \tau) F_3^*(\mathbf{x}_a, \mathbf{y}_a, t') \Delta \dot{u}_3(\mathbf{x}_a, \tau) \dot{F}_3^*(\mathbf{x}_a, \mathbf{y}_a, t')] dS \ d\tau \qquad (2.16)$$

where

$$K_3^*(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{\mu}{2\pi c_2} \frac{n_\beta(\mathbf{y}_a)\mathbf{x}_\beta}{r} \frac{t'}{L^3} \frac{r}{L} H(t' \ r/c_2) \\ \dot{K}_3^*(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{\mu}{2\pi c_2} \frac{n_\beta(\mathbf{y}_a)\mathbf{x}_\beta}{r} \frac{1}{L} H(t' \ r/c_2) \qquad (2.17)$$

and

$$F_3^*(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{\mu}{2\pi c_2} \frac{(t' \ r/c_2)}{rL^3} \left( 1 - \frac{n_\alpha(\mathbf{x}_a)\mathbf{x}_a}{r} \frac{n_\beta(\mathbf{y}_a)\mathbf{x}_\beta}{r} \right) \\ \times \left( 1 + \frac{3r}{c_2t' + r} \right) H(t' \ r/c_2), \\ \dot{F}_3^*(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{\mu}{2\pi c_2} \frac{1}{rL} \left( 1 - \frac{n_\alpha(\mathbf{x}_a)\mathbf{x}_a}{r} \frac{n_\beta(\mathbf{y}_a)\mathbf{x}_\beta}{r} \right) \\ \times \left( 1 + \frac{2r}{c_2t' + r} \right) H(t' \ r/c_2), \\ F_3^*(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{\mu}{2\pi c_2} \frac{n_\alpha(\mathbf{x}_a)\mathbf{x}_a}{r} \frac{n_\beta(\mathbf{y}_a)\mathbf{x}_\beta}{r} \frac{1}{c_2L} \qquad (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 $n_1 = 0, n_2 = 1, x_2 = 0$ and $\partial V = \emptyset$. Therefore the integral equation that related the crack opening and the traction is presented in the form of

$$p_3(\mathbf{y}_a, t) = \int_{\mathfrak{S}, \partial V} [\Delta u_3(\mathbf{x}_a, \tau) F_3^*(\mathbf{x}_a, \mathbf{y}_a, t') \Delta \dot{u}_3(\mathbf{x}_a, \tau) \dot{F}_3^*(\mathbf{x}_a, \mathbf{y}_a, t')] dS \ d\tau \qquad (2.19)$$

with

$$F_3^*(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{\mu}{2\pi c_2 r} \frac{(t' \ r/c_2)}{L^3}, \quad \dot{F}_3^*(\mathbf{x}_a, \mathbf{y}_a, t') = \frac{\mu}{2\pi c_2 r} \frac{1}{L} \qquad (2.20)$$

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

## 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 $\Omega$ which is divided into $N$ boundary elements with $Q$ nodes of interpolation per element.
Standard Lagrange interpolation polynomials have the form of

$$ \phi_q(\xi) = \frac{(\xi \xi_0)(\xi \xi_1)\dots(\xi \xi_{q-1})(\xi \xi_{q+1})\dots(\xi \xi_Q)}{(\xi_q \xi_0)(\xi_q \xi_1)\dots(\xi_q \xi_{q-1})(\xi_q \xi_{q+1})\dots(\xi_q \xi_Q)} \quad (3.1) $$

For piece wise constant approximation

$$ \phi_q(\mathbf{x}_\alpha) = \begin{cases} 1 & \forall \mathbf{x}_\alpha \in \Omega_n, \\ 0 & \forall \mathbf{x}_\alpha \notin \Omega_n. \end{cases} \quad (3.2) $$

The observation time interval $\mathfrak{I} = [t_0, T]$ is divided into $M$ time steps with $P$ nodes of interpolation per time step

For piece wise constant approximation

$$ \psi_0(t) = H(t \ t_m + \Delta t) \ H(t \ t_m) \quad (3.3) $$

where $H(t)$ represents the Heaviside unit step function.

For piece wise linear approximation

$$ \psi_1(t) = \frac{1}{\Delta t}(R(t \ t_m + \Delta t) \ 2R(t \ t_m) + R(t \ t_m \ \Delta t)) \quad (3.4) $$

where $R(t)$ represents the unit ramp function, defined as $R(t) = t$ if $t > 0$ and $0$ otherwise.

Now the crack opening, its velocity and traction on the boundary element $\Omega_n$ and the time step $[t_{m-1}, t_m]$ can be presented in the form of

$$ \Delta u_3(\mathbf{x}_\alpha, \tau) = \sum_{p}^{P} \sum_{1 q}^{Q} \Delta u_3(\mathbf{x}_\alpha^q, t^m) \phi_q(\mathbf{x}) \psi_p(\tau), \quad \mathbf{x}_\alpha \in \Omega_n, \quad \tau \in [t_{m-1}, t_m], $$

$$ \Delta \dot{u}_3(\mathbf{x}_\alpha, \tau) = \sum_{p}^{P} \sum_{1 q}^{Q} \Delta u_3(\mathbf{x}_\alpha^q, t^m) \phi_q(\mathbf{x}) \dot{\psi}_p(\tau), \quad \mathbf{x}_\alpha \in \Omega_n, \quad \tau \in [t_{m-1}, t_m], $$

$$ p_3(\mathbf{x}_\alpha, \tau) = \sum_{p}^{P} \sum_{1 q}^{Q} p_3(\mathbf{x}_\alpha^q, t^m) \phi_q(\mathbf{x}) \psi_p(\tau), \quad \mathbf{x}_\alpha \in \Omega_n, \quad \tau \in [t_{m-1}, t_m] \quad (3.5) $$

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

$$ \begin{aligned} p_3(\mathbf{y}_r, t_M) &= \sum_{n}^{Q} \sum_{1 q}^{Q} \sum_{1 m}^{M} \sum_{0 p}^{P} \left( \begin{aligned} & F_3^{nq}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) \Delta u_3(\mathbf{x}_n, t_m) \\ & \dot{F}_3^{nq}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) \Delta \dot{u}_3(\mathbf{x}_n, t_m) \end{aligned} \right) \\ & \qquad \left( \begin{aligned} & F_3^{nq}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) \Delta u_3(\mathbf{x}_n, t_m) \\ & \dot{F}_3^{nq}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) \Delta \dot{u}_3(\mathbf{x}_n, t_m) \end{aligned} \right) \end{aligned} \quad (3.6) $$

where

$$ F_3^{q,p}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) = \int_{t_{m-1}}^{t_m} \int_{\Omega_n} F_3^{nq}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) \tau \phi_q(\xi) dS \psi_p(\tau) d\tau, $$

$$ \dot{F}_3^{q,p}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) = \int_{t_{m-1}}^{t_m} \int_{\Omega_n} \dot{F}_3^{nq}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) \tau \phi_q(\xi) dS \dot{\psi}_p(\tau) d\tau, \quad (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

$$ F_3^{q,p}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) = \int_{\Omega_n} \int_{t_{m-1}}^{t_m} F_3^{nq}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) \tau \psi_p(\tau) d\tau \phi_q(\xi) dS $$

$$ \dot{F}_3^{q,p}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) = \int_{\Omega_n} \int_{t_{m-1}}^{t_m} \dot{F}_3^{nq}(\mathbf{y}_r, \mathbf{x}_n, t_M \quad t_m) \tau \dot{\psi}_p(\tau) d\tau \phi_q(\xi) dS \quad (3.8) $$

In the case of piece wise constant time approximation inter polation polynomial and its derivative with respect to time are

$$ \psi_0(t) = H(t \ t_{m+1}) \ H(t \ t_m), \quad \dot{\psi}_0(t) = \delta(t \ t_{m+1}) \ \delta(t \ t_m) \quad (3.9) $$

Substituting them and kernels (2.20) in (3.8) we obtain integral that has to be evaluated in the form of

$$ F_3^{M,m}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, t_M) = \frac{\mu}{2\pi c_2 r} \int_{t_{m-1}}^{t_m} \left( \frac{(t_M \tau) r/c_2}{((t_M \tau)^2 (r/c_2)^2)^{3/2}} + \frac{\delta(\tau \ t_{m+1}) \delta(\tau \ t_m)}{\sqrt{(t_M \tau)^2 (r/c_2)^2}} \right) H(t_M \tau \ r/c_2) d\tau \quad (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_{m-1} \le r/c_2$, the effect of perturbation has not arrived yet at the point $\mathbf{x}_\alpha$

$$ F_3^{M,m}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, t^M) = 0 \quad (3.11) $$

2. For $t_M \ t_m \le r/c_2 < t_M \ t_{m-1}$, the effect of the part of the perturbation has arrived at $\mathbf{x}_\alpha$ but the part of the perturbation produced for $r$ such that $r/c_2 > t_m \ \tau$ has not arrived at $\mathbf{x}_\alpha$ yet

$$ F_3^{M,m}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, t_M) = \frac{\mu}{2\pi r^2} \frac{(t_M \ t_{m-1})}{L_{m-1}} \quad (3.12) $$

Here

$$ L_m = \sqrt{(t_M \ t_m)^2 (r/c_2)^2} $$

3. For $t_n \ t_m < r/c_2$, the effect of part of the perturbation has arrived at $\mathbf{x}_\alpha$ but the part of the perturbation produced for $r$ such that $r/c_2 > t_M \ \tau$ has not arrived at $\mathbf{x}_\alpha$ yet

$$ F_3^{M,m}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, t_M) = \frac{\mu}{2\pi r^2} \left( \frac{(t_M \ t_{m-1})}{L_{m-1}} \frac{(t_M \ t_m)}{L_m} \right) \quad (3.13) $$

Now we proceed with space integration. Considering (3.12) and (3.13) we have to evaluate the following integral:

$$ \int_{-\Delta_n}^{\Delta_n} F_3^{M,m}(x_1, y_1^r, t_M) dx_1 = \frac{\mu}{2\pi} \int_{-\Delta_n}^{\Delta_n} \frac{1}{(x_1 y_1^r)^2} \frac{c_2(t_M t_{m-1})}{L^m} dx_1 \quad (3.14) $$

In the case if $r \ne n$ integral is regular and can be evaluated numerically or analytically without any problems. Analytical evaluation is given as

$$ F_3^{M,m}(r^{n,r}, t_M) = \frac{\mu}{\pi (r^{n,r})^2 + \Delta_n^2} \left( \frac{L_{m-1}^+}{t_M} \frac{L_{m-1}^-}{t_{m-1}} - \frac{L_{m-1}^-}{t_M} \frac{L_{m-1}^+}{t_{m-1}} \right) \quad (3.15) $$

Here

$$ r^{n,r} = |x_1^n - y_1^n|, $$

and

$$ L_m^{\pm} = \sqrt{(t_M \ t_{m-1})^2 ((r^{n,r} \pm \Delta_n)/c_2)^2} $$

In the case where $r = n \ x_1 \in [y_1^n, \Delta_n, y_1^n + \Delta_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

$$ \begin{aligned} F.P. \int_{-\Delta_n}^{\Delta_n} F_3^{M,m}(x_1, y_1^r, t_M) dx_1 \\ = \frac{\mu}{2\pi} F.P. \int_{-\Delta_n}^{\Delta_n} \frac{\phi(x_1)}{(x_1 y_1^r)^2} dx_1, \quad x_1 \in [\Delta_n, \Delta_n] \end{aligned} \quad (3.16) $$
where

$$ \phi(x_1) = \frac{c_2(t_M \quad t_{m-1})}{L^m} $$

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

$$ F_3^{M,m}(r^{n,n}, t_M) = \frac{2\mu}{\pi \Delta_n} \sqrt{1} \left( \frac{\Delta_n}{c_2(t_M \quad t_{m-1})} \right)^2 \quad (3.17) $$

In the case of piece wise linear time approximation interpolation polynomial and its derivative with respect to time are

$$ \begin{aligned} \psi_1(t) &= \frac{1}{\Delta t}(R(t \quad t_m + \Delta t) \quad 2R(t \quad t_m) + R(t \quad t_m \quad \Delta t)), \\ \psi_1(t) &= \frac{1}{\Delta t}(H(t \quad t_m + \Delta t) \quad 2H(t \quad t_m) + H(t \quad t_m \quad \Delta t)) \end{aligned} \quad (3.18) $$

Substituting them and kernels (2.20) in (3.8) we obtain an integral that has to be evaluated in the form of

$$ \begin{aligned} F_3^{M,m}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, t_M) = & \frac{\mu}{2\pi c_2 r} \int_{t_m}^{t_m} \frac{(\tau \quad t_{m-1})}{\Delta t} \\ & \times \left( \frac{(t_M \quad \tau) \quad r/c_2}{((t_M \quad \tau)^2 \quad (r/c_2)^2)^{3/2}} + \frac{\delta(\tau \quad t_{m+1}) \quad \delta(\tau \quad t_m)}{\sqrt{(t_M \quad \tau)^2 \quad (r/c_2)^2}} \right) \\ & \times H(t_M \quad \tau \quad r/c_2) d\tau + \frac{\mu}{2\pi c_2 r} \int_{t_m}^{t_m+1} \frac{(t_{m+1} \quad \tau)}{\Delta t} \\ & \times \left( \frac{(t_M \quad \tau) \quad r/c_2}{((t_M \quad \tau)^2 \quad (r/c_2)^2)^{3/2}} + \frac{\delta(\tau \quad t_m + \Delta t) \quad \delta(\tau \quad t_m)}{\sqrt{(t_M \quad \tau)^2 \quad (r/c_2)^2}} \right) \\ & \times H(t_M \quad \tau \quad r/c_2) d\tau \end{aligned} \quad (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 \quad t_{m-1} \le r/c_2, \Delta t(n \quad m+1) \le r/c_2$, ($t_n \quad t_{m-1} = \Delta t(n \quad m+1)$) the effect of the perturbation has not yet arrived at point $\mathbf{x}_\alpha$

$$ F_3^n(\mathbf{x}_\alpha, \mathbf{y}_\alpha, t^n) = 0 \quad (3.20) $$

2. For $t_n \quad t_m \le r/c_2 < t_n \quad t_{m-1}$, the effect of the part of the perturbation has arrived at $\mathbf{x}_\alpha$ but the part of the perturbation produced for $r$ such as $r/c_2 > t_n \quad \tau$ has not arrived at $\mathbf{x}_\alpha$ yet

$$ F_3^{M,m}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, t_M) = \frac{\mu}{2\pi r^2 \Delta t} L_{m-1} \quad (3.21) $$

3. For $t_n \quad t_{m+1} < r/c_2 \le t_n \quad t_m$, the effect of part of the perturbation has arrived at $\mathbf{x}_\alpha$ but the part of the perturbation produced for $r$ such as $r/c_2 > t_n \quad \tau$ has not arrived at $\mathbf{x}_\alpha$ yet

$$ F_3^{M,m}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, t_n) = \frac{\mu}{2\pi r^2 \Delta t} (L_{m-1} \quad 2L_m) \quad (3.22) $$

4. For $r/c_2 \le t_n \quad t_{m+1}$, the effect of part of the perturbation has arrived at $\mathbf{x}_\alpha$ but the part of the perturbation produced for $r$ such as $r/c_2 > t_n \quad \tau$ has not arrived at $\mathbf{x}_\alpha$ yet

$$ F_3^{M,m}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, t_n) = \frac{\mu}{2\pi r^2 \Delta t} (L_{m-1} \quad 2L_m + L_{m+1}) \quad (3.23) $$

Now we proceed with space integration. When considering (3.21) (3.23) we have to evaluate the following integral:

$$ \int_{-\Delta_n}^{\Delta_n} F_3^{M,m}(x_1, y_1^r, t_M) dx_1 $$

$$ = \frac{\mu}{2\pi c_2 \Delta t} \int_{-\Delta_n}^{\Delta_n} \frac{1}{(x_1 \quad y_1^r)^2} \sqrt{c_2^2(t_M \quad t_{m-1})^2} (x_1 \quad y_1^r)^2 dx_1 \quad (3.24) $$

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

$$ F_3^{M,m}(r^{n,r}, t_M) = \frac{\mu}{2\pi \Delta t} \left( \frac{L_{m-1}^{-r}}{r^{n,r} \Delta_n} + \arctan \left( \frac{r^{n,r} \Delta_n}{c_2 L_{m-1}^{r}} \right) \right) \frac{L_{m-1}^{+r}}{r^{n,r} + \Delta_n} \arctan \left( \frac{r^{n,r} + \Delta_n}{c_2 L_{m-1}^{+r}} \right) \quad (3.25) $$

In the case that $r=n \quad x_1 \in [y_1^r \Delta_n, y_1^r + \Delta_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

$$ F_3^{M,m}(r^{n,n}, t_M) = \frac{\mu}{2\pi \Delta_n} \frac{(t_M \quad t_{m-1})}{\Delta t} \sqrt{1} \frac{\Delta_n}{c_2(t_M \quad t_{m-1})} \quad (3.26) $$

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]

$$ \mathbf{p}^M = \sum_{m=0}^{M} \mathbf{A}^{M,m} \cdot \Delta \mathbf{u}^m \quad (3.27) $$

where Mis a number of time steps and N is a number of BE.

$$ \begin{gathered} \mathbf{p}^M = \begin{vmatrix} p_3(\mathbf{y}_1, t_M) \\ \vdots \\ p_3(\mathbf{y}_N, t_M) \end{vmatrix}, \\ \mathbf{A}^{M,m} = \begin{vmatrix} F_3^{M,m}(\mathbf{y}_1, \mathbf{x}_1, t_M) & \dots & F_3^{M,m}(\mathbf{y}_1, \mathbf{x}_N, t_M) \\ \vdots & \ddots & \vdots \\ F_3^{M,m}(\mathbf{y}_N, \mathbf{x}_1, t_M) & \dots & F_3^{M,m}(\mathbf{y}_N, \mathbf{x}_N, t_M) \end{vmatrix}, \\ \Delta \mathbf{u}^m = \begin{vmatrix} \Delta u_3(\mathbf{y}_1, t_M) \\ \vdots \\ \Delta u_3(\mathbf{y}_N, t_M) \end{vmatrix} \end{gathered} \quad (3.28) $$

Now, taking into account the initial conditions from Eq. (3.27) the explicit time stepping algorithm follows

$$ \Delta \mathbf{u}^M = (\mathbf{A}^0)^{-1} \left( \mathbf{p}^M \sum_{m=1}^{M-1} \mathbf{A}^{M-m+1} \cdot \Delta \mathbf{u}^m \right), \quad (3.29) $$

here $(\mathbf{A}^0)^{-1}$ is the inverse of the system matrix $\mathbf{A}^0$ 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 $q_3^0(x_1, t)$, $\forall \mathbf{x} \in \Omega, \forall t \in \mathfrak{F}$ is assigned;

* the problem without constrains is solved and the unknown quantities on the contact surfaces $\Delta u_3(x_1, t)$ are defined;

* the normal and tangential components of the vector of the contact forces are corrected to satisfy the unilateral restrictions

$$ q_3^{k+1}(x_1, t) = P_\tau[q_3^k(x_1, t) \quad \rho_\tau \partial_t \Delta u_3^{k+1}(x_1, t)], \quad (3.30) $$

where

$$ P_\tau[q_3] = \begin{cases} q_3(x_3, t) & \text{if } |q_3| \le k_\tau q_n(x_3, t) \\ k_\tau q_n(x_3, t)q_3/|q_3| & \text{if } |q_3| > k_\tau q_n(x_3, t) \end{cases} \quad (3.31) $$
is the operator of the orthogonal projection onto set $|q_3| \le k_\tau q_n(x_3, t)$, coefficient $\rho_\tau$ 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 elasto dynamic contact problems have been presented in [57]. Particularly we find that coefficient $\rho_\tau$ has a significant influence on algorithm convergence. In the case that $\rho_\tau$ 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 $\omega$ propagating in the plane $R^2$ 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 $\mathbf{x}_\alpha$, multiplied by factor $e^{-i\omega 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 $\omega$

$$ 
\begin{aligned}
u_3(\mathbf{x}_\alpha, t) &= \operatorname{Re} \left\{ \sum_{-\infty}^{\infty} u_3^k(\mathbf{x}_\alpha) e^{i\omega_k t} \right\}, \\
\sigma_{3\beta}(\mathbf{x}_\alpha, t) &= \operatorname{Re} \left\{ \sum_{-\infty}^{\infty} \sigma_{3\beta}^k(\mathbf{x}_\alpha) e^{i\omega_k t} \right\}
\end{aligned}
\quad (4.1) $$

where

$$ 
\begin{aligned}
u_3^k(\mathbf{x}_\alpha) &= \frac{\omega_k}{2\pi} \int_0^T u_3(\mathbf{x}_\alpha, t) e^{i\omega_k t} dt, \\
\sigma_{3\beta}^k(\mathbf{x}_\alpha) &= \frac{\omega_k}{2\pi} \int_0^T \sigma_{3\beta}(\mathbf{x}_\alpha, t) e^{i\omega_k t} dt
\end{aligned}
\quad (4.2) $$

In the same way, traction $p_3(x_1, t)$ on the crack edges and their opening $\Delta u_3(x_1, t)$ may be expanded into Fourier series

$$ 
\begin{aligned}
p_3(x_1, t) &= \operatorname{Re} \left\{ \sum_{-\infty}^{\infty} p_3^k(x_1) e^{i\omega_k t} \right\}, \\
\Delta u_3(x_1, t) &= \operatorname{Re} \left\{ \sum_{-\infty}^{\infty} \Delta u_3^k(x_1) e^{i\omega_k t} \right\}
\end{aligned}
\quad (4.3) $$

where

$$ 
\begin{aligned}
p_3^k(x_1) &= \frac{\omega_k}{2\pi} \int_0^T p_3(x_1, t) e^{i\omega_k t} dt, \\
\Delta u_3^k(x_1) &= \frac{\omega_k}{2\pi} \int_0^T \Delta u_3(x_1, t) e^{i\omega_k t} dt
\end{aligned}
\quad (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

from of

$$ 
\begin{aligned}
& \partial_\rho \partial_\beta u_3^k + \frac{\omega_k^2}{c_2^2} u_3^k + \frac{\rho_k}{\mu} b_3^k = 0 \quad \text{for } k=0, \pm 1, \pm 2, \dots, \infty \quad \forall \mathbf{x}_\alpha \in \Omega
\end{aligned}
\quad (4.5) $$

In [20,21,23, 24] it was shown that the Fourier series expansions of the displacement discontinuity $\Delta u_3^k(\mathbf{x}_\alpha)$ and traction $p_3^k(\mathbf{x}_\alpha)$ are related by the boundary integral equations (BIE) in the form of

$$ 
\begin{aligned}
& p_3^k(\mathbf{x}_\alpha) = F.P. \int_{\Omega} F_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) \Delta u_3^k(\mathbf{y}_\alpha) d\Omega \\
& \text{for } k=0, \pm 1, \pm 2, \dots, \pm \infty \quad \forall \mathbf{x}_\alpha \in \Omega
\end{aligned}
\quad (4.6) $$

The kernels $F_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k)$ may be obtained from the fundamental solutions $U_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_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(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) = \frac{i}{4\mu} H_0^{(1)}(l_2^k), \quad l_2^k = r\omega_k/c_2 \quad (4.7) $$

where $H_0^1(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 $\mathbf{x}_\alpha$ and $\mathbf{y}_\alpha$ respectively to the fundamental solutions $U_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k)$ we obtain

$$ F_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) = \mu^2 \partial_n \partial_n U_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) \quad (4.8) $$

where $\partial_n = n_\alpha \partial_\alpha$ is the normal derivative.

The normal derivatives of the fundamental solutions $U_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k)$ are calculated by the equations

$$ 
\begin{aligned}
\partial_n U_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) &= \frac{n_\alpha x_\alpha dU_3(r, \omega_k)}{r dr}, \\
\partial_n \partial_n U_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha) &= \frac{n_\alpha n_\beta}{r} \left( \delta_{\alpha\beta} \frac{x_\alpha x_\beta}{r^2} \right) \frac{dU_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha)}{dr} \\
&\quad + \frac{n_\alpha x_\alpha n_\beta x_\beta d^2 U_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha)}{r dr^2}
\end{aligned}
\quad (4.9) $$

Taking into account that in the case of the plane finite crack shown in Fig. 1 vectors of outward normal are $n_1 = 0$, $n_2 = 1$ and $x_2 = 0$, the normal derivative for the fundamental solution $U_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k)$ have the form of

$$ \partial_n \partial_n U_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) = \frac{dU_3(r, \omega_k)}{r dr} \quad (4.10) $$

Using the formula for the derivative of the Bessel function $H_0^1(z)$ [1]

$$ \frac{d}{dr} H_0^{(1)}\left(\frac{\omega r}{c_2}\right) = \frac{\omega}{c_2} H_1^{(1)}\left(\frac{\omega r}{c_2}\right) \quad (4.11) $$

we obtain the kernel $F_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k)$. For our specific case it has the form of

$$ F_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) = \mu \frac{i\omega_k H_1^{(1)}(l_1^k)}{r c_2} \quad (4.12) $$

The kernel is a complex value function. It may be represented in the form of

$$ F_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) = F_3^{\text{Re}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) + i F_3^{\text{Im}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) \quad (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

$$ H_\nu^{(1)} = J_\nu(z) + iY_\nu(z) \quad (4.14) $$

Now the real and imaginary parts of the fundamental solution have the form of

$$ F_3^{\text{Re}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) = \mu \frac{\omega_k}{r c_2} J_1(l_2), \quad F_3^{\text{Im}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) = \mu \frac{\omega_k}{r c_2} J_1(l_2) \quad (4.15) $$
Let us consider in detail the structure of the kernel $F_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k)$ and analyze its behavior when $\mathbf{x}_\alpha \to \mathbf{y}_\alpha$. The Bessel functions of the first and second kind and the first order may be represented by the series expansion [1]

$$ 
\begin{aligned}
J_1(z) &= \frac{z}{2} \left[ 1 + \frac{1}{1!2!} \left(\frac{z}{2}\right)^2 + \frac{1}{2!3!} \left(\frac{z}{2}\right)^4 \dots \right] \\
Y_1(z) &= \frac{2}{\pi} \left( \ln(z/2) + \gamma J_1(z) + \frac{z}{2\pi} \left[ 1 + \left(1+1/2\right) \frac{1}{1!2!} \left(\frac{z}{2}\right)^2 + \dots \right] \right) \frac{2}{\pi z}
\end{aligned} 
\quad (4.16) 
$$

From these representations it follows that $J_1(z) \to z$ and $Y_1(z) \to 1/z$ for $z \to 0$ and therefore for

$$ F_3^{\text{Re}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) \to \frac{1}{r^2}, \quad F_3^{\text{Im}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) \to 0 \quad \text{for } \mathbf{x}_\alpha \to \mathbf{y}_\alpha \quad (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 $p_3^k(\mathbf{x}_\alpha)$ and $\Delta u_3^k(\mathbf{x}_\alpha)$. Separating real and imaginary parts

$$ 
\begin{aligned}
p_3^k(\mathbf{x}_\alpha) &= p_3^{\text{Re}}(\mathbf{x}_\alpha, k) + ip_3^{\text{Im}}(\mathbf{x}_\alpha, k), \\
\Delta u_3^k(\mathbf{x}_\alpha) &= \Delta u_3^{\text{Re}}(\mathbf{x}_\alpha, k) + i\Delta u_3^{\text{Im}}(\mathbf{x}_\alpha, k)
\end{aligned} 
\quad (4.18) 
$$

and using representation (4.13) for the kernels $F_3(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k)$ we obtain the system of integral equations for each $k$ in the form of

$$ 
\begin{aligned}
p_3^{\text{Re}}(\mathbf{y}_\alpha, k) &= F.P. \int_{\Omega} F_3^{\text{Re}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) \Delta u_3^{\text{Re}}(\mathbf{x}_\alpha, k) d\Omega \\
& \quad + \int_{\Omega} F_3^{\text{Im}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) \Delta u_3^{\text{Im}}(\mathbf{x}_\alpha, k) d\Omega, \\
p_3^{\text{Im}}(\mathbf{y}_\alpha, k) &= \int_{\Omega} F_3^{\text{Im}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) \Delta u_3^{\text{Re}}(\mathbf{x}_\alpha, k) d\Omega \\
& \quad + F.P. \int_{\Omega} F_3^{\text{Re}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k) \Delta u_3^{\text{Im}}(\mathbf{x}_\alpha, k) d\Omega
\end{aligned} 
\quad (4.19) 
$$

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

$$ 
\begin{aligned}
p_3^{\text{Re}}(\mathbf{y}_r, \omega_k) &= \sum_{n=1}^{N} [F_3^{\text{Re}}(\mathbf{x}_n, \mathbf{y}_r, \omega_k) \Delta u_3^{\text{Re}}(\mathbf{x}_n, k) \\
& \quad F_3^{\text{Im}}(\mathbf{x}_n, \mathbf{y}_r, \omega_k) \Delta u_3^{\text{Im}}(\mathbf{x}_n, k)], \\
p_3^{\text{Im}}(\mathbf{y}_r, k) &= \sum_{n=1}^{N} [F_3^{\text{Im}}(\mathbf{x}_n, \mathbf{y}_r, \omega_k) \Delta u_3^{\text{Re}}(\mathbf{x}_n, k) \\
& \quad F_3^{\text{Re}}(\mathbf{x}_n, \mathbf{y}_r, \omega_k) \Delta u_3^{\text{Im}}(\mathbf{x}_n, k)]
\end{aligned} 
\quad (4.20) 
$$

Their coefficients are defined by the equations

$$ 
\begin{aligned}
F_k^{\text{Re}}(\mathbf{x}_n, \mathbf{y}_r, \Delta_n) &= \int_{2\Delta_n} F_3^{\text{Re}}(\mathbf{x}, \mathbf{y}_r, \omega_k) d\Omega, \\
F_k^{\text{Im}}(\mathbf{x}_n, \mathbf{y}_r, \Delta_n) &= \int_{2\Delta_n} F_3^{\text{Im}}(\mathbf{x}, \mathbf{y}_r, \omega_k) d\Omega
\end{aligned} 
\quad (4.21) 
$$

Let us consider in detail how to calculate these coefficients. If points $\mathbf{x}_\alpha$ and $\mathbf{y}_\alpha$ 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 $\mathbf{x}_\alpha$ and $\mathbf{y}_\alpha$, belong to the same boundary elements some integrals in (4.21) are hypersingular. For calculation of these hypersingular integrals above mentioned regularization method

will be used. For that purpose we present the kernel $F_3^{\text{Re}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, \omega_k)$ in the form of

$$ F_k^{\text{Re}}(\mathbf{x}^n, \mathbf{y}^r, \Delta_n) = \int_{2\Delta_n} [F_3^{\text{Re}}(\mathbf{x}, \mathbf{y}^r, k) F_3^{\text{St}}(\mathbf{x}, \mathbf{y}^r)] dx + F.P. \int_{2\Delta_n} F_3^{\text{St}}(\mathbf{x}, \mathbf{y}^r) dx \quad (4.22) $$

Here $F_3^{\text{St}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha) = 2\mu/\pi r^2$ is the fundamental solution for the elastostatic problem, extracted from $F_3^{\text{Re}}(\mathbf{x}_\alpha, \mathbf{y}_\alpha, k)$ by considering the series expansion (4.16) and the asymptotic analysis for the case of $\omega_k \to 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

$$ F_3^{\text{St}}(\mathbf{x}^n, \mathbf{y}^r, \Delta_n) = \frac{4\mu}{\pi} \frac{\Delta_n}{(\mathbf{x}^n \cdot \mathbf{y}^r)^2} \Delta_n^2 \quad (4.23) $$

The residual 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])

$$ 
\begin{aligned}
J_1^*(m, n, c_\alpha) &= \int_{-\Delta_n}^{\Delta_n} (\omega/c_\alpha r) J_1(l_\alpha^*) d\beta \\
&= (\omega/2c_\alpha)^2 \left[ 2\Delta_n + \sum_{k=0}^{\infty} \left( 1 \right)^{k \gamma_k^{mn} (\Delta_n, c_\alpha)} \frac{(k!) (k+1)!}{(k+1)!} \right] \\
Y_1^*(m, n, c_\alpha) &= \int_{-\Delta_n}^{\Delta_n} (\omega/c_\alpha r) Y_1(l_\alpha^*) d\beta = \gamma J_1^*(m, n, c_\alpha) \quad (\omega/2c_\alpha)^2 \\
& \quad \times \left( 2\Delta_n + \sum_{k=0}^{\infty} \left( 1 \right)^{k \gamma_k^{mn} (\Delta_n, c_\alpha)} \frac{(k!) (\psi(k+1)+\psi(k+2))}{(k+1)!} \right) \\
& \quad + 2(\omega/2c_\alpha)^2 \sum_{k=0}^{\infty} \frac{(1)^k}{k! (k+1)!} \\
& \quad \times (\eta_k^{mn} (\Delta_n, c_\alpha) \quad \gamma_k^{mn} (\Delta_n, c_\alpha) / (2k+1))
\end{aligned} 
\quad (4.24) 
$$

where $r^* = |y^m x^n + \beta|$, $r_{+\Delta} = |y^m x^n + \Delta_n|$, $r_{-\Delta} = |y^m x^n - \Delta_n|$ and

$$ 
\begin{aligned}
\gamma_k^{mn}(\Delta_n, c_\alpha) &= \int_{-\Delta_n}^{\Delta_n} (r^*\omega/2c_\alpha)^{2k} d\beta = \frac{(\omega/2c_\alpha)^{2k}}{2k+1} (r_{+\Delta}^{2k+1} + r_{-\Delta}^{2k+1}), \\
\eta_k^{mn}(\Delta_n, c_\alpha) &= \frac{(\omega/2c_\alpha)^{2k}}{2k+1} (r_{+\Delta}^{2k+1} \ln(r_{+\Delta}\omega/2c_\alpha) + r_{-\Delta}^{2k+1} \ln(r_{-\Delta}\omega/2c_\alpha))
\end{aligned} 
\quad (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 $\omega_k$ being low or high. The system BEM equation (4.20) does not lose stability with increasing $\omega_k$.

For convenience and compactness the BEM equation (4.20) can be presented in the form of

$$ \mathbf{p}^k = \mathbf{A}^k \cdot \Delta \mathbf{u}^k \quad \text{for } k=0, \pm 1, \pm 2, \dots, \infty \quad (4.26) $$

where

$$ 
\begin{bmatrix} p_3^{\text{Re}}(\mathbf{y}_1, \omega_k) \\ \vdots \\ p_3^{\text{Re}}(\mathbf{y}_N, \omega_k) \\ p_3^{\text{Im}}(\mathbf{y}_1, \omega_k) \\ \vdots \\ p_3^{\text{Im}}(\mathbf{y}_N, \omega_k) \end{bmatrix}, \quad \Delta \mathbf{u}^k = \begin{bmatrix} \Delta u_3^{\text{Re}}(\mathbf{y}_1, \omega_k) \\ \vdots \\ \Delta u_3^{\text{Re}}(\mathbf{y}_N, \omega_k) \\ \Delta u_3^{\text{Im}}(\mathbf{y}_1, \omega_k) \\ \vdots \\ \Delta u_3^{\text{Im}}(\mathbf{y}_N, \omega_k) \end{bmatrix} \quad (4.27) 
$$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Crack opening in TD for impulse load.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Crack opening in the TD for different number of BE.</figcaption></figure>

$$A^k \begin{bmatrix} F_3^{\text{Re}}(\mathbf{x}_1 - \mathbf{y}_1, \omega_k) & \dots & F_3^{\text{Re}}(\mathbf{x}_N - \mathbf{y}_1, \omega_k) & F_3^{\text{Im}}(\mathbf{x}_1 - \mathbf{y}_1, \omega_k) & \dots & F_3^{\text{Im}}(\mathbf{x}_N - \mathbf{y}_1, \omega_k) \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ F_3^{\text{Re}}(\mathbf{x}_1 - \mathbf{y}_N, \omega_k) & \dots & F_3^{\text{Re}}(\mathbf{x}_N - \mathbf{y}_N, \omega_k) & F_3^{\text{Im}}(\mathbf{x}_1 - \mathbf{y}_N, \omega_k) & \dots & F_3^{\text{Im}}(\mathbf{x}_N - \mathbf{y}_N, \omega_k) \\ F_3^{\text{Im}}(\mathbf{x}_1 - \mathbf{y}_1, \omega_k) & \dots & F_3^{\text{Im}}(\mathbf{x}_N - \mathbf{y}_1, \omega_k) & F_3^{\text{Re}}(\mathbf{x}_1 - \mathbf{y}_1, \omega_k) & \dots & F_3^{\text{Re}}(\mathbf{x}_N - \mathbf{y}_1, \omega_k) \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ F_3^{\text{Im}}(\mathbf{x}_1 - \mathbf{y}_N, \omega_k) & \dots & F_3^{\text{Im}}(\mathbf{x}_N - \mathbf{y}_N, \omega_k) & F_3^{\text{Re}}(\mathbf{x}_1 - \mathbf{y}_N, \omega_k) & \dots & F_3^{\text{Re}}(\mathbf{x}_N - \mathbf{y}_N, \omega_k) \end{bmatrix} \quad (4.28)$$

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 algorithm for the solution of the problem in the FD includes the following steps:

* the initial distribution of the contact forces $q_3^0(x_1, t)$ in the TD is assigned;

* the Fourier coefficients $q_3(x_1, \omega_k)$ in the FD are calculated using formula (4.4);

* the FD BEM equations (4.20) and (4.26) are solved and $\Delta u_3(x_1, \omega_k)$ is calculated in the FD;

* the TD value $\Delta u_3^n(\mathbf{x}, t)$ is calculated using formula (4.3);

* the tangential components of the vector of contact forces are corrected to satisfy the unilateral restrictions

$$q_3^{k+1}(x_1, t) = P_{\tau}[q_3^k(x_1, t) \quad \rho_{\tau} \partial_t \Delta u_3^{k+1}(x_1, t)] \quad (4.29)$$

where

$$P_{\tau}[q_3] = \begin{cases} q_3(x_3, t) & \text{if } q_3 \le k_{\tau}q_n(x_3, t) \\ k_{\tau}q_n(x_3, t)q_3/q_3 & \text{if } |q_3| > k_{\tau}q_n(x_3, t) \end{cases} \quad (4.30)$$
is the operator of the orthogonal projection onto set $|q_3| \le k_r q_n(x_3, t)$, coefficient $\rho_r$ 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 $\mu = 1$, $c_2 = 1$ and $\nu = 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

$$u_3(x_1) = \frac{4p_3l}{\mu} \sqrt{1} (x_1/l)^2 \quad (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

<figure><img src="image_5.png" /><figcaption>Fig. 4. Frictional contact force and crack opening for $k_2$ 0.25 and $k_r$ 0.2.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Frictional contact force and crack opening for $k_2$ 0.25 and $k_r$ 0.4.</figcaption></figure>
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 by 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.

$q_n = 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 $\nu = 0.25$, specific density $\rho = 7800$ kg/m$^3$.

The crack opening and frictional contact force in the middle point of the crack ($x_1 = 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

<figure><img src="image_7.png" /><figcaption>Fig. 6. Frictional contact force and crack opening for $k_2$ 0.75 and $k_3$ 0.2.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Frictional contact force and crack opening for $k_2$ 0.75 and $k_3$ 0.4.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Frictional contact force for $k_2$ 0.25 and $k_r$ 0.2.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Crack opening for $k_2$ 0.25 and $k_r$ 0.2.</figcaption></figure>

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 friction 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

$$u_3 = \frac{K_{III}}{\mu} \sqrt{\frac{2\epsilon}{\pi}} \sin(\theta/2) \quad (5.2)$$

Using limit transition for $\theta = \pi$ we obtain

$$K_{III} = \lim_{\epsilon \to 0} \frac{\mu \sqrt{\pi}}{\sqrt{2\epsilon}} u_3(l, \epsilon, t) \quad (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_r = 0.2$ and $k_r = 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
<figure><img src="image_11.png" /><figcaption>Fig. 10. Frictional contact force for $k_2$ 0.25 and $k_r$ 0.4.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Crack opening for $k_2$ 0.25 and $k_r$ 0.4.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Number of iterations at each time step.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Stress intensity factor against wave number.</figcaption></figure>
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 formula tion are much better. For example convergence of the algorithm in the TD formulation is more sensitive to the mesh refinement and coefficient $ρ_s$ 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 2D, homogeneous, isotropic and linear elastic solids by con sidering 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 com puted 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 elasto dynamic 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 mixed-mode dynamic crack analysis. Computer and Structures 1996;59(6):1021-31.

[14] Freund LB. 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] Garcia-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 two-dimensional 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.
[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 Techiccheckie 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.

[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.