# Solution of the two-dimensional scalar wave equation by the time-domain boundary element method: Lagrange truncation strategy in time integration

Jose Antonio Marques Carrer, Webe João Mansur

► To cite this version:

Jose Antonio Marques Carrer, Webe João Mansur. Solution of the two-dimensional scalar wave equation by the time-domain boundary element method: Lagrange truncation strategy in time integration. Structural Engineering and Mechanics, 2006, 23 (3), pp.263-278. 10.12989/sem.2006.23.3.263. hal-01352433

HAL Id: hal-01352433

https://hal.science/hal-01352433

Submitted on 7 Aug 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
# Solution of the two-dimensional scalar wave equation by the time-domain boundary element method: Lagrange truncation strategy in time integration

J. A. M. Carrer

Programa de Pós-Graduação em Métodos Numéricos em Engenharia, Universidade Federal do Paraná,
Caixa Postal 19011, CEP 81531-990, Curitiba, PR, Brasil

W. J. Mansur

Programa de Engenharia Civil, COPPE/UFRJ, Universidade Federal do Rio de Janeiro,
Caixa Postal 68506, CEP 21945-970, Rio de Janeiro, Brasil

**Abstract.** This work presents a time-truncation scheme, based on the Lagrange interpolation polynomial, for the solution of the two-dimensional scalar wave problem by the time-domain boundary element method. The aim is to reduce the number of stored matrices, due to the convolution integral performed from the initial time to the current time, and to keep a compromise between computational economy and efficiency and the numerical accuracy. In order to verify the accuracy of the proposed formulation, three examples are presented and discussed at the end of the article.

**Keywords:** truncation strategy; boundary element method; scalar wave equation; TD-BEM.

## 1. Introduction

The Boundary Element Method (BEM) has been applied to solve time-dependent problems quite successfully, as demonstrated by the several works, dealing with different approaches, published during the last years. For general purposes, these approaches can be classified according to the nature of the fundamental solution adopted. The use of time-dependent fundamental solution originates time-domain formulations (TD-BEM). TD-BEM formulations, beside providing good representation of causality and time response jumps and, consequently, leading to very accurate results, fulfill the radiation condition, which makes them suitable for infinite domain analysis, e.g. (Mansur 1983, Dominguez 1993, Mansur *et al.* 1998). The use of static fundamental solution, on the other hand, originate two formulations, classified according to the maintenance or not of the inertial domain integral in the BEM equations as: D-BEM, that keeps the domain integral (D means domain) in the equations (Carrer and Telles 1992, Hatzigeorgiou and Beskos 2001), and DR-BEM
(DR means Dual Reciprocity) that, by means of suitable interpolation functions, substitutes the
domain integral by boundary integrals (Kontoni and Beskos 1993, Partridge et al. 1992).
Formulations based on frequency and Laplace domains are also available, e.g., (Manolis 1983, de
Lacerda et al. 1996). More recently, a formulation based on the Operational Quadrature Method
appeared in the literature (Gaul and Schanz 1999, Schanz 2001).

This work is concerned with the solution of 2D scalar wave problem and is based on the TD-BEM formulation. The aim here is to reduce the storage computational cost due to the evaluation of the convolution integral that appear in the formulation. In other words, the aim is to reduce the number of assembled matrices, necessary to take into account the time history contribution, and to preserve the accuracy of the standard TD-BEM formulation. It seems that the best strategy to achieve this goal is to truncate the time integration, as previously presented (Demirel and Wang 1987, Mansur and de Lima-Silva 1992, Soares Jr. and Mansur 2004): the present work is concerned with this topic. Here, the whole time interval of analysis ($0 \le t \le t_n$) is divided into two parts, both constituted by time steps, $\Delta t$, of equal size: the first one is restricted to the interval $t_k \le t \le t_{k+1}$, where $t_k$ is a specific value of time. In this first interval, the time integration is effectively done, i.e., the matrices are assembled in the standard way. In the second interval, defined by $0 \le t \le t_k$, discrete values of time are chosen and a Lagrange polynomial, passing through these discrete values of time, is constructed. Proceeding in this way, only matrices at these specific values of time are appropriately assembled, and matrices corresponding to intermediate values are computed by interpolation. The numerical computation of the convolution integral is carried out in this way. In order to simplify the nomenclature, the first interval, constituted by the last $n_{INT}$ time intervals, will be referred to as *integration interval*; the second interval, by its turn, will be referred to as *interpolation interval*. Naturally, the length of the so-called *integration interval* and the degree of the polynomial in the *interpolation interval* are problem dependent parameters; incorrect choices of these parameters can lead to not reliable results.

Linear boundary elements and linear triangular cells were employed, respectively, to approximate
the boundary and the part of the domain with non-homogeneous initial conditions (note that the
solution of problems with non-homogeneous initial conditions does not present any difficulty: the
domain integrals, related to the initial conditions are computed as in the standard TD-BEM
formulation). Linear and constant time variation were assumed, respectively, for the potential and its
normal derivative (flux) and, as usual in TD-BEM formulations, time integration was carried out
analytically.

Three examples are presented and discussed at the end of the article, in order to verify the
applicability of the proposed interpolation scheme.

## 2. TD-BEM formulation

The TD-BEM equations can be written by employing the kernel regularization procedure (Mansur 1983) or the concept of finite part of integrals (FPI) (Hadamard 1952). If, as it is usual in TD-BEM formulations, time integration is performed analytically, the resulting time integrated kernels from both representations are the same (Mansur and Carrer 1993), that is, both representations are entirely equivalent. The use of the Hadamard's concept, however, leads to more compact expressions and, it is the authors' opinion, provides an elegant representation of the equations involved in the analysis. A brief summary of the equations is given below.
<figure><img src="image_2.png" /><figcaption>Fig. 1 Internal angle for computation of c (ξ)</figcaption></figure>

## 2.1 Basic integral equation

Time-domain integral representation of the 2-D scalar wave propagation problem is written as (Carrer and Mansur 1996):

$$ 
\begin{aligned}
4 \pi c(\xi) u(\xi, t) = & \int_{\Gamma} \int_{0}^{t^{+}} u^{*}(X, t; \xi, \tau) p(X, \tau) d\tau d\Gamma(X) - \int_{\Gamma} \frac{\partial r}{\partial n} \int_{0}^{t^{+}} u_{r}^{*}(X, t; \xi, \tau) u(X, \tau) d\tau d\Gamma(X) - \\
& \frac{1}{c} \int_{\Gamma} \frac{\partial r}{\partial n} u_{o}^{*}(X, t; \xi) u_{o}(X) d\Gamma(X) + \frac{1}{c^2} \int_{\Omega} u_{o}^{*}(X, t; \xi) v_{o}(X) d\Omega(X) + \\
& \frac{1}{c} \int_{\Omega} \left[ \frac{u_{o}^{*}(X, t; \xi)}{r} - b_{o}^{*}(X, t; \xi) \right] u_{o}(X) d\Omega(X) + \frac{1}{c} \int_{\Omega} u_{o}^{*}(X, t; \xi) \frac{\partial u_{o}(X)}{\partial r} d\Omega(X)
\end{aligned} 
\quad (1) $$

In Eq. (1), $\Gamma$ is the boundary; $\Omega$ is the domain, or the part of the domain, that presents non-homogeneous initial conditions; the coefficient $c(\xi)$ assumes the same values of the static case, i.e., it is equal to 1 ($\xi \in \Omega$) or ($\alpha/2\pi$) ($\xi \in \Gamma$ and $\alpha$ is the internal angle depicted in Fig. 1); and the subscript 'o' means that $\tau=0$.

The fundamental solution, $u^*(X, t; \xi, \tau)$, that corresponds to the effect of a source represented by an impulse at $t = \tau$ located at $X = \xi$ propagating with velocity equal to $c$, has the following expression:

$$ u^*(X, t; \xi, \tau) = U^*(X, t; \xi, \tau) H[c(t - \tau) - r] \quad (2) $$

where:

$$ U^*(X, t; \xi, \tau) = \frac{2c}{\sqrt{c^2(t-\tau)^2 - r^2}} \quad (3) $$

In expression (2), $H[c(t - \tau) - r]$ stands for the Heaviside function ($r = r(X; \xi)$) is the distance between the field ($X$) and the source ($\xi$) points).

The functions $b^*(X, t; \xi, \tau)$ and $u_r^*(X, t; \xi, \tau)$, in Eq. (1), are given by:

$$ b^*(X, t; \xi, \tau) = B^*(X, t; \xi, \tau) H[c(t - \tau) - r] \quad (4) $$

$$ u_r^*(X, t; \xi, \tau) = U_r^*(X, t; \xi, \tau) H[c(t - \tau) - r] \quad (5) $$
where:

$$B^*(X, t; \xi, \tau) = \frac{2c[c(t-\tau)-r]}{[c^2(t-\tau)^2 - r^2]^{3/2}} \quad (6)$$

and

$$U_r^*(X, t; \xi, \tau) = \frac{2cr}{[c^2(t-\tau)^2 - r^2]^{3/2}} \quad (7)$$

The symbol $\mathfrak{F}$ on the second term on the right-hand-side of Eq. (1) stands for the FPI (Hadamard 1952):

$$\mathfrak{F}_0^t u_r^*(X, t; \xi, \tau)u(X, \tau)d\tau = \lim_{\tau \to t-r/c} \left\{ \int_0^\tau U_r^*(X, t; \xi, \tau)u(X, \tau)d\tau - \frac{1}{c}U^*(X, t; \xi, \tau)u(X, \tau) \right\} \quad (8)$$

## 2.2 Space derivative boundary integral equation for internal points

The derivative of Eq. (1) with respect to a generic direction $m(\xi)$, if $\xi \in \Omega$, can be written as (Carrer and Mansur 1996):

$$\begin{aligned} 4\pi \frac{\partial u(\xi, t)}{\partial m(\xi)} = & -\int_{\Gamma} \mathfrak{F}_0^t u_r^*(X, t; \xi, \tau)p(X, \tau)d\tau(\mathbf{r}^o \cdot \mathbf{m}^o)d\Gamma(X) \\ & -\int_{\Gamma} \mathfrak{F}_0^t u_r^*(X, t; \xi, \tau)u(X, \tau)d\tau\frac{1}{r}[(\mathbf{r}^o \cdot \mathbf{m}^o)(\mathbf{r}^o \cdot \mathbf{n}^o) - (\mathbf{m}^o \cdot \mathbf{n}^o)]d\Gamma(X) \\ & +\int_{\Gamma} \frac{\partial}{\partial r}(\mathfrak{F}_0^t u_r^*(X, t; \xi, \tau)u(X, \tau)d\tau)(\mathbf{r}^o \cdot \mathbf{m}^o)(\mathbf{r}^o \cdot \mathbf{n}^o)d\Gamma(X) \\ & -\frac{1}{c}\int_{\Gamma} u_o^*(X, t; \xi)u_o(X)\frac{1}{r}[(\mathbf{r}^o \cdot \mathbf{m}^o)(\mathbf{r}^o \cdot \mathbf{n}^o) - (\mathbf{m}^o \cdot \mathbf{n}^o)]d\Gamma(X) \\ & +\frac{1}{c}\int_{\Gamma} \frac{\partial u_o^*(X, t; \xi)}{\partial r}u_o(X)(\mathbf{r}^o \cdot \mathbf{m}^o)(\mathbf{r}^o \cdot \mathbf{n}^o)d\Gamma(X) + \frac{\partial}{\partial m(\xi)}(I_{\Omega, v_o} + I_{\Omega, u_o}) \end{aligned} \quad (9)$$

The FPI in the first term on the right-hand-side of Eq. (9) is interpreted as:

$$\mathfrak{F}_0^t u_r^*(X, t; \xi, \tau)p(X, \tau)d\tau = \lim_{\tau \to t-r/c} \left\{ \int_0^\tau U_r^*(X, t; \xi, \tau)p(X, \tau)d\tau + U^*(X, t; \xi, \tau)p(X, \tau) \right\} \quad (10)$$

The derivative with respect to $r$ of the FPI indicated on the third term on the right-hand-side of Eq. (9) is defined as follows:

$$\frac{\partial}{\partial r}(\mathfrak{F}_0^t u_r^*(X, t; \xi, \tau)u(X, \tau)d\tau) = \lim_{\tau \to t-r/c} \left\{ \int_0^\tau \frac{\partial U_r^*(X, t; \xi, \tau)}{\partial r}u(X, \tau)d\tau - \frac{2}{c}U_r^*(X, t; \xi, \tau)u(X, \tau) - \frac{1}{c}U^*(X, t; \xi, \tau)\frac{\partial u(X, \tau)}{\partial \tau} \right. \\ \left. + \frac{1}{c^2}\left[U^*(X, t; \xi, \tau)\frac{\partial u(X, \tau)}{\partial \tau} + U_\tau^*(X, t; \xi, \tau)u(X, \tau)\right] \right\} \quad (11)$$
where:

$$ U_{\tau}^{*}(X, t; \xi, \tau) = \frac{2c^{3}(t-\tau)}{[c^{2}(t-\tau)^{2} - r^{2}]^{3/2}} \quad (12) $$

## 2.3 Time derivative boundary integral equation for internal points

The derivative of Eq. (1) with respect to time, if $\xi \in \Omega$, can be written as (Carrer and Mansur 1996):

$$ 
\begin{aligned}
4\pi \frac{\partial u(\xi, t)}{\partial t} = & \int_{\Gamma} \oint_{0}^{t} u_{t}^{*}(X, t; \xi, \tau) p(X, \tau) d\tau d\Gamma(X) \\
& - \int_{\Gamma} \frac{\partial}{\partial t} \left( \oint_{0}^{t} u_{r}^{*}(X, t; \xi, \tau) u(X, \tau) d\tau \right) \frac{\partial r}{\partial n} d\Gamma(X) \\
& - \frac{1}{c} \int_{\Gamma} \frac{\partial r}{\partial n} \frac{\partial u_{o}^{*}(X, t; \xi)}{\partial t} u_{o}(X) d\Gamma(X) \\
& + \frac{\partial}{\partial t} (I_{\Omega, v_o} + I_{\Omega, u_o})
\end{aligned}
\quad (13) $$

The FPI in the first term on the right-hand-side of Eq. (13) is interpreted as:

$$ \oint_{0}^{t} u_{t}^{*}(X, t; \xi, \tau) p(X, \tau) d\tau = \lim_{\tau \to t-r/c} \left\{ \int_{0}^{\tau} U_{t}^{*}(X, t; \xi, \tau) p(X, \tau) d\tau + U^{*}(X, t; \xi, \tau) p(X, \tau) \right\} \quad (14) $$

The function $u_t^*(X, t; \xi, \tau)$ in expression (14) is given by:

$$ u_{t}^{*}(X, t; \xi, \tau) = U_{t}^{*}(X, t; \xi, \tau) H[c(t - \tau) - r] \quad (15) $$

where:

$$ U_{t}^{*}(X, t; \xi, \tau) = - \frac{2c^{3}(t-\tau)}{[c^{2}(t-\tau)^{2} - r^{2}]^{3/2}} \quad (16) $$

The time derivative of the FPI indicated on the second term on the right-hand-side of Eq. (13) is defined as follows:

$$ \lim_{\tau \to t-r/c} \left\{ \int_{0}^{\tau} \frac{\partial U_{r}^{*}(X, t; \xi, \tau)}{\partial t} u(X, \tau) d\tau + U_{r}^{*}(X, t; \xi, \tau) u(X, \tau) - \frac{1}{c} U^{*}(X, t; \xi, \tau) \frac{\partial u(X, \tau)}{\partial \tau} \right\} \quad (17) $$

The space and time derivatives of the initial conditions domain integrals were only indicated in Eqs. (9) and (13). For a detailed discussion concerning domain integration, the reader is referred to (Mansur 1983, Carrer and Mansur 1996).
### 3. Numerical procedure

For the numerical solution of Eq. (1) the boundary is approximated by linear boundary elements, and the domain (or the part of it in which non-homogeneous initial conditions appear) is approximated by linear triangular cells. Time approximation assumes linear variation for the potential and constant variation for the potential normal derivative (flux), along the time steps $\Delta t$ in which the overall time of analysis, $0 \le t \le t_n$, is divided. Time integrals are computed analytically and the remaining kernels, in the boundary integrals, are computed numerically by the Gaussian quadrature.

The application of the discretized version of the boundary integral equation to all boundary nodes produces a system of algebraic equations, which can be written according to:

$$ C u^n + \sum_{m=0}^{n} H^{nm} u^m = \sum_{m=0}^{n} G^{nm} p^m + F^n \quad (18) $$

In Eq. (18), diagonal matrix C is constituted by the $c(\xi)$ coefficients, matrices $H^{nm}$ and $G^{nm}$ result from the spatial integration of the time-integrated kernels related, respectively, to $u_r^*(X, t; \xi, \tau)$ and to $u^*(X, t; \xi, \tau)$, and the vector $F^n$ contains the initial conditions contributions (see Eq. (1)). Note, additionally, that the subscripts $n$ and $m$ stand for the time $t_n$ (final time) and $t_m$ (previous times), respectively.

After imposing the boundary conditions, the system of equations represented by Eq. (18) can be solved for the boundary unknowns.

For a more detailed discussion concerning these matters, the reader is referred to (Mansur 1983, Dominguez 1993, Carrer and Mansur 1996).

Note that equations similar to Eq. (18) are obtained for the solution of Eqs. (9) and (13) and that, in these equations, one has $C = I$.

### 4. Lagrange interpolation polynomial for time truncation

According to Eq. (1), or to its corresponding discretized version given by Eq. (18), it is necessary to take into account the contribution of the history, i.e., it is necessary to take into account the contribution of all responses previous to $t_n$ to obtain the response at time $t_n$. If very large values of $n$ are required, i.e., if late time results are required, and if there is no memory available in the computer, the use of time truncation procedures becomes desirable and justified. Naturally, such a procedure introduces approximations in overall time integration: these approximations depend on the length of the intervals in which the time integration is effectively done and the interpolation takes place. The question that arises is: What is the best strategy for an efficient time truncation? It is the authors' opinion that the response to this question is not conclusive, yet. Truncation procedures, for infinite domain applications, were first reported by Demirel and Wang (1987): in this work, it is assumed that there is a time, say $t_s$, before which the contributions to the response at a late time $t_n, t_n \gg t_s$, can be disregarded; in this way, the convolution integrals in Eqs. (1), (9) and (13) start at $\tau=t_s$ instead of at $\tau=0$. Another truncation scheme, which applies quite well to bounded domains analyses, was presented by (Mansur and de Lima-Silva 1992): now the early time contributions are not disregarded, but the time integrals are computed in a simplified way in the interval $[0, t_s]$. Another work, based on multi-linear and Chebyshev-Lagrange polynomials, was
<figure><img src="image_3.png" /><figcaption>Fig. 2 Definition of interpolation and integration intervals</figcaption></figure>

recently presented (Soares Jr. and Mansur 2004). In order to provide topics for discussion, this work presents a time truncation strategy based on the use of Lagrange interpolation polynomials. The basic idea can be outlined as follows:

i) initially, the final time of analysis, say $t_n$, is divided in time intervals $\Delta t$, such that $t_n = n\Delta t$;

ii) the number of time intervals in which the time integrations will be carried out in the standard way, $n_{INT}$, is then defined (or specified as a input datum in the computer program). Note that $n_{INT}$ defines an interval named here as *integration interval* and that the convolution, diversely from the standard TD-BEM formulation, is computed from an intermediate value of time to the current (final) time;

iii) in the remaining interval, called *interpolation interval*, lgr discrete values of time are selected and the corresponding matrices are appropriately assembled. Then, the matrices corresponding to intermediate values of time are computed by employing a Lagrange interpolation polynomial. Therefore, if the interpolation interval is constituted by the first $k$ time intervals $\Delta t$, $0 \le t \le k\Delta t = t_k$, Eq. (18) can be written as follows (the same is valid for the discretized versions of Eqs. (9) and (13)):

$$ C\mathbf{u}^n + \sum_{m=k}^{n} \mathbf{H}^{nm}\mathbf{u}^m + \sum_{m=0}^{k-1} \mathbf{L}_{h}^{nm}\mathbf{u}^m = \sum_{m=k}^{n} \mathbf{G}^{nm}\mathbf{p}^m + \sum_{m=0}^{k-1} \mathbf{L}_{g}^{nm}\mathbf{p}^m + \mathbf{F}^n \quad (19) $$

In Eq. (19), the interpolated matrices associated to matrices **H** and **G** are denoted by $L_h$ and $L_g$, respectively. It is important to mention that necessarily $t_o$ and $t_k$ must belong to the set of selected discrete time values. As a matter of fact, $t_o$ is always the first value and $t_k$ is always the last value of the *interpolation interval*; intermediate values can be equally spaced between them, or not. An illustration of the scheme is shown in Fig. 2.

A generic interpolated matrix, say **B**, assumed to be a function of time, can be represented as:

$$ \mathbf{B}(t) = \sum_{m=0}^{lgr} \mathbf{B}^{nm} L_m(t) \quad (20) $$

in which $\mathbf{B}^{nm}$ represent the matrices computed appropriately at the *lgr* selected discrete time values, for $t = t_n$. The Lagrange interpolation polynomial can be defined according to:

$$ L_m(t) = \frac{(t - \bar{t}_o)(t - \bar{t}_1)(t - \bar{t}_2)\dots(t - \bar{t}_{m-1})(t - \bar{t}_{m+1})\dots(t - \bar{t}_{lgr})}{(\bar{t}_m - \bar{t}_o)(\bar{t}_m - \bar{t}_1)(\bar{t}_m - \bar{t}_2)\dots(\bar{t}_m - \bar{t}_{m-1})(\bar{t}_m - \bar{t}_{m+1})\dots(\bar{t}_m - \bar{t}_{lgr})}; \quad \bar{t}_o \le t \le \bar{t}_{lgr} \quad (21) $$

and has the property:

$$ L_m(\bar{t}_m) = 1 \quad \text{and} \quad L_m(\bar{t}_j) = 0 \quad \text{if} \quad j \neq m \quad (22) $$
It is important to point out that $\bar{t}_j$ does not represent $t_j = j\Delta t$ but, instead of it, $\bar{t}_j$ represents the j-th selected discrete time value and that, according to what was previously mentioned, $\bar{t}_o = t_o$ and $\bar{t}_{lgr} = t_k$.

Another aspect to be mentioned is the generality of the proposed procedure: it can easily be applied to alternative TD-BEM formulations (Yu et al. 1998, Carrer and Mansur 2002).

# 5. Examples

In the examples presented in this work, reference will be made to the dimensionless $\beta$ parameter:

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

in which $\ell$ is the boundary element length. The choice of the value of the $\beta$ parameter is a problem dependent task: as a general rule, small and large values, inside the interval $0 < \beta < 1$, must be avoided.

The following notation will be employed in the examples: $n$ represents the total number of time intervals and $n_{INT}$ represents the number of time intervals that constitute the integration interval; besides, equally spaced time values were adopted to construct the Lagrange polynomial in the interpolation interval. Along the discussion, the proposed formulation, for simplicity, will be referred to as Lagrange formulation.

Additionally, in the first and in the second example, $E$ is the Young's modulus.

## 5.1 One-dimensional rod under compression (waveguide)

This example consists of a one-dimensional rod fixed at one extreme and free at the other, that is subjected to a compression load suddenly applied at $t=0$ and kept constant during the analysis, see Fig. 3. The material is such that $c=1$. The boundary discretization employed 24 elements, as depicted in Fig. 4. The number of time intervals is $n=320$, and the time interval length was defined for $\beta=0.6$. Results furnished by the standard TD-BEM formulation, corresponding to the potential at node A(a, b/2) and to the flux at node B(0, b/2), are presented in Figs. 5 and 6, respectively, and were included to demonstrate how the proposed interpolation scheme can produce reliable results,

<figure><img src="image_4.png" /><figcaption>Fig. 3 One-dimensional rod: geometry and loading definition</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 One-dimensional rod: boundary discretization and selected nodes</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5 One-dimensional rod: potential at boundary node A(a, b/2): standard TD-BEM</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 7 One-dimensional rod: potential at boundary node A(a, b/2): Lagrange formulation with $n = 320$, $n_{INT} = 40$, and 14th order polynomial</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6 One-dimensional rod: flux at boundary node B(0, b/2): standard TD-BEM</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8 One-dimensional rod: flux at boundary node B(0, b/2): Lagrange formulation with $n = 320$, $n_{INT} = 40$, and 14th order polynomial</figcaption></figure>

even at a reduced storage cost: the results presented here were obtained by adopting $n_{INT} = 40$ and a Lagrange interpolation polynomial of 14th order. In the interpolation interval, the discrete time values, say $\bar{t}_j$, are equally spaced, i.e., $j = 0, 20, 40, ..., 280$, always bearing in mind that $\bar{t}_j = j\Delta t$. It is important to mention that only $2(41+14) = 110$ matrices were assembled and stored, instead of the $(321+320) = 641$ matrices required by the standard TD-BEM formulation; in other words,
<figure><img src="image_10.png" /><figcaption>Fig. 9 One-dimensional rod: potential at point C(a/2, b/2): Lagrange formulation with $n = 320$, $n_{INT} = 40$, and 14th order polynomial</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10 One-dimensional rod: potential space derivative $\partial u/\partial x$ at point C(a/2, b/2): Lagrange formulation with $n = 320$, $n_{INT} = 40$, and 14th order polynomial</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11 One-dimensional rod: potential time derivative $\partial u/\partial t$ at point C(a/2, b/2): Lagrange formulation with $n = 320$, $n_{INT} = 40$, and 14th order polynomial</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12 One-dimensional rod: potential at boundary node A(a, b/2): Lagrange formulation with $n = 320$, $n_{INT} = 40$, and 7th order polynomial</figcaption></figure>

the proposed scheme reduces the number of the stored matrices to 17.2% of the total number required by the standard TD-BEM. Results for the boundary nodes A(a, b/2) and B(0, b/2) are presented in Figs. 7 and 8, respectively. Results for the potential, and its space and time derivatives
at internal point C(a/2, b/2) are presented in Figs. 9, 10 and 11. Note that, from Fig. 5 to Fig. 11, the BEM results are always compared with the corresponding analytical solution. The oscillations around the analytical solution, in Figs. 8, 10 and 11, appear even in the results from the standard TD-BEM formulation, and better results can be obtained only with the use of a more refined boundary element mesh. It is important to mention that several analyses have been performed by the authors, aiming at finding a pattern for the adoption of the best length for the *integration interval* and the order of the interpolation polynomial. The conclusion is that the choice of these parameters depends on the problem; thus, the experience of the analyst plays an important role in this matter. To illustrate this, another analysis was included; see Fig. 12, in which results not so accurate (when compared to those from Fig. 7) were achieved for the potential at node A(a, b/2): although the integration interval was kept the same, $n_{INT} = 40$, a poor interpolation polynomial of 7th order was adopted: the values $\bar{t}_j$ are equally spaced but, now, one has $j = 0, 40, 80, 120..., 280$. Now, the number of matrices assembled and stored is equal to $2(41+7)=96$, which means that only 15% of the total number of matrices are effectively stored.

## 5.2 Circular cavity in an infinite medium

This example consists of a circular cavity, in an infinite medium, subjected to an internal pressure suddenly applied at $t=0$ and kept constant during the time, see Fig. 13. The boundary discretization, shown in Fig. 14, employed 24 elements. The number of time intervals is $n=380$, and the time interval length was defined for $\beta=0.645$. The results presented here were obtained by adopting $n_{INT}=80$ and a Lagrange interpolation polynomial of 15th order, with discrete time values, $\bar{t}_j$, equally spaced, i.e., $j = 0, 20, 40, ..., 300$. In this example, the advantage of the TD-BEM formulation becomes more evident: only the cavity boundary needs to be approximated and the results at any internal point, no matter how far it is from the cavity centre, can be computed without any domain discretization. The results furnished by the proposed formulation are compared with those furnished by the standard TD-BEM formulation in Fig. 15, for the potential at node A(R, 0), and in Figs. 16 to 18 for the potential and its space and time derivatives at internal point B(2R,0). The results related to the potential, in Figs. 15 and 16, are in very good agreement. The results related to space and time derivatives, furnished by both formulations, present oscillations around

<figure><img src="image_14.png" /><figcaption>Fig. 13 Circular cavity: geometry and loading definition</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14 Circular cavity: boundary discretization and selected nodes</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 15 Circular cavity: potential at boundary node A(R, 0): Lagrange formulation with $n = 380$, $n_{INT} = 80$, and 15th order polynomial</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16 Circular cavity: potential at point B(2R, 0): Lagrange formulation with $n = 380$, $n_{INT} = 80$, and 15th order polynomial</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17 Circular cavity: potential space derivative $\partial u/\partial x$ at point B(2R, 0): Lagrange formulation with $n = 380$, $n_{INT} = 80$, and 15th order polynomial</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18 Circular cavity: potential time derivative $\partial u/\partial t$ at point B(2R, 0): Lagrange formulation with $n = 380$, $n_{INT} = 80$, and 15th order polynomial</figcaption></figure>

their asymptotic values: in the results from the Lagrange formulation, the oscillation is more pronounced, but are still in agreement with the corresponding ones from the standard formulation. In this example, the number of matrices required by the proposed scheme is given by $2(81 + 15) = 192$; comparing this value with the number of matrices required by the standard TD-
BEM formulation, given by $(381 + 380) = 761$, one can verify that the analysis was performed by employing only 25.2% of the matrices originally required.

### 5.3 Square membrane under prescribed initial velocity

The square membrane depicted in Fig. 19, with initial velocity field $v_o = c$ prescribed over the sub-domain $\Omega_o$ and with zero displacements prescribed all over the boundary, is analysed in this

<figure><img src="image_20.png" /><figcaption>Fig. 19 Square membrane: geometry and loading definition</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20 Square membrane: boundary and domain discretization and selected nodes</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 21 Square membrane: flux at boundary node $A(a, a/2)$: Lagrange formulation with $n = 200$, $n_{INT} = 60$, and 14th order polynomial</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 22 Square membrane: potential at point $B(4a/5, a/2)$: Lagrange formulation with $n = 200$, $n_{INT} = 60$, and 14th order polynomial</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. 23 Square membrane: potential space derivative $\partial u/\partial x$ at point B(4a/5, a/2): Lagrange formulation with $n = 200$, $n_{INT} = 60$, and 14th order polynomial</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 24 Square membrane: potential time derivative $\partial u/\partial t$ at point B(4a/5, a/2): Lagrange formulation with $n = 200$, $n_{INT} = 60$, and 14th order polynomial</figcaption></figure>

example. The boundary discretization employed 32 elements and $\Omega_0$ was divided into four cells, as shown in Fig. 20. As previously pointed out (Mansur 1983), the adoption of $\beta < 0.6$ does not introduce any great amount of noise in the BEM results. For this reason, the value $\beta = 0.2$ was adopted. Besides, as the number of discrete values increases, a better picture of the results can be inferred from the numerical results. Fig. 21 presents the results for the flux at boundary node A(a, a/2). Potential, space and time derivatives results, for internal point B(4a/5, a/2), are presented in Figs. 22, 23 and 24, respectively. This analysis was carried out by adopting $n = 200$, $n_{INT} = 60$, and an interpolation polynomial of 14th order, with discrete time values, $\bar{t}_j$, equally spaced, i.e., $j = 0, 10, 20, ..., 140$. In this example, the number of matrices required by the proposed scheme is equal to $2(61+14)=150$ whereas the number of matrices required by the standard TD-BEM formulation is $(201+200)=401$; consequently, the analysis was performed by employing only 37.4% of the matrices originally required. Finally, one can conclude that the numerical results agree quite well with the analytical solution (Morse and Ingard 1968).

## 6. Conclusions

In this work, a strategy for the solution of the 2-D scalar wave propagation problem, by the TD-BEM formulation, is developed. The aim is to reduce the computational costs from the assemblage and from the storage of the matrices related to the time-history contributions to the results at a specific value of time. These matrices are due to the convolution integral presented in the basic TD-BEM formulation. Computational costs are reduced by partially computing the convolution integral, i.e., the time integration is no longer performed from $t_o$ to $t_n$ but, instead of it, from some value $t_k$ to $t_n$; the interval $[t_k, t_n]$ is designated *integration interval*, meaning that the matrices are appropriately
computed there (BEM matrices). In the remaining interval, $[t_o, t_k]$, designated *interpolation interval*, the matrices are computed by interpolation. To do so, a Lagrange polynomial is constructed by selecting discrete *lgr* values of time in the interval $[t_o, t_k]:$ BEM matrices are computed for these *lgr* values and, finally, matrices corresponding to values of time different from the selected ones are computed by interpolation. It is important to mention that the Lagrange interpolation formulation was also adopted for the computation of space and time derivatives of the potential at internal points and for the analysis of problems with non-homogeneous initial conditions. For the chosen examples, the numerical results can be considered good. As the experience plays an important role in the choice of the new parameters, i.e., the *integration interval* length and the order of the interpolation polynomial, in the absence of a study concerning error estimation, two practical recommendations are suggested by the authors: i) the use of interpolation polynomials of order greater than or equal to 10, with equally spaced time values; ii) the ratio between $n_{INT}$ and *n* must belong to the interval: $0.15 \le n_{INT}/n \le 0.30$. Note that the proposed procedure can easily be extended to elastodynamics, as well.

## References

Carrer, J.A.M. and Mansur, W.J. (1996), “Time-domain BEM analysis for the 2D scalar wave equation: Initial conditions contributions to space and time derivatives”, *Int. J. Numer. Meth. Eng.*, **39**, 2167-2188.

Carrer, J.A.M. and Mansur, W.J. (2002), “Time-dependent fundamental solution generated by a not impulsive source in the boundary element method analysis of the 2D scalar wave equation”, *Commun. Numer. Meth. Eng.*, **18**, 277-285.

Carrer, J.A.M. and Telles, J.C.F. (1992), “A boundary element formulation to solve transient dynamic elastoplastic problems”, *Comput. Struct.*, **45**, 707-713.

Demirel, V. and Wang, S. (1987), “Efficient boundary element method for two-dimensional transient wave propagation problems”, *Appl. Mathematical Modelling*, **11**, 411-416.

Dominguez, J. (1993), *Boundary Elements in Dynamics*, Computational Mechanics Publications, Southampton, Boston.

Gaul, L. and Schanz, M. (1999), “A comparative study of three boundary element approaches to calculate the transient response of viscoelastic solids with unbounded domains”, *Comput. Method Appl. Mech. Eng.*, **179**, 111-123.

Hadamard, J. (1952), *Lectures on Cauchy's Problem in Linear Partial Differential Equations*, Dover Publications, New York.

Hatzigeorgiou, G.D. and Beskos, D.E. (2001), “Transient dynamic response of 3-D elastoplastic structures by the D/BEM”, *Proc. of the XXIII International Conf. on the Boundary Element Method*, (eds. D.E. Beskos, C.A. Brebbia, J.T. Katsikadelis, G.D. Manolis), Lemnos, Greece.

Kontoni, D.P.N. and Beskos, D.E. (1993), “Transient dynamic elastoplastic analysis by the dual reciprocity BEM”, *Engineering Analysis with Boundary Elements*, **12**, 1-16.

de Lacerda, L.A., Wrobel, L.C. and Mansur, W.J. (1996), “A boundary integral formulation for two-dimensional acoustic radiation in a subsonic uniform flow”, *J. of the Acoustical Society of America*, **100**, 98-107.

Manolis, G.D. (1983), “A comparative study on three boundary element method approaches to problems in elastodynamics”, *Int. J. Numer. Meth. Eng.*, **19**, 73-91.

Mansur, W.J. (1983), “A time-stepping technique to solve wave propagation problems using the boundary element method”, Ph.D. Thesis, University of Southampton, England.

Mansur, J.W. and Carrer, J.A.M. (1993), “Two-dimensional transient BEM analysis for the scalar wave equation: Kernels”, *Engineering Analysis with Boundary Elements*, **12**, 283-288.

Mansur, W.J., Carrer, J.A.M. and Siqueira, E.F.N. (1998), “Time discontinuous linear traction approximation in time domain BEM scalar wave propagation analysis”, *Int. J. Numer. Meth. Eng.*, **42**, 667-683.
Mansur, W.J. and de Lima-Silva, W. (1992), "Efficient time truncation in two-dimensional BEM analysis of transient wave propagation problems", *Earthq. Eng. Struct. Dyn.*, **21**, 51-63.

Morse, P.M. and Ingard, K.V. (1968), *Theoretical Acoustics*, McGraw-Hill, London.

Partridge, P.W., Brebbia, C.A. and Wrobel, L.C. (1992), *The Dual Reciprocity Boundary Element Method*, Computational Mechanics Publications, Southampton, Boston.

Schanz, M. (2001), *Wave Propagation in Viscoelastic and Poroelastic Continua*, Lectures Notes in Applied Mechanics, Vol. 2, Springer-Verlag, Berlin, Heidelberg.

Soares Jr., D. and Mansur, W.J. (2004), "Compression of time generated matrices in two-dimensional time-domain elastodynamic BEM analysis", *Int. J. Numer. Meth. Eng.*, **61**, 1209-1218.

Yu, G., Mansur, W.J., Carrer, J.A.M. and Gong, L. (1998), "A linear θ method applied to 2D time-domain BEM", *Commun. Numer. Meth. Eng.*, **14**, 1171-1179.