Alternative time-marching schemes for elastodynamic analysis with the domain boundary element method formulation

This work presents alternative time-marching schemes for performing elastodynamic analysis by the Boundary Element Method. The use of the static fundamental solution and the maintenance of the domain integral associated to the accelerations characterize the formulation employed in this work. It is called D-BEM, D meaning domain. Time response is obtained by employing step-by-step time-marching procedures similar to those adopted in the Finite Element Method. Among all integration procedures, Houbolt scheme became the most popular used to march in time with D-BEM formulation, in spite of the presence of a high numerical damping. In order to improve the integration, this work presents alternative schemes that can be used to perform elastodynamic analysis by the BEM with a better damping control. In order to verify the accuracy of the proposed scheme, three examples are presented and discussed at the end of this work.


Introduction
Transient elastodynamic analysis by the BEM is an area of research that offers a wide range of possibilities for the solution of the problem. This means that various methodologies are available and a lot of research work can be done.
Time-domain BEM formulations (TD-BEM), e.g., Mansur (1983), Dominguez (1993), Mansur et al. (1998), provide good representation of causality and time response jumps and lead to very accurate results. Besides, the fulfillment of the radiation condition makes them suitable for infinite domain analysis. Formulations that use the static fundamental solution can be 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, e.g., Carrer and Telles (1994) and Hatzigeorgiou and Beskos (2001), and DR-BEM (DR means Dual Reciprocity), that by means of suitable interpolation functions substitutes the domain integral by boundary integrals, e.g., Kontoni and Beskos (1993) and Partridge et al. (1992). Formulations based on frequency and Laplace domains are also available, e.g., Manolis (1983) and de Lacerda et al. (1996). More recently, a formulation based on the Operational Quadrature Method appeared in the literature, e.g., Gaul and Schanz (1999) and Schanz (2001).
In this article, attention will be devoted to the twodimensional D-BEM formulation. It is expected that some conclusions found here will also be valid for the DR-BEM formulation. The discussion carried out here is concerned with the use of alternative time-marching schemes for the D-BEM formulation. The Houbolt scheme, e.g., Houbolt (1974), Bathe (1996), and Weaver and Johnston (1983), until now has been the most popular algorithm for performing transient analysis by the D-and DR-BEM formulations. Probably, the damping properties, presented by the Houbolt scheme in FEM analysis, are benefic in BEM analysis, justifying the wide range of applications that use this technique. The search for other schemes, however, is an important task and some work related to the use of alternative time-marching solution procedures is discussed here. Two schemes are proposed in this work: initially, the BEM version of the Hilbert-a procedure, Hilbert et al. (1977), is presented; afterwards, according to the time-marching procedure adopted, one has the schemes named Houbolt-a and Newmark-a. In order to demonstrate their efficiency and accuracy, three examples are presented and discussed at the end of the article.
It is important to notice that the D-BEM (or the DR-BEM) approach presented here follows a strategy similar to that of the FEM Galerkin formulation, i.e., a time independent weighting function (the static Green function in the case developed here) is employed to establish the basic weighted residual sentence. Thus the inertial domain integral remains and step-by-step time marching algorithms, similar to those employed in the FEM, apply. Advantages and drawbacks of methods that employ domain discretization arise in the D-BEM and DR-BEM approaches: they are quite adequate to consider non-linear effects, e.g., Carrer and Telles (1994), Hatzigeorgiou and Beskos (2001), Kontoni and Beskos (1993), etc., but are not appropriate for infinite domain analyses. On the other hand, the TD-BEM approach is quite suitable for infinite domain analyses, as the radiation condition is automatically fulfilled, and is quite robust for linear problems. The best strategy for those engaged in doing dynamic linear and non-linear analyses of infinite domain wave propagation problems is to have a computer code where both formulations are available. Such a code should include routines containing the coupling of TD-BEM and D-BEM, so that the former can work as an efficient transmitting boundary for the later: in this way one will have an efficient and robust code for dynamic linear and non-linear bounded and unbounded domains.

D-BEM Formulation and time-marching schemes
Before presenting the BEM formulation, it is important to mention that the Sommerfeld radiation condition is not fulfilled by the integral equation obtained using the static fundamental solution as weighting function and, as a consequence, infinite and semi-infinite domains cannot be treated by boundary discretizations only. However, in spite of this unfavorable aspect, the use of the static fundamental solution leads to a very simple approach. Besides, this formulation presents some advantage in elastoplastic analysis, since the part of the domain where inelastic variables are expected to occur needs to be discretized, see for instance Telles (1983), Carrer and Telles (1992), and Beskos (1995. Domain discretization can be avoided if an alternative approach, known by dual reciprocity (DR-BEM), is employed (the reader is referred to Partridge et al. (1992) for additional details on the DR-BEM).
The starting D-BEM integral equation is written as follows: In Eq. (1) G is the boundary and W is the domain of the problem, X and n are the field and source points, and q is the constant mass density. The fundamental solution u Ã ik n; X ð Þcorresponds to a displacement at the field point X in k direction associated to a unity load applied at the source point n in i direction. The same interpretation is valid for the traction p Ã ik ðn; XÞ In plane strain problems, one has: where r = r(n, X) is the distance between X and n and d ik is the Kronecker delta. In expression (2) G is the shear modulus and m is the Poisson coefficient (note that in the section Numerical Applications reference will be made to the parameter E = G/2(1+m) (Young modulus)); in expression (3) n stands for the direction of the outward normal vector (n k is its component in the k direction). Additionally: Plane strain expressions are valid for plane stress if the Poisson coefficient m is replaced by " is obtained by geometric considerations. The reader is referred to Hartmann (1980) for additional details concerning its direct computation. Note that for internal points one has C ik (n) = d ik .
In the present work, linear elements were employed in the boundary discretization, and linear triangular cells were employed in the domain discretization. In order to solve the problem, the acceleration in Eq. (1) is approximated by recurrence formulae involving displacements at different time levels (including the displacement at the present time). As a consequence, the unknowns to be determined at a given time instant, say t n+1 , are not only the displacements and the tractions on the boundary, but also the displacements at the internal points. The application of the discretized version of Eq. (1) to all boundary nodes and internal points produces an enlarged system of equations, schematically written as (see Carrer and Telles (1994)): In Eq. (5) the superscripts b and d correspond to the boundary 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 corresponds to the position of the field point. The identity matrix I is related to the coefficients C ik = d ik of the internal points.
Equation (5) can be written in a more compact way as: The time integration scheme proposed in this work was inspired in the Hilbert-a method, e.g., Weaver and Jonhston (1987), and Hilber et al. (1977), and starts writing a D-BEM version of this scheme as follows: Subsequently, the well known algorithms by Houbolt (1974) and Newmark (1959) are alternatively used in order to approximate the accelerations. The Houbolt scheme, Houbolt (1974), adopts the following approximations for the velocity and acceleration: In the Newmark scheme, Newmark (1959), one has: The stability and accuracy of the Newmark scheme depend on the correct choice of the parameters b and c. If one takes b = 1/6 and c = 1/2, expressions (10) and (11) correspond to the linear acceleration method. Stability, however, is not guaranteed; as a matter of fact, only unreliable results (not presented here) were obtained when b = 1/6 is adopted. If b = 1/4 and c = 1/2, expressions (10) and (11) correspond to the average acceleration method and the scheme is unconditionally stable, Weaver and Jonhston (1987). These are the values adopted in this work. Substituting the approximation for the acceleration given by Eq. (9) in Eq. (7), one has the scheme designated here as Houbolt-a; analogously, the Newmark-a scheme is obtained after substituting Eq. (11) in Eq. (7). After performing all algebraic operations, a general expression can be written as follows (note that vector g n on the righthand-side of the system of equations contains the contribution corresponding to the previous instants of time): It is important to point out that an adequate choice of the time-step plays a fundamental role in the analysis. A dimensionless variable was adopted in order to provide a measure of the time-step Dt (see Mansur (1983)). It is defined here as: where l is the smallest boundary element length and c d is the primary wave propagation velocity: A correct choice for the value of the parameter b Dt is conditioned by the experience of the researcher. Diversely from the FEM formulation, where the definition of the critical time-steps can be done appropriately, e.g., Bathe (1996) and Weaver and Jonhston (1987), in BEM formulations, as can be observed, for instance, in the works by Mansur (1983), Partridge et al. (1992), Kontoni and Beskos (1993) and Hatzigeorgiou and Beskos (2001), there is no specific studies concerning this matter; for this reason, the experience plays a fundamental role in the adoption of the time-step. The simplest and straightforward recommendation states that values of b Dt greater than unity must be avoided.
As mentioned previously at the Introduction to this work, the damping properties, presented by the Houbolt scheme in FEM analyses, are benefic in BEM analyses. The use of the Newmark scheme, on the other hand, tends to produce amplification in the numerical results, leading to unreliable results. The introduction of the a parameter, as demonstrated in the examples included in this work, provided a better control of the damping and rendered possible the analyses with the Newmark scheme. It is verified that a > 0 increases the damping and that a < 0 reduces the damping. Consequently, Newmark-a and Houbolt-a schemes must employ, respectively, positive and negative values of a.  Fig. 2. In order to simulate the one-dimensional problem, the Poisson coefficient is taken as null. For the other material parameters, the following values were adopted: E = 100. 0 and q = 1.0.
It was verified empirically that negative values of the a parameter reduce the damping and that positive values Fig. 1. One dimensional rod: geometry and loading definition tend to increase the damping. As the Houbolt scheme presents a high numerical damping, it seems that its use with negative values of a is the most appropriate combination. However, care should be taken in order to avoid amplification of the numerical solution. The correctness of this conclusion can be verified in Figs. 3 and 4 in which the results obtained with a = )0.05 and corresponding to the displacement component u x at boundary node A(a, a/4) and to the traction component p x at boundary node B(0, a/4) are compared with the corresponding analytical solutions and with the solutions corresponding to a = 0. In this analysis, it was adopted b Dt ¼ 1=3. It can be observed that the Houbolt-a results related to the displacement component u x are better than those of the standard Houbolt (Fig. 3); on the other hand, the Houbolt-a produces a more pronounced jump than the standard Houbolt for the traction component p x (Fig. 4). In order to demonstrate the correctness of the choice of the parameters a and b Dt , two additional analyses were carried out for different values of the b Dt parameter: the first analysis employed b Dt ¼ 1=6 and, for the second analysis, b Dt ¼ 1=2 was adopted. Due to the difficulty in the representation of the jumps that occur along the time, only the results related to the traction component p x will be discussed next (naturally, if the results related to p x can be considered acceptable, the same conclusion will be valid to the results related to u x ). The results depicted in Fig. 5 present a large level of oscillation and can not be considered acceptable, demonstrating that b Dt ¼ 1=6 is not a good choice. The same conclusion is valid for the use of b Dt ¼ 1=2: the results presented in Fig. 6 are less accurate that those presented in Fig. 4.
Newmark scheme, on the other hand, tends to become unstable in the D-BEM formulation, producing unreliable results (see Fig. 8). For this reason, and bearing in mind that positive values of the a parameter produce an increasing of the damping, the combination of the Newmark scheme with positive values of a (a = 0.10) was successfully accomplished, as can be verified from Fig. 7 (u x at A(a, a/4)) and Fig. 8 (p x at B(0, a/4)). The results related to u x are in better agreement with the analytical solution than those furnished by the standard Newmark scheme, as shown in Fig. 7. Besides, the introduction of the damping was responsible for the reliable results     Fig. 8 and related to p x . In this analysis, b Dt ¼ 2=3 was adopted, because greater values of the timestep tend to produce more damping. The same discussion concerning the choice of the b Dt parameter was carried out again. Figs. 9 and 10 present the results related to the traction component p x obtained with b Dt ¼ 1=3 and b Dt ¼ 1, respectively. Keeping in mind the criterion of a better representation of the jumps, the results in Fig. 9 are less accurate than those in Fig. 8, showing that small values of Dt are not adequate when employing the Newmark-a scheme. Although from the results in Fig. 10 a good picture of the analytical response can be inferred, these results present a more pronounced jump than those presented in Fig. 8, demonstrating that the best choice is b Dt ¼ 2=3. In Figs. 3 and 7, u st is the horizontal static displacement; u st = )0.12) and, in Figs. 4, 5, 6 and 8, 9, 10 p st is the static reaction horizontal component.

Deep beam simply supported
This plane stress example consists of a simply supported deep beam, depicted in Fig. 11, submitted to a suddenly applied uniform load kept constant in time. FEM and BEM analyses were carried out in order to verify the accuracy of the results furnished by the latter. FEM analysis employed linear quadrilateral finite elements and the Newmark integration scheme with Dt = 0.10. FEM mesh is depicted in Fig. 12. BEM mesh is the same already depicted in Fig. 2. Only half of the beam was discretized, and the essential boundary conditions were taken into account as follows: the displacements on the left side were restricted in the vertical direction; on the right side the displacements were restricted in the horizontal direction. Proceeding this way, the symmetry of the problem was simulated appropriately. The material parameters are: E = 100.0, m ¼ 1=3 and q = 1.5. The results corresponding to the vertical displacement component u y at boundary node A(0,0), obtained with b Dt ¼ 1=3, are presented in Figs. 13-14. In this example, the results furnished by the Houbolt scheme (standard version or with a = )0.05) are quite stable, presenting a good agreement with the FEM results. The results furnished by the Newmark-a scheme with a = 0.20, on the other hand, are remarkably better than those furnished by the standard Newmark scheme, as can be seen in Fig. 14; besides, they present a good agreement with the FEM ones. In Figs. 13-14, u st is the static vertical displacement, computed according to the classical beam theory; u st = )0.024.

Circular cavity under internal pressure
This example consists of a circular cavity of radius R subjected to an internal pressure suddenly applied and kept constant in time, as depicted in Fig. 15. Boundary and domain discretizations are depicted in Fig. 16 (note that the symmetry of the problem was automatically considered by the integration of reflected elements and cells, see Brebbia et al. (1984)). It is important to notice that the cell mesh was extended far away from the cavity in order to avoid the undesirable reflection effects produced by the  artificial outer boundary (as mentioned before, this is the main limitation of the D-BEM formulation; for this reason, in this example the outer boundary is a circle whose radius is equal to 9.4 times R). In order to simulate the plane stress analysis carried out by Chow and Koenig (1966), assumed here as the analytical one (m = 0.3 and E = 100.0), the following parameters have been adopted in both BEM plane strain analyses: m = 0.2308 and E = 94.6769.
Results corresponding to the circumferential and radial stress components at internal point A(2.02R, 0) are shown, respectively, in Figs. 17 and 18 for the Houbolt-a scheme with a = )0.10, and in Figs. 19 and 20 for the Newmark-a scheme with a = 0.20. In both analyses b Dt ¼ 3=4 was adopted. Note again the benefic influence of the a parameter in the analysis carried out with the Newmark scheme, producing the reliable results in Figs. 19 and 20. On the other hand, the results furnished by the Houbolt-a

Conclusions
This work is concerned with presenting some improvements, related to the time-marching process, for the D-BEM formulation. Essentially, the proposed procedures consist in the association between Newmark and Houbolt schemes with the Hilbert-a scheme. Better control of the numerical damping (obtained by taking a < 0 in the Houbolt scheme) and of the spurious oscillations (obtained by taking a > 0 in the Newmark scheme) is the    produced, validating the proposed formulation and stimulating the continuity of this research work, which means its application to the DR-BEM and to 3-D elasto-dynamics, as well as to non-linear problems. Another topic that deserves attention is concerned with the use of time variable values of the a parameter. Finally, it is important to mention that, according to the authors' point of view, some theoretical studies are still required when performing dynamic analysis with the BEM: these studies are related with the determination of the time-step length and with the determination of the a parameter and become necessary to confirm the empirical values adopted so far.