# Stationary Dynamic Stress Solutions for a Rectangular Load Applied within a 3D Viscoelastic Isotropic Full-Space

Authors: E. Romanini; J. Labaki; E. Mesquita; R. C. Silva
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This paper presents stress influence functions for uniformly distributed, time-harmonic rectangular loads within a three-dimensional, viscoelastic, isotropic full-space. The coupled differential equations relating displacements and stresses in the full-space are solved through double Fourier integral transforms in the wave number domain, in which they can be solved algebraically. The final stress fields are expressed in terms of double indefinite integrals arising from the Fourier transforms. The paper presents numerical schemes with which to integrate these functions accurately. The article presents numerical validation of the synthesized stress kernels and their behavior for high frequencies and large distances from the excitation source. The influence of damping ratio on the dynamic results is also investigated. This article is complementary to previous results of the authors in which the corresponding displacement solutions were derived. Stress influence functions, together with their displacement counterparts, are a fundamental part of many numerical methods of discretization such the boundary element method.


---

## Body

## 1. Introduction

The displacements and stress fields of an elastic full-space are coupled according to Navier-Cauchy’s differential equation. Solution to these equations may be obtained through integral transforms, in which the coupled equations, written in physical space coordinates, are transformed into the wave number space where they can be solved algebraically. Numerous authors have used this technique to obtain stress and displacement fields through Fourier [1, 2], Radon [3–5], and Hankel transforms [6], which are selected according to the coordinate system or type of loading involved in a particular form of the Navier-Cauchy equations. For an extensive literature review on this topic, please refer to Mesquita, Romanini, and Labaki [7].Influence functions—the response of a medium to nonsingular loadings—can be used directly to model loading problems such as pavements, foundations, and footings [8], as well as wave propagation through unbounded media [9], while adequately accounting for Sommerfeld’s radiation condition [10]. Moreover, they can be used as auxiliary states within the framework of the boundary element method (BEM) to model more sophisticated problems [11, 12].The authors of the present article have previously derived a solution for displacement influence functions for rectangular loads within a viscoelastic full-space [7]. The intricacies of the numerical computation of their final indefinite integral expressions were discussed by the authors in Labaki, Romanini, and Mesquita [13]. Despite their contribution to the understanding of the response of elastic media, solutions for more practical, arbitrarily shaped soil-foundation interaction problems through the BEM require the corresponding stress influence functions to be derived. This is the goal of the present article. This work presents stress influence functions for uniformly distributed rectangular loads within three-dimensional, viscoelastic, isotropic full-spaces. An analytical procedure involving double Fourier transforms is used to uncouple the differential equations of motion so that they can be solved algebraically. Closed-form solutions for arbitrary time-harmonic rectangular loadings are presented in the transformed domain, in terms of indefinite double integrals. The article discusses the characteristics of the corresponding integrands and numerical schemes with which to integrate the final equations. Original numerical results are shown for selected representative parameters.


## 2. Problem Statement

Consider an unbounded, three-dimensional, isotropic, viscoelastic full-space, with complex-valued Lamé’s constantsμ and λ, damping coefficient η, and mass density ρ, described in terms of a rectangular coordinate system (x,y,z). A hysteretic damping model is considered for this medium [14], and the complex-valued Lamé’s constants relate to their real counterparts μr and λr according to the elastic-viscoelastic correspondence principle [15]: μ=μr(1+iη) and λ=λr(1+iη), in which i=√-1. Vertical or horizontal time-harmonic loads of circular frequency ω are applied within a rectangular surface of sides 2A and 2B in the x- and y-directions, respectively, on the xy plate (z=0) (Figure 1). The loads are uniformly distributed within this area.Figure 1
External vertical and horizontal loadings applied within the full-space.
(a) (b) (c)In the absence of body forces, the displacement fieldU=U(x) due to these loadings in the frequency domain is described by(1)μ∇2U+λ+μ∇∇•U=-ω2ρUThe corresponding stress field σ is obtained from U according to(2)σij=λ∇•Uδij+μUi,j+Uj,iThe present problem consists of deriving a solution for σ.


## 3. Solution Strategy

The vector displacement fieldui(i=x,y,z) of a continuous, linear-elastic, unbounded full-space may be assumed to be continuously differentiable and vanish at spatial infinity, such that it can be resolved into the linear superposition of a curl-free vector field Δ and a divergence-free vector field Ω, according to the Helmholtz decomposition (Eringen and Suhubi, 1979)(3)ui=-1kL2Δ,i+2kS2eimnΩn,m,in which(4)kL2=ω2ρλ+2μ,kS2=ω2ρμare dilatational and distortional wave numbers, respectively. The corresponding stress field becomes, in terms of this decomposition,(5)σijμ=δij1-2n2n2Δ-2kL2Δ,ij+2kS2eiklΩl,kj+ejklΩk,li,in which δij is the Kroenecker delta and n2=k2L/k2S. Trial solutions may be established for Δ and Ω as(6)Δ1=A1kL2exp⁡αLz+iβx+γy,(7)Ωj1=Bj1kS2exp⁡αSz+iβx+γy,(8)Δ2=A2kL2exp⁡-αLz+iβx+γy,(9)Ωj2=Bj2kS2exp⁡-αSz+iβx+γy,in which A(m) and Bn(m)(m=1,2;n=1,2,3) are arbitrary functions corresponding to specific boundary-value problems; i.e., they will be determined by the boundary conditions of a specific problem. The index m indicates if the trial solution corresponds to Δ and Ωj within domain m=1(-∞<z≤0) or domain m=2(0≤z<+∞). It is necessary to split the trial functions into two expressions, corresponding to the two domains, m=1(-∞<z≤0) and domain m=2(0≤z<+∞), so that the arguments αLz and αSz of the exponential terms vanish at infinity, as required by Sommerfeld’s radiation condition [10].The conditions thatΔ(m) and Bn(m) must be, respectively, curl-free and divergence-free yield(10)αL,S2=β2+γ2-kL,S2,(11)B31,2=∓iαSβB1+γB2.The substitution of (6) to (9) into (3) and (5) yields the expressions of displacement and stress fields in the wave number domain β,γ for media 1 and 2 in terms of the trial functions:(12)uX1,2=-A1,2iβe±αLz±2αSB11,2βγ+B21,2γ2-αS2e±αSzeiβx+γy,(13)uY1,2=-A1,2iγe±αLz±2αSB11,2-β2+αS2-B21,2βγe±αSzeiβx+γy,(14)uZ1,2=∓A1,2αLe±αLz-2iB11,2γ-B21,2βe±αSzeiβx+γy,(15)σXX1,2=μA1,2β2-γ2-αS2+2αL2e±αLz±4iβαSB11,2βγ+B21,2γ2-αS2e±αSzeiβx+γy(16)σXY1,2=2μA1,2βγe±αLz±iαSB11,2γ2β+βαS2-β3+B21,2γ3-γαS2-β2γe±αSzeiβx+γy(17)σXZ1,2=2μ∓iA1,2βαLe±αLz+2B11,2βγ+B21,2γ2-αS2-β2e±αSzeiβx+γy(18)σYY1,2=μA1,2γ2-β2-αS2+2αL2e±αLz±4iγαSB11,2αS2-β2-B21,2βγe±αSzeiβx+γy(19)σYZ1,2=2μ∓iA1,2γαLe±αLz+B11,2αS2+γ2-β2-2B21,2βγe±αSzeiβx+γy,(20)σZZ1,2=μ-A1,2γ2+β2+αS2e±αLz∓4iαSB11,2γ-B21,2βe±αSzeiβx+γy.The six remaining arbitrary functions A(m) and Bn(m)(m,n=1,2) in (12) to (20) arise upon establishing both kinematic compatibilities,(21)uj1x,y,z=0=uj2x,y,z=0,and equilibrium conditions(22)σZj1x,y,z=0-σZj2x,y,z=0=pjβ,γ,at the loaded interface between media 1 and 2 (z=0). The solution of these six algebraic equations summarized in (21) and (22) (j=x,y,z) results in A(m) and Bn(m) corresponding to the full-space boundary-value problem, in which(23)pjβ,γ=p-jβ,γeiβx+γyis the time-harmonic external loads of circular frequency ω and amplitude p-j at the interface between media 1 and 2, presented in the wavenumber domain β,γ. For this boundary-value problem, the solution for the six trial functions is given by(24)A1,2=12iβp-Xβ,γαLμ∗kS∗2+12γp-Yβ,γαLμ∗kS∗2±12p-Zβ,γμ∗kS∗2,(25)B11,2=∓14p-Yβ,γμ∗kS∗2+14iγp-Zβ,γμ∗αSkS∗2,(26)B21,2=±14p-Xβ,γμ∗kS∗2-14iβp-Zβ,γμ∗αSkS∗2.


### 3.1. Distributed Loads

The nonsingular loading case proposed in this article may be obtained by integrating the external point loadpj(β,γ) within the rectangular patch (-A≤x≤A, -B≤y≤B, z=0) where it is applied. An expression for this loading in the wavenumber domain has been derived by Mesquita, Romanini, and Labaki [7]:(27)p¯jβ,γ=-4Pjβγsin⁡βAsin⁡γB,j=x,y,z,where Pj is the spatially constant amplitude of this uniformly distributed load.


### 3.2. Stress Fields in the Physical Domain

The substitution of (24) to (26) into (15) to (20) yields the stress fields in terms of the wavenumbers kS and kL. Their physical domain counterparts are obtained upon applying an inverse double Fourier transformation to those expressions:(28)σXX=PXμπ2∫0∞∫0∞αSβ2-γ2+2αL2-αS2e-αLz+2αLγ2-αS2e-αSzαSαLsβsβxdkβsγcγykγdkγ+PYμπ2∫0∞∫0∞αSβ2-γ2+2αL2-αS2e-αLz-2αLβ2e-αSzαSαLsβcβxkβdkβsγsγydkγ+zzPZμπ2∫0∞∫0∞β2-γ2+2αL2-αS2e-αLz+2β2e-αSzαSαLsβcβxkβdkβsγcγydkγ(29)σXY=PXμπ2∫0∞∫0∞2β2αSe-αLz-αL2β2-kS2e-αSzαSαLsβcβxkβdkβsγsγydkγ+PYμπ2∫0∞∫0∞2γ2αSe-αLz-αL2γ2-kS2e-αSzαSαLsβsβxdkβsγcγykγdkγ-2zzPZμπ2∫0∞∫0∞αSαLe-αLz-αSαLe-αSzαSαLsβsβxdkβsγsγydkγ(30)σXZ=zzPXμπ2∫0∞∫0∞2β2αSαLe-αLz+αSαLγ2-β2-αS2e-αSzαSαLsβcβxkβdkβsγcγykγdkγ+2zzPYμπ2∫0∞∫0∞αSαLe-αLz-αSαLe-αSzαSαLsβsβxdkβsγsγydkγ-PZμπ2∫0∞∫0∞2αSαL2e-αLz-αLγ2+β2+αS2e-αSzαSαLsβsβxdkβsγcγykγdkγ(31)σYY=-PXμπ2∫0∞∫0∞αSβ2-γ2-2αL2+αS2e-αLz+2αLγ2e-αSzαSαLsβsβxdkβsγcγykγdkγ+PYμπ2∫0∞∫0∞αSβ2-γ2-2αL2+αS2e-αLz-2αLβ2-αS2e-αSzαSαLsβcβxkβdkβsγsγydkγ-zzPZμπ2∫0∞∫0∞αSαLβ2-γ2-2αL2+αS2e-αLz+2αSαLγ2e-αSzαSαLsβcβxkβdkβsγcγykγdkγ(32)σYZ=2zzPXμπ2∫0∞∫0∞αSαLe-αLz-αSαLe-αSzαSαLsβsβxdkβsγsγydkγ+zzPYμπ2∫0∞∫0∞2γ2αSαLe-αLz-αSαLγ2-β2+αS2e-αSzαSαLsβcβxkβdkβsγcγykγdkγ-PZμπ2∫0∞∫0∞2αSαL2e-αLz-αLγ2+β2+αS2e-αSzαSαLsβcβxkβdkβsγsγydkγ(33)σZZ=-PXμπ2∫0∞∫0∞αSβ2+γ2+αL2e-αLz-2αS2αLe-αSzαSαLsβsβxdkβsγcγykγdkγ-PYμπ2∫0∞∫0∞αSβ2+γ2+αS2e-αLz-2αS2αLe-αSzαSαLsβcβxkβdkβsγsγydkγ-zzPZμπ2∫0∞∫0∞αSαLβ2+γ2+αS2e-αLz-2αSαLβ2+γ2e-αSzαSαLsβcβxkβdkβsγcγykγdkγin which(34)β=a0Akβ,γ=a0Akγ,(35)kS∗2=a0A211+iμ,kL∗2=a0A2n21+iμ,(36)sβ=sin⁡a0kβ,sβx=sina0kβxA,sγ=sina0b0kγ,sγy=sina0b0kγyB,(37)cβ=cosa0kβ,cβx=cosa0kβxA,cγ=cosa0b0kγ,cγy=cosa0b0kγyB,(38)n2=kL2kS2,(39)αL2=β2+γ2-kL2,αS2=β2+γ2-kS2,(40)a0=ωAμρ,in which a0 is a normalized frequency of excitation. Note in (28) to (33) that unified expressions encompassing both media 1 and 2 can be obtained by adjusting the expressions by a factor z/|z|.


### 3.3. Numerical Evaluation of the Stress Fields

In order to obtain the stress fields of the full-space in the physical domain, the double indefinite integrals expressing these fields (see (28) to (33)) must be evaluated numerically. The evaluation of these integrals requires a careful selection of appropriate numerical schemes. The integrands of the stress fields are characterized by the superposition of a singular kernel and an oscillatory-decaying kernel.The singular kernel presents two singularities corresponding to the pressure and shear waves in the full-space (Labaki, Mesquita, and Romanini, 2012), which are confined within a predictable region. A monotonic decay of the singular kernel is observed beyond this region. In order to deal with these singularities, a small damping factorη has been incorporated to the elastic constants according to the elastic-viscoelastic correspondence principle [15]. This strategy smooths out the singularities in that kernel and enables the integral within the singular region to be evaluated with ordinary adaptive Gaussian quadratures (Piessens, Doncker-Kapenga, and Überhuber, 1983). The resulting integral then corresponds to a viscoelastic full-space, rather than an elastic one. This is a reasonable model for many applications of these influence functions.Beyond the singular region, the behavior of the integrand is dominated by the oscillatory-decaying kernel. In this second region, the integrand oscillates about zero and decays in amplitude indefinitely. A method to approximate indefinite integrals of such functions has been proposed by Longman [16] in terms of Euler’s expansion of convergent series. The present implementation uses Longman’s method to deal with the oscillatory-decaying component of the integrands. For a more complete description of the integration scheme and the division of the integrand into singular and oscillatory-decaying regions, please refer to (Labaki, Mesquita, and Romanini, 2012).


## 3.1. Distributed Loads

The nonsingular loading case proposed in this article may be obtained by integrating the external point loadpj(β,γ) within the rectangular patch (-A≤x≤A, -B≤y≤B, z=0) where it is applied. An expression for this loading in the wavenumber domain has been derived by Mesquita, Romanini, and Labaki [7]:(27)p¯jβ,γ=-4Pjβγsin⁡βAsin⁡γB,j=x,y,z,where Pj is the spatially constant amplitude of this uniformly distributed load.


## 3.2. Stress Fields in the Physical Domain

The substitution of (24) to (26) into (15) to (20) yields the stress fields in terms of the wavenumbers kS and kL. Their physical domain counterparts are obtained upon applying an inverse double Fourier transformation to those expressions:(28)σXX=PXμπ2∫0∞∫0∞αSβ2-γ2+2αL2-αS2e-αLz+2αLγ2-αS2e-αSzαSαLsβsβxdkβsγcγykγdkγ+PYμπ2∫0∞∫0∞αSβ2-γ2+2αL2-αS2e-αLz-2αLβ2e-αSzαSαLsβcβxkβdkβsγsγydkγ+zzPZμπ2∫0∞∫0∞β2-γ2+2αL2-αS2e-αLz+2β2e-αSzαSαLsβcβxkβdkβsγcγydkγ(29)σXY=PXμπ2∫0∞∫0∞2β2αSe-αLz-αL2β2-kS2e-αSzαSαLsβcβxkβdkβsγsγydkγ+PYμπ2∫0∞∫0∞2γ2αSe-αLz-αL2γ2-kS2e-αSzαSαLsβsβxdkβsγcγykγdkγ-2zzPZμπ2∫0∞∫0∞αSαLe-αLz-αSαLe-αSzαSαLsβsβxdkβsγsγydkγ(30)σXZ=zzPXμπ2∫0∞∫0∞2β2αSαLe-αLz+αSαLγ2-β2-αS2e-αSzαSαLsβcβxkβdkβsγcγykγdkγ+2zzPYμπ2∫0∞∫0∞αSαLe-αLz-αSαLe-αSzαSαLsβsβxdkβsγsγydkγ-PZμπ2∫0∞∫0∞2αSαL2e-αLz-αLγ2+β2+αS2e-αSzαSαLsβsβxdkβsγcγykγdkγ(31)σYY=-PXμπ2∫0∞∫0∞αSβ2-γ2-2αL2+αS2e-αLz+2αLγ2e-αSzαSαLsβsβxdkβsγcγykγdkγ+PYμπ2∫0∞∫0∞αSβ2-γ2-2αL2+αS2e-αLz-2αLβ2-αS2e-αSzαSαLsβcβxkβdkβsγsγydkγ-zzPZμπ2∫0∞∫0∞αSαLβ2-γ2-2αL2+αS2e-αLz+2αSαLγ2e-αSzαSαLsβcβxkβdkβsγcγykγdkγ(32)σYZ=2zzPXμπ2∫0∞∫0∞αSαLe-αLz-αSαLe-αSzαSαLsβsβxdkβsγsγydkγ+zzPYμπ2∫0∞∫0∞2γ2αSαLe-αLz-αSαLγ2-β2+αS2e-αSzαSαLsβcβxkβdkβsγcγykγdkγ-PZμπ2∫0∞∫0∞2αSαL2e-αLz-αLγ2+β2+αS2e-αSzαSαLsβcβxkβdkβsγsγydkγ(33)σZZ=-PXμπ2∫0∞∫0∞αSβ2+γ2+αL2e-αLz-2αS2αLe-αSzαSαLsβsβxdkβsγcγykγdkγ-PYμπ2∫0∞∫0∞αSβ2+γ2+αS2e-αLz-2αS2αLe-αSzαSαLsβcβxkβdkβsγsγydkγ-zzPZμπ2∫0∞∫0∞αSαLβ2+γ2+αS2e-αLz-2αSαLβ2+γ2e-αSzαSαLsβcβxkβdkβsγcγykγdkγin which(34)β=a0Akβ,γ=a0Akγ,(35)kS∗2=a0A211+iμ,kL∗2=a0A2n21+iμ,(36)sβ=sin⁡a0kβ,sβx=sina0kβxA,sγ=sina0b0kγ,sγy=sina0b0kγyB,(37)cβ=cosa0kβ,cβx=cosa0kβxA,cγ=cosa0b0kγ,cγy=cosa0b0kγyB,(38)n2=kL2kS2,(39)αL2=β2+γ2-kL2,αS2=β2+γ2-kS2,(40)a0=ωAμρ,in which a0 is a normalized frequency of excitation. Note in (28) to (33) that unified expressions encompassing both media 1 and 2 can be obtained by adjusting the expressions by a factor z/|z|.


## 3.3. Numerical Evaluation of the Stress Fields

In order to obtain the stress fields of the full-space in the physical domain, the double indefinite integrals expressing these fields (see (28) to (33)) must be evaluated numerically. The evaluation of these integrals requires a careful selection of appropriate numerical schemes. The integrands of the stress fields are characterized by the superposition of a singular kernel and an oscillatory-decaying kernel.The singular kernel presents two singularities corresponding to the pressure and shear waves in the full-space (Labaki, Mesquita, and Romanini, 2012), which are confined within a predictable region. A monotonic decay of the singular kernel is observed beyond this region. In order to deal with these singularities, a small damping factorη has been incorporated to the elastic constants according to the elastic-viscoelastic correspondence principle [15]. This strategy smooths out the singularities in that kernel and enables the integral within the singular region to be evaluated with ordinary adaptive Gaussian quadratures (Piessens, Doncker-Kapenga, and Überhuber, 1983). The resulting integral then corresponds to a viscoelastic full-space, rather than an elastic one. This is a reasonable model for many applications of these influence functions.Beyond the singular region, the behavior of the integrand is dominated by the oscillatory-decaying kernel. In this second region, the integrand oscillates about zero and decays in amplitude indefinitely. A method to approximate indefinite integrals of such functions has been proposed by Longman [16] in terms of Euler’s expansion of convergent series. The present implementation uses Longman’s method to deal with the oscillatory-decaying component of the integrands. For a more complete description of the integration scheme and the division of the integrand into singular and oscillatory-decaying regions, please refer to (Labaki, Mesquita, and Romanini, 2012).


## 4. Validation and Numerical Results

Figures2–4 show comparisons between the results from the present implementation and cases from the literature. Figure 2 considers the case of a static (a0=0), concentrated load within an elastic, three-dimensional full-space. The case of concentrated load can be approached with the present implementation by making A and B small. In these results, A=B=0.01. The results consider a Poisson ratio ν=0.25, coordinates y=0.3 and z=0.5, and show selected stress components along a line 0.1≤x≤0.5. The results agree with the classical Kelvin solution for this problem [17].Figure 2
Comparison of the present solution for static, concentrated loads.
(a) (b)Figure 3
Comparison of the present solution for dynamic loads: vertical loads.
(a) (b)Figure 4
Comparison of the present solution for dynamic loads: horizontal loads.
(a) (b)The reference solution for the case of dynamic, distributed loads has been presented by Barros and Mesquita [18], which considers a two-dimensional full-space. The two-dimensional problem can be approached with the present implementation by making the ratio B/A large. In these results, A=1 and B=50. The results consider ν=0.25, x=0.8, y=0, and z=0.5, which correspond to a point on the plane x-z, for the sake of comparison with the 2D solution. The results are compared for a frequency range 0.5≤a0≤2.5.


### 4.1. Boundary Conditions

The aim of this section is to show that the present implementation satisfies the boundary conditions of the problem (see (22)). Figure 5 shows two representative components of the stress fields near the loaded surface of the full-space (|x|≤A; |y|≤B; z=0). When a unit load is applied at that surface, these components must satisfy σZZZ=σXZX=∓0.5 for z=0±; |x|≤A; |y|≤B (under the loaded surface), and σZZZ=σXZX=0 otherwise. The stress fields near the loaded surface are particularly difficult to compute. Equations (28) to (33) show that the decaying component of the integrands vanish for very small values of z. Longman’s series extrapolation integration method, currently being used to compute the oscillatory-decaying portion of that integrand, is poorly suited to deal with slowly convergent series. Nevertheless, Figure 5 shows that the stress components tend monotonically to the boundary conditions as the depth z approaches the loaded surface. For these results, -2≤x≤2, y=0.5, A=B=1, and a0=0.Figure 5
Verification of the boundary conditions in the present implementation.
(a) (b)


### 4.2. Influence of the Damping Factor

Figures6–9 show how the behavior of selected stress components is affected by different damping ratios. In these results, A=B=1, ν=0.25, y=z=0.5, a0=0.5, and 0≤x≤10.Figure 6
Influence of the damping factor: vertical loads.
(a) (b)Figure 7
Influence of the damping factor: vertical loads.
(a) (b)Figure 8
Influence of the damping factor: horizontal loads.
(a) (b)Figure 9
Influence of the damping factor: horizontal loads.
(a) (b)Figures6–9 show that the effect of the damping coefficient in the real part of the stress components is negligible, while a significant influence is observed in the imaginary parts. The way in which damping ratios affect these imaginary components at various distances from the excitation source is different for each of the stress components shown. A small damping coefficient is introduced in the elastic constants in the present implementation in order to facilitate the numerical integration of (28) to (33). It is important to keep in mind the effect that different values of this coefficient will have on the solution.Additionally, Figures10 and 11 present the effect of damping ratios in the stress fields for various frequencies of excitation. For this analysis, x=y=z=2, A=B=1, and ν=0.25. The results show that the stress amplitude decreases more quickly for larger damping coefficients and frequencies of excitation. This increase in amplitude decay with larger damping coefficients is physically consistent [14].Figure 10
Stress solutions distant from the excitation source: vertical loads.
(a) (b)Figure 11
Stress solutions distant from the excitation source: horizontal loads.
(a) (b)


### 4.3. Long-Distance Behavior

One of the applications of influence functions such as the one presented in this article is as fundamental solutions in the analysis of problems through the boundary element method (BEM). In the solution of typical problems with the BEM, these influence functions will have to be evaluated at points that are far from the excitation source. This section presents representative cases in which the stress fields are computed in such situations. Figures12 and 13 show selected stress components for three different frequencies of excitation, and y=z=0.5, A=B=1, η=0.05, ν=0.25, and 5≤x≤50. These qualitative results show that the present implementation is capable of computing stress fields at large distances from the excitation source.Figure 12
Stress solutions distant from the excitation source: vertical loads.
(a) (b)Figure 13
Stress solutions distant from the excitation source: horizontal loads.
(a) (b)


### 4.4. Higher-Frequency Behavior

Even though this work presents the time-harmonic solution of stress fields, the corresponding transient response may be obtained from these solutions through Fourier or Laplace transforms. In order for these transforms to be accomplished successfully, it is necessary that the time-harmonic solutions be evaluated at considerably high frequencies. The smallest time step in the resulting transient solution is inversely proportional to the highest frequency in the original time-harmonic signal. Figures14 and 15 show that the present solution is capable of computing time-harmonic stress fields for considerably high frequencies of excitation. In these results, x=y=z=2.0, A=B=1, ν=0.25, η=0.01, and 0≤a0≤20.Figure 14
Higher-frequency behavior: vertical loads.Figure 15
Higher-frequency behavior: horizontal loads.


## 4.1. Boundary Conditions

The aim of this section is to show that the present implementation satisfies the boundary conditions of the problem (see (22)). Figure 5 shows two representative components of the stress fields near the loaded surface of the full-space (|x|≤A; |y|≤B; z=0). When a unit load is applied at that surface, these components must satisfy σZZZ=σXZX=∓0.5 for z=0±; |x|≤A; |y|≤B (under the loaded surface), and σZZZ=σXZX=0 otherwise. The stress fields near the loaded surface are particularly difficult to compute. Equations (28) to (33) show that the decaying component of the integrands vanish for very small values of z. Longman’s series extrapolation integration method, currently being used to compute the oscillatory-decaying portion of that integrand, is poorly suited to deal with slowly convergent series. Nevertheless, Figure 5 shows that the stress components tend monotonically to the boundary conditions as the depth z approaches the loaded surface. For these results, -2≤x≤2, y=0.5, A=B=1, and a0=0.Figure 5
Verification of the boundary conditions in the present implementation.
(a) (b)


## 4.2. Influence of the Damping Factor

Figures6–9 show how the behavior of selected stress components is affected by different damping ratios. In these results, A=B=1, ν=0.25, y=z=0.5, a0=0.5, and 0≤x≤10.Figure 6
Influence of the damping factor: vertical loads.
(a) (b)Figure 7
Influence of the damping factor: vertical loads.
(a) (b)Figure 8
Influence of the damping factor: horizontal loads.
(a) (b)Figure 9
Influence of the damping factor: horizontal loads.
(a) (b)Figures6–9 show that the effect of the damping coefficient in the real part of the stress components is negligible, while a significant influence is observed in the imaginary parts. The way in which damping ratios affect these imaginary components at various distances from the excitation source is different for each of the stress components shown. A small damping coefficient is introduced in the elastic constants in the present implementation in order to facilitate the numerical integration of (28) to (33). It is important to keep in mind the effect that different values of this coefficient will have on the solution.Additionally, Figures10 and 11 present the effect of damping ratios in the stress fields for various frequencies of excitation. For this analysis, x=y=z=2, A=B=1, and ν=0.25. The results show that the stress amplitude decreases more quickly for larger damping coefficients and frequencies of excitation. This increase in amplitude decay with larger damping coefficients is physically consistent [14].Figure 10
Stress solutions distant from the excitation source: vertical loads.
(a) (b)Figure 11
Stress solutions distant from the excitation source: horizontal loads.
(a) (b)


## 4.3. Long-Distance Behavior

One of the applications of influence functions such as the one presented in this article is as fundamental solutions in the analysis of problems through the boundary element method (BEM). In the solution of typical problems with the BEM, these influence functions will have to be evaluated at points that are far from the excitation source. This section presents representative cases in which the stress fields are computed in such situations. Figures12 and 13 show selected stress components for three different frequencies of excitation, and y=z=0.5, A=B=1, η=0.05, ν=0.25, and 5≤x≤50. These qualitative results show that the present implementation is capable of computing stress fields at large distances from the excitation source.Figure 12
Stress solutions distant from the excitation source: vertical loads.
(a) (b)Figure 13
Stress solutions distant from the excitation source: horizontal loads.
(a) (b)


## 4.4. Higher-Frequency Behavior

Even though this work presents the time-harmonic solution of stress fields, the corresponding transient response may be obtained from these solutions through Fourier or Laplace transforms. In order for these transforms to be accomplished successfully, it is necessary that the time-harmonic solutions be evaluated at considerably high frequencies. The smallest time step in the resulting transient solution is inversely proportional to the highest frequency in the original time-harmonic signal. Figures14 and 15 show that the present solution is capable of computing time-harmonic stress fields for considerably high frequencies of excitation. In these results, x=y=z=2.0, A=B=1, ν=0.25, η=0.01, and 0≤a0≤20.Figure 14
Higher-frequency behavior: vertical loads.Figure 15
Higher-frequency behavior: horizontal loads.


## 5. Concluding Remarks

This paper presented a derivation of stress fields within homogeneous, three-dimensional, viscoelastic, isotropic full-spaces due to time-harmonic, uniformly distributed rectangular loads. The derivation is obtained upon the decomposition of the displacement fields in the full-space in terms of a curl-free and a divergence-free vector fields, and their subsequent double Fourier space transform. The resulting transformation yielded algebraic expression for the displacement and stress fields within the full-space, in which the boundary conditions were applied. The resulting stress fields were expressed in terms of double indefinite integrals that must be computed numerically for each analysis. The results showed that the present implementation produces physically consistent results and is robust for a wide range of parameters. The present stress field solution can be used together with the displacement solutions presented previously by the authors in this journal to solve a variety of elastodynamic problems through the boundary element and other meshless methods.



---

*Source: 4738498-2019-03-26.xml*