# Boundary element method formulations for the solution of the scalar wave equation in one-dimensional problems

Jose Antonio Marques Carrer, Vivianne Luiza Costa

► To cite this version:

Jose Antonio Marques Carrer, Vivianne Luiza Costa. Boundary element method formulations for the solution of the scalar wave equation in one-dimensional problems. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2015, 37 (3), pp.959-971. 10.1007/s40430-014-0226-z. hal-01581476

HAL Id: hal-01581476

https://hal.science/hal-01581476

Submitted on 4 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Boundary element method formulations for the solution of the scalar wave equation in one-dimensional problems

J. A. M. Carrer · V. L. Costa

**Abstract** This work is concerned with the development of two boundary element method (BEM) formulations for the solution of one-dimensional scalar wave propagation problems. The first formulation is called TD-BEM, TD meaning time-domain, as it employs a time-dependent fundamental solution. The second formulation is called D-BEM, D meaning domain, and employs the fundamental solution from the static problem. The Houbolt and the Newmark methods are employed for the time-marching in the D-BEM approach. Two examples, constituted of five analyses, are included.

**Keywords** Boundary element method · Scalar wave · Houbolt method · Newmark method

## 1 Introduction

Problems in one dimension have the merit of giving the researcher experience and confidence to face more elaborated analysis. Most of the difficulties that appear in two- and three-dimensional formulations are simpler to be handled in one-dimensional formulations, turning the study attractive and justifying, in this way, the attention devoted to the solution of one-dimensional problems.

In this work, two BEM formulations are presented for the solution of scalar wave propagation problems in one

dimension. Before presenting the BEM formulations, it is important to note that if a finite interval, say $[0, L]$, is the problem domain, then the boundary is constituted only by the nodes at $x = 0$ and at $x = L$. The first formulation is named TD-BEM, TD meaning time-domain, and employs a time-dependent fundamental solution, e.g. [4, 9, 17, 22]. This formulation is very elegant from the mathematical point of view—the concept of causality appears quite naturally—and, as all the concepts used in two- and three-dimensional analysis already appear here, gives the researcher experience with this kind of BEM formulations. Following the recommendation of Mansur [17] linear and constant time variations are assumed, respectively, for the potential and its normal derivative. The domain, or part of it, is discretized only when the problem presents non-homogeneous initial conditions. The domain discretization employs linear cells. Note that, in the absence of initial conditions, only time approximations are required, as the boundary consists only of the two nodes at $x = 0$ and at $x = L$. The second formulation, named D-BEM, where D means domain, employs a static, or steady-state, fundamental solution instead of a time-dependent one, e.g. [5, 12, 21]. As a consequence, a domain integral containing the second-order time derivative of the potential, or the acceleration, appears in the BEM basic integral equation: this domain integral is responsible for the designation D-BEM, as the solution of the problem requires the discretization of the entire domain, diversely from the TD-BEM formulation, where the domain discretization is required only when non-homogeneous initial conditions appear. Linear cells are employed in the domain discretization. The acceleration is approximated by the Houbolt method [13], which has been extensively used, during the last years, not only in D-BEM formulations but also in the so-called DR-BEM formulations. In the DR formulations,

<footnote>J. A. M. Carrer · V. L. Costa
PPGMNE: Programa de Pós-Graduação e m Métodos Numéricos em Engenharia, Universidade Federal do Paraná, Caixa Postal 19011, Curitiba, PR CEP 81531-990, Brazil
e-mail: carrer@ufpr.br</footnote>
DR meaning dual reciprocity, the domain integral is transformed into boundary integrals by suitable interpolation functions, see [1, 14, 15, 25, 26]. The Newmark method can be cited as well. As it is well known, the Newmark method [19] has been widely used in finite element method (FEM) formulations presenting, over the Houbolt method, a better control of the stability and accuracy, according to the values of the parameters $\beta$ and $\gamma$, see [2, 8]. For this reason, it was chosen to be implemented in the D-BEM formulation developed in this work. It is important to mention that although other time-marching schemes were recently proposed, see [5, 7, 23], the search for acceleration approximations is a task that still deserves attention. The presence of non-homogeneous initial conditions is taken into account by following the procedure presented by Carrer et al. [6]. For a more complete discussion concerning dynamic analysis by the BEM, the reader is referred to Beskos [3] and Manolis and Polyzos [16].

It is important to state that the main novelty of the present paper falls in the development of the TD-BEM formulation for one-dimensional problems. Another task that also deserves mention is the discussion concerning the use of different time-marching schemes in the D-BEM formulation. Finally, it is the authors' opinion that the possibility of development of different BEM formulations for time-domain problems brings, among many questions, the following: which formulation is the best for the problem at hand? The authors hope to have thrown some light on the subject.

Finally, two examples are presented and discussed at the end of the article. The first example contains three analyses. The second example contains two analyses with time variable boundary condition.

## 2 BEM formulation

The scalar wave equation for the one-dimensional problem is given by:

$$ \frac{\partial^2 u}{\partial x^2} = \frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} \qquad (1) $$

where the potential function $u(x, t)$ can represent, for instance, the transversal displacement of a string or the longitudinal displacement of a bar, $c$ is the wave propagation velocity, and $t$ is the time. The domain $\Omega$ is restricted to the interval $[0, L]$; consequently, the boundary $\Gamma$ is represented only by the nodes at $x = 0$ and at $x = L$.

The boundary conditions to Eq. (1), for $t \ge 0$ are the essential (or Dirichlet) in the $\Gamma_u$ part of the boundary and

the natural (or Neumann) in the $\Gamma_p$ part of the boundary. They are written as:

$$ \begin{aligned} u(x, t) &= \bar{u}(x, t) \quad \text{in } \Gamma_u \\ p(x, t) &= \frac{du(x, t)}{dn} = \bar{p}(x, t) \quad \text{in } \Gamma_p \end{aligned} \qquad (2) $$

Note that:

1. in the natural boundary condition, $n(x)$ stands for the direction of the unit outward vector normal to the boundary;

2. for general purposes $\Gamma = \Gamma_u \cup \Gamma_p$;

3. $p(x, t) = \frac{du(x, t)}{dn}$ is proportional to the bar support reaction;

The initial conditions, over the domain $\Omega$, are:

$$ \begin{aligned} u(x, 0) &= u_0(x) \\ \dot{u}(x, 0) &= v_0(x) = \left. \frac{\partial u}{\partial t} \right|_{t=0} \end{aligned} \qquad (3) $$

Anticipating, at this point, a discussion concerning the BEM formulations, one must note that in Eq. (1) it is assumed that the material does not present variable density and, consequently, the velocity $c$ is constant. In this way, the BEM formulations, as presented in the subsequent sections, are limited to the solution of problems constituted of only one type of homogeneous material. This limitation, however, can be overcome with the use of the technique of the sub-regions, in which different parts of the domain with different properties can be taken into account. Each different part constitutes a sub-region and the common interfaces between different sub-regions generate additional equations that can easily be incorporated into the basic ones. The additional equations come from the compatibility and equilibrium conditions. The interested reader is referred to Santiago et al. [24] for additional information concerning this matter. Problems presenting materials with continuously varying density (or velocity) can be simulated approximately with the use of sub-regions. In other words, increasing the number of sub-regions a better representation for problems with continuously varying properties is achieved.

### 2.1 The TD-BEM formulation

By adopting as the weighting function the fundamental solution that satisfies Eq. (4) below:

$$ \frac{\partial^2 u^*(x, t, \xi, \tau)}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2 u^*(x, t, \xi, \tau)}{\partial \tau^2} = \delta(x - \xi)\delta(t - \tau) \qquad (4) $$

a residual statement can be applied to Eq. (1), thus generating the basic BEM equation written below as:
$$
\begin{equation}
\begin{aligned}
u(\xi, t) = & \int_0^t p^*(x, t, \xi, \tau) u(x, \tau) d\tau|_{x=L} + \int_0^t p^*(x, t, \xi, \tau) u(x, \tau) d\tau|_{x=0} \\
& - \int_0^t u^*(x, t, \xi, \tau) p(x, \tau) d\tau|_{x=L} - \int_0^t u^*(x, t, \xi, \tau) p(x, \tau) d\tau|_{x=0} \\
& - \frac{1}{c^2} \int_0^L u^*(x, t, \xi, 0) v_0(x) dx + \frac{1}{c^2} \int_0^L \frac{\partial u^*(x, t, \xi, 0)}{\partial \tau} |_{\tau=0} u_0(x) dx
\end{aligned}
\tag{5}
\end{equation}
$$

The fundamental solution, see Graff [10], represents the
effect, at the field point x, of a unit impulse applied at the
source point ξ at τ = t and is given by:

$$
u^*(x, t, \xi, \tau) = -\frac{c}{2} H[c(t - \tau) - r] \quad (6)
$$

where $H[...]$ is the Heaviside function and $r = lx - \xi l$ is
the distance between the field and source points.

The fundamental solution presents the following
properties:

1. causality: $u^*(x,t,\xi,\tau) = 0$ if $c(t-\tau) < lx - \xi l$;

2. time translation: $u^*(x, t, \xi, \tau) = u^*(x, t + t_1, \xi, \tau + t_1)$;

3. reciprocity: $u^*(x,t,\xi,\tau) = u^*(\xi,-\tau,x,-t)$

The function $p^*(x, t, \xi, \tau)$ that appears in Eq. (5), represents the normal derivative of the fundamental solution. Its expression is given by:

$$
p^*(x, t, \xi, \tau) = \frac{du^*(x, t, \xi, \tau)}{dn(x)} = \frac{c}{2} \delta[r - c(t - \tau)] \quad (7)
$$

Although the limits of integration of the domain inte-
grals in Eq. (5) correspond to the entire domain, the inte-
gration is effectively carried out only in the part of the
domain that presents non-homogeneous initial conditions.

After substituting Eq. (7) in Eq. (5), the Dirac delta
properties can be applied, see [11, 18], resulting in:

$$
\begin{equation}
\begin{split}
u(\xi, t) ={}& \frac{1}{2} u(L; t - \frac{L-\xi}{c}) + \frac{1}{2} u(0; t - \frac{\xi}{c}) \\
& - \int_0^t u^*(x, t, \xi, \tau) p(x, \tau) d\tau|_{x=L} - \int_0^t u^*(x, t, \xi, \tau) p(x, \tau) d\tau|_{x=0} \\
& - \frac{1}{c^2} \int_0^L u^*(x, t, \xi, 0) v_0(x) dx + \frac{1}{c^2} \int_0^L \frac{\partial u^*(x, t, \xi, 0)}{\partial \tau} \Big|_{\tau=0} u_0(x) dx
\end{split}
\tag{8}
\end{equation}
$$

The domain integration is analyzed in details in what follows. The domain is approximated by linear cells $\Omega_k = [x_k, x_{k+1}]$ of the same length, say $\Delta L = x_{k+1} - x_k$, with $0 \le x_k$ and $x_{k+1} \le L$. Inside each cell, a linear variation is assumed for the initial conditions and one can write, for $x_k \le x \le x_{k+1}$:

$$
v_0(x) = \Phi_I(x)v_0(x_k) + \Phi_F(x)v_0(x_{k+1}) \quad (9)
$$

and

$$
u_0(x) = \Phi_I(x)u_0(x_k) + \Phi_F(x)u_0(x_{k+1}) \quad (10)
$$

where:

$$
\begin{equation}
\begin{aligned}
\Phi_I(x) &= \frac{x_{k+1} - x}{\Delta L} \\
\Phi_F(x) &= \frac{x - x_k}{\Delta L}
\end{aligned}
\tag{11}
\end{equation}
$$

The first domain integral is denoted as:

$$
I_v(x) = \frac{1}{c^2} \int_0^L u^*(x, t, \xi, 0) v_0(x) dx \quad (12)
$$

Two relative positions of the source point with respect to
the cell are considered: in the first, $\xi \le x_k$; in the second,
$\xi \ge x_{k+1}$. With respect to the first position, the three cases
listed below are analyzed:

1. $\xi + ct < x_k$;

2. $x_k < \xi + ct < x_{k+1}$;

3. $\xi + ct > x_{k+1}$;

For the first case, due to the causality property, $I_v = 0$.
For the second and third cases, the resulting integration can
be represented as:

$$
I_v = \frac{1}{c^2} \int_0^L u^*(x, t, \xi, 0) v_0(x) dx = \frac{1}{c^2} \int_{x_k}^{x_f} u^*(x, t, \xi, 0) v_0(x) dx \quad (13)
$$

where, for the second case, $x_f = ct + \xi$ and, for the third
case, $x_f = x_{k+1}$. After substituting (6) and (10) in Eq. (13),
it can be rewritten as:

$$
I_v = - \frac{1}{2c\Delta L} \left[ \left( x_{k+1}x - \frac{x^2}{2} \right) \left( \frac{x^2}{2} - x_kx \right) \right]_{x_k}^{x_f} \left\{ \frac{v_0(x_k)}{v_0(x_{k+1})} \right\} \quad (14)
$$

when $\xi \ge x_{k+1}$, a similar procedure can be followed in
order to compute the integral, that is, three cases are
identified and listed below:

1 $x_{k+1} < \xi - ct$;

2 $x_k < \xi - ct < x_{k+1}$;

3 $\xi - ct < x_k$;

For the first case, due to the causality property, $I_v = 0$.
For the second and third cases, the resulting integration can
be represented as:
$$I_v = \frac{1}{c^2} \int_0^L u^*(x, t, \xi, 0) v_0(x) dx = \frac{1}{c^2} \int_{x_f}^{x_{k+1}} u^*(x, t, \xi, 0) v_0(x) dx \quad (15)$$

where, for the second case, $x_f = \xi - ct$ and, for the third case, $x_f = x_k$. After substituting (7) and (9) in Eq. (15), it can be rewritten as:

$$I_v = -\frac{1}{2c\Delta L} \left[ \left( x_{k+1}x - \frac{x^2}{2} \right) \left( \frac{x^2}{2} - x_kx \right) \right]_{x_f}^{x_{k+1}} \left\{ \begin{matrix} v_0(x_k) \\ v_0(x_{k+1}) \end{matrix} \right\} \quad (16)$$

The second domain integral is written as:

$$I_u(x) = \frac{1}{c^2} \int_0^L \left. \frac{\partial u^*(x, t, \xi, \tau)}{\partial \tau} \right|_{\tau=0} u_0(x) dx \quad (17)$$

where:

$$\left. \frac{\partial u^*(x, t, \xi, \tau)}{\partial \tau} \right|_{\tau=0} = \frac{c^2}{2} \delta(r - ct) \quad (18)$$

Following the previous reasoning, for $\xi < x_k$, one can consider two situations. In the first, if $\xi + ct < x_k$ or if $\xi + ct > x_{k+1}$, then $I_u = 0$. In the second situation, one has $x_k \le \xi + ct \le x_{k+1}$. Additionally, the distance between field and source points is given by: $r = x - \xi$, allowing Eq. (18) be written as:

$$\left. \frac{\partial u^*(x, t, \xi, \tau)}{\partial \tau} \right|_{\tau=0} = \frac{c^2}{2} \delta[x - (\xi + ct)] \quad (19)$$

Finally, after applying the delta properties:

$$I_u = \frac{1}{2} u(\xi + ct) \quad (20)$$

By taking into account the linear approximation of $u_0(x)$ inside each cell, one finally has:

$$\begin{aligned} I_u &= \frac{1}{2} u(\xi + ct) \\ &= \frac{1}{2} \left[ \left( \frac{(x_{k+1} - \xi) - ct}{\Delta L} \right) \left( \frac{(\xi - x_k) + ct}{\Delta L} \right) \right] \left\{ \begin{matrix} u_0(x_k) \\ u_0(x_{k+1}) \end{matrix} \right\} \end{aligned} \quad (21)$$

For $\xi > x_{k+1}$, $I_u = 0$ if $\xi - ct < x_k$ or if $\xi - ct > x_{k+1}$. If $x_k \le \xi - ct \le x_{k+1}$, then $r = \xi - x$ and Eq. (18) is written as:

$$\left. \frac{\partial u * (x, t, \xi, \tau)}{\partial \tau} \right|_{\tau=0} = \frac{c^2}{2} \delta[x - (\xi - ct)] \quad (22)$$

which results in:

$$I_u = \frac{1}{2} u(\xi - ct) \quad (23)$$

or:

$$\begin{aligned} I_u &= \frac{1}{2} u(\xi - ct) \\ &= \frac{1}{2} \left[ \left( \frac{(x_{k+1} - \xi) + ct}{\Delta L} \right) \left( \frac{(\xi - x_k) - ct}{\Delta L} \right) \right] \left\{ \begin{matrix} u_0(x_k) \\ u_0(x_{k+1}) \end{matrix} \right\} \end{aligned} \quad (24)$$

After carrying out all the integrations, the resulting system of equations can be written as:

$$H^{nm} u^n = G^{nm} p^n + \sum_{m=0}^{n-1} G^{nm} p^m - H^{nm} u^m + U^n u_0 + V^n v_0 \quad (25)$$

In Eq. (25), matrices $H^{nm}$, $G^{nm}$ and $H^{nm}$, $G^{nm}$ have dimension $2 \times 2$, as only boundary variables are associated with them.

After the solution of the system, potential at internal points can be easily determined.

## 2.2 The D-BEM formulation

The fundamental solution employed in the development of the D-BEM formulation is the solution of the equation below:

$$\frac{d^2 u^*(\xi, x)}{dx^2} = \delta(x - \xi) \quad (26)$$

Now, the function $u^*(\xi, x)$ represents the effect, at the field point $x$, of a unit source, $\delta(x - \xi)$, applied at the source point $\xi$. The expression for $u^*(\xi, x)$ can be found, for instance, in [27], and is given by:

$$u^*(\xi, x) = \frac{|x - \xi|}{2} = \frac{r}{2} \quad (27)$$

Note that when $x = \xi$ the fundamental solution is null, diversely from the two- and three-dimensional ones, that present singularities of the type $\ln(r)$ and $1/r$, respectively.

Although the time variable does not appear in the fundamental solution, it can be used in a weighted residuals statement to find an approximate solution to Eq. (1), in this way generating the basic integral equation of the D-BEM formulation. This equation is written as follows:
$$
\begin{equation}
\begin{aligned}
u(\xi, t) ={}& p^*(\xi, x)u(x, t)|_{x=L} + p^*(\xi, x)u(x, t)|_{x=0} \\
& - u^*(\xi, x)p(x, t)|_{x=L} - u^*(\xi, x)p(x, t)|_{x=0} \\
& + \frac{1}{c^2} \int_0^L u^*(\xi, x)\ddot{u}(x, t)dx
\end{aligned}
\tag{28}
\end{equation}
$$

Note that Eq. (28) is valid when $\xi \in \Gamma$ and when $\xi \in \Omega$;
in other words, the coefficient $c(\xi)$ that appears in 2D and
3D BEM formulations is the same, independently of the
position of the source point, and is equal to 1.

The expression for $p^*(\xi, x)$ is given by:

$$
p^*(\xi, x) = \frac{1}{2} \tag{29}
$$

The solution for Eq. (28) requires the entire domain discretization, which is carried out by employing, as in the TD-BEM formulation, linear cells. The acceleration is assumed to vary linearly inside each cell. Once the spatial discretization has been carried out, the resulting matrices can be assembled, thus generating an enlarged system of equations, typical from D-BEM formulations, [6], written below:

$$
\begin{bmatrix} H^{bb} & 0 \\ H^{db} & I \end{bmatrix} \begin{Bmatrix} u_{n+1}^b \\ u_{n+1}^d \end{Bmatrix} = \begin{bmatrix} G^{bb} \\ G^{db} \end{bmatrix} \begin{Bmatrix} p_{n+1}^b \end{Bmatrix} \\
+ \frac{1}{c^2} \begin{bmatrix} M^{bb} & M^{bd} \\ M^{db} & M^{dd} \end{bmatrix} \begin{Bmatrix} \ddot{u}_{n+1}^b \\ \ddot{u}_{n+1}^d \end{Bmatrix} \quad (30)
$$

In Eq. (30), the subscript (n + 1) represents the time
t<sub>n+1</sub> = (n + 1)Δt, n ≥ 0, where Δt is the selected time
interval. The superscripts b and d correspond to the
boundary nodes and to the domain internal points,
respectively. In the sub-matrices, the first superscript cor-
responds to the source point and the second superscript, to
the field point. The boundary unknowns are the potential at
Γ<sub>p</sub> and the flux at Γ<sub>u</sub> (as usual in BEM formulations) and
the potential at the internal points. It is important to men-
tion that the assemblage of such an enlarged system of
equations is necessary because the domain integral relates
boundary values to domain values. Additionally, I is the
identity matrix associated with the internal points.

The sub-matrices $H^{bb}$, $H^{db}$, $G^{bb}$, $G^{db}$ are represented schematically below. It was assumed that the domain discretization employed N cells with the same length, say $\Delta L$; consequently, the internal points are numbered from 3 to $(N+1)$, i.e. $x_k = (k - 2)\Delta L$, $k = 3,4,...,(N+1)$:

$$
H^{bb} = \begin{bmatrix} \frac{1}{2} & -\frac{1}{2} \\ -\frac{1}{2} & \frac{1}{2} \end{bmatrix} \qquad (31)
$$

$$
G^{bb} = \begin{bmatrix} 0 & -\frac{L}{2} \\ -\frac{L}{2} & 0 \end{bmatrix} \qquad (32)
$$

$$
H^{\text{db}} = \begin{bmatrix} -\frac{1}{2} & -\frac{1}{2} \\ -\frac{1}{2} & -\frac{1}{2} \\ -\frac{1}{2} & -\frac{1}{2} \\ \vdots & \vdots \\ -\frac{1}{2} & -\frac{1}{2} \end{bmatrix} \qquad (33)
$$

$$
G^{\text{db}} = \begin{bmatrix} -\frac{\Delta L}{2} & -\frac{\Delta L}{2} \\ -\frac{2\Delta L}{2} & -\frac{L-\Delta L}{2} \\ \vdots & \vdots \\ -\frac{(s-2)\Delta L}{2} & -\frac{L-(s-2)\Delta L}{2} \\ \vdots & \vdots \\ -\frac{(N-2)\Delta L}{2} & -\frac{(N-1)\Delta L}{2} \\ -\frac{(N-1)\Delta L}{2} & -\frac{\Delta L}{2} \end{bmatrix} \quad (34)
$$

In Eq. (34), the letter *s* represents the position of the
source point: ξ = x_s = (s - 2)ΔL, s = 3,4,...,(N + 1).

The same reasoning followed in the computation of the
coefficients of the domain integrals of the TD-BEM for-
mulation is adopted for the computation of the coefficients
in sub-matrices $M^{bb}$, $M^{bd}$, $M^{db}$, $M^{dd}$. Thus, considering a
generic cell $\Omega_k = [x_k, x_{k+1}]$ and a generic source point $\xi_s$,
$s = 1,2,...,(N+1)$, with $\xi_s \le x_k$, a typical domain integral
can be represented as:

$$
\begin{align}
I_{\Omega} &= \frac{1}{c^2} \int_{x_k}^{x_{k+1}} u^*(\xi, x) \ddot{u}(x, t_{n+1}) dx \nonumber \\
&= \frac{1}{2c\Delta L} [M_I \quad M_F] \begin{Bmatrix} \ddot{u}(x_k, t_{n+1}) \\ \ddot{u}(x_{k+1}, t_{n+1}) \end{Bmatrix} \tag{35}
\end{align}
$$

where:

$$
M_I = - \frac{(x_{k+1} - \xi)^3 - (x_k - \xi)^3}{3} + \frac{(x_{k+1} - \xi)}{2} \left[ (x_{k+1} - \xi)^2 - (x_k - \xi)^2 \right] \quad (36)
$$

and

$$
M_F = \frac{(x_{k+1} - \xi)^3 - (x_k - \xi)^2}{3} + \frac{(3 - x_k)}{2} \left[ (x_{k+1} - \xi)^2 - (x_k - \xi)^2 \right] \quad (37)
$$

For the situation in which $\xi_s \ge x_{k+1}$, the same result is
obtained, but multiplied by $(-1)$.
### 2.2.1 Time-marching scheme

The Houbolt method, [13], is obtained by cubic Lagrange interpolation of $u = u(t)$ from time $t_{n-2} = (n-2)\Delta t$ to time $t_{n+1} = (n+1)\Delta t$. Exact time differentiation gives the following approximations to velocity and acceleration:

$$ \dot{u}_{n+1} = \frac{1}{6\Delta t} [11u_{n+1} - 18u_n + 9u_{n-1} - 2u_{n-2}] \quad (38) $$

$$ \ddot{u}_{n+1} = \frac{1}{\Delta t^2} [2u_{n+1} - 5u_n + 4u_{n-1} - u_{n-2}] \quad (39) $$

After substituting (39) in (30), time-march can start, according to the equation below:

$$ \begin{bmatrix} (c\Delta t)^2 H^{bb} - 2M^{bb} & -2M^{bd} \\ (c\Delta t)^2 H^{db} - 2M^{db} & (c\Delta t)^2 I - 2M^{dd} \end{bmatrix} \begin{Bmatrix} u_{n+1}^b \\ u_{n+1}^d \end{Bmatrix} = \begin{bmatrix} (c\Delta t)^2 G^{bb} \\ (c\Delta t)^2 G^{db} \end{bmatrix} \{p_{n+1}^b\} + \frac{1}{c^2} \begin{bmatrix} M^{bb} & M^{bd} \\ M^{db} & M^{dd} \end{bmatrix} \times \begin{Bmatrix} -5u_n^b + 4u_{n-1}^b - u_{n-2}^b \\ -5u_n^d + 4u_{n-1}^d - u_{n-2}^d \end{Bmatrix} \quad (40) $$

As mentioned before, the Houbolt method has been widely used with both D-BEM and DR-BEM formulations, see Carrer and Mansur [6] and [1].

For the Newmark method, [19], the approximations are:

$$ \dot{u}_{n+1} = \frac{\gamma}{\beta\Delta t} [u_{n+1} - u_n] + \frac{\beta - \gamma}{\beta} \dot{u}_n - \frac{\gamma - 2\beta}{2\beta} \Delta t \ddot{u}_n \quad (41) $$

$$ \ddot{u}_{n+1} = \frac{1}{\beta\Delta t^2} [u_{n+1} - u_n] - \frac{1}{\beta\Delta t}\dot{u}_n - \frac{1-2\beta}{2\beta}\ddot{u}_n \quad (42) $$

The stability and accuracy of the Newmark method depend on the correct choice of the parameters $\beta$ and $\gamma$. Some remarks concerning the Newmark method are useful: (1) the method is unstable for $\gamma < 1/2$; (2) unconditional stability occurs when $2\beta \ge \gamma \ge 1/2$; (3) taking $\gamma > 1/2$ introduces artificial damping but reduces the accuracy of the method to first order. If one takes $\beta = 1/6$ and $\gamma = 1/2$, expressions (41) and (42) correspond to the linear acceleration method. For $\beta = 1/4$ and $\gamma = 1/2$, expressions (41) and (42) correspond to the average acceleration method. These observations can be found in FEM text books, see [2, 8].

After substituting (42) in (30), one has:

$$ \begin{bmatrix} \beta(c\Delta t)^2 H^{bb} - M^{bb} & -M^{bd} \\ \beta(c\Delta t)^2 H^{db} - M^{db} & \beta(c\Delta t)^2 I - M^{dd} \end{bmatrix} \begin{Bmatrix} u_{n+1}^b \\ u_{n+1}^d \end{Bmatrix} = \begin{bmatrix} \beta(c\Delta t)^2 G^{bb} \\ \beta(c\Delta t)^2 G^{db} \end{bmatrix} \{p_{n+1}^b\} - \frac{1}{c^2} \begin{bmatrix} M^{bb} & M^{bd} \\ M^{db} & M^{dd} \end{bmatrix} \times \begin{Bmatrix} u_n^b + \Delta t \dot{u}_n^b + \frac{1-2\beta}{2} \Delta t^2 \ddot{u}_n^b \\ u_n^d + \Delta t \dot{u}_n^d + \frac{1-2\beta}{2} \Delta t^2 \ddot{u}_n^d \end{Bmatrix} \quad (43) $$

Equations (40) and (43) are written in a concise manner as follows:

$$ \bar{H}u_{n+1} = \bar{G}p_{n+1} + g_n \quad (44) $$

in which the contributions of the previous instants of time are stored in vector $g_n$.

As usual in time domain BEM formulations, a dimensionless variable, say $\beta_{\Delta t}$ (Courant limit), is adopted in order to provide a measure of the time-step $\Delta t$, [17]. This variable is defined according to:

$$ \beta_{\Delta t} = \frac{c\Delta t}{\Delta L} \quad (45) $$

For the one-dimensional problem discussed here as mentioned before, $\Delta L$ is the length of the smallest cell used in the domain discretization.

## 3 Numerical results

Some of the problems presented in this section present non-homogeneous initial conditions. The computation of the initial conditions contributions, in the case of the D-BEM formulation, is carried out by following the approach presented by [6]. The procedure is summarized below:

### 1. Houbolt method:

when $n = 0$ in expression (40), one is faced with the values $u_{-2}$ and $u_{-1}$. These values are related with the initial conditions as follows:

$$ u_{-2} = u_0 - 2\Delta t \dot{u}_0 \quad (46) $$

$$ u_{-1} = u_0 - \Delta t \dot{u}_0 \quad (47) $$

### 2. Newmark method:

Equation (43), when $n = 0$, is written as:

$$ \begin{bmatrix} (c\Delta t)^2 H^{bb} - 2M^{bb} & -2M^{bd} \\ (c\Delta t)^2 H^{db} - 2M^{db} & (c\Delta t)^2 I - 2M^{dd} \end{bmatrix} \begin{Bmatrix} u_1^b \\ u_1^d \end{Bmatrix} = \begin{bmatrix} (c\Delta t)^2 G^{bb} \\ (c\Delta t)^2 G^{db} \end{bmatrix} \{p_1^b\} - \begin{bmatrix} 2M^{bb} & 2M^{bd} \\ 2M^{db} & 2M^{dd} \end{bmatrix} \times \begin{Bmatrix} u_0^b + \Delta t \dot{u}_0^b \\ u_0^d + \Delta t \dot{u}_0^d \end{Bmatrix} \quad (48) $$

With the purpose of obtaining accurate results for both the longitudinal displacement and the support reaction of the bar analyzed, several attempts have been made aiming at the determination of the best values of the parameter $\beta_{\Delta t}$ for each method. It was observed that the support reaction jumps in time are best characterized
by increasing the values of $β_{Δt}$; in other words, the oscillations that appear in the results concerned with the support reaction, in the vicinity of points of discontinuity, are attenuated by increasing the values of $β_{Δt}$; the counterpart of this attenuation, however, can be an excessive damping of the results related to the displacement. For the Houbolt method, the best choice falls on $β_{Δt} = 0.32$; for the Newmark method, the best choice is $β_{Δt} = 0.60$. The use of the Newmark method brings another concern: the choice of adequate values of the parameters $β$ and $γ$. The conclusion at which the authors arrived is that accurate results can be achieved for $u(x,t)$ and $p(x,t)$ when $β = 0.3$ and $γ = 0.56$ (note that the damping introduced by taking $γ = 0.56$ is benefic for the support reaction results, see Figs. 4, 8 and 12. In accordance with which was mentioned previously, the adoption of smaller values for the parameter $γ$ is followed by an improvement in the oscillations around the points of discontinuity in the support reaction response, whereas the use of larger values introduce an excessive amount of damping in the displacements. As a final remark, although the aim here is to obtain results for the displacement and for the support reaction with the same level of accuracy, the choice of the values of $β_{Δt}$, $β$ and $γ$ is conditioned, mainly, by the accuracy of the support reaction results.

The D-BEM analysis employed 32 cells of the same length in the domain discretization.

The results produced by the TD-BEM formulation present a remarkable accuracy when compared with the corresponding analytical solutions and, for this reason, were always taken as the basic ones for comparison purposes. It is important to mention that in this formulation there is no lower limit in the choice of the time-step length, diversely from the two-dimensional case, see [17]. However, in order to provide a good description of the results, the same time-step adopted in the D-BEM formulation with the Houbolt method was selected.

By adopting a bar fixed at the left end, at $x = 0$, and free at the right end, at $x = L$, the results are always presented in the following sequence: (1) displacement at boundary node $x = L$ and (2) support reaction at boundary node $x = 0$ from TD-BEM and D-BEM (Houbolt method) analyses; (3) displacement at boundary node $x = L$ and (4) support reaction at boundary node $x = 0$ from TD-BEM and D-BEM (Newmark method) analyses. The length is $L = 2$ m and the material is the steel; consequently, $c = 5,000$ m/s.

In spite of selecting a specific material and length for the bar, in all figures the results are presented in dimensionless form.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Bar under a heaviside-type forcing function: TD-BEM and D-BEM (Houbolt) analysis, displacement at $x = L$</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2 Bar under a heaviside-type forcing function: TD-BEM and D-BEM (Houbolt) analysis, support reaction at $x = 0$</figcaption></figure>

### 3.1 Bar under a Heaviside-type forcing function and under initial conditions prescribed over the entire domain

In this example, the selected bar is subjected to a Heaviside-type forcing function applied instantaneously at $t = 0$ and kept constant in time, that is, $\bar{p}(L, t) = PH(t - 0)$. The initial conditions:
<figure><img src="image_4.png" /><figcaption>Fig. 3 Bar under a heaviside-type forcing function: TD-BEM and D-BEM (Newmark) analysis, displacement at $x = L$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5 Bar under initial displacement field over the entire domain: TD-BEM and D-BEM (Houbolt) analysis, displacement at $x = L$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4 Bar under a heaviside-type forcing function: TD-BEM and D-BEM (Newmark) analysis, support reaction at $x = 0$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6 Bar under initial displacement field over the entire domain: TD-BEM and D-BEM (Houbolt) analysis, support reaction at $x = 0$</figcaption></figure>

$$u_0(x) = Ux, \quad 0 \le x \le L \qquad (49)$$

$$\dot{u}_0(x) = V, \quad 0 \le x \le L \qquad (50)$$

are also taken into account.

The analytical solution is given by:

$$u(x,t) = P x + \frac{8L}{\pi^2} \sum_{n=1}^{\infty} \mu_n \sin \left[ \frac{(2n-1)\pi x}{2L} \right] \qquad (51)$$

where
<figure><img src="image_8.png" /><figcaption>Fig. 7 Bar under initial displacement field over the entire domain: TD-BEM and D-BEM (Newmark) analysis, displacement at $x = L$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8 Bar under initial displacement field over the entire domain: TD-BEM and D-BEM (Newmark) analysis, support reaction at $x = 0$</figcaption></figure>

$$ \mu_n = (P - U) \frac{(-1)^n}{(2n-1)^2} \cos \left[ \frac{(2n-1)\pi ct}{2L} \right] + \frac{V}{c} \sin \left[ \frac{(2n-1)\pi ct}{2L} \right] \quad (52) $$

This example can be divided in three ones, according to the values of the constants $P$, $U$ and $V$. In the first analysis, $U = V = 0$ and the results are presented in Figs. 1, 2, 3, 4.

<figure><img src="image_10.png" /><figcaption>Fig. 9 Bar under initial velocity field over the entire domain: TD-BEM and D-BEM (Houbolt) analysis, displacement at $x = L$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10 Bar under initial velocity field over the entire domain: TD-BEM and D-BEM (Houbolt) analysis, support reaction at $x = 0$</figcaption></figure>

The results in Figs. 5, 6, 7, 8 correspond to the second analysis, in which $P = V = 0$.

Finally, in the third analysis one has $P = U = 0$ with the results presented in Figs. 9, 10, 11, 12.

Note that in the last two analyses the TD-BEM required domain discretization, which was carried out with only one cell.
<figure><img src="image_12.png" /><figcaption>Fig. 11 Bar under initial velocity field over the entire domain: TD-BEM and D-BEM (Newmark) analysis, displacement at $x = L$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12 Bar under initial velocity field over the entire domain: TD-BEM and D-BEM (Newmark) analysis, support reaction at $x = 0$</figcaption></figure>

Equations (46), (47) and (48) are employed for the computation of the initial conditions' contributions in the second and third analyses, leading to reliable results and, consequently, demonstrating their usefulness. In all analyses, results related to the displacement, computed by employing the Houbolt method, present less numerical damping than the corresponding ones computed with the Newmark method. Results from the TD-BEM formulation

<figure><img src="image_14.png" /><figcaption>Fig. 13 Bar with sinusoidal boundary condition $\bar{p}(L,t) = P \sin(\psi t)$, $\psi = 2,000$ Hz: TD-BEM and D-BEM (Houbolt) analysis, displacement at $x = L$</figcaption></figure>

do not present numerical damping, even for values of ct/L greater than those presented in the figures of this section. Results related to the support reaction present oscillations around the points of discontinuity for the Houbolt method. For the Newmark method, the introduction of numerical damping was beneficial, at least at the beginning of the analysis, as a noticeable damping will deteriorate the results for greater values of ct/L. The points of discontinuity, or the jumps, are very well represented by the TD-BEM results. As a matter of fact, the TD-BEM results present a delay equal to $\Delta t$ relatively to the analytical solution (that is another reason for choosing small values of $\Delta t$) but do not present the Gibbs phenomenon.

## 3.2 Bar with a time variable boundary condition

The bar, now, is subjected to a time variable sinusoidal boundary condition, written as: $\bar{p}(L,t) = P \sin(\psi t)$.

By assuming homogeneous initial conditions, the analytical solution to this problem is given by:

$$ u(x,t) = U \left\{ \frac{x}{L} \sin(\psi t) + \sum_{n=1}^{\infty} \frac{2\psi}{\phi_n^2 - \psi^2} \frac{(-1)^n}{n\pi} \times [\phi_n \sin(\phi_n t) - \psi \sin(\psi t)] \sin\left(\frac{n\pi x}{L}\right) \right\} \quad (53) $$

where

$$ \phi_n = \frac{nc\pi}{2L} \quad (54) $$
<figure><img src="image_15.png" /><figcaption>**Fig. 14** Bar with sinusoidal boundary condition $\bar{p}(L,t) = P \sin(\psi t)$, $\psi = 2,000$ Hz: TD-BEM and D-BEM (Houbolt) analysis, support reaction at $x = 0$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15** Bar with sinusoidal boundary condition $\bar{p}(L,t) = P \sin(\psi t)$, $\psi = 2,000$ Hz: TD-BEM and D-BEM (Newmark) analysis, displacement at $x = L$</figcaption></figure>

Two analysis were carried out. In the first, $\psi = 2,000$ Hz. The results related to the Houbolt method are presented in Figs. 13 and 14. Newmark results are presented in Figs. 15 and 16.

In the second analysis, $\psi = 3,927$ Hz was considered. Note that this value coincides with the first natural

<figure><img src="image_17.png" /><figcaption>**Fig. 16** Bar with sinusoidal boundary condition $\bar{p}(L,t) = P \sin(\psi t)$, $\psi = 2,000$ Hz: TD-BEM and D-BEM (Newmark) analysis, support reaction at $x = 0$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 17** Bar with sinusoidal boundary condition $\bar{p}(L,t) = P \sin(\psi t)$, $\psi = 3,927$ Hz: TD-BEM and D-BEM (Houbolt) analysis, displacement at $x = L$</figcaption></figure>

frequency, see Rao [20], and, then, the phenomenon of resonance occurs. The results are presented in Figs. 17 and 18 for the Houbolt method and in Figs. 19 and 20 for the Newmark method. A good agreement is observed in all figures between the D-BEM and the TD-BEM results.
<figure><img src="image_19.png" /><figcaption>**Fig. 18** Bar with sinusoidal boundary condition $\bar{p}(L,t) = P \sin(\psi t)$, $\psi = 3,927$ Hz: TD-BEM and D-BEM (Houbolt) analysis, support reaction at $x = 0$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 19** Bar with sinusoidal boundary condition $\bar{p}(L,t) = P \sin(\psi t)$, $\psi = 3,927$ Hz: TD-BEM and D-BEM (Newmark) analysis, displacement at $x = L$</figcaption></figure>

## 4 Conclusions

The D-BEM formulation presented in this work, for the solution of scalar wave propagation problems, proved to be very efficient, leading to accurate examples. Naturally, experience still plays an important role in the selection of the time interval, in other words, in the selection of the

<figure><img src="image_21.png" /><figcaption>**Fig. 20** Bar with sinusoidal boundary condition $\bar{p}(L,t) = P \sin(\psi t)$, $\psi = 3,927$ Hz: TD-BEM and D-BEM (Newmark) analysis, support reaction at $x = 0$</figcaption></figure>

parameter $\beta_{\Delta t}$, for both the Houbolt and the Newmark methods, and in the choice of the parameters $\beta$ and $\gamma$ in the Newmark method. In this way, this work intends to provide some useful values of these parameters to those interested in this attractive area of research. The possibility of solving problems with non-homogeneous initial conditions is another task that deserves attention and the expressions presented with this aim proved to be correct (as a matter of fact, such expressions were already successfully employed in the solution of two-dimensional problems). The TD-BEM formulation, on the other hand, produced reliable results and fully justifies what was written at the beginning of this work, that is, it gives experience and confidence to face more elaborate problems.

## References

1. Agnantiaris JP, Polyzos D, Beskos DE (1998) Three-dimensional structural vibration analysis by the dual reciprocity BEM. Comput Mech 21:372-381

2. Bathe KJ (1996) Finite element procedures. Prentice Hall Inc, New Jersey

3. Beskos DE (1997) Boundary elements in dynamic analysis: part II (1986-1996). Appl Mech Rev 50:149-197

4. Carrer JAM, Mansur WJ (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 Methods Eng 18:277-285

5. Carrer JAM, Mansur WJ (2004) Alternative time-marching schemes for elastodynamic analysis with the domain boundary element method formulation. Comput Mech 34:387-399
6. Carrer JAM, Mansur WJ, Vanzuit RJ (2009) Scalar wave equation by the boundary element method: A D-BEM approach with non-homogeneous initial conditions. Comput Mech 44:31–44

7. Chien CC, Chen YH, Chuang CC (2003) Dual reciprocity BEM analysis of 2D transient elastodynamic problems by time-discontinuous galerkin FEM. Eng Anal Bound Elem 27:611–624

8. Cook RD, Malkus DS, Plesha ME (1989) Concepts and applications of finite element analysis. Wiley, New York

9. Dominguez J (1993) Boundary elements in dynamics. Computational Mechanics, Southampton

10. Graff KF (1991) Wave motion in elastic solids. Dover Publications, Inc., New York

11. Greenberg MD (1971) Application of green's functions in engineering and sciences. Prentice-Hall Inc, Eaglewood Cliffs

12. Hatzigeorgiou GD, Beskos DE (2002) Dynamic elastoplastic analysis of 3D structures by the domain/boundary element method. Comput Struct 80:339–347

13. Houbolt JC (1950) A recurrence matrix solution for the dynamic response of elastic aircraft. J Aeronaut Sci 17:540–550

14. Ochiai Y, Sladek V, Sladek J (2005) Transient heat conduction analysis by triple-reciprocity boundary element method. Eng Anal Bound Elem 30:194–204

15. Loeffler CF, Neves FP, Oliveira PC (2010) Dual reciprocity boundary element formulation applied to non-linear Darcian diffusive-advective problems. In: Proceedings of the international conference on boundary elements XI, Berlin

16. Manolis G, Polyzos D (2009) Recent advances in boundary element methods: a volume to honor professor dimitri Beskos. Springer, Berlin

17. Mansur WJ (1983). A time-stepping technique to solve wave propagation problems using the boundary element method. Southampton University, Ph. D Tesis

18. Mei CC (1995) Mathematical analysis in engineering. Cambridge University Press, Cambridge

19. Newmark NM (1959) A method of computation for structural Dynamics. ASCE J Eng Mech Div 85:67–94

20. Rao SS (1995) Mechanical vibrations, 3rd edn. Addison-Wesley Publishing Company, Inc., USA

21. Soares D Jr, Telles JCF, Mansur WJ (2006) Boundary elements with equilibrium satisfaction—a consistent formulation for dynamic problems considering non-linear effects. Int J Numer Meth Eng 65:701–713

22. Soares D Jr, Vinagre MP (2008) Numerical computational of electromagnetic fields by the time-domain boundary element method and the complex variable method. Comput Model Eng Sci 25:1–8

23. Souza LA, Carrer JAM, Martins CJ (2004) A fourth order finite difference method applied to elastodynamics: finite element and boundary element formulations. Struct Eng Mech 17:735–749

24. Santiago JAF, Telles JCF, Valentin VA (1999) An optimized block matrix manipulation for boundary elements with subregions. Adv Eng Softw 30:701–713

25. Tanaka M, Kurokawa K, Matsumoto T (2008) A time-stepping DRBEM for transient heatconduction in anisotropic solids. Eng Anal Bound Elem 32:1046–1053

26. Vera-Tudela CAR, Telles JCF (2005) A numerical green's function and dual reciprocity BEM method to solve elastodynamic crack problems. Eng Anal Bound Elem 29:204–209

27. Vladimirov VS (1979) Generalized functions in mathematical physics. Nauka, Moscow