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

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.


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

Operational quadrature method
The basic steps of the OQM [11 -13] are outlined in what follows. The convolution integral is written below where the inner ' p ' denotes convolution, and u p ðtÞ represents, in BEM notation, the fundamental solution to the problem. The solution to this problem can be represented as: In expression (2), v n represents the quadrature weights (or integration weights) that constitute the coefficients of the power series that approximate the Laplace transformû p ðsÞ of the fundamental solution u p ðtÞ; that iŝ where s ¼ gðzÞ=Dt and z is a complex variable. It is also required that lû p ðsÞl ! 0 when lsl ! 1 for RðsÞ $ q (q is a real positive number). The coefficients of the series in expression (3) are furnished by the Cauchy integral formula given below: where r is the radius of a circle in the domain of analyticity of the functionû p ðsÞ: If a polar coordinate system is adopted, the integral presented in expression (4) is approximated by the trapezoidal rule with L equal steps (2p/L ) [12]. The function gð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 and then The function gðzÞ clearly characterises the multistep method and must be AðaÞ-stable with positive angle a, stable in a neighbourhood of infinity, strongly zero-stable, and consistent of order p. If an error d is assumed in the computation ofû p ðsÞ in Eq. (4), the choice L ¼ N and r N ¼ ffiffi d p leads to an error in v n of order Oð ffiffi d p Þ:

Boundary element formulation
The scalar wave equation is written as [18] 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 G ¼ G u < G p ): 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 ›uðx; tÞ=›tl t¼0 ¼ 0 in V < G (V stands for the domain).
The TD-BEM equation that corresponds to Eq. (7) can be written as [1] 4pcðjÞuðj where u p ðx; t; j; tÞ is the fundamental solution, p p ðx; t; j; tÞ ¼ ›u p ðx; t; j; tÞ=›n; and the coefficient cðjÞ is the same as that of the static problem. Following Eq. (2), the convolutions implicitly indicated in Eq. (10) are written as ðt n 0 u p ðx; t; j; tÞpðx; tÞdt ¼ X n k¼0 g j n2k ðx; j; DtÞp j k ðxÞ; where it was already assumed the boundary discretised in J elements G j . Boundary variables uðx; tÞ and pðx; tÞ; for each element, are written as: The weights g n and 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 where F 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 r N ¼ ffiffi d p with an error d ¼ 10 210 : The function gðzÞ of order p, was taken as [11 -13]: The fundamental solution for 2D problems in the Laplace domain, and its normal derivative, in Eq. (13), are given by [18] u p ðr; sÞ ¼ 2K 0 ðsrÞ;p p ðr; sÞ ¼ 22sK 1 ðsrÞ ›r ›n ð15Þ where r ¼ lx 2 jl is the distance between the source (j ) 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 j i can be written as: After writing Eq. (16) for all nodes of the discretised boundary, the following representation arises: In Eq. (17), C is a diagonal matrix that contains the coefficients cðj i Þ and n and k correspond to the discrete times t n ¼ nDt and t k ¼ kDt; respectively. In order to solve the problem numerically, after imposing the boundary conditions, one has the following general expression for t n ¼ nDt : 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:

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 b ¼ cðDtÞ=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 gðzÞ was taken as the first and second order polynomial functions and, according to Eq. (14), can be written as: for first order polynomial ð20Þ

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 2 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 b ¼ 0:6 and b ¼ 0:3: It is observed that the use of a small value for b produced a more accurate result diversely  from the TD-BEM formulation [1] for which the best choice was b ¼ 0: 6: For the results in Fig. 2b, b ¼ 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 b ¼ 0:3 and for b ¼ 0:054; respectively.

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, b ¼ 0:6 was adopted. Fig. 5b presents the same comparison between the first order polynomial and the TD-BEM formulation and, again, b ¼ 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.

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 gð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 b 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.