# Convolution Quadrature Time-domain Boundary Element Method and Acceleration by the Fast Multipole Method in 2-D Viscoelastic Wave Propagation

Takahiro Saitoh, Sohichi Hirose, Takuo Fukui

► To cite this version:

Takahiro Saitoh, Sohichi Hirose, Takuo Fukui. Convolution Quadrature Time-domain Boundary Element Method and Acceleration by the Fast Multipole Method in 2-D Viscoelastic Wave Propagation. Theoretical and Applied Mechanics Japan, 2009, 57, pp.385-393. 10.11345/nctam.57.385. hal-01581489

HAL Id: hal-01581489

https://hal.science/hal-01581489

Submitted on 4 Sep 2017

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
# Convolution Quadrature Time-domain Boundary Element Method and Acceleration by the Fast Multipole Method in 2-D Viscoelastic Wave Propagation

Takahiro SAITOH*, Sohichi HIROSE** and Takuo FUKUI*

This paper presents a new time-domain boundary element method (BEM) in 2-D viscoelastic wave propagation. The conventional time-domain BEM approach cannot be used in general, since it is difficult to obtain closed time-domain fundamental solutions in viscoelastic wave propagation. To overcome the difficulty, in this paper, the convolution quadrature method (CQM) developed by Lubich is applied to 2-D viscoelastic wave propagation. In the proposed method, the convolution integral is numerically approximated by quadrature formulas, whose weights are computed by using the Laplace transform of the fundamental solution in 2-D viscoelastodynamics. In addition, the fast multipole method (FMM) is applied to improve computational efficiency.

## 1. INTRODUCTION

Since the boundary element method (BEM) is known as a suitable numerical approach for wave analysis, time-domain transient problems have been solved by many researchers using the BEM formulated by Mansur and Brebbia(1), and Hirose(2). However, it is known that the use of a direct time-domain BEM sometimes causes numerical instability in time-stepping procedures.

Recently, the convolution quadrature boundary element method (CQ-BEM) has been researched by Schanz(3) and Abreu et al.(4) for transient analysis, following the proposal of the convolution quadrature method (CQM) by Lubich(5) in 1988. In the formulation of the CQ-BEM, convolution integrals in the time-domain boundary element formulation are numerically approximated by the CQM. Application of the CQM to the time-domain BEM not only improves the numerical stability of the time-domain BEM, but also expands its range of engineering applications, because the CQ-BEM utilizes the Laplace-domain, and not time-domain, fundamental solutions derived from the governing equations of a problem. In fact, the conventional time-domain BEM approach cannot be used for viscoelastic problems because no time-domain fundamental solution is known. Therefore, using the CQ-BEM is particularly helpful for viscoelastic problems.

Thus, the CQ-BEM is elegant. However, it shares a critical disadvantage with the conventional time-domain BEM: these methods cannot solve large-scale problems because of low computational efficiency. If this shortcoming is resolved, the CQ-BEM should become a more powerful technique for transient analysis of wave propagation problems.

<footnote>* Department of Nuclear Power and Safety Energy Engineering, University of Fukui, Fukui</footnote>

<footnote>** Department of Mechanical and Environmental Informatics, Tokyo Institute of Technology, Tokyo</footnote>
To overcome the difficulty, in this paper we propose an efficient CQ-BEM using the fast multipole method (FMM)⁽⁶⁾⁽⁷⁾ developed by Rokhlin and Greengard. In the following, after a review of the CQM, the conventional time-domain BEM and CQ-BEM formulations for 2-D viscoelastic wave propagation are presented. The CQ-BEM using the FMM is then discussed. Finally, numerical examples are shown to validate the proposed method.

## 2. CONVOLUTION QUADRATURE METHOD (CQM)

In this section, the convolution quadrature method (CQM) is described briefly. The convolution quadrature method (CQM), first proposed by Lubich, approximates the convolution $f * g(t)$ by a discrete convolution using the Laplace transform of the time-dependent function $f(t - \tau)$. In general, the convolution integral is defined as follows:

$$f * g(t) = \int_{0}^{t} f(t - \tau)g(\tau)d\tau, \quad t \ge 0 \qquad (1)$$

where * denotes the convolution. The convolution integral defined by Eq.(1) is approximated by the CQM as follows:

$$f * g(n\Delta t) \approx \sum_{j=0}^{n} \omega_{n-j}(\Delta t)g(j\Delta t), \quad n = 0, 1, \dots, N \qquad (2)$$

where time $t$ was divided into $N$ equal steps $\Delta t$. Moreover, $\omega_j(\Delta t)$ denotes the quadrature weights, which are determined by the coefficients of the following power series with complex variable $z$, namely,

$$F\left(\frac{\delta(\zeta)}{\Delta t}\right) = \sum_{n=0}^{\infty} \omega_n(\Delta t) z^n. \qquad (3)$$

In Eq.(3), $F$ is the Laplace transform of the time-dependent function $f$ and $\delta(\zeta)$ is the quotient of the generating polynomials of a linear multistep method. The power series defined in Eq.(3) can be calculated by Cauchy's integral formula. Considering a polar coordinate transformation, Cauchy's integral is approximated by a trapezoidal rule with $L$ equal steps $2\pi/L$ as follows:

$$\omega_n(\Delta t) = \frac{1}{2\pi i} \int_{|\zeta|=R} F\left(\frac{\delta(\zeta)}{\Delta t}\right) \zeta^{-n-1} d\zeta \approx \frac{R^{-n}}{L} \sum_{l=0}^{L-1} F\left(\frac{\delta(\zeta_l)}{\Delta t}\right) e^{-\frac{2\pi i n l}{L}} \qquad (4)$$

where $i$ is the imaginary unit and $\zeta_l$ is given by $\zeta_l = Re^{2\pi i l/L}$. In addition, $R$ is the radius of a circle in the domain of analyticity of $F$.

## 3. CQ-BEM FORMULATION FOR 2-D VISCOELASTIC WAVE PROPAGATION

We consider 2-D viscoelastic wave scattering by scatterer $\bar{D}$ in an exterior viscoelastic medium $D$. When the incident wave $u^{in}$ hits the boundary surface $S$ of scatterer $\bar{D}$, scattered waves are generated by the interaction between scatterer $\bar{D}$ and the incident wave. Assuming zero initial conditions, i.e., $u_i(x, t = 0) = 0$ and $\partial u_i(x, t = 0)/\partial t = 0$, the governing equations and boundary conditions are written as

$$\mu(t) * \dot{u}_{i,jj}(\mathbf{x}, t) + \left( K(t) + \frac{1}{3}\mu(t) \right) * \dot{u}_{j,ij}(\mathbf{x}, t) = \rho \ddot{u}_i(\mathbf{x}, t) \quad \text{in } D \qquad (5)$$

$$u_i = \bar{u}_i \quad \text{on } S_1, \quad t_i = \bar{t}_i \quad \text{on } S_2, \quad S_2 = S \setminus S_1. \qquad (6)$$
In Eq.(5) and Eq.(6), $\mu(t)$ and $K(t)$ are the shear relaxation and bulk relaxation functions, respectively, which satisfy the axiom of nonretroactivity, i.e., $\mu(t) = K(t) = 0, -\infty < t < 0$. In addition, $\rho$ is the density of elastic medium $D$. $u_i(x, t)$ and $t_i(x, t)$ are the displacement and traction components of the displacement and traction vector in Cartesian coordinate system at $x$ and time $t$, respectively. Moreover, $\bar{u}_i$ and $\bar{t}_i$ are the given boundary conditions, $(\cdot)_i$ denotes space derivative, and $(\cdot)$ indicates the partial derivative with respect to time $t$. Summation over repeated subscripts is implied and indices range from 1 to 2 in this paper.

In general, the time-domain fundamental solution $U_{ij}(x, y, t)$ and its double layer kernel $T_{ij}(x, y, t)$ cannot be obtained analytically for a linear viscoelastic problem as mentioned above. However, assuming the existence of an explicit form of the time-domain fundamental solutions, the time-domain boundary integral equation in 2-D viscoelastodynamics can be expressed by

$$C_{ij}(\mathbf{x})u_i(\mathbf{x}, t) = u_i^{\text{in}}(\mathbf{x}, t) + \int_S U_{ij}(\mathbf{x}, \mathbf{y}, t) * t_j(\mathbf{y}, t)dS_y - \int_S T_{ij}(\mathbf{x}, \mathbf{y}, t) * u_j(\mathbf{y}, t)dS_y \quad (7)$$

where $C_{ij}$ is the free term. Normally, the boundary integral equation is discretized by using the appropriate interpolation functions for the unknown values and solved by a time-stepping algorithm. The boundary integral equation (7), however, cannot be solved using such a scheme because there is no explicit time-domain fundamental solution. Therefore, the CQM is applied to Eq.(7) to overcome the difficulty.

If we discretize the boundary surface $S$ into $M$ elements by a piecewise constant approximation of the unknown displacement $u_i$ and traction $t_i$, taking the limit of $x \in D \to x \in S$, and finally using CQM for the convolutions in Eq.(7), we arrive at the following discretized boundary integral equation for time increment $\Delta t$ and $n$ steps as follows:

$$\frac{1}{2} u_i(\mathbf{x}, n\Delta t) = u_i^{\text{in}}(\mathbf{x}, n\Delta t) + \sum_{\alpha=1}^{M} \sum_{k=1}^{n} \left[ A_{ij}^{n-k}(\mathbf{x}, \mathbf{y}^{\alpha}) t_j^{\alpha}(k\Delta t) - B_{ij}^{n-k}(\mathbf{x}, \mathbf{y}^{\alpha}) u_j^{\alpha}(k\Delta t) \right]. \quad (8)$$

Here, $A_i^m$ and $B_i^m$ are influence functions defined by

$$A_{ij}^{m}(\mathbf{x}, \mathbf{y}) = \frac{\mathcal{R}^{-m}}{L} \sum_{l=0}^{L-1} \int_{S} \hat{U}_{ij}(\mathbf{x}, \mathbf{y}, s_l) e^{-\frac{2\pi i m l}{L}} dS_y \quad (9)$$

$$B_{ij}^{m}(\mathbf{x}, \mathbf{y}) = \frac{\mathcal{R}^{-m}}{L} \sum_{l=0}^{L-1} \int_{S} \hat{T}_{ij}(\mathbf{x}, \mathbf{y}, s_l) e^{-\frac{2\pi i m l}{L}} dS_y \quad (10)$$

where $s_l$ is given by $s_l = \delta(\zeta_l)/\Delta t$. In addition, $\hat{U}_{ij}(x, y, s)$ and $\hat{T}_{ij}(x, y, s)$ are the Laplace-domain fundamental solutions, which can be analytically obtained as follows:

$$\hat{U}_{ij}(\mathbf{x}, \mathbf{y}, s) = \frac{1}{2\pi\mu^*(s)} \left[ K_0(s_T r)\delta_{ij} - \frac{1}{s_T^2}\{K_0(s_T r) - K_0(s_L r)\}_{,ij} \right] \quad (11)$$

$$\hat{T}_{ij}(\mathbf{x}, \mathbf{y}, s) = n_j(\mathbf{y})\rho(c_L^{*2} - 2c_T^{*2})\hat{U}_{ik,j}(\mathbf{x}, \mathbf{y}, s) + \rho c_T^{*2} (\hat{U}_{ij,k}(\mathbf{x}, \mathbf{y}, s) + \hat{U}_{ik,j}(\mathbf{x}, \mathbf{y}, s)) n_k(\mathbf{y}) \quad (12)$$

where $c_T^*$ and $c_L^*$ are the transverse and longitudinal wave velocities given by $c_T^* = \sqrt{\mu^*(s)/\rho}$ and $c_L^* = \sqrt{(K^*(s) + (4/3)\mu^*(s))/\rho}$, respectively. In addition, $\mu^*(s)$ and $K^*(s)$ are defined by $\mu^*(s) = s\hat{\mu}(s)$ and $K^*(s) = s\hat{K}(s)$ using the Laplace-transformed relaxation functions $\hat{\mu}(s)$ and $\hat{K}(s)$, respectively. Also, $\delta_{ij}$ is the Kronecker delta, $r$ is given by $r = |x - y|$ and $K_n$ is the modified Bessel function of the second kind. Note that $s_L$ and $s_T$ are defined by $s_L = \delta(\zeta)/(c_L^*\Delta t)$ and $s_T = \delta(\zeta)/(c_T^*\Delta t)$ due to the
<figure><img src="image_2.png" /><figcaption>Figure 1. Multipole expansion and Graf's addition theorem.</figcaption></figure>

simple expression in Eq.(11). The calculations of Eq.(9) are identical to the discrete Fourier transform. Therefore, Eq.(9) can be evaluated rapidly by means of the fast Fourier transform (FFT) algorithm. Arranging Eq.(8) according to the boundary conditions, we obtain

$$ 
\begin{aligned}
& \frac{1}{2} u_i(\mathbf{x}, n\Delta t) + \sum_{\alpha=1}^{M} [B_{ij}^{0}(\mathbf{x}, \mathbf{y}^{\alpha})u_j^{\alpha}(n\Delta t) - A_{ij}^{0}(\mathbf{x}, \mathbf{y}^{\alpha})t_j^{\alpha}(n\Delta t)] \\
& = u_i^{\text{in}}(\mathbf{x}, n\Delta t) + \sum_{\alpha=1}^{M} \sum_{k=1}^{n-1} [A_{ij}^{n-k}(\mathbf{x}, \mathbf{y}^{\alpha})t_j^{\alpha}(k\Delta t) - B_{ij}^{n-k}(\mathbf{x}, \mathbf{y}^{\alpha})u_j^{\alpha}(k\Delta t)].
\end{aligned} 
\quad (13) $$

For the n-th time step, all the quantities on the right-hand side of Eq.(13) are known. Therefore, the unknown values $u_i^{\alpha}(n\Delta t)$ and $t_i^{\alpha}(n\Delta t)$ in the n-th time step can be obtained by solving Eq.(13).

Unfortunately, we cannot solve a large-scale problem with a large number of boundary elements $M$ and time steps $N$ by the CQ-BEM because the required computational complexity and memory become $O(LM^2N)$ and $O(M^2N)$ in Eq.(13), respectively. Therefore, the CQ-BEM is accelerated by the FMM in this research.

# 4. APPLICATION OF FMM TO CQ-BEM

The FMM, developed by Rokhlin and Greengard$^{(6)(7)}$, is a technique for reducing computational time and memory requirements for a large-scale problem. Since the FMM and its application to the BEM have been described in detail in other published papers, we will summarize only the essential formulas here.

## 4.1 MULTIPOLE EXPANSION AND LOCAL EXPANSION

To apply the FMM to the CQ-BEM, we need to separate the variables of the fundamental solutions (11) and (12). For example, the fundamental solution in Eq.(11) is transformed into the equation:

$$ \hat{U}_{ij}(\mathbf{x}, \mathbf{y}, s) = \frac{1}{\mu^*(s)s_T^2} [\Phi_{,i}^U + e_{3ij}\Psi_{,j}^U]. \quad (14) $$

Here $\Phi^U$ and $\Psi^U$ are the displacement potentials, which are defined by

$$ \Phi^U = \frac{1}{2\pi} K_0(s_L |\mathbf{x} - \mathbf{y}|)_k, \quad \Psi^U = e_{3kl} \frac{1}{2\pi} K_0(s_T |\mathbf{x} - \mathbf{y}|)_l. \quad (15) $$
Considering a point $\mathbf{y}_0$ near a source point $\mathbf{y}$ as shown in Fig.1, and using Graf's addition theorem for Eq.(15), we can obtain the multipole expansions as follows:

$$ \Phi^U(\mathbf{x}) = \frac{1}{2\pi} \sum_{n=-\infty}^{\infty} M_n^U K_n(s_L r) e^{in\theta}, \quad \Psi^U(\mathbf{x}) = \frac{1}{2\pi} \sum_{n=-\infty}^{\infty} N_n^U K_n(s_T r) e^{in\theta} \quad (16) $$

where $r$ and $\theta$ are the polar coordinate components of a field point $\mathbf{x}$ originated at $\mathbf{y}_0$ as shown in Fig.1. In Eq.(16), the coefficients $M_n^U$ and $N_n^U$ do not depend on the field point $\mathbf{x}$. These coefficients are called the multipole moment on $\mathbf{y}_0$, and are given by

$$ M_n^U = -\frac{\partial}{\partial y_k} [I_n(s_L \rho) e^{-in\phi}], \quad N_n^U = -e_{3kl} \frac{\partial}{\partial y_l} [I_n(s_T \rho) e^{-in\phi}]. \quad (17) $$

In Eq.(17), $I_n$ shows the modified Bessel function of the first kind. Similarly, multipole moments $M_n^T$ and $N_n^T$ for the double layer kernel defined in Eq.(12) are derived, as well as $M_n^U$ and $N_n^U$ for the fundamental solution defined in Eq.(11).

In addition, local expansions are also obtained corresponding to the multipole expansions defined in Eq.(16) as follows:

$$ \Phi(\mathbf{x}) = \sum_{n=-\infty}^{\infty} L_n I_n(s_L r) e^{-in\theta}, \quad \Psi(\mathbf{x}) = \sum_{n=-\infty}^{\infty} H_n I_n(s_T r) e^{-in\theta}. \quad (18) $$

In Eq.(18), $r$ and $\theta$ are the polar coordinate components of a field point $\mathbf{x}$ originated at the local expansion point $\mathbf{x}_0$. In addition, $L_n$ and $H_n$ are local expansion coefficients corresponding to multipole moments $M_n$ and $N_n$, respectively.

## 4.2 TRANSLATION FORMULAS FOR FMM

Once the multipole moments are obtained, we can quickly evaluate the matrix-vector products of the discretized boundary integral equation (13) using the fast multipole algorithm. The translation formulas for the fast multipole algorithm, M2M, M2L, and L2L, are also derived from Graf's addition theorem, and are defined as follows:

$$ M'_n = \sum_{m=-\infty}^{\infty} M_m I_{n-m}(s_L \rho) e^{-i(n-m)\phi}, \quad N'_n = \sum_{m=-\infty}^{\infty} N_m I_{n-m}(s_T \rho) e^{-i(n-m)\phi} \quad (\text{M2M}) \quad (19) $$

$$ L_n = \frac{1}{2\pi} \sum_{m=-\infty}^{\infty} (-1)^m M_m K_{n+m}(s_L \rho) e^{i(n+m)\phi}, \quad H_n = \frac{1}{2\pi} \sum_{m=-\infty}^{\infty} (-1)^m N_m K_{n+m}(s_T \rho) e^{i(n+m)\phi} \quad (\text{M2L}) \quad (20) $$

$$ L'_n = \sum_{m=-\infty}^{\infty} (-1)^{n-m} L_m I_{n-m}(s_L \rho) e^{i(n-m)\phi}, \quad H'_n = \sum_{m=-\infty}^{\infty} (-1)^{n-m} H_m I_{n-m}(s_T \rho) e^{i(n-m)\phi} \quad (\text{L2L}) \quad (21) $$

where $M'_n$, $N'_n$, and $L'_n$, $H'_n$ are new multipole and local expansion coefficients obtained by each translation formula. Moreover, In Eqs.(19)-(21), $\rho$ and $\phi$ are the polar coordinate components of an original multipole and local point based on a new multipole and local expansion point.

However, there is still one problem to resolve before we implement the CQ-BEM using the FMM. We note that the modified Bessel functions $I_n$ and $K_n$ used in the translation formulas tend exponentially to infinity for a large argument $z$ as follows:

$$ I_n(z) \sim e^{\Re[z]}, \quad K_n(z) \sim e^{-z}. \quad (22) $$
<figure><img src="image_3.png" /><figcaption>Figure 2. Viscoelastic wave scattering by a cavity.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 4. Vertical total displacement $u_2/u_0$ at point B on the cavity.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3. Vertical total displacement $u_2/u_0$ at point A on the cavity.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Vertical total displacement $u_2/u_0$ at point C on the cavity.</figcaption></figure>

This fact sometimes causes instability of the translation formulas (19)-(21) when cell size is large. To overcome this problem, we introduce a scaling of the multipole moments (17) using the scaled modified Bessel functions $\tilde{I}_n(z) = e^{-z}I_n(z)$ and $\tilde{K}_n(z) = e^z K_n(z)$.

# 5. NUMERICAL EXAMPLES

In the following practical applications, the three-element standard linear model(8) is considered as a typical elastic model.

## 5.1 VISCOELASTIC WAVE SCATTERING BY A CAVITY

We first consider the problem of an incident plane SV-wave scattered by a cavity with radius $a$ as shown in Fig.2. In this analysis, we use a Ricker wavelet modified for a viscoelastic wave, as the incident wave. The modified Ricker wavelet $u_i^{\text{in}}(x, t)$ is defined by a Fourier transform of the original Ricker wavelet as

$$u_i^{\text{in}}(x, t) = u_0 \delta_{i2} \mathcal{F}^{-1} \left[ 2t_0 \sqrt{\pi} \left( \frac{wt_0}{2} \right)^2 e^{-\left( \frac{wt_0}{2} \right)^2} e^{i(k_T x_1 + wt_s)} \right] \quad (23)$$

where $t_0$ corresponds to the peak frequency $2\pi/t_0$, $t_s$ corresponds to the minimum peak in the time domain and $u_0$ is the reference displacement. The difference from the elastic wave is that the wavenumber
<figure><img src="image_7.png" /><figcaption>Figure 6. Multiple scattering model.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Time variations of $u_1/u_0$ of the incident viscoelastic wave calculated by Eq.(24) at $x_1/a = 0.0, 5.0, 10.0$ and $15.0$.</figcaption></figure>

$k_T$ is complex and depends on the angular frequency $\omega$.

Figs.3-5 show the time histories of the displacements $u_2/u_0$ at points A, B, and C on the cavity as shown in Fig.2. Solutions using two different time increments, $c_{L0}\Delta t/a = 0.125$ and $c_{L0}\Delta t/a = 0.0625$, an elastic solution, and a numerical solution obtained by Fukui et. al. (9) are shown in Figs.3-5 for comparison. We solved this problem solely through the CQ-BEM without using the FMM. The boundary of the cavity is discretized into 72 boundary elements using piecewise constant approximation, and we assumed a traction-free condition: $t_i = 0$. The CQM parameters, $N$ and $L$, are given by $N = L = 128$. In addition, $\mathcal{R}$ is assumed to be $\mathcal{R} = 0.91398170$ ($\epsilon = 10^{-10}$). The material parameters of the three-element standard linear model are $\mu_R/\mu_0 = 0.5$, $K/\mu_0 = 1.0$, $\tau_\sigma = T_0$, and $\tau_\epsilon = 0.5T_0$, where $\mu_R$ is the relaxed shear modulus, $\mu_0$ is the initial shear modulus, $\tau_\sigma$ is the stress relaxation time, and $\tau_\epsilon$ is the strain relaxation time. Moreover, the time $T_0$ is defined by $T_0 = 2a/c_{L0}$, which is the time that an elastic wave with initial wave velocity $c_{L0} = \sqrt{(K + (4/3)\mu_0)/\rho}$ passes through a cavity with radius $a$ as shown in Fig.2. We set time $t/T_0 = 0.0$ at the moment the incident wave hits point A on the boundary of the cavity in Fig.2.

We can see that the results obtained by the CQ-BEM are in good agreement with the analytical numerical results calculated by Fukui et al.

## 5.2 VISCOELASTIC WAVE SCATTERING BY MANY CAVITIES

Next, viscoelastic wave scattering by many cavities, as shown in Fig.6, is solved using the CQ-BEM accelerated by the FMM. In this analysis, we use a desktop computer with an Intel(R) Xeon(R) 1.60GHz CPU 5310 and 2GB of memory. The components of an incident viscoelastic plane wave with initial wave velocity $c_{L0}$ are defined as follows:

$$u_i^{\text{in}}(\mathbf{x}, t) = u_0 \delta_{i1} \mathcal{F}^{-1} \left[ \frac{\omega_0^2 (1 - e^{i\omega T_0})}{2i\omega(\omega^2 - \omega_0^2)} e^{ik_L x_1} \right] \quad (24)$$

where $\omega_0 = 2\pi/T_0$. Fig.7 shows the time variations of $u_1/u_0$ of the incident viscoelastic waves calculated by Eq.(24) at $x_1/a = 0, 5, 10$, and $15$. The parameters of the three-element standard model are $\mu_R/\mu_0 = 0.85$, $K/\mu_0 = 5/3$, $\tau_\sigma = 0.5T_0$, and $\tau_\epsilon = 17T_0/40$. For comparison, incident elastic waves with wave velocity $c_{L0}$ corresponding to the incident viscoelastic waves are plotted by solid lines in Fig.7. We can see that the wave-forms of the incident viscoelastic wave decrease while the wave is
<figure><img src="image_9.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c)</figcaption></figure>

<caption>Figure 8. Time histories of the displacements $u_1/u_0$ Figure 9. Time histories of the displacements $u_1/u_0$</caption>

(a) $t/T_0 = 1.25$, (b) $t/T_0 = 6.25$, (c) $t/T_0 = 11.25$.

<figure><img src="image_12.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c)</figcaption></figure>

around the cavities in the elastic case. (a)

$t/T_0 = 1.25$, (b) $t/T_0 = 6.25$, (c) $t/T_0 = 11.25$.

propagating.

Fig.8-(a), (b), and (c) show the time histories of displacement $u_1/u_0$ around cavities at $t/T_0 = 1.25$,
$t/T_0 = 6.25$, and $t/T_0 = 11.25$, respectively. The results for the elastic case are also shown for comparison in Fig.9-(a), (b), and (c). We set the time $t/T_0 = 0.0$ at the moment the incident wave hits the cavities in the far left column in Fig.6. The parameters for the CQM are taken as $N = L = 256$, $\rho = 0.95609320$ ($\epsilon = 10^{-10}$). The time increment is $c_{L0} \Delta t/a = 0.125$, the number of DOF in each time step is 8192, and the truncation number for the summation of multipole expansion and local expansion such as Eq.(16) and Eq.(18) is chosen as $p = \pm 10$. This problem cannot be solved by the use of the CQ-BEM alone because of memory restrictions. Therefore, the FMM is applied to accelerate the calculation of the matrix-vector products of the discretized boundary integral equation and to conserve memory. In addition, OpenMP with eight threads is used to parallelize this analysis. We can see that scattered waves are generated by the interaction of the incident wave with each cavity. The effect of attenuation from the viscoelastic model can be seen by comparing Fig.8 and Fig.9. Thus, the CQ-BEM using the FMM is very effective for large-scale problems.

## 6. CONCLUSION

In this paper, the CQ-BEM using the FMM was developed for 2-D viscoelastic wave propagation. In the proposed method, the convolution integrals in a time-domain boundary integral equation are discretized using the CQM, and the Laplace-domain fundamental solutions are used to calculate influence functions. Therefore, the proposed method has the potential to tackle problems where no time-domain fundamental solutions are available, e.g., poroelastodynamics or 3-D viscoelastodynamics. Moreover, the FMM was applied to accelerate the calculations of matrix-vector products for the retarded potential and to reduce the memory required in the CQ-BEM. As numerical examples, scattering problems of incident waves by cavities were demonstrated to validate the proposed method. In the near future, we will extend the CQ-BEM using the FMM to 3-D viscoelastic wave propagation.

## REFERENCES

1) Mansur, W. J. and Brebbia, C. A., "Transient elastodynamics using a time-stepping technique", In: Boundary Elements, C. A. Brebbia, T. Futagami and M. Tanaka (Eds), (1983), pp. 677-698.

2) Hirose, S., "Boundary Integral equation method for transient analysis of 3-D cavities and inclusions", *Engineering analysis with Boundary Elements*, 8-3, (1991), pp. 146-153.

3) Schanz, M., "Wave propagation in viscoelastic and poroelastic continua", A Boundary Element Approach, Springer, (2001).

4) Abreu, A. I. and Carrer, J. A. M. and Mansur, W. J., "Scalar wave propagation in 2D: a BEM formulation based on the operational quadrature method", *Engineering analysis with Boundary Elements*, 27, (2003), pp. 101-105.

5) Lubich, C., "Convolution quadrature and discretized operational calculus I", *Numer. Math.*, 52, (1988), pp. 129-145.

6) Rokhlin, V., "Rapid solution of integral equations of classical potential theory", *J. Comput. Phys.*, Vol.60, (1985), pp. 187-207.

7) Greengard, L. and Rokhlin, V., "A fast algorithm for particle simulations", *J. Comput. Phys.*, Vol.73, (1987), pp. 325-348.

8) Fung, Y. C., "Foundation of Solid Mechanics", Prentice-Hall, (1965).

9) Kobayashi, S. et. al., "Wave analysis and boundary element methods", Kyoto University Press (in Japanese).