Efficiency improvement of the frequency-domain BEM for rapid transient elastodynamic analysis

Efficiency improvement of the frequency-domain BEM for rapid transient elastodynamic analysis

Jinyou Xiao J. Xiao College of Astronautics, Northwestern Polytechnical University, Xi’an 710072, China
22email: xiaojy@nwpu.edu.cnW. Ye Department of Mechanical Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong
44email: mewye@ust.hk
   Wenjing Ye J. Xiao College of Astronautics, Northwestern Polytechnical University, Xi’an 710072, China
22email: xiaojy@nwpu.edu.cnW. Ye Department of Mechanical Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong
44email: mewye@ust.hk
   Lihua Wen J. Xiao College of Astronautics, Northwestern Polytechnical University, Xi’an 710072, China
22email: xiaojy@nwpu.edu.cnW. Ye Department of Mechanical Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong
44email: mewye@ust.hk
Received: date / Accepted: date

The frequency-domain fast boundary element method (BEM) combined with the exponential window technique leads to an efficient yet simple method for elastodynamic analysis. In this paper, the efficiency of this method is further enhanced by three strategies. Firstly, we propose to use exponential window with large damping parameter to improve the conditioning of the BEM matrices. Secondly, the frequency domain windowing technique is introduced to alleviate the severe Gibbs oscillations in time-domain responses caused by large damping parameters. Thirdly, a solution extrapolation scheme is applied to obtain better initial guesses for solving the sequential linear systems in the frequency domain. Numerical results of three typical examples with the problem size up to 0.7 million unknowns clearly show that the first and third strategies can significantly reduce the computational time. The second strategy can effectively eliminate the Gibbs oscillations and result in accurate time-domain responses.

Elastodynamics Fast boundary element method Wave propagation Exponential window method Frequency domain windowing

1 Introduction

Transient analysis is encountered in almost every field in science and engineering. For most real-world problems, numerical simulation is the only viable approach for obtaining accurate solutions. The boundary element method (BEM) has been well-recognized as a powerful numerical method to treat wave propagation problems. Its reduced dimensionality, high accuracy and its ability in capturing rapid transitions and steep gradients of the fields have made this method particularly suitable for problems with complex geometry, for example, dynamic analysis of porous solids and fracture analysis. In view of the increasingly complex systems that we have encountered nowadays, the development of efficient BEM for large-scale wave propagation simulations has become indispensable.

According to the different approximation solution strategies in time space, the BEM for treating elastodynamic wave propagation problems generally follows two approaches, namely, time-domain approaches and frequency-domain approaches (see e.g. the reviews by Beskos Bde97 () and Costabel CM04 ()). Time-domain approaches can be further classified into time-stepping methods and the space-time integral equation method. In these methods the physical problems are directly solved in the real time domain, thus one can observe the phenomenon as it evolves. However, such methods require an adequate choice of the time step size. An improper time step could lead to instability or numerical damping. For recent development of time-domain methods, see e.g. SA97 (); BS12 (); ADF12 ().

In the frequency-domain approach, the frequency-domain boundary integral equation is established and solved at a series of sampling frequencies. The results are then transformed back into time domain by inverse Laplace transform which can be numerically evaluated by the discretized Fourier transform (DFT) AM87 (); PTB05 (); PGS10 (). This approach is attractive for transient analysis because it is: (1) stable in obtaining time-domain solutions, (2) simple in implementation and (3) suitable for parallel solution. In addition, one notes that the frequency-domain problems by itself already has a broad range of applications such as harmonic wave propagation in solids Bde97 (), air damping on micro resonators DY04 (), heat transfer SST12 (), to name a few.

There are two issues need to be dealt with when computing the transient response using the inverse Fourier transform. The first issue is regarding to the discretization and truncation errors introduced by the DFT MR08 (). The discretization in DFT causes wrap-aroud (or aliasing) error. For damped systems this error can be removed by adding trailing zeroes to damp out the free vibration before the end of the period. However, this scheme is inefficient for low damped systems, and more importantly, it does not work for non-dissipative systems such as the elastic solids considered in this work. To circumvent this problem, the exponential window method (EWM) KR92 (); Humar02 () was employed in the BEM transient elastodynamic analysis PGSG11 (); XYCZ12 (). This method adds an artificial damping to the system to effectively damp out the free vibrations. By properly choosing the damping parameter, denoted by , in the EWM, the response period can be arbitrarily set independent of the actual time-domain response. In addition, it has been found that a large can improve the conditioning of the BEM coefficient matrix. However one cannot arbitrarily increase because it also amplifies the truncation error of DFT and therefore causes unacceptable Gibbs oscillations in the later time response. In this paper, a frequency-domain windowing technique is introduced to eliminate the Gibbs oscillations, and consequently accurate time response can be obtained in all simulation periods. The frequency-domain windowing technique in combine with the EWM leads to a simple but efficient modified Fourier transform approach for BEM elastodynamic transient analysis.

The second issue concerns the computational efficiency. In the frequency-domain approach the transient problem is transformed into independent time-harmonic problems. The computational cost grows almost linearly with . For many practical problems sampling frequencies are often required to obtain accurate time-domain solutions. This implies a huge computational cost because solving one frequency-domain eqation in a three-dimension domain using classic BEM is already computationally very costly. To resolve this issue, a series of work has been recently published on the development of efficient integral equation approaches for frequency-domain elastodynamics. Examples include, but are not limited to, fast multipole SBD08 (); TC09 (); CBS08 (), -matrix BA10 () and precorrected-FFT (pFFT) accelerated BEM YZY10 (), as well as efficient BEM with adaptive cross approximation MS08 (). In this paper, the pFFT technique is employed to accelerate the frequency-domain BEM analysis XYCZ12 () due to its ease in implementation and high computational efficiency.

A key factor that affects the computational cost of the frequency domain approach is the number of iterations required to solve each linear system. In order to achieve further speedup, a least square extrapolation method is proposed to obtain better initial guesses for the iterative solution procedure.

The basic setting of the frequency-domain approach for BEM transient analysis is briefly recalled in Section 2. The modified Fourier transform method for transient elastodynamic anaylsis is proposed in Section 3. In Section 4 the precorrected-FFT and the extrapolation method for initial guess are presented. The performance of the methods put forward in this paper is tested by three examples in Section 5, and a summary is drawn in Section 6.

2 Frequency domain BEM for transient analysis

2.1 Frequency domain boundary element method

The first step in the frequency-domain approach for transient elastodynamic analysis is to transform the time domain into the frequency domain via Fourier transform. Let denote the region of space occupied by a three-dimensional elastic solid with isotropic constitutive properties defined by Lamé constants and , Poisson’s ratio and mass density . In a mixed boundary value problem, the boundary consists of a displacement boundary and a traction boundary such that . The displacement governing equation for linear elasticity in the frequency domain is formulated as


where, is the circular frequency; is the displacement in the frequency domain; and denote the Nabla and the Laplace operators, respectively; is the traction operator defining the stress-strain relation based on Hooke’s law; and are prescribed boundary displacement and traction on and , respectively.

The boundary integral equation corresponding to equation (1) is formulated as


where, (P.V.) indicates a Cauchy principal value (CPV) of the singular integral; the free-term is equal to for a smooth boundary at ; and denote the th components of the displacement and traction generated at by a unit point force applied at along the direction , given by




is the unit normal to at directed outwards of ;

and are the respective wavenumbers of S and P elastic waves, so that

The boundary integral equation (2) is numerically solved using a surface discretization with triangular boundary elements. Piecewise-constant approximation is employed in this work. A standard collocation approach yields a discretized system of the BIE (2) of the form


where matrices and with dimensions of by , , consist of boundary integrals of kernel functions and on each element respectively; the -vectors and contain the nodal displacement and traction at each element. Note that the contribution of the free-term has been absorbed into the matrix . By enforcing the boundary conditions, the linear system (6) becomes


where, the -vector collects the unknown nodal displacement and traction components and can be obtained by solving the above linear system. All the matrix and vectors in (7) are functions of frequency . The computational cost for an iterative solver is . This can be greatly reduced by using the pFFT method in section 4.1.

2.2 Conventional frequency-domain approach

Equation (1) or (2) is solved at a series of discrete frequencies. The resulting displacements and tractions are then transformed back to the time domain by using the inverse DFT to obtain the time-domain response.

Let and be the circular frequency and time resolutions, respectively; and be the time period of the transient response. Given the number of sampling points , one has the basic relations and . Note that the Nyquist frequency should be chosen such that responses corresponding to frequencies higher than are insignificant and can thus be discarded.

The procedure for obtaining the transient responses using the frequency-domain approach can be summarized as:

  1. Determine a frequency resolution and the number of sampling points . needs to be small enough to minimize the loss of useful frequency information.

  2. Sample boundary condition , which could be the given boundary displacement or the traction, and perform DFT to obtain the frequency-domain boundary conditions at frequencies ,

  3. Conduct analysis for the first frequencies, that is, , to obtain the frequency-domain response for the first frequencies. The responses corresponding to the last can be determined utilizing the conjugate symmetric property of the DFT as,

    where, and denotes the complex conjugate of .

  4. Perform IDFT for to obtain time-domain response , i.e., the desired boundary displacement and/or traction,


Discretization and truncation errors arise when the forward and inverse Fourier transforms are numerically evaluated using DFT and IDFT in (8) and (9). The discretization error is also known as aliasing error. To eliminate it, the time-domain signal is required to be zero for . The truncation error is due to the truncation of the frequency spectrum. It causes high frequency oscillations near the discontinuities of the signals, known as Gibbs oscillations.

In next sections, techniques for reducing the two errors are introduced. Consequently, a general and efficient modified Fourier transform (MFT) method for elastodynamic transient analysis is obtained. The MFT method itself is not new. It was first developed for analyzing power system transients by Mullineux’s group about four decades ago; see MR08 () and the references therein. In this paper, it is applied in elastodynamic transient analysis.

3 Modified Fourier transform method

In the MFT method, two data processing techniques, one with respect to the time-domain data and another to the frequency-domain data, are introduced to reduce the discretization and truncation errors associated with the DFT and IDFT.

3.1 Discretization error & time-domain damping

As aforementioned, the periodic nature of the DFT requires that the time-domain response tends to zero at the end of a period. This however often cannot be satisfied in the dynamic analysis using the frequency domain approach. A common way to circumvent this problem for a dissipated system is to add trailing zeroes to prolong the duration so that the free vibration can be damped out before the end of the period. Such an approach is clearly not applicable for a non-dissipated system. Another more general method suitable for all systems is the exponential window method (EWM) KR92 (); Humar02 (). In this method, a complex frequency shift, amounting to an artificial damping, is introduced to efficiently damp out the free vibrations within the selected period. The actual response is obtained after a proper scaling of the damped response.

Recently, the EWM was used in 2D and 3D BEM elastodynamic transient analysis PGSG11 (); XYCZ12 (). Accurate time domain response can be obtained by a proper choice of the damping parameter . To be specific, let be the displacement in the time domain and its Fourier transform is in (1). In applying the EWM, the displacement is scaled by an exponential window function to obtain a damped displacement , i.e.,


Since is a bounded function, the damped response tends to be zero for large . The Fourier transform of is , meaning that can be obtained by simply replacing the frequency in by , where . To obtain the transient response, one needs to solve the problem at a series of complex frequencies . Once the damped response is obtained, the actual response can be obtained as


The choice of the damping parameter determines the performance of the EWM. On one hand, it cannot be too large since the factor in (11) acts as an amplifier which magnifies the cumulative errors due to the truncation of the frequency space which will be discussed in section 3.2 and BEM error. On the other hand, given the factor is used to damp the time-domain response, it could be supposed that a high value for is required to sufficiently eliminate the aliasing error. A compromise between the two ends leading to the following rule-of-thumb for the selection of KR92 ()


This rule was used in the work described in PGSG11 (); XYCZ12 (). It was found in XYCZ12 () that in BEM analysis a choice of the coefficient between could often lead to good time-domain results.

In addition to reducing the aliasing error, another important effect of lies in that a positive can improve the conditioning of the system matrices in (6) and (7). This is a nature result of the Gershgorin circle theorem, because the exponential function reduces the absolute values of the non-diagonal entries of the matrices, thus making the eigenvalues more clustered. Hence, the damping factor actually acts as a perconditioner (called artifical damping perconditioner, ADP, in this paper) whose performance can be enhanced by increasing . As such, a larger is always desirable.

Matrix preconditioning is a crucial problem in BEM analysis. This is particularly true for high frequency and multi-domain problems, in which the coefficient matrices tend to be much more ill-conditioned. Various preconditioners exist (e.g., INYN12 (); CSB12 ()) and can be used. However, it is well known that the performance of many existing preconditioners is generally problem-dependent and there always is a trade-off between their performance and computational cost. On the contrary, the ADP is simple yet efficient, because its implementation requires no extra computation and memory. Whereas, it can largely increase the convergence speed of the iterative solvers; see Section 5.

Meanwhile, the ADP can be used in combination with other preconditioners. For instance, in this work a block diagonal preconditioner is used together with the time-domain damping technique described in this section to solve (7).

3.2 Truncation error & frequency-domain windowing

It is observed that choosing in BEM analysis can cause unacceptable Gibbs oscillations in time-domain responses at later times XYCZ12 (). Here, a frequency-domain data windowing technique is introduced to alleviate the the Gibbs oscillations. Consequently, a larger value of can be set. This implies a considerable improvement in the conditioning of the BEM matrices when compared with ; see evidences in Section 5.

The implementation of this technique is rather simple MR08 (). Before the damped frequency-domain response is transformed back into the time domain using IDFT, it is first multiplied by a suitable data window ; that is, the function is transformed in lieu of itself. The function of the frequency-domain windowing is to reduce the Gibbs oscillations by averaging out the oscillations within one period.

There exist a number of window functions, see e.g., MR08 (); CC95 (); CC98 (). As far as BEM elastodynamic analysis is concerned, we find that the Hanning and Blackman windows often yield satisfactory results. The discretized window functions are given by


3.3 A summary

To summarize, the time-domain damping in conjunction with the frequency-domain windowing techniques leads to a MFT method for elastodynamic transient analysis. The basic implementation procedure of the MFT method is similar to that described in section 2.2. The distinction lies in that in the MFT method the sampling frequencies are instead of , and the DFT formula (8) and IDFT formula (9) are replaced by





4 Rapid solution of linear systems

The main computational cost for the frequency-domain approach is the solution of linear systems (7) for a sequence of sampling frequencies , i.e., the parameterized linear systems


where, , and are defined analogously. In this section, two methods are described to considerably reduce the computational cost, namely, a pFFT method for accelerating the BEM matrix-vector multiplication, and a frequency extrapolation method for obtaining good initial guesses and thus reducing the total number of iterations in solving the sequential systems.

4.1 PFFT accelerated elastodynamic BEM

The computational cost of a matrix-vector multiplication in traditional BEM is , which could be prohibitively expensive for large-scale elastodynamic simulations. Within the last two decades, significant progress has been made in the development of accelerated BEM techniques. With the accelerated techniques, the computational complexity, including both the CPU time and memory usage, can be considerably reduced, therefore enabling BEM-based large-scale computation.

The first accelerated technique used in elastodynamic BEM is the FMM F98 (), which is further developed by several groups; see, e.g., TC09 (); CBS08 (); TNK03 (). In BA10 () the -matrix method is applied in frequency-domain elastodynamic BEM. Recently the pFFT method was used to accelerate elastodynamic BEM analysis by the present authors XYCZ12 (). The pFFT method is very easy to implement and to be parallelized, so a high computational efficiency can be achieved with ease. For more information about the pFFT technique and the details of its implementation, the readers are referred to XYCZ12 ().

Theoretically, for problem domains with small surface-to-volume ratios such as porous solids, optimal order for CPU time and for memory can be achieved by pFFT, which is comparable to the FMM and the -matrix method. For other problems, the CPU time of pFFT is between and . In these cases, the FMM and the -matrix method would be advantageous. However numerically, due to the small overhead, the actually computational cost of pFFT could be less than that of FMM and the -matrix method for problems with up to order ; see XYCZ12 () for comparisons.

4.2 Solution extrapolation method

Iteratively solving each linear system in the sequence (15) needs an initial guess, denoted by . Typically, is independently chosen as zero or one vectors. Here a least square extrapolation method is proposed to determine a better initial guess for the solution of based on the calculated solutions at previous frequencies in the sequence; this approach is termed as the solution extrapolation method (SEM).

The motivation is inspired by the fact that solution depend on the frequency . If is a smooth function of and that the frequency step is sufficiently small, it should be possible to obtain by extrapolation from the previous solutions; that is,


where, , . Let . Then the coefficient vector can be obtained by minimizing the square error


By computing the QR-decomposition , one gets the least square solution of (17)

Thus, the initial guess is can be computed as


The main computational cost for evaluating is the computation of matrix-vector products .

Optimal value for the number of preceding solutions in (16) clearly depends on the smoothness of the solutions and thus is problem dependent. For the cases studied in this work, it is found that is optimal. Generally, a large does not improve the accuracy of the initial guess, because the solution is often not smooth over a large frquency range and the frequency step may not be so small as required.

5 Numerical examples

In this section, three examples are selected to validate and demonstrate the performance of the MFT method in Section 3 and the SEM for initial guess in Section 4.2. The first example is a classical benchmark for transient elastodynamic analysis: a prismatic rod subject to a sudden applied uniform pressure at one end. The analytic solution is available and is employed to evaluate the accuracy and the efficiency of the method. The second case is an aluminum plate subject to an impact loading. This example is chosen to test the method for problems with a time-dependent arbitrary loading. Finally to test the ability of the method for solving large-scale problems, a unit cube containing 64 spherical cavities subject to a sudden applied uniform pressure is simulated.

In the following simulations, constant elements are used. The pFFT method in XYCZ12 () is used to accelerate the BEM analysis. All the linear systems are preconditioned using the block diagonal preconditioner and solved by the generalized minimal residual method (GMRES) SS86 () with convergence tolerance . Calculations are performed on a workstation with a Xeon 3.0GHz CPU and 30GB RAM.

5.1 Example 1: prismatic rod

As the first test case, a prismatic rod with dimensions of as sketched in Figure 1 is simulated. The rod is fixed at its left end and subject to a Heaviside traction = at its right end, where N/m and is the Heaviside function. By setting the Poisson’s ratio of the material to zero, the 3-D problem is reduced to a 1-D problem, of which the analytical solution is available, see XYCZ12 (). The Young’s modulus and the density of the material are set to be Pa and kg/m, respectively.

Figure 1: Prismatic rod subject to a step traction pulse

5.1.1 Effects of the ADP and the SEM

The boundary of the prismatic rod is partitioned into 12800 triangular elements. The response period and the number of samplings are set to be s and , respectively. To test the effects of the damping coefficient and the SEM on accelerating the convergence of the linear system solver, simulations are carried out under the three cases specified by the first three columns of Table 1. For cases without using the SEM, the initial guesses are set to be zero vectors.

Cases SEM Total its Total time (h)
Case 1 2 29237 16.0
Case 2 4 21758 14.4
Case 3 4 ✓, 17562 12.1
Table 1: Three test cases and their CPU times
Figure 2: Total number of iterations in solving the linear systems to the sampling frequencies

The number of iterations for solving the sequential linear systems (15) corresponding to the first 129 sampling frequencies are plotted in Figure 2. The total number of iterations and the total CPU time for solving all the 129 systems in each of the three cases are listed in the last two columns of Table 1. One can see that by using a larger damping with , a 25% reduction in the total number of iterations can be achieved when compared with the case of . By further using the SEM, the reduction is up to about 40%. However, the reduction of the total CPU time is not as considerable as the number of iterations. This is because the matrix-vector multiplication in this problem is relatively fast and the preprocessing time (i.e., time for computing the near-field interations) occupies a large portion.

5.1.2 Effects of the frequency windowing technique

As described in Section 3.2, the frequency windowing technique can be employed to eliminate the Gibbs oscillations in the time-domain response. Here this is verified using the results of case 3, i.e. , in Table 1. Our numerical experience shows that both the Blackman and Hanning windows obtains almost the same results. Thus only the results of Hanning window are presented here.

Figure 3: Displacement response at the free end of the rod

(a) Total time response

(b) Local magnification

Figure 4: Traction response at the fixed end of the rod

Figure 3 and 4 demonstrate both the analytical and simulated displacements at the free end and the tractions at the fixed end as functions of time; the non-dimensional time (, with and being the wave speed and length of the rod, respectively) is employed in the plots. It is observed that, before using the frequency windowing technique, both the displacement and traction responses show drastic oscillations in later times. The magnitude of the traction oscillations is about 4 times larger than the largest traction value; see Figure 4 (a). But after the frequency windowing technique is used, the red curves obviously show that the oscillations are effectively removed and the resulting responses agree very well with the analytical solutions in all the simulation period; see Figure 4 (b).

5.1.3 Long time behavior of the method

It is well known that time-domain BEM analysis often suffers from either strong numerical damping or instabilities when the calculation time period is long. Here, the long-time behavior of the present frequency domain method is studied. The simulation conditions are the same as the case 3 in Section 5.1.1; but the response period and the number of samplings are set to be s and , respectively. Thus, the corresponding time step is s and 257 frequency domain BEM analyses are needed.

The computed traction history is shown in Figure 5. Clearly, the proposed method shows a very nice behavior in the entire simulation period. Computed traction responses of similar problems can be found in Figure 12 of BMS12 () and Figure 9 of ADF12 (). It seems that the total numbers of time steps required in their work are much larger than that in the present simulation. For example, the number of time steps for obtaining Figure 9 of ADF12 () is about , whereas in the present method only 257 frequency steps corresponding to 257 time steps are needed.

Figure 5: A long time traction response at the fixed end of the rod

5.2 Example 2: a plate with a hole subject to an impact loading

The performance of the present method is further demonstrated by a more realistic problem; i.e., an aluminum plate with a hole subject to a vertical impact force as illustrated in Figure 6. The material properties of the aluminum plate are: GPa, kg/m and . For more details of the simulation conditions and the load history, the readers are referred to Section 4.2 of XYCZ12 (). In the present computation, a mesh with 14072 triangular elements is used, s, . The damping coefficient is set to be 4, the SEM with is used. The Hanning window is employed for frequency windowing.

Figure 6: A plate with a hole (unit: mm)

The responses of displacement and strain in the load-direction at points S (see Figure 6) are computed using the present method and compared with the FEM and/or experiment results in Section 4.2 of XYCZ12 (); see Figure 7. By using the frequency windowing technique the Gibbs oscillations in both displacement and strain responses in the later times are effectively removed. The windowed results coincide very well with the FEM results. The agreement of the computed strain with the experimental measurement is also quite good overall.

The total CPU time for this simulation is about 5.4 hours. The memory usage is 840 MB, slightly more than the 723 MB in XYCZ12 (), because here we use GMRES solver which is more memory consuming than the IDR(s) solver used in XYCZ12 ().

(a) Displacement

(b) Strain

Figure 7: The displacement and strain responses at point S

5.3 Example 3: a cube with 64 spherical cavities

To demonstrate the performance of the presented method for solving large-scale elastodynamic wave propagation problems, a model consisting of a solid cube with 64 identical spherical cavities was built and simulated; see Figure 8. The non-overlapped spherical cavities are randomly distributed inside the cube and their radius is 0.072 m. The cube is fixed at the bottom surface () and is subject to a uniform step traction at the upper surface (). The load is the same as that in Section 5.1. All other surfaces are traction free. The material is aluminum and the properties are the same as those defined in Section 5.2. This example was also used in XYCZ12 (), but here a finer mesh with elements is used; hence, the total degrees of freedoms is nearly 0.7 million. The response period and the sampling size are set to be s and , respectively. Damping coefficient and the SEM with are used. The Hanning window is employed in frequency windowing.

Figure 8: A unit cube containing 64 spherical cavities subjecting to a step load

The simulation takes about 101 hours, and consumes 21 GB memory. The CPU time still seems to be quite large. Further reduction of computational expense can be achieved by employing a more effective pre-conditioning method and by parallel computing. The time responses of the displacement and traction in -direction at four boundary points are plotted in Figure 9. Compared with the 1D solutions of the solid rod presented in Section 5.1, the oscillation patterns of the displacement and traction seem reasonable.

(a) Displacement

(b) Traction

Figure 9: The displacement and traction responses for a few selected points

6 Summary

Fast solvers are indispensable for large-scale wave propagation analysis. Recently, a fast frequency domain pFFT BEM incorporated with the exponential window technique was developed for solving transient elastodynamic problems XYCZ12 (). It has been shown that this method is stable and general in the sense that: (1) the method can be used to obtain accurate time responses for an arbitrarily chosen time period; (2) the method is applicable to solving problems with damping, either small or large, as long as the damping term is linear and the resulting frequency-domain equation has an equivalent integral formulation

In this work, the computational efficiency of the method in XYCZ12 () is further improved via three approaches:

  1. Introducing large damping via the exponent window function. It can be proved that the damping actually plays the role of preconditioner, denoted as the artificial damping preconditioner (ADP) in this paper. Applying a large damping can in general improve the conditioning of the BEM matrices, thus reducing the total number of iterations required to achieve convergence. Numerical results show that the reduction is around 25% for compared with .

  2. Using frequency windowing technique in Fourier inversion. Exponential function with a large damping coefficient would cause unacceptable oscillations in the later time responses. The oscillations stem from the Gibbs phenomenon which is inherent in expanding non-periodic functions using Fourier series. The frequency domain windowing technique is an effective way to abate the Gibbs phenomenon and is employed in the Fourier inversion of the frequency domain BEM results. It is found that the well-known Blackman and Hanning windows perform equally well in removing the Gibbs oscillations in the BEM solutions. The performance of the Hanning window is demonstrated by numerical results in Section 5. The time-domain damping in conjunction with the frequency-domain windowing techniques results in a modified Fourier transform method for elastodynamic wave propagation simulation.

  3. Using solution extrapolation method (SEM) for obtaining good initial guess. In this method, the initial guess for the current linear system is obtained from the solutions of the linear systems corresponding to several previous sampling frequencies via extrapolation. The aim is to reduce the number of iterations for iterative solvers. Numerical results indicate that the reduction can be up to 15% and is generally problem dependent.

The effectiveness of the above strategies is demonstrated by three examples. The largest model has nearly 0.7 million DOFs, and is simulated on a workstation with a Xeon 3.0GHz CPU and 30GB RAM.

JX and LW were supported by the NSFC under Grant 11102154 and 11074201, and the New Teacher Fund for Doctor Station from the Chinese Ministry of Education under Grant 20106102120009. WY was supported by Hong Kong Research Grants Council under Competitive Earmarked Research Grant 621411.


  • (1) Beskos DE. Boundary element methods in dynamic analysis: Part II (1986-1996). Applied Mechanics Reviews 1997; 50(3):149–197.
  • (2) Costabel M. Time-dependent problems with the boundary integral method, in: Ency- clopedia of Computational Mechanics, Stein E, de Borst R., Hughes TJR (eds), Wiley, 2004.
  • (3) Schanz M, Antes H. Application of ‘operational quadrature methods’ in time domain boundary element methods. Meccanica 1997; 32(3):179–186.
  • (4) Banjai L, Schanz M. Wave Propagation Problems Treated with Convolution Quadrature and BEM. Lecture Notes in Applied and Computational Mechanics 2012; 63:145–184.
  • (5) Aimi A, Diligenti M, Frangi A, Guardasoni C. A stable 3D energetic Galerkin BEM approach for wave propagation interior problems. Engineering Analysis with Boundary Elements 2012; 36(3):1756–1765.
  • (6) Ahmad S, Manolis GD. Dynamic analysis of 3-D structures by a transformed boundary element method. Comput Mech 1987; 2:185–196.
  • (7) Polyzos D, Tsepoura KG, Beskos DE. Transient dynamic analysis of 3-D gradient elastic solids by BEM. Computers and Structures 2005; 83:782–792.
  • (8) Phan AV, Gray LJ, Salvadori A. Transient analysis of the dynamic stress intensity factors using SGBEM for frequency-domain elastodynamics. Computer Methods in Applied Mechanics and Engineering 2010; 199(45-48):3039–3050.
  • (9) Ding J, Ye W. A fast integral approach for drag force calculation due to oscillatory slip stokes flows. International Journal for Numerical Methods in Engineering 2004; 60(9):1535–1567.
  • (10) Simões I, Simões N, Tadeu A, Reis M, Vasconcellos CAB, Mansur WJ. Experimental validation of a frequency domain BEM model to study 2D and 3D heat transfer by conduction. Engineering Analysis with Boundary Elements 2012; 36(11):1686–1698.
  • (11) Moreno P, Ramirez A. Implementation of the Numerical Laplace Transform: A Review. IEEE Transactions on Power Delivery 2008; 23(4):2599–2609.
  • (12) Kausel E, Roësset JM. Frequency domain analysis of undamped systems. Journal of Engineering Mechanics 1992; 118(4):721–734.
  • (13) Humar JL. Dynamics of structures. 2nd edition. Balkema, 2002.
  • (14) Phan AV, Guduru V, Salvadori A, Gray LJ. Frequency domain analysis by the exponential window method and SGBEM for elastodynamics. Computational Mechanics 2011; 48(5):615–630.
  • (15) Xiao J, Ye W, Cai Y, Zhang J. Precorrected FFT accelerated BEM for large-scale transient elastodynamic analysis using frequency-domain approach. International Journal for Numerical Methods in Engineering 2012; 90(1):116–134.
  • (16) Sanz JA, Bonnet M, Domínguez J. Fast multipole method applied to 3-D frequency domain elastodynamics. Engineering Analysis with Boundary Elements 2008; 32:787–795.
  • (17) Tong MS, Chew WC. Multilevel fast multipole algorithm for elastic wave scattering by large three-dimensional objects. Journal of Computational Physics 2009; 228(1):921–932.
  • (18) Chaillat S, Bonnet M, Semblat JF. A multi-level fast multipole BEM for 3-D elastodynamics in the frequency domain. Comput Methods Appl Mech Eng 2008; 197:4233–4249.
  • (19) Benedetti I, Aliabadi MH. A fast hierarchical dual boundary element method for three-dimensional elastodynamic crack problems. International Journal for Numerical Methods in Engineering 2010; 84(9):1038–1067.
  • (20) Yan ZY, Zhang J, Ye W. Rapid solution of 3-D oscillatory elastodynamics using the pFFT accelerated BEM. Engineering Analysis with Boundary Elements 2010; 34(11):956–962.
  • (21) Messner M, Schanz M. Adaptive cross approximation in an elastodynamic boundary element formulation. PAMM 2008; 8:10309–10310.
  • (22) Isakari H, Niino K, Yoshikawa H, Nishimura N. Calderon’s preconditioning for periodic fast multipole method for elastodynamics in 3D. International Journal for Numerical Methods in Engineering 2012; 90(4):484–505.
  • (23) Chaillat S, Semblat J, Bonnet M. A preconditioned 3-D multi-region fast multipole solver for seismic wave propagation in complex geometries. Commun Comput Phys 2012; 11(2):594–609.
  • (24) Chen LY, Chen JT, Hong HK, Chen CH. Application of Cesaro mean and the L-curve for the deconvolution problem. Soil Dynamics and Earthquake Engineering 1995; 14(5):361–373.
  • (25) Chen JT, Chen CH. Analytical study and numerical experiments for Laplace equation with over specified boundary conditions. Appl. Math. Modelling 1998; 22:703–725.
  • (26) Fujiwara H. The fast multipole method for the integral equations of seismic scattering problems. Geophysical Journal International 1998; 133(3):773–782.
  • (27) Takahashi T, Nishimura N, Kobayashi S. A fast BIEM for three-dimensional elastodynamics in time domain. Engineering Analysis with Boundary Elements 2003; 27:491–506.
  • (28) Saad Y, Schultz MH. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Sci and Stat Comput 1986; 7(3):856–869.
  • (29) Banjai L, Messner M, Schanz M. Runge-Kutta convolution quadrature for the Boundary Element Method. Comput. Methods Appl. Mech. Engrg. 2012; 245-246(3):90–101.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description