In part I of the present article the formulation for a dynamic stationary semianalytical solution for a spatially constant load applied over a rectangular surface within a viscoelastic isotropic full-space has been presented. The solution is obtained within the frame of a double Fourier integral transform. These inverse integral transforms must be evaluated numerically. In the present paper, the technique to evaluate numerically the inverse double Fourier integrals is described. The procedure is validated, and a number of original displacement results for the stationary loading case are reported.

1. Introduction

In an accompanying paper [1], closed-form solutions, in the wave number domain, for displacements within a viscoelastic homogeneous isotropic full-space excited by rectangular loadings are furnished. Expressions for the solutions in the original physical domain were obtained with the use of the inverse double Fourier integral transform. The resulting expressions have to be integrated numerically. The present second part of the paper describes the strategy used to perform numerically the inverses of the double Fourier integral transforms. The strategy is validated by comparisons with known special cases. A number of original numerical results for displacements within the three-dimensional full-space under spatially rectangular loads of constant amplitude, harmonically varying in time are furnished.

2. Numerical Integration

2.1. Viscoelasticity and Singular Integrands

In the first part of this paper [1], a set of general expressions for the displacement response of a 3D isotropic viscoelastic full-space were obtained. Equation (2.1) shows a typical expression that needs to be evaluated. It is a typical expression resulting from the inversion of the double Fourier integral transform. In this case it describes the vertical displacement of the full-space due to a vertical load applied over a rectangular area with dimensions (see Appendix for notation), :

Equation (2.1) presents two improper integrations, to be performed over the wave number variables and . An outer improper integration is performed over the variable . For every value of , an improper inner integration over the variable must be evaluated. Now the structure of (2.1) is analyzed.

The integrand of (2.1) represents, without loss of generality, the behavior of the integrands of all the remaining displacement terms furnished in part I [1]. These integrands are composed by a kernel multiplied by oscillating functions and . For the elastic case, these kernels, here represented by in (2.2), present two singularities, which correspond to the dilatational and shear waves which can propagate in the isotropic full space.

Figures 1(a) and 1(b) show the behavior of the integrand kernel as a function of the integration variable for the following set of parameters: Lamé constant and Poisson’s ratio , loading half-width , nondimensional frequency , and damping coefficient ranging from (elastic case) to . It can be shown that, for values of the integration variable larger than 1, the absolute value of the kernel function decays monotonically.

The singularities corresponding to the shear and dilatational waves can be observed at and , in Figures 1(a) and 1(b), respectively. It is also observed that these singularities are smoothed when internal damping is introduced. In the present work viscoelastic behavior is assumed and the elastic constants are made complex, according to the elastic-viscoelastic correspondence principle [1, 2]. As a consequence, the singularities in the kernels are smoothed and no integration over a singularity must be evaluated [3].

Now consider the improper inner integration over , described in (2.1). The difference between the behavior of the full integrand and that of its kernel is the oscillating nature of the former. Figure 2 shows the behavior of the complete integrand for the inner evaluation given in (2.1). As can be seen, the integrand presents an oscillatory and decaying nature. The singularity remains at the points and .

Figures 1 and 2 characterize the line of inner improper integration over for the value of the outer improper integration variable. The integrands present the same behavior with respect to the variable . In Figure 3, the behavior of the kernel is shown for the plane of double-integration variables and . Notice that, along the planes and , the integration kernels reproduce the behavior shown in Figure 1.

2.2. Numerical Inversion of the Fourier Integrals

Based on the behavior of the integrands shown in the previous section, the following methodology of integration is devised. The plane of integration is divided into 3 regions as shown in Figure 4.

Region I is the region in which the integrands may present a singularity. Mathematically, this region is determined by the ranges and . Actually it can be shown that the singularities in the plane are to be found at the locations determined by the equations and (see Figure 3). In the present implementation the values have been chosen. It can also be shown that the oscillatory character of the integrand increases for larger values of the distance , and , as well as with the increase of the nondimensional frequency parameter .

In this first and limited integration region, the integrand presents a quasi-singular character, superposed by an oscillating behavior. An adaptive 2D Gaussian quadrature scheme is applied to evaluate the integrands. The number of integration points is systematically increased until a degree of convergence is obtained.

In Region II the integration is proper in one direction and improper in the other direction . There is no singular or quasi-singular behavior of the integrands in this region. There is an oscillatory behavior in the range and a monotonic decreasing oscillatory behavior in the direction of the improper integration . In the direction, standard 1D Gaussian quadrature is employed. For the direction, in which the integrand is oscillating and monotonically decaying, the scheme proposed by Longman [3] for improper integrals is applied. This is a very efficient integration strategy. Recently the authors of this paper developed a scheme to implement Longman’s strategy on general-purpose graphics processing units (GPGPUs) [4]. The implementation of a GPGPU was able to increase the integration performance by almost 900 times, making it feasible to evaluate this inverse double integral transforms on commodity computers.

In Region III , the integrands oscillate and decay monotonically in both directions, and Longman’s integration strategy is used in these two directions.

3. Validation

The methodology described in the previous section was used to determine the displacements within the full space due to a dynamic load uniformly distributed over an area of sides . For the synthesized solutions, the origin of the coordinate system is at the center of the loading area. The solutions depend on many parameters such as the coordinate of the point at which the solutions is calculated, , and , on the dimensionless frequency , on the constitutive parameters and , on the dimensions of the loading area and , and also on the internal damping coefficient .

The present solution can be partially validated using the limiting case of a static solution due to a concentrated point load . Kane [5] presents a classical closed-form Green’s function for concentrated static loads in the interior of the 3D isotropic full space, known as Kelvin’s solution. To simulate the case of a concentrated static load, the nondimensional frequency was set to and the half-widths of the loading were set to . The other parameters are , , and . Figures 5(a) to 5(i) show the results obtained by the Kelvin solution and by the present approach (2-Fourier). The 9 components of the displacement tensor are determined. The solutions were obtained along the line of the full space given by and .

The corresponding Kelvin solution for the dynamic load case is presented by Kitahara [6]. This dynamic Kelvin solution for concentrated loads has been integrated over an area with sides using standard Gaussian quadrature in order to obtain results that could be used to validate the present solution for distributed dynamic loads.

Figure 6 shows the complete set of 9 dynamic displacement components obtained by integrating the solution presented by Kitahara [6] and the solution obtained by the present implementation. In the graphics the results are indicated, respectively, as (Kelvin) and (2-Fourier). The graphics show the real and imaginary parts of the components of displacement along a line in the interior of the full space given by and . The other parameters are , , , , and .

In Figure 7, the frequency-dependent behavior of the displacement solution is shown. The real and imaginary parts of the components and are determined for the point within the full space with coordinates . The dimensionless frequency varies from to . The other parameters are , and .

Now consider an increase in the loading dimension . If is large enough, the present 3D solution should approach the plane strain condition at the plane . Barros and De Mesquita Neto [7] determined the displacement solutions for 2D full spaces under distributed loads of spatially constant amplitude. Figures 8(a) and 8(b) show the displacements components and obtained for the 2D plane strain case [7] compared to the present 3D formulation for the case in which the ratio . The other parameters of the problem are (plane ), , and . It can be seen that the 3D solution tends to the 2D case for increasing values of the ratio.

The difference observed in Figure 8 between the two solutions comes from the fact that a 3D solution is being compared with one in which a plane strain condition is established (2D case).

The results for the static and the dynamic solutions determined by the present approach compare very well with the case of the Kelvin solution, with the integrated Kitahara’s closed-form dynamic solution, and with Barros’ 2D plane strain solution. It is understood that these comparisons, static and dynamic, for all components of the displacement tensor validate the present approach and implementation.

4. Numerical Results

In this section a series of numerical studies are presented which allow assessing the influence of the distinct parameters on the dynamic displacement solutions.

4.1. Influence of the Damping

Part I of this paper described how a model of hysteretic damping was introduced in the formulation according to the elastic-viscoelastic correspondence principle [2]. Figure 9 shows the influence of distinct damping factors on the amplitude of the displacement components. The vertical axis is shown in log scale. The other parameters of the problem are and ; , and .

It can be observed from Figure 9 that an increase in the internal damping coefficient will cause the amplitude to decay faster along the coordinate . This increase in amplitude decay of propagating waves in viscoelastic unbounded solids is in accordance with the literature [8].

The frequency behavior of the solution is similar. Higher damping coefficients will stiffen the domain and decrease the displacement amplitude compared to the ones with smaller damping coefficients. This behavior can be seen in Figures 10(a) and 10(b) for the components and . The other parameters used in Figure 10 are , and .

4.2. Long-Distance Behavior

For the present formulation to be used, for instance, as an auxiliary state for the indirect formulation of the boundary element method (BEM) [9], the displacement solutions must be able to be determined at fairly large distances from the loading area. This section shows a few examples of this case. Figure 11 shows the real part of two displacement components at a line along the axis with a low frequency . The parameters of this problem are and .

Figure 12 shows the real part of two displacement components at a line along the axis for the case of a higher frequency, . The parameters of this problem are: , and .

Analogously, Figure 13 shows the real part of two displacement components at a line along the axis for an even higher frequency, . The parameters of this problem are , and .

These results indicate that the present implementation is capable of determining displacement components at fairly large distances from the loading source, for low and relative high frequencies.

4.3. Higher-Frequency Behavior

Frequency-dependent solutions may be used to obtain transient solutions via Laplace or Fourier integral transforms, by relating the frequency parameter with the time variable. It is also well known that there is a relation between the highest frequency in which a solution may be determined and the smallest time step of the transient solutions [10]. Therefore, it is important to verify if the present implementation is able to determine displacement solutions at larger frequencies.

Figure 14 shows the frequency behavior for two displacement components for a dimensionless frequency parameter varying from . The absolute value of the displacement amplitude is shown. The components and were calculated for the parameters , and .

Figures 15(a) and 15(b) show, exemplarily, the frequency behavior of the components and for frequencies up to . The parameters are , and .

The presented results indicate that the synthesized and implemented frequency domains solutions may be used to obtain transient solutions via, for instance, the FFT algorithm. To exemplify this statement, consider a continuum with a shear wave velocity of  m/s and a load half-width  m. for these parameters, the solutions shown in Figures 14 and 15 would lead to transient solutions with minimum time steps of and , respectively. These time steps are sufficiently small for many transient applications.

4.4. Behavior of the Solution for Nearly Incompressible Continua

Figure 16 shows the frequency solution for the displacement components and for the case in which Poisson’s ratio approaches the incompressibility limit . Results for three distinct values of , namely, , and , are shown in the figure. The other parameters are , and . These results show that the present implementation remains stable while going to the incompressibility limit and may be used to determine solutions for the case of nearly incompressible media.

5. Concluding Remarks

In the first part of the present paper, the displacement response of a three-dimensional isotropic viscoelastic full space under stationary dynamic loading was synthesized using the double Fourier integral transform. Expressions for the displacements in the original physical domain were furnished as inverse double integral transforms. In the present paper, the strategy applied to evaluate numerically the inverse double Fourier integral transforms is described. The obtained results were validated by comparison with a static and a dynamic concentrated load solutions. All the 9 components of the displacement tensor have been presented and validated. Furthermore, the influence of the damping coefficient on the solutions was investigated. It has been shown that, for the present 3D solution, an increase in the ratio of the loading area is able to recover the 2D plane strain solution. The behavior of the displacement solutions at larger distances from the loading area have been presented, showing the stability and accuracy of the implemented numerical scheme. The behavior of the solutions for high frequencies is also investigated. Stable numerical results for very high frequencies were obtained, indicating that the present solutions may be used in conjunction with an FFT algorithm to obtain transient solutions. It is expected that the determined displacement solutions may be used as auxiliary states in formulations such as the indirect boundary element method.


In (2.1), the following notation is adopted:


The research leading to this paper has been supported by Fapesp, CNPq, Capes, and the UFMS-Foundation. This is gratefully acknowledged.