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

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.


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

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 2A × 2B on the plane x-y are furnished and numerically evaluated.Figures 1 a to 1 c show the cases of normal and tangential loads applied to a rectangular surface.The problem is governed by the Navier-Cauchy differential equations which, in terms of the displacement components U i U i x, y, z, ω i x, y, z and in the absence of body forces, may be expressed as 27 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 : In the present work the internal damping coefficients are assumed to be equal, η Iλ η Iμ η.The real part of the complex coefficient is taken to be constant and of unit value, η Rλ η Rμ 1. Considering the Kroenecker Delta δ ij , for a linear isotropic continua the components of the stress tensor σ ij can be expressed in terms of the displacement components U i :

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 Defining the complex velocities of dilatational c * L and distortional c * S waves and complex wave numbers k * L and k * S such that

3.7
Equation 3.6 may be written as Equation 3.6 illustrates well the very complex nature of the displacement equations of motion.Considering the vector identities It is possible to transform 3.6 into two decoupled wave equations by taking the divergent and the curl of that expression, leading to Equations 3.10 and 3.11 shows that the displacement field U i is composed of a rotation free scalar dilatation field Δ and a dilatation-free rotational field Ω i .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 : The stress field can also be expressed in terms of dilatation and rotation parts by substituting 3.8 in 2.3 :

3.13
In 3.13 , the following parameter is used:

Trial Solutions
To proceed with the solution, the unbounded full space is divided in two domains, −∞ < z ≤ 0 and 0 < z < ∞, respectively.For each domain a trial solution Ansatz function is given.
For domain 1, 15 And for domain 2,

3.18
The superscripts 1 and 2 in 3.15 to 3.18 are used to denote quantities associated with the domains bounded by −∞ < z ≤ 0 and 0 < z < ∞, respectively.The 8 constants present at the trial solutions 3.15 to 3.18 , namely, A m and B m j with j 1, 2, 3 and m 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, A m / 0 and B m j / 0 with j 1, 2, 3 and m 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 α L and α S must be a positive value.The variables β and γ in the trial solutions are wave numbers corresponding, respectively, to the x and y Cartesian directions.
The additional criterion given in 3.12 , which establishes that the vector field Ω has zero divergence Ω i,i 0 , allows the three components of the displacement field U i to be uniquely determined from the four components of Δ and Ω n n 1, 2, 3 29 .Applying 3.12 in 3.16 and 3.18 yields

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 Analogously, 3.17 and 3.18 are substituted into 3.8 to obtain the displacement field of the domain 2:

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 : S e α S z e i βx γy , 3.27 2 βγ e α S z e i βx γy , 3.30 8 Mathematical Problems in Engineering 2 βγ e α S z e i βx γy , 3.31 2 β e α S z e i βx γy .

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 A m and B m j with j 1, 2, 3 and m 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 z 0 is prescribed.Using 3.21 to 3.26 continuity conditions gives rise to the following set of equations: e i βx γy 0, 4.1 e i βx γy 0, 4.2 e i βx γy 0.

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

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 x-y z 0 see Figure 1 a , given by the following expression in the Fourier wave number domain: The boundary conditions prescribe stress continuity at the interface in the tangential directions x and y, and a jump or discontinuity in the normal direction z due to the applied normal load p Z β, γ , leading to the expressions: x, y, z 0 0.

4.9
The solution of 4.9 leads to

4.12
The displacement field in the Fourier transformed wave number domain due to the vertical loading p Z β, γ is obtained by substituting 4.10 to 4.12 into 3.21 to 3.26 : z e i βx γy .

4.13
In this paper, U m ik indicates the displacement of a point at medium m in the direction i due to a loading in the direction k.

The Second BVP-Displacements Field due to the Tangential Load in Y 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 x-y z 0 see Figure 1 b , given by the following expression in the Fourier wave number domain: The boundary conditions prescribe stress continuity at the interface in the normal direction z and in the tangential directions x as well as a jump or discontinuity in the y direction due to the applied normal load p Y β, γ , leading to the expressions: The solution of the system 4.16 leads to 0.

4.17
The displacement field in the Fourier transformed wave number domain due to the transversal loading p Y β, γ is obtained by substituting 4.17 into 3.21 to 3.26 : e −α L z − e −α S z e i βx γy . 4.18

The Third BVP-Displacements Field due to the Tangential Load in X 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 x-y z 0 see Figure 1 c , given by the following expression in the Fourier wave number domain: The boundary conditions prescribe stress continuity at the interface in the normal direction z and in the tangential directions y as well as a jump or discontinuity in the x direction due to the applied normal load p X β, γ , leading to the expressions:

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, A

4.21
The solution of the system 4.21 leads to

Mathematical Problems in Engineering
The displacement field in the Fourier transformed wave number domain due to the transversal loading p X β, γ is obtained by substituting 4.22 into 3.21 to 3.26 : e −α L z − e −α S z e i βx γy . 4.23

Displacements Field in the Transformed Space
Equations 4.13 , 4.18 , and 4.23 express the component of displacement U m ik with distinct expressions for each media m 1 and m 2. In the following equations, a single expression for the components U ik is given for both media through a small adaptation of the coordinate z.In these equations, U ik indicates the displacement of a point of coordinates x, y, z of the full space in the direction i i x, y, z due to a loading applied in the direction k k x, y, z .
For the first problem with normal loading p Z : |z| e i βx γy .

4.24
For the second problem with tangential loading p Y : e −α L |z| − e −α S |z| e i βx γy .

4.25
And for the third problem with tangential loading p X :

Uniformly Distributed Traction Loadings in the Wave Number Domain
Consider the vertical harmonic load p Z x, y, z depicted in Figure 1 a .Let this load have a uniform intensity p Z distributed over a rectangular surface with dimensions 2A × 2B, placed at the interface plane x-y z 0 .In the original physical domains x, y, z this load can be written as

4.27
The loading can be transformed to the wave number domain p Z β, γ using the double integral Fourier transform: 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 p i β, γ i x, y, z 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.

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 x, y, z .
For the first BVP with normal excitation p

4.31
For the second BVP with tangential excitation p

4.32
For the third BVP with tangential excitation p

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.37
In 4.35 to 4.37 , the following notation is adopted:

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 2A × 2B in the plane x-y z 0 of the full space.These expressions must be evaluated numerically.

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.

Figure 1 :
Figure 1: Traction loadings for the three boundary value problems.

Equations 4 . 2 m 1 , 2 ,
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, can be obtained:

4 . 15 These traction boundary conditions 4 . 2 m 1 , 2 ,
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, can be obtained:
sin A 0 k β , s βx sin A 0 k β x A , s γ sin A 0 b 0 k γ , s γy sin 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 .
y, z 0 e −iβx dx e −iγy dy, 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 x, y, z .