A step-by-step approach in the time-domain BEM formulation for the scalar wave equation

. This article is concerned with the presentation of a time-domain BEM approach applied to the solution of the scalar wave equation for 2D problems. The basic idea is quite simple: the basic variables of the problem at time tn (potential and flux) are computed with the results related to the potential and to its time derivative at time tn- 1 playing the role of "initial conditions". This time-marching scheme needs the computation of the potential and its time derivative at all boundary nodes and internal points, as well as the entire discretization of the domain. The convolution integrals of the standard time domain BEM formulation, however, are not computed; the matrices assembled, only at the initial time interval, are those related to the potential, flux and to the potential time derivative. Two examples are presented and discussed at the end of the article, in order to verify the accuracy and potentialities of the proposed formulation.


Introduction
Time-domain analysis by the boundary element method (BEM) is a very attractive and challenging area of research; as a consequence a lot of works have been presented so far, enriching the range of applications of the BEM and the literature concerning numerical methods in engineering. For general purposes, the different BEM approaches can be classified into two main groups: the first group, that employs the static fundamental solutions, and the second group, that employs time-dependent fundamental solutions; transformed-domain (frequency and Laplace) approaches can also be included in this group. A brief discussion concerned with different BEM approaches is carried out in what follows. For a more complete discussion concerning dynamic analysis by the BEM, the reader is referred to (Beskos 1997(Beskos , 2003. equation (potential integral equation), thus generating the velocity integral equation. As the potential and the velocity at time tn--1 play the role of "initial conditions" for the potential, for the flux and for the velocity at time tm the domain integrals, that in the TD-BEM formulation are restricted only to the part of the domain with non-homogeneous initial conditions, are now extended to the entire domain. Besides, by following an approach similar to that from the D-BEM formulation, e.g. Mansur 2004, Hatzigeorgiou andBeskos 2002), potential and velocity integral equations are applied at boundary nodes and internal points simultaneously, thus generating an enlarged system of equations. After solving this system of equations one has the values of the potential and of the flux at boundary nodes and the values of the potential at internal points. The values of the velocity are computed directly, since the boundary variables are already known. One question can be formulated at this moment: can this approach be used for infinite domain analysis? If results at late times are required, the answer is no: in this case the interruption of the domain discretization generates a fictitious boundary and, consequently, the influence of waves reflected there tends to invalidate the results. Although the response to the question is negative, it is the authors' opinion that the present study demonstrates the immense possibilities of research in the TD-BEM. Also, the two examples presented at end of this work demonstrate the potentialities of the proposed formulation and encourage future research in this area.

Standard TD-BEM integral equations
In what follows, a brief summary of the TD-BEM formulation is presented. For a more complete discussion concerning this topic, the reader is referred to (Mansur 1983, Dominguez 1993).

1 Potential integral equation
The potential integral equation associated to the TD-BEM formulation of the scalar wave problem, by making use of the concept of fmite part of integrals, FPI, (Hadamard 1952) can be written as Mansur 1996, 2006) (1) The following notation was employed in Eq. (1 ): [' is the boundary; n is the domain, or the part of the domain, that presents non-homogeneous initial conditions uo(X) and/or vo(X) = O'u(X, t) I ; the Ot t= 0 variable p(X, t) = du(X, t)ldn(X), where n(X) is the unit outward normal to the boundary at X, is the flux; and the coefficient c( ~ on the left-hand-side of Eq. (1) assumes the same values of the static case, e.g. (Mansur 1983).
The fundamental solution, u * (X, t; ~' r) , that corresponds to the effect of a source represented by an impulse at t = r located at X= ,; propagating with velocity equal to c, has the following expression The symbol of on the second term on the right-hand-side of Eq. (1) stands for the FPI, defmed by Hadamard (1952) 4 1 +u; (x, t;,;, r)u(X, r)dr = lim { rru: (x, t;,;, r)u(X, r)dr-lu*(x, t;,;, r)u(X, r)} (5) Jo r--+t-rlc Jo c

2 Velocity integral equation
The integral equation associated to the potential time derivative can be written as Mansur 1996, 2006) 4Jrc(,;)O'u~ t) = 1 1:+ u; (x, t;,;, r)p(X, r)drdf'(X) -1! (fa+ u;(x, t;,;, r)u(X, r)dr)! df'(X) _l r ar au;cx, t;,;) u (X)df'(X) + a( 1 + 1 ) where, for simplicity and because the resulting expressions do not require any special treatment in the domain integration, the time derivative of the domain integrals related to the non-homogeneous initial conditions are indicated by the last term on the right-hand-side of Eq. (6). The FPI in the first term on the right-hand-side of Eq. (6) is interpreted as The function u 1 (X, t; ,;, r) in expression (7) is given by The time derivative of the FPI indicated on the second term on the right-hand-side of Eq. (6) is defined as

Numerical solution
For the numerical solution of the problem, linear and constant time variations were assumed, respectively, for the potential and for the flux. Following the usual time-domain BEM approach, e.g. (Mansur 1983, Dominguez 1993, time integrals were computed analytically. It is important to mention that time integrals can also be computed numerically; in this case, one should employ the Kutt quadrature formulae (Kutt 1975) to evaluate the fmite part integrals that appear in Eq. (1) and in Eq. (6). For the boundary approximation, isoparametric linear elements are employed; the domain is approximated by triangular linear cells, in which the initial conditions vary linearly, e.g. .
After space and time discretization has been accomplished, the application of Eq.
(1) to all boundary nodes and internal points simultaneously (in the next section this procedure becomes clearer) generates the system of equations below In Eq. (10) the subscript n stands for the time tn and the subscript m, in the summation symbol, stands for the times tm, m = 0, 1, 2, ... , (n-1), previous to tn; the superscripts band d correspond to the boundary and to the domain (internal points), respectively. In the sub-matrices, the first and second superscripts correspond, respectively, to the positions of source and field points. In the first matrix on the left-hand-side of Eq. (10), the C matrix is related to the coefficients c( ~) of the boundary nodes whereas the identity matrix, I = H'ld, is related to the coefficients c( ~) = 1 of the internal points. Alternatively, Eq. (10) can be written in a more compact way as follows After the solution of the system of equations represented by Eq. (11 ), the velocity v(~, t) = O'u(~, t)!Ot at boundary nodes and internal points is computed by employing Eq. (6); in matrix form, one has a system of equations, similar to that represented by Eq. (10), given below In compact form, Eq. (12) can be written as follows Since matrix C is a diagonal matrix, Eq. (13) can easily be solved for Vw Note that the over bared matrices in Eqs. (12), (13) correspond to the matrices that appear in Eqs. (10-11).

Step-by-step procedure
The solution procedure employed in this article is based mainly in the time translation property, e.g. (Mansur 1983), which can be written as (14) where t 1 is any time increment.
(1) and (6), enables one to solve the time-dependent problem by considering the results related to u,_ 1 and v,_ 1 as "initial conditions" in the computation of the results at time tn-That is, in the computation of Un and Vm the time integration is restricted only to the interval [0, ~t]. For the approach proposed here, the following versions of Eq. (11) and Eq. (13) can be written (16)   (17) The values of the potentials and of the velocities, treated as "initial conditions", are stored in the vectors Un-I and v,_ 1 • Note that proceeding this way, the entire domain needs to be discretized. This feature generates the main limitation of the proposed approach, i.e., the proposed approach is suitable for fmite domain applications; infinite domains analyses are feasible only if early time responses are requested.

Examples
Two examples are presented next. For each example three analyses were carried out, aiming at verifying the accuracy and convergence of the numerical results to the available analytical solutions (Stephenson 1970). Comparison is also made between the results obtained from the present formulation with those obtained from the standard TD-BEM formulation. From the starting mesh of the first analysis, constituted by 48 boundary elements and 256 internal cells, see Fig. 1, subsequent meshes were obtained by duplicating the number of boundary elements and quadrupling the number of internal cells. In this way, the second analysis was carried out with 96 boundary elements and 1024 cells; the third analysis, with 192 elements and 4096 cells. As usual in time-domain BEM analyses, the dimensionless parameter fJ = c/)..t/1, which relates the wave velocity (c), the smallest element length (/), and the time interval /).t (assumed constant), is used to estimate the time-step length. It is important to mention that, for the proposed approach, the optimum time-step length is obtained by taking fJ = 0.5. One must observe that this is not the value recommended by (Mansur 1983) for the standard TD-BEM formulation, i.e., fJ = 0.6.

1 One-dimensional rod under a Heaviside-type forcing function
This is the classical example of a one-dimensional rod under a Heaviside-type forcing function, applied instantaneously at t = 0 and kept constant in time.
Boundary conditions and geometry are depicted in Fig. 2. The analytical solution to this problem can be found in (Stephenson 1970) and is given by u(x,t) =fa[:!+ 82f (-1f 2sin((2n-1)nx)cos((2n-1)Jrct)] (18) E a 7r n = 1 (2n _ 1)  1) and (6) (represented under the summation symbol in Eqs. (11) and (13)), are stored to compute the time-history contributions to the results at the current time, tn; for this reason, the standard TD-BEM analyses presented in this work were carried out by employing mesh 1. An alternative to avoid the storage of a large number of matrices could be the use of a truncation strategy, such as that presented by Carrer and Mansur (2006). Note that a small amount of numerical damping is also observed in the standard TD-BEM results. In order to verify the stability, the results related to the potential at boundary node A(a, b/2) are compared with the analytical solution, Eq. (18), for a larger observation time, see Fig. (6).

2 One-dimensional rod under sinusoidal initial conditions
In this example, the same bar of the previous example, now fixed at both ends, is analysed when subjected to initial conditions of the type In analysis (a) The analytical solution to this problem can be found in (Stephenson 1970) u(x, t) = b 0 cos(tret!a)sin(nx/a)

In analysis (b) one has
The analytical solution to this problem can be found in (Stephenson 1970) u The three meshes adopted to perform analyses (a) and (b) are the same of the previous example. By considering, initially, analysis (a), results related to the potential and to velocity at point C of coordinates (a/2, b/2) are depicted, respectively, in Figs. 7 and 8. Results related to the flux at point A(a, b/2) are depicted in Fig. 9. Fig. 10   Standard TD-BEM analyses were carried out with fJ = 0.62 and the corresponding results were included in Figs. 7-9 and 11-13. In this example, the results related to the potential at the boundary node C(a/2, b/2) are compared with the analytical solutions, given by Eqs. (19) and (20), in Figs. 15 and 16, respectively, for a larger observation time, in order to verify the stability.
One can conclude from Figs. 6, 15 and 16 that the results from the present formulation are quite stable; however, the reduction of the numerical damping is a task that must deserve special attention.

Conclusions
In this work, a modified version of the standard TD-BEM formulation is presented. The main characteristic of the proposed formulation and, perhaps, what could turn it very attractive, is that the convolution integrals are no longer evaluated (as a matter of fact, time integration is restricted to only one interval). Now, the results at a given value of the time, say tm are computed with the results related to the potential and to its time derivative at the time tn-1 assumed as "initial conditions". Note that such a version of the TD-BEM formulation, in addition to the basic integral equation, called here potential integral equation, also makes use of the integral equation corresponding to the velocity. Both integral equations are applied simultaneously to all boundary nodes and internal points, in order to obtain the results that will be used as "initial conditions" for the next time-step. This is the main disadvantage of the proposed formulation: as the entire domain needs to be discretized in order to compute the contribution of the "initial conditions" to the results at the time tm the analysis of infmite domain problems becomes almost impossible or, at least, very expensive in the computational point of view, if late time results are required. It is important to point out that, although the numerical results are in good agreement with the analytical ones, special attention must be devoted to the reduction of the numerical damping observed in the analyses presented here: if this problem is overtaken, the procedure proposed in this work can easily be extended to elastodynamics.