Stationary Dynamic Displacement Solutions for a Rectangular Load Applied within a 3D Viscoelastic Isotropic Full Space—Part I: Formulation

Abstract

A dynamic stationary semianalytical solution for a spatially constant load applied over a rectangular surface within a viscoelastic isotropic full space is presented. The solution is obtained within the frame of a double Fourier integral transform. Closed-form solutions for general loadings within the full space are furnished in the transformed wave number domain. Expressions for three boundary value problems, associated to a normal and two tangential rectangular loadings in the original physical space, are given in terms of a double inverse Fourier integral. These inverse integral transforms must be evaluated numerically. In the second part of the present paper a strategy to evaluate these integrals is described, the procedure validated and a number of original results are reported.

1. Introduction

The dynamics of unbounded domains is characterized by outgoing and nonreflected waves which carry energy away from the perturbation source. The withdrawal of energy in the form of nonreflected outgoing waves, radiating waves, introduces in the unbounded medium a damping mechanism called geometric damping. The mathematical expression governing this condition is known as Sommerfeld radiation condition (SRC) [1, 2]. Many relevant geomechanical problems deal with the modeling of engineering structures interacting with the soil. A proper modeling of the dynamic soil-structure interaction problem (DSSI) must take the Sommerfeld radiation condition into account.

This has been recognized by early researchers like Reissner [3], who developed a semianalytical solution for a harmonic vibrating circular load applied at the surface of an elastic half-space. Later on, Lysmer [4] used Reissner’s solution within a superposition scheme to describe the interaction of rigid circular surface foundation interacting with the soil. Based on the semianalytical solution of a rectangular load applied at the surface of the half-space Wong and Luco [5] devised a scheme to model the interaction of surface foundation of arbitrary shape also with the elastic soil.

The synthesis of a semianalytical solution for a rectangular load of constant amplitude over a 3D viscoelastic-half space was determined by Gaul in his dissertation [6]. The work of Gaul used a double Fourier integral transform as the mathematical framework. The final solution was obtained through the numerical realization of inverse integral transforms. Gaul applied the obtained solution to analyze surface foundation interacting with the homogeneous viscoelastic half-space [7]. Half-space semianalytical solutions were used by Dasgupta [8] in conjunction with the Finite Element Method to model the interaction of rigid foundations embedded in the elastic half-space. A hybrid approach for 3D analysis based on the distributed load solutions has been proposed by Mita and Luco [9].

The solution of a 2D transverse isotropic full space subjected to a concentrated load and also to a spatially distributed load of constant amplitude and finite width under harmonic time behavior was determined by Rajapakse and Wang [10, 11]. The axisymmetric case was also solved by Wang [12] and used within the context of the Indirect Boundary Element Method (IBEM) to describe the dynamic response of buried structures.

Distributed load solutions for a transverse isotropic media with inclined principal axis were furnished by Barros and De Mesquita Neto [13]. Similar solutions for spatially constant distributed loads of finite width applied inside a 2D transverse isotropic were used by Barros and De Mesquita Neto [14] to analyzed the structure-soil-structure interaction with the IBEM. The solution for a distributed load with a singularity at one extremity was determined by Barros and De Mesquita Neto [15]. This solution was used to obtain a very efficient and accurate scheme to analyze the interaction of 2D rigid foundations with a half-space. A semianalytical solution for linear porous saturated full space was synthesized by Senjuntichai and Rajapakse [16].

In the previously cited work the mathematical framework to obtain these semianalytical solutions was the single and double Fourier integral transforms for the 2D and 3D cases, respectively. For the axisymmetric case the Hankel transform has been used.

In the last years the Radon transform has been applied to synthesize fundamental (point load) solutions for anisotropic full spaces [17, 18]. Using the Fourier and the Radon integral transforms, semianalytical solutions for concentrated and distributed 3D harmonic loads on the half-space surface were obtained, respectively, by Mesquita et al. [19] and by Adolph et al. [20].

Distributed solutions for the 2D static case and their application to the Boundary Element Method have been reported by Crouch and Starfield [21]. A series of static solutions for distributed loads in transversely isotropic half-spaces have been listed by Wang and Liao [22]. Static closed-form solutions for transversely isotropic 3D half-spaces for the case of buried rectangular and triangular loads, with nonconstant spatial load variation, were synthesized by Wang and his coworkers [23, 24].

Distributed load solutions for full spaces or half-spaces can be used directly as elements of analysis in geomechanical engineering problems [25] or as auxiliary states, auxiliary solutions, to model more general problems together with a numerical scheme such as the Boundary Element Method [14, 21, 26].

In the present paper a dynamic stationary semianalytical solution for a spatially constant load applied over a rectangular surface within a viscoelastic isotropic full space is presented. The solution is obtained within the frame of a double Fourier integral transform. Close form solutions for general loadings within the full space are furnished in the transformed wave number domain. Expressions for the solution in the original physical domain for the case of a rectangular load of constant amplitude are furnished as inverse double Fourier integrals. In an accompanying paper the strategy used to perform the numerical integration of these expressions is addressed. The validation and a series of dynamic results are presented.

2. Problem Statement

The boundary value problem to be solved consists of a 3D isotropic full space subjected to a harmonic load. The presented formulation allows any spatially constant load distribution to be considered. Solution expressions for concentrated loads (Green’s functions) are indicated but not evaluated. The displacement solutions for a spatially uniformly distributed normal and tangential load over a rectangular surface with dimensions 2𝐴×2𝐵 on the plane 𝑥-𝑦 are furnished and numerically evaluated. Figures 1(a) to 1(c) show the cases of normal and tangential loads applied to a rectangular surface.

(a) first BVP load 𝑝𝑍

(b) second BVP load 𝑝𝑌

(c) third BVP load 𝑝𝑋

Figure 1 

Traction loadings for the three boundary value problems.

The problem is governed by the Navier-Cauchy differential equations which, in terms of the displacement components 𝑈_𝑖=𝑈_𝑖(𝑥,𝑦,𝑧,𝜔)(𝑖=𝑥,𝑦,𝑧) and in the absence of body forces, may be expressed as [27]𝜇^∗𝑈_{𝑖,𝑗𝑗}+𝜆^∗+𝜇^∗𝑈_{𝑘,𝑘𝑖}=−𝜔²𝜌𝑈_𝑖. (2.1)

In (2.1), 𝜌 is the density of the medium and 𝜔 is the circular frequency and 𝜇^∗ and 𝜆^∗ are complex Lamé parameters, containing the elastic constitutive parameters of the medium 𝜆, 𝜇, and a viscoelastic model represented by frequency-dependent internal damping coefficients [28]:𝜆^∗=𝜆𝜂_𝑅𝜆+𝑖𝜂_𝐼𝜆(𝜔,𝜇_∗)=𝜇𝜂_𝑅𝜇+𝑖𝜂_𝐼𝜇(𝜔). (2.2)

In the present work the internal damping coefficients are assumed to be equal, 𝜂_𝐼𝜆=𝜂_𝐼𝜇=𝜂. The real part of the complex coefficient is taken to be constant and of unit value, 𝜂_𝑅𝜆=𝜂_𝑅𝜇=1. Considering the Kroenecker Delta 𝛿_𝑖𝑗, for a linear isotropic continua the components of the stress tensor 𝜎_𝑖𝑗 can be expressed in terms of the displacement components 𝑈_𝑖:𝜎_𝑖𝑗=𝜇^∗𝑈_𝑖,𝑗+𝑈_𝑗,𝑖+𝛿_𝑖𝑗𝜆^∗𝑈_𝑘,𝑘. (2.3)

3. Solution Strategy

The solution strategy adopted in this paper is based on the classical Helmhotz decomposition of a vector field into a scalar dilatation and a vector rotational field [29].

With the aid of the vector identity∇²𝐔=∇(∇•𝐔)−∇×(∇×𝐔) (3.1) which may also be written in indicial notation𝑈_{𝑖,𝑗𝑗}=𝑈_{𝑗,𝑗𝑖}−𝑒_𝑖𝑗𝑘𝑒_𝑘𝑙𝑚𝑈_{𝑚,𝑙𝑗}. (3.2) Equation (2.1) may be written as𝜇^∗𝑈_{𝑗,𝑗𝑖}−𝑒_𝑖𝑗𝑘𝑒_𝑘𝑙𝑚𝑈_{𝑚,𝑙𝑗}+𝜆^∗+𝜇^∗𝑈_{𝑘,𝑘𝑖}=−𝜔²𝜌𝑈_𝑖. (3.3)

Now using the following definition for the rotation components Ω_𝑛:2Ω_𝑛𝑒_𝑛̂=𝑒_𝑛𝑘𝑙𝑈_𝑙,𝑘𝑒_𝑛̂=∇×𝐔, (3.4) and for the dilatation partΔ=𝑈_𝑗,𝑗=∇•𝐔. (3.5) Equation (3.3) can be recast into𝜆^∗+2𝜇^∗𝜔²𝜌Δ_,𝑖−𝜇^∗𝜔²𝜌2𝑒_𝑖𝑚𝑛Ω_𝑛,𝑚+𝑈_𝑖=0. (3.6)

Defining the complex velocities of dilatational 𝑐^∗_𝐿 and distortional 𝑐^∗_𝑆 waves and complex wave numbers 𝑘^∗_𝐿 and 𝑘^∗_𝑆 such that𝑘^∗_𝐿²=𝜔²𝑐^∗_𝐿²,𝑐^∗_𝐿²=𝜆^∗+2𝜇^∗𝜌,𝑘^∗_𝑆²=𝜔²𝑐^∗_𝑆²,𝑐^∗_𝑆²=𝜇^∗𝜌. (3.7) Equation (3.6) may be written as𝑈_𝑖1=−𝑘^∗_𝐿²Δ_,𝑖+2𝑘^∗_𝑆²𝑒_𝑖𝑚𝑛Ω_𝑛,𝑚=0. (3.8) Equation (3.6) illustrates well the very complex nature of the displacement equations of motion. Considering the vector identities𝑒_𝑖𝑗𝑘Δ_𝑘,𝑗𝜕=𝟎,𝜕𝑥_𝑖𝑒_𝑖𝑚𝑛Ω_𝑛,𝑚=0. (3.9) It is possible to transform (3.6) into two decoupled wave equations by taking the divergent and the curl of that expression, leading to1Δ+𝑘^∗_𝐿²Δ_,𝑖𝑖Ω=0,_𝑖+1𝑘^∗_𝑆²Ω_{𝑖,𝑗𝑗}=0. (3.10)(3.11) Equations (3.10) and (3.11) shows that the displacement field 𝑈_𝑖 is composed of a rotation free scalar dilatation field Δ and a dilatation-free rotational field Ω_𝑖. This decomposition of the equations of motion (2.1) into two independent wave equations is known as the Helmholtz decomposition of a vector field. The decomposition is complete under the restriction that the rotation field must be divergent-free [29]:Ω_𝑖,𝑖=0. (3.12)

The stress field can also be expressed in terms of dilatation and rotation parts by substituting (3.8) in (2.3):𝜎_𝑖𝑗=𝜇^∗𝛿_𝑖𝑗1−2𝑛²𝑛²2Δ−𝑘^∗_𝐿²Δ_,𝑖𝑗+2𝑘^∗_𝑆²𝑒_𝑖𝑘𝑙Ω_{𝑙,𝑘𝑗}+𝑒_𝑗𝑙𝑘Ω_{𝑘,𝑙𝑖}. (3.13)

In (3.13), the following parameter is used:𝑛²=𝑘^∗_𝐿²𝑘^∗_𝑆². (3.14)

Trial Solutions
To proceed with the solution, the unbounded full space is divided in two domains, −∞<𝑧≤0 and 0<𝑧<+∞, respectively. For each domain a trial solution (Ansatz function) is given.
For domain 1,Δ⁽¹⁾=𝐴⁽¹⁾𝑘^∗_𝐿²exp𝛼_𝐿,Ω𝑧+𝑖𝛽𝑥+𝛾𝑦_𝑗(1)=𝐵_𝑗⁽¹⁾𝑘^∗_𝑆²exp𝛼_𝑆𝑧+𝑖𝛽𝑥+𝛾𝑦. (3.15)(3.16)
And for domain 2,Δ⁽²⁾=𝐴⁽²⁾𝑘^∗_𝐿²exp−𝛼_𝐿,Ω𝑧+𝑖𝛽𝑥+𝛾𝑦_𝑗(2)=𝐵_𝑗⁽²⁾𝑘^∗_𝑆²exp−𝛼_𝑆𝑧+𝑖𝛽𝑥+𝛾𝑦. (3.17)(3.18)
The superscripts (1) and (2) in (3.15) to (3.18) are used to denote quantities associated with the domains bounded by −∞<𝑧≤0 and 0<𝑧<+∞, respectively. The 8 constants present at the trial solutions (3.15) to (3.18), namely, 𝐴^(𝑚) and 𝐵_𝑗^(𝑚) with (𝑗=1,2,3) and (𝑚=1,2) are to be determined by the boundary conditions prescribed at the interface of the two domains.
If the trial functions (3.15) to (3.18) are substituted into (3.10) and (3.11) then, for a nontrivial solution of the problem, 𝐴^(𝑚)≠0 and 𝐵_𝑗^(𝑚)≠0 with (𝑗=1,2,3) and (𝑚=1,2), it is necessary that the following conditions are satisfied:𝛼²_𝐿=𝛽²+𝛾²−𝑘^∗_𝐿²,𝛼²_𝑆=𝛽²+𝛾²−𝑘^∗_𝑆². (3.19)
The trial solutions (3.15) to (3.18) must also satisfy the Sommerfeld radiation condition, which states that the waves at points far from the energy source must be outgoing and not incoming waves [1]. In the context of the trial functions (3.15) to (3.18), the radiation condition implies that the real part of the coefficients 𝛼_𝐿 and 𝛼_𝑆 must be a positive value. The variables 𝛽 and 𝛾 in the trial solutions are wave numbers corresponding, respectively, to the 𝑥 and 𝑦 Cartesian directions.
The additional criterion given in (3.12), which establishes that the vector field Ω has zero divergence (Ω_𝑖,𝑖=0), allows the three components of the displacement field 𝑈_𝑖 to be uniquely determined from the four components of Δ and Ω_𝑛(𝑛=1,2,3) [29]. Applying (3.12) in (3.16) and (3.18) yields𝐵₃⁽¹⁾=−𝑖𝛼_𝑆𝛽𝐵₁+𝛾𝐵₂,𝐵₃⁽²⁾=𝑖𝛼_𝑆𝛽𝐵₁+𝛾𝐵₂. (3.20)

Displacement Solutions
The expressions for the displacement field of the domain 1 are obtained when (3.15) and (3.16) are substituted into (3.8) 𝑈_𝑋⁽¹⁾=−𝐴⁽¹⁾𝑖𝛽exp𝛼_𝐿𝑧+2𝛼_𝑆𝐵₁⁽¹⁾𝛽𝛾+𝐵₂⁽¹⁾𝛾²−𝛼²_𝑆exp𝛼_𝑆𝑧𝑖,𝑈exp𝛽𝑥+𝛾𝑦_𝑌(1)=−𝐴⁽¹⁾𝑖𝛾exp𝛼_𝐿𝑧+2𝛼_𝑆𝐵₁⁽¹⁾−𝛽²+𝛼²_𝑆−𝐵₂⁽¹⁾𝛽𝛾exp𝛼_𝑆𝑧𝑖𝑈exp𝛽𝑥+𝛾𝑦,_𝑍(1)=−𝐴⁽¹⁾𝛼_𝐿exp𝛼_𝐿𝑧−2𝑖𝐵₁⁽¹⁾𝛾−𝐵₂⁽¹⁾𝛽exp𝛼_𝑆𝑧exp𝑖𝛽𝑥+𝛾𝑦. (3.21)(3.22)(3.23)
Analogously, (3.17) and (3.18) are substituted into (3.8) to obtain the displacement field of the domain 2:𝑈_𝑋⁽²⁾=−𝐴⁽²⁾𝑖𝛽exp−𝛼_𝐿𝑧−2𝛼_𝑆𝐵₁⁽²⁾𝛽𝛾−𝐵₂⁽²⁾−𝛾²+𝛼²_𝑆exp−𝛼_𝑆𝑧𝑖,𝑈exp𝛽𝑥+𝛾𝑦_𝑌(2)=−𝐴⁽²⁾𝑖𝛾exp−𝛼_𝐿𝑧+2𝛼_𝑆𝐵₁⁽²⁾𝛽²−𝛼²_𝑆+𝐵₂⁽²⁾𝛽𝛾exp−𝛼_𝑆𝑧𝑖𝑈exp𝛽𝑥+𝛾𝑦,_𝑍(2)=𝐴⁽²⁾𝛼_𝐿exp−𝛼_𝐿𝑧−2𝑖𝐵₁⁽²⁾𝛾−𝐵₂⁽²⁾𝛽exp−𝛼_𝑆𝑧exp𝑖𝛽𝑥+𝛾𝑦. (3.24)(3.25)(3.26)

Stress Solutions
The expressions for the stress field of the domain 1 are obtained when (3.15) and (3.16) are substituted into (3.13): 𝜎⁽¹⁾_𝑋𝑋=𝜇^∗𝐴⁽¹⁾𝛽²−𝛾²−𝛼²_𝑆+2𝛼²_𝐿𝑒^𝛼𝐿𝑧+4𝑖𝛽𝛼_𝑆𝐵₁⁽¹⁾𝛽𝛾+𝐵₂⁽¹⁾𝛾²−𝛼²_𝑆𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝜎⁽¹⁾_𝑋𝑌=2𝜇^∗𝐴⁽¹⁾𝛽𝛾𝑒^𝛼𝐿𝑧+𝑖𝛼_𝑆𝐵₁⁽¹⁾𝛾²𝛽+𝛽𝛼²_𝑆−𝛽³+𝐵₂⁽¹⁾𝛾³−𝛾𝛼²_𝑆−𝛽²𝛾𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}𝜎_(1)𝑋𝑍,=2𝜇^∗−𝑖𝐴⁽¹⁾𝛽𝛼_𝐿𝑒^𝛼𝐿𝑧+2𝐵₁⁽¹⁾𝛽𝛾+𝐵₂⁽¹⁾𝛾²−𝛼²_𝑆−𝛽²𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}𝜎_(1)𝑌𝑌,=𝜇^∗𝐴⁽¹⁾𝛾²−𝛽²−𝛼²_𝑆+2𝛼²_𝐿𝑒^𝛼𝐿𝑧+4𝑖𝛾𝛼_𝑆𝐵₁⁽¹⁾𝛼²_𝑆−𝛽²−𝐵₂⁽¹⁾𝛽𝛾𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}𝜎_(1)𝑌𝑍,=2𝜇^∗−𝑖𝐴⁽¹⁾𝛾𝛼_𝐿𝑒^𝛼𝐿𝑧+𝐵₁⁽¹⁾𝛼²_𝑆+𝛾²−𝛽²−2𝐵₂⁽¹⁾𝛽𝛾𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}𝜎_(1)𝑍𝑍,=𝜇^∗−𝐴⁽¹⁾𝛾²+𝛽²+𝛼²_𝑆𝑒^𝛼𝐿𝑧−4𝑖𝛼_𝑆𝐵₁⁽¹⁾𝛾−𝐵₂⁽¹⁾𝛽𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}. (3.27)(3.28)(3.29)(3.30)(3.31)(3.32)
Analogously, (3.17) and (3.21) are substituted into (3.13) to obtain the stress field of the domain 2:𝜎⁽²⁾_𝑋𝑋=𝜇^∗𝐴⁽²⁾𝛽²−𝛾²−𝛼²_𝑆+2𝛼²_𝐿𝑒^{−𝛼𝐿𝑧}−4𝑖𝛽𝛼_𝑆𝐵₁⁽²⁾𝛽𝛾+𝐵₂⁽²⁾𝛾²−𝛼²_𝑆𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝜎⁽²⁾_𝑋𝑌=2𝜇^∗𝐴⁽²⁾𝛽𝛾𝑒^{−𝛼𝐿𝑧}−𝑖𝛼_𝑆𝐵₁⁽²⁾𝛾²𝛽+𝛽𝛼²_𝑆−𝛽³+𝐵₂⁽²⁾𝛾³−𝛾𝛼²_𝑆−𝛽²𝛾𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}𝜎_(2)𝑋𝑍,=2𝜇^∗𝑖𝐴⁽¹⁾𝛽𝛼_𝐿𝑒^{−𝛼𝐿𝑧}+2𝐵₁⁽²⁾𝛽𝛾+𝐵₂⁽²⁾𝛾²−𝛼²_𝑆−𝛽²𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}𝜎_(2)𝑌𝑌,=𝜇^∗𝐴⁽²⁾𝛾²−𝛽²−𝛼²_𝑆+2𝛼²_𝐿𝑒^{−𝛼𝐿𝑧}−4𝑖𝛾𝛼_𝑆𝐵₁⁽²⁾𝛼²_𝑆−𝛽²−𝐵₂⁽²⁾𝛽𝛾𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}𝜎_(2)𝑌𝑍,=2𝜇^∗𝑖𝐴⁽²⁾𝛾𝛼_𝐿𝑒^{−𝛼𝐿𝑧}+𝐵₁⁽²⁾𝛼²_𝑆+𝛾²−𝛽²−2𝐵₂⁽²⁾𝛽𝛾𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}𝜎_(2)𝑍𝑍,=𝜇^∗−𝐴⁽²⁾𝛾²+𝛽²+𝛼²_𝑆𝑒^{−𝛼𝐿𝑧}+4𝑖𝛼_𝑆𝐵₁⁽²⁾𝛾−𝐵₂⁽²⁾𝛽𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}. (3.33)(3.34)(3.35)(3.36)(3.37)(3.38)

4. Boundary-Value Problem

The displacement and stress solutions given in (3.21) to (3.38) can be regarded as solutions in the wave number domain (𝛽,𝛾). It can be shown that these solutions written in terms of the trial functions (3.15) to (3.18) are tantamount to solutions in the transformed Fourier wave number domain (𝛽,𝛾) [12, 13, 26]. There are 6 independent constants to be determined from the 6 boundary conditions at the domains’ interface. The 8 constants 𝐴^(𝑚) and 𝐵_𝑗^(𝑚) with (𝑗=1,2,3) and (𝑚=1,2) that are subjected to the 2 restrictions given in (3.20).

Since (3.21) to (3.38) are in the Fourier wave number domain (𝛽,𝛾), the boundary conditions at the interface at the two medium must also be given in this wave number domain. There are two kinds of boundary conditions to be prescribed at the interface, namely, displacement and tractions conditions.

Displacement Boundary Conditions
Regardless of the traction boundary conditions, displacement continuity at the domain interfaces (𝑧=0) is prescribed. Using (3.21) to (3.26) continuity conditions gives rise to the following set of equations: 𝑈_𝑋⁽¹⁾𝑥,𝑦,𝑧=0−𝑈_𝑋⁽²⁾=𝑥,𝑦,𝑧=0−𝑖𝛽𝐴⁽¹⁾+𝑖𝛽𝐴⁽²⁾+2𝛽𝛾𝛼_𝑆𝐵₁⁽¹⁾+2𝛾²−𝛼²_𝑆𝛼_𝑆𝐵₂⁽¹⁾+2𝛽𝛾𝛼_𝑆𝐵₁⁽²⁾+2𝛾²−𝛼²_𝑆𝛼_𝑆𝐵₂⁽²⁾𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}𝑈=0,_𝑌(1)𝑥,𝑦,𝑧=0−𝑈_𝑌⁽²⁾=𝑥,𝑦,𝑧=0−𝑖𝛾𝐴⁽¹⁾+𝑖𝛾𝐴⁽²⁾−2𝛽²−𝛼²_𝑆𝛼_𝑆𝐵₁⁽¹⁾−2𝛽𝛾𝛼_𝑆𝐵₂⁽¹⁾−2𝛽²−𝛼²_𝑆𝛼_𝑆𝐵₁⁽²⁾−2𝛽𝛾𝛼_𝑆𝐵₂⁽²⁾𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}𝑈=0,_𝑍(1)𝑥,𝑦,𝑧=0−𝑈_𝑍⁽²⁾=𝑥,𝑦,𝑧=0−𝛼_𝐿𝐴⁽¹⁾−𝛼_𝐿𝐴⁽²⁾−2𝑖𝛾𝐵₁⁽¹⁾+2𝑖𝛽𝐵₂⁽¹⁾+2𝑖𝛾𝐵₁⁽²⁾−2𝑖𝛽𝐵₂⁽²⁾𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}=0. (4.1)(4.2)(4.3)

Traction Boundary Conditions
The stress boundary conditions vary according to whether the loading is vertical (applied in the direction of 𝑧) or transversal (applied in the directions of 𝑥 or 𝑦), see Figure 1. These distinct boundary conditions will give rise to three boundary value problems. The following subsections investigate each of these cases.

4.1. The First BVP-Displacements Field due to Vertical Load in the Wave Number Domain

Consider the case in which the infinite medium is subjected to a harmonic vertical load, applied in its plane 𝑥-𝑦(𝑧=0) (see Figure 1(a)), given by the following expression in the Fourier wave number domain:𝑝_𝑍𝛽,𝛾=𝑝_𝑍𝛽,𝛾exp𝑖𝛽𝑥+𝛾𝑦. (4.4)

The boundary conditions prescribe stress continuity at the interface in the tangential directions 𝑥 and 𝑦, and a jump or discontinuity in the normal direction 𝑧 due to the applied normal load 𝑝_𝑍(𝛽,𝛾), leading to the expressions: 𝜎⁽¹⁾_𝑍𝑍𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑍𝑥,𝑦,𝑧=0=𝑝_𝑍,𝜎𝛽,𝛾_(1)𝑍𝑋𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑋𝜎𝑥,𝑦,𝑧=0=0,_(1)𝑍𝑌𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑌𝑥,𝑦,𝑧=0=0. (4.5)

These boundary conditions are applied to (3.27) to (3.38) to obtain:𝜎⁽¹⁾_𝑍𝑍𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑍𝑥,𝑦,𝑧=0=𝜇^∗−𝛾²+𝛽²+𝛼²_𝑆𝐴⁽¹⁾−4𝑖𝛼_𝑆𝛾𝐵₁⁽¹⁾+4𝑖𝛼_𝑆𝛽𝐵₂⁽¹⁾+𝛾²+𝛽²+𝛼²_𝑆𝐴⁽²⁾−4𝑖𝛼_𝑆𝛾𝐵₁⁽²⁾+4𝑖𝛼_𝑆𝛽𝐵₂⁽²⁾𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}=𝑝_𝑍𝛽,𝛾𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝜎⁽¹⁾_𝑍𝑋𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑋𝑥,𝑦,𝑧=0=2𝜇^∗−𝑖𝛼_𝐿𝛽𝐴⁽¹⁾+2𝛽𝛾𝐵₁⁽¹⁾+𝛾²−𝛽²−𝛼²_𝑆𝐵₂⁽¹⁾−𝑖𝛼_𝐿𝛽𝐴⁽²⁾−2𝛽𝛾𝐵₁⁽²⁾−𝛾²−𝛽²−𝛼²_𝑆𝐵₂⁽²⁾𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}𝜎=0,_(1)𝑍𝑌𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑌𝑥,𝑦,𝑧=0=2𝜇^∗−𝑖𝛼_𝐿𝛾𝐴⁽¹⁾+𝛾²−𝛽²+𝛼²_𝑆𝐵₁⁽¹⁾−2𝛽𝛾𝐵₂⁽¹⁾−𝑖𝛼_𝐿𝛾𝐴⁽²⁾−𝛾²−𝛽²+𝛼²_𝑆𝐵₁⁽²⁾+2𝛽𝛾𝐵₂⁽²⁾𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}=0. (4.6)(4.7)(4.8)

Equations (4.1) to (4.3) and (4.6) to (4.8) form an algebraic system of six equations, from which the six unknowns of the problem, 𝐴₁^(𝑚), 𝐵₁^(𝑚), and 𝐵₂^(𝑚)(𝑚=1,2), can be obtained: −𝑖𝛽𝐴⁽¹⁾+𝑖𝛽𝐴⁽²⁾+2𝛽𝛾𝛼_𝑆𝐵₁⁽¹⁾+2𝛽𝛾𝛼_𝑆𝐵₁⁽²⁾+2𝛾²−𝛼²_𝑆𝛼_𝑆𝐵₂⁽¹⁾+2𝛾²−𝛼²_𝑆𝛼_𝑆𝐵₂⁽²⁾=0,−𝑖𝛾𝐴⁽¹⁾+𝑖𝛾𝐴⁽²⁾−2𝛽²−𝛼²_𝑆𝛼_𝑆𝐵₁⁽¹⁾−2𝛽²−𝛼²_𝑆𝛼_𝑆𝐵₁⁽²⁾−2𝛽𝛾𝛼_𝑆𝐵₂⁽¹⁾−2𝛽𝛾𝛼_𝑆𝐵₂⁽²⁾=0,−𝛼_𝐿𝐴⁽¹⁾−𝛼_𝐿𝐴⁽²⁾−2𝑖𝛾𝐵₁⁽¹⁾+2𝑖𝛾𝐵₁⁽²⁾+2𝑖𝛽𝐵₂⁽¹⁾−2𝑖𝛽𝐵₂⁽²⁾−𝛾=0,₂+𝛽²+𝛼²_𝑆𝐴⁽¹⁾+𝛾²+𝛽²+𝛼²_𝑆𝐴⁽²⁾−4𝑖𝛼_𝑆𝛾𝐵₁⁽¹⁾−4𝑖𝛼_𝑆𝛾𝐵₁⁽²⁾+4𝑖𝛼_𝑆𝛽𝐵₂⁽¹⁾+4𝑖𝛼_𝑆𝛽𝐵₂⁽²⁾=𝑝_𝑍𝜇^∗,−𝑖𝛼_𝐿𝛽𝐴⁽¹⁾−𝑖𝛼_𝐿𝛽𝐴⁽²⁾+2𝛽𝛾𝐵₁⁽¹⁾−2𝛽𝛾𝐵₁⁽²⁾+𝛾²−𝛽²−𝛼²_𝑆𝐵₂⁽¹⁾−𝛾²−𝛽²−𝛼²_𝑆𝐵₂⁽²⁾=0,−𝑖𝛼_𝐿𝛾𝐴⁽¹⁾−𝑖𝛼_𝐿𝛾𝐴⁽²⁾+𝛾²−𝛽²+𝛼²_𝑆𝐵₁⁽¹⁾−𝛾²−𝛽²+𝛼²_𝑆𝐵₁⁽²⁾−2𝛽𝛾𝐵₂⁽¹⁾+2𝛽𝛾𝐵₂⁽²⁾=0. (4.9)

The solution of (4.9) leads to𝐴⁽¹⁾=−𝐴⁽²⁾=12𝑝_𝑍𝛽,𝛾𝜇^∗𝑘^∗_𝑆²,𝐵₁⁽¹⁾=𝐵₁⁽²⁾=14𝑖𝛾𝑝_𝑍𝛽,𝛾𝜇^∗𝛼_𝑆𝑘^∗_𝑆²𝐵₂₍₁₎,=𝐵₂⁽²⁾1=−4𝑖𝛽𝑝_𝑍𝛽,𝛾𝜇^∗𝛼_𝑆𝑘^∗_𝑆². (4.10)(4.11)(4.12)

The displacement field in the Fourier transformed wave number domain due to the vertical loading 𝑝_𝑍(𝛽,𝛾) is obtained by substituting (4.10) to (4.12) into (3.21) to (3.26):𝑈⁽¹⁾_𝑋𝑍1=−2𝑖𝛽𝑝_𝑍𝛽,𝛾𝜇^∗𝑘^∗_𝑆²𝑒^𝛼𝐿𝑧−𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽¹⁾_𝑌𝑍1=−2𝑖𝛾𝑝_𝑍𝛽,𝛾𝜇^∗𝑘^∗_𝑆²𝑒^𝛼𝐿𝑧−𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽¹⁾_𝑍𝑍1=−2𝑝_𝑍𝛽,𝛾𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝛼_𝐿𝑒^𝛼𝐿𝑧−𝛾²+𝛽²𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽²⁾_𝑋𝑍=12𝑖𝛽𝑝_𝑍𝛽,𝛾𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿𝑧}−𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽²⁾_𝑌𝑍=12𝑖𝛾𝑝_𝑍𝛽,𝛾𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿𝑧}−𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽²⁾_𝑍𝑍1=−2𝑝_𝑍𝛽,𝛾𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝛼_𝐿𝑒^{−𝛼𝐿𝑧}−𝛾²+𝛽²𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}. (4.13)

In this paper, 𝑈^(𝑚)_𝑖𝑘 indicates the displacement of a point at medium 𝑚 in the direction 𝑖 due to a loading in the direction 𝑘.

4.2. The Second BVP-Displacements Field due to the Tangential Load in 𝑌 Direction in the Wave Number Domain

Consider the case in which the infinite medium is subjected to a harmonic vertical load, applied in its plane 𝑥-𝑦(𝑧=0) (see Figure 1(b)), given by the following expression in the Fourier wave number domain:𝑝_𝑌𝛽,𝛾=𝑝_𝑌𝛽,𝛾exp𝑖𝛽𝑥+𝛾𝑦. (4.14)

The boundary conditions prescribe stress continuity at the interface in the normal direction 𝑧 and in the tangential directions 𝑥 as well as a jump or discontinuity in the 𝑦 direction due to the applied normal load 𝑝_𝑌(𝛽,𝛾), leading to the expressions:𝜎⁽¹⁾_𝑍𝑍𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑍𝜎𝑥,𝑦,𝑧=0=0,_(1)𝑍𝑋𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑋𝜎𝑥,𝑦,𝑧=0=0,_(1)𝑍𝑌𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑌𝑥,𝑦,𝑧=0=𝑝_𝑌𝛽,𝛾. (4.15)

These traction boundary conditions (4.15) substituted in (3.27) to (3.38) together with continuity equations (4.1) to (4.3) form an algebraic system of six equations, from which the six unknowns of the problem, 𝐴₁^(𝑚), 𝐵₁^(𝑚), and 𝐵₂^(𝑚)(𝑚=1,2), can be obtained: −𝑖𝛽𝐴⁽¹⁾+𝑖𝛽𝐴⁽²⁾+2𝛽𝛾𝛼_𝑆𝐵₁⁽¹⁾+2𝛽𝛾𝛼_𝑆𝐵₁⁽²⁾+2𝛾²−𝛼²_𝑆𝛼_𝑆𝐵₂⁽¹⁾+2𝛾²−𝛼²_𝑆𝛼_𝑆𝐵₂⁽²⁾=0,−𝑖𝛾𝐴⁽¹⁾+𝑖𝛾𝐴⁽²⁾−2𝛽²−𝛼²_𝑆𝛼_𝑆𝐵₁⁽¹⁾−2𝛽²−𝛼²_𝑆𝛼_𝑆𝐵₁⁽²⁾−2𝛽𝛾𝛼_𝑆𝐵₂⁽¹⁾−2𝛽𝛾𝛼_𝑆𝐵₂⁽²⁾=0,−𝛼_𝐿𝐴⁽¹⁾−𝛼_𝐿𝐴⁽²⁾−2𝑖𝛾𝐵₁⁽¹⁾+2𝑖𝛾𝐵₁⁽²⁾+2𝑖𝛽𝐵₂⁽¹⁾−2𝑖𝛽𝐵₂⁽²⁾−𝛾=0,₂+𝛽²+𝛼²_𝑆𝐴⁽¹⁾+𝛾²+𝛽²+𝛼²_𝑆𝐴⁽²⁾−4𝑖𝛼_𝑆𝛾𝐵₁⁽¹⁾−4𝑖𝛼_𝑆𝛾𝐵₁⁽²⁾+4𝑖𝛼_𝑆𝛽𝐵₂⁽¹⁾+4𝑖𝛼_𝑆𝛽𝐵₂⁽²⁾=0,−𝑖𝛼_𝐿𝛽𝐴⁽¹⁾−𝑖𝛼_𝐿𝛽𝐴⁽²⁾+2𝛽𝛾𝐵₁⁽¹⁾−2𝛽𝛾𝐵₁⁽²⁾+𝛾²−𝛽²−𝛼²_𝑆𝐵₂⁽¹⁾−𝛾²−𝛽²−𝛼²_𝑆𝐵₂⁽²⁾=0,−𝑖𝛼_𝐿𝛾𝐴⁽¹⁾−𝑖𝛼_𝐿𝛾𝐴⁽²⁾+𝛾²−𝛽²+𝛼²_𝑆𝐵₁⁽¹⁾−𝛾²−𝛽²+𝛼²_𝑆𝐵₁⁽²⁾−2𝛽𝛾𝐵₂⁽¹⁾+2𝛽𝛾𝐵₂⁽²⁾=𝑝_𝑌2𝜇^∗. (4.16) The solution of the system (4.16) leads to𝐴⁽¹⁾=𝐴⁽²⁾=12𝛾𝑝_𝑌𝛽,𝛾𝛼_𝐿𝜇^∗𝑘^∗_𝑆²,𝐵₁⁽¹⁾=−𝐵₁⁽²⁾1=−4𝑝_𝑌𝛽,𝛾𝜇^∗𝑘^∗_𝑆²,𝐵₂⁽¹⁾=𝐵₂⁽²⁾=0. (4.17)

The displacement field in the Fourier transformed wave number domain due to the transversal loading 𝑝_𝑌(𝛽,𝛾) is obtained by substituting (4.17) into (3.21) to (3.26):𝑈⁽¹⁾_𝑋𝑌=12𝛾𝛽𝑝_𝑌𝛽,𝛾𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝑒^𝛼𝐿𝑧−𝛼_𝐿𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽¹⁾_𝑌𝑌=12𝑝_𝑌𝛽,𝛾𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛾²𝛼_𝑆𝑒^𝛼𝐿𝑧+𝛼_𝐿𝛽²−𝛼²_𝑆𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽¹⁾_𝑍𝑌1=−2𝑖𝛾𝑝_𝑌𝛽,𝛾𝜇^∗𝑘^∗_𝑆²𝑒^𝛼𝐿𝑧−𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽²⁾_𝑋𝑌=12𝛾𝛽𝑝_𝑌𝛽,𝛾𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝑒^{−𝛼𝐿𝑧}−𝛼_𝐿𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽²⁾_𝑌𝑌=12𝑝_𝑌𝛽,𝛾𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛾²𝛼_𝑆𝑒^{−𝛼𝐿𝑧}+𝛼_𝐿𝛽²−𝛼²_𝑆𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽²⁾_𝑍𝑌=12𝑖𝛾𝑝_𝑌𝛽,𝛾𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿𝑧}−𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}. (4.18)

4.3. The Third BVP-Displacements Field due to the Tangential Load in 𝑋 Direction in the Wave Number Domain

Consider the case in which the infinite medium is subjected to a harmonic vertical load, applied in its plane 𝑥-𝑦(𝑧=0) (see Figure 1(c)), given by the following expression in the Fourier wave number domain:𝑝_𝑋𝛽,𝛾=𝑝_𝑋𝛽,𝛾exp𝑖𝛽𝑥+𝛾𝑦. (4.19)

The boundary conditions prescribe stress continuity at the interface in the normal direction 𝑧 and in the tangential directions 𝑦 as well as a jump or discontinuity in the 𝑥 direction due to the applied normal load 𝑝_𝑋(𝛽,𝛾), leading to the expressions:𝜎⁽¹⁾_𝑍𝑍𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑍𝜎𝑥,𝑦,𝑧=0=0,_(1)𝑍𝑋𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑋𝑥,𝑦,𝑧=0=𝑝_𝑋,𝜎𝛽,𝛾_(1)𝑍𝑌𝑥,𝑦,𝑧=0−𝜎⁽²⁾_𝑍𝑌𝑥,𝑦,𝑧=0=0. (4.20)

These traction boundary conditions (4.20) substituted in (3.27) to (3.38) together with continuity equations (4.1) to (4.3) form an algebraic system of six equations, from which the six unknowns of the problem, 𝐴₁^(𝑚), 𝐵₁^(𝑚), and 𝐵₂^(𝑚)(𝑚=1,2), can be obtained:−𝑖𝛽𝐴⁽¹⁾+𝑖𝛽𝐴⁽²⁾+2𝛽𝛾𝛼_𝑆𝐵₁⁽¹⁾+2𝛽𝛾𝛼_𝑆𝐵₁⁽²⁾+2𝛾²−𝛼²_𝑆𝛼_𝑆𝐵₂⁽¹⁾+2𝛾²−𝛼²_𝑆𝛼_𝑆𝐵₂⁽²⁾=0,−𝑖𝛾𝐴⁽¹⁾+𝑖𝛾𝐴⁽²⁾−2𝛽²−𝛼²_𝑆𝛼_𝑆𝐵₁⁽¹⁾−2𝛽²−𝛼²_𝑆𝛼_𝑆𝐵₁⁽²⁾−2𝛽𝛾𝛼_𝑆𝐵₂⁽¹⁾−2𝛽𝛾𝛼_𝑆𝐵₂⁽²⁾=0,−𝛼_𝐿𝐴⁽¹⁾−𝛼_𝐿𝐴⁽²⁾−2𝑖𝛾𝐵₁⁽¹⁾+2𝑖𝛾𝐵₁⁽²⁾+2𝑖𝛽𝐵₂⁽¹⁾−2𝑖𝛽𝐵₂⁽²⁾−𝛾=0,₂+𝛽²+𝛼²_𝑆𝐴⁽¹⁾+𝛾²+𝛽²+𝛼²_𝑆𝐴⁽²⁾−4𝑖𝛼_𝑆𝛾𝐵₁⁽¹⁾−4𝑖𝛼_𝑆𝛾𝐵₁⁽²⁾+4𝑖𝛼_𝑆𝛽𝐵₂⁽¹⁾+4𝑖𝛼_𝑆𝛽𝐵₂⁽²⁾=0,−𝑖𝛼_𝐿𝛽𝐴⁽¹⁾−𝑖𝛼_𝐿𝛽𝐴⁽²⁾+2𝛽𝛾𝐵₁⁽¹⁾−2𝛽𝛾𝐵₁⁽²⁾+𝛾²−𝛽²−𝛼²_𝑆𝐵₂⁽¹⁾−𝛾²−𝛽²−𝛼²_𝑆𝐵₂⁽²⁾=𝑝_𝑋2𝜇^∗,−𝑖𝛼_𝐿𝛾𝐴⁽¹⁾−𝑖𝛼_𝐿𝛾𝐴⁽²⁾+𝛾²−𝛽²+𝛼²_𝑆𝐵₁⁽¹⁾−𝛾²−𝛽²+𝛼²_𝑆𝐵₁⁽²⁾−2𝛽𝛾𝐵₂⁽¹⁾+2𝛽𝛾𝐵₂⁽²⁾=0. (4.21) The solution of the system (4.21) leads to𝐴⁽¹⁾=𝐴⁽²⁾=12𝑖𝛽𝑝_𝑋𝛽,𝛾𝛼_𝐿𝜇^∗𝑘^∗_𝑆²,𝐵₁⁽¹⁾=𝐵₁⁽²⁾𝐵=0,₂₍₁₎=−𝐵₂⁽²⁾=14𝑝_𝑋𝛽,𝛾𝜇^∗𝑘^∗_𝑆². (4.22)

The displacement field in the Fourier transformed wave number domain due to the transversal loading 𝑝_𝑋(𝛽,𝛾) is obtained by substituting (4.22) into (3.21) to (3.26):𝑈⁽¹⁾_𝑋𝑋=12𝑝_𝑋𝛽,𝛾𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛽²𝛼_𝑆𝑒^𝛼𝐿𝑧+𝛼_𝐿𝛾²−𝛼²_𝑆𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽¹⁾_𝑌𝑋1=+2𝛽𝛾𝑝_𝑋𝛽,𝛾𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝑒^𝛼𝐿𝑧−𝛼_𝐿𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽¹⁾_𝑍𝑋1=−2𝑖𝛽𝑝_𝑋𝛽,𝛾𝜇^∗𝑘^∗_𝑆²𝑒^𝛼𝐿𝑧−𝑒^𝛼𝑆𝑧𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽²⁾_𝑋𝑋=12𝑝_𝑋𝛽,𝛾𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛽²𝛼_𝑆𝑒^{−𝛼𝐿𝑧}+𝛼_𝐿𝛾²−𝛼²_𝑆𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽²⁾_𝑌𝑋1=+2𝛽𝛾𝑝_𝑋𝛽,𝛾𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝑒^{−𝛼𝐿𝑧}−𝛼_𝐿𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈⁽²⁾_𝑍𝑋=12𝑖𝛽𝑝_𝑋𝛽,𝛾𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿𝑧}−𝑒^{−𝛼𝑆𝑧}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}. (4.23)

4.4. Displacements Field in the Transformed Space

Equations (4.13), (4.18), and (4.23) express the component of displacement 𝑈^(𝑚)_𝑖𝑘 with distinct expressions for each media 𝑚=1 and 𝑚=2. In the following equations, a single expression for the components 𝑈_𝑖𝑘 is given for both media through a small adaptation of the coordinate 𝑧. In these equations, 𝑈_𝑖𝑘 indicates the displacement of a point of coordinates (𝑥,𝑦,𝑧) of the full space in the direction 𝑖(𝑖=𝑥,𝑦,𝑧) due to a loading applied in the direction 𝑘(𝑘=𝑥,𝑦,𝑧).

For the first problem with normal loading 𝑝_𝑍:𝑈_𝑋𝑍=𝑧𝛽,𝛾|𝑧|𝑖𝛽𝑝_𝑍𝛽,𝛾2𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈_𝑌𝑍=𝑧𝛽,𝛾|𝑧|𝑖𝛾𝑝_𝑍𝛽,𝛾2𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈_𝑍𝑍𝛽,𝛾=−𝑝_𝑍𝛽,𝛾2𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝛼_𝐿𝑒^{−𝛼𝐿|𝑧|}−𝛾²+𝛽²𝑒^{−𝛼𝑆|𝑧|}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}. (4.24) For the second problem with tangential loading 𝑝_𝑌:𝑈_𝑋𝑌𝛽,𝛾=𝛽𝛾𝑝_𝑌𝛽,𝛾2𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝑒^{−𝛼𝐿|𝑧|}−𝛼_𝐿𝑒^{−𝛼𝑆|𝑧|}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈_𝑌𝑌𝛽,𝛾=𝑝_𝑌𝛽,𝛾2𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛾²𝛼_𝑆𝑒^{−𝛼𝐿|𝑧|}+𝛼_𝐿𝛽²−𝛼²_𝑆𝑒^{−𝛼𝑆|𝑧|}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈_𝑍𝑌=𝑧𝛽,𝛾|𝑧|𝑖𝛾𝑝_𝑌𝛽,𝛾2𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}. (4.25) And for the third problem with tangential loading 𝑝_𝑋:𝑈_𝑋𝑋𝛽,𝛾=𝑝_𝑋𝛽,𝛾2𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛽²𝛼_𝑆𝑒^{−𝛼𝐿|𝑧|}+𝛼_𝐿𝛾²−𝛼²_𝑆𝑒^{−𝛼𝑆|𝑧|}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈_𝑌𝑋𝛽,𝛾=𝛽𝛾𝑝_𝑋𝛽,𝛾2𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝑒^{−𝛼𝐿|𝑧|}−𝛼_𝐿𝑒^{−𝛼𝑆|𝑧|}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)},𝑈_𝑍𝑋=𝑧𝛽,𝛾|𝑧|𝑖𝛽𝑝_𝑋𝛽,𝛾2𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑒^{𝑖(𝛽𝑥+𝛾𝑦)}. (4.26) Equations (4.24) to (4.26) furnish displacement solutions for the full space for the three boundary value problems in the wave number domain for any loading traction function 𝑝_𝑖(𝛽,𝛾)(𝑖=𝑥,𝑦,𝑧).

4.5. Uniformly Distributed Traction Loadings in the Wave Number Domain

Consider the vertical harmonic load 𝑝_𝑍(𝑥,𝑦,𝑧) depicted in Figure 1(a). Let this load have a uniform intensity 𝑝_𝑍 distributed over a rectangular surface with dimensions 2𝐴×2𝐵, placed at the interface plane 𝑥-𝑦(𝑧=0). In the original physical domains (𝑥,𝑦,𝑧) this load can be written as𝑝_𝑍𝑥,𝑦,𝑧=0=−𝑝_𝑍,∈|𝑥|𝑦∈|𝑥|𝑦𝑥,𝑦,𝑧=0≤𝐴∪≤𝐵,0,𝑥,𝑦,𝑧=0>𝐴∪>𝐵. (4.27)

The loading can be transformed to the wave number domain 𝑝_𝑍(𝛽,𝛾) using the double integral Fourier transform:𝑝_𝑍𝛽,𝛾=^∞_−∞^∞_−∞𝑝_𝑍𝑥,𝑦,𝑧=0𝑒^{−𝑖𝛽𝑥}𝑑𝑥𝑒^{−𝑖𝛾𝑦}𝑑𝑦,𝑝_𝑍𝛽,𝛾=−𝑝_𝑍^𝐵_−𝐵^𝐴_−𝐴𝑒^{−𝑖𝛽𝑥}𝑑𝑥𝑒^{−𝑖𝛾𝑦}𝑑𝑦 (4.28) and finallly𝑝_𝑍𝛽,𝛾=−4𝑃_𝑍𝛽𝛾sin𝛽𝐴sin𝛾𝐵. (4.29)

Analogously, the tangential loads in the wave number domain may be expressed as𝑝_𝑌𝛽,𝛾=−4𝑃_𝑌,𝛽𝛾sin𝛽𝐴sin𝛾𝐵𝑝_𝑋𝛽,𝛾=−4𝑃_𝑋.𝛽𝛾sin𝛽𝐴sin𝛾𝐵 (4.30)

Substituting (4.29), (4.30), and (83), respectively into expressions (4.24), (4.25), and (4.26) results in the displacement solutions for the three boundary value problems in the transformed wave number domain. The next step is to recover the solution in the original physical Cartesian domain (𝑥,𝑦,𝑧).

It should be noted that the solutions presented in (4.24), (4.25), and (4.26) for the displacements of the three boundary value problems in the wave number domain are quite general, in the sense that any loading function 𝑝_𝑖(𝛽,𝛾)(𝑖=𝑥,𝑦,𝑧) may be used as a traction condition at the interface. If the wave number integral transform of Dirac’s Delta distribution is used as the loading function, the dynamic stationary Green’s function, the full space fundamental solution is obtained.

4.6. Displacement Fields in the Original Cartesian Space

Applying a double inverse Fourier integral to (4.24) to (4.26), considering the traction loadings given in (4.29) to (83), leads to the displacement solutions in the original physical space (𝑥,𝑦,𝑧).

For the first BVP with normal excitation 𝑝_𝑍𝑈_𝑋𝑍1𝑥,𝑦,𝑧=−4𝜋^2∞_−∞𝑏𝑧𝛽,𝛾|𝑧|𝑖2𝑃_𝑍𝛾𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑈𝑑𝛽𝑑𝛾,_𝑌𝑍1𝑥,𝑦,𝑧=−4𝜋^2∞_−∞𝑏𝑧𝛽,𝛾|𝑧|𝑖2𝑃_𝑍𝛽𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑈𝑑𝛽𝑑𝛾,_𝑍𝑍1𝑥,𝑦,𝑧=+4𝜋^2∞_−∞𝑏𝛽,𝛾2𝑃_𝑍𝛽𝛾𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝛼_𝐿𝑒^{−𝛼𝐿|𝑧|}−𝛾²+𝛽²𝑒^{−𝛼𝑆|𝑧|}𝑑𝛽𝑑𝛾. (4.31)

For the second BVP with tangential excitation 𝑝_𝑌𝑈_𝑋𝑌1𝑥,𝑦,𝑧=−4𝜋^2∞_−∞𝑏𝛽,𝛾2𝑃_𝑌𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝑒^{−𝛼𝐿|𝑧|}−𝛼_𝐿𝑒^{−𝛼𝑆|𝑧|}𝑈𝑑𝛽𝑑𝛾,_𝑌𝑌1𝑥,𝑦,𝑧=−4𝜋^2∞_−∞𝑏𝛽,𝛾2𝑃_𝑌𝛽𝛾𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛾²𝛼_𝑆𝑒^{−𝛼𝐿|𝑧|}+𝛼_𝐿𝛽²−𝛼²_𝑆𝑒^{−𝛼𝑆|𝑧|}𝑈𝑑𝛽𝑑𝛾,_𝑍𝑌1𝑥,𝑦,𝑧=−4𝜋^2∞_−∞𝑏𝑧𝛽,𝛾|𝑧|𝑖2𝑃_𝑌𝛽𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑑𝛽𝑑𝛾. (4.32)

For the third BVP with tangential excitation 𝑝_𝑋𝑈_𝑋𝑋1𝑥,𝑦,𝑧=−4𝜋^2∞_−∞𝑏𝛽,𝛾2𝑃_𝑋𝛽𝛾𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛽²𝛼_𝑆𝑒^{−𝛼𝐿|𝑧|}+𝛼_𝐿𝛾²−𝛼²_𝑆𝑒^{−𝛼𝑆|𝑧|}𝑈𝑑𝛽𝑑𝛾,_𝑌𝑋1𝑥,𝑦,𝑧=−4𝜋^2∞_−∞𝑏𝛽,𝛾2𝑃_𝑋𝛼_𝐿𝛼_𝑆𝜇^∗𝑘^∗_𝑆²𝛼_𝑆𝑒^{−𝛼𝐿|𝑧|}−𝛼_𝐿𝑒^{−𝛼𝑆|𝑧|}𝑈𝑑𝛽𝑑𝛾,_𝑍𝑋1𝑥,𝑦,𝑧=−4𝜋^2∞_−∞𝑏𝑧𝛽,𝛾|𝑧|𝑖2𝑃_𝑋𝛾𝜇^∗𝑘^∗_𝑆²𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑑𝛽𝑑𝛾. (4.33)

In (4.31) to (4.33), the following notation has been adopted:𝑏𝛽,𝛾=sin𝛽𝐴sin𝛾𝐵exp𝑖𝛽𝑥+𝛾𝑦. (4.34)

The improper integrals from (4.31) to (4.33) can be simplified according to the behavior of their integrands—whether they are odd or even functions. For the first problem this simplification leads to𝑈_𝑋𝑍=𝑧𝑥,𝑦,𝑧|𝑧|𝐷_𝑁𝑃_𝑍^∞₀^∞₀𝛽𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑠_𝛽𝑠_𝛽𝑥𝐴₀𝑘_𝛽𝑑𝑘_𝛽𝑠_𝛾𝑐_𝛾𝑦𝐴₀𝑏₀𝑘_𝛾𝑑𝑘_𝛾,𝑈_𝑌𝑍=𝑧𝑥,𝑦,𝑧|𝑧|𝐷_𝑁𝑃_𝑍^∞₀^∞₀𝛾𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑠_𝛽𝑐_𝛽𝑥𝐴₀𝑘_𝛽𝑑𝑘_𝛽𝑠_𝛾𝑠_𝛾𝑦𝐴₀𝑏₀𝑑𝑘_𝛾,𝑈_𝑍𝑍𝑥,𝑦,𝑧=𝐷_𝑁𝑃_𝑍^∞₀^∞₀𝛼_𝑆𝛼_𝐿𝑒^{−𝛼𝐿|𝑧|}−𝛽²+𝛾²𝑒^{−𝛼𝑆|𝑧|}𝛼_𝑆𝑠_𝛽𝑐_𝛽𝑥𝐴₀𝑘_𝛽𝑑𝑘_𝛽𝑠_𝛾𝑐_𝛾𝑦𝐴₀𝑏₀𝑘_𝛾𝑑𝑘_𝛾. (4.35) For the second BVP with tangential excitation 𝑝_𝑌𝑈_𝑋𝑌𝑥,𝑦,𝑧=−𝐷_𝑁𝑃_𝑌^∞₀^∞₀𝛼𝛽𝛾^𝑆𝑒^{−𝛼𝐿|𝑧|}−𝛼_𝐿𝑒^{−𝛼𝑆|𝑧|}𝛼_𝑆𝛼_𝐿𝑠_𝛽𝑠_𝛽𝑥𝐴₀𝑘_𝛽𝑑𝑘_𝛽𝑠_𝛾𝑠_𝛾𝑦𝐴₀𝑏₀𝑘_𝛾𝑑𝑘_𝛾,𝑈_𝑌𝑌𝑥,𝑦,𝑧=−𝐷_𝑁𝑃_𝑌^∞₀^∞₀𝛾²𝛼_𝑆𝑒^{−𝛼𝐿|𝑧|}+𝛼_𝐿𝛽²−𝛼²_𝑆𝑒^{−𝛼𝑆|𝑧|}𝛼_𝑆𝛼_𝐿𝑠_𝛽𝑐_𝛽𝑥𝐴₀𝑘_𝛽𝑑𝑘_𝛽𝑠_𝛾𝑐_𝛾𝑦𝐴₀𝑏₀𝑘_𝛾𝑑𝑘_𝛾,𝑈_𝑍𝑌𝑥,𝑦,𝑧=𝐷_𝑁𝑃_𝑌𝑧|𝑧|^∞₀^∞₀𝛾𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑠_𝛽𝑐_𝛽𝑥𝐴₀𝑘_𝛽𝑑𝑘_𝛽𝑠_𝛾𝑠_𝛾𝑦𝐴₀𝑏₀𝑘_𝛾𝑑𝑘_𝛾. (4.36) For the third BVP with tangential excitation 𝑝_𝑋𝑈_𝑋𝑋𝑥,𝑦,𝑧=−𝐷_𝑁𝑃_𝑋^∞₀^∞₀𝛽²𝛼_𝑆𝑒^{−𝛼𝐿|𝑧|}+𝛾²−𝛼²_𝑆𝛼_𝐿𝑒^{−𝛼𝑆|𝑧|}𝛼_𝑆𝛼_𝐿𝑠_𝛽𝑐_𝛽𝑥𝐴₀𝑘_𝛽𝑑𝑘_𝛽𝑠_𝛾𝑐_𝛾𝑦𝐴₀𝑏₀𝑘_𝛾𝑑𝑘_𝛾,𝑈_𝑌𝑋𝑥,𝑦,𝑧=−𝐷_𝑁𝑃_𝑋^∞₀^∞₀𝛼𝛽𝛾^𝑆𝑒^{−𝛼𝐿|𝑧|}−𝛼_𝐿𝑒^{−𝛼𝑆|𝑧|}𝛼_𝑆𝛼_𝐿𝑠_𝛽𝑠_𝛽𝑥𝐴₀𝑘_𝛽𝑑𝑘_𝛽𝑠_𝛾𝑠_𝛾𝑦𝐴₀𝑏₀𝑘_𝛾𝑑𝑘_𝛾,𝑈_𝑍𝑋𝑥,𝑦,𝑧=𝐷_𝑁𝑃_𝑋𝑧|𝑧|^∞₀^∞₀𝛽𝑒^{−𝛼𝐿|𝑧|}−𝑒^{−𝛼𝑆|𝑧|}𝑠_𝛽𝑠_𝛽𝑥𝐴₀𝑘_𝛽𝑑𝑘_𝛽𝑠_𝛾𝑐_𝛾𝑦𝐴₀𝑏₀𝑘_𝛾𝑑𝑘_𝛾. (4.37)

In (4.35) to (4.37), the following notation is adopted:𝐷_𝑁=2𝐴𝐵𝜋²𝜇^∗𝜂_𝑅𝑆+𝑖𝜂_𝐼𝑆,𝑏₀=𝐵𝐴,𝐴𝛽=⁰𝐴𝑘_𝛽,𝐴𝛾=⁰𝐴𝑘_𝛾,𝑘^∗_𝑆²=𝐴₀𝐴²1𝜂_𝑅𝑆+𝑖𝜂_𝐼𝑆,𝑘^∗_𝐿²=𝐴₀𝐴²𝑛²𝜂_𝑅𝐿+𝑖𝜂_𝐼𝐿,𝑠_𝛽=sin𝐴₀𝑘_𝛽,𝑠_𝛽𝑥=sin𝐴₀𝑘_𝛽𝑥𝐴,𝑠_𝛾=sin𝐴₀𝑏₀𝑘_𝛾,𝑠_𝛾𝑦=sin𝐴₀𝑏₀𝑘_𝛾𝑦𝐵,𝑐_𝛽=cos𝐴₀𝑘_𝛽,𝑐_𝛽𝑥=cos𝐴₀𝑘_𝛽𝑥𝐴,𝑐_𝛾=cos𝐴₀𝑏₀𝑘_𝛾,𝑐_𝛾𝑦=cos𝐴₀𝑏₀𝑘_𝛾𝑦𝐵. (4.38)

Equations (4.35) to (4.37) represent the final set of expressions of the displacement field of a three-dimensional isotropic, viscoelastic full space subjected to transversal and vertical harmonic loadings uniformly distributed over a rectangular surface 2𝐴×2𝐵 in the plane 𝑥-𝑦(𝑧=0) of the full space. These expressions must be evaluated numerically.

5. Concluding Remarks

In the present paper, general expressions for displacement response of a three-dimensional isotropic viscoelastic full space under stationary dynamic loading were obtained. A particular solution for the case of a rectangular load of constant spatial amplitude applied at the interior of the full space is given in the transformed wave number domains. Expressions for three boundary value problems, associated to a normal and two tangential loadings in the original physical space, are given in terms of a double inverse Fourier integral. These inverse integral transforms must be evaluated numerically. In the second part of the present paper a strategy to evaluate these integrals is described, the procedure validated and a number of original results are reported.

Acknowledgments

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