1. Introduction

The problem of multi-crack interaction is of important interest as the interaction often results in an increase of the stress intensity factors (SIFs) which reduce the fracture strength of the structure. Although the interaction between multiple cracks under stactic loading conditions has attracted a large amount of attention from numerous research groups, e.g., ^[1,2,3,4,5] , there has been just a limited number of work devoted to the problem of multi-crack interaction under dynamic loading conditions, e.g., ^{[6,7,8,9,10,11]} . At issue is the lack of studies using versatile numerical methods, such as finite element method (FEM) and boundary element method (BEM), for investigating multi-crack dynamic interactions while these numerical methods are known for effectively dealing with a variety of practical crack problems (see, e.g., ^[12,13]). The objective of this work is to show that the dynamic interaction between multiple cracks in infinite domains can be accurately and effectively analyzed using a versatile numerical method called the symmetric-Galerkin boundary element method (SGBEM). This work studied cracks whose surfaces are symmetrically impact-loaded with a Heaviside pressure such as those on the inner surface of a large pressure vessel under sudden pressurization. The interaction between the cracks was assessed by investigating the dynamic SIFs (DSIFs) as the cracks approach each other.

The BEM is a powerful alternative to the well-established FEM. The key feature of the BEM is that only the boundary of the domain is discretized. For problems of wave scattering from a finite body in an infinite domain, such as those considered in this work, this means that artificial truncation of the domain is not required. For fracture analysis the important implications are that the singular stress field ahead of the crack is not approximated, and that remeshing a propagating crack is straightforward. A variant of the BEM employing a Galerkin approximation, SGBEM (e.g., ^[14]), possesses several important advantages for fracture applications: (a) the method employs displacement boundary integral equation (BIE) on the boundary part where displacement is prescribed and traction BIE on the boundary part where traction is known. As the name implies, this results in a symmetric coefficient matrix, and this remains true for fracture problems providing that the unknowns on the crack faces are the jumps in displacement; (b) the presence of both displacement and traction BIE enables fracture problems to be solved without artificial sub-domains which were employed in early BEM analyses (see ^[15] for example); and (c) unlike collocation methods, there is no smoothness requirement on the displacement (e.g., ^[16]) in order to evaluate the hypersingular integral; thus, standard continuous elements can be employed. The Galerkin approach can therefore easily exploit the highly effective modified quarterpoint (MQP) quadratic element ^[17] to accurately capture the crack tip behavior. As a result, the MQP element was employed to compute the DSIFs in the present work.

Dynamic analysis can be formulated either in the time domain (e.g., ^[18]) or frequency domain (e.g., ^[19]). In the frequency domain, analyses can be performed in the Fourier or Laplace spaces. If time solutions are needed following a frequency-domain analysis, frequency-time conversion can be done using fast Fourier transform (FFT, e.g., ^[20]) or Laplace transform (e.g., ^[19]) depending on the type of space used in the frequency-domain analysis. The dynamic interaction analysis in this work was conducted using the SGBEM in the Laplace domain (Laplace-SGBEM) in conjunction with the MQP element as this technique has shown to be accurate and effective in computing the DSIFs ^[19].

---

2. Symmetric Boundary Integral Formulations for Elastodynamics in the Laplace Domain

First, consider the body in Figure 1 without the crack. In the Laplace domain, the BIE for a source point P interior to a domain is given by

$\overline {\mathcal U} (P,s) \equiv {\bar u_k}(P,s) - \int_\Gamma {\left[{{{\bar U}_{kj}}(P,Q,s){{\bar t}_j}(Q,s) - {{\bar T}_{kj}}(P,Q,s){\mkern 1mu} {{\bar u}_j}(Q,s)} \right]} {\text{ d}}Q = 0$

(1)

where $Q$ is a filed point, $s$ is a Laplace parameter, ${\bar u}_j$ and ${\bar t}_j$ denote the transformed displacement and traction vectors, respectively, and $\Gamma$ is the boundary of the 2-D domain under consideration.

For $P$ off the boundary, the kernel tensors ${\bar U}_{kj}$ and ${\bar T}_{kj}$ are not singular and it is permissible to differentiate Eq. (1) with respect to $P$, yielding the transformed displacement gradient equation. It can be shown that the limits of the integrals in these equations as $P$ approaches the boundary exist. From now on, for $P \in \Gamma$, the BIEs are understood in this limiting sense. Substitution of the transformed displacement gradients into Hooke's law and then Cauchy's relation results in the following equation for surface traction:

$\overline {\mathcal T} (P,s) \equiv {\bar t_k}(P,s) - {n_\ell }(P)\int_\Gamma {\left[{{{\bar D}_{kj\ell }}(P,Q,s){\mkern 1mu} {{\bar t}_j}(Q,s) - {{\bar S}_{kj\ell }}(P,Q,s){{\bar u}_j}(Q,s)} \right]} {\text{d}}Q = 0$

(2)

In Eq. (2), $n_{\ell}$ is the outward normal vector to the boundary $\Gamma$. Expressions for the elastodynamic kernel tensors $\bar U_{kj}, \bar T_{kj}, \bar D_{kj\ell}$ and $\bar S_{kj\ell}$ in Eqs. (1) and (2) can be found in, e.g., ^[21].

The traction equation (2) is essential for treating crack geometries, and in the symmetric-Galerkin approach it is this equation that is employed on the crack surface.

The Galerkin boundary integral formulation is obtained by taking the shape functions $\psi_m$ employed in approximating the boundary tractions and displacements as weighting functions for Eqs. (1) and (2). For SGBEM, the displacement BIE (1) is employed on the boundary part $\Gamma_u$ where displacements are specified, while the traction BIE (2) is used on the boundary part $\Gamma_t$ where tractions are prescribed,

$\int_{{\Gamma _u}} {{\psi _m}} (P){\mkern 1mu} \bar {\mathcal U}(P,s){\mkern 1mu} {\text{d}}P = 0$

(3)

$\int_{{\Gamma _t}} {{\psi _m}} (P)\bar {\mathcal T}(P,s){\mkern 1mu} {\text{d}}P = 0$

(4)

As mentioned above, this arrangement results in a symmetric coefficient matrix. These formulas remain the same for fracture analysis, with the proviso that the unknowns on the crack are now the transformed displacement jump $\Delta {\bar u}_k$, and thus only one crack face is required to be discretized. In this work, a numerical implementation of Eqs. (3) and (4) is carried out with the standard quadratic element. Employing the parameter space $\xi\in\lbrack0, ~1\rbrack$, and defining $\xi_1=0$, $\xi_2=1/2$ and $\xi_3=1$, the quadratic shape functions are defined by $\psi_{\ell}(\xi_m) = \delta_{\ell m}$ and hence

$\psi_1(t) = (1-\xi)(1-2\xi)$

(5)

$\psi_2(t) = 4\xi(1-\xi)$

(6)

$\psi_3(t) = \xi(2\xi-1)$

(7)

One of the advantages of the frequency-domain analysis is that Eqs. (1) and (2) have a similar form as those in elastostatics. Details on the treatment of singular integrals in Eqs. (3) and (4) can be found in, for example, Ref. ^[19].

Microcracks in structures, such as those considered in this work, are preferred to be modeled as cracks in infinite/unbounded domains. Unlike some other numerical methods, the SGBEM can easily handle this type of problems without the need of using large artificial boundaries. The reader is referred to, for example, Ref. ^[22], for more details on SGBEM analysis of cracks in unbounded domains.

---

3. Fracture Analysis by the SGBEM for Laplace-Space Elastodynamics

Frequency domain analysis of cracks, such as the analysis of the DSIFs and dynamic T-stress, using the SGBEM in the Laplace space was introduced in ^[19]. In this section, the formulations for this type of analysis are briefly reviewed.

---

3.1. Formulations

If a crack of boundary $\Gamma_c$ is added to the domain of boundary $\Gamma = \Gamma_u \cup \Gamma_t$ considered in the previous section, the new total boundary becomes $\stackrel{*}{\Gamma} = \Gamma \cup \Gamma_c$ (see Figure 1). The crack is composed of two symmetrically loaded surfaces $\Gamma^+_c$ and $\Gamma^-_c$ which are initially coincident. Let $\stackrel{*}{\Gamma}_t = \Gamma_t + \Gamma^+_c$. In this case, the displacement and traction BIEs are written as

$\mathop {\mathcal U}\limits^* (P,s) = \bar {\mathcal U}(P,s) + \int_{\Gamma _c^ + } {{{\bar T}_{kj}}} (P,Q,s){\mkern 1mu} \Delta {\bar u_j}(Q,s){\mkern 1mu} {\text{ d}}Q = 0$

(8)

$\mathop {\mathcal T}\limits^* (P,s) = \bar {\mathcal T}(P,s) + n_\ell ^ + (P)\int_{\Gamma _c^ + } {{{\bar S}_{kj\ell }}} (P,Q,s){\mkern 1mu} \Delta {\bar u_j}(Q,s){\mkern 1mu} {\text{ d}}Q{\text{ }} = 0$

(9)

where $n_{\ell}^+$ is the outward normal vector to $\Gamma_c^+$, and as the transformed displacement jump vector $\Delta {\bar u}_j$ across the crack surfaces is used as the unknown on the crack, only one crack surface, e.g., $\Gamma_c^+$, needs to be discretized. It is well known that the traction BIE (9) is essential for treating crack geometries.

The use of $\Delta {\bar u}_j$ as the unknown on the crack as mentioned above is needed for obtaining a symmetric coefficient matrix. The symmetric-Galerkin formulation is given by

$\int_{{\Gamma _u}} {{\psi _m}} (P){\mkern 1mu} \mathop {\mathcal U}\limits^* (P,s){\text{d}}P = 0$

(10)

$\int_{{{\mathop \Gamma \limits^* }_t}} {{\psi _m}} (P){\mkern 1mu} \mathop {\mathcal T}\limits^* (P,s){\mkern 1mu} {\text{d}}P = 0$

(11)

---

3.2. Dynamic Stress Intensity Factors

For stationary cracks in elastic materials (as those considered in this work), the DSIFs under plane strain situations, and thereby their transforms can be determined from the asymptotic expansion for the displacement field in the vicinity of a crack tip as follows:

$\begin{gathered} {{\bar K}_{\text{I}}}(s) = \frac{\mu }{{4(1 - \nu )}}\mathop {\lim }\limits_{r \to 0} \sqrt {\frac{{2\pi }}{r}} n\Delta {{\bar u}_n}(s) \hfill \\ {{\bar K}_{{\text{II}}}}(s) = \frac{\mu }{{4(1 - \nu )}}\mathop {\lim }\limits_{r \to 0} \sqrt {\frac{{2\pi }}{r}} n\Delta {{\bar u}_t}(s) \hfill \\ \end{gathered}$

(12)

where $\mu$ and $\nu$ are the shear modulus and Poison's ratio of the material, $\Delta {\bar u}_n$ and $\Delta {\bar u}_t$ are, respectively, the normal and tangential components of the transformed crack displacement jump vector, and $r$ is the distance to the crack tip.

In both finite and boundary element modeling of discrete cracks, the standard approach consists of incorporating the quarter-point (QP) element ^[23,24] to improve the accuracy of the SIF calculations (e.g., ^[25,26]). However, as discussed in ^[27], the standard QP element in general fails to satisfy a constraint that exists in the series expansion of the crack opening displacement at the tip of an arbitrary crack geometry (see also ^[28]). This has led to the development of an improved MQP element ^[17]. It was demonstrated in ^[17,29] that the accuracy of the computed SIFs can be significantly improved by incorporating this MQP element into the SGBEM. Thus, the MQP element is employed to determine the DSIFs in this work.

By using the MQP shape functions in Eq. (12), the transformed DSIFs are accurately given by

$\begin{gathered} {{\bar K}_{\text{I}}}(s) = \frac{\mu }{{12(1 - \nu )}}\sqrt {\frac{{2\pi }}{L}} \left( {8\Delta \bar u_n^{(2)} - \Delta \bar u_n^{(3)}} \right) \hfill \\ {{\bar K}_{{\text{II}}}}(s) = \frac{\mu }{{12(1 - \nu )}}\sqrt {\frac{{2\pi }}{L}} \left( {8\Delta \bar u_t^{(2)} - \Delta \bar u_t^{(3)}} \right) \hfill \\ \end{gathered}$

(13)

where $L$ is the distance between the two end-nodes of the crack-tip element, and the superscripts (2) and (3) denote its quarter-point node and non-tip end-node, respectively.

---

4. Frequency Domain Analysis using SGBEM and Fast Laplace Inverse Transform

Let $\bar{P}(s)$ and $\bar{F}(s)$ respectively be the Laplace transforms of load $P(t)$ and time response $F(t)$, where $t$ is time, from a frequency analyis of a dynamic system. It is known that the output $\bar{F}(s)$ is related to the input $\bar{P}(s)$ as

$\bar{F}(s) = \bar{H}(s)\,\bar{P}(s)$

(14)

where $\bar{H}(s)$ is called the frequency response which is the response of the system due to a unit harmonic load $e^{ist}$ where $i$ is the imaginary unit.

Figure 2 depicts a model for obtaining time responses of a system from frequency response analysis. First, the SGBEM is employed to obtain the frequency response $\bar{H}(s)$ of the system under the unit harmonic load ($e^{ist}$). In the meantime, the Laplace transform (LT) is used to convert $P(t)$ to $\bar{P}(s)$. Relation (14) is then applied to obtain the dynamic response $\bar{F}(s)$. Finally, the fast Laplace inverse transform (FLIT) proposed by Durbin ^[30] is utilized to transform $\bar{F}(s)$ into the desired time history $F(t)$.

A procedure for obtaining time responses by using the FLIT can be summarized as follows:

(a) Determine a value for the duration of analysis $T_f$. Note that the frequency resolution $\Delta f$ ($f=\omega /(2\pi)$) is related to $T_f$ as $\Delta f = 1/T_f$. The frequency resolution should be small enough to minimize the loss of frequency information.

(b) Determine a value for the real part $a_o$ of the Laplace parameters such that

$a_oT_f = 5 \ldots 10$

(15)

(c) Perform SGBEM analysis for a series of ($N_o+1$) Laplace parameters, namely $s_k = a_o + i\frac{2\pi k}{T_f}$ where $k=0 \ldots N_o$ to obtain the frequency response $H(s_k)$. Note that the very first sample ($k=0$) of the series corresponds to the static sample and $N_o$ does not have to be a power-of-two number as in the cases of FFT;

(e)] Perform the Laplace transform for the loading function $P(t)$ to obtain

$\bar P(s) = \int_0^\infty P (t){\mkern 1mu} {e^{ - st}}{\mkern 1mu} {\text{d}}t$

(16)

(f) Calculate $\bar{F}(s_k)=\bar{H}(s_k)\, \bar{P}(s_k)$.

(g)] Perform the FLIT for $\bar{F}(s_k)$ to obtain the desired time history as follows ^[30]:

$F(t_j) = \frac{2}{T_f}\,e^{a_o j\Delta t}\left[-\frac{1}{2}\Re\{\bar{F}(a_o)\} + \Re\left\{\sum_{k=0}^{N_o}\bar{F}(s_k)\left(\cos kj\frac{2\pi}{N} + i\sin kj\frac{2\pi}{N}\right)\right\}\right]$

(17)

In this equation, $t_j = j\Delta t = j\frac{T_f}{N}$ where $j = 0 \ldots N-1$, $N = k_\text{if} N_o$ and $k_\text{if} = 1, 2, \ldots$. Here, $k_\text{if} > 1$ is used as an interpolation factor to increase the resolution of the time-history curve by reducing the time increment $\Delta t$.

---

5. Numerical Examples

In this section, three numerical examples of the dynamic interaction between multiple cracks are investigated via the DSIFs computed using the Laplace-SGBEM technique outlined in Sections 2 through 4. All the cracks are situated in infinite domains. The crack surfaces are symmetrically loaded by a pressure $p(t)=p_oH(t)$ (see Figure 3) where $p_o$ is the magnitude of the pressure and $H(t)$ is the Heaviside step function numerically defined as (see also Figure 4 where $T_f$ is the final time of analysis).

$\begin{gathered} H(t) = 0{\text{if }}t < 0 \hfill \\ H(t) = 1/2{\text{if }}t = 0 \hfill \\ H(t) = 1{\text{if }}t > 0 \hfill \\ \end{gathered}$

(18)

The material properties of the domains are $\mu=76.923$ GPa, $\nu=0.3$ and mass density $\rho=5, 000$ kg/m$^3$. In this work, the analysis of the dynamic interaction between the cracks is based on the normalized mode-I and mode-II DSIFs defined as

${F_{\text{I}}} = \frac{{{K_{\text{I}}}}}{{{p_o}\sqrt {\pi a} }}, \qquad {F_{{\text{II}}}} = \frac{{{K_{{\text{II}}}}}}{{{p_o}\sqrt {\pi a} }}$

(19)

which are functions of normalized time $c_\text{s}t/a$ where $a$ is the half length of the crack under consideration and $c_\text{s}=\sqrt{\mu/\rho}$ is the shear wave velocity of the material.

For the SGBEM analysis of all the three problems examined in this section, 10 equal-length quadratic elements were used to discretize each of the cracks. Among these quadratic elements, MQP elements were employed as the crack-tip elements for the purpose of obtaining highly accurate DSIFs.

For the FLIT employed in these three problems, a convergence study from the frequency-domain analysis resulted in the following selection: $c_\text{s}T_f/a=120$, $a_oT_f=5$ and $N_o=100$.

---

5.1. Three Collinear Cracks

The first example deals with the dynamic interaction of three collinear cracks, each of length $2a$ and separated by a gap $c$ as shown in Figure 5. As the cracks are collinear and perpendicular to the applied load, the sliding mode (mode II) is absent in this case ($F_\text{II}=0$).

Figures 6(a) and 6(b) show the histories of the normalized mode-I DSIF for two cases of the gap between the cracks, namely, $c/a=1$ and $c/a=0.2$. Due to the symmetry of the problem, the mode-I DSIF results at crack tips $A_1$, $A_2$ and $A_3$ are the same as at tips $B_3$, $B_2$ and $B_1$, respectively. In these figures, $F_\text{I}^\text{S}$ denotes the corresponding static mode-I SIF resulted from the static pressure $p=p_o$.

It can be seen from Figures 6(a) and 6(b) that the curves of the mode-I DSIF attain their peak within a normalized time of 10 after the crack surfaces are suddenly loaded. The largest peak occurs at the tips of the middle crack while the lowest peak occurs at crack tips $A_1$ and $B_3$. Also, as the gap between the cracks (represented by ratio $c/a$) decreases, the maximum values of $F_\text{I}$ at the tips increase as expected.

Subsequent travels of the diffracted waves between the crack tips result in decaying oscillations of $F_\text{I}$ about $F_\text{I}^\text{S}$.

---

5.2. Three Parallel Cracks

This example considers three parallel cracks of the same length $2a$ (see Figure 7). The cracks are aligned in the direction of the applied load and the distance between two neighboring cracks is denoted as $c$.

Normalized mode-I DSIF histories are shown in Figures 8(a) and 8(b) for $c/a=1$ and $c/a=0.2$, respectively. While mode-II is absent in the case of the collinear cracks studied in the previous example, it does exist for the top and bottom cracks in this example (see Figures 9(a) and 9(b)) due to the presence of the middle crack which causes shearing deformations of the top and bottom cracks.

For both fracture modes, the maximum (absolute) value of the DSIFs occur at the top and bottom cracks. As the cracks approach each other, the maximum absolute value of the mode-II DSIF for the top and bottom cracks increases while the maximum mode-I DSIFs of all the cracks decreases. The reduction of $F_\text{I}$ can be explained by the fact that, as the distance $c$ between the cracks is small enough, the opening of any crack is obstructed by the opening of the neighboring crack(s). Also, the decrease of $F_\text{I}$ for the middle crack is more pronounced as the opening of this crack is hindered by the opening of both the top and bottom cracks. For this example, the decaying oscillatory behavior of both $F_\text{I}$ and $F_\text{II}$ about their corresponding static value $F_\text{I}^\text{S}$ and $F_\text{I}^\text{S}$ can mainly be attributed to the subsequent travels of the reflected waves between the two neighboring cracks.

---

5.3. Four Radial Cracks

In this last example, four equally spaced radial cracks (two horizontal cracks of length $2a$ and two vertical cracks of length $2b=a$) as depicted in Figure 10 were studied. Figures 11(a) and 11(b) show the effect of $c/a$ on $F_\text{I}$ for the horizontal cracks while Figures 12(a) and 12(b) show the said effect on $F_\text{I}$ for the vertical cracks. Due to the symmetry of the problem, there are only four distinct histories of $F_\text{I}$ for a given value of $c/a$ at the four pairs of crack tips ($A_1, A_2$), ($A_3, A_4$), ($B_1, B_2$) and ($B_3, B_4$).

As in the previous examples, all the DSIF curves exhibit the decaying oscillatory behavior about their corresponding static SIF value. As the gap $c$ is not small enough, the maximum values of $F_\text{I}$ occur at tips $A$'s where the corresponding static values are also larger than those at tips $B$'s. However, the peaks of $F_\text{I}$ as well as their corresponding static values, switch from tips $A$'s to tips $B$'s as the ratio $c/a$ decreases. This is expected as the interaction effect is more pronounced at crack tips that are close enough to each other.

---

6. Conclusion

A study of dynamic interaction between multiple cracks whose surfaces are symmetrically impactloaded in infinite domains using the SGBEM for elastodynamics in the Laplace space was presented in this work. Four major concluding remarks can be drawn from this investigation as follows:

(a) Although there is no similar work in the literature which could be used to validate the numerical results reported in the present study, the numerical esults for the DSIFs obtained here can be justified as they oscillate about their corresponding static value as time elapses. This DSIF behavior can also be found in other studies where cracks are situated in infinite domains, e.g., ^[6,9,10,11].

(b) For cracks loaded with pressure in the form of the Heaviside step function, the peak of mode- I (and mode-II, if any) DSIF for a given crack tip occurs almost immediately after the load is applied.

(c) Due to crack-tip interaction, the maximum (absolute) value of the DSIFs at a crack tip can increase or decrease as it is approached by a neighboring crack tip.

(d) Following the sudden loading on the crack surfaces, subsequent travels of the diffracted waves back and forth between the crack tips (and the reflected waves between the cracks, if any) result in damped-like oscillatory behavior of the DSIFs about their corresponding static value.

---

Conflict of Interest

The author declares no conflict of interest in this paper.

---