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

Three different boundary element methods (BEM) for transient dynamic crack analysis in twodimensional (2-D), homogeneous, anisotropic and linear elastic solids are presented. Hypersingular traction boundary integral equations (BIEs) in frequencydomain, 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 Laplacetransform 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. 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 Ch. Zhang(B) Department of Civil Engineering, University of Siegen, Paul-Bonatz-Str. 9-11, 57076 Siegen, Germany e-mail: c.zhang@uni-siegen.de


A comparative study of three BEM for transient dynamic crack analysis of 2-D anisotropic solids
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 twodimensional (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-openingdisplacements (CODs) are the fundamental unknown 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 Laplacedomain [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 Laplacedomain 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.

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] Hooke's law the initial conditions and the traction boundary conditions on the crack-faces In Eqs. (1)-(4), σ ij and u i denote the stress and the displacement components, C ijkl is the fourth order elasticity tensor, p i (x, t) is the traction vector, n j is the outward unit normal vector, c = + c + − c represents the upper and the lower crack-faces, ρ is the mass density, and p * i (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 where x = (x 1 , x 2 ) and ξ = (ξ 1 , ξ 2 ) represent the field and the source points, p * ij are the traction fundamental solutions, u j (x, t) are the CODs defined by and an asterisk * denotes Riemann convolution which is defined by The traction fundamental solution p * ij is related to the displacement fundamental solution u * ij by Substituting Eq. (5) into Hooke's law (2) we obtain an integral representation formula for the traction components as where and n k (ξ ) represents the outward unit normal vector to the boundary at the collocation point. By taking the limiting process ξ → + c , the following time-domain traction BIEs are obtained where = 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).

Elastodynamic fundamental solutions
Unlike in the isotropic case, elastodynamic fundamental solutions for homogeneous, anisotropic and linear elastic solids cannot be given in closed-forms and they are much more complicated. For 2-D case, fundamental 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 solutions 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] where η = (η 1 , η 2 ) is the wave propagation vector, ω is the circular frequency, γ m ij = V im V jm (no sum over m) is the projection operator, with V im and ρc 2 m being the eigenvectors and the eigenvalues of the Christoffel tensor ij ij (η 1 , η 2 ) − ρc 2 m δ ij V jm = 0, (no sum over m), (13) in which δ ij denotes the Kronecker delta, c m represents the phase velocities of elastic waves, and the Christoffel tensor is defined by The function (ζ ) is given by where ci and si denote the cosine and the sine integrals which are defined by with ζ being real and − indicating the Cauchy principal value integral. Two-dimensional Laplace-domain displacement fundamental 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 where s is the Laplace-transform parameter, with z being complex and Ei being the exponential integral defined by It can be easily shown that both functions and have a logarithmic singularity. To deal with this singularity it is advantageous to split the fundamental solutions into a singular static part and a regular dynamic part as [27,28] 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 (22) while the regular dynamic part can be written as with the regular continuous functions The static displacement fundamental solution (22) can be reduced to the following explicit expression [7][8][9] where are the complex counterpart of the integration and the collocation points, the matrices P jm , Q mi and the complex constants µ 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 where In Eqs. (33)-(38), (ζ ) = d (ζ )/dζ , (z) = d (z)/dz, and the auxiliary functions L jm , S m ij , T m ij and U m ij are given in the Appendix.

Numerical implementation of the BEM
To solve the hypersingular TBIEs (11) as well as their counterparts in the frequency-domain and the Laplacetransformed 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 1 -continuity requirement of the CODS in the hypersingular TBIEs.

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 0 m ) 2 when the integration point coincides with the collocation point. After the discretization, the hypersingular integrals to be computed have the following expression where e is the boundary element under consideration, φ denotes the quadratic shape-function and n(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 the Jacobian of this transformation is then given by where the following relations have been used (see Fig. 1). By using Eqs. (40) and (41), the hypersingular integral (39) can be rewritten as Considering the shape function φ as a function of the complex variable χ m , and using the first two terms of its Taylor-series expansion at χ m = 0, i.e., the integral I m can be recast into 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

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 Laplacetransform 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 where f (s) is the Laplace-transform of f (t), and 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.

Time-domain BEM
In the time-domain BEM, the Riemann convolution integral is approximated by the quadrature formula of Lubich [11,12], which is given by where the time-interval t is divided into N equal time-steps t, and the weights ω n ( t) are defined by In Eq. (51), -g(·) is the Laplace-transform of g(t), ζ m = re 2πim/N , and -r = 1/2N , with being the numerical error in computing the Laplace-transform g(·).
It should remarked here that a backward differential formula (BDF) of the second order is used for δ(ζ m ), and the error in computing ω n ( t) is of the order O( √ ) [11,12]. In this analysis, different -values between 10 −6 and 10 −12 have been tested to verify the effects of on the numerical accuracy. Our numerical tests have shown that in this -range there are no notable influences of 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 timedomain 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 A(s m ), can be computed by using the following equation According to Eq. (51), the system matrix at the nth timestep A n can be obtained as where s m = δ(ζ m )/ t. Finally, a system of linear algebraic equations for the discrete CODs can be obtained as By invoking the zero initial conditions (3), Eq. (54) leads to the following explicit time-stepping scheme [29,30,33] for computing the unknown CODs at the nth time-step.
In Eq. (55), A 0 −1 is the inverse of the system matrix A 0 at the time-step n = 0.

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 √ 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] where r and θ 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 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 SFIs and the CODs where

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α 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α/30, while the first three elements from the crack-center have an angle 0.50α, 0.30α and 0.15α, 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  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.

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 The selected material constants correspond to a Graphite-epoxy composite [33]. The same BEM mesh as in the isotropic case has been applied. For a tensile impact crack-face loading, the normalized mode-I and mode-II dynamic SIFs versus the dimensionless time t c 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 = √ C 66 /ρ, and a timestep c T 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.

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 opening angle 2α as shown in Fig. 3. A tensile impact crack-face loading in the radial direction in the form of 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 σ 0 √ π a, where σ 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 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 α = 15 • and α = 30 • . In the case of α = 60 • , 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  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 nonuniform 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 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 t ≤ l ≤ 2c L t [6], where c L is the longitudinal wave velocity, l is the element-size, and 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 timedomain BEM is less sensitive to the used time-steps, 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   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 Laplacetransform 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 Riemman 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 cracktips 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 timedomain BEM using the Lubich's quadrature, -approximately 2N frequencies in the frequencydomain BEM using the fast Fourier inverse transform, and -n × N Laplace-transform parameters in the Laplacedomain 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.
The auxiliary functions S m ij , T m ij and U m ij in Eqs. (34)-(38) are given by S m ij = C jklr γ m ir n k (x)η l , T m ij = (C irk1 + µ m C irk2 ) n r (ξ )L jm Q mk , (no sum over m), U m ij = −C irkl S m kj n r (ξ )η l .