Stationary Dynamic Displacement Solutions for a Rectangular Load Applied within a 3 D Viscoelastic Isotropic Full Space — Part II : Implementation , Validation , and Numerical Results

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.


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 fullspace under spatially rectangular loads of constant amplitude, harmonically varying in time are furnished.

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 2A × 2B see Appendix for notation , U ZZ : Equation 2.1 presents two improper integrations, to be performed over the wave number variables k β and k γ .An outer improper integration is performed over the variable k γ .For every value of k γ , an improper inner integration over the variable k β 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 U ij i, j x, y, z furnished in part I 1 .These integrands are composed by a kernel K ij i, j x, y, z multiplied by oscillating functions c i , c ij , s i and s ij i β, γ , j x, y .For the elastic case, these kernels, here represented by K ZZ 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 K ZZ as a function of the integration variable k β for the following set of parameters: Lamé constant μ 1 and Poisson's ratio ν 0.25, loading half-width A B 1, nondimensional frequency A 0 1, and damping coefficient ranging from η 0 elastic case to η 0.2.It can be shown that, for values of the integration variable k β 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 k β 1 and k β ≈ 0.55, 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 η / 0 and the elastic constants are made complex, according to the elastic-viscoelastic correspondence principle 1, 2 .As a consequence, the singularities in the kernels K ij are smoothed and no integration over a singularity must be evaluated 3 .Now consider the improper inner integration over k β , described in 2.1 .The difference between the behavior of the full integrand and that of its kernel K ZZ 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 k β ≈ 0.55 and k β 1.
Figures 1 and 2 characterize the line of inner improper integration over k β for the value k γ 0 of the outer improper integration variable.The integrands present the same behavior  with respect to the variable k γ .In Figure 3, the behavior of the kernel K ZZ is shown for the plane of double-integration variables k β and k γ .Notice that, along the planes k β 0 and k γ 0, the integration kernels reproduce the behavior shown in Figure 1.

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 k β , k γ 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 0 ≤ k β ≤ k β0 and 0 ≤ k γ ≤ k γ0 .Actually it can be shown that the singularities in the plane k β , k γ are to be found at the locations determined by the equations k 2 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 0 ≤ k γ ≤ k γ0 and improper in the other direction k β0 ≤ k β < ∞ .There is no singular or quasi-singular behavior of the integrands in this region.There is an oscillatory behavior in the range 0 ≤ k γ ≤ k γ0 and a monotonic decreasing oscillatory behavior in the direction of the improper integration k β0 ≤ k β < ∞ .In the k γ direction, standard 1D Gaussian quadrature is employed.For the k β 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 k γ0 ≤ k γ ≤ ∞; k β0 ≤ k β ≤ ∞ , the integrands oscillate and decay monotonically in both directions, and Longman's integration strategy is used in these two directions.

Validation
The methodology described in the previous section was used to determine the displacements U ij i, j x, y, z within the full space due to a dynamic load uniformly distributed over an area of sides 2A × 2B.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, x, y, and z, on the dimensionless frequency A 0 , on the constitutive parameters μ and λ, on the dimensions of the loading area A and B, and also on the internal damping coefficient η.
The present solution can be partially validated using the limiting case of a static solution A 0 → 0 due to a concentrated point load A → 0, B → 0 .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 A 0 0.001 and the half-widths of the loading were set to A B 0.01.The other parameters are ν 0.25, μ 1, and η 0.01.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 U ij i, j x, y, z are determined.The solutions were obtained along the line of the full space given by y z 0.5 and −2 ≤ x ≤ 2.
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 2A × 2B 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 U ij i, j x, y, z 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 y z 0.5 and −2 ≤ x ≤ 2. The other parameters are ν 0.25, μ 1, η 0.01, A B 1, and A 0 5.In Figure 7, the frequency-dependent behavior of the displacement solution is shown.The real and imaginary parts of the components U XY and U ZZ are determined for the point within the full space with coordinates x y z 0.5.The dimensionless frequency varies from A 0 0 to A 0 4. The other parameters are μ 1.0, ν 0.25, η 0.01, and A B 1.
Now consider an increase in the loading dimension B. If B is large enough, the present 3D solution should approach the plane strain condition at the plane x-z, y 0 .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 U XX and U ZX obtained for the 2D plane strain case 7 compared to the present 3D formulation for the case in which the ratio B/A 10.The other parameters of the problem are −4 ≤ x ≤ 4, y 0 plane x-z , z 0.5, A 0 0.1, A 1, ν 0.25, μ 1, and η 0.1.It can be seen that the 3D solution tends to the 2D case for increasing values of the A/B 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 U ij i, j x, y, z validate the present approach and implementation.

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.

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  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 x.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

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 x axis with a low frequency A 0 0.1.The parameters of this problem are 0 ≤ x ≤ 50, y z 0.5, A B 1, ν 0.25, μ 1 and η 0.05.  Figure 12 shows the real part of two displacement components at a line along the x axis for the case of a higher frequency, A 0 1.0.The parameters of this problem are: 0 ≤ x ≤ 40, y z 0.5, A B 1, ν 0.25, μ 1, and η 0.05.
Analogously, Figure 13 shows the real part of two displacement components at a line along the z axis for an even higher frequency, A 0 2.0.The parameters of this problem are 0 ≤ z ≤ 50, x y 2.0, A B 1, ν 0.25, μ 1, and η 0.05.
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.

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 0 < A 0 < 20.The absolute value of the displacement amplitude is shown.The components U XZ and U Y Y were calculated for the parameters x y z 2.0, A B 1, ν 0.25, μ 1, and η 0.01.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 c S 250 m/s and a load half-width A 1 m. for these parameters, the solutions shown in Figures 14 and 15 would lead to transient solutions with minimum time steps of Δt min 1.26 • 10 −3 s and Δt min 0.5 • 10 −3 s, respectively.These time steps are sufficiently small for many transient applications.

Behavior of the Solution for Nearly Incompressible Continua
Figure 16 shows the frequency solution for the displacement components U XY A 0 and U Y Y A 0 for the case in which Poisson's ratio approaches the incompressibility limit

Appendix
In 2.1 , the following notation is adopted: A.1

Figure 1 :
Figure 1: a Behavior of the typical integrand kernel K ZZ as a function of the normalized wave number k β in the range 0 ≤ k β ≤ 1.5, k γ 0. b Detail of the singularity due to the dilatational wavefront.In the range 0.5 ≤ k β ≤ 0.65, k γ 0.

Figure 2 :
Figure 2: Behavior of the typical integrand with its oscillatory part, as a function of the normalized wave number k β in the range 0 ≤ k β ≤ 4, for k γ 0, A 0 3.

Figure 4 :
Figure 4: Regions of integration for the double inverse Fourier transform.

3 .
In the present implementation the values k β0 k γ0 1.2 have been chosen.It can also be shown that the oscillatory character of the integrand increases for larger values of the distance x, y, and z, as well as with the increase of the nondimensional frequency parameter A 0 .

Figure 5 :
Figure 5: Static displacements solutions due to a concentrated point load within the full-space along the line −2 ≤ x ≤ 2, y z 0.5.

Figure 8 :
Figure 8: Reproducing the 2D case by increasing the ratio B/A.

Figure 9 :
Figure 9: Effect of the damping factor on the amplitude of the displacement components.

Figure 10 :
Figure 10: Effect of the damping factor on the frequency behavior of the viscoelastic media.

Figure 11 :
Figure 11: Displacements solutions far from the origin for a low frequency A 0 0.1.

Figure 12 :Figure 13 :
Figure 12: Displacements solutions far from the origin for a higher frequency, A 0 1.0.

Figure 14 :
Figure 14: Frequency-dependent behavior of the displacement solutions.

Figure 15 :
Figure 15: Behavior of the displacement solutions for higher frequencies.
sen A 0 k β , s βx sen A 0 k β x A , s γ sen A 0 b 0 k γ , s γy sen A 0 b 0 k γ y B , c β cos A 0 k β , c βx cos A 0 k β x A , c γ cos A 0 b 0 k γ , c γy cos A 0 b 0 k γ y B .