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

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.


Introduction
The development of BEM formulations for solving timedependent 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 timedomain) , 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

Scalar wave equation by the boundary element method: a D-BEM approach with non-homogeneous initial conditions
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 timemarching 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 β and γ , 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.

D-BEM formulation
The scalar wave equation for 2D problems over a domain limited by a boundary is given by: 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 ≥ 0, are given by: u(X, t) = u(X, t) on u (essential or Dirichlet boundary condition) 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 = u ∪ p The initial conditions, over the domain , are: In the BEM literature, the solution of the Poisson equation: where δ(X − ξ) 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 ξ . The fundamental solution is: 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 * (ξ, 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: The coefficient c(ξ ) is computed according to: where α is the angle depicted in Fig. 1. The function p * (ξ ; X ) is the normal derivative of the fundamental solution: 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 α Ω Γ ξ Fig. 1 Internal angle for the computation of c(ξ ) 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: In Eq. 9, in order to simplify the notation, the subscript (n + 1) represents the time t n+1 = (n + 1) t, 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 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(ξ ) = 1 of the internal points.
From Eq. 9, the boundary unknowns are the potential at p and the flux at 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 M bd in Eq. 9.
For the Newmark family methods, Newmark [21], the approximations are: The stability and accuracy of the Newmark method depend on the correct choice of the parameters β and γ . Some remarks concerning the Newmark method are useful: (i) the method is unstable for γ < 1/2; (ii) unconditional stability occurs when 2β ≥ γ ≥ 1/2, (iii) taking γ > 1/2 introduces artificial damping but reduces the accuracy of the method to first order. If one takes β = 1/6 and γ = 1/2, expressions (13) and (14) correspond to the linear acceleration method. For β = 1/4 and γ = 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 β and γ not usual in FEM applications.
After substituting (14) in (9), one has: Equations 12 and 15 can be represented in a concise manner as: (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 β t , was adopted in order to provide a measure of the time-step t, see Mansur [3], Carrer and Mansur [10]: where is the length of the smallest element used in the boundary discretization.

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

Houbolt method
In the Houbolt method, the computation of the velocity and accelerations at time t n+1 = (n + 1) t requires the knowledge of the values of u from time t n−2 = (n − 2) t to time t n+1 = (n + 1) t. At the beginning of the analysis, i.e., at the beginning of the time-marching process, n = 0 and, consequently, the values u −2 and u −1 must be computed appropriately in order to provide a good start of the analysis. For the determination of u −1 , initiallyu 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: Solving Eq. 18 for u 1 : One can also assume thatu 0 can be computed by employing a central finite difference formula, which gives: Solving Eq. 20 for u 1 : From (19) and (21) one has u −1 , as a function ofu 0 and u 0 : For the determination of u −2 , a similar procedure is followed: nowu −1 is initially computed by employing the forward and the backward finite difference formulae at t = − t and assuming that the finite difference results are equal, that is: Solving Eq. 23 for u −2 : Substituting (22) in (24) one obtains the value of u −2 as a function ofu 0 and u 0 :

Newmark method
At the beginning of the time-marching process one can write, by taking n = 0 in expressions (13) and (14): In expressions (26) and (27),ü 0 is the only unknown; in order to determine this value, one can assume that: After solving (28) foru 1 and substituting the resulting expression in (26), one has: The substitution of (29) in (27) produces the following expression: Bearing Eq. 30 in mind, Eq. 15 at the first time step (n = 0) can be particularized and written according to: 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: Particularizing Eqs. (32) and (33) for n = 0, one can arrive at the expression for u −1 : Note thatü 0 can be computed from Eq. 9, rewritten concisely below for t = 0: It is important to note that the values of p 0 , in Eq. 35, are directly computed from the analytical expressions of the initial condition u 0 . Finally,ü 0 is given by: 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 M a singular matrix and, consequently,ü 0 can not be computed from Eq. 36. If the analysis produces a non-singular matrix M, as occurs in the third example, Eq. 36 can be used.
For the use of the Houbolt method it becomes necessary to compute u 1 and u 2 ; from Eq. 32 one has: The computation of u 2 with the aid of Eq. 32 (by taking n = 1) requires the knowledge ofü 1 and, consequently, the knowledge of p 1 . This observation becomes clearer by writing an equation equivalent to Eq. 36: As 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: and Eq. 15 can be directly used for n ≥ 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.

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 β and γ . The choice of an adequate value of the parameter β 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.

One-dimensional bar with a time variable boundary condition
This example consists of a one-dimensional bar defined over 0 ≤ x ≤ a, fixed at x = 0 and subjected to a time variable sinusoidal boundary condition at x = a, i.e.
By assuming homogeneous initial conditions, the analytical solution to this problem is given by, see Stephenson [23]: where: In the analysis presented here, λ = π 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 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 β t = 1/3. The results for the potential at point C(a/2, a/4)  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 β and γ : for practical purposes, the best choice falls on β = 0.30 and on γ ≥ 0.52. This choice was conditioned by the accuracy in the results related to the flux at node B(0, a/4); here γ = 0.55 was adopted. Best results were achieved for β 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.  The second and the third analyses are carried out by assuming, respectively, initial conditions of the type: For the general case, the analytical solution is given by, see Stephenson [23]: 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   This example deals with a square membrane subjected to the initial conditions over the domain 0 ≤ x ≤ a, 0 ≤ y ≤ a, see Fig. 19:  The general analytical solution to this problem, for a rectangular membrane with dimensions a and b, according to Kreyszig [24], is:  where: 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  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 β = 0.30 and γ = 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 γ parameter was adopted.  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,   The analytical solution to this problem is given by, see Kreyszig [24]:  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

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