Transient wave propagation phenomena at viscoelastic half-spaces under distributed surface loadings

This article analyzes the transient wave propagation phenomena that take place at 2D viscoelastic half-spaces subjected to spatially distributed surface loadings and to distinct temporal excitations. It starts with a fairly detailed review of the existing strategies to describe transient analysis for elastic and viscoelastic continua by means of the Boundary Element Method (BEM). The review explores the possibilities and limitations of the existing transient BEM procedures to describe dynamic analysis of unbounded viscoelastic domains. It proceeds to explain the strategy used by the authors of this article to synthesize numerically fundamental solutions or auxiliary states that allow an accurate analysis of transient wave propagation phenomena at the surface of viscoelastic half-spaces. In particular, segments with spatially constant and linear stress distributions over a halfspace surface are considered. The solution for the superposition of constant and discontinuous adjacent elements as well as linear and continuous stress distributions is addressed. The influence of the temporal excitation type and duration on the transient response is investigated. The present study is based on the numerical solution of stress boundary value problems of (visco)elastodynamics. In a first stage, the solution is obtained in the frequency domain. A numerical integration strategy allows the stationary solutions to be determined for very high frequencies. The transient solutions are obtained, in a second stage, by applying the Fast Fourier Transform (FFT) algorithm to the previously synthesized frequency domain solutions. Viscoelastic effects are taken into account by means of the elastic-viscoelastic correspondence principle. By analyzing the transient solution of the stress boundary value problems, it is possible to show that from every surface stress discontinuity three wave fronts are generated. (continued on next page...)


INTRODUCTION
In this introduction the main developments regarding the transient analysis of elastic and viscoelastic, bounded and unbounded domains by the Boundary Element Method is reviewed.
In spite of the many attempts to develop Finite Element Methodologies (FEM) able to perform linear dynamic analysis of unbounded elastic domains [10,11,25,49] the study of such domains requiring the satisfaction of the Sommerfeld radiation condition [61] is, "par excellence", the realm of the Boundary Element Method [8,9,26].The most striking advantages of the BEM concerning the dynamic analysis of unbounded domains is the natural fulfillment of the Sommerfeld radiation condition and the reduction of the problem dimension by one.To take full advantage of these properties the BEM must be formulated using an auxiliary state (a fundamental solution, a Green's function or similar state) that satisfies simultaneously the operator under analysis and the Sommerfeld radiation condition.An account of the main currently available strategies to perform the dynamic analysis of unbounded domains may be found in Mesquita and Pavanello [48].
The first publication on transient Boundary Element schemes for scalar wave propagation problems are due to Mansur and Brebbia [35,36].For two-dimensional elastic problems, a time stepping procedure was introduced by Mansur [34] as well as by Mansur and Brebbia [37].Three-dimensional transient elastic BE formulation was first presented by Karabalis and Beskos [29].In this article the authors used Stokes solution with a Heaviside time impulse function together with constant spatial and temporal discretizations.
As early as 1985, Antes has also produced a time-stepping algorithm for the 2D elastodynamic analysis [1].Antes and his co-workers continued the development of transient analysis and also produced strategies to coupled BEM-FEM problems [62,63].The introduction of a Dirac's Delta as a time excitation function in 3-D transient elastic analysis was due to Banerjee et al. (1986), [2].
From this point on, several variants for the transient analysis were introduced, aiming to improve accuracy and stability of the time domain procedures.A 3-D non-singular transient formulation in space and time based on Stokes fundamental solution and presenting a Heaviside window as a temporal excitation function was given by Coda and Venturini in 1995, [14].A smoother polynomial temporal excitation function was introduced by Coda and Venturini in 1996, [15].A multi-domain strategy for 3D transient problems able to couple bounded and unbounded domains described, respectively, by the BEM and FEM was presented by Coda, Venturini and Aliabadi in 1999, [17].
Boundary element analysis of transient viscoelastic and non-linear plastic phenomena was reported by Coda and Venturini [16].They used a static fundamental solution together with domain integration, performed by means of cells, to include inertia effects as well as non-elastic phenomena, i.e., viscoelasticity and plasticity.The most evident limitation of this procedure is that it is limited to bounded domains.It cannot model unbounded domains or fulfill the Sommerfeld radiation condition.Rizos and his co-workers followed a different path.Rizos and Karabalis [52] developed a B-spline fundamental solution for the transient 3D elastic isotropic problem.The numerical solution of the resulting BIE was based on higher order temporal and spatial discretizations.Using the B-Spline fundamental solution, or auxiliary state, Rizos and his co-workers developed strategies to coupled the BEM with the FEM or with rigid bodies to allow for transient analysis of 3D elastic bonded and unbounded domains [50,51,53,54].Carrer and Mansur, [12] determined stresses and velocities in 2D transient dynamic analysis.
All previously described transient formulations are time stepping algorithms.Within each new time step a new system matrix is generated and must be stored.For analysis consisting of large time periods with many time steps these procedures demand high storage capacities.
The previous paragraphs described significant moments regarding the development of transient boundary element methodologies mainly for an ideally elastic media.Next the research efforts to include viscoelastic phenomena into the transient BEM analysis is addressed.It is well known that many materials exhibit viscoelastic properties which play an important role in the continuum dynamic response.The soil [27,33], concrete [60] and polymers [28] are typical examples of these materials.
To perform a transient viscoelastic analysis by the boundary element method, in which only the boundary of the domain would need to be discretized and in which the Sommerfeld radiation condition would be satisfied, a viscoelastic transient fundamental solution or equivalent auxiliary state would be necessary.Unfortunately, there are many viscoelastic material models and the continuum also presents a large spectrum of viscoelastic properties [20].So it is not possible to synthesize "the one transient viscoelastic fundamental solution".For every viscoelastic model, or for every measured set of viscoelastic parameters, a new fundamental solution would be necessary.This complexity of models and material behavior is what makes the solution of transient viscoelastic problems by the BEM a challenging problem.
Although there is no "single fundamental solution" for transient viscoelastic analysis, researchers have spent a lot of effort to develop strategies to face this challenge.There are two main approaches to model linear viscoelastic phenomena.The first approach is to use differential operators to describe the viscoelastic behavior [13,20].In this case the viscoelastic model is usually given in the frequency domain.The second approach is to make use of the so-called relaxation functions [13,20].The relaxation functions are viscoelastic properties given in the time domain.A connection between the description by differential relations and by the relaxation functions can be established [33].
One classical approach to solve transient problems is to use the Laplace transform.The time domain equations are transformed into the frequency domain.In this transformed domain, the elastic-viscoelastic correspondence principle [13] may be used to exchange the elastic constitutive parameters by its viscoelastic counterparts.The inverse Laplace transform, back to the original time domain, will yield the transient solution of the viscoelastic problem.This has been used by many authors to yield transient solutions, see for instance [64].There are at least two limitations in this approach that must be mentioned.The first one is the difficulty to include time varying boundary conditions in the analysis [38].The second is that the numerical algorithms available to perform the inverse Laplace domain are not able to furnish very small time steps and sometimes fail to give a very detailed description of the transient phenomena.
In 1992, Gaul, Schanz and Fiedler reported an article in which the transient viscoelastic analysis was performed by the Laplace transform method and by the use of elastic-viscoelastic correspondence principle, [23].The modeling capabilities of the differential operator approach was enhanced by including fractional order derivatives in the model.This allowed a better matching of the viscoelastic properties with smaller number of parameters.They reported a time domain viscoelastic fundamental solution for the Maxwell model [23].This is an important result, but it is well known that the Maxwell model is not appropriate to model solids [20].
Another alternative was proposed by Gaul and Schanz in 1997 [22].The analysis employed the elastic-viscoelastic correspondence principle with fractional differentiation to describe the constitutive equations.But the implementation of the viscoelasticity was provided within each time step, in the Laplace domain.The results for distinct inverse algorithms, including Talbot and Durbin and Crump, were presented.
A new approach to describe transient problems by the BEM was introduced by Schanz and Antes.They developed the Operational Quatrature Method proposed by Lubich [31,32] to be incorporated into the realm of the BEM for transient analysis.They named their approach Convolution Quadrature Method (CQM).The first work was concerned with the transient response of elastic problems [57].At the same year the work was extended to include viscoelastic effects [58].This very promising strategy does not require a time domain fundamental solution and avoids very sensitive fundamental solutions, which may yield not very stable results.This approach has been extended to distinct class of problems such as poroelasticity [56].Later, Schanz applied successfully the CQM to quasi-static viscoelastic problems [59].
A different research line was pursued by Coda and his co-workers.Using constitutive differential equations for viscoelastic models Coda and his co-workers presented a quasi-stationary formulation based on cell discretization for the Kelvin [43] and for the Boltzmann [42] viscoelastic models.This formulation leads to a matrix differential equation, the solution of which does not require the storage of a new system matrix at every time step.It also allowed the inclusion of time varying boundary conditions.But since domain discretization through cells was required, no unbounded domains could be modeled.
In a later work they improved their formulation to eliminate the domain integrals, resulting in a quasi-static viscoelastic formulation for Kelvin and Boltzmann models.The twodimensional case was treated in [41] and the 3-D case was discussed in the paper of 2003 [42].Through the use of differential viscoelastic relations they did avoid the use of relaxation functions and the convolutional solution procedure.This improved procedure made possible the analysis of unbounded domains.Here the resulting matrix differential equation could be solved by a time marching scheme without the need to store system matrices at every time step.They also presented a new Finite Element formulation for viscoelasticity [39,40] and a coupling procedure for BEM-FEM [41].
The works by Coda and his co-workers have many advantages, as has already been mentioned.The main disadvantages are related to the quasi-static character of the problems that can be handled, when analyzing unbounded domains.That means no inertia terms can be considered, hence no wave propagation problems may be described.The other aspect to be commented is that for every viscoelastic model a new formulation must be written.The authors show the steps to be followed for other viscoelastic models, but the work must still be done.
During the last 20 years the Brazilian research group, to which most of the authors of this article belong, have been involved in the development of Green's functions and auxiliary distributed load solutions for full-spaces, half-spaces, and layered media.Initially frequency domain solutions for 2-D and 3-D viscoelastic half-spaces were developed [44].A solution for a layer on rigid bottom or over a half-space was developed by Romanini [55].Green's functions for horizontally multi-layered 3D soil profiles were presented by Barros, [7].Green's functions for transverse isotropic full-spaces and half were synthesized by Barros and Mesquita, [5].
The synthesized solutions were concentrated load solutions, known as Green's functions, and non-singular load solutions [4] as well as singular-ended distributed load solutions, [5].These frequency domain non-singular solutions were incorporated in a direct version of the BEM [6] as well as in indirect BE formulations [3].
The common mathematical basis for the previously mentioned synthesized frequency domain solutions was the Fourier integral transform.Viscoelastic effects were taken into account by the elastic-viscoelastic principle, which for the stationary, frequency domain, case is rather straight forward.Efficient numerical Fourier inversion strategies were developed leading to frequency domain solutions which can be accurately obtained at very large distances or at very high frequencies [30], [18].
The main features of this semi-numerical Fourier integral approach can be listed.First, the numerically synthesized transient Green's Functions or auxiliary states were originally obtained in the Frequency Domain, so any linear viscoelastic model, which has a frequency domain representation, may be included in the analysis.This includes experimentally determined viscoelastic constitutive parameters.
Second, the synthesized transient solutions are dynamical, that means that all inertia effects are already included in the solution.This type of solution may be incorporated in a boundary integral scheme which does not require domains discretizations and neither cell integrations.The resulting transient BE scheme is a convolutional scheme in which time varying boundary conditions may be incorporated without any difficulty.
Third, the reported methodology is able to determine a viscoelastic isotropic and anisotropic transient solutions, which could not be obtained by the other reviewed strategies.The larger limitation of this approach is that the synthesis of frequency domain solutions may be computationally expensive.But the solutions may well be determined using distributed computation on a cluster of machines or on combined CPU-GPU systems [18,19].The synthesized 2D solutions were successfully used to determine short time and long time response of layered unbounded viscoelastic media [46].
Based on the previous developments, the idea of the present article is to show that the viscoelastic frequency domain solutions numerically synthesized by the authors can describe transient phenomena at very small time-steps and can accurately depict the wave propagation phenomena, including sharp wave fronts.The viscoelastic transient response obtained can be correlated to the arrival of the distinct wave front existing at the transient half-space response.For the best of the authors' knowledge such a detailed recovery of the wave propagation phenomena was not reported in the reviewed literature.
The analysis performed in this article is based on the transient solution of two stress boundary value problems of visco-elastodynamics. Spatially distributed constant and linear traction loads at the half-space surface are considered.The spatially distributed loads may by applied as a Dirac delta temporal loading or through Heaviside type time windows with various durations.
The stress boundary value problems are solved in the frequency domain with the aid of the Fourier integral transform.A numerical integration strategy is developed, which allows half-space Frequency Response Functions (FRFs) to be synthesized for very high frequencies.
The application of the FFT algorithm on these high frequency solutions will yield time domain solutions with very small time steps.These numerical solutions are capable of recovering, in great detail, the wave propagation phenomena that take place at the half-space surface and at its vicinity.
It will be shown that the transient response depends strongly on the number of load elements applied at the half-space surface as well as on the type of temporal excitation.Dirac delta excitation type will result in displacement oscillations every time a wave front, generated by a stress discontinuity, impinges the observation point.These phenomena occur for both constant and discontinuous but also for linear and continuous load elements.The number of jumps or oscillations at the observation point is related to the number of surface stress discontinuities.
In particular the article shows that for half-spaces excited by spatially constant stress distributions applied at time intervals of finite duration, the transient response will depend largely on the loading duration.If the loading time span is smaller than the time required for all wave fronts to overcome the distance between the observation point and the farthest surface stress discontinuity, then the transient response will still present some oscillations.A smoother transient behavior is obtained for loads applied at time windows of larger durations.The relative importance of spatial discretization and temporal excitation on the transient response will be discussed in detail.

Equations of motion for the viscoelastic continuum
The equations of motion governing the transient behavior of a 2D plane strain homogeneous, isotropic and viscoelastic continuum expressed in terms of the displacement components u i (t) (i=x,z ) are given by [13]: In equation ( 1), the relaxation functions for the isotropic continuum are given by µ(t) and λ(t), whereas the continuum density is ρ and b i represents the body forces acting in the domain.Neglecting the body forces, b i =0, and assuming a harmonic time dependence u i (x, t) = ūi (x) e iωt , the frequency domain counterpart of equation ( 1) can be written as [13]: In equation ( 2), the variable ω is the circular frequency, λ * (ω) and µ * (ω) are the complex and frequency dependent Lamé constants containing the linear viscoelastic model of the continuum.For linear viscoelastic models the complex Lamé contants may be expressed as: and The storage and loss moduli of the continuum, defined in equations ( 3) and ( 4), are, respectively, λ * S (ω), µ * S (ω)and λ * L (ω), µ * L (ω).The damping factors η λ (ω)and η µ (ω) shown in equations ( 3) and ( 4) are defined as:

Boundary conditions
In this article, two sets of boundary conditions will be used to synthesize the necessary auxiliary solutions.Both auxiliary problems are stress boundary value problems in which vertical, spatially constant and a spatially linear traction distributions are prescribed at a half-space surface.The surface stress boundary conditions are show in figure 1.
The mathematical formulation for the linear traction distribution of width 2a applied at the half-space surface Γ S is: In equation ( 6), t z1 and t z2 represent the vertical tractions amplitude, respectively, at points x 1 and x 2 .The constant distribution is included with t z1 = t z2 .

DISPLACEMENT SOLUTION IN THE FREQUENCY DOMAIN
The described stress boundary value problems have been solved in the frequency domain with the aid of the Fourier integral transform.The final solution step involves a numerical evaluation of the inverse Fourier transform with respect to the pair (x,β), in which x represents the space and β the wave number domain variables.For the case of the vertical displacement component u zz (x, z, ω), originated due to a vertical surface traction distribution t zi (x, z = 0, ω) (i = 1, 2), the solution is given by the following expression [4] 2 ) ( βa cos(βa)−sin(βa) )] exp(iβx)dβ (7) A detailed expression for the kernel H zz (β, ω) of equation ( 7) can be found in Appendix 1.An important part of this research is the development of an integration strategy for equation (7) based on the Longman [30] algorithm for improper integration combined with an adaptive quadrature scheme [21].The strategy is able to perform accurate integrations of equations with structures resembling that of (7) for very high values of the circular frequency ω.This integration strategy may be dramatically accelerated by the use of graphics hardware [18].
A typical stationary displacement solution u zz is shown in figures 2 and 3.For these solutions, the half-space was loaded with a spatially constant stress distribution t z1 = t z2 = 1 Figure 2a shows the real part of the displacement solution u zz for low values of the dimensionless frequency A 0 .Figure 2b depicts the same solution for very high frequency values.The low frequency behavior of the imaginary part of u zz can be seen in figure 3a. Figure 3b shows the high frequency content of the same imaginary part in log scale.The low frequency content of these solutions have already been validated in [4].Figures 2 and 3 demonstrate the capability of the developed integration strategy to determine the half-space solutions at very high frequencies.It should also be noted that for most stationary civil engineering applications the analysis is limited to low frequencies, A 0 <10 [24].The high frequency results are determined to ensure transient responses with very short time steps by the FFT algorithm [45,47].

Superposition of linear elements in the frequency domain
The stationary displacement at point x i , u zzi (ω) = u zz (x i , ω), due to a series of n linear elements acting simultaneously at the half-space surface, as shown in figure 4, can be determined by superposition of the solution given in equation ( 7) In expression (8), t zi (i = 1, n) are the traction amplitudes the load points with coordinates x i (i = 1, n).The coordinate of the center of the i-th loading element is given by x 0i = (x i + x i+1 )/2.The coefficients I k zz (x − x 0i , z, ω)with (k = 1, 2) and (i = 1, n − 1) can be regarded as influence functions that are to be determined numerically from equation (7).

TRANSIENT SOLUTIONS
The transient solutions are obtained by applying the Fast Fourier Transform (FFT) algorithm with respect to the pair (ω, t) or (A 0 , t) to the previously synthesized frequency domain solutions.Prior to being submitted to the FFT algorithm, the signal of the stationary solution u zz (ω) must undergo a processing as described in [45].The signal processing involves the determination of a cut-off frequency, the construction of a smooth transition to vanishing values above the cut-off frequency, the addition of zeros at the high end of the spectrum as well as the elimination of the rigid body component present in the 2D half-space solution.
If the applied external traction excitation remains with constant unit amplitude for all circular frequencies, t zj (ω) = 1, then the stationary solution measured at point x i due to an excitation applied at point x j is the frequency response function (FRF) of the continuum for this stress boundary value problem, u zz (x i , x j , t zj = 1, ω) = H zz (x i , x j , ω).It corresponds to Latin American Journal of Solids and Structures 9(2012) 453 -473 the spectrum of the transient continuum response due to a temporal Dirac delta excitation of unit amplitude.In order to obtain transient solutions for a general temporal excitation f z (t), it is first necessary to obtain the spectrum or the Fourier transform of the forcing function, F z (ω).The product of the frequency response function H zz (ω) with the excitation spectrum F z (ω) will provide the frequency displacement solution to the general transient loading, u zz (ω) = H zz (ω)F z (ω).Finally, the application of the FFT algorithm to the frequency domain signal u zz (ω) will lead to the transient solution u zz (t).It will be shown that this procedure results in accurate solutions which are able to capture all the essential features of the underlying transient wave propagation phenomena.
One important issue of the Discrete Fourier Transform (DFT) algorithm is the relation ∆A 0 = 2π/T max between the frequency step ∆A 0 = A 0k+1 − A 0k (k = 0, n max ) and the maximum time reachable by the transient process, T max .Equally important is the relation ∆t = 2π/A 0 max , connecting the time step ∆t and the maximum sampling frequency A 0 max .This implies that to reproduce transient phenomena with very small time steps ∆t the maximum sampling frequency A 0 max must be correspondently high.
In this section, the transient solution of several stress boundary value problems (SBVP) will be investigated.The temporal excitations to be considered are the Dirac delta and time windows of various finite durations.The spatial traction distributions at the half-space surface are the constant element, the linear element as well as superposition of constant and linear elements.
With the values of the parameters a the velocity of the compression (c P ), shear (c S ) and Rayleigh (c R ) waves of the elastic continuum are, respectively: Considering a temporal discretization with time step ∆t, the numerical realization of the Dirac delta δ N (t) applied to the time instant t=t a is implemented according to the following simplifying rule:

Transient response to a Dirac Delta. Constant element
The frequency domain solutions shown in figures 2 and 3 were submitted to the signal processing treatment described above and to the FFT algorithm.Initially, a Dirac temporal loading applied at t a = 0 [s] is considered.Figure 5 shows the long term transient vertical displacement response u zz (t) for the point with coordinates x = z = 0.After an abrupt initial displacement, the solution decreases monotonically towards the zero value of the initial undisturbed solution.
Although this solution is consistent with the expected behavior, in the initial time instants there are complex phenomena taking place which deserve a more detailed analysis.Figure 6 shows in detail the initial time instants of the solution due to a constant spatial loading.It corresponds to the initial time steps of the solution depicted in figure 5.It can be seen that, at the very initial time steps, there is a massive displacement which takes place

Transient response to a Dirac delta. Superposition of constant elements
To substantiate this analysis, two more numerical investigations are conducted.The first is the superposition of three constant elements of equal traction amplitude, t z1 =t z2 =t z3 =t z4 =1 [N].For this example, the element half-width is a = 1/3 [m].This case is depicted in figure 7. The coordinates of the element edges are given by x i (i = 1, 4).The central point of the elements are designated by x 0j (j = 1, 3).The results, presented so far, point out that the procedure developed to obtain the transient displacement solution of half-space problems subjected to stress loadings is very accurate, produces very small time steps, and allows to recover the complex wave propagation phenomena generated by a temporal Delta loading.Furthermore, it also permits to show that under a distributed spatial traction and a temporal Dirac delta loading, the maximum displacement amplitude at any surface point is generated by the passing of the Rayleigh waves.The arrival of the Rayleigh waves at a given observation point causes a sudden and large increase in the displacement amplitude.

Transient response to a time loading of finite duration. Spatially constant stress distribution
In this section, the transient half-space displacement response due to a time window of finite duration is analyzed.The surface stress boundary conditions have the following time dependency: In equation (10), H(t) represents the Heaviside or unit step function.The time window is characterized by its duration, ∆t w = t f − t i , where the initial and final window times are given by t i and t f , respectively.
A single element of width 2a = 2 [m] and constant spatial surface stress distribution t z1 = t z2 = 1 [N] is subjected to temporal excitations of finite duration.Figure 12 shows the transient displacement response for the point at the loading center, u zz (x = z = 0, t), for time window of distinct durations: ∆t w1 = 0.464 [s], ∆t w2 = 0.928 [s], ∆t w3 = 1.392 [s] and ∆t w4 = 2, 321 [s].For these cases the initial time is t i = 0 [s].The response due to a Dirac delta time excitation is also shown.Due to the fact that larger time excitation durations introduce larger momentum in the half-space and, consequently, produce larger response amplitudes, the results shown in figure 12 are normalized by the maximum amplitude of each response, u zz (t, ∆t w )/u zz max (∆t w ).
Figure 12 shows that the temporal windows of finite duration produce displacement patterns which are very different from those of the temporal Dirac delta excitation.During the time that the excitation is applied, t i < t < t f , a massive displacement takes place with the largest velocity of the medium, c P .This massive displacement overrides the displacement fluctuations originally caused by the arrival of the wave fronts generated at the stress discontinuities.There are still some small fluctuations that can be associated to the arrival of the wave fronts, but the strong amplitude peaks stemming from the passage of the Rayleigh wave are no longer present.
In analyzing figure 12 is important to notice that if the time window is smaller then the time required for a given wave front to arrive at the observation point, this and the subsequent wave fronts cause some abrupt changes or oscillations in the transient response.On the other hand, if the time window duration ∆t w is larger than the time required for the last wave front to impinge the observation station, the solution is smoothed by the integrated displacement response stemming from the Heaviside type excitation.

CONCLUDING REMARKS
In the present work a review of existing strategies to solve transient elastic and viscoelastic problems for bounded and unbounded domains with the Boundary Element Method is given.Furthermore, transient stress boundary value problems of (visco)elastodynamics have been numerically solved.Spatially constant and linear vertical stress distributions applied at the surface of an isotropic 2D half-plane have been considered.Initially the problem is solved in the Latin American Journal of Solids and Structures 9(2012) 453 -473 frequency domain.Accurate Frequency Response Functions (FRFs) for the stress boundary value problem are synthesized for very large frequencies.These FRFs are used in conjunction with the Fast Fourier Transform (FFT) algorithm to yield transient half-space displacement solutions to temporal excitations ranging from a Dirac delta to time windows of various finite durations.The presented procedure is able to depict, accurately, the complex wave propagation phenomena that take place in the transient half-space solution.
Considering a temporal Dirac delta excitation, it is shown that 3 wave fronts are generated at those surface points where stress discontinuities are present.These wave fronts can be associated to compression, shear and surface Rayleigh waves.By solving and analyzing spatially constant and linear surface stress distributions and also a superposition of these loadings, it has been possible to correlate the transient displacements at any observation point on the half-space surface with the arrival of the mentioned wave fronts.The half-space displacement solution due to a temporal Dirac delta excitation presents various oscillations.These oscillations can be correlated to the arrival of the wave fronts at the considered observation point.The amplitude and the time instants of the oscillations depend on the half-space surface discretization scheme.For both constant and linear elements the larger oscillations will occur when the Rayleigh wave fronts impinge the observation point.
When constant temporal forces are applied within time windows of finite durations, the resulting displacement pattern is completely different from the responses due to a Dirac time loading.If the time window is smaller then the time required for any wave front to arrive at the observation point, this and the subsequent wave fronts cause oscillations in the transient response.If the time window duration is larger than the time required for the last wave front to impinge the observation station the solution is smoothed by the integrated displacement response stemming from the finite duration time window.
The accuracy of the transient solution strategy obtained in the very simple superposition schemes do suggest that these numerically synthesized solutions may be used in more general problems within the context of the Boundary Element Method.

APPENDIX: EXPRESSION FOR THE KERNEL H ZZ (β)
Considering the wave number β relating the circular frequency ω and the continuum wave velocity c, β= ω/c, the integration kernel H zz (β) of equation ( 7) is given by: with (12) And

Figure 1
Figure 1 Stress boundary conditions at the half-space surface [N].The solution was determined at the central point of the loading surface with coordinates x = z = 0.0 The other parameters are: load half-width a = 1 [m], density ρ = 1 [kg/m 3 ], Poisson ratio ν = 0.33 and the Lamé constant µ * S (ω) = µ = 1 [N/m 2 ].The case of constant Latin American Journal of Solids and Structures 9(2012) 453 -473 hysteretic viscoelastic model with damping factor η µ = η λ = 0.01 was considered.The real and imaginary parts of the displacement solution are expressed as function of dimensionless frequency parameter A 0 = ω a/c S in which c S = (µ/ρ) 1/2 is the shear wave velocity of the elastic continuum.

Figure 2 Figure 3
Figure 2 Real part of the stationary displacement solution uzz-a) low frequency behavior, -b) high frequency behavior

Figure 4
Figure 4 Scheme for the superposition of linear solutions

Figure 6
Figure 6 Short term displacement response uzz(x = z = 0, t)and arrival of wave fronts

Figure 7
Figure 7 Superposition of 3 constant elements of equal traction amplitude Figure 8 shows the transient response measured at the origin of the coordinate system x = x 02 .The wave fronts, starting at the edges x = x 1 and x = x 4 , have to travel the distance |x 1 − x 02 | = |x 4 − x 02 | = 3a to reach the observation point x = x 02 .The time required for the distinct waves (P,S,R) to propagate the distances 3a are t P (3a) = 0.5 [s], t S(3a) = 1.0 [s], and t R(3a) = 1.07 [s], respectively.These time instants are marked in figure 8.The solution pattern of the former example is again recognizable, see fig.6.After the initial displacement with velocity c P , the solution stays in a constant plateau until the remaining two wave fronts (S,R) impinge the observation point.

Figure 9 Figure 8 Figure 9
Figure 8 Half-space transient response at point x 02 due to 3 equal constant loads.Dirac delta excitation

LatinFigure 10
Figure 10 Three linear loading elements at the half-space surface

Figure 12
Figure 12 Transient response uzz(t, ∆tw) due to distinct windows excitations with distinct durations ∆tw.