# Scalar wave propagation in 2D: a BEM formulation based on the operational quadrature method

Ana Ibis Abreu, Jose Antonio Marques Carrer, Webe João Mansur

► To cite this version:

Ana Ibis Abreu, Jose Antonio Marques Carrer, Webe João Mansur. Scalar wave propagation in 2D: a BEM formulation based on the operational quadrature method. Engineering Analysis with Boundary Elements, 2003, 27 (2), pp.101-105. 10.1016/S0955-7997(02)00087-5. hal-01580922

HAL Id: hal-01580922

https://hal.science/hal-01580922

Submitted on 3 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
# Scalar wave propagation in 2D: a BEM formulation based on the operational quadrature method

A.I. Abreu, J.A.M. Carrer, W.J. Mansur

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

This work presents a boundary element method formulation for the analysis of scalar wave propagation problems. The formulation presented here employs the so-called operational quadrature method, by means of which the convolution integral, presented in time-domain BEM formulations, is substituted by a quadrature formula, whose weights are computed by using the Laplace transform of the fundamental solution and a linear multistep method. Two examples are presented at the end of the article with the aim of validating the formulation.

Keywords: Boundary element; Wave propagation; Dynamic

## 1. Introduction

Due to the importance and wide range of application of the analyses concerning scalar wave propagation and elastodynamics, several BEM formulations were developed during the last years and, for the purposes of the present work, they will be classified in time-domain and transformed-domain formulations.

Time-domain formulations (TD-BEM) [1-3] provide good representation of the causality and lead to very accurate results. Besides, the fulfilment of the radiation condition makes them suitable for infinite domain analysis. Transformed-domain formulations [4-7], on the other hand, initially solve the problem in a transformed-domain (frequency domain or Laplace domain) and determine the answer in the time-domain (when desired) by means of suitable transformation formulae [8-10].

The BEM formulation developed in this work for the solution of problems governed by the scalar wave equation employs the fundamental solution in the Laplace domain and provides direct solution to the problem in the time-domain. This characteristic of the formulation turns it very attractive. The approach is based on the operational quadrature method (OQM) [11-13], and will be referred

to along the text as OQM-BEM. The formulation developed in the present work has already been successfully employed in 3D elastodynamic BEM formulations [14] and in 3D viscoelastic BEM formulations [15,16]. Poroelastic applications can be found in the recent book by Schanz [17].

In Section 2, the basic steps of the OQM will be presented, as well as its extension to the BEM. In order to solve the problem, linear boundary elements are employed. Two examples are analysed, with the aim of validating the extension of the formulation to the BEM and of verifying the accuracy of the results that it can produce.

## 2. Operational quadrature method

The basic steps of the OQM [11-13] are outlined in what follows. The convolution integral is written below

$$y(t) = \int_{0}^{t} u^{*}(t - \tau)f(\tau)d\tau = u^{*}(t) * f(t) \quad (1)$$

where the inner '*' denotes convolution, and $u^{*}(t)$ represents, in BEM notation, the fundamental solution to the problem. The solution to this problem can be represented as:

$$y(n\Delta t) = \sum_{k=0}^{n} \omega_{n-k}(\Delta t)f(k\Delta t), \quad n = 0, 1, \dots, N \quad (2)$$
In expression (2), $\omega_n$ represents the quadrature weights (or integration weights) that constitute the coefficients of the power series that approximate the Laplace transform $\hat{u}^*(s)$ of the fundamental solution $u^*(t)$, that is

$$ \hat{u}^*\left(\frac{\gamma(z)}{\Delta t}\right) = \sum_{n=0}^{\infty} \omega_n(\Delta t) z^n \quad (3) $$

where $s = \gamma(z)/\Delta t$ and $z$ is a complex variable. It is also required that $|\hat{u}^*(s)| \to 0$ when $|s| \to \infty$ for $\Re(s) \ge q$ ($q$ is a real positive number).

The coefficients of the series in expression (3) are furnished by the Cauchy integral formula given below:

$$ \begin{aligned} \omega_n(\Delta t) &= \frac{1}{2\pi i} \int_{|z|=\rho} \hat{u}^*\left(\frac{\gamma(z)}{\Delta t}\right) z^{-n-1} \, dz \\ &= \frac{\rho^{-n}}{L} \sum_{l=0}^{L-1} \hat{u}^* \left( \frac{\gamma(\rho e^{i2\pi L/L})}{\Delta t} \right) e^{-inL\pi/L} \end{aligned} \quad (4) $$

where $\rho$ is the radius of a circle in the domain of analyticity of the function $\hat{u}^*(s)$. If a polar coordinate system is adopted, the integral presented in expression (4) is approximated by the trapezoidal rule with $L$ equal steps $(2\pi/L)$ [12].

The function $\gamma(z)$, previously utilized in Eqs. (3) and (4), is the quotient of the polynomials generated by a linear multistep method: if this method is employed for approximating a certain function, say $x(t)$, which, by its turn, is the solution of the first order differential equation $dx(t)/dt = sx(t) + f(t)$ (with $x(0) = 0$), one has

$$ x(t) \approx \sum_{j=0}^{k} \alpha_j x_{n-j} = \Delta t \sum_{j=0}^{k} \beta_j (sx_{n-j} + f((n-j)\Delta t)) \quad (5) $$

and then

$$ \gamma(z) = \frac{\alpha_0 + \cdots + \alpha_k z^k}{\beta_0 + \cdots + \beta_k z^k} \quad (6) $$

The function $\gamma(z)$ clearly characterises the multistep method and must be $A(\alpha)$-stable with positive angle $\alpha$, stable in a neighbourhood of infinity, strongly zero-stable, and consistent of order $p$. If an error $\delta$ is assumed in the computation of $\hat{u}^*(s)$ in Eq. (4), the choice $L = N$ and $\rho^N = \sqrt{\delta}$ leads to an error in $\omega_n$ of order $O(\sqrt{\delta})$.

### 3. Boundary element formulation

The scalar wave equation is written as [18]

$$ \nabla^2 u(x, t) - \frac{1}{c^2} \frac{\partial^2 u(x, t)}{\partial t^2} = 0 \quad (7) $$

where $u(x, t)$ represents the potential, $c$ is the propagation velocity, and $t$ is the time. The boundary conditions are given by (note that $\Gamma = \Gamma_u \cup \Gamma_p$):

$$ u(x, t) = \bar{u}(x, t), \quad x \in \Gamma_u \quad (8) $$

$$ p(x, t) = \frac{\partial u(x, t)}{\partial \mathbf{n}} = \bar{p}(x, t), \quad x \in \Gamma_p \quad (9) $$

In the above expression, $p(x, t)$ is the flux. Additionally, for the formulation developed here, one has zero initial conditions, i.e. $u(x, 0) = 0$ and $\partial u(x, t)/\partial t|_{t=0} = 0$ in $\Omega \cup \Gamma$ ($\Omega$ stands for the domain).

The TD-BEM equation that corresponds to Eq. (7) can be written as [1]

$$ \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} \int_{0}^{t^{+}} p^{*}(x, t; \xi, \tau)u(x, \tau)d\tau d\Gamma(x) \end{aligned} \quad (10) $$

where $u^*(x, t; \xi, \tau)$ is the fundamental solution, $p^*(x, t; \xi, \tau) = \partial u^*(x, t; \xi, \tau)/\partial \mathbf{n}$, and the coefficient $c(\xi)$ is the same as that of the static problem. Following Eq. (2), the convolutions implicitly indicated in Eq. (10) are written as

$$ \int_{0}^{t_n} u^*(x, t; \xi, \tau)p(x, \tau)d\tau = \sum_{k=0}^{n} \mathbf{g}_{n-k}^{j}(x, \xi, \Delta t)\mathbf{p}_{k}^{j}(x), \quad (11) $$

$n = 0, 1, \dots, N$

$$ \int_{0}^{t_n} p^*(x, t; \xi, \tau)u(x, \tau)d\tau = \sum_{k=0}^{n} \mathbf{h}_{n-k}^{j}(x, \xi, \Delta t)\mathbf{u}_{k}^{j}(x), $$

$n = 0, 1, \dots, N$

where it was already assumed the boundary discretised in $J$ elements $\Gamma_j$. Boundary variables $u(x, t)$ and $p(x, t)$, for each element, are written as:

$$ \mathbf{p}_k^j(x) = \mathbf{p}^j(x, k\Delta t), \quad \mathbf{u}_k^j(x) = \mathbf{u}^j(x, k\Delta t) \quad (12) $$

The weights $\mathbf{g}_n$ and $\mathbf{h}_n$ in Eq. (11) are computed by following a similar pattern as that presented in Eq. (4). The resulting expressions for the OQM-BEM formulation are given below

$$ \mathbf{g}_n^j(x, \xi, \Delta t) = \frac{\rho^{-n}}{L} \sum_{l=0}^{L-1} \int_{\Gamma_l} \hat{u}^* r(x, \xi); \frac{\gamma(\rho e^{i2\pi L/L})}{\Delta t} \times \mathbf{\Phi}_j(x) d\Gamma(x) e^{-inL\pi/L} \quad (13) $$

$$ \mathbf{h}_n^j(x, \xi, \Delta t) = \frac{\rho^{-n}}{L} \sum_{l=0}^{L-1} \int_{\Gamma_l} \hat{p}^* r(x, \xi); \frac{\gamma(\rho e^{i2\pi L/L})}{\Delta t} \times \mathbf{\Phi}_j(x) d\Gamma(x) e^{-inL\pi/L} \quad (13) $$

where $\mathbf{\Phi}_j(x)$ represents the interpolation functions utilized in the boundary discretisation (in fact, linear interpolation functions were employed).

In the analyses performed in the present work, it was adopted $L = N$ and $\rho^N = \sqrt{\delta}$ with an error $\delta = 10^{-10}$. The function $\gamma(z)$ of order $p$, was taken as [11–13]:

$$ \gamma(z) = \sum_{n=1}^{p} \frac{1}{n} (1-z)^n \quad (14) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. One-dimensional rod analysis: (a) geometry and boundary conditions; (b) boundary discretisation and selected boundary nodes and internal points.</figcaption></figure>

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

The fundamental solution for 2D problems in the Laplace domain, and its normal derivative, in Eq. (13), are given by [18]

$$ \hat{u}^*(r, s) = 2K_0 sr, \quad \hat{p}^*(r, s) = -2sK_1 sr \frac{\partial r}{\partial \mathbf{n}} \quad (15) $$

where $r = |x - \xi|$ is the distance between the source $(\xi)$ and field $(x)$ points, and $K_0 sr$ is the modified Bessel function of order zero and second type and $K_1 sr$ is the modified Bessel function of first order and second type.

By employing Eqs. (11) and (13), the discrete version of Eq. (10) for the source boundary node $\xi_i$ can be written as:

$$ c(\xi_i)u(\xi_i, t_n) = \sum_{j=1}^{J} \sum_{k=0}^{n} \mathbf{g}_{n-k}^{j}(r, \Delta t) \mathbf{p}_{k}^{j}(x) \\ - \sum_{j=1}^{J} \sum_{k=0}^{n} \mathbf{h}_{n-k}^{j}(r, \Delta t) \mathbf{u}_{k}^{j}(x) \quad (16) $$

After writing Eq. (16) for all nodes of the discretised boundary, the following representation arises:

$$ C\mathbf{u}^n = \sum_{k=0}^{n} G^{n-k} \mathbf{p}^k - \sum_{k=0}^{n} H^{n-k} \mathbf{u}^k \quad (17) $$

In Eq. (17), $C$ is a diagonal matrix that contains the coefficients $c(\xi_i)$ and $n$ and $k$ correspond to the discrete times $t_n = n\Delta t$ and $t_k = k\Delta t$, respectively.

In order to solve the problem numerically, after imposing the boundary conditions, one has the following general

expression for $t_n = n\Delta t$:

$$ A^0 y^n = f^n + \sum_{k=0}^{n-1} f^k \quad (18) $$

In expression (18), unknowns and contribution of the boundary conditions at time $t_n$ are stored, respectively, in vectors $y^n$ and $f^n$. Vector $f^k$ contains the contribution of the history and is given by:

$$ f^k = G^{n-k} p^k - H^{n-k} u^k \quad (19) $$

## 4. Examples

Two numerical applications are presented next in order to verify the accuracy of the numerical results provided by the formulation presented in this work. In all analyses, reference will be made to the dimensionless variable $\beta = c(\Delta t)/l$ that gives a measure of the time-step length to be adopted to perform the analysis ($l$ is the smaller boundary element length).

The function $\gamma(z)$ was taken as the first and second order polynomial functions and, according to Eq. (14), can be written as:

$$ \gamma(z) = 1 - z \quad \text{for first order polynomial} \quad (20) $$

$$ \gamma(z) = \frac{3}{2} - 2z + \frac{1}{2}z^2 \quad \text{for second order polynomial} \quad (21) $$

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

<figure><img src="image_5.png" /><figcaption>Fig. 2. Flux at boundary node $E(a, b/2)$ for the one-dimensional rod analysis: (a) second order polynomial with $\beta = 0.6$ and $\beta = 0.3$; (b) first order polynomial with $\beta = 0.054$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 3. Potential at points A(0, b/2), B(a/4, b/2), C(a/2, b/2) and D(3a/4, b/2) for the one-dimensional rod analysis: (a) second order polynomial with $\beta = 0.3$; (b) first order polynomial with $\beta = 0.054$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Circular cavity analysis: (a): geometry and boundary conditions; (b): boundary discretisation and selected boundary node and internal points.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Potential at points A(R, 0), B(2R, 0) and C(3R, 0) for the circular cavity: (a) second order polynomial with $\beta = 0.6$; (b) first order polynomial and $\beta = 0.6$.</figcaption></figure>

### 4.1. One-dimensional rod under a Heaviside-type forcing function

This first example consists of a one-dimensional rod under a Heaviside-type forcing function applied at $t = 0$ and kept constant from this time onwards, $p(0, 0) = P/EH(t - 0)$ ($E$ is longitudinal elasticity modulus) according to Fig. 1a. The boundary element mesh depicted in Fig. 1b also presents the selected boundary nodes and

internal points at which the computed numerical responses will be compared with the analytical ones.

Flux time history at the boundary node $E(a, b/2)$ is depicted in Fig. 2a and b for the second and first order polynomials, according to expressions (20) and (21), respectively.

Fig. 2a presents the results of two analyses, carried out with $\beta = 0.6$ and $\beta = 0.3$. It is observed that the use of a small value for $\beta$ produced a more accurate result diversely
from the TD-BEM formulation [1] for which the best choice was $\beta = 0.6$.

For the results in Fig. 2b, $\beta = 0.054$ was adopted. Fig. 3a and b show the potential time-histories at boundary node A(0, b/2) and at internal points B(a/4, b/2), C(a/2, b/2) and D(3a/4, b/2), for $\beta = 0.3$ and for $\beta = 0.054$, respectively.

## 4.2. Circular cavity under a Heaviside-type forcing function

This example consists of a circular cavity with radius $R$ in an unbounded medium, under a internal pressure applied at $t = 0$ and kept constant from this time onwards, see Fig. 4a. The boundary element mesh is shown in Fig. 4b.

In Fig. 5a, the potential time-histories at boundary node A(R, 0) and internal points B(2R, 0) and C(3R, 0) furnished by the second order polynomial, are compared with the solution provided by the TD-BEM formulation [1]. For both formulations, $\beta = 0.6$ was adopted. Fig. 5b presents the same comparison between the first order polynomial and the TD-BEM formulation and, again, $\beta = 0.6$ was adopted. The results in Fig. 5a and b are in good agreement with the results provided by the TD-BEM formulation and, probably due to the absence of reflections, a large time-step could be adopted for the first order polynomial.

## 5. Conclusions

The formulation developed in this work, based on the OQM, shows an elegant way of solving the wave propagation problem, provided one has the corresponding fundamental solution in the Laplace domain. It also seems to be very promising, due to the quite simple extension to a large class of problems, such as 2D elastodynamic.

In order to assess the accuracy of the formulation, the first and second order polynomials have been adopted to represent the function $\gamma(z)$ that characterises the multistep method. In the first example, it was shown the dependence of the results on the correct choice of the time-step length, and also the convergence of them, as the dimensionless parameter $\beta$ reaches its optimal value. The second order polynomial was adopted at this point of the study. Then, the first order polynomial was adopted, also providing reliable results. Due to the absence of reflected waves, the cavity analysis presents certain independence on the choice of the time-step length. For this reason, the same time-step was adopted for both polynomials.

Improvements to the present formulation are presently being investigated by the authors and these include

the extension to 2D elastodynamic and poroelastodynamic problems.

## References

[1] Mansur WJ. A time-stepping technique to solve wave propagation problems using the boundary element method. PhD Thesis. University of Southampton, Southampton, UK; 1983.

[2] Dominguez J. Boundary elements in dynamics. Southampton: Computational Mechanics Publications; 1993.

[3] Antes H, von Estorff O. On FEM-BEM coupling for fluid-structure interaction analyses in time-domain. Int J Numer Meth Engng 1991; 31:1151-68.

[4] Cruse TA, Rizzo FJ. A direct formulation and numerical solution of the general transient elastodynamic problem. I. J Math Anal Appl 1968;22:244-59.

[5] Cruse TA. A direct formulation and numerical solution of the general transient elastodynamic problem. II. J Math Anal Appl 1968;22: 341-55.

[6] Manolis GD. A comparative study on three boundary element method approaches to problems in elastodynamics. Int J Numer Meth Engng 1983;19:73-91.

[7] Manolis GD, Beskos DE. Dynamic stress concentration studies by boundary integrals and Laplace transform. Int J Numer Meth Engng 1981;17:573-99.

[8] Dubner H, Abate J. Numerical investigation of Laplace transforms by relating them to the finite Fourier cosine transform. J Assoc Comput Mach 1983;15:115-23.

[9] Durbin F. Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate's method. Comput J 1974;17: 371-6.

[10] Narayanan GV, Beskos DE. Numerical operational methods for time-dependent linear problems. Int J Numer Meth Engng 1982;18: 1829-54.

[11] Lubich C. Convolution quadrature and discretized operational calculus. I. Numer Math 1988;52:129-45.

[12] Lubich C. Convolution quadrature and discretized operational calculus. II. Numer Math 1988;52:413-25.

[13] Lubich C. On the multistep time discretization of linear initial-boundary value problems and their boundary integral equations. Numer Math 1994;67:365-89.

[14] Schanz M, Antes H. Application of 'Operational Quadrature Methods' in time domain boundary element methods. Meccanica 1997;32:179-86.

[15] Schanz M, Antes H. New visco- and elastodynamic time domain boundary element formulation. Comput Mech 1997;20:452-9.

[16] Gaul L, Schanz M. A comparative study of three boundary element approaches to calculate the transient response of viscoelastic solids with unbounded domains. Comput Meth Appl Mech Engng 1999;179: 111-23.

[17] Schanz M. Wave propagation in viscoelastic and poroelastic continua. Lectures notes in applied mechanics, vol. 2. Berlin: Springer; 2001.

[18] Morse PM, Feshbach H. Methods of theoretical physics. New York: McGraw-Hill; 1953.