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

HAL
open science

# Scalar wave equation by the boundary element method: a D-BEM approach with non-homogeneous initial conditions

J. A. M. Carrer, W. J. Mansur, R. J. Vanzuit

► To cite this version:

J. A. M. Carrer, W. J. Mansur, R. J. Vanzuit. Scalar wave equation by the boundary element method: a D-BEM approach with non-homogeneous initial conditions. Computational Mechanics, 2009, 44 (1), pp.31. 10.1007/s00466-008-0353-4 . hal-01352250

HAL Id: hal-01352250

https://hal.science/hal-01352250

Submitted on 6 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Scalar wave equation by the boundary element method: a D-BEM approach with non-homogeneous initial conditions

J. A. M. Carrer

W. J. Mansur

PPGMNE: Programa de Pós-Graduação em Métodos Numéricos em Engenharia, Universidade Federal do Paraná, Caixa Postal 19011, Curitiba, PR CEP 81531-990, Brazil
carrer@ufpr.br, carrer@mat.ufpr.br

R. J. Vanzuit

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

**Abstract** This work is concerned with the development of a D-BEM approach to the solution of 2D scalar wave propagation problems. The time-marching process can be accomplished with the use of the Houbolt method, as usual, or with the use of the Newmark method. Special attention was devoted to the development of a procedure that allows for the computation of the initial conditions contributions. In order to verify the applicability of the Newmark method and also the correctness of the expressions concerned with the computation of the initial conditions contributions, four examples are presented and the D-BEM results are compared with the corresponding analytical solutions.

**Keywords** Boundary element method ·

Scalar wave equation · Initial conditions · D-BEM

## 1 Introduction

The development of BEM formulations for solving time-dependent problems is a very attractive area of research, as demonstrated by the great number of formulations that appeared during the last years, enriching the BEM literature concerning this matter. With the purpose of attesting the contribution of this work, a brief discussion concerning the BEM application to time-dependent problems is carried out

in what follows. For a very complete discussion concerning dynamic analysis by the BEM, the reader is referred to Beskos [1,2].

Bearing in mind that the integral equations can be obtained by means of weighted residuals statement, the solution of time-dependent problems can be accomplished with the use of time-dependent fundamental solutions: in this case, BEM formulations are denominated TD-BEM (TD meaning time-domain), e.g. Mansur [3], Dominguez [4], Carrer and Mansur [5]. TD-BEM formulations are very elegant, from the mathematical point of view, and the fulfilment of the radiation condition makes them suitable for infinite domain analyses. Also, good representations of causality and of time response jumps lead to very accurate results. These favorable 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, based on truncation in the time integration, were presented by Demirel and Wang [6], Mansur and de Lima-Silva [7], Soares and Mansur [8], Carrer and Mansur [9].

Alternatively, one can use static fundamental solutions instead of time-dependent fundamental solutions. In this case, the BEM basic integral equation presents a domain integral with the kernel constituted by the fundamental solution multiplied by the second order time derivative of the potential. In the so-called D-BEM formulations (D meaning domain) this domain integral is kept in the BEM equations, e.g. Carrer and Mansur [10], Hatzigeorgiou and Beskos [11]. On the other hand, by means of suitable interpolation functions the domain integral can be transformed into boundary integrals, generating the so-called DR-BEM formulations (DR meaning double reciprocity), e.g. Kontoni and Beskos [12], Partridge et al. [13], Agnantiaris et al. [14,15]. Although
the treatment given to the domain integral is completely different, one common feature of the D-BEM and DR-BEM formulations is that the time variable does not appear explicitly in the integral equations. As a consequence, an adequate approximation to the acceleration is of fundamental importance if reliable results are to be found; this requirement was fulfilled by the Houbolt method, Houbolt [16], successfully employed with the BEM. Although alternative time-marching schemes were recently proposed by Carrer and Mansur [10], Souza et al. [17], Chien et al. [18], the search for other approximations is a task that still deserves attention. Among the various schemes presented in the Finite Element Method (FEM) literature, e.g. Bathe [19], Cook et al. [20], the Newmark family of methods can be cited. As it is well known, the Newmark method, Newmark [21] has been widely used in 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 Bathe [19], Cook et al. [20]. For this reason, it was chosen to be implemented in the D-BEM formulation presented in this work. This is the first matter to be discussed here. The second task, that deserves special attention, is concerned with the solution of problems with non-homogeneous initial conditions: a general approach is presented to solve this kind of problems by employing both the Houbolt and the Newmark time-marching schemes.

Four examples are presented and discussed at the end of the article, with the aim of validating the proposed formulation. Additionally, it is expected that the conclusions from this work remain valid for the DR-BEM formulation.

## 2 D-BEM formulation

The scalar wave equation for 2D problems over a domain $\Omega$ limited by a boundary $\Gamma$ is given by:

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

where $u(X, t)$, that is generically referred to as potential function, can represent the transversal displacements of a membrane, $c$ is the wave propagation velocity, $t$ is the time and $X$ represents the point of coordinates $(x, y)$.

The boundary conditions to Eq. 1, for $t \ge 0$, are given by:

$$ \begin{aligned} u(X, t) &= \bar{u}(X, t) \quad \text{on } \Gamma_u \\ (\text{essential or Dirichlet boundary condition}) \\ p(X, t) &= \frac{du(X, t)}{dn(X)} = \bar{p}(X, t) \quad \text{on } \Gamma_p \\ (\text{natural or Neumann boundary condition}) \end{aligned} \quad (2) $$

Note that:

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

(ii) for general purposes $\Gamma = \Gamma_u \cup \Gamma_p$

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

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

In the BEM literature, the solution of the Poisson equation:

$$ \frac{\partial^2 u^*(\xi, X)}{\partial x^2} + \frac{\partial^2 u^*(\xi, X)}{\partial y^2} = -\delta(X - \xi), \quad (4) $$

where $\delta(X - \xi)$ is the Dirac delta function, is the so-called fundamental solution and represents the effect, in the field point $X$, of a Dirac function applied at the source point $\xi$. The fundamental solution is:

$$ u^*(\xi, X) = \frac{1}{2\pi} \ln\left(\frac{1}{r}\right), \quad (5) $$

where $r$ is the distance between the field and the source points.

Note that although the time variable does not appear in the fundamental solution $u^*(\xi, X)$, it can be used in a weighted residuals statement for finding an approximate solution to Eq. 1, in this way generating the basic integral equation of the D-BEM formulation. This equation can be written as follows:

$$ \begin{aligned} c(\xi)u(\xi, t) = & \int_{\Gamma} u^*(\xi, X)p(X, t) \, d\Gamma(X) \\ & - \int_{\Gamma} p^*(\xi, X)u(X, t) \, d\Gamma(X) \\ & - \frac{1}{c^2} \int_{\Omega} u^*(\xi, X)\ddot{u}(X, t) \, d\Omega(X). \end{aligned} \quad (6) $$

The coefficient $c(\xi)$ is computed according to:

$$ c(\xi) = \frac{\alpha}{2\pi}, \quad (7) $$

where $\alpha$ is the angle depicted in Fig. 1. The function $p^*(\xi; X)$ is the normal derivative of the fundamental solution:

$$ p^*(\xi, X) = \frac{du^*(\xi, X)}{dn} = \frac{du^*(\xi, X)}{dr} \frac{dr}{dn}. \quad (8) $$

In order to solve Eq. 6, boundary and domain discretizations must be carried out and an approximation to the acceleration must be adopted. In the present work, the discretization of the boundary is accomplished by linear boundary elements; the discretization of the domain, by triangular linear cells. The reader is referred to Mansur [3], and Carrer
<figure><img src="image_3.png" /><figcaption>Fig. 1 Internal angle for the computation of c(ξ)</figcaption></figure>

and Mansur [22] for further details concerning this matter. Once the spatial discretization has been carried out, the resulting matrices can be assembled, thus generating an enlarged system of equations, written below:

$$ 
\begin{bmatrix} \mathbf{H}^{bb} & \mathbf{0} \\ \mathbf{H}^{db} & \mathbf{1} \end{bmatrix} \left\{ \begin{matrix} \mathbf{u}_{n+1}^b \\ \mathbf{u}_{n+1}^d \end{matrix} \right\} = \begin{bmatrix} \mathbf{G}^{bb} \\ \mathbf{G}^{db} \end{bmatrix} \left\{ \mathbf{p}_{n+1}^b \right\} - \frac{1}{c^2} \begin{bmatrix} \mathbf{M}^{bb} & \mathbf{M}^{bd} \\ \mathbf{M}^{db} & \mathbf{M}^{dd} \end{bmatrix} \left\{ \begin{matrix} \ddot{\mathbf{u}}_{n+1}^b \\ \ddot{\mathbf{u}}_{n+1}^d \end{matrix} \right\}. \quad (9) $$

In Eq. 9, in order to simplify the notation, the subscript $(n+1)$ represents the time $t_{n+1} = (n+1)\Delta t$, where $\Delta 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 corresponds to the position of the source point and the second superscript, to the position of the field point. The identity matrix is related to the coefficients $c(\xi) = 1$ of the internal points.

From Eq. 9, the boundary unknowns are the potential at $\Gamma_p$ and the flux at $\Gamma_u$ (as usual in BEM formulations) and the potential at the internal points. It is important to mention that the assemblage of such an enlarged system of equations is necessary because the domain integral relates boundary values to domain values. This remark is confirmed by matrix $\mathbf{M}^{bd}$ in Eq. 9.

## 2.1 Time-marching scheme

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

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

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

After substituting (11) in (9), the time-marching scheme can start:

$$ 
\begin{aligned}
& \begin{bmatrix} (c\Delta t)^2 \mathbf{H}^{bb} + 2\mathbf{M}^{bb} & 2\mathbf{M}^{bd} \\ (c\Delta t)^2 \mathbf{H}^{db} + 2\mathbf{M}^{db} & (c\Delta t)^2 \mathbf{I} + 2\mathbf{M}^{dd} \end{bmatrix} \begin{bmatrix} \mathbf{u}_{n+1}^b \\ \mathbf{u}_{n+1}^d \end{bmatrix} \\
& = \begin{bmatrix} (c\Delta t)^2 \mathbf{G}^{bb} \\ (c\Delta t)^2 \mathbf{G}^{db} \end{bmatrix} \begin{bmatrix} \mathbf{p}_{n+1}^b \\ \mathbf{p}_{n+1}^d \end{bmatrix} \\
& \quad - \begin{bmatrix} \mathbf{M}^{bb} & \mathbf{M}^{bd} \\ \mathbf{M}^{db} & \mathbf{M}^{dd} \end{bmatrix} \begin{bmatrix} -5\mathbf{u}_n^b + 4\mathbf{u}_{n-1}^b - \mathbf{u}_{n-2}^b \\ -5\mathbf{u}_n^d + 4\mathbf{u}_{n-1}^d - \mathbf{u}_{n-2}^d \end{bmatrix}
\end{aligned} \quad (12) $$

As mentioned before, the Houbolt method has been widely used with both D-BEM and DR-BEM formulations, e.g. Carrer and Mansur [10], Hatzigeorgiou and Beskos [11], Kontoni and Beskos [12], Partridge et al. [13], Agnantiaris et al. [14, 15].

For the Newmark family methods, Newmark [21], the approximations are:

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

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

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: (i) the method is unstable for $\gamma < 1/2$; (ii) unconditional stability occurs when $2\beta \ge \gamma \ge 1/2$, (iii) 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 (13) and (14) correspond to the linear acceleration method. For $\beta = 1/4$ and $\gamma = 1/2$, expressions (13) and (14) correspond to the average acceleration method. These observations can be found in FEM text books, e.g. Bathe [19], Cook et al. [20]. For the D-BEM formulation developed here, however, reliable results were obtained only with the adoption of values for $\beta$ and $\gamma$ not usual in FEM applications.

After substituting (14) in (9), one has:

$$ 
\begin{aligned}
& \begin{bmatrix} (\beta(c\Delta t)^2 \mathbf{H}^{bb} + \mathbf{M}^{bb}) & \mathbf{M}^{bd} \\ (\beta(c\Delta t)^2 \mathbf{H}^{db} + \mathbf{M}^{db}) & (\beta(c\Delta t)^2 \mathbf{I} + \mathbf{M}^{dd}) \end{bmatrix} \begin{bmatrix} \mathbf{u}_{n+1}^b \\ \mathbf{u}_{n+1}^d \end{bmatrix} \\
& = \begin{bmatrix} \beta(c\Delta t)^2 \mathbf{G}^{bb} \\ \beta(c\Delta t)^2 \mathbf{G}^{db} \end{bmatrix} \left\{ \mathbf{p}_{n+1}^b \right\} + \begin{bmatrix} \mathbf{M}^{bb} & \mathbf{M}^{bd} \\ \mathbf{M}^{db} & \mathbf{M}^{dd} \end{bmatrix} \\
& \quad \times \left\{ \mathbf{u}_n^b + \Delta t \dot{\mathbf{u}}_n^b + \frac{(1-2\beta)}{2} \Delta t^2 \ddot{\mathbf{u}}_n^b \right. \\
& \quad \left. \mathbf{u}_n^d + \Delta t \dot{\mathbf{u}}_n^d + \frac{(1-2\beta)}{2} \Delta t^2 \ddot{\mathbf{u}}_n^d \right\}
\end{aligned} \quad (15) $$

Equations 12 and 15 can be represented in a concise manner as:

$$ \bar{\mathbf{H}} \mathbf{u}_{n+1} = \bar{\mathbf{G}} \mathbf{p}_{n+1} + \mathbf{g}_n \quad (16) $$
in which the contributions of previous instants of time are stored in vector $g_n$.

It is important to point out that an adequate choice of the time-step plays a fundamental role in the analysis. A dimensionless variable, say $\beta_{\Delta t}$, was adopted in order to provide a measure of the time-step $\Delta t$, see Mansur [3], Carrer and Mansur [10]:

$$ \beta_{\Delta t} = \frac{c \Delta t}{\ell} \qquad (17) $$

where $\ell$ is the length of the smallest element used in the boundary discretization.

### 3 Initial conditions contributions

In what follows, a description of the procedures adopted to take into account non-homogeneous initial conditions is given in details.

#### 3.1 Houbolt method

In the Houbolt method, the computation of the velocity and accelerations at time $t_{n+1} = (n + 1)\Delta t$ requires the knowledge of the values of $\mathbf{u}$ from time $t_{n-2} = (n - 2)\Delta t$ to time $t_{n+1} = (n + 1)\Delta t$. At the beginning of the analysis, i.e., at the beginning of the time-marching process, $n = 0$ and, consequently, the values $\mathbf{u}_{-2}$ and $\mathbf{u}_{-1}$ must be computed appropriately in order to provide a good start of the analysis.

For the determination of $\mathbf{u}_{-1}$, initially $\dot{\mathbf{u}}_0$ is computed by employing the forward and the backward finite difference formulae at $t = 0$ and by assuming that the finite difference expressions are equal, that is:

$$ \dot{\mathbf{u}}_0 = \frac{\mathbf{u}_1 - \mathbf{u}_0}{\Delta t} = \frac{\mathbf{u}_0 - \mathbf{u}_{-1}}{\Delta t} \qquad (18) $$

Solving Eq. 18 for $\mathbf{u}_1$:

$$ \mathbf{u}_1 = 2\mathbf{u}_0 - \mathbf{u}_{-1} \qquad (19) $$

One can also assume that $\dot{\mathbf{u}}_0$ can be computed by employing a central finite difference formula, which gives:

$$ \dot{\mathbf{u}}_0 = \frac{\mathbf{u}_1 - \mathbf{u}_{-1}}{2\Delta t} \qquad (20) $$

Solving Eq. 20 for $\mathbf{u}_1$:

$$ \mathbf{u}_1 = 2\Delta t \dot{\mathbf{u}}_0 + \mathbf{u}_{-1} \qquad (21) $$

From (19) and (21) one has $\mathbf{u}_{-1}$, as a function of $\dot{\mathbf{u}}_0$ and $\mathbf{u}_0$:

$$ \mathbf{u}_{-1} = \mathbf{u}_0 - \Delta t \dot{\mathbf{u}}_0 \qquad (22) $$

For the determination of $\mathbf{u}_{-2}$, a similar procedure is followed: now $\dot{\mathbf{u}}_{-1}$ is initially computed by employing the forward and the backward finite difference formulae at

$t = -\Delta t$ and assuming that the finite difference results are equal, that is:

$$ \dot{\mathbf{u}}_{-1} = \frac{\mathbf{u}_0 - \mathbf{u}_{-1}}{\Delta t} = \frac{\mathbf{u}_{-1} - \mathbf{u}_{-2}}{\Delta t} \qquad (23) $$

Solving Eq. 23 for $\mathbf{u}_{-2}$:

$$ \mathbf{u}_{-2} = 2\mathbf{u}_{-1} - \mathbf{u}_0 \qquad (24) $$

Substituting (22) in (24) one obtains the value of $\mathbf{u}_{-2}$ as a function of $\dot{\mathbf{u}}_0$ and $\mathbf{u}_0$:

$$ \mathbf{u}_{-2} = \mathbf{u}_0 - 2\Delta t \dot{\mathbf{u}}_0 \qquad (25) $$

#### 3.2 Newmark method

At the beginning of the time-marching process one can write, by taking $n = 0$ in expressions (13) and (14):

$$ \dot{\mathbf{u}}_1 = \frac{\gamma}{\beta \Delta t} [\mathbf{u}_1 - \mathbf{u}_0] + \frac{(\beta - \gamma)}{\beta} \dot{\mathbf{u}}_0 - \frac{(\gamma - 2\beta)}{2\beta} \Delta t \ddot{\mathbf{u}}_0 \qquad (26) $$

$$ \ddot{\mathbf{u}}_1 = \frac{1}{\beta \Delta t^2} [\mathbf{u}_1 - \mathbf{u}_0] - \frac{1}{\beta \Delta t} \dot{\mathbf{u}}_0 - \frac{(1-2\beta)}{2\beta} \ddot{\mathbf{u}}_0 \qquad (27) $$

In expressions (26) and (27), $\ddot{\mathbf{u}}_0$ is the only unknown; in order to determine this value, one can assume that:

$$ \ddot{\mathbf{u}}_0 = \frac{\dot{\mathbf{u}}_1 - \dot{\mathbf{u}}_0}{\Delta t} \qquad (28) $$

After solving (28) for $\dot{\mathbf{u}}_1$ and substituting the resulting expression in (26), one has:

$$ \ddot{\mathbf{u}}_0 = \frac{2}{\Delta t^2} [\mathbf{u}_1 - \mathbf{u}_0] - \frac{2}{\Delta t} \dot{\mathbf{u}}_0 \qquad (29) $$

The substitution of (29) in (27) produces the following expression:

$$ \ddot{\mathbf{u}}_1 = \frac{2}{\Delta t^2} [\mathbf{u}_1 - \mathbf{u}_0] - \frac{2}{\Delta t} \dot{\mathbf{u}}_0 \qquad (30) $$

Bearing Eq. 30 in mind, Eq. 15 at the first time step ($n = 0$) can be particularized and written according to:

$$ \begin{bmatrix} ((c \Delta t)^2 \mathbf{H}^{bb} + 2 \mathbf{M}^{bb}) & 2 \mathbf{M}^{bd} \\ ((c \Delta t)^2 \mathbf{H}^{db} + 2 \mathbf{M}^{db}) & ((c \Delta t)^2 \mathbf{I} + 2 \mathbf{M}^{dd}) \end{bmatrix} \begin{Bmatrix} \mathbf{u}_1^b \\ \mathbf{u}_1^d \end{Bmatrix} = \begin{bmatrix} (c \Delta t)^2 \mathbf{G}^{bb} \\ (c \Delta t)^2 \mathbf{G}^{db} \end{bmatrix} \begin{Bmatrix} \mathbf{p}_1^b \\ \mathbf{p}_1^d \end{Bmatrix} + \begin{bmatrix} 2 \mathbf{M}^{bb} & 2 \mathbf{M}^{bd} \\ 2 \mathbf{M}^{db} & 2 \mathbf{M}^{dd} \end{bmatrix} \begin{Bmatrix} \mathbf{u}_0^b + \Delta t \dot{\mathbf{u}}_0^b \\ \mathbf{u}_0^d + \Delta t \dot{\mathbf{u}}_0^d \end{Bmatrix} \qquad (31) $$

Before proceeding to the next section, a discussion concerning the procedures presented here, summarized by Eqs. 22, 25 and 29, seems to be necessary and can start with the following question: why a classical starting procedure, such as the central difference method, was not employed, as
suggested by Bathe [19]? To answer this question, and also justify the validity of the proposed procedures, it is important to recall how the central difference method is applied as a starting procedure. In this manner, if the central difference method is adopted, it is assumed that:

$$ \ddot{\mathbf{u}}_n = \frac{\mathbf{u}_{n+1} - 2\mathbf{u}_n + \mathbf{u}_{n-1}}{\Delta t^2} \quad (32) $$

and

$$ \dot{\mathbf{u}}_n = \frac{\mathbf{u}_{n+1} - \mathbf{u}_{n-1}}{2\Delta t} \quad (33) $$

Particularizing Eqs. (32) and (33) for $n = 0$, one can arrive at the expression for $\mathbf{u}_{-1}$:

$$ \mathbf{u}_{-1} = \mathbf{u}_0 - \Delta t \dot{\mathbf{u}}_0 + \frac{\Delta t^2}{2} \ddot{\mathbf{u}}_0 \quad (34) $$

Note that $\ddot{\mathbf{u}}_0$ can be computed from Eq. 9, rewritten concisely below for $t = 0$:

$$ \mathbf{H}\mathbf{u}_0 = \mathbf{G}\mathbf{p}_0 - \mathbf{M}\ddot{\mathbf{u}}_0 \quad (35) $$

It is important to note that the values of $\mathbf{p}_0$, in Eq. 35, are directly computed from the analytical expressions of the initial condition $\mathbf{u}_0$. Finally, $\ddot{\mathbf{u}}_0$ is given by:

$$ \ddot{\mathbf{u}}_0 = \mathbf{M}^{-1}[\mathbf{G}\mathbf{p}_0 - \mathbf{H}\mathbf{u}_0] \quad (36) $$

Here appears the main difficulty for employing this starting procedure with the D-BEM: the use of double nodes in the boundary discretization (as is done in Examples 1 and 2) or in the domain (as will be explained next, in Example 4) makes $\mathbf{M}$ a singular matrix and, consequently, $\ddot{\mathbf{u}}_0$ can not be computed from Eq. 36. If the analysis produces a non-singular matrix $\mathbf{M}$, as occurs in the third example, Eq. 36 can be used.

For the use of the Houbolt method it becomes necessary to compute $\mathbf{u}_1$ and $\mathbf{u}_2$; from Eq. 32 one has:

$$ \mathbf{u}_1 = 2\mathbf{u}_0 - \mathbf{u}_{-1} + \Delta t^2 \ddot{\mathbf{u}}_0 \quad (37) $$

The computation of $\mathbf{u}_2$ with the aid of Eq. 32 (by taking $n = 1$) requires the knowledge of $\ddot{\mathbf{u}}_1$ and, consequently, the knowledge of $\mathbf{p}_1$. This observation becomes clearer by writing an equation equivalent to Eq. 36:

$$ \ddot{\mathbf{u}}_1 = \mathbf{M}^{-1}[\mathbf{G}\mathbf{p}_1 - \mathbf{H}\mathbf{u}_1] \quad (38) $$

As $\mathbf{p}_1$ is no longer determined directly, with the Houbolt method this starting procedure fails.

On the other hand, if the Newmark method is chosen, see Eq. 14, one can write:

$$ \ddot{\mathbf{u}}_1 = \frac{1}{\beta \Delta t^2} [\mathbf{u}_1 - \mathbf{u}_0] - \frac{1}{\beta \Delta t} \dot{\mathbf{u}}_0 - \frac{(1-2\beta)}{2\beta} \ddot{\mathbf{u}}_0 \quad (39) $$

and Eq. 15 can be directly used for $n \ge 0$. A comparison between the results provided by the central difference starting

procedure and the procedure developed here for the Newmark method is carried out in the third example of the next section.

## 4 Numerical results

The main concern of the first example and of the first analysis of the second example is to verify the applicability of the Newmark method, Newmark [21], in the D-BEM formulation; in other words, to determine a possible range of variation of the parameters $\beta$ and $\gamma$. The choice of an adequate value of the parameter $\beta_{\Delta t}$ is also discussed. The BEM results are always compared with the corresponding analytical solution, computed following the procedures described by Stephenson [23] and Kreyszig [24]. In order to simplify the notation, analyses carried out by employing the Houbolt or the Newmark methods are referred to as Houbolt or Newmark analyses.

### 4.1 One-dimensional bar with a time variable boundary condition

This example consists of a one-dimensional bar defined over $0 \le x \le a$, fixed at $x = 0$ and subjected to a time variable sinusoidal boundary condition at $x = a$, i.e.

$$ \bar{u}(0, t) = 0 \quad (40) $$

$$ \bar{u}(a, t) = U \sin \lambda t \quad (41) $$

By assuming homogeneous initial conditions, the analytical solution to this problem is given by, see Stephenson [23]:

$$ u(x, t) = U \left[ \frac{x}{a} \sin \lambda t + \sum_{r=1}^{\infty} w_r(t) \sin \left( \frac{r\pi x}{a} \right) \right] \quad (42) $$

where:

$$ w_r(t) = \frac{2\lambda}{(\phi^2 - \lambda^2)} \frac{(-1)^r}{r\pi} [\phi \sin \phi t - \lambda \sin \lambda t] \quad \text{with} \quad \phi = \frac{cr\pi}{a} \quad (43) $$

In the analysis presented here, $\lambda = \pi/24$. A particular feature of this analysis is to be mentioned: as the essential boundary condition at $x = a$ is known (see expression (41)), the expression of the acceleration is easily computed and substitutes the approximations given by (11) and (14) in the systems of Eqs. 12 and 15. The contribution of the known acceleration term is then added directly to the independent vector $\mathbf{g}_n$ in Eq. 16.

Treated as a 2D problem, the original bar is represented by a rectangular domain of length equal to a and height equal to $a/2$. The discretization, with double nodes at the corners, employed 48 boundary elements and 256 internal cells, see Fig. 2.

For the Houbolt analysis, best results were achieved for $\beta_{\Delta t} = 1/3$. The results for the potential at point $C(a/2, a/4)$
<figure><img src="image_4.png" /><figcaption>Fig. 2 One-dimensional bar: boundary and domain discretization</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3 One-dimensional bar with sinusoidal boundary condition: Houbolt analysis, potential at $C(a/2, a/4)$</figcaption></figure>

are presented in Fig. 3; the results for the flux at node $B(0, a/4)$ are presented in Fig. 4.

Bearing in mind that the use of the Newmark method requires the introduction of damping, Carrer and Mansur [10], several analyses were carried out before achieving useful values of the parameters $\beta$ and $\gamma$: for practical purposes, the best choice falls on $\beta = 0.30$ and on $\gamma \ge 0.52$. This choice was conditioned by the accuracy in the results related to the flux at node $B(0, a/4)$; here $\gamma = 0.55$ was adopted. Best results were achieved for $\beta_{\Delta t} = 2/3$. The results for the potential at point $C(a/2, a/4)$ are presented in Fig. 5; the results for the flux at node $B(0, a/4)$ are presented in Fig. 6. As before, good agreement is observed between D-BEM and analytical results, thus confirming the usefulness of the Houbolt method and demonstrating that the Newmark method is suitable for D-BEM analyses.

<figure><img src="image_6.png" /><figcaption>Fig. 4 One-dimensional bar with sinusoidal boundary condition: Houbolt analysis, flux at $B(0, a/4)$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5 One-dimensional bar with sinusoidal boundary condition: Newmark analysis, potential at $C(a/2, a/4)$</figcaption></figure>

## 4.2 One-dimensional rod under a Heaviside-type forcing function and under initial conditions prescribed over the entire domain

This is the classical example of a one-dimensional rod fixed at one side ($x = 0$) and free at the other ($x = a$). Three analyses were performed by adopting the same mesh of the previous example (see Fig. 2): in the first analysis a Heaviside-type
<figure><img src="image_8.png" /><figcaption>Fig. 6 One-dimensional bar with sinusoidal boundary condition: Newmark analysis, flux at $B(0, a/4)$</figcaption></figure>

forcing function $\bar{p}(a) = PH(t-0)$, applied instantaneously at $t=0$ and kept constant in time, is imposed to the free side. The second and the third analyses are carried out by assuming, respectively, initial conditions of the type: $u_0(x) = Ux$, $0 \le x \le a$, and $\dot{u}_0(x) = V$, $0 \le x \le a$.

For the general case, the analytical solution is given by, see Stephenson [23]:

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

Expression (44) can be particularized to each one of the above mentioned analysis by taking the constants $P$, or $U$, or $V$ null.

Houbolt results for the potential at node $A(a, a/4)$ are shown in Figs. 7, 8 and 9. Results for the flux at node $B(0, a/4)$ are shown in Figs. 10, 11 and 12. The time interval was selected by taking $\beta_{\Delta t} = 1/3$.

Newmark results were obtained by choosing $\beta = 0.30$ and $\gamma = 0.60$, as more damping was required for obtaining accurate results for the flux. Results for the potential at node $A(a, a/4)$ are shown in Figs. 13, 14 and 15. Results for the flux at node $B(0, a/4)$ are shown in Figs. 16, 17 and 18. Here, the time interval was selected by taking $\beta_{\Delta t} = 2/3$.

The same conclusions of the previous example are valid, i.e., the good agreement between the analytical and the D-BEM results demonstrates that the Newmark method is suitable for D-BEM analyses, depending on the correct

<figure><img src="image_9.png" /><figcaption>Fig. 7 One-dimensional bar under a Heaviside-type forcing function: Houbolt analysis, potential at $A(a, a/4)$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8 One-dimensional bar under initial displacement field over the entire domain: Houbolt analysis, potential at $A(a, a/4)$</figcaption></figure>

choice of the parameters $\beta$ and $\gamma$. Although the results in Figs. 13, 14 and 15 present more damping than those in Figs. 7, 8 and 9, the presence of damping proved to be positive, leading to the reliable results in Figs. 16, 17 and 18. Further, the second and third analyses show that expressions (22), (25) and (29) are capable of taking into account the contributions of the initial conditions.
<figure><img src="image_11.png" /><figcaption>Fig. 9 One-dimensional bar under initial velocity field over the entire domain: Houbolt analysis, potential at A(a, a/4)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10 One-dimensional bar under a Heaviside-type forcing function: Houbolt analysis, flux at B(0, a/4)</figcaption></figure>

### 4.3 Square membrane under prescribed initial displacement over the entire domain

This example deals with a square membrane subjected to the initial conditions over the domain $0 \le x \le a$, $0 \le y \le a$, see Fig. 19:

$$u_0(x, y) = U x(x - a)y(y - a) \quad (45)$$

$$\dot{u}_0(x, y) = 0 \quad (46)$$

<figure><img src="image_13.png" /><figcaption>Fig. 11 One-dimensional bar under initial displacement field over the entire domain: Houbolt analysis, flux at B(0, a/4)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12 One-dimensional bar under initial velocity field over the entire domain: Houbolt analysis, flux at B(0, a/4)</figcaption></figure>

The general analytical solution to this problem, for a rectangular membrane with dimensions a and b, according to Kreyszig [24], is:

$$\begin{aligned} u(x, y, t) = & \frac{64Ua^2b^2}{\pi^6} \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} \frac{1}{m^3n^3} \sin\left(\frac{m\pi x}{a}\right) \\ & \times \sin\left(\frac{n\pi y}{b}\right) \cos(\lambda_{mn}\pi ct) \end{aligned} \quad (47)$$
<figure><img src="image_15.png" /><figcaption>Fig. 13 One-dimensional bar under a Heaviside-type forcing function: Newmark analysis, potential at A(a, a/4)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14 One-dimensional bar under initial displacement field over the entire domain: Newmark analysis, potential at A(a, a/4)</figcaption></figure>

where:

$$ \lambda_{mn} = \sqrt{\frac{m^2}{a^2} + \frac{n^2}{b^2}} \qquad (48) $$

In this analysis, the boundary discretization employed 80 elements, without double nodes at the corners due to the symmetry of the problem, and the square domain, 800 cells, see Fig. 20.

<figure><img src="image_17.png" /><figcaption>Fig. 15 One-dimensional bar under initial velocity field over the entire domain: Newmark analysis, potential at A(a, a/4)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16 One-dimensional bar under a Heaviside-type forcing function: Newmark analysis, flux at B(0, a/4)</figcaption></figure>

The results corresponding to the displacement at point A(a/2, a/2) and to the support reaction at node B(a, a/2), from the Houbolt analyses, are shown in Figs. 21 and 22, respectively.

Newmark analyses were carried out by adopting $\beta = 0.30$ and $\gamma = 0.52$ and the results are shown in Figs. 23 and 24.

In this example the presence of a large amount of damping is not so important as it was in the previous ones; for this reason, a smaller value to the $\gamma$ parameter was adopted.
<figure><img src="image_19.png" /><figcaption>**Fig. 17** One-dimensional bar under initial displacement field over the entire domain: Newmark analysis, flux at $B(0, a/4)$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 18** One-dimensional bar under initial velocity field over the entire domain: Newmark analysis, flux at $B(0, a/4)$</figcaption></figure>

Houbolt and Newmark analyses were carried out with the same time interval, computed from $\beta_{\Delta t} = 3/10$.

The good agreement between D-BEM and analytical results, depicted in Figs. 21, 22, 23 and 24, was already expected, as in this example time jumps do not appear.

Before finishing this example, a comparison between the results provided by the Newmark method, as proposed in this work, with the results provided by the central difference method (as a starting method), is presented in Fig. 25,

<figure><img src="image_21.png" /><figcaption>**Fig. 19** Square membrane under initial displacement field over the entire domain</figcaption></figure>

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

for the displacement at point $A(a/2, a/2)$, and in Fig. 26, for the support reaction at node $B(a, a/2)$. In this comparison, the time interval was adopted by taking $\beta_{\Delta t} = 3/10$. The good agreement observed between both the numerical results attests the applicability of the proposed formulation to D-BEM analyses.

4.4 Square membrane under prescribed initial velocity over part of the domain

The square membrane depicted in Fig. 27, with an initial velocity field $\dot{u}_0(x, y) = V$ prescribed over the sub-domain
<figure><img src="image_23.png" /><figcaption>Fig. 21 Square membrane under initial displacement field over the entire domain: displacement at point A(a/2, a/2): Houbolt analysis</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 22 Square membrane under initial displacement field over the entire domain: support reaction at node B(a, a/2): Houbolt analysis</figcaption></figure>

$\Omega_0$ and with zero displacements prescribed all over the domain, is analysed in this example.

The analytical solution to this problem is given by, see Kreyszig [24]:

$$u(x, y, t) = \frac{4V}{c\pi^3} \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} \frac{1}{mn\lambda_{mn}} \sin\left(\frac{m\pi x}{a}\right) \\ \times \sin\left(\frac{n\pi y}{a}\right) G_{mn} \quad (49)$$

<figure><img src="image_25.png" /><figcaption>Fig. 23 Square membrane under initial displacement field over the entire domain: displacement at point A(a/2, a/2): Newmark analysis</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 24 Square membrane under initial displacement field over the entire domain: support reaction at node B(a, a/2): Newmark analysis</figcaption></figure>

where:

$$\lambda_{mn} = \frac{\sqrt{m^2 + n^2}}{a} \quad (50)$$

and

$$G_{mn} = \left( \cos\left(\frac{3n\pi}{5}\right) - \cos\left(\frac{2n\pi}{5}\right) \right) \\ \times \left( \cos\left(\frac{3m\pi}{5}\right) - \cos\left(\frac{2m\pi}{5}\right) \right) \sin(\lambda_{mn}\pi ct) \quad (51)$$
<figure><img src="image_27.png" /><figcaption>**Fig. 25** Square membrane under initial displacement field over the entire domain: displacement at point A(a/2, a/2): Newmark analysis and central difference method</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 26** Square membrane under initial displacement field over the entire domain: support reaction at node B(a, a/2): Newmark analysis and central difference method</figcaption></figure>

This example was analysed previously with the use of the TD-BEM formulation, e.g. Mansur [3] and Carrer and Mansur [9] and, more recently, with the use of a BEM formulation based on the convolution quadrature method, see Abreu et al. [25]. In the D-BEM formulation presented here, care should be taken when considering the initial conditions contributions. Because linear interpolation is adopted for the initial conditions in the internal cells, the use of double-nodes

<figure><img src="image_29.png" /><figcaption>**Fig. 27** Square membrane under prescribed initial velocity over part of the domain</figcaption></figure>

in the boundary of Ω₀ became necessary in order to avoid the spreading of the initial conditions to the cells in the neighbouring of Ω₀. Aiming at providing a good representation of the jumps in the support reaction response, a more refined mesh, constituted of 160 elements and 3,200 cells, not shown here, was employed. The results of the Houbolt analyses, corresponding to the displacement at point A(a/2, a/2) and to the support reaction at node B(a, a/2), are depicted in Figs. 28 and 29, respectively. These results were obtained by adopting β<sub>Δt</sub> = 3/10. The Newmark analyses was carried out with β = 0.30, γ = 0.52 and β<sub>Δt</sub> = 3/5; the results are depicted in Figs. 30 and 31. The D-BEM responses are in good agreement with the analytical ones, although it must be pointed out that the TD-BEM formulation provides a better representation of the jumps in the support reaction response. Nevertheless, the good responses furnished by the D-BEM formulation demonstrate its applicability to the solution of this kind of problems.

## 5 Conclusions

The D-BEM formulation is a very promising approach and, consequently, it is expected that some research work concerning its development will be done in the next years. With the purpose of contributing with the development of the D-BEM formulation, this work in concerned with two subjects: the first one is the search for alternative methods to perform the march in time, as the Houbolt method has been, during the last years, the only one successfully employed in the D-BEM formulation. For this reason, the Newmark
<figure><img src="image_30.png" /><figcaption>**Fig. 28** Square membrane under prescribed initial velocity over part of the domain: displacement at point A(a/2, a/2): Houbolt analysis</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Fig. 29** Square membrane under prescribed initial velocity over part of the domain: support reaction at node B(a, a/2): Houbolt analysis</figcaption></figure>

method was adopted to carry out this discussion and, on the basis of the examples presented here, it seems that this goal was achieved, as stable and accurate results were found. The second subject is concerned with the solution of problems with non-homogeneous initial conditions; in other words, it is concerned with finding general expressions to the terms $\mathbf{u}_{-2}$ and $\mathbf{u}_{-1}$, in the Houbolt method, and to the term $\mathbf{u}_0$, in the Newmark method. The results presented here demonstrate the validity of the expressions found and encourage the

<figure><img src="image_32.png" /><figcaption>**Fig. 30** Square membrane under prescribed initial velocity over part of the domain: displacement at point A(a/2, a/2): Newmark analysis</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Fig. 31** Square membrane under prescribed initial velocity over part of the domain: support reaction at node B(a, a/2): Newmark analysis</figcaption></figure>

extension of the present work to elastodynamics, as well as to the DR-BEM formulation and also to the Finite Element Method.

## References

1. Beskos DE (1997) Boundary elements in dynamic analysis: part II (1986–1996). Appl Mech Rev 50:149–197
2. Beskos DE (2003) Dynamic analysis of structures and structural systems. In: Beskos DE, Maier G (eds) Boundary element advances in solid mechanics. CISM, Udine

3. Mansur WJ (1983) A Time-stepping technique to solve wave propagation problems using the boundary element method. Ph.D. Thesis, University of Southampton

4. Dominguez J (1993) Boundary elements in dynamics. Computational Mechanics Publications, Southampton, Boston

5. 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

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

7. Mansur WJ, deLima-Silva W (1992) Efficient time truncation in two-dimensional BEM analysis of transient wave propagation problems. Earthquake Eng Struct Dyn 21:51-63

8. Soares Jr D, Mansur WJ (2004) Compression of time generated matrices in two-dimensional time-domain elastodynamic BEM analysis. Int J Numer Methods Eng 61:1209-1218

9. Carrer JAM, Mansur WJ (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

10. Carrer JAM, Mansur WJ (2004) Alternative time-marching schemes for elastodynamic analysis with the domain boundary element method formulation. Comput Mech 34:387-399

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

12. Kontoni DPN, Beskos DE (1993) Transient dynamic elastoplastic analysis by the dual reciprocity BEM. Eng Anal Bound Elements 12:1-16

13. Partridge PW, Brebbia CA, Wrobel LC (1992) The dual reciprocity boundary element method. Computational Mechanics Publications, Southampton, Boston

14. Agnantiaris JP, Polyzos D, Beskos DE (1996) Some studies on dual reciprocity BEM for elastodynamic analysis. Comput Mech 17:270-277

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

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

17. 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

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

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

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

21. Newmark NM (1959) A method of computation for structural dynamics. ASCE J Eng Mech Div 85:67-94

22. Carrer JAM, Mansur WJ (1996) Time-Domain BEM analysis for the 2D scalar wave equation: initial conditions contributions to space and time derivatives. Int J Numer Methods Eng 39:2169-2188

23. Stephenson G (1970) An introduction to partial differential equations for science students. Longman, London

24. Kreyszig E (1999) Advanced engineering mathematics, 8th edn. Wiley, New York

25. Abreu AI, Mansur WJ, Carrer JAM (2006) Initial conditions contributions in a BEM formulation based on the convolution quadrature method. Int J Numer Methods Eng 67:417-434