Acessibilidade / Reportar erro

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

Abstract

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 in uence 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. The displacement velocity of these wave fronts can be associated to compression, shear and Rayleigh waves. It is shown that the half-space transient displacement solutions present abrupt jumps or oscillations which can be correlated to the arrival of these wave fronts at the observation point. Such a detailed analysis connecting half-space transient responses to the wave propagation fronts in viscoelastic half-spaces have not been reported in the reviewed literature.

half-space; transient wave propagation; elastodynamics; viscoelasticity


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

E MesquitaI,* * Author email: euclides@fem.unicamp.br ; H AntesII; LH ThomazoIII; M AdolphIV

IDepartment of Computational Mechanics - DMC/FEM,State University at Campinas - UNICAMP,C.P. 6122 Campinas, Brazil

IIInstitut fr Angewandte Mechanik, Technische Universitt Braunschweig, Spielmannstr. 11, D-38106 Braunschweig, Germany

IIIDepartment of Computational Mechanics - DMC/FEM, State University at Campinas - UNICAMP, C.P. 6122 Campinas, Brazil

IVDepartment of Computational Mechanics - DMC/FEM, State University at Campinas - UNICAMP, C.P. 6122 Campinas, Brazil

ABSTRACT

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 in uence 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. The displacement velocity of these wave fronts can be associated to compression, shear and Rayleigh waves. It is shown that the half-space transient displacement solutions present abrupt jumps or oscillations which can be correlated to the arrival of these wave fronts at the observation point. Such a detailed analysis connecting half-space transient responses to the wave propagation fronts in viscoelastic half-spaces have not been reported in the reviewed literature.

Keywords: half-space, transient wave propagation, elastodynamics, viscoelasticity

1 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, 48] 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 [50].

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 two-dimensional case was treated in[41] and the 3-D case was discussed in the paper of 200342 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.

2 STATEMENT OF THE PROBLEM

2.1 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 ui (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 bi represents the body forces acting in the domain. Neglecting the body forces, bi=0, and assuming a harmonic time dependence ui (x, t) = i (x)eiω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, (ω), (ω) and (ω), (ω). The damping factor ηλ (ω) and ηµ (ω) shown in equations (3) and (4) are defined as:

2.2 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), tz1 and tz2 represent the vertical tractions amplitude, respectively, at points x1 and x2. The constant distribution is included with tz1 = tz2 .

3 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 zz (x, z, ω), originated due to a vertical surface traction distribution zi (x,z = 0,ω) (i = 1,2), the solution is given by the following expression [4]

A detailed expression for the kernel

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

zz is shown in figures 2 and 3. For these solutions, the half-space was loaded with a spatially constant stress distribution tz1 = tz2 = 1 [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/m3], Poisson ratio ν = 0.33 and the Lamé constant (ω) = µ = 1 [N/m2 ]. The case of constant 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 A0 = ω a/cS in which cS = (µ/ρ)1/2 is the shear wave velocity of the elastic continuum.

Figure 2a shows the real part of the displacement solution

zz for low values of the dimensionless frequency A0. Figure 2b depicts the same solution for very high frequency values. The low frequency behavior of the imaginary part of 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, A0 < 10 [24]. The high frequency results are determined to ensure transient responses with very short time steps by the FFT algorithm[45, 47].

3.1 Superposition of linear elements in the frequency domain

The stationary displacement at point xi, zzi (ω) = zz (xi, ω), 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), tzi (i = 1,n) are the traction amplitudes the load points with coordinates xi (i = 1,n). The coordinate of the center of the i-th loading element is given by x0i = (xi + xi + 1 )/2. The coefficients (x - x0i , 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).

4 TRANSIENT SOLUTIONS

The transient solutions are obtained by applying the Fast Fourier Transform (FFT) algorithm with respect to the pair (ω, t) or (A0, t) to the previously synthesized frequency domain solutions. Prior to being submitted to the FFT algorithm, the signal of the stationary solution 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,

zj (ω) = 1, then the stationary solution measured at point xi due to an excitation applied at point xj is the frequency response function (FRF) of the continuum for this stress boundary value problem, zz (xi, xj, zj = 1, ω) = zz (xi, xj, ω). It corresponds to 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 fz (t), it is first necessary to obtain the spectrum or the Fourier transform of the forcing function, z (ω). The product of the frequency response function zz (ω) with the excitation spectrum z (ω) will provide the frequency displacement solution to the general transient loading, zz (ω) = zz (ω)z (ω). Finally, the application of the FFT algorithm to the frequency domain signal zz (ω) will lead to the transient solution uzz (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 ΔA0 = 2π/Tmax between the frequency step ΔA0 = A0k + 1 - A0k (k = 0, nmax) and the maximum time reachable by the transient process, Tmax. Equally important is the relation Δt = 2π/A0max, connecting the time step Δt and the maximum sampling frequency A0max. This implies that to reproduce transient phenomena with very small time steps Δt the maximum sampling frequency A0max 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 = 1 [m], (ω) = µ = 1 [N/mm2 ], ν = 0.33 and ρ = 1 [kg/m3 ] the velocity of the compression (cP), shear (cS) and Rayleigh (cR) waves of the elastic continuum are, respectively: cP = 2 [m/s], cS = 1 [m/s] and cR ≈ 0.932 [m/s].

Considering a temporal discretization with time step Δt, the numerical realization of the Dirac delta δN (t) applied to the time instant t = ta is implemented according to the following simplifying rule:

4.1 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 ta = 0 [s] is considered. Figure 5 shows the long term transient vertical displacement response uzz (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 with the largest velocity of the continuum, cP. After a very short time t < 0.1 [s], the solution reaches a plateau and is kept constant up to the time t ≈ 0.5 [s]. This is precisely the time required by the compression wave to cover the distance from the loading edges x = |a| to the center point x = 0, tP(a) = 0.5 [s]. Then the solution reaches and stays at a new plateau until other disturbances arrive. These new disturbances are, respectively, the shear cS and the Rayleigh cR wave fronts also originated at the traction discontinuous loading edges, x = |a|. The arrival times of these wave fronts, respectively, tS(a) = 1.0 [s] and tR(a)≈ 1.07 [s], are marked in figure 6. It can be seen that the maximum displacement takes place when the Rayleigh waves, stemming from the traction discontinuities at the edges x = |a|, arrive and are superposed at the observation point x = z = 0. After the passage of the Rayleigh waves the solution decays monotonically.

4.2 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, tz1 = tz2 = tz3 = tz4 = 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 xi (i = 1,4). The central point of the elements are designated by x0j (j = 1,3).

The traction discontinuities are located at the points x = x1 and x = x4 . At these points the wave fronts are generated. Figure 8 shows the transient response measured at the origin of the coordinate system x = x02 . The wave fronts, starting at the edges x = x1 and x = x4 , have to travel the distance |x1 - x02 | = |x4 - x02 | = 3a to reach the observation point x = x02 . The time required for the distinct waves (P,S,R) to propagate the distances 3a are tP(3a) = 0.5 [s], tS(3a) = 1.0 [s], and tR(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 cP, the solution stays in a constant plateau until the remaining two wave fronts (S,R) impinge the observation point.

Figure 9 shows the transient response at point x01 , which has two distinct distances from the loading edges: |x1 - x01 | = a and |x4 - x01 | = 5a. The time required for the waves (P,S,R) to propagate the distances a and 5a are tP(a) = 0.167 [s], tS(a) = 0.333 [s], tR(a) = 0.358 [s], tP(3a) = 0.83 [s], tS(3a) = 1.67 [s] and tR(3a) = 1.79 [s], respectively. These time instants are again marked in figure 9. In this picture the maxima displacements stemming from the passage of the two Rayleigh wave fronts at the observation point x01 are clearly recognizable.

Next, the response to a spatial linear loading, as shown in figure 10, is analyzed. For this case the largest load intensity, at point x4 , has the amplitude tz4 = 1[N]. For this boundary condition, the only stress discontinuity is at point x = x4 . Accordingly, only from this point the wave fronts should be generated. Figure 11 depicts the transient response uzz (x1, t) at point x = x1. The distance from the load discontinuity to the observation point is |x1 - x4 | = 6a = 2m. The time required for the different wave fronts to overcome this distance is tP(6a) = 1.0 [s], tS(6a) = 2.0 [s], and tR(6a) = 2.15 [s], respectively. These time instants are indicated in figure 11. As expected, the only wave fronts noticeable in the transient response are those corresponding to the time required by the waves to travel the distance |x1 - x4 | = 6a. In this solution the passing of the Rayleigh wave front is again vividly observed.

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.

4.3 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, Δtw = tf - ti, where the initial and final window times are given by ti and tf, respectively.

A single element of width 2a = 2 [m] and constant spatial surface stress distribution tz1 = tz2 = 1 [N] is subjected to temporal excitations of finite duration. Figure 12 shows the transient displacement response for the point at the loading center, uzz (x = z = 0,t), for time window of distinct durations: Δtw1 = 0.464 [s], Δtw2 = 0.928 [s], Δtw3 = 1.392 [s] and Δtw4 = 2,321 [s]. For these cases the initial time is ti = 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, uzz (ttw)/uzzmaxtw).

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, ti < t < tf, a massive displacement takes place with the largest velocity of the medium, cP. 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 Δtw 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.

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

Acknowledgments The research leading to this article has been supported by the brazilian funding agencies Capes, Cnpq and Fapesp. This is gratefully acknowledged.

Received 16 Apr 2012

In revised form 16 May 2012

6 APPENDIX 1: EXPRESSION FOR THE KERNEL Hzz ( β )

Considering the wave number β relating the circular frequency ω and the continuum wave velocity c, β= ω/c, the integration kernel Hzz (β) of equation (7) is given by:

with

And

  • [1] H. Antes. A boundary element procedure for transient wave propagations in two-dimensional isotropic elastic media. Finite Elements Analysis and Design, 1:313322, 1985.
  • [2] P.K. Banerjee, S. Ahmad, and G.D. Manolis. Transient elastodynamic analysis of 3-d problems by boundary element method. Earthquake Engineering and Structural Dynamics, 14:933949, 1986.
  • [3] P.L.A. Barros and E. Mesquita. Singular and higher order elements for accurate interface stresses in dynamic soil structure interaction. In Proceedings ASCE 12th Engineering Mechanics Division Conference, pages 146149, La Jolla, California, May 17-20 1998.
  • [4] P.L.A. Barros and E. Mesquita. Elastodynamic green's functions for orthotropic plane strain continua with inclined axis of symmetry. International Journal for Solids and Structures, 36:47674788, 1999.
  • [5] P.L.A. Barros and E. Mesquita. Singular-ended spline interpolation for 2d boundary element analysis. International Journal for Numerical Methods in Engineering, 47(5):951967, 2000.
  • [6] P.L.A. Barros and E. Mesquita. On the dynamic interaction and cross-interaction of 2d rigid structures with or thotropic elastic media possessing general principal axes orientation. Meccanica, 36(4):367378, 2001.
  • [7] R.M.B Barros. Funações de Green e de Influência para meios visco-elásticos transversalmente isotrópicos no domínio da frequância (Green and influence functions for transverse isotropic media in the frequency domain). PhD thesis, University of Campinas, Unicamp, 2001.
  • [8] D.E. Beskos. Boundary element methods in dynamic analysis. Applied Mechanics Reviews, 40(1):0123, 1987.
  • [9] D.E. Beskos. Boundary element methods in dynamic analysis: Part ii (1986-1996). Applied Mechanics Reviews, 50(3):149197, 1997.
  • [10] P. Bettess. Infinite Elements Penshaw Press.
  • [11] P. Bettess. Infinite elements. International Journal for Numerical Methods in Engineering, 13(1):6364, 1977.
  • [12] J.A.M. Carrer and W.J. Mansur. Stress and velocity in 2D transient elastodynamic analysis by the boundary element method. Engineering Analysis with Boundary Elements, 23(3):233245, 1999.
  • [13] R.M. Christensen. Theory of Viscoelasticity. Academic Press, NY, 2nd. ed. edition, 1982.
  • [14] H.B. Coda and W.S. Venturini. Non-singular time-stepping bem for transient elastodynamic analysis. Engineering Analysis with Boundary Elements, 15:1118, 1995.
  • [15] H.B. Coda and W.S. Venturini. Further improvements on three dimensional transient bem elastodynamic analysis. Engineering Analysis with Boundary Elements, 17:231243, 1996.
  • [16] H.B. Coda and W.S. Venturini. Dynamic non-linear stress analysis by the mass matrix bem. Engineering Analysis with Boundary Elements, 24:623632, 2000.
  • [17] H.B. Coda, W.S. Venturini, and M.H. Aliabadi. A general 3D BEM/FEM coupling applied to elastodynamic con tinua/frame structures interaction analysis. International Journal for Numerical Methods in Engineering, 46(5):695712, 1999.
  • [18] E.Mesquita, J. Labaki, and L.O.S. Ferreira. An implementation of the longman's integration method on graphics hardware. CMES - Computer Modeling in Engeneering and Sciences, 51(2):143168, 2009.
  • [19] E.Mesquita, J. Labaki, and L.O.S. Ferreira. Constant boundary elements on graphics hardware: a gpu-cpu comple mentary implementation. RBCM, Revista Brasileira de Ciências Mecânicas, 3(4):475482, 2011.
  • [20] W.N. Findley, J.S. Lai, and K. Onaran. Creep and Relaxation of Nonlinear Viscoelastic Materials: With an Intro duction to Linear Viscoelasticity Dover Publications, 1989.
  • [21] G.E. Forsythe. Computer methods for mathematical computations Englewood Cliffs Prentice Hall, New Jersey, 1977.
  • [22] L. Gaul and M. Schanz. Boundary element calculation of transient response of viscoelastic solids based on inverse transformationn. Meccanica, 32:171178, 1997.
  • [23] L. Gaul, M. Schanz, and C. Fiedler. Viscoelastic formulations of BEM in time and frequency domain. Engineering Analysis with Boundary Elements, 10:137141, 1992.
  • [24] G. Gazetas. Analysis of machine foundation vibrations: state of the art. Soil Dynamics and Earthquake Engineering, 2:242, 1983.
  • [25] D. Givoli. Numerical methods for problems in infinite domains. Elsevier, 1992.
  • [26] W.S. Hall and G. Olivetto. Boundary Element Methods for Soil-Structure Interaction Kluwer Academic Publishers, Dordrecht, 1992.
  • [27] B.O. Hardin. The nature of damping in sand. Journal of Soil Mechanics and Foundation Engineering, 91:6397, 1965.
  • [28] Y. Huang, S.L. Crouch, and S.G. Mogilevskaya. A time domain direct boundary integral method for a viscoelastic plane with circular holes and elastic inclusions. Engineering Analysis with Boundary Elements, pages 113, 2006. article in press.
  • [29] D.L. Karabalis and D.E. Beskos. Dynamic response of 3-d rigid surface foundations by time domain boundary element method. Earthquake Engineering and Structural Dynamics, 12:7393, 1984.
  • [30] I.M. Longman. Note on method for computing infinite integrals of oscillatory functions. Cambridge Philosophical Society, 2:764768, 1956.
  • [31] C. Lubich. Convolution quadrature and discretized operational calculus. i. Numerische Mathematik, 52:129145, 1988.
  • [32] C. Lubich. Convolution quadrature and discretized operational calculus. ii. Numerische Mathematik, 52:413425, 1988.
  • [33] N. Makris and J. Zhang. Time-domain viscoelastic analysis of earth structures. Earthquake Engineering and Struc 02 tural Dynamics, 29:745768, 2000.
  • [34] W.J. Mansur. A time-stepping technique to solve wave propagation problems using the boundary element method PhD thesis, University of Southampton, 1983.
  • [35] W.J. Mansur and C.A. Brebbia. Formulation of the boundary element method for transient problems governed by the scalar wave equation. Applied Mathematical Modelling, 6:307311, August 1982.
  • [36] W.J. Mansur and C.A. Brebbia. Numerical implementation of the boundary element method for transient problems goverved by the scalar wave equation. Applied Mathematical Modelling, 6:299306, August 1982.
  • [37] W.J. Mansur and C.A. Brebbia. Transient elastodynamics using a time-stepping technique. In T. Futagami & M. Tanaka C.A.Brebbia, editor, In Proc. 5th Int. Sem. BEM Engng, pages 677698, Berlin, 1983. Springer Verlag.
  • [38] A.D. Mesquita and H.B. Coda. An alternative time integration procedure for boltzmann viscoelasticity: A bem approach. Computers and Structures, 79:14871496, 2001.
  • [39] A.D.Mesquita and H.B. Coda. Alternative viscoelastic procedure for finite elements. Applied Mathematical Modelling, 26:501516, 2002.
  • [40] A.D. Mesquita and H.B. Coda. Boundary integral equation method for general viscoelastic analysis. International Journal for Solids and Structures, 39:26432664, 2002.
  • [41] A.D. Mesquita and H.B. Coda. New methodology for the treatment of two dimensional viscoelastic coupling problems. Computer Methods in Applied Mechanics and Engineering, 192:19111927, 2003.
  • [42] A.D. Mesquita and H.B. Coda. A simple kelvin and boltzmann viscoelastic analysis fo three-dimensional solids by the boundary element method. Engineering Analysis with Boundary Elements, 27:885895, 2003.
  • [43] A.D. Mesquita, H.B. Coda, and W.S. Venturini. An alternative time marching process for viscoelastic analysis by bem and fem. International Journal for Numerical Methods in Engineering, 51:11571173, 2001.
  • [44] E. Mesquita. Zur Dynamischen Wechselwirkung von Fundamenten auf dem Viskoelastischen Halbraum. (On the dynamic interaction of foundations over half spaces) VDI Fortschrittberichte, Reihe 11: Schwingungstechnik, NR , VDI-Verlag Duesseldorf, ISBN 3-18-142011-5, 1989.
  • [45] E. Mesquita, M. Adolph, P.L.A. Barros, and E. Romanini. Transient green and influence functions for plane strain visco-elastic half-spaces. Proc. IABEM SYMPOSIUM 2002, The University of Texas at Austin, Texas, May 28-31 2002.
  • [46] E. Mesquita, M. Adolph, P.L.A. Barros, and E. Romanini. Transient green's functions and distributed load solutions for plane strain, transversely isotropic viscoelastic layers. Latin American Journal of Solids and Structures, 1:75100, 2003.
  • [47] E. Mesquita, M. Adolph, and E. Romanini. Transient response of three dimensional rigid structures interacting with visco-elastic soil profiles. In EM2002, the 15th Engineering Mechanics Division Conference of the ASCE, pages 18, Columbia University, NY, June 2-5 2002. proceedings in CD-ROM.
  • [48] E. Mesquita and R. Pavanello. Numerical strategies for the dynamics of unbounded domains. Computational and Applied Mathematics, 24(1):0126, 2005.
  • [49] S.H. Park. A posteriori evaluation of wave reflection for adaptive analysis of wave propagation in unbounded domains. Computer Methods in Applied Mechanics and Engineering, 193(45-47):49474959, 2009.
  • [50] D.C. Rizos. A rigid surface boundary element for soil-structure interaction analysis in the direct time domain. Computational Mechanics, 26:582591, 2000.
  • [51] D.C. Rizos and D.L. Karabalis. An advanced direct time domain bem formulation for general 3-d elastodynamic problems. Computational Mechancis, 15:249269, 1994.
  • [52] D.C. Rizos and D.L. Karabalis. A time domain bem for 3-d elastodynamic analysis using the b-spline fundamentalsolutions. Computational Mechancis, 22:108115, 1998.
  • [53] D.C. Rizos and Z. Wang. Coupled bem-fem solutions for direct time domain soil-structure interaction analysis. Engineering Analysis with Boundary Elements, 26:877888, 2002.
  • [54] D.C. Rizos and Z.Y.Wang. Applications of staggered bem-fem solutions to soil structure interaction. In In Proceedings of the 13th ASCE Engineering Mechanics Division Specialty Conference, 1999.
  • [55] E. Romanini. Síntese de funções de influência e Green para o tratamento da interaçõo dinâmica solo-estrutura através de equaçõees integrais de contorno. (Synthesis of Green and influence functions to analyse dynamic soil structure interaction by boundary integral equations) PhD thesis, University of Campinas, Unicamp (in Portuguese), 1995.
  • [56] M. Schanz. Wave Propagation in Viscoelastic and Poroelastic Continua - A Boundary Element Approach Springer Verlag, Berlin, 2001.
  • [57] M. Schanz and H. Antes. Application of 'operational quadrature methods' in time domain boundary element methods. Meccanica, 32:179186, 1997.
  • [58] M. Schanz and H. Antes. A new visco- and elastodynamic time domain boundary element formulation. Computational Mechanics, 20:452459, 1997.
  • [59] M. Schanz, H. Antes, and T. Räberg. Convolution quadrature bem for quasi-static visco- and poroelastic continua. Computers and Structures, 83:673684, 2005.
  • [60] B. Senzale, P.W. Partridge, and G.J. Creus. General boundary elements solution for ageing viscoelastic structures. International Journal for Numerical Methods in Engineering, 50:14551468, 2001.
  • [61] A. Sommerfeld. Partial Differential Equations Academic Press, New York, 1949.
  • [62] O. von Estorff. Coupling of bem and fem in the time domain - some remarks on ist applicability and efficiency. Computers and Structures, 44:325337, 1992.
  • [63] O. von Estorff and H. Antes. On fem-bem coupling for fluid-structure interaction analysis in the time domain. International Journal for Numerical Methods in Engineering, 31:11511168, 1991.
  • [64] Y. Wang and R.K.N.D. Rajapakse. Transient fundamental solutions for a transversely isotropic elastic half-space. Proc. Royal Soc. London A, 442:505531, 1993.
  • *
    Author email:
  • Publication Dates

    • Publication in this collection
      21 Jan 2013
    • Date of issue
      Aug 2012

    History

    • Received
      16 Apr 2012
    • Accepted
      16 May 2012
    Individual owner www.lajss.org - São Paulo - SP - Brazil
    E-mail: lajsssecretary@gmsie.usp.br