Solution of the two-dimensional scalar wave equation by the time-domain boundary element method: Lagrange truncation strategy in time integration

This work presents a time-truncation scheme, based on the Lagrange interpolation polynomial, for the solution of the two-dimensional scalar wave problem by the time-domain boundary element method. The aim is to reduce the number of stored matrices, due to the convolution integral performed from the initial time to the current time, and to keep a compromise between computational economy and efficiency and the numerical accuracy. In order to verify the accuracy of the proposed formulation, three examples are presented and discussed at the end of the article.


Introduction
The Boundary Element Method (BEM) has been applied to solve time-dependent problems quite successfully, as demonstrated by the several works, dealing with different approaches, published during the last years. For general purposes, these approaches can be classified according to the nature of the fundamental solution adopted. The use of time-dependent fundamental solution originates time-domain formulations (TD-BEM). TD-BEM formulations, beside providing good representation of causality and time response jumps and, consequently, leading to very accurate results, fulfill the radiation condition, which makes them suitable for infinite domain analysis, e.g. (Mansur 1983, Dominguez 1993. The use of static fundamental solution, on the other hand, originate two formulations, classified according to the maintenance or not of the inertial domain integral in the BEM equations as: D-BEM, that keeps the domain integral (D means domain) in the equations Telles 1992, Hatzigeorgiou andBeskos 2001), and DR-BEM (DR means Dual Reciprocity) that, by means of suitable interpolation functions, substitutes the domain integral by boundary integrals (Kontoni andBeskos 1993, Partridge et al. 1992). Formulations based on frequency and Laplace domains are also available, e.g., (Manolis 1983, de Lacerda et al. 1996. More recently, a formulation based on the Operational Quadrature Method appeared in the literature (Gaul andSchanz 1999, Schanz 2001).
This work is concerned with the solution of 2D scalar wave problem and is based on the TD-BEM formulation. The aim here is to reduce the storage computational cost due to the evaluation of the convolution integral that appear in the formulation. In other words, the aim is to reduce the number of assembled matrices, necessary to take into account the time history contribution, and to preserve the accuracy of the standard TD-BEM formulation. It seems that the best strategy to achieve this goal is to truncate the time integration, as previously presented (Demirel and Wang 1987, Mansur and de Lima-Silva 1992, Soares Jr. and Mansur 2004: the present work is concerned with this topic. Here, the whole time interval of analysis (0 :s; t :s; tn) is divided into two parts, both constituted by time steps, At, of equal size: the first one is restricted to the interval tk :s; t :s; tn, where tk is a specific value of time. In this first interval, the time integration is effectively done, i.e., the matrices are assembled in the standard way. In the second interval, defmed by 0 :s; t :s; tb discrete values of time are chosen and a Lagrange polynomial, passing through these discrete values of time, is constructed. Proceeding in this way, only matrices at these specific values of time are appropriately assembled, and matrices corresponding to intermediate values are computed by interpolation. The numerical computation of the convolution integral is carried out in this way. In order to simplifY the nomenclature, the first interval, constituted by the last niNr time intervals, will be referred to as integration interval; the second interval, by its tum, will be referred to as interpolation interval. Naturally, the length of the so-called integration interval and the degree of the polynomial in the interpolation interval are problem dependent parameters; incorrect choices of these parameters can lead to not reliable results.
Linear boundary elements and linear triangular cells were employed, respectively, to approximate the boundary and the part of the domain with non-homogeneous initial conditions (note that the solution of problems with non-homogeneous initial conditions does not present any difficulty: the domain integrals, related to the initial conditions are computed as in the standard TD-BEM formulation). Linear and constant time variation were assumed, respectively, for the potential and its normal derivative (flux) and, as usual in TD-BEM formulations, time integration was carried out analytically.
Three examples are presented and discussed at the end of the article, in order to verifY the applicability of the proposed interpolation scheme.

TD-BEM formulation
The TD-BEM equations can be written by employing the kernel regularization procedure (Mansur 1983) or the concept of fmite part of integrals (FPI) (Hadamard 1952). If, as it is usual in TD-BEM formulations, time integration is performed analytically, the resulting time integrated kernels from both representations are the same (Mansur and Carrer 1993), that is, both representations are entirely equivalent. The use of the Hadamard's concept, however, leads to more compact expressions and, it is the authors' opinion, provides an elegant representation of the equations involved in the analysis. A brief summary of the equations is given below.

Basic integral equation
Time-domain integral representation of the 2-D scalar wave propagation problem is written as (Carrer and Mansur 1996): In Eq. (1), r is the boundary; 0 is the domain, or the part of the domain, that presents nonhomogeneous initial conditions; the coefficient c( ~) assumes the same values of the static case, i.e., it is equal to 1 ( ~ E 0) or ( a/2~r) ( ~ E r and a is the internal angle depicted in Fig. 1 ); and the subscript 'o' means that r= 0.
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: where: • 2c U (X, t;~, r) = --;:::::::;.=:::::: In expression (2), H[c(t-r)r] stands for the Heaviside function (r = r(X; ~) is the distance between the field (X) and the source(~) points). where: and u*(x t·c; r) = 2 cr r ' ' ' The symbolry on the second term on the right-hand-side ofEq.

Space derivative boundary integral equation for internal points
The derivative of Eq. (1) with respect to a generic direction m(~), if~ E n, can be written as (Carrer and Mansur 1996): .lr 0 The FPI in the first term on the right-hand-side of Eq. (9) is interpreted as: The derivative with respect to r of the FPI indicated on the third term on the right-hand-side of Eq. (9) is defmed as follows: where:

Time derivative boundary integral equation for internal points
The derivative of Eq. (1) with respect to time, if~ E n, can be written as (Carrer and Mansur 1996): The FPI in the first term on the right-hand-side of Eq. (13) is interpreted as: (14) is given by: where: The time derivative of the FPI indicated on the second term on the right-hand-side of Eq. (13) is defmed as follows: The space and time derivatives of the initial conditions domain integrals were only indicated in Eqs. (9) and (13). For a detailed discussion concerning domain integration, the reader is referred to (Mansur 1983, Carrer and.

Numerical procedure
For the numerical solution of Eq. (1) the boundary is approximated by linear boundary elements, and the domain (or the part of it in which non-homogeneous initial conditions appear) is approximated by linear triangular cells. Time approximation assumes linear variation for the potential and constant variation for the potential normal derivative (flux), along the time steps At in which the overall time of analysis, 0 :s; t :s; tm is divided. Time integrals are computed analytically and the remaining kernels, in the boundary integrals, are computed numerically by the Gaussian quadrature.
The application of the discretized version of the boundary integral equation to all boundary nodes produces a system of algebraic equations, which can be written according to: In Eq. (18), diagonal matrix C is constituted by the c(~) coefficients, matrices Hnm and Gnm result from the spatial integration of the time-integrated kernels related, respectively, to u;(.x, t;~, r) and to u*{.X, t;~, r), and the vector Fn contains the initial conditions contributions (see Eq. (1)). Note, additionally, that the subscripts n and m stand for the time tn (fmal time) and tm (previous times), respectively. After imposing the boundary conditions, the system of equations represented by Eq. (18) can be solved for the boundary unknowns.
For a more detailed discussion concerning these matters, the reader is referred to (Mansur 1983, Dominguez 1993, Carrer and Mansur 1996. Note that equations similar to Eq. (18) are obtained for the solution of Eqs. (9) and (13) and that, in these equations, one has C = I.

Lagrange interpolation polynomial for time truncation
According to Eq. (1), or to its corresponding discretized version given by Eq. (18), it is necessary to take into account the contribution of the history, i.e., it is necessary to take into account the contribution of all responses previous to tn to obtain the response at time tn. If very large values of n are required, i.e., if late time results are required, and if there is no memory available in the computer, the use of time truncation procedures becomes desirable and justified. Naturally, such a procedure introduces approximations in overall time integration: these approximations depend on the length of the intervals in which the time integration is effectively done and the interpolation takes place. The question that arises is: What is the best strategy for an efficient time truncation? It is the authors' opinion that the response to this question is not conclusive, yet. Truncation procedures, for infmite domain applications, were first reported by Demirel and Wang (1987): in this work, it is assumed that there is a time, say t., before which the contributions to the response at a late time tm tn >> t., can be disregarded; in this way, the convolution integrals in Eqs. (1), (9) and (13) start at r= ts instead of at r= 0. Another truncation scheme, which applies quite well to bounded domains analyses, was presented by (Mansur and de Lima-Silva 1992) recently presented (Soares Jr. and Mansur 2004). In order to provide topics for discussion, this work presents a time truncation strategy based on the use of Lagrange interpolation polynomials. The basic idea can be outlined as follows: i) initially, the fmal time of analysis, say tn, is divided in time intervals At, such that tn = nAt; ii) the number of time intervals in which the time integrations will be carried out in the standard way, num is then defmed (or specified as a input datum in the computer program  (19), the interpolated matrices associated to matrices Hand G are denoted by Lh and Lg, respectively. It is important to mention that necessarily t 0 and tk must belong to the set of selected discrete time values. As a matter of fact, t 0 is always the first value and tk is always the last value of the interpolation interval; intermediate values can be equally spaced between them, or not. An illustration of the scheme is shown in Fig. 2.
A generic interpolated matrix, say B, assumed to be a function of time, can be represented as: lgr in which Bnm represent the matrices computed appropriately at the lgr selected discrete time values, for t = tn. The Lagrange interpolation polynomial can be defined according to: and has the property: It is important to point out that t 1 does not represent ~=}At but, instead of it, t 1 represents the j-th selected discrete time value and that, according to what was previously mentioned, t 0 = t 0 and tlgr = tk.
Another aspect to be mentioned is the generality of the proposed procedure: it can easily be applied to alternative TD-BEM formulations (Yu et al. 1998, Carrer andMansur 2002).

Examples
In the examples presented in this work, reference will be made to the dimensionless fJ parameter: in which C is the boundary element length. The choice of the value of the fJ parameter is a problem dependent task: as a general rule, small and large values, inside the interval 0 < fJ < 1, must be avoided.
The following notation will be employed in the examples: n represents the total number of time intervals and niNr represents the number of time intervals that constitute the integration interval; besides, equally spaced time values were adopted to construct the Lagrange polynomial in the interpolation interval. Along the discussion, the proposed formulation, for simplicity, will be referred to as Lagrange formulation.
Additionally, in the first and in the second example, E is the Young's modulus.

One-dimensional rod under compression (waveguide)
This example consists of a one-dimensional rod fixed at one extreme and free at the other, that is subjected to a compression load suddenly applied at t = 0 and kept constant during the analysis, see Fig. 3. The material is such that c = 1. The boundary discretization employed 24 elements, as depicted in Fig. 4. The number of time intervals is n = 320, and the time interval length was defmed for fJ= 0.6. Results furnished by the standard TD-BEM formulation, corresponding to the potential at node A(a, b/2) and to the flux at node B(O, b/2), are presented in Figs. 5 and 6, respectively, and were included to demonstrate how the proposed interpolation scheme can produce reliable results,  Figs. 7 and 8, respectively. Results for the potential, and its space and time derivatives at internal point C(a/2, b/2) are presented in Figs. 9, 10 and 11. Note that, from Fig. 5 to Fig. 11, the BEM results are always compared with the corresponding analytical solution. The oscillations around the analytical solution, in Figs. 8, 10 and 11, appear even in the results from the standard TD-BEM formulation, and better results can be obtained only with the use of a more refmed boundary element mesh. It is important to mention that several analyses have been performed by the authors, aiming at fmding a pattern for the adoption of the best length for the integration interval and the order of the interpolation polynomial. The conclusion is that the choice of these parameters depends on the problem; thus, the experience of the analyst plays an important role in this matter. To illustrate this, another analysis was included; see Fig. 12, in which results not so accurate (when compared to those from Fig. 7) were achieved for the potential at node A(a, b/2): although the integration interval was kept the same, niNr = 40, a poor interpolation polynomial of 7th order was adopted: the values t 1 are equally spaced but, now, one has j = 0, 40, 80, 120 ... , 280. Now, the number of matrices assembled and stored is equal to 2 (41 + 7) = 96, which means that only 15% of the total number of matrices are effectively stored.

Circular cavity in an infinite medium
This example consists of a circular cavity, in an infmite medium, subjected to an internal pressure suddenly applied at t = 0 and kept constant during the time, see their asymptotic values: in the results from the Lagrange formulation, the oscillation is more pronounced, but are still in agreement with the corresponding ones from the standard formulation.
In this example, the number of matrices required by the proposed scheme is given by 2(81 + 15) = 192; comparing this value with the number of matrices required by the standard TD-BEM formulation, given by (381 + 380) = 761, one can verify that the analysis was performed by employing only 25.2% of the matrices originally required.

Square membrane under prescribed initial velocity
The square membrane depicted in Fig. 19, with initial velocity field V 0 = c prescribed over the sub-domain 0, and with zero displacements prescribed all over the boundary, is analysed in this u=o a u==o a/5l~ u=o  shown in Fig. 20. As previously pointed out (Mansur 1983), the adoption of fJ< 0.6 does not introduce any great amount of noise in the BEM results. For this reason, the value fJ= 0.2 was adopted. Besides, as the number of discrete values increases, a better picture of the results can be inferred from the numerical results. Fig. 21 presents the results for the flux at boundary node A(a, a/2). Potential, space and time derivatives results, for internal point B(4a/5, a/2), are presented in Figs. 22, 23 and 24, respectively. This analysis was carried out by adopting n = 200, niNr = 60, and an interpolation polynomial of 14th order, with discrete time values, t 1 , equally spaced, i.e., j = 0, 10, 20, ... , 140. In this example, the number of matrices required by the proposed scheme is equal to 2 ( 61 + 14) = 15 0 whereas the number of matrices required by the standard TD-BEM formulation is (201 + 200) = 401; consequently, the analysis was performed by employing only 37.4% of the matrices originally required. Finally, one can conclude that the numerical results agree quite well with the analytical solution (Morse and Ingard 1968).

Conclusions
In this work, a strategy for the solution of the 2-D scalar wave propagation problem, by the TD-BEM formulation, is developed. The aim is to reduce the computational costs from the assemblage and from the storage of the matrices related to the time-history contributions to the results at a specific value of time. These matrices are due to the convolution integral presented in the basic TD-BEM formulation. Computational costs are reduced by partially computing the convolution integral, i.e., the time integration is no longer performed from t 0 to tn but, instead of it, from some value tk to tn; the interval [tk, tn] is designated integration interval, meaning that the matrices are appropriately computed there (BEM matrices). In the remaining interval, [t 0 , tk], designated interpolation interval, the matrices are computed by interpolation. To do so, a Lagrange polynomial is constructed by selecting discrete lgr values of time in the interval [t 0 , tk]: BEM matrices are computed for these lgr values and, fmally, matrices corresponding to values of time different from the selected ones are computed by interpolation. It is important to mention that the Lagrange interpolation formulation was also adopted for the computation of space and time derivatives of the potential at internal points and for the analysis of problems with non-homogeneous initial conditions. For the chosen examples, the numerical results can be considered good. As the experience plays an important role in the choice of the new parameters, i.e., the integration interval length and the order of the interpolation polynomial, in the absence of a study concerning error estimation, two practical recommendations are suggested by the authors: i) the use of interpolation polynomials of order greater than or equal to 10, with equally spaced time values; ii) the ratio between niNr and n must belong to the interval: 0.15 :s; niNrln :s; 0.30. Note that the proposed procedure can easily be extended to elastodynamics, as well.