# An efficient stabilized boundary element formulation for 2D time-domain acoustics and elastodynamics

Delfim Soares, Webe João Mansur

► To cite this version:

Delfim Soares, Webe João Mansur. An efficient stabilized boundary element formulation for 2D time-domain acoustics and elastodynamics. Computational Mechanics, 2007, 40 (2), 10.1007/s00466-006-0104-3. hal-01365306

HAL Id: hal-01365306

https://hal.science/hal-01365306

Submitted on 13 Sep 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
# An efficient stabilized boundary element formulation for 2D time-domain acoustics and elastodynamics

D. Soares Jr. · W. J. Mansur

**Abstract** The present paper describes a procedure that improves efficiency, stability and reduces artificial energy dissipation of the standard time-domain direct boundary element method (BEM) for acoustics and elastodynamics. Basically, the developed procedure modifies the boundary element convolution-related vector, being very easy to implement into existing codes. A stabilization parameter is introduced into the recent-in-time convolution operations and the operations related to the distant-in-time convolution contributions are approximated by matrix interpolations. As it is shown in the numerical examples presented at the end of the text, the proposed formulation substantially reduces the BEM computational cost, as well as its numerical instabilities.

**Keywords** Time-domain BEM · Acoustics · Elastodynamics · Time convolution · Stability

## 1 Introduction

The very first works concerning time-domain boundary element methods (TDBEM) applied to acoustic and el-

astodynamic problems reported numerical instabilities [6,8,10–12], which can become quite severe for late time.

There are nowadays many publications that report robust techniques to improve TDBEM stability. In the approach developed by Frangi [9], space and time shape functions are not independent as in standard TDBEM approaches: a modified space–time approximation, which represents more accurately the causality of the phenomenon, is employed. Yu et al. [22,23] achieved more stable TDBEM algorithms by introducing, respectively, extra time and space (Galerkin approach) weighting integrals. Three other quite simple and effective schemes presently available are the linear $\theta$ method [14,20,21,24], the $\varepsilon$ scheme and the half-step scheme [1,15]. Many other stabilization schemes can be found in the literature; some being more complex than those mentioned above (see [4,5]) or else, dedicated to minimize one specific source of instability [3,13].

The present paper describes a new procedure that improves stability and artificial energy dissipation of TDBEM for acoustic and elastodynamic modelling, besides reducing drastically CPU time and storage area. The developed procedure ($\alpha$-$\delta$ method) utilizes: (1) an $\alpha$ stabilization parameter to modify the recent-in-time convolution operations, by replacing current values by weighted averages; (2) matrix interpolations to approximate the distant-in-time convolution contributions (the $\delta$ parameter controls the size of the recent/distant-in-time regions). The algorithm hereby developed, as one will notice, is quite efficient and stable.

The $\alpha$-$\delta$ TDBEM presented here associates stability and efficiency, and it is expected to be also more robust than previously reported TDBEM transmitting boundaries (as for instance in a BEM–FEM coupling context [16,17,24]).

<footnote>D. Soares Jr. · W. J. Mansur (✉)
Civil Engineering Department, COPPE/UFRJ,
CXP 68506, CEP 21945-970 Rio de Janeiro, RJ, Brazil
e-mail: webe@coc.ufrj.br</footnote>

<footnote>D. Soares Jr.
Structural Engineering Department,
UFJF, Cidade Universitária,
CEP 36036-330 Juiz de Fora, MG, Brazil</footnote>

<footnote>D. Soares Jr.
e-mail: delfim@coc.ufrj.br</footnote>
## 2 Standard boundary element formulation

The time-domain integral representation for acoustic and elastodynamic wave propagation, in the absence of initial conditions and internal sources/forces, can be written as described by Eqs. (1a) and (1b), respectively

$$ 
\begin{aligned}
c(\xi)p(\xi, t) &= \int_0^{t_+} \int_0^{t_+} p * (X, t; \xi, \tau) q(X, \tau) d\Gamma(X) d\tau \\ 
&\quad -FP \int_0^{t_+} \int_0^{t_+} q * (X, t; \xi, \tau) p(X, \tau) d\Gamma(X) d\tau
\end{aligned} 
\quad (1a) $$

$$ 
\begin{aligned}
c_{ik}(\xi) u_k(\xi, t) &= \int_0^{t_+} \int_0^{t_+} u_{ik} * (X, t; \xi, \tau) \tau_k(X, \tau) d\Gamma(X) d\tau \\ 
&\quad -FP \int_0^{t_+} \int_0^{t_+} \tau_{ik} * (X, t; \xi, \tau) u_k(X, \tau) d\Gamma(X) d\tau
\end{aligned} 
\quad (1b) $$

where $c(\xi)$ and $c_{ik}(\xi)$ are geometrical parameters; $p * (X, t; \xi, \tau)$ represents the potential acoustic fundamental solution at field point $X$, corresponding to a unit impulsive source applied at point $\xi$ at time $\tau$; $u_{ik} * (X, t; \xi, \tau)$ represents the displacement elastodynamic fundamental solution component in the $k$ direction at field point $X$, corresponding to a unit impulsive force applied on the $i$ direction at point $\xi$ at time $\tau$; $q * (X, t; \xi, \tau)$ and $\tau_{ik} * (X, t; \xi, \tau)$ are the flux and traction related fundamental solutions, respectively. FP stands for the finite part of the integral.

The boundary integral equations (1) can be treated numerically [2,10], resulting in a system of equations that can be expressed in matrix form. Adopting $H^n$ and $G^n$ as being the matrices computed at the current time-step $n$, related to the flux/traction fundamental solution and the potential/displacement fundamental solution, respectively, the system of equations (2) is obtained

$$ \mathbf{C}\mathbf{V}_{\mathrm{H}}^{n}=\sum_{l=1}^{n}\left(\mathbf{G}^{n-l+1} \mathbf{V}_{\mathrm{G}}^{l}-\mathbf{H}^{n-l+1} \mathbf{V}_{\mathrm{H}}^{l}\right) \quad (2) $$

where, in the case of acoustics [Eq. (1a)], $\mathbf{V}_H$ and $\mathbf{V}_G$ entries are nodal values of the boundary potentials ($p(X, \tau)$) and fluxes ($q(X, \tau)$), respectively. On the other hand, if an elastodynamic problem is being considered [Eq. (1b)], $\mathbf{V}_H$ and $\mathbf{V}_G$ are the vectors in which the nodal boundary displacement ($u_k(X, \tau)$) and traction ($\tau_k(X, \tau)$) components are assembled, respectively.

After introducing the boundary conditions in Eqs. (2), the following final expressions are obtained

$$ \mathbf{A} \mathbf{X}^n = \mathbf{B} \mathbf{Y}^n + \mathbf{R}^n \quad (3a) $$

$$ \mathbf{R}^n = \sum_{l=1}^{n-1} \left( \mathbf{G}^{n-l+1} \mathbf{V}_G^l - \mathbf{H}^{n-l+1} \mathbf{V}_H^l \right) \quad (3b) $$

where, as usual in time-domain BEM, the entries of $\mathbf{X}^n$ are unknown variables at boundary nodes at discrete time $t_n$, while the entries of vector $\mathbf{Y}^n$ are the according known nodal values. Matrices $\mathbf{A}$ and $\mathbf{B}$ are obtained from the combination of matrices $\mathbf{C}+\mathbf{H}^1$ and $\mathbf{G}^1$, taking into account, as well, the prescribed boundary conditions. $\mathbf{R}^n$ is the vector related to the convolution process of the BEM; it represents the complete history up to $t_{n-1}$. The time stepping procedure indicated by Eq. (3) requires convolutions ($l = 1, 2, \dots, n-1$) to be carried out for $n = 1, 2, \dots, N$, where the final time is $t_N = N\Delta t$, $\Delta t$ being the time step.

Further details on the implementation of standard time-domain BEM algorithms can be found, for instance, in Mansur [10], Carrer and Mansur [2], Dominguez [7] etc.

## 3 Efficient stabilized boundary element formulation

In order to achieve an efficient stabilized boundary element formulation, the present work modifies the BEM convolution vector [Eq. (3b)], introducing a stabilization parameter into the recent-in-time convolution operations and a time-truncation procedure to compute the distant-in-time convolution contributions. No modification is introduced in the evaluation of the standard BEM matrices $\mathbf{H}^n$ and $\mathbf{G}^n$, only their manipulation along the convolution process is modified. Thus, the proposed formulation is easy to implement into standard BEM existing codes.

The algorithm that describes the developed boundary element formulation is represented by Eqs. (4) below (which are analogous to the standard BEM Eqs. (3)). Details about the elaboration of the terms in Eqs. (4) are presented next, in the sub-sections that follow.

$$ \bar{\mathbf{A}} \mathbf{X}^n = \bar{\mathbf{B}} \mathbf{Y}^n + \bar{\mathbf{R}}^n + \hat{\mathbf{R}}^n \quad (4a) $$

$$ \bar{\mathbf{R}}^n = \sum_{l=n-L+1}^{n-1} \left\{ \mathbf{G}^{n-l+1} \sum_{k=-1}^{1-\delta_{l,n-1}} \left[ J(\alpha, k) \mathbf{V}_G^{l+k} \right] - \mathbf{H}^{n-l+1} \sum_{k=-1}^{1-\delta_{l,n-1}} \left[ J(\alpha, k) \mathbf{V}_H^{l+k} \right] \right\} \quad (4b) $$
$$ \hat{\mathbf{R}}^n = \sum_{k=1}^{m} \left\{ \mathbf{G}_k \sum_{l=1}^{n-L} \left[ I(n-l+1, k) \mathbf{V}_G^l \right] - \mathbf{H}_k \sum_{l=1}^{n-L} \left[ I(n-l+1, k) \mathbf{V}_H^l \right] \right\} \quad (4c) $$

In Eq. (4), the entries of $\mathbf{X}^n$ and $\mathbf{Y}^n$ are, once more, unknown and known nodal values, respectively, at boundary nodes at discrete time $t_n$. Matrices $\tilde{\mathbf{A}}$ and $\tilde{\mathbf{B}}$ are obtained from the combination of matrices $\mathbf{C} + \mathbf{H}^1 + J(\alpha, 1)\mathbf{H}^2$ and $\mathbf{G}^1 + J(\alpha, 1)\mathbf{G}^2$, taking into account the prescribed boundary conditions of the problem ($\alpha$ is a stabilization parameter). $\hat{\mathbf{R}}^n$ is the vector related to the recent-in-time convolution process (it represents the time history from $t_{n-L+1}$ up to $t_{n-1}$, where $L$ is a parameter which defines the truncation point of the BEM convolution process) and $\hat{\mathbf{R}}^n$ is a truncated BEM convolution vector. $J(\cdot)$ and $I(\cdot)$ are scalar functions; they define stabilization coefficients and interpolation parameters, respectively. $\delta_{i,j} = 0$ if $i \neq j$ and $\delta_{i,j} = 1$ if $i = j$. Further details about matrices $\tilde{\mathbf{A}}$ and $\tilde{\mathbf{B}}$ and vector $\hat{\mathbf{R}}^n$ are given in Sub-sect. 3.1; details about vector $\hat{\mathbf{R}}^n$ are given in Sub-sect. 3.2.

## 3.1 Recent-in-time convolution contributions

For the recent-in-time [i.e., $n-L+1 \le l \le n$ in Eq. (3b)] convolution operations, an $\alpha$ parameter, which linearly combines three consecutive time-step results, is introduced for the stabilization of the standard boundary element formulation. The combination of consecutive results has been used in several applications (adoption of a relaxation parameter in iterative processes, for instance [16,17]), showing good performance in what concerns the stabilization of an oscillatory tendency. The basic idea here is to smooth the convolution operations, trying to avoid any possible numerical instability. Several algorithms may be developed to smooth these operations based on the combination of consecutive time-step results; the present work adopts a simple linear combination of three time-consecutive terms, as described by

$$ \bar{\mathbf{V}}^l = \left( \frac{1}{2}\alpha \right) \left( \mathbf{V}^{l-1} + \mathbf{V}^{l+1} \right) + (1 - \alpha) \mathbf{V}^l \quad (5) $$

By replacing $\mathbf{V}^l$ by $\bar{\mathbf{V}}^l$ in Eq. (3b), one may obtain Eq. (4b), as well as matrices $\tilde{\mathbf{A}}$ and $\tilde{\mathbf{B}}$. Taking into account Eq. (5), the coefficient function $J(\cdot)$, presented in Eq. (4b), is given by

$$ J(\alpha, k) = \frac{1}{2}\alpha k^2 + (1 - k^2)(1 - \alpha) \quad (6) $$

As it will be illustrated in the Sect. 4, this algorithm is very effective, reducing considerably typical TDBEM numerical oscillations, giving as a result a stable algorithm. Not only this formulation completely eliminates or substantially reduces usual instabilities by adopting $0.5 \le \alpha \le 1.0$, but also it is extremely simple and easy to implement into standard BEM existing codes.

## 3.2 Distant-in-time convolution contributions

This section discussion is concerned with a time-truncation approach where operations related to the distant-in-time [i.e., $l \le n-L$ in Eq. (3b)] contribution of the convolution process are effectuated approximately. This approximation consists of evaluating the matrices $\mathbf{H}^{n-l+1}$ and $\mathbf{G}^{n-l+1}$ by interpolation, based on a few m-calculated matrices $\mathbf{H}_k$ and $\mathbf{G}_k$ ($k = 1, 2, \dots, m$), computed at appropriate discrete times $T_k$ ($T_L \le T_k \le t_N$, where $T_L$ is the time limit after which the approximations take place). The integer $m$ is an input parameter to the approximated convolution analysis and it indicates the number of key time-points ($T_k$) to be used in the interpolation procedure (within the interlude of these key time-points, the interpolations will occur).

In order to obtain matrices $\mathbf{H}^{n-l+1}$ and $\mathbf{G}^{n-l+1}$ as function of $m$ previously selected matrices $\mathbf{H}_k$ and $\mathbf{G}_k$, an interpolation parameter $I(\cdot)$ is introduced into Eq. (3b), originating Eq. (4c). There are several interpolation procedures that might be chosen to specify $I(\cdot)$; for a detailed discussion, the reader is referred to Soares and Mansur [18]. The present work adopts a multi-linear interpolation algorithm; i.e., linear interpolation between the $m$ interpolation time-points is adopted. Taking into account the multi-linear interpolation procedure, the $I(\cdot)$ parameter is given by

$$ I(j, k) = \frac{j\Delta t - T_{k-1}}{T_k - T_{k-1}} [H(j\Delta t - T_{k-1})H(T_k - j\Delta t)](1 - \delta_{k,1}) + \frac{T_{k+1} - j\Delta t}{T_{k+1} - T_k}[H(j\Delta t - T_k)H(T_{k+1} - j\Delta t)](1 - \delta_{k,m}) \quad (7) $$

where $H(\cdot)$ is the Heaviside function. In the present work, the discrete times $T_k$, which are the key-times to the interpolation approximation, are adopted as

$$ T_k = T_L + (T_N - T_L) \left( \frac{k-1}{m-1} \right)^{m'} \quad (8) $$

and the time limit $T_L$ is adopted as

$$ T_L = \Delta t N^{(1-\delta m)} + \frac{\bar{r}}{c} \quad (9) $$
where $\delta$ is a precision control parameter; $\bar{r}$ is the average distance between boundary nodes and $c$ is the (primary) wave velocity. The calculus of the time limit $T_L$ [Eq. (9)] is based on the behaviour of the function $f(t) = 1/t$, which is a simple function, similar in most aspects to the time related kernels of the models being considered.

As it can be seen in Eq. (8), the values of the discrete times $T_k$ are not necessarily equally spaced once the input parameter $\eta$ is introduced. This scheme is adopted since better approximations are obtained if more interpolation points are chosen in the neighbourhood of $T_L$, region in which kernels gradient and amplitude are more significant. The discrete times $T_k$ become equally spaced if $\eta = 0$.

It is important to observe that the kernels that are here being interpolated decrease monotonically with time increase; thus, linearly interpolated values are always greater than real ones. This property is beneficial: as it is well known, standard TDBEM algorithms usually introduce some spurious numerical damping; as a consequence, the multi-linear interpolation approximation acts as a subtle dispersion-correction algorithm (besides drastically reducing the BEM computational cost).

## 4 Numerical applications

Time-step length plays an important role in the time-domain BEM analyses [10]. A measure of the time-step length is computed according to the following expression

$$ \beta = \frac{c\Delta t}{\ell} \qquad (10) $$

where $c$ is the (primary) wave velocity and $\ell$ is the boundary element length.

In the present section, two numerical applications are considered. The first application (finite-domain problem) is concerned with an acoustic rod; the discussion presents numerical results obtained taking into account different discretization ($\beta$) and stabilization ($\alpha$) parameters. In the second application (infinite-domain problem), an elastodynamic problem, composed by two close circular cavities, is focused. In both applications, typical numerical instabilities occur for standard TDBEM analyses, according to the parameter $\beta$ adopted. In fact, taking into account numerical aspects, these problems are quite sensitive, due to the multiple wave reflections that occur.

As one can observe in the following applications, the proposed methodology deals quite well with the

usual time-domain BEM limitations, namely: numerical instabilities and excessive oscillation, energy dispersion (numerical damping) and high computational cost.

### 4.1 Acoustic rod

This example [10] consists of an acoustic rectangular domain, submitted to a time Heaviside type flux applied at one of its end and prescribed null potentials applied at the opposite one. Figure 1a shows the boundary conditions ($f(t) = H(t)$) and geometry ($a = 4.0 \text{ m}$ and $b = 2.0 \text{ m}$) of the body. Figure 1b shows the boundary discretization: the symmetry of the problem was considered and only 24 linear boundary elements of equal size ($\ell = 0.25 \text{ m}$) were employed (an interesting feature of the boundary element formulation is that symmetric bodies under symmetric loads can be analysed without discretization of the symmetry axes; this can be accomplished by an automatic condensation process, which integrates over reflected elements and performs the assemblage of the final coefficients in a reduced matrix [19]). The solution of this 2D plane wave problem is in fact independent of the coordinate in the vertical direction, i.e., the plane wave front propagates in the longitudinal direction. Thus, its solution is the same as that of longitudinal waves propagating along a rod.

The wave propagation velocity along the rod is $c = 2,500 \text{ m/s}$. Different time discretizations have been adopted for the analyses. The time history results concerning potentials at points A and B and fluxes at point C are depicted in Figs. 2,3,4 and 5, for different discretization (namely: $\beta = 0.2, \beta = 0.3, \beta = 0.4$ and $\beta = 0.5$) and stabilization (namely: $\alpha = 0.0, \alpha = 0.3, \alpha = 0.5$ and $\alpha = 0.7$) parameters. As one may observe, the proposed formulation is quite stable (one should note that the standard BEM is described by $\alpha = 0.0$).

In the above-described analyses, $\delta = 0.00$ was adopted, i.e., no truncation on the convolution process was considered. Taking into account the approximations presented in Sub-sect. 3.2 for the distant-in-time convolution contributions ($\delta = 0.04, m = 10$ and $\eta = 0.5$)

<figure><img src="image_2.png" /><figcaption>Fig. 1 Acoustic rod: a geometry, boundary conditions and selected points; b boundary elements discretization</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Potential and flux time history results for $\beta = 0.2$: **a** $\alpha = 0.0$ (standard BEM); **b** $\alpha = 0.3$; **c** $\alpha = 0.5$; **d** $\alpha = 0.7$</figcaption></figure>

were adopted), some time history results are presented in Fig. 6, for $\alpha = 0.5$ and for four different $\beta$-values, namely: $\beta = 0.25$, $\beta = 0.5$, $\beta = 0.75$ and $\beta = 1.00$. As one may notice, the results are once more stable and no perceptible difference, concerning accuracy, may be observed among the truncated ($\delta = 0.04$) and the complete ($\delta = 0.00$) convolution analyses. However, concerning efficiency aspects, the difference among the analyses is quite outstanding. The relations among the truncated and the complete time convolution analyses are described in Table 1, considering two computational-cost

aspects: total CPU time and storage area. As one can observe, the developed formulation may easily reduce more than 80% the time-domain BEM computational cost (one should observe that in Table 1 the relation “modified approach cost/standard approach cost” is presented, in percentage).

In Fig. 7, another interesting feature of the presented formulation is depicted. As it has been mentioned, the multi-linear interpolation approximation acts as a subtle dispersion-correction algorithm, once standard time-domain BEM analyses are usually associated with some
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Potential and flux time history results for $\beta = 0.3$; **a** $\alpha = 0.0$ (standard BEM); **b** $\alpha = 0.3$; **c** $\alpha = 0.5$; **d** $\alpha = 0.7$</figcaption></figure>

numerical damping. By increasing the error parameter ($\delta = 0.06$), which controls the time limit ($T_L$) of the distant-in-time convolution contributions, the kernels are over-estimated and perceptible differences among the truncated and complete convolution results may be observed. However, if not excessive, this error is beneficial: it corrects the BEM numerical damping tendency. A clear exemplification of the above aspect is depicted in Fig. 7, where the results of the truncated analysis may be considered more accurate than the ones related to the complete convolution analysis.

## 4.2 Elastodynamic cavities

This example consists of two circular cavities located close to each other in an unbounded elastodynamic medium. One cavity is submitted to a time Heaviside type pressure-load, whereas the other remains unloaded. Figure 8a shows the boundary conditions ($f(t) = H(t)$) and geometric configuration ($R = 1.0$ m) of the problem. Figure 8b shows the cavities boundary discretization: 16 linear boundary elements of equal size ($\ell = 0.3902$ m) were adopted to model each cavity.
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Potential and flux time history results for $\beta = 0.4$; **a** $\alpha = 0.0$ (standard BEM); **b** $\alpha = 0.3$; **c** $\alpha = 0.5$; **d** $\alpha = 0.7$</figcaption></figure>

The physical properties of the medium are: $E = 70.0 \text{ N/m}^2$ (Young Modulus); $\nu = 0.2$ (Poisson rate); $\rho = 2.7 \text{ kg/m}^3$ (mass density). Two different time-step lengths were adopted for the analyses, namely: $\Delta t = 0.01 \text{ s}$ ($\beta = 0.1376$) and $\Delta t = 0.01455 \text{ s}$ ($\beta = 0.2001$).

The time history results concerning normal displacements at points A and B are presented in Fig. 9, for different discretization ($\beta = 0.1376$ and $\beta = 0.2001$) and stabilization ($\alpha = 0.0, \alpha = 0.5, \alpha = 0.75$ and $\alpha = 1.00$) parameters. As one may once again observe, the pro-

posed formulation considerably reduces or completely removes the usual BEM numerical instabilities.

Taking into account interpolation approximations for the distant-in-time convolution contributions ($\delta = 0.03, m = 10$ and $\eta = 0.5$ were adopted), some time history results are presented in Fig. 10, for $\alpha = 0.75$ and for $\beta = 0.1376$ and $\beta = 0.2001$. The computational cost relations between the truncated ($\delta = 0.03$) and the complete ($\delta = 0.00$) time convolution analyses are described in Table 1. As it may be observed, the proposed formu-
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Potential and flux time history results for $\beta = 0.5$: **a** $\alpha = 0.0$ (standard BEM); **b** $\alpha = 0.3$; **c** $\alpha = 0.5$; **d** $\alpha = 0.7$</figcaption></figure>

lation once again substantially reduced the time-domain BEM computational cost, without introducing any perceptible accuracy difference.

## 5 Conclusions

Presently it is possible to find in the specialized literature many different schemes, which can either eliminate, or else reduce to quite acceptable levels the instability observed in many TDBEM applications concerning acoustics and elastodynamics.

Here, a new scheme has been presented, the $\alpha-\delta$ scheme, which completely eliminates or substantially reduces instability, decreases artificial (numerical) energy dissipation and improves computational efficiency. The procedure can be employed either for 2D or 3D applications; in the former, the described compression scheme can be used to drastically reduce CPU time and storage area.

Two examples were presented where a parametric study concerning the $\alpha, \delta$ and $\beta$ analysis parameters was carried out, the conclusions being that quite a wide
<figure><img src="image_7.png" /><figcaption>**Table 1** Computational cost relations between the truncated and the complete time convolution analyses</figcaption></figure>

<table><thead><tr><th>Problem</th><th>β</th><th>CPU time (%)</th><th>Storage area (%)</th></tr></thead><tbody><tr><td rowspan="4">Acoustics δ = 0.04</td><td>0.25</td><td>15</td><td>8</td></tr><tr><td>0.50</td><td>18</td><td>11</td></tr><tr><td>0.75</td><td>23</td><td>13</td></tr><tr><td>1.00</td><td>24</td><td>14</td></tr><tr><td rowspan="2">Elastodynamics δ = 0.03</td><td>0.1376</td><td>20</td><td>13</td></tr><tr><td>0.2001</td><td>26</td><td>15</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 8** Elastodynamic cavities: a geometry, boundary conditions and selected points; b boundary elements discretization</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Potential time history results ($\alpha = 0.5$ and $\beta = 2.0$) at point A for $\delta = 0.00$ (complete convolution), $\delta = 0.04$ and $\delta = 0.06$ (truncated convolution)</figcaption></figure>

interval of variation of these parameters lead to stable numerical results; many thousands of time steps could be used even in the most severe case (first example), which is prone to numerical instabilities. In fact, the proposed methodology efficiently managed to overcome the main drawbacks of TDBEM formulations: the introduction of the $\alpha$ parameter is extremely efficient dealing with the TDBEM instability problem and its awkward dependency on the space-time discretization parameter $\beta$; moreover, the introduction of the $\delta$ parameter overcomes the high computational cost of the time convolution process, introducing controlled errors (it can be as small as the user wants), which are beneficial to correct the numerical damping tendency of the TDBEM formulation.

The present work smoothed the recent-time convolution operations by a simple weighted linear combination
<figure><img src="image_10.png" /><figcaption>**Fig. 10** Displacement time history results ($\alpha = 0.75$) at point A for $\delta = 0.00$ (complete convolution) and $\delta = 0.03$ (truncated convolution): **a** $\beta = 0.1376$; **b** $\beta = 0.2001$</figcaption></figure>

<figure><img src="image_11.png" /></figure>
of three time-consecutive elements; however, many
other equally simple derived alternatives are possible.
The authors suggest as future work the development of
such simple alternative schemes, and a comparison with
published algorithms (see Sect. 1) in order to identify
the optimum (one or more). Another point to be inves-
tigated concerns coupling of more than one algorithm,
specially time weighting and Galerking with the present
procedure, linear θ, ε scheme, half-step etc.

References

1. Birgisson B, Siebrits E, Peirce AP (1999) Elastodynamic direct boundary element methods with enhanced numerical stability properties. Int J Numer Meth Eng 46:871–888

2. Carrer JAM, Mansur WJ (1993) Two-dimensional transient BEM analysis for the scalar wave equation: kernels. Eng Anal Bound Elem 12:283–288

3. Carrer JAM, Mansur WJ (2000) Time discontinuos linear trac-
tion approximation in time-domain BEM: 2D elastodynamics.
Int J Numer Meth Eng 49:833–848

4. Carrer JAM, Mansur WJ (2002) Time-dependent fundamen-
tal solution generated by a not impulsive source in the bound-
ary element method analysis of the 2D scalar wave equation.
Commun Numer Meth Eng 18:277–285

5. Coda HB (1993) Tri-dimensional transient analysis of struc-
tures by the combination of boundary and finite elements (in
portuguese) DSc Thesis, USP/São Carlos

6. Cole DM, Kosloff DD, Minster JB (1978) Numerical bound-
ary integral-equation method for elastodynamics. 1. Bull Seis-
mol Soc Am 68:1331–1357

7. Dominguez J (1993) Boundary elements in dynamics. Com-
putational Mechanics Publications, Southampton, Boston

8. Friedman MB, Shaw R (1962) Diffraction of pulse by cylindri-
cal obstacles of arbitrary cross section. Trans ASME J Appl
Mech 25:40–46

9. Frangi A (2000) Causal shape functions in the time domain boundary element method. Comput Mech 25:533–541

10. Mansur WJ (1983) A time-stepping technique to solve wave propagation problems using the boundary element method. PhD Thesis, University of Southampton

11. Mansur WJ, Brebbia CA (1982) Formulation of the bound-
ary element method for transient problems governed by the
scalar wave-equation. Appl Math Model 6:307–311

12. Mansur WJ, Brebbia CA (1982) Numerical implementation of the boundary element method for 2 dimensional transient scalar wave-propagation problems. Appl Math Model 6:299–306

13. Mansur WJ, Carrer JAM, Siqueira EFN (1998) Time discon-
tinuous linear traction approximation in time-domain BEM
scalar wave propagation analysis. Int J Numer Meth Eng
42:667–683

14. Mansur WJ, Yu G, Carrer JAM, Lie ST (2000) The θ scheme for time-domain BEM/FEM coupling applied to the 2-D scalar wave equation. Commun Numer Meth Eng 16:439–448

15. Siebrits E, Peirce AP (1997) Implementation and applica-
tion of elastodynamic boundary element discretization with
improved stability properties. Eng Comput 14:669–695

16. Soares D Jr, von Estorff O, Mansur WJ (2004) Iterative cou-
pling of BEM and FEM for nonlinear dynamic analyses. Com-
put Mech 34:67–73

17. Soares D Jr, von Estorff O, Mansur WJ (2005) Efficient nonlinear solid-fluid interaction analysis by an iterative BEM/FEM coupling. Int J Numer Meth Eng 64:1416–1431

18. Soares D Jr, Mansur WJ (2004) Compression of time gener-
ated matrices in two-dimensional time-domain elastodynamic
BEM analysis. Int J Numer Meth Eng 61:1209–1218

19. Telles JCF (1983) The boundary element method applied to inelastic problems. Springer, Berlin Heidelberg New York

20. Yu G, Mansur WJ, Carrer JAM (1999) The linear θ method for 2D elastodynamic BE analysis. Comput Mech 24:82–89

21. Yu G, Mansur WJ, Carrer JAM, Gong L (1998) A linear θ method applied to 2D time-domain BEM. Commun Numer Meth Eng 14:1171–1179

22. Yu G, Mansur WJ, Carrer JAM, Gong L (1998) Time weight-
ing in time domain BEM. Eng Anal Bound Elem 22:175–181

23. Yu G, Mansur WJ, Carrer JAM, Gong L (2000) Stability of Galerkin and collocation time domain boundary element methods as applied to the scalar wave equation. Comput Struct 74:495–506

24. Yu G, Mansur WJ, Carrer JAM, Lie ST (2001) A more stable scheme for BEM/FEM coupling applied to two-dimensional elastodynamics. Comput Struct 79:811–823