# A comparative study of three BEM for transient dynamic crack analysis of 2-D anisotropic solids

Felipe García-Sánchez, Chuanzeng Zhang

► To cite this version:

Felipe García-Sánchez, Chuanzeng Zhang. A comparative study of three BEM for transient dynamic crack analysis of 2-D anisotropic solids. Computational Mechanics, 2007, 40 (4), pp.753 - 769. 10.1007/s00466-006-0137-7. http-01382198

HAL Id: hal-01382198

https://hal.science/hal-01382198

Submitted on 19 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
# A comparative study of three BEM for transient dynamic crack analysis of 2-D anisotropic solids

Felipe García-Sánchez · Chuanzeng Zhang

**Abstract** Three different boundary element methods (BEM) for transient dynamic crack analysis in two-dimensional (2-D), homogeneous, anisotropic and linear elastic solids are presented. Hypersingular traction boundary integral equations (BIEs) in frequency-domain, Laplace-domain and time-domain with the corresponding elastodynamic fundamental solutions are applied for this purpose. In the frequency-domain and the Laplace-domain BEM, numerical solutions are first obtained in the transformed domain for discrete frequency or Laplace-transform parameters. Time-dependent results are subsequently obtained by means of the inverse Fourier-transform and the inverse Laplace-transform algorithm of Stehfest. In the time-domain BEM, the quadrature formula of Lubich is adopted to approximate the arising convolution integrals in the time-domain BIEs. Hypersingular integrals involved in the traction BIEs are computed through a regularization process that converts the hypersingular integrals to regular integrals, which can be computed numerically, and singular integrals which can be integrated analytically. Numerical results for the dynamic stress intensity factors are presented and discussed for a finite crack in an infinite domain subjected to an impact crack-face loading.

**Keywords** Hypersingular traction BEM · Frequency-domain BEM · Laplace-domain BEM · Time-domain BEM · 2-D anisotropic solids

## 1 Introduction

Dynamic crack analysis has many important applications in engineering sciences such as in fracture and damage mechanics, quantitative non-destructive material testing, geophysics and geomechanics. Since analytical solutions for dynamic crack problems can be obtained only for very simple cases, numerical methods have to be applied in general to solve the arising initial boundary-value problems. Among many available numerical methods, the boundary element method (BEM) is an accurate and efficient numerical tool for dynamic crack analysis, at least for homogeneous and linear elastic solids.

From the mathematical points of view, a crack in two-dimensional (2-D) elastic solids is a line with two coincident faces, which leads to a degeneration of the classical displacement BEM formulation over both crack-faces. This problem can be avoided by using the dual BEM, where the displacement boundary integral equations (DBIEs) are used over one of the crack-faces while the traction boundary integral equations (TBIEs) are applied to other crack-face. TBIEs can be obtained by the partial differentiation of the DBIEs and the subsequent application of the Hooke's law. Another remedy to overcome the degeneration of the DBIEs for crack analysis is the use of the hypersingular TBIEs only on one of the crack-faces, where the crack-opening-displacements (CODs) are the fundamental unknowns.

<footnote>F. García-Sánchez
Departamento de Ingeniería Civil, de Materiales y Fabricación, Universidad de Málaga,
Plaza El Ejido S/N 29013 Málaga, Spain
e-mail: fgsanchez@uma.es</footnote>

<footnote>Ch. Zhang(✉)
Department of Civil Engineering, University of Siegen,
Paul-Bonatz-Str. 9-11, 57076 Siegen, Germany
e-mail: c.zhang@uni-siegen.de</footnote>
quantities. For cracked solids with finite boundaries, a combined use of both the DBIEs and the TBIEs is promising. In this case, the TBIEs are applied to one of the crack-faces while the DBIEs are used on the external boundaries of the cracked solids. An overview of different possibilities for crack analysis in cracked elastic solids with bounded domains by BEM can be found in [31]. The numerically computed CODs can be used in a data post-processing to obtain the stress intensity factors (SIFs), which are the most important crack-tip characterizing parameters in linear elastic fracture mechanics.

Over the past years, three different BEM formulations, namely, the frequency-domain [5,9], the Laplace-domain [3] and the time-domain [2,24,25,29,33] BEM, are often applied to transient elastodynamic crack analysis. To analyze their accuracy and efficiency, a comparative study is performed in this paper. Hypersingular TBIEs are applied for this purpose. A collocation method for the spatial discretization of the hypersingular BIEs is adopted. Hypersingular integrals are dealt with by a regularization technique based on a suitable variable change [8,9]. In the frequency-domain and the Laplace-domain BEM, hypersingular BIEs in the transformed domain are first solved numerically for discrete frequency and Laplace-transform parameters. To obtain the time-dependent solutions, fast Fourier inverse transform and Stehfest's Laplace-inversion algorithm [23] are applied in the frequency-domain and the Laplace-domain BEM, respectively. The time-domain BEM uses a convolution quadrature formula of Lubich [11,12] for approximating the arising convolution integrals and it leads directly to time-dependent solutions. Numerical examples for computing transient elastodynamic SIFs in homogeneous and linear elastic solids of general anisotropy are presented to compare the accuracy and the efficiency of the three different BEM formulations.

## 2 Problem statement and elastodynamic BIEs

We consider a finite crack in a 2-D, homogeneous, linear elastic and anisotropic solid. The deformation of the cracked solid is assumed to be in a state of either generalized plane strain or generalized plane stress. In the absence of body forces, the cracked anisotropic solid satisfies the equations of motion [1]

$$ \sigma_{ij,j} = \rho \ddot{u}_i, \quad (1) $$

Hooke's law

$$ \sigma_{ij} = C_{ijkl} u_{k,l}, \quad (2) $$

the initial conditions

$$ u_i(\mathbf{x}, t) = \dot{u}_i(\mathbf{x}, t) = 0, \quad \text{for } t=0, \quad (3) $$

and the traction boundary conditions on the crack-faces

$$ p_i(\mathbf{x}, t) = \sigma_{ij}(\mathbf{x}, t)n_j(\mathbf{x}) = p_i^*(\mathbf{x}, t), \quad \mathbf{x} \in \Gamma_c. \quad (4) $$

In Eqs. (1)–(4), $\sigma_{ij}$ and $u_i$ denote the stress and the displacement components, $C_{ijkl}$ is the fourth order elasticity tensor, $p_i(\mathbf{x}, t)$ is the traction vector, $n_j$ is the outward unit normal vector, $\Gamma_c = \Gamma_c^+ + \Gamma_c^-$ represents the upper and the lower crack-faces, $\rho$ is the mass density, and $p_i^*(\mathbf{x}, t)$ is the prescribed crack-face loading, respectively. Also, a comma after a quantity represents partial derivatives with respect to spatial variables, and superscript dots stand for the time differentiations of the quantity. Unless otherwise stated, the conventional summation rule over double indices is implied, and the indices $i$ and $j$ take the values 1 and 2.

The displacements can be represented by the following boundary integral

$$ u_i(\xi, t) = - \int_{\Gamma_c^+} p_{ij}^*(\xi, \mathbf{x}, t) * \Delta u_j(\mathbf{x}, t) d\Gamma, \quad (5) $$

where $\mathbf{x} = (x_1, x_2)$ and $\xi = (\xi_1, \xi_2)$ represent the field and the source points, $p_{ij}^*$ are the traction fundamental solutions, $\Delta u_j(\mathbf{x}, t)$ are the CODs defined by

$$ \Delta u_i(\mathbf{x}, t) = u_i(\mathbf{x} \in \Gamma_c^+, t) - u_i(\mathbf{x} \in \Gamma_c^-, t), \quad (6) $$

and an asterisk * denotes Riemann convolution which is defined by

$$ g(\mathbf{x}, t) * h(\mathbf{x}, t) = \int_{0}^{t} g(\mathbf{x}, t - \tau)h(\mathbf{x}, \tau)d\tau. \quad (7) $$

The traction fundamental solution $p_{ij}^*$ is related to the displacement fundamental solution $u_{ij}^*$ by

$$ p_{ij}^* = C_{jkln} \frac{\partial u_{il}^*}{\partial x_n} n_k(\mathbf{x}). \quad (8) $$

Substituting Eq. (5) into Hooke's law (2) we obtain an integral representation formula for the traction components as

$$ p_j(\xi) = - \int_{\Gamma_c^+} s_{ij}^*(\xi, \mathbf{x}, t) * \Delta u_j(\mathbf{x}, t) d\Gamma, \quad (9) $$

where

$$ s_{ij}^* = C_{ikln} \frac{\partial p_{lj}^*}{\partial \xi_n} n_k(\xi), \quad (10) $$

and $n_k(\xi)$ represents the outward unit normal vector to the boundary at the collocation point. By taking the
limiting process $\xi \to \Gamma_c^+$, the following time-domain
traction BIEs are obtained

$$
p_j(\xi) = - \int_{\Gamma_c^+} s_{ij}^*(\xi, \mathbf{x}, t) * \Delta u_j(\mathbf{x}, t) d\Gamma, \quad (11)
$$

where $\int \cdot \int$ stands for the finite-part integral of Hadamard.

Frequency-domain and Laplace-domain TBIEs can
be obtained directly by applying the Fourier-transform
and the Laplace-transform to Eq. (11).

**3 Elastodynamic fundamental solutions**

Unlike in the isotropic case, elastodynamic fundamen-
tal solutions for homogeneous, anisotropic and linear
elastic solids cannot be given in closed-forms and they
are much more complicated. For 2-D case, fundamen-
tal solutions can be represented by a line-integral over
a unit circle, while for 3-D case they can be given as
a surface-integral over a unit sphere. In the present
analysis, two different elastodynamic fundamental solu-
tions are adopted, namely, the frequency-domain and
the Laplace-domain dynamic fundamental solutions.

Two-dimensional frequency-domain displacement
fundamental solution has the following expression [27]

$$
u_{ij}^*(\xi, \mathbf{x}, \omega) = \frac{1}{8\pi^2} \int_{|\eta|=1} \sum_{m=1}^{2} \frac{\gamma_{ij}^m}{\rho c_m^2} \Phi(\omega|\eta \cdot (\mathbf{x} - \xi)|/c_m) dS_{\eta}, \quad (12)
$$

where $\eta = (\eta_1, \eta_2)$ is the wave propagation vector, $\omega$ is
the circular frequency, $\gamma_{ij}^m = V_{im}V_{jm}$ (no sum over $m$)
is the projection operator, with $V_{im}$ and $\rho c_m^2$ being the
eigenvectors and the eigenvalues of the Christoffel ten-
sor $\Gamma_{ij}$

$$
[\Gamma_{ij}(\eta_1, \eta_2) - \rho c_m^2 \delta_{ij}] V_{jm} = 0, \quad (\text{no sum over } m), \quad (13)
$$

in which $\delta_{ij}$ denotes the Kronecker delta, $c_m$ represents
the phase velocities of elastic waves, and the Christoffel
tensor is defined by

$$
\Gamma_{ij}(\eta_1, \eta_2) = C_{ikjl} \eta_k \eta_l. \tag{14}
$$

The function $\Phi(\zeta)$ is given by

$$
\Phi(\zeta) = i\pi e^{i\zeta} - 2[\cos(\zeta) \operatorname{ci}(\zeta) + \sin(\zeta) \operatorname{si}(\zeta)], \quad (15)
$$

where ci and si denote the cosine and the sine integrals
which are defined by

$$
\mathrm{ci}(\zeta) = - \int_{\zeta}^{\infty} \frac{\cos t}{t} \mathrm{d}t; \quad \mathrm{si}(\zeta) = - \int_{\zeta}^{\infty} \frac{\sin t}{t} \mathrm{d}t, \qquad (16)
$$

with $\zeta$ being real and $f$ indicating the Cauchy principal
value integral.

Two-dimensional Laplace-domain displacement fun-
damental solution has been recently published in [28] for
piezoelectric materials. From [28], 2-D Laplace-domain
displacement fundamental solution for homogeneous,
anisotropic and linear elastic solids can be obtained as

$$
u_{ij}^*(\xi, \mathbf{x}, s) = \frac{1}{8\pi^2} \int_{|\eta|=1} \sum_{m=1}^{2} \frac{\gamma_{ij}^m}{\rho c_m^2} \Psi(s|\eta \cdot (\mathbf{x} - \xi)|/c_m) dS_{\eta}, \quad (17)
$$

where *s* is the Laplace-transform parameter,

$$
\Psi(z) = -[e^{-z} \mathrm{Ei}(z) + e^{z} \mathrm{Ei}(-z)], \quad (18)
$$

with z being complex and Ei being the exponential integral defined by

$$
\mathrm{Ei}(z) = - \int_{-1}^{\infty} \frac{e^{-zt}}{t} \mathrm{d}t; \quad \mathrm{Ei}(-z) = - \int_{1}^{\infty} \frac{e^{-zt}}{t} \mathrm{d}t. \quad (19)
$$

It can be easily shown that both functions $\Phi$ and $\Psi$
have a logarithmic singularity. To deal with this singular-
ity it is advantageous to split the fundamental solutions
into a singular static part and a regular dynamic part as
[27,28]

$$
u_{ij}^*(\xi, \mathbf{x}, \omega) = u_{ij}^{R*}(\xi, \mathbf{x}, \omega) + u_{ij}^{S*}(\xi, \mathbf{x}), \quad (20)
$$

$$
u_{ij}^*(\xi, \mathbf{x}, s) = u_{ij}^{R*}(\xi, \mathbf{x}, s) + u_{ij}^{S*}(\xi, \mathbf{x}), \quad (21)
$$

where the superscripts *R* and *S* denote the regular
dynamic part and the singular static part, respectively.
The singular static part is independent of *ω* and *s* and
has the following expression

$$
u_{ij}^{S*}(\xi, \mathbf{x}) = -\frac{1}{4\pi^2} \int_{|\eta|=1} \sum_{m=1}^{2} \frac{\gamma_{ij}^m}{\rho c_m^2} \log |\eta \cdot (\mathbf{x} - \xi)| dS_{\eta}, \quad (22)
$$

while the regular dynamic part can be written as

$$
u_{ij}^{R*}(\xi, \mathbf{x}, \omega) = \frac{1}{8\pi^2} \int_{|\eta|=1} \sum_{m=1}^{2} \frac{\gamma_{ij}^m}{\rho c_m^2} \Phi^R(\omega/c_m, |\eta \cdot (\mathbf{x}-\xi)|) dS_{\eta}, \quad (23)
$$

$$
u_{ij}^{R*}(\xi, \mathbf{x}, s) = \frac{1}{8\pi^2} \int_{|\eta|=1} \sum_{m=1}^{2} \frac{\gamma_{ij}^m}{\rho c_m^2} \Psi^R(s/c_m, |\eta \cdot (\mathbf{x}-\xi)|) dS_{\eta}, \quad (24)
$$

with the regular continuous functions

$$
\Phi^R(x, y) = \Phi(xy) + 2 \log y, \quad (25)
$$

$$
\Psi^R(x, y) = \Psi(xy) + 2 \log y. \quad (26)
$$
The static displacement fundamental solution (22) can be reduced to the following explicit expression [7–9]

$$u_{ij}^{S*}(\xi, \mathbf{x}) = -\frac{1}{\pi} \mathrm{Re} \left[ \sum_{m=1}^{2} P_{jm} Q_{mi} \ln(z_m - z_m^0) - \ln(i + \mu_m) \right], \quad (27)$$

where

$$z_m = x_1 + \mu_m x_2, \quad z_m^0 = \xi_1 + \mu_m \xi_2 \quad (28)$$

are the complex counterpart of the integration and the collocation points, the matrices $P_{jm}$, $Q_{mi}$ and the complex constants $\mu_m$ are given in the Appendix.

By substituting Eqs. (20)–(24) into Eqs. (8) and (10), the corresponding traction and higher-order traction fundamental solutions $p_{ij}^*$ and $s_{ij}^*$ can be obtained as

$$p_{ij}^{*}(\xi, \mathbf{x}, \omega) = p_{ij}^{R*}(\xi, \mathbf{x}, \omega) + p_{ij}^{S*}(\xi, \mathbf{x}), \quad (29)$$

$$p_{ij}^{*}(\xi, \mathbf{x}, s) = p_{ij}^{R*}(\xi, \mathbf{x}, s) + p_{ij}^{S*}(\xi, \mathbf{x}), \quad (30)$$

$$s_{ij}^{*}(\xi, \mathbf{x}, \omega) = s_{ij}^{R*}(\xi, \mathbf{x}, \omega) + s_{ij}^{S*}(\xi, \mathbf{x}), \quad (31)$$

$$s_{ij}^{*}(\xi, \mathbf{x}, s) = s_{ij}^{R*}(\xi, \mathbf{x}, s) + s_{ij}^{S*}(\xi, \mathbf{x}), \quad (32)$$

where

$$p_{ij}^{S*}(\xi, \mathbf{x}) = \frac{1}{\pi} \mathrm{Re} \left[ \sum_{m=1}^{2} L_{jm} Q_{mi} \frac{\mu_m n_1(\mathbf{x}) - n_2(\mathbf{x})}{z_m - z_m^0} \right], \quad (33)$$

$$p_{ij}^{R*}(\xi, \mathbf{x}, \omega) = \frac{1}{8\pi^2} \int_{|\eta|=1} \sum_{m=1}^{2} \frac{S_{ij}^m}{\rho c_m^2} \cdot \frac{\omega}{c_m} \cdot \Phi'(\omega|\eta \cdot (\mathbf{x}-\xi)|/c_m) \times \mathrm{sign}[\eta \cdot (\mathbf{x}-\xi)] \mathrm{d}S_{\eta}, \quad (34)$$

$$p_{ij}^{R*}(\xi, \mathbf{x}, s) = \frac{1}{8\pi^2} \int_{|\eta|=1} \sum_{m=1}^{2} \frac{S_{ij}^m}{\rho c_m^2} \cdot \frac{s}{c_m} \cdot \Psi'(s|\eta \cdot (\mathbf{x}-\xi)|/c_m) \times \mathrm{sign}[\eta \cdot (\mathbf{x}-\xi)] \mathrm{d}S_{\eta}, \quad (35)$$

$$s_{ij}^{S*}(\xi, \mathbf{x}) = \frac{1}{\pi} \mathrm{Re} \left[ \sum_{m=1}^{2} T_{ij}^{m} \frac{\mu_{m} n_{1}(\mathbf{x}) - n_{2}(\mathbf{x})}{(z_{m} - z_{m}^{0})^{2}} \right], \quad (36)$$

$$s_{ij}^{R*}(\xi, \mathbf{x}, \omega) = \frac{1}{8\pi^2} \int_{|\eta|=1} \sum_{m=1}^{2} \frac{U_{ij}^m}{\rho c_m^2} \left(\frac{\omega}{c_m}\right)^2 \times \Phi(\omega|\eta \cdot (\mathbf{x}-\xi)|/c_m) \mathrm{d}S_{\eta}, \quad (37)$$

$$s_{ij}^{R*}(\xi, \mathbf{x}, s) = \frac{1}{8\pi^2} \int_{|\eta|=1} \sum_{m=1}^{2} \frac{U_{ij}^m}{\rho c_m^2} \left(\frac{s}{c_m}\right)^2 \times \Psi(s|\eta \cdot (\mathbf{x}-\xi)|/c_m) \mathrm{d}S_{\eta}. \quad (38)$$

In Eqs. (33)–(38), $\Phi'(\xi) = d\Phi(\xi)/d\zeta$, $\Psi'(z) = d\Psi(z)/dz$, and the auxiliary functions $L_{jm}$, $S_{ij}^m$, $T_{ij}^m$ and $U_{ij}^m$ are given in the Appendix.

## 4 Numerical implementation of the BEM

To solve the hypersingular TBIEs (11) as well as their counterparts in the frequency-domain and the Laplace-transformed domain, a collocation method for the spatial discretization by using quadratic elements is developed. Discontinuous quadratic elements are adopted over the crack-face. It should be mentioned here that the use of discontinuous elements is necessary in order to fulfill the C¹-continuity requirement of the CODS in the hypersingular TBIEs.

### 4.1 Treatment of hypersingular integrals

As mentioned previously, the TBIEs (11) and their counterparts in the transformed domain involve hypersingular integrals of the type $1/(z_m - z_m^0)^2$ when the integration point coincides with the collocation point. After the discretization, the hypersingular integrals to be computed have the following expression

$$I_m = \int_{\Gamma_e} \frac{\mu_m n_1(\mathbf{x}) - n_2(\mathbf{x})}{(z_m - z_m^0)^2} \phi \mathrm{d}\Gamma, \quad (39)$$

where $\Gamma_e$ is the boundary element under consideration, $\phi$ denotes the quadratic shape-function and $\mathbf{n}(\mathbf{x})$ is the outward unit normal vector to the boundary. By means of a suitable variable change [8,9], the hypersingular integral (39) can be regularized as described in the following.

Introducing a complex distance between the collocation and the integration points as a new variable

$$\chi_m = z_m - z_m^0 = (x_1 - \xi_1) + \mu_m (x_2 - \xi_2), \quad (40)$$

the Jacobian of this transformation is then given by

$$\frac{d\chi_m}{d\Gamma} = \frac{d\chi_m}{dx_1} \frac{dx_1}{d\Gamma} + \frac{d\chi_m}{dx_2} \frac{dx_2}{d\Gamma} = \mu_m n_1(\mathbf{x}) - n_2(\mathbf{x}), \quad (41)$$

where the following relations

$$\frac{dx_1}{d\Gamma} = \cos(\theta) = -n_2(\mathbf{x}), \quad \frac{dx_2}{d\Gamma} = \sin(\theta) = n_1(\mathbf{x}) \quad (42)$$

have been used (see Fig. 1).

By using Eqs. (40) and (41), the hypersingular integral (39) can be rewritten as

$$I_m = \int_{\Gamma_e} \frac{1}{\chi_m^2} \phi \mathrm{d}\chi_m. \quad (43)$$
<figure><img src="image_2.png" /><figcaption>Fig. 1 Differential boundary element geometry around the integration point with the outward unit normal vector</figcaption></figure>

Considering the shape function $\phi$ as a function of the complex variable $\chi_m$, and using the first two terms of its Taylor-series expansion at $\chi_m = 0$, i.e.,

$$ \phi(\chi_m) = \phi(\chi_m = 0) + \left. \frac{d\phi}{d\chi_m} \right|_{\chi_m=0} \chi_m + O(\chi_m^2) \\ = \phi_0 + \phi'_0 \chi_m + O(\chi_m^2), \quad (44) $$

the integral $I_m$ can be recast into

$$ I_m = \int_{\Gamma_e} \frac{\phi - (\phi_0 + \phi'_0 \chi_m)}{\chi_m^2} d\chi_m + \phi_0 \int_{\Gamma_e} \frac{1}{\chi_m^2} d\chi_m \\ + \phi'_0 \int_{\Gamma_e} \frac{1}{\chi_m} d\chi_m. \quad (45) $$

The first integral in Eq. (45) is regular and can be computed numerically by using standard Gaussian quadrature formula. The second and the third integral are hypersingular and strongly singular, but they can be evaluated analytically as

$$ \int_{\Gamma_e} \frac{1}{\chi_m^2} d\chi_m = - \left. \frac{1}{\chi_m} \right|_{\Gamma_e}, \\ \int_{\Gamma_e} \frac{1}{\chi_m} d\chi_m = \log(\chi_m) \bigg|_{\Gamma_e}. \quad (46) $$

## 4.2 Frequency-domain and Laplace-domain BEM

In the frequency-domain and the Laplace-domain BEM, the boundary value problem is first solved numerically for discrete values of the frequency and the Laplace-transform parameters. Subsequently, the corresponding time-domain solutions are obtained by using the fast Fourier inverse transform in the frequency-domain BEM and the inverse Laplace-transform algorithm of Stehfest [23] in the Laplace-domain BEM.

According to the Stehfest's inversion algorithm [23] a time-dependent function $f(t)$ can be approximated by

$$ f(t) = \frac{\ln 2}{t} \sum_{n=1}^{N} v_n \hat{f}(n \ln 2/t), \quad (47) $$

where $\hat{f}(s)$ is the Laplace-transform of $f(t)$, and

$$ v_n = (-1)^{n+N/2} \sum_{k=(n+1)/2}^{\min(n,N/2)} \frac{(2k)! k^{N/2}}{(N/2 - k)! k!(k-1)!(n-k)!(2k-n)!} \quad (48) $$

Stehfest suggested to use a single precision arithmetic and $N = 10$ for the truncation limit in order to obtain accurate results. In this work, $N=12$ and $24$ with a double precision arithmetic are used.

In the course of this study, we have also tested the Durbin's method for the inverse Laplace-transform, which has been preferred and suggested in some previous investigations on dynamic problems [10,14]. Our own experiences have confirmed that the Durbin's method is generally more accurate but also more complicated and computationally much more expensive than the Stehfest's algorithm. The Durbin's method uses a complex Laplace-transform parameter, while the Stehfest's algorithm applies a real Laplace-transform parameter and is easier to implement. For the dynamic crack problems considered in the present study, no significant changes in the numerical results have been noted by using these two different inversion methods. For this reason, the Stehfest's inversion algorithm is applied in our numerical examples, which will be presented in Sect. 6.

## 4.3 Time-domain BEM

In the time-domain BEM, the Riemann convolution integral

$$ f(t) = g(t) * h(t) = \int_{0}^{t} g(t - \tau)h(\tau)d\tau \quad (49) $$

is approximated by the quadrature formula of Lubich [11,12], which is given by

$$ f(n \Delta t) = \sum_{j=0}^{n} \bar{\omega}_{n-j}(\Delta t)h(j \Delta t), \quad n = 0, 1, 2, \dots, N, \quad (50) $$

where the time-interval $t$ is divided into $N$ equal time-steps $\Delta t$, and the weights $\bar{\omega}_n(\Delta t)$ are defined by

$$ \bar{\omega}_n(\Delta t) = \frac{r^{-n}}{N} \sum_{m=0}^{N-1} \hat{g}\left(\frac{\delta(\zeta_m)}{\Delta t}\right) e^{-2\pi i m N/N}. \quad (51) $$
In Eq. (51),

- $\hat{g}(\cdot)$ is the Laplace-transform of $g(t)$,

- $\delta(\zeta_m) = \sum_{j=1}^{2} (1 - \zeta_m)^j / j$,

- $\zeta_m = re^{2\pi im/N}$, and

- $r = \epsilon^{1/2N}$, with $\epsilon$ being the numerical error in computing the Laplace-transform $\hat{g}(\cdot)$.

It should be remarked here that a backward differential formula (BDF) of the second order is used for $\delta(\zeta_m)$, and the error in computing $\bar{\omega}_n(\Delta t)$ is of the order $O(\sqrt{\epsilon})$ [11,12]. In this analysis, different $\epsilon$-values between $10^{-6}$ and $10^{-12}$ have been tested to verify the effects of $\epsilon$ on the numerical accuracy. Our numerical tests have shown that in this $\epsilon$-range there are no notable influences of $\epsilon$ on the numerical results.

In contrast to the conventional time-domain BEM [24,25,31], the present method uses the Laplace-domain instead of the time-domain elastodynamic fundamental solutions. This is advantageous for cases where time-domain dynamic fundamental solutions are not available but their Laplace-transforms can be obtained. As representative examples we just mention the viscoelastic and the dynamic poroelastic problems, which have been investigated in details by Schanz [20], Gaul and Schanz [10] and Schanz [21] using the quadrature formula of Lubich [11,12].

For the present crack problem in an infinite anisotropic domain, the Laplace-domain system matrix, denoted by $\hat{\mathbf{A}}(s_m)$, can be computed by using the following equation

$$ \hat{\mathbf{A}}(s_m) = \sum_{e=1}^{E} \int_{\Gamma_e} \hat{s}_{ij}^* \left( \xi, \mathbf{x}, \frac{\zeta_m}{\Delta t} \right) \phi(\mathbf{x}) \, d\Gamma \quad (52) $$

According to Eq. (51), the system matrix at the nth time-step $\mathbf{A}^n$ can be obtained as

$$ \mathbf{A}^n = \frac{r^{-n}}{N} \sum_{m=0}^{N-1} \hat{\mathbf{A}}(s_m) e^{-2\pi im/N}, \quad (53) $$

where $s_m = \delta(\zeta_m)/\Delta t$.

Finally, a system of linear algebraic equations for the discrete CODs can be obtained as

$$ \sum_{j=0}^{n} \mathbf{A}^{n-j} \Delta \mathbf{u}^{j} = \mathbf{p}^{j}. \quad (54) $$

By invoking the zero initial conditions (3), Eq. (54) leads to the following explicit time-stepping scheme [29, 30,33]

$$ \Delta \mathbf{u}^n = (\mathbf{A}^0)^{-1} \left( \mathbf{p}^n - \sum_{j=1}^{n-1} \mathbf{A}^{n-j} \Delta \mathbf{u}^j \right) \quad (55) $$

for computing the unknown CODs at the nth time-step. In Eq. (55), $(\mathbf{A}^0)^{-1}$ is the inverse of the system matrix $\mathbf{A}^0$ at the time-step $n=0$.

## 5 Computation of SIFs

In all three BEM as presented in previous sections, straight quarter-point elements are adopted at the crack-tips in order to capture the local $\sqrt{r}$-behavior of the CODs near the crack-tips. This allows us to compute the dynamic SIFs very efficiently and accurately.

In the vicinity of the crack-tip, the displacement field has the following asymptotic expressions [22]

$$ u_1 = \sqrt{\frac{2r}{\pi}} \left\{ K_I \mathrm{Re} \left[ \frac{1}{\mu_1 - \mu_2} \left( \mu_1 p_2 \sqrt{\cos \theta + \mu_2 \sin \theta} - \mu_2 p_1 \sqrt{\cos \theta + \mu_1 \sin \theta} \right) \right] + K_{II} \mathrm{Re} \left[ \frac{1}{\mu_1 - \mu_2} \left( p_2 \sqrt{\cos \theta + \mu_2 \sin \theta} - p_1 \sqrt{\cos \theta + \mu_1 \sin \theta} \right) \right] \right\}, \quad (56) $$

$$ u_2 = \sqrt{\frac{2r}{\pi}} \left\{ K_I \mathrm{Re} \left[ \frac{1}{\mu_1 - \mu_2} \left( \mu_1 q_2 \sqrt{\cos \theta + \mu_2 \sin \theta} - \mu_2 q_1 \sqrt{\cos \theta + \mu_1 \sin \theta} \right) \right] + K_{II} \mathrm{Re} \left[ \frac{1}{\mu_1 - \mu_2} \left( q_2 \sqrt{\cos \theta + \mu_2 \sin \theta} - q_1 \sqrt{\cos \theta + \mu_1 \sin \theta} \right) \right] \right\}, \quad (57) $$

where $r$ and $\theta$ are polar coordinates with the origin at the crack-tip, $K_I$ and $K_{II}$ are the mode-I and the mode-II SIFs, and

$$ \begin{aligned} p_m &= b_{11}\mu_m^2 + b_{12} - b_{16}\mu_m, \\ q_m &= (b_{12}\mu_m^2 + b_{22} - b_{26}\mu_m)/\mu_m. \end{aligned} \quad (58) $$

In Eq. (58), $b_{ij}$ ($i, j = 1, 2, 6$) are the materials compliance coefficients.

By substituting Eqs. (56) and (57) into Eq. (6), we obtain the following relationship between the dynamic SIFs and the CODs

$$ \left\{ \begin{matrix} K_I \\ K_{II} \end{matrix} \right\} = \sqrt{\frac{\pi}{8d}} \cdot \frac{1}{\Delta} \left[ \begin{matrix} H_{11} & H_{12} \\ H_{21} & H_{22} \end{matrix} \right] \left\{ \begin{matrix} \Delta u_1 \\ \Delta u_2 \end{matrix} \right\}, \quad (59) $$
where

$$H_{ij} = \begin{bmatrix} \operatorname{Im}\left(\frac{q_1-q_2}{\mu_1-\mu_2}\right) & \operatorname{Im}\left(\frac{p_2-p_1}{\mu_1-\mu_2}\right) \\ \operatorname{Im}\left(\frac{\mu_1 q_2 - \mu_2 q_1}{\mu_1 - \mu_2}\right) & \operatorname{Im}\left(\frac{\mu_2 p_1 - \mu_1 p_2}{\mu_1 - \mu_2}\right) \end{bmatrix}, \quad (60)$$

$$\Delta = H_{11}H_{22} - H_{12}H_{21}. \quad (61)$$

In Eq. (59), $\Delta u_1$ and $\Delta u_2$ are the CODs at the collocation point closest to the crack-tip either in the frequency-domain, or in the Laplace-transformed domain or in the time-domain depending on the applied BEM, and $d$ is the distance between this collocation point and the crack-tip.

Using Eq. (59), dynamic SIFs can be obtained directly once the CODs have been computed numerically by the BEM. In the frequency-domain and the Laplace-domain BEM, an inverse transform is subsequently applied to obtain time-dependent SIFs. In contrast, time-domain BEM yields directly time-dependent SIFs.

## 6 Numerical results

To compare the three different BEM, both straight and curved cracks are analyzed. The straight crack has a finite length $2a$, and is considered to be either in an isotropic or a fully anisotropic unbounded plane subjected to a tensile or a shear impact crack-face loading (see Fig. 2). As an example for curved cracks, we consider a circular arc-shaped crack with a radius $r$ and opening angle $2\alpha$ embedded in a fully anisotropic plane, which is subjected to a radial impact crack-face loading as shown in Fig. 3. In all cases, the accuracy and the efficiency of the three BEM are compared. Furthermore the stability of the time-domain BEM is studied by using different time-steps.

Numerical calculations have been carried out by using ten elements as depicted in Fig. 4. For the straight crack, the size-ratio between two consecutive elements is constant, which is equal to 2 between the central and the crack-tip element. For the curved crack, the straight crack-tip elements have a length of $r\alpha/30$, while the first three elements from the crack-center have an angle $0.50\alpha$, $0.30\alpha$ and $0.15\alpha$, respectively. The size of the fourth element from the crack-center is determined by the remaining size of the crack. For convenience of presentation and discussion, the following abbreviations are introduced:

- FD = frequency-domain BEM,

- LD = Laplace-domain BEM,

- TD = Time-domain BEM,

- FFT = Fast Fourier Transform,

- LP = Laplace-transform parameter.

<figure><img src="image_3.png" /><figcaption>Fig. 2 A straight crack in an unbounded domain subjected to a tensile impact loading (left) and a shear impact loading (right)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3 A circular arc-shaped crack in an unbounded domain subjected to a tensile impact loading in radial direction</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 BEM meshes used for straight and arc-shaped (opening angle = $2\alpha$) cracks</figcaption></figure>

### 6.1 Isotropic solids

In order to compare our numerical results with available analytical solutions we first consider a straight finite crack of length $2a$ in an unbounded, homogeneous, isotropic and linearly elastic solid. The crack is subjected to a tensile impact crack-face loading $\sigma_{22}(t) = \sigma_0 H(t)$
<figure><img src="image_6.png" /><figcaption>Fig. 5 Normalized dynamic SIFs for a crack in an isotropic solid subjected to a tensile and a shear impact crack-face loading</figcaption></figure>

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

or a shear impact crack-face loading $\sigma_{12}(t) = \tau_0 H(t)$, where $\sigma_0$ and $\tau_0$ are the loading amplitudes and $H(t)$ is the Heaviside step function. Poisson's ratio is taken as 1/4.

Numerical results are presented and compared with the analytical solutions [26] in Figs. 5 and 6. Here, the dynamic SIFs are normalized by their corresponding static values $K_I^{st} = \sigma_0 \sqrt{\pi a}$ and $K_{II}^{st} = \tau_0 \sqrt{\pi a}$, and the dimensionless time $tc_L/a$ is used for convenience, where $c_L$ is the velocity of the longitudinal wave. In the time-domain BEM, a time-step $c_T \Delta t = a/20$ is selected. As we can see in Figs. 5 and 6, the numerical results obtained by all three BEM show very good agreement with the analytical solutions, though some small discrepancies occur near the peak values of the normalized dynamic SIFs.

To study the stability of the time-domain BEM, additional calculations for three different time-steps have been performed. The corresponding numerical results are presented in Fig. 6, which imply that the present time-domain BEM is quite insensitive to the selected time-steps. This is an important advantage over the classical time-domain BEM [24,25,31], which suffers from the stability problem.

## 6.2 General anisotropic solids

### 6.2.1 Straight crack

Next, we consider a straight finite crack of length $2a$ in an unbounded, homogeneous and linear elastic solid with a general anisotropy. The following elastic constants have

<figure><img src="image_8.png" /><figcaption>Fig. 6 Influences of time-steps on the normalized dynamic SIF for a crack in an isotropic solid subjected to a tensile impact crack-face loading</figcaption></figure>

been chosen in the numerical calculations

$$ \begin{bmatrix} C_{11} & C_{12} & C_{16} \\ C_{22} & C_{26} \\ \text{sym.} & C_{66} \end{bmatrix} = \begin{bmatrix} 95.46 & 28.93 & 44.67 \\ 25.91 & 15.56 \\ \text{sym.} & 32.68 \end{bmatrix} \text{ GPa}, \qquad (62) $$

The selected material constants correspond to a Graphite-epoxy composite [33]. The same BEM mesh as in the isotropic case has been applied.
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Normalized mode-I and mode-II dynamic SIFs for a crack in an anisotropic solid subjected to a tensile impact loading</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8** Normalized mode-I and mode-II dynamic SIFs for a crack in an anisotropic solid subjected to a shear impact loading</figcaption></figure>

For a tensile impact crack-face loading, the normalized mode-I and mode-II dynamic SIFs versus the dimensionless time $tc_T/a$ are presented in Fig. 7, while the corresponding numerical results for a shear impact crack-face loading are given in Fig. 8. Here, the shear wave velocity is defined by $c_T = \sqrt{C_{66}/\rho}$, and a time-step $c_T \Delta t = a/20$ is used in the time-domain BEM. The results are compared with that obtained by Zhang [33], who used a Galerkin-method for the spatial discretization.

The good agreement between our numerical results and that of Zhang [33] as shown in Figs. 7 and 8

verifies the accuracy of the three implemented BEM. As in the isotropic case, our time-domain BEM is quite stable.

### 6.2.2 Circular arc-shaped crack

Now we consider a circular arc-shaped crack embedded in an unbounded anisotropic domain with the same material constants as used in the preceding case. The crack is defined by the radius $r$ of the circle and the
<figure><img src="image_11.png" /><figcaption>**Fig. 9** Normalized mode-I and mode-II dynamic SIFs for an arc-shaped crack with $\alpha = 15^\circ$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10** Normalized mode-I and mode-II dynamic SIFs for an arc-shaped crack with $\alpha = 30^\circ$</figcaption></figure>

opening angle $2\alpha$ as shown in Fig. 3. A tensile impact crack-face loading in the radial direction in the form of $\sigma_r = \sigma_0 H(t)$ is considered.

The numerical results are presented as normalized mode-I and mode-II dynamic SIFs versus dimensionless time, which is defined as in the preceding case for a straight crack. The dynamic SIFs are normalized by $\sigma_0\sqrt{\pi a}$, where $\sigma_0$ is the loading amplitude and $2a$ is the arc-chord length of the crack.

A comparison of the numerical results obtained by the three BEM is shown in Figs. 9, 10, and 11 for the two crack-tips. Also here, a time-step $c_T \Delta t = a/20$ is used in the time-domain BEM. As expected, the differences

between the dynamic SIFs at both crack-tips increase as the curvature of the arc-shaped crack increases. The numerical results show that, except for the mode-II SIF at the crack-tip B, the time at which the peaks of the SIFs are achieved does not vary substantially with the curvature in the cases of $\alpha = 15^\circ$ and $\alpha = 30^\circ$. In the case of $\alpha = 60^\circ$, the peak $K_I$-factor at the crack-tip A is shifted to a smaller time instant. The results obtained by the time-domain BEM and the frequency-domain BEM agree very well, while the results obtained by the Laplace-domain BEM using 12 Laplace-transform parameters show some differences. The main difference lies in the tendency that the Laplace-domain BEM
<figure><img src="image_13.png" /><figcaption>**Fig. 11** Normalized mode-I and mode-II dynamic SIFs for an arc-shaped crack with $\alpha = 60^{\circ}$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12** Normalized mode-I and mode-II dynamic SIFs for an arc-shaped crack with $\alpha = 15^{\circ}$ by using three different time-steps</figcaption></figure>

with 12 Laplace-transform parameters cannot predict the peaks accurately.

In order to investigate the influence of the used number of the Laplace-transform parameters on the accuracy of the Laplace-domain BEM, numerical calculations have also been carried out for 24 Laplace-transform parameters. The results show that the agreement with other two BEM are improved now in the first part of the plot, i.e., for a dimensionless time below about 5.5. However, in the large time-range, oscillation of the numerical results appears which indicates an instability of the method.

The stability of the time-domain BEM is tested here again for three opening angles of the arc-shaped crack. The numerical results for the normalized dynamic SIFs are presented in Figs. 12, 13, and 14 versus the dimensionless time $tc_T/a$. For three quite different time-steps selected, no essential differences in the numerical results are observed, which confirms again that the present time-domain BEM is quite insensitive to the used time-steps. Moreover, the curvature of the crack-faces does not affect the stability behavior of the present time-domain BEM, at least for the investigated cases in this analysis.
<figure><img src="image_15.png" /><figcaption>Fig. 13 Normalized mode-I and mode-II dynamic SIFs for an arc-shaped crack with $\alpha = 30^\circ$ by using three different time-steps</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14 Normalized mode-I and mode-II dynamic SIFs for an arc-shaped crack with $\alpha = 60^\circ$ by using three different time-steps</figcaption></figure>

Mesh-sensitivity of the dynamic BEM is a critical issue for dynamic problems. For a straight crack in an anisotropic solid subjected to a tensile impact loading, a mesh-sensitivity study of the three presented BEM is carried out for different numbers and distributions of the boundary elements. The used material constants are given in Eq. (62). The corresponding numerical results for the normalized dynamic SIFs are presented in Figs. 15, 16, 17, 18, 19, 20, and 21. From this study the following conclusions can be drawn:

- The three BEM presented in this paper are quite insensitive to the used element-number.

- Ten quadratic elements are sufficient to obtain accurate numerical results.

- A uniform mesh is also adequate for the dynamic crack problem under consideration, although a non-uniform mesh is applied in most part of the analysis.

In addition, the following comments to the mesh-sensitivity of the dynamic BEM should be made:

- In the classical time-domain BEM, the spatial discretization and the time discretization are not independent of each other. This implies that the mesh-sensitivity and the time-step sensitivity cannot
<figure><img src="image_17.png" /><figcaption>**Fig. 15** Mesh-sensitivity of the frequency-domain BEM (uniform mesh)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 16** Mesh-sensitivity of the frequency-domain BEM</figcaption></figure>

be separated in the classical time-domain BEM. A proper choice of the element-size depends on the choice of the time-step, which is important to the quality and the stability of the BEM. A too small element-size (or too small time-step) may lead to an instability of the numerical scheme, while a too large element-size (or too large time-step) may cause a physically unrealistic large numerical damping of the results. A brief review on the subject can be found in reference [32]. To the best knowledge of the authors, yet there are no rigorous mathematical stability proofs for the classical time-domain BEM based on collocation methods both for spatial

and time discretizations. However, many previous including our own numerical experiences show that for isotropic elastic solids stable and accurate numerical results can be obtained by using $c_L \Delta t \le \Delta l \le 2c_L \Delta t$ [6], where $c_L$ is the longitudinal wave velocity, $\Delta l$ is the element-size, and $\Delta t$ is the used time-step. The situation for anisotropic elastic solids becomes even more tangled due to the directional dependence of the wave velocities. The present time-domain BEM by using Lubich's quadrature is in any case less sensitive to the used mesh-size in comparison to the classical time-domain BEM, since the present time-domain BEM is less sensitive to the used time-steps,
<figure><img src="image_19.png" /><figcaption>**Fig. 17** Mesh-sensitivity of the Laplace-domain BEM (uniform mesh)</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 18** Mesh-sensitivity of the Laplace-domain BEM</figcaption></figure>

which implies implicitly a less sensitivity of the method to the spatial mesh-size.

- In the frequency-domain BEM for time-harmonic analysis, the mesh-sensitivity is governed by the ratio of the nodal spacing and the wave-length. To reliably capture the sinusoidal waveform by using polynomial shape functions in the frequency-domain BEM, sufficient nodes within the wave-length are required. A commonly quoted and accepted rule of thumb recommends 8–10 nodes per wave-length (i.e., four or five quadratic elements) are needed. For large scale problems and high frequencies (short wave-lengths) this may cause some serious difficulties due to the

substantially large memory, storage and computing time, which may exceed the available computer resources. To overcome this difficulty, several advanced methods have been proposed in recent years. Here, we just mention the fast multipole method (FMM) [15, 18, 19], the wave boundary element method (WBEM) or the wave basis functions method [4, 16, 17] based on the partition of unity method (PUM) [13]. A review of some advanced computational methods for wave simulation in high frequency range has been presented by Bettess [4]. For the present transient dynamic crack analysis, numerical calculations for different frequencies are
<figure><img src="image_21.png" /><figcaption>**Fig. 19** Mesh-sensitivity of the time-domain BEM (uniform mesh)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 20** Mesh-sensitivity of the time-domain BEM</figcaption></figure>

needed. Therefore, it is impractical to use different element-numbers according to the above mentioned rule of thumb. Our numerical tests (see Figs. 15 and 16) show that a fixed element-number can still yield accurate results for transient dynamic crack analysis.

A mesh-sensitivity of the Laplace-domain BEM is also expected, which depends on the value of the Laplace-transform parameter. However to the best of the author's knowledge, yet no detailed analysis on the subject can be found in literature, which remains a future research topic. For the present transient dynamic crack analysis, it is sufficient to ignore

the effects of the dependence of the mesh on the Laplace-transform parameter (see Figs. 17 and 18).

## 7 Conclusions

In this paper, three different traction BEM for transient dynamic crack analysis in 2-D, homogeneous, generally anisotropic and linear elastic solids are presented and compared. Hypersingular TBIEs are applied in all three BEM. The arising hypersingular integrals in the TBIEs
<figure><img src="image_23.png" /><figcaption>**Fig. 21** Comparison of the computational costs of the three BEM ($t_L=19$ minutes $32$ seconds)</figcaption></figure>

are computed via a regularization technique based on a suitable variable change. This allows us to recast a hypersingular integral into a regular integral plus a strongly singular integral and another new hypersingular integral. Regular integrals are computed numerically by using standard Gaussian quadrature formula, while the strongly singular and the new hypersingular integrals are evaluated analytically.

In the frequency-domain and the Laplace-domain BEM, the boundary value problem is first solved in the transformed domain for discrete frequency and Laplace-transform parameters. Then, time-domain results are obtained by inverse Fourier-transform and inverse Laplace-transform. In the time-domain BEM, the quadrature formula of Lubich is applied for approximating the Riemann convolution integrals, which requires only Laplace-domain fundamental solutions instead of time-domain fundamental solutions. At crack-tips, quarter-point elements are adopted in all three BEM to describe the local behavior of the CODs at the crack-tips properly. Dynamic SIFs are obtained directly from the numerically computed CODs.

The implemented BEM presented in this paper are general and can be used for straight and curved cracks. Numerical examples for computing dynamic SIFs show that the implemented BEM are very accurate and robust, and they are valid even for isotropic material properties as a special case of the anisotropic materials. It is well known that the isotropic case is a mathematically degenerated case of the general anisotropic

BEM formulation, which may cause some numerical difficulties.

Regarding the computational efficiency or cost, the time-domain BEM is the fastest method while the Laplace-domain BEM with the Stehfest's inversion algorithm is the most expensive one, see Fig. 21. In order to obtain the numerical results for N time-steps, the system matrices have to be computed for

- $N/2 + 1$ Laplace-transform parameters in the time-domain BEM using the Lubich's quadrature,

- approximately $2N$ frequencies in the frequency-domain BEM using the fast Fourier inverse transform, and

- $n \times N$ Laplace-transform parameters in the Laplace-domain BEM using the Stehfest's inversion algorithm, with $n$ being the truncation limit ($n=12$ or $24$ in this analysis).

On the other hand, the Laplace-domain BEM with the Stehfest's inversion algorithm have two important advantages: it is easy to implement and it uses only a real Laplace-transform parameter.

Though the present crack analysis is shown in this paper for unbounded anisotropic solids, the extension of the three BEM to transient dynamic crack analysis in bounded anisotropic solids is straight-forward and the corresponding results will be reported in future.

**Acknowledgements** This work is supported by the German Research Foundation (DFG) under the project number ZH 15/5-1 and ZH 15/5-2, which is gratefully acknowledged.

## Appendix: Auxiliary functions

In Eq. (27), $\mu_m$ are the complex roots of the following characteristic equation

$$ \det[\Gamma_{ij}(1, \mu_m)] = 0. \qquad (63) $$

The matrices $P$ and $Q$ are determined by

$$ \Gamma(1, \mu_m)P = 0, \quad \text{or} \quad \Gamma_{ij}(1, \mu_m)P_{jm} = 0, $$

$$ Q = P^{-1} (B^{-1} + \bar{B}^{-1})^{-1}, \qquad (64) $$

where

$$ B = iPL^{-1}, $$

$$ \begin{aligned} L &= L_{im} = (C_{2ij1} + C_{2ij2}\mu_m) P_{jm} \\ &= -\frac{1}{\mu_m} (C_{1ij1} + C_{1ij2}\mu_m) P_{jm}, \end{aligned} \qquad (65) $$
The auxiliary functions $S_{ij}^m$, $T_{ij}^m$ and $U_{ij}^m$ in Eqs. (34)–(38) are given by

$$S_{ij}^m = C_{jklr} \gamma_{lr}^m n_k(\mathbf{x}) \eta_l, \quad (66)$$

$$T_{ij}^m = (C_{irk1} + \mu_m C_{irk2}) n_r(\xi) L_{jm} Q_{mk}, \quad (67)$$
(no sum over m),

$$U_{ij}^m = -C_{irkl} S_{kj}^m n_r(\xi) \eta_l. \quad (68)$$

## References

1. Achenbach JD (1980) Wave propagation in elastic solids. North Holland, Amsterdam

2. Albuquerque EL, Sollero P, Aliabadi MH (2002) The boundary element method applied to time dependent problems in anisotropic materials. Int J Solids Struct 39:1405–1422

3. Albuquerque EL, Sollero P, Fedelinski P (2003) Dual reciprocity boundary element method in Laplace domain applied to anisotropic dynamic problems. Comput Struct 81:1703–1713

4. Bettess P (2004) Short wave scattering: problems and techniques. Philos Trans R Soc Lond A 362:421–443

5. Dineva P, Rangelov T, Gross D (2005) BIEM for 2D steady-state problems in cracked anisotropic materials. Eng Anal Bound Elem 29:689–698

6. Domínguez J (1993) Boundary elements in dynamics. Computational Mechanics Publications, Southampton, UK

7. Eshelby JD, Read WT, Shockley W (1953) Anisotropic elasticity with applications to dislocation theory. Acta Metall 1:251–259

8. García F, Sáez A, Domínguez J (2004) Traction boundary elements for cracks in anisotropic solids. Eng Anal Bound Elem 28:667–676

9. García-Sánchez F, Sáez A, Domínguez J (2006) Two-dimensional time-harmonic BEM for cracked anisotropic solids. Eng Anal Bound Elem 30:88–99

10. Gaul L, Schanz M (1999) A comparative study of three boundary element approaches to calculate the transient response of viscoelastic solids with unbounded domains. Comput Methods Appl Mech Eng 179:111–123

11. Lubich C (1988) Convolution quadrature and discretized operational calculus, Part I. Numer Math 52:129–145

12. Lubich C (1988) Convolution quadrature and discretized operational calculus, Part II. Numer Math 52:413–425

13. Melenk JM, Babuska I (1996) The partition of unity finite element method. Basic theory and applications. Comput Methods Appl Mech Eng 136:289–314

14. Narayanan GV, Beskos DE (1982) Numerical operational methods for time-dependent linear problems. Int J Numer Methods Eng 18:1829–1854

15. Nishimura N (2002) Fast multipole accelerated boundary integral equation methods. Appl Mech Rev 55:299–324

16. Perrey-Debain E, Trevelyan J, Bettess P (2003) P-wave and S-wave decomposition in boundary integral equation for plane elastodynamic problems. Commun Numer Methods Eng 19:945–958

17. Perrey-Debain E, Trevelyan J, Bettess P (2004) Wave boundary elements: a theoretical overview presenting applications in scattering of short waves. Eng Anal Bound Elem 28:131–141

18. Rokhlin V (1985) Rapid solution of integral equations of classical potential theory. J Comput Phys 60:187–207

19. Rokhlin V (1990) Rapid solution of integral equations of scattering theory in two dimensions. Comput Phys 86:414–439

20. Schanz M (1999) A boundary element formulation in time domain for viscoelastic solids. Commun Numer Methods Eng 15:799–809

21. Schanz M (2001) Wave propagation in viscoelastic and poroelastic continua. Springer, Berlin Heidelberg New York

22. Sih GC, Paris PC, Irwin GR (1965) On cracks in rectilinearly anisotropic bodies. Int J Fract 1:189–203

23. Stehfest H (1996) Commun. Algorithm 368: numerical inversion of Laplace transform; an overview and recent developments. Comput Methods Appl Mech Eng 139:3–47

24. Tan A, Hirose S, Zhang Ch, Wang C-Y (2005) A time-domain BEM for transient wave scattering analysis by a crack in anisotropic solids. Eng Anal Bound Elem 29:610–623

25. Tan A, Hirose S, Zhang Ch (2005) A time-domain collocation-Galerkin BEM for transient dynamic crack analysis in anisotropic solids. Eng Anal Bound Elem 29:1025–1038

26. Thau SA, Lu TH (1971) Transient stress intensity factors for a finite crack in an elastic solid caused by a dilatational wave. Int J Solids Struct 7:731–750

27. Wang C-Y, Achenbach JD (1994) Elastodynamic fundamental solutions for anisotropic solids. Geophys J Int 118:384–392

28. Wang C-Y, Zhang Ch (2005) 3-D and 2-D dynamic Green's functions and time-domain BIEs for piezoelectric solids. Eng Anal Bound Elem 29:454–465

29. Zhang Ch (2000) Transient elastodynamic antiplane crack analysis of anisotropic solids. Int J Solids Struct 37:6107–6130

30. Zhang Ch (2002) A 2-D time-domain BIEM for dynamic analysis of cracked orthotropic solids. Comp Model Eng Sci 3:381–398

31. Zhang Ch (2002) A 2-D hypersingular time-domain traction BEM for transient elastodynamic crack analysis. Wave Motion 35:17–40

32. Zhang Ch, Savaidis A (2003) 3-D transient dynamic crack analysis by a novel time-domain BEM. Comput Model Eng Sci 4:603–618

33. Zhang Ch (2005) Transient dynamic crack analysis in anisotropic solids. In: Ivankovic A, Aliabadi MH (eds) Crack dynamics Chap. 3. WIT Press, Southampton, Boston, pp 71–120