An efficient stabilized boundary element formulation for 2D time-domain acoustics and elastodynamics

The present paper describes a procedure that improves efﬁciency, stability and reduces artiﬁcial energy dissipation of the standard time-domain direct boundary element method (BEM) for acoustics and elastodynamics. Basically, the developed procedure modiﬁes the boundary element convolution-related vector, being very easy to implement into existing codes. A stabilization parameter is introduced into the recent-in-time convolution operations and the operations related to the distant-in-time convolution contributions areapproximatedbymatrixinterpolations.Asitisshown in the numerical examples presented at the end of the text, the proposed formulation substantially reduces the BEM computational cost, as well as its numerical instabilities.

There are nowadays many publications that report robust techniques to improve TDBEM stability. In the approach developed by Frangi [9], space and time shape functions are not independent as in standard TDBEM approaches: a modified space-time approximation, which represents more accurately the causality of the phenomenon, is employed. Yu et al. [22,23] achieved more stable TDBEM algorithms by introducing, respectively, extra time and space (Galerkin approach) weighting integrals. Three other quite simple and effective schemes presently available are the linear θ method [14,20,21,24], the ε scheme and the half-step scheme [1,15]. Many other stabilization schemes can be found in the literature; some being more complex than those mentioned above (see [4,5]) or else, dedicated to minimize one specific source of instability [3,13].
The present paper describes a new procedure that improves stability and artificial energy dissipation of TDBEM for acoustic and elastodynamic modelling, besides reducing drastically CPU time and storage area. The developed procedure (α-δ method) utilizes: (1) an α stabilization parameter to modify the recent-in-time convolution operations, by replacing current values by weighted averages; (2) matrix interpolations to approximate the distant-in-time convolution contributions (the δ parameter controls the size of the recent/distant-intime regions). The algorithm hereby developed, as one will notice, is quite efficient and stable.
The α-δ TDBEM presented here associates stability and efficiency, and it is expected to be also more robust than previously reported TDBEM transmitting boundaries (as for instance in a BEM-FEM coupling context [16,17,24]).

Standard boundary element formulation
The time-domain integral representation for acoustic and elastodynamic wave propagation, in the absence of initial conditions and internal sources/forces, can be written as described by Eqs. (1a) and (1b), respectively where c(ξ ) and c ik (ξ ) are geometrical parameters; p * (X, t; ξ , τ ) represents the potential acoustic fundamental solution at field point X, corresponding to a unit impulsive source applied at point ξ at time τ ; u ik * (X, t; ξ , τ ) represents the displacement elastodynamic fundamental solution component in the k direction at field point X, corresponding to a unit impulsive force applied on the i direction at point ξ at time τ ; q * (X, t; ξ , τ ) and τ ik * (X, t; ξ , τ ) are the flux and traction related fundamental solutions, respectively. FP stands for the finite part of the integral. The boundary integral equations (1) can be treated numerically [2,10], resulting in a system of equations that can be expressed in matrix form. Adopting H n and G n as being the matrices computed at the current timestep n, related to the flux/traction fundamental solution and the potential/displacement fundamental solution, respectively, the system of equations (2) is obtained where, in the case of acoustics [Eq. (1a)], V H and V G entries are nodal values of the boundary potentials (p(X, τ )) and fluxes (q(X, τ )), respectively. On the other hand, if an elastodynamic problem is being considered [Eq, (1b)], V H and V G are the vectors in which the nodal boundary displacement (u k (X, τ )) and traction (τ k (X, τ )) components are assembled, respectively.
After introducing the boundary conditions in Eqs. (2), the following final expressions are obtained where, as usual in time-domain BEM, the entries of X n are unknown variables at boundary nodes at discrete time t n , while the entries of vector Y n are the according known nodal values. Matrices A and B are obtained from the combination of matrices C + H 1 and G 1 , taking into account, as well, the prescribed boundary conditions. R n is the vector related to the convolution process of the BEM; it represents the complete history up to t n−1 . The time stepping procedure indicated by Eq. (3) requires convolutions (l = 1, 2, . . . , n − 1) to be carried out for n = 1, 2, . . . , N, where the final time is t N = N t, t being the time step. Further details on the implementation of standard time-domain BEM algorithms can be found, for instance, in Mansur [10], Carrer and Mansur [2], Dominguez [7] etc.

Efficient stabilized boundary element formulation
In order to achieve an efficient stabilized boundary element formulation, the present work modifies the BEM convolution vector [Eq. (3b)], introducing a stabilization parameter into the recent-in-time convolution operations and a time-truncation procedure to compute the distant-in-time convolution contributions. No modification is introduced in the evaluation of the standard BEM matrices H n and G n , only their manipulation along the convolution process is modified. Thus, the proposed formulation is easy to implement into standard BEM existing codes.
The algorithm that describes the developed boundary element formulation is represented by Eqs. (4) below (which are analogous to the standard BEM Eqs. (3)). Details about the elaboration of the terms in Eqs. (4) are presented next, in the sub-sections that follow.
In Eq. (4), the entries of X n and Y n are, once more, unknown and known nodal values, respectively, at boundary nodes at discrete time t n . MatricesĀ andB are obtained from the combination of matrices C + H 1 + J(α, 1)H 2 and G 1 + J(α, 1)G 2 , taking into account the prescribed boundary conditions of the problem (α is a stabilization parameter).R n is the vector related to the recent-in-time convolution process (it represents the time history from t n−L+1 up to t n−1 , where L is a parameter which defines the truncation point of the BEM convolution process) andR n is a truncated BEM convolution vector. J(·) and I(·) are scalar functions; they define stabilization coefficients and interpolation parameters, respectively. δ i,j represents the Kronecker delta, i.e., Further details about matricesĀ andB and vectorR n are given in Sub-sect. 3.1; details about vectorR n are given in Sub-sect. 3.2.

Recent-in-time convolution contributions
For the recent-in-time [i.e., n−L+1 ≤ l ≤ n in Eq. (3b)] convolution operations, an α parameter, which linearly combines three consecutive time-step results, is introduced for the stabilization of the standard boundary element formulation. The combination of consecutives results has been used in several applications (adoption of a relaxation parameter in iterative processes, for instance [16,17]), showing good performance in what concerns the stabilization of an oscillatory tendency. The basic idea here is to smooth the convolution operations, trying to avoid any possible numerical instability. Several algorithms may be developed to smooth these operations based on the combination of consecutive time-step results; the present work adopts a simple linear combination of three time-consecutive terms, as described bȳ By replacing V l byV l in Eq. (3b), one may obtain Eq. (4b), as well as matricesĀ andB. Taking into account Eq. (5), the coefficient function J(·), presented in Eq. (4b), is given by As it will be illustrated in the Sect. 4, this algorithm is very effective, reducing considerably typical TDBEM numerical oscillations, giving as a result a stable algorithm. Not only this formulation completely eliminates or substantially reduces usual instabilities by adopting 0.5 ≤ α ≤ 1.0, but also it is extremely simple and easy to implement into standard BEM existing codes.

Distant-in-time convolution contributions
This section discussion is concerned with a time-truncation approach where operations related to the distant-in-time [i.e., l ≤ n − L in Eq. (3b)] contribution of the convolution process are effectuated approximately. This approximation consists of evaluating the matrices H n−l+1 and G n−l+1 by interpolation, based on a few m- where T L is the time limit after which the approximations take place). The integer m is an input parameter to the approximated convolution analysis and it indicates the number of key time-points (T k ) to be used in the interpolation procedure (within the interlude of these key time-points, the interpolations will occur).
In order to obtain matrices H n−l+1 and G n−l+1 as function of m previously selected matrices H k and G k , an interpolation parameter I(·) is introduced into Eq. (3b), originating Eq. (4c). There are several interpolation procedures that might be chosen to specify I(·); for a detailed discussion, the reader is referred to Soares and Mansur [18]. The present work adopts a multilinear interpolation algorithm; i.e., linear interpolation between the m interpolation time-points is adopted. Taking into account the multi-linear interpolation procedure, the I(·) parameter is given by where H(·) is the Heaviside function. In the present work, the discrete times T k , which are the key-times to the interpolation approximation, are adopted as and the time limit T L is adopted as where δ is a precision control parameter;r is the average distance between boundary nodes and c is the (primary) wave velocity. The calculus of the time limit T L [Eq. (9)] is based on the behaviour of the function f (t) = 1/t, which is a simple function, similar in most aspects to the time related kernels of the models being considered.
As it can be seen in Eq. (8), the values of the discrete times T k are not necessarily equally spaced once the input parameter η is introduced. This scheme is adopted since better approximations are obtained if more interpolation points are chosen in the neighbourhood of T L , region in which kernels gradient and amplitude are more significant. The discrete times T k become equally spaced if η = 0.
It is important to observe that the kernels that are here being interpolated decrease monotonically with time increase; thus, linearly interpolated values are always greater than real ones. This property is beneficial: as it is well known, standard TDBEM algorithms usually introduce some spurious numerical damping; as a consequence, the multi-linear interpolation approximation acts as a subtle dispersion-correction algorithm (besides drastically reducing the BEM computational cost).

Numerical applications
Time-step length plays an important role in the timedomain BEM analyses [10]. A measure of the time-step length is computed according to the following expression where c is the (primary) wave velocity and is the boundary element length. In the present section, two numerical applications are considered. The first application (finite-domain problem) is concerned with an acoustic rod; the discussion presents numerical results obtained taking into account different discretization (β) and stabilization (α) parameters. In the second application (infinite-domain problem), an elastodynamic problem, composed by two close circular cavities, is focused. In both applications, typical numerical instabilities occur for standard TDBEM analyses, according to the parameter β adopted. In fact, taking into account numerical aspects, these problems are quite sensitive, due to the multiple wave reflections that occur.
As one can observe in the following applications, the proposed methodology deals quite well with the usual time-domain BEM limitations, namely: numerical instabilities and excessive oscillation, energy dispersion (numerical damping) and high computational cost.

Acoustic rod
This example [10] consists of an acoustic rectangular domain, submitted to a time Heaviside type flux applied at one of its end and prescribed null potentials applied at the opposite one. Figure 1a shows the boundary conditions (f (t) = H(t)) and geometry (a = 4.0 m and b = 2.0 m) of the body. Figure 1b shows the boundary discretization: the symmetry of the problem was considered and only 24 linear boundary elements of equal size ( = 0.25 m) were employed (an interesting feature of the boundary element formulation is that symmetric bodies under symmetric loads can be analysed without discretization of the symmetry axes; this can be accomplished by an automatic condensation process, which integrates over reflected elements and performs the assemblage of the final coefficients in a reduced matrix [19]). The solution of this 2D plane wave problem is in fact independent of the coordinate in the vertical direction, i.e., the plane wave front propagates in the longitudinal direction. Thus, its solution is the same as that of longitudinal waves propagating along a rod.
The wave propagation velocity along the rod is c = 2, 500 m/s. Different time discretizations have been adopted for the analyses. The time history results concerning potentials at points A and B and fluxes at point C are depicted in Figs. 2,3,4 and 5, for different discretization (namely: β = 0.2, β = 0.3, β = 0.4 and β = 0.5) and stabilization (namely: α = 0.0, α = 0.3, α = 0.5 and α = 0.7) parameters. As one may observe, the proposed formulation is quite stable (one should note that the standard BEM is described by α = 0.0).
In the above-described analyses, δ = 0.00 was adopted, i.e., no truncation on the convolution process was considered. Taking into account the approximations presented in Sub-sect.  were adopted), some time history results are presented in Fig. 6, for α = 0.5 and for four different β-values, namely:β = 0.25, β = 0.5, β = 0.75 and β = 1.00. As one may notice, the results are once more stable and no perceptible difference, concerning accuracy, may be observed among the truncated (δ = 0.04) and the complete (δ = 0.00) convolution analyses. However, concerning efficiency aspects, the difference among the analyses is quite outstanding. The relations among the truncated and the complete time convolution analyses are described in Table 1, considering two computational-cost aspects: total CPU time and storage area. As one can observe, the developed formulation may easily reduce more than 80% the time-domain BEM computational cost (one should observe that in Table 1 the relation "modified approach cost/standard approach cost" is presented, in percentage). In Fig. 7, another interesting feature of the presented formulation is depicted. As it has been mentioned, the multi-linear interpolation approximation acts as a subtle dispersion-correction algorithm, once standard timedomain BEM analyses are usually associated with some numerical damping. By increasing the error parameter (δ = 0.06), which controls the time limit (T L ) of the distant-in-time convolution contributions, the kernels are over-estimated and perceptible differences among the truncated and complete convolution results may be observed. However, if not excessive, this error is beneficial: it corrects the BEM numerical damping tendency. A clear exemplification of the above aspect is depicted in Fig. 7, where the results of the truncated analysis may be considered more accurate than the ones related to the complete convolution analysis.

Elastodynamic cavities
This example consists of two circular cavities located close to each other in an unbounded elastodynamic medium. One cavity is submitted to a time Heaviside type pressure-load, whereas the other remains unloaded. Figure 8a shows the boundary conditions (f (t) = H(t)) and geometric configuration (R = 1.0 m) of the problem. Figure 8b shows the cavities boundary discretization: 16 linear boundary elements of equal size ( = 0.3902 m) were adopted to model each cavity. The time history results concerning normal displacements at points A and B are presented in Fig. 9, for different discretization (β = 0.1376 and β = 0.2001) and stabilization (α = 0.0, α = 0.5, α = 0.75 and α = 1.00) parameters. As one may once again observe, the pro-posed formulation considerably reduces or completely removes the usual BEM numerical instabilities.
Taking into account interpolation approximations for the distant-in-time convolution contributions (δ = 0.03, m = 10 and η = 0.5 were adopted), some time history results are presented in Fig. 10, for α = 0.75 and for β = 0.1376 and β = 0.2001. The computational cost relations between the truncated (δ = 0.03) and the complete (δ = 0.00) time convolution analyses are described in Table 1. As it may be observed, the proposed formu-

Conclusions
Presently it is possible to find in the specialized literature many different schemes, which can either eliminate, or else reduce to quite acceptable levels the instability observed in many TDBEM applications concerning acous-tics and elastodynamics. Here, a new scheme has been presented, the α-δ scheme, which completely eliminates or substantially reduces instability, decreases artificial (numerical) energy dissipation and improves computational efficiency. The procedure can be employed either for 2D or 3D applications; in the former, the described compression scheme can be used to drastically reduce CPU time and storage area.
Two examples were presented where a parametric study concerning the α, δ and β analysis parameters was carried out, the conclusions being that quite a wide   interval of variation of these parameters lead to stable numerical results; many thousands of time steps could be used even in the most severe case (first example), which is prone to numerical instabilities. In fact, the proposed methodology efficiently managed to overcome the main drawbacks of TDBEM formulations: the introduction of the α parameter is extremely efficient dealing with the TDBEM instability problem and its awkward dependency on the space-time discretization parameter β; moreover, the introduction of the δ parameter overcomes the high computational cost of the time convolution process, introducing controlled errors (it can be as small as the user wants), which are beneficial to correct the numerical damping tendency of the TDBEM formulation. The present work smoothed the recent-time convolution operations by a simple weighted linear combination of three time-consecutive elements; however, many other equally simple derived alternatives are possible. The authors suggest as future work the development of such simple alternative schemes, and a comparison with published algorithms (see Sect. 1) in order to identify the optimum (one or more). Another point to be investigated concerns coupling of more than one algorithm, specially time weighting and Galerking with the present procedure, linear θ , ε scheme, half-step etc.