# A step-by-step approach in the time-domain BEM formulation for the scalar wave equation

Jose Antonio Marques Carrer, Webe João Mansur

► To cite this version:

Jose Antonio Marques Carrer, Webe João Mansur. A step-by-step approach in the time-domain BEM formulation for the scalar wave equation. Structural Engineering and Mechanics, 2007, 27 (6), pp.683-696. 10.12989/sem.2007.27.6.683. hal-01352431

HAL Id: hal-01352431

https://hal.science/hal-01352431

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
# A step-by-step approach in the time-domain BEM formulation for the scalar wave equation

J. A. M. Carrer

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

W. J. Mansur

*Programa de Engenharia Civil, COPPE/UFRJ, Caixa Postal 68506, CEP 21945-970, Rio de Janeiro, RJ, Brasil*

**Abstract.** This article is concerned with the presentation of a time-domain BEM approach applied to the solution of the scalar wave equation for 2D problems. The basic idea is quite simple: the basic variables of the problem at time $t_n$ (potential and flux) are computed with the results related to the potential and to its time derivative at time $t_{n-1}$ playing the role of “initial conditions”. This time-marching scheme needs the computation of the potential and its time derivative at all boundary nodes and internal points, as well as the entire discretization of the domain. The convolution integrals of the standard time-domain BEM formulation, however, are not computed; the matrices assembled, only at the initial time interval, are those related to the potential, flux and to the potential time derivative. Two examples are presented and discussed at the end of the article, in order to verify the accuracy and potentialities of the proposed formulation.

**Keywords:** time-domain BEM; scalar wave equation; step-by-step BEM analysis.

## 1. Introduction

Time-domain analysis by the boundary element method (BEM) is a very attractive and challenging area of research; as a consequence a lot of works have been presented so far, enriching the range of applications of the BEM and the literature concerning numerical methods in engineering. For general purposes, the different BEM approaches can be classified into two main groups: the first group, that employs the static fundamental solutions, and the second group, that employs time-dependent fundamental solutions; transformed-domain (frequency and Laplace) approaches can also be included in this group. A brief discussion concerned with different BEM approaches is carried out in what follows. For a more complete discussion concerning dynamic analysis by the BEM, the reader is referred to (Beskos 1997, 2003).
The use of static fundamental solutions generated two research lines: one that keeps the domain integrals in the equations of the method, designated D-BEM (D meaning domain), e.g. (Carrer and Mansur 2004, Hatzigeorgiou and Beskos 2002), and another one that, by the use of suitable interpolation functions transforms the domain integrals into boundary integrals, designated DR-BEM (DR meaning double reciprocity), e.g. (Kontoni and Beskos 1993, Partridge et al. 1992). As the time is not a variable, i.e., time does not appear explicitly in the integral equations, an approximation to the acceleration is required in order to perform the march in time. As long as the authors know, the Houbolt scheme (Houbolt 1950) was the only approximation successfully employed with the BEM. Alternative time-marching schemes were recently proposed by (Carrer and Mansur 2004, Souza et al. 2004).

Transformed-domain BEM formulations also have been the subject of intense research work and are mainly related either to the frequency-domain approach, (Dominguez 1993, Gaul and Wenzel 2002, Mansur et al. 2004), or to the Laplace-domain approach, e.g. (Gaul and Schanz 1999). When employing these formulations, the problem is initially solved in the transformed-domain for a suitable number of discrete values of the transform parameter and, if required, the solution in the time-domain can be obtained by means of adequate inverse transformation procedures, e.g. (Dubner and Abate 1968, Durbin 1974) for Laplace transforms, and the well known DFT (or FFT) algorithms for Fourier transforms. Another BEM formulation, based on the works by (Lubich 1988a, 1988b), is called *operational quadrature method* or *convolution quadrature method*, e.g. (Schanz and Antes 1997, Schanz 2001, Abreu *et al.* 2003, Antes *et al.* 2004, Abreu *et al.* 2006). The main characteristics of this formulation, that can be applied to problems where time-domain fundamental solutions are not available, are: i) the use of fundamental solutions in the Laplace transformed domain and ii) the numerical approximation of the time convolutions, presented in the time-domain BEM equations, by a quadrature formula based on a linear multistep method, which provides direct solution in the time-domain.

The formulations that employ time-dependent fundamental solutions, designated by TD-BEM (TD meaning time-domain), are very elegant from the mathematical point of view, e.g. (Mansur 1983, Dominguez 1993, Carrer and Mansur 2002). The fulfilment of the radiation condition makes the TD-BEM formulations suitable for infinite domain analyses. Together with this favourable characteristic, the good representation of causality and time response jumps leads to very accurate results. These advantageous characteristics, however, are counterbalanced by the high computational effort required to compute the time convolution integrals that appear in the TD-BEM integral equations. In order to overtake this difficulty, some works concerned with the reduction of the computational cost, by truncation in the time integration, were presented (Demirel and Wang 1987, Mansur and de Lima-Silva 1992, Soares and Mansur 2004, Carrer and Mansur 2006). Focusing the attention in this aspect of the research, this work is concerned with the solution of scalar wave equation for 2D problems by employing a modified version of the TD-BEM formulation. The approach proposed in this work is quite simple and is based on the computation of the basic variables of the problem, i.e. the potential and its normal derivative (flux), at time $t_n$ with the potential and its time derivative at time $t_{n-1}$ treated as the “initial conditions”. Following this approach and taking into account the time translation property of the kernels, e.g. (Mansur 1983), the computation of the convolution integrals, which appear in the standard TD-BEM formulation, is not required. As a matter of fact, as demonstrated along the text, time integration is restricted only to the discrete interval $[0, \Delta t]$. For the development of the proposed approach, it becomes necessary to calculate the derivative, with respect to time, of the potential from the basic BEM integral
equation (potential integral equation), thus generating the velocity integral equation. As the potential and the velocity at time $t_{n-1}$ play the role of “initial conditions” for the potential, for the flux and for the velocity at time $t_n$, the domain integrals, that in the TD-BEM formulation are restricted only to the part of the domain with non-homogeneous initial conditions, are now extended to the entire domain. Besides, by following an approach similar to that from the D-BEM formulation, e.g. (Carrer and Mansur 2004, Hatzigeorgiou and Beskos 2002), potential and velocity integral equations are applied at boundary nodes and internal points simultaneously, thus generating an enlarged system of equations. After solving this system of equations one has the values of the potential and of the flux at boundary nodes and the values of the potential at internal points. The values of the velocity are computed directly, since the boundary variables are already known. One question can be formulated at this moment: can this approach be used for infinite domain analysis? If results at late times are required, the answer is no: in this case the interruption of the domain discretization generates a fictitious boundary and, consequently, the influence of waves reflected there tends to invalidate the results. Although the response to the question is negative, it is the authors' opinion that the present study demonstrates the immense possibilities of research in the TD-BEM. Also, the two examples presented at end of this work demonstrate the potentialities of the proposed formulation and encourage future research in this area.

## 2. Standard TD-BEM integral equations

In what follows, a brief summary of the TD-BEM formulation is presented. For a more complete discussion concerning this topic, the reader is referred to (Mansur 1983, Dominguez 1993).

### 2.1 Potential integral equation

The potential integral equation associated to the TD-BEM formulation of the scalar wave problem, by making use of the concept of finite part of integrals, FPI, (Hadamard 1952) can be written as (Carrer and Mansur 1996, 2006)

$$ 
\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} \Big|_{\Gamma} 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) $$

The following notation was employed in Eq. (1): $\Gamma$ is the boundary; $\Omega$ is the domain, or the part of the domain, that presents non-homogeneous initial conditions $u_o(X)$ and/or $v_o(X) = \frac{\partial u(X, t)}{\partial t} \bigg|_{t=0}$; the variable $p(X, t) = du(X, t)/dn(X)$, where $n(X)$ is the unit outward normal to the boundary at $X$, is the flux; and the coefficient $c(\xi)$ on the left-hand-side of Eq. (1) assumes the same values of the static case, e.g. (Mansur 1983).

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 \text{with} \quad U^*(X, t; \xi, \tau) = \frac{2c}{\sqrt{c^2(t-\tau)^2 - r^2}} \quad (2) $$

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 \text{with} \quad B^*(X, t; \xi, \tau) = \frac{2c[c(t-\tau)-r]}{[c^2(t-\tau)^2 - r^2]^{3/2}} \quad (3) $$

and

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

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

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

## 2.2 Velocity integral equation

The integral equation associated to the potential time derivative can be written as (Carrer and Mansur 1996, 2006)

$$ 4\pi c(\xi) \frac{\partial u(\xi, t)}{\partial t} = \int_{\Gamma} \mathcal{F} \int_{0}^{t^+} u_t^*(X, t; \xi, \tau) p(X, \tau) d\tau d\Gamma(X) - \int_{\Gamma} \frac{\partial}{\partial t} \left( \mathcal{F} \int_{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}) \quad (6) $$

where, for simplicity and because the resulting expressions do not require any special treatment in the domain integration, the time derivative of the domain integrals related to the non-homogeneous initial conditions are indicated by the last term on the right-hand-side of Eq. (6).

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

$$ \mathcal{F} \int_0^t u_t^*(X, t; \xi, \tau) p(X, \tau) d\tau = \lim_{\tau \to t-r/c} \left\{ \int_0^t U_t^*(X, t; \xi, \tau) p(X, \tau) d\tau + U^*(X, t; \xi, \tau) p(X, \tau) \right\} \quad (7) $$

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

$$ u_t^*(X, t; \xi, \tau) = U_t^*(X, t; \xi, \tau)H[c(t-\tau)-r] \quad \text{with} \quad U_t^*(X, t; \xi, \tau) = \frac{2c^3(t-\tau)}{[c^2(t-\tau)^2 - r^2]^{3/2}} \quad (8) $$

The time derivative of the FPI indicated on the second term on the right-hand-side of Eq. (6) is defined as
$$ \frac{\partial}{\partial t} \left( \oint_0^t u_r^*(X, t; \xi, \tau) u(X, \tau) d\tau \right) = $$

$$ \lim_{\tau \to t-r/c} \left\{ \int_0^\tau \frac{\partial U_r^*(X, t; \xi, \tau)}{\partial \tau} 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 (9) $$

### 3. Numerical solution

For the numerical solution of the problem, linear and constant time variations were assumed, respectively, for the potential and for the flux. Following the usual time-domain BEM approach, e.g. (Mansur 1983, Dominguez 1993), time integrals were computed analytically. It is important to mention that time integrals can also be computed numerically; in this case, one should employ the Kutt quadrature formulae (Kutt 1975) to evaluate the finite part integrals that appear in Eq. (1) and in Eq. (6). For the boundary approximation, isoparametric linear elements are employed; the domain is approximated by triangular linear cells, in which the initial conditions vary linearly, e.g. (Carrer and Mansur 2004).

After space and time discretization has been accomplished, the application of Eq. (1) to all boundary nodes and internal points simultaneously (in the next section this procedure becomes clearer) generates the system of equations below

$$ \begin{bmatrix} \mathbf{C} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \left\{ \begin{matrix} \mathbf{u}_n^b \\ \mathbf{u}_n^d \end{matrix} \right\} + \sum_{m=0}^{n} \begin{bmatrix} \mathbf{H}_{nm}^{bb} & \mathbf{0} \\ \mathbf{H}_{nm}^{db} & \mathbf{I} \end{bmatrix} \left\{ \begin{matrix} \mathbf{u}_m^b \\ \mathbf{u}_m^d \end{matrix} \right\} = \sum_{m=0}^{n} \begin{bmatrix} \mathbf{G}_{nm}^{bb} \\ \mathbf{G}_{nm}^{db} \end{bmatrix} \left\{ \mathbf{p}_m^b \right\} + \begin{bmatrix} \mathbf{U}_n^{bb} & \mathbf{U}_n^{bd} \\ \mathbf{U}_n^{db} & \mathbf{U}_n^{dd} \end{bmatrix} \left\{ \begin{matrix} \mathbf{u}_o^b \\ \mathbf{u}_o^d \end{matrix} \right\} + \begin{bmatrix} \mathbf{V}_n^{bb} & \mathbf{V}_n^{bd} \\ \mathbf{V}_n^{db} & \mathbf{V}_n^{dd} \end{bmatrix} \left\{ \begin{matrix} \mathbf{v}_o^b \\ \mathbf{v}_o^d \end{matrix} \right\} \quad (10) $$

In Eq. (10) the subscript *n* stands for the time *t*<sub>*n*</sub> and the subscript *m*, in the summation symbol, stands for the times *t*<sub>*m*</sub>, *m* = 0, 1, 2, ..., (*n* − 1), previous to *t*<sub>*n*</sub>; the superscripts *b* and *d* correspond to the boundary and to the domain (internal points), respectively. In the sub-matrices, the first and second superscripts correspond, respectively, to the positions of source and field points. In the first matrix on the left-hand-side of Eq. (10), the **C** matrix is related to the coefficients *c*(*ξ*) of the boundary nodes whereas the identity matrix, **I** = **H**<sup>dd</sup>, is related to the coefficients *c*(*ξ*) = 1 of the internal points. Alternatively, Eq. (10) can be written in a more compact way as follows

$$ \mathbf{C} \mathbf{u}_n + \sum_{m=0}^{n} \mathbf{H}_{nm} \mathbf{u}_m = \sum_{m=0}^{n} \mathbf{G}_{nm} \mathbf{p}_m + \mathbf{U}_n \mathbf{u}_o + \mathbf{V}_n \mathbf{v}_o \quad (11) $$

After the solution of the system of equations represented by Eq. (11), the velocity $v(\xi, t) = \partial u(\xi, t)/\partial t$ at boundary nodes and internal points is computed by employing Eq. (6); in matrix form, one has a system of equations, similar to that represented by Eq. (10), given below

$$ \begin{bmatrix} \mathbf{C} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \left\{ \begin{matrix} \mathbf{v}_n^b \\ \mathbf{v}_n^d \end{matrix} \right\} = \sum_{m=0}^{n} \left[ \frac{\bar{\mathbf{G}}_{nm}^{bb}}{\bar{\mathbf{G}}_{nm}^{db}} \right] \left\{ \mathbf{p}_m^b \right\} - \sum_{m=0}^{n} \left[ \frac{\bar{\mathbf{H}}_{nm}^{bb}}{\bar{\mathbf{H}}_{nm}^{db}} \mathbf{0} \right] \left\{ \mathbf{u}_m^b \right\} + \left[ \frac{\bar{\mathbf{H}}_{nm}^{db}}{\bar{\mathbf{H}}_{nm}^{bb}} \mathbf{I} \right] \left\{ \mathbf{u}_m^d \right\} $$
$$ \begin{bmatrix} \bar{\mathbf{U}}_n^{bb} & \bar{\mathbf{U}}_n^{bd} \\ \bar{\mathbf{U}}_n^{db} & \bar{\mathbf{U}}_n^{dd} \end{bmatrix} \begin{Bmatrix} \mathbf{u}_o^b \\ \mathbf{u}_o^d \end{Bmatrix} + \begin{bmatrix} \bar{\mathbf{V}}_n^{bb} & \bar{\mathbf{V}}_n^{bd} \\ \bar{\mathbf{V}}_n^{db} & \bar{\mathbf{V}}_n^{dd} \end{bmatrix} \begin{Bmatrix} \mathbf{v}_o^b \\ \mathbf{v}_o^d \end{Bmatrix} \quad (12) $$

In compact form, Eq. (12) can be written as follows

$$ \mathbf{C}\mathbf{v}_n = \sum_{m=0}^{n} \bar{\mathbf{G}}_{nm}\mathbf{p}_m - \sum_{m=0}^{n} \bar{\mathbf{H}}_{nm}\mathbf{u}_m + \bar{\mathbf{U}}_n\mathbf{u}_o + \bar{\mathbf{V}}_n\mathbf{v}_o \quad (13) $$

Since matrix C is a diagonal matrix, Eq. (13) can easily be solved for $v_n$.

Note that the over bared matrices in Eqs. (12), (13) correspond to the matrices that appear in Eqs. (10-11).

## 4. Step-by-step procedure

The solution procedure employed in this article is based mainly in the time translation property, e.g. (Mansur 1983), which can be written as

$$ u^*(X, t; \xi, \tau) = u^*(X, t + t_1; \xi, \tau + t_1) \quad (14) $$

where $t_1$ is any time increment.

Due to the time-translation property, one can write (note that, in the expression below, $k^*(X, t; \xi, \tau)$ represents a generic kernel and $q(X, \tau)$, a boundary variable)

$$ \int_{0}^{\Delta t} k^*(X, t; \xi, \tau)q(X, \tau)d\tau = \int_{\Delta t}^{2\Delta t} k^*(X, t; \xi, \tau)q(X, \tau)d\tau = \dots = \int_{(n-1)\Delta t}^{n\Delta t} k^*(X, t; \xi, \tau)q(X, \tau)d\tau \quad (15) $$

The application of the translation property, as illustrated by Eq. (14), together with Eqs. (1) and (6), enables one to solve the time-dependent problem by considering the results related to $\mathbf{u}_{n-1}$ and $\mathbf{v}_{n-1}$ as “initial conditions” in the computation of the results at time $t_n$. That is, in the computation of $\mathbf{u}_n$ and $\mathbf{v}_n$, the time integration is restricted only to the interval $[0, \Delta t]$. For the approach proposed here, the following versions of Eq. (11) and Eq. (13) can be written

$$ \mathbf{C}\mathbf{u}_n + \mathbf{H}_{nn}\mathbf{u}_n = \mathbf{G}_{nm}\mathbf{p}_n + \mathbf{U}_n\mathbf{u}_{n-1} + \mathbf{V}_n\mathbf{v}_{n-1} \quad (16) $$

$$ \mathbf{C}\mathbf{v}_n = \bar{\mathbf{G}}_{nn}\mathbf{p}_n - \bar{\mathbf{H}}_{nn}\mathbf{u}_n + \bar{\mathbf{U}}_n\mathbf{u}_{n-1} + \bar{\mathbf{V}}_n\mathbf{v}_{n-1} \quad (17) $$

The values of the potentials and of the velocities, treated as “initial conditions”, are stored in the vectors $\mathbf{u}_{n-1}$ and $\mathbf{v}_{n-1}$. Note that proceeding this way, the entire domain needs to be discretized. This feature generates the main limitation of the proposed approach, i.e., the proposed approach is suitable for finite domain applications; infinite domains analyses are feasible only if early time responses are requested.

## 5. Examples

Two examples are presented next. For each example three analyses were carried out, aiming at
<figure><img src="image_2.png" /><figcaption>Fig. 1 Boundary and domain discretization: Mesh 1</figcaption></figure>

verifying the accuracy and convergence of the numerical results to the available analytical solutions (Stephenson 1970). Comparison is also made between the results obtained from the present formulation with those obtained from the standard TD-BEM formulation. From the starting mesh of the first analysis, constituted by 48 boundary elements and 256 internal cells, see Fig. 1, subsequent meshes were obtained by duplicating the number of boundary elements and quadrupling the number of internal cells. In this way, the second analysis was carried out with 96 boundary elements and 1024 cells; the third analysis, with 192 elements and 4096 cells. As usual in time-domain BEM analyses, the dimensionless parameter $\beta = c\Delta t/l$, which relates the wave velocity (c), the smallest element length ($l$), and the time interval $\Delta t$ (assumed constant), is used to estimate the time-step length. It is important to mention that, for the proposed approach, the optimum time-step length is obtained by taking $\beta = 0.5$. One must observe that this is not the value recommended by (Mansur 1983) for the standard TD-BEM formulation, i.e., $\beta = 0.6$.

## 5.1 One-dimensional rod under a Heaviside-type forcing function

This is the classical example of a one-dimensional rod under a Heaviside-type forcing function, applied instantaneously at $t = 0$ and kept constant in time.

Boundary conditions and geometry are depicted in Fig. 2.

The analytical solution to this problem can be found in (Stephenson 1970) and is given by

$$u(x, t) = \frac{P}{E} a \left[ \frac{x}{a} + \frac{8}{\pi^2} \sum_{n=1}^{\infty} \frac{(-1)^n}{(2n-1)^2} \sin\left(\frac{(2n-1)\pi x}{2a}\right) \cos\left(\frac{(2n-1)\pi ct}{2a}\right) \right] \quad (18)$$

<figure><img src="image_3.png" /><figcaption>Fig. 2 One-dimensional rod under a Heaviside-type forcing function: geometry and boundary conditions</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3 One-dimensional rod under a Heaviside-type forcing function: Potential at node A(a, b/2)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 One-dimensional rod under a Heaviside-type forcing function: Velocity at node A(a, b/2)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5 One-dimensional rod under a Heaviside-type forcing function: Flux at node B(0, b/2)</figcaption></figure>

Results related to the potential at the boundary node A(a, b/2) are compared with the analytical solution, Eq. (18), in Fig. 3. The results related to the velocity at boundary node A(a, b/2) and to the flux at boundary node B(0, b/2) are compared with the corresponding analytical solutions (obtained from Eq. (18)) in Figs. 4 and 5, respectively. A numerical damping, which diminishes with the refinement of the mesh, is observed in the results from the three analyses. For a more complete discussion, results obtained from the standard TD-BEM formulation were also included in Figs. 3-5. Relatively to the TD-BEM formulation it is important to mention that the matrices generated by the
<figure><img src="image_7.png" /><figcaption>Fig. 6 One-dimensional rod under a Heaviside-type forcing function: potential at node A(a, b/2) for a larger observation time</figcaption></figure>

convolution integrals in Eqs. (1) and (6) (represented under the summation symbol in Eqs. (11) and (13)), are stored to compute the time-history contributions to the results at the current time, $t_n$; for this reason, the standard TD-BEM analyses presented in this work were carried out by employing mesh 1. An alternative to avoid the storage of a large number of matrices could be the use of a truncation strategy, such as that presented by Carrer and Mansur (2006). Note that a small amount of numerical damping is also observed in the standard TD-BEM results.

In order to verify the stability, the results related to the potential at boundary node A(a, b/2) are compared with the analytical solution, Eq. (18), for a larger observation time, see Fig. (6).

## 5.2 One-dimensional rod under sinusoidal initial conditions

In this example, the same bar of the previous example, now fixed at both ends, is analysed when subjected to initial conditions of the type

In analysis (a)

$$u_o = b_o \sin(\pi x/a) \quad \text{and} \quad v_o = 0$$

The analytical solution to this problem can be found in (Stephenson 1970)

$$u(x, t) = b_o \cos(\pi ct/a) \sin(\pi x/a) \quad (19)$$

In analysis (b) one has

$$v_o = d_o \sin(\pi x/a) \quad \text{and} \quad u_o = 0$$

The analytical solution to this problem can be found in (Stephenson 1970)

$$u(x, t) = d_o \frac{a}{\pi c} \sin(\pi ct/a) \sin(\pi x/a) \quad (20)$$
<figure><img src="image_8.png" /><figcaption>Fig. 7 One-dimensional rod under sinusoidal initial conditions: Potential at point C(a/2, b/2), analysis (a)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8 One-dimensional rod under sinusoidal initial conditions: Velocity at point C(a/2, b/2), analysis (a)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9 One-dimensional rod under sinusoidal initial conditions: Flux at node A(a, b/2), analysis (a)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10 One-dimensional rod under sinusoidal initial conditions: Potential along the rod for different values of time, analysis (a)</figcaption></figure>

The three meshes adopted to perform analyses (a) and (b) are the same of the previous example. By considering, initially, analysis (a), results related to the potential and to velocity at point C of coordinates (a/2, b/2) are depicted, respectively, in Figs. 7 and 8. Results related to the flux at point A(a, b/2) are depicted in Fig. 9. Fig. 10 presents the configurations of the bar for different values of
<figure><img src="image_12.png" /><figcaption>Fig. 11 One-dimensional rod under sinusoidal initial conditions: Potential at point C(a/2, b/2), analysis (b)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12 One-dimensional rod under sinusoidal initial conditions: Velocity at point C(a/2, b/2), analysis (b)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13 One-dimensional rod under sinusoidal initial conditions: Flux at node A(a, b/2), analysis (b)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14 One-dimensional rod under sinusoidal initial conditions: Potential along the rod for different values of time, analysis (b)</figcaption></figure>

time, represented by the parameter ct/a. Results corresponding to analysis (b) are presented in Figs. 11 to 14, in the same sequence they were presented for analysis (a). Similarly to the first example, one can observe that the mesh refinement is followed by a reduction in the numerical
<figure><img src="image_16.png" /><figcaption>Fig. 15 One-dimensional rod under sinusoidal initial conditions: potential at point C(a/2, b/2), analysis (a), for a larger observation time</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16 One-dimensional rod under sinusoidal initial conditions: potential at point C(a/2, b/2), analysis (b), for a larger observation time</figcaption></figure>

damping; consequently, a better agreement between numerical and analytical results is achieved. Standard TD-BEM analyses were carried out with $\beta = 0.62$ and the corresponding results were included in Figs. 7-9 and 11-13.

In this example, the results related to the potential at the boundary node C(a/2, b/2) are compared with the analytical solutions, given by Eqs. (19) and (20), in Figs. 15 and 16, respectively, for a larger observation time, in order to verify the stability.

One can conclude from Figs. 6, 15 and 16 that the results from the present formulation are quite stable; however, the reduction of the numerical damping is a task that must deserve special attention.

## 6. Conclusions

In this work, a modified version of the standard TD-BEM formulation is presented. The main characteristic of the proposed formulation and, perhaps, what could turn it very attractive, is that the convolution integrals are no longer evaluated (as a matter of fact, time integration is restricted to only one interval). Now, the results at a given value of the time, say $t_n$, are computed with the results related to the potential and to its time derivative at the time $t_{n-1}$ assumed as “initial conditions”. Note that such a version of the TD-BEM formulation, in addition to the basic integral equation, called here potential integral equation, also makes use of the integral equation corresponding to the velocity. Both integral equations are applied simultaneously to all boundary nodes and internal points, in order to obtain the results that will be used as “initial conditions” for the next time-step. This is the main disadvantage of the proposed formulation: as the entire domain needs to be discretized in order to compute the contribution of the “initial conditions” to the results at the time $t_n$, the analysis of infinite domain problems becomes almost impossible or, at least, very
expensive in the computational point of view, if late time results are required. It is important to point out that, although the numerical results are in good agreement with the analytical ones, special attention must be devoted to the reduction of the numerical damping observed in the analyses presented here: if this problem is overtaken, the procedure proposed in this work can easily be extended to elastodynamics.

## References

Abreu, A.I., Carrer, J.A.M. and Mansur, W.J. (2003), “Scalar wave propagation in 2D: A BEM formulation based on the operational quadrature method”, *Eng. Anal. Bound. Elem.*, **27**, 101-105.

Abreu, A.I., Mansur, W.J. and Carrer, J.A.M. (2006), “Initial conditions contribution in a BEM formulation based on the convolution quadrature method”, *Int. J. Numer. Meth. Eng.*, **67**, 417-434.

Antes, H., Schanz, M. and Alvermann, S. (2004), “Dynamic analyses of plane frames by integral equations for bars and timoshenko beams”, *J. Sound Vib.*, **276**, 807-836.

Beskos, D.E. (1997), “Boundary elements in dynamic analysis: Part II (1986-1996)”, *Appl. Mecha. Rev.*, **50**, 149-197.

Beskos, D.E. (2003), *Dynamic Analysis of Structures and Structural Systems* (in Boundary Element Advances in Solid Mechanics, D. Beskos and G. Maier, editors) CISM, Udine.

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**, 2169-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. En.*, **18**, 277-285.

Carrer, J.A.M. and Mansur, W.J. (2004), “Alternative time-marching schemes for elastodynamic analysis with the domain boundary element method formulation”, *Comput. Mech.*, **34**, 387-399.

Carrer, J.A.M. and Mansur, W.J. (2006), “Solution of the two-dimensional scalar wave equation by the time-domain boundary element method: Lagrange truncation strategy in time integration”, *Struct. Eng. Mech.*, **23**, 263-278.

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

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

Dubner, H. and Abate, J. (1968), “Numerical inversion of Laplace transforms by relating them to the Finite fourier cosine transform”, *J. Assoc. Comput. Machinery*, **15**, 115-123.

Durbin, F. (1974), “Numerical inversion of Laplace transforms: An efficient improvement to Dubner and Abate’s method”, *Comput. J.*, **17**, 371-376.

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. Meth. Appl. Mech. Eng.*, **179**, 111-123.

Gaul, L. and Wenzel, W. (2002), “A coupled symmetric BE-FE method for acoustic fluid-structure interaction”, *Eng. Analysis with Boundary Elements*; **26**, 629-636.

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

Hatzigeorgiou, G.D. and Beskos, D.E. (2002), “Dynamic elastoplastic analysis of 3-D structures by the domain/boundary element method”, *Comput. Struct.*, **80**, 339-347.

Houbolt, J.C. (1950), “A recurrence matrix solution for the dynamic response of elastic aircraft”, *J. Aeronautical Sci.*, **17**, 540-550.

Kontoni, D.P.N. and Beskos, D.E. (1993), “Transient dynamic elastoplastic analysis by the dual reciprocity BEM”, *Eng. Anal. Bound. Elem.*, **12**, 1-16.

Kutt, H.R. (1975), “The numerical evaluation of principal value integrals by finite-part integration”, *Numer.*
Math., **24**, 205-210.

Lubich, C. (1988a), “Convolution quadrature and discretized operational calculus I”, *Numer. Math.*, **52**, 129-145.

Lubich, C. (1988b), “Convolution quadrature and discretized operational calculus II”, *Numer. Math.*, **52**, 413-425.

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, W.J. and deLima-Silva, W. (1992), “Efficient time truncation in two-dimensional BEM analysis of transient wave propagation problems”, *Earthq. Eng. Struct. Dyn.*, **21**, 51-63.

Mansur, W.J., Abreu, A.I. and Carrer, J.A.M. (2004), “Initial conditions contribution in frequency-domain analysis”, *Comput. Model. Eng. Sci.*, **6**, 31-42.

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), “Application of 3D time domain boundary element formulation to wave propagation in poroelastic solids”, *Eng. Anal. Bound. Elem.*, **25**, 363-376.

Schanz, M. and Antes, H. (1997), “Application of ‘Operational Quadrature Methods’ in time domain boundary element methods”, *Meccanica*, **32**, 179-186.

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.

Souza, L.A., Carrer, J.A.M. and Martins, C.J. (2004), “A fourth order finite difference method applied to elastodynamics: Finite element and boundary element formulations”, *Struct. Eng. Mech.*, **17**, 735-749.

Stephenson, G. (1970), *An Introduction to Partial Differential Equations for Science Students*, Longman.