Three-dimensional transient half-space dynamics using the dual reciprocity boundary element method

ARTICLE IN PRESS

Engineering Analysis with Boundary Elements 32 (2008) 597–618 www.elsevier.com/locate/enganabound

Three-dimensional transient half-space dynamics using the dual reciprocity boundary element method Andrej Tosecky´a, Yvona Kolekova´c, Gu¨nther Schmida,, Valerij Kalinchukb,1 a

Faculty of Civil Engineering, Ruhr University Bochum, Universita¨tsstrasse 150, 44801 Bochum, Germany Institute of Mechanics and Applied Mathematics, University of Rostov, Stachky ave 200/1, 344090 Rostov on Don, Russian Federation c Faculty of Civil Engineering, Slovak Technical University, Radlinske´ho 11, 81368 Bratislava, Slovakia

b

Received 5 December 2005; accepted 11 October 2007 Available online 18 December 2007

Abstract The dual reciprocity boundary element method (DRBEM) is studied thoroughly in the present paper. To the best knowledge of the authors, the DRBEM has never been applied to 3D half-space dynamics previously. In the present paper, the mathematical derivation of the method is presented, with stress placed on peculiarities of the method when applied to transient half-space dynamics. It has been found that semi-infinite domains (half-space) are more difficult to simulate, since the truncation of the discretization of the half-space surface by boundary elements results in excessive spurious wave reflection on the border of the discretized region. Mathematical derivation of the method is followed by its numerical implementation. Wave propagation due to various kinds of loading in the time-domain is studied afterward, aimed at the validation of the method. The main advantage of the DRBEM over its counterparts, used to model infinite and semiinfinite domains (such as classical BEM formulation, integral transformations, thin layer method), is that it produces time- and frequencyindependent matrices (mass and stiffness matrix), by preserving a boundary-only discretization (no internal nodes necessary). This makes the method very attractive, since this feature is very close to the common engineering understanding and the final equation of motion has a similar form like the one known from the finite element method. Moreover, the formulation allows for seamless incorporation of nonhomogeneous initial conditions, i.e. non-zero initial displacements and velocities and surface tractions can be prescribed. r 2007 Elsevier Ltd. All rights reserved. Keywords: Dual reciprocity BEM; Semi-infinite domains; Half-space dynamics; Wave propagation; Moving loads

1. Introduction The most frequently used numerical method in the common civil engineering practice is the finite element method (FEM). However, there are cases where problems are encountered where the use of the classical FEM is cumbersome and inadequate. One such case is the modeling of the building subgrade, which, by its nature is a semi-infinite medium. Unboundedness of the subgrade disables any seamless application of the FEM which always requires domain discretization. Of course, only a limited portion of the semi-infinite medium can be discretized and Corresponding author. Tel.: +49 173 5443146.

E-mail addresses: [email protected], [email protected] (G. Schmid). 1 Present address: Southern Scientific Centre of Russian Academia of Science, Chekhov st. 41, 344006 Rostov on Don, Russia. 0955-7997/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2007.10.004

the region of interest must be truncated at some distance. Truncation of the region leads to serious errors in dynamic analysis, where waves reflected from the artificially created boundary swamp the results to such an extent that these become unacceptable. Some solutions to this problem can be found in the work of Wolf and Song [1]. In subgrade modeling other methods are firmly established and regarded as reliable and robust. Three groups of methods are to be mentioned: (1) Classical boundary element method (BEM). Its application to transient half-space dynamics requires the knowledge of the fundamental solution of the wave equation. This solution is known, so that a straightforward derivation of a pure boundary-only formulation for a general linear elastodynamic problem is possible. However, there are two things to frown upon. Firstly,

ARTICLE IN PRESS 598

A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

since the elastodynamic fundamental solution is itself time-dependent the formulation leads to a convolution integral. This feature complicates the numerical implementation of the classical BEM seriously, for it demands a new set of system matrices to be assembled for each single time-step of the analysis. CPU time and memory storage requirements are prone to become prohibitive. Secondly, the fundamental solution of elastodynamics assumes homogeneous initial conditions, what in turn means, that only homogeneous initial conditions can be prescribed over the domain of interest if no domain discretization is employed. An interested reader is referred to [2] to find more on this subject. (2) Methods based on integral transformations (IT) follow the idea of transforming the governing differential equation together with the given boundary conditions into an another space (Fourier, Laplace, wavenumber space). It is usually much easier to find the sought solution in the transformed space. The major difficulty resides in the back-transformation of the solution into the original space. An exact back-transformation is only rarely possible, so numerical techniques must come to aid. Numerical effort at this back-transformation is usually an enormously time-consuming operation. This method is reviewed in detail in [3], together with suggestions on how to accelerate its applications. (3) The thin layer method (TLM) enables the effective modeling of a stratified half-space in the frequencydomain. The method combines spatial discretization by FEM in vertical direction together with an analytical formulation of the wave propagation into infinity in horizontal direction realized by solving the problem in the wavenumber–frequency domain. The theoretical background of the method can be found in [4–6]. Understanding the methods listed above and their related techniques requires a considerable knowledge of mathematics and mechanics. Their formulation differs significantly from the formulation of the FEM, with which most design engineers are familiar. It would be therefore desirable to come up with a method enabling the modeling of a semi-infinite media without having to discretize its volume and, at the same time, resulting in an equation of motion in a similar form as the one known from the FEM, i.e.: M u€ ðtÞ þ C u_ ðtÞ þ KuðtÞ ¼ f ðtÞ,

(1)

where M, C and K are mass, damping and stiffness matrix, respectively. Vectors u€ ðtÞ, u_ ðtÞ, uðtÞ and f ðtÞ contain timedependent accelerations, velocities, displacements and external forces, respectively. One method possibly fulfilling the demands stated above is the dual reciprocity boundary element method (DRBEM) [7]. Its formulation yields an equation of motion very similar to (1) where the definition of internal nodes is not necessary for a sufficiently accurate solution, though if done so, an improved quality of results can be expected.

To the best knowledge of the authors, the majority of applications of the DRBEM to transient dynamics is limited to closed domains. Applications to transient solid dynamics of infinite and semi-infinite domains could not be found in the literature at all. There are some texts on DRBEM applications to elastodynamic problems in bounded domains [8]. The application of the DRBEM for exterior and interior domains for a Stokes flow is discussed in [9]. As will be shown later in the text, if the domain under consideration is an unbounded domain, DRBEM formulation has to fulfill the regularity conditions at infinity. It requires special treatment which is different from the DRBEM applicable to closed domains. Let us commence now with a concise mathematical derivation of the DRBEM for a general elastodynamic problem. A detailed treatise on the specialties of infinite and semi-infinite domains will follow. 2. Derivation of DRBEM formulation for a general elastodynamic problem Starting point of the derivation of the boundary integral equation for elastodynamics using the dual reciprocity approach is the equilibrium equation of a homogeneous linear elastic solid. This equation in term of stresses is written as sjk;k ðx; tÞ þ rðbj ðx; tÞ  u€ j ðx; tÞÞ ¼ 0;

j; k ¼ 1; 2; 3,

(2)

where sjk;k ðx; tÞ is a stress tensor component, bj ðx; tÞ is a component of the body force vector (force per unit mass) and r denotes material density. Variables x and t denote space- and time-dependency of the functions, respectively. Application of the method of weighted residuals to Eq. (2) leads to the weak form of the equation of motion: Z ðsjk;k ðx; tÞ þ rðbj ðx; tÞ  u€ j ðx; tÞÞÞu^ j ðxÞ dO ¼ 0 (3) O

or term by term as Z Z sjk;k ðx; tÞu^ j ðxÞ dO þ rbj ðx; tÞu^ j ðxÞ dO O O Z ru€ j ðx; tÞu^ j ðxÞ dO ¼ 0 

ð4Þ

O

with u^ j ðxÞ being a weighting function (note that the chosen weighting function is time-independent). At this moment it is worth noting that the classical BEM formulation uses the ^ tÞ as weighting fundamental solution of elastodynamics uðx; function. This assumption leads to the elimination of the domain integrals involved. On the other hand, the dual reciprocity approach uses the time-independent fundamental solution of elastostatics as the weighting function. An instantaneous merit of this assumption is the disappearance of the convolution product in the formulation and its replacement by a dot product. The price to be paid for this simplification is that one volume integral persists in the formulation. Its straightforward reduction to a boundary integral is no longer possible as it is by the classical BEM

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

formulation. Integrating Eq. (4) by parts, applying Gauss’s theorem, Betti’s reciprocal work theorem and for the sake of simplicity assuming zero prescribed body forces over the domain O, i.e. bj ðx; tÞ ¼ 0,2 we obtain Z Z cij ðnÞuj ðn; tÞ þ t^ij ðxÞuj ðx; tÞ dG ¼ tj ðx; tÞu^ ij ðxÞ dG G G Z ru€ j ðx; tÞu^ ij ðxÞ dO,  O

599

leads to the following result: Z ru€ j ðx; tÞu^ ij ðxÞ dO O

r

N X

 Z n a€ ðtÞ cij ðnÞcjm ðnÞ  t^ij ðxÞcnjm ðxÞ dG n

G

n¼1

Z þ

G

 u^ ij ðxÞZnjm ðxÞ dG

ð9Þ

ð5Þ where n denotes the so-called source point, cij are the boundary coefficients, t^ij ðxÞ the surface tractions corresponding to the fundamental state and tj ðx; tÞ are the surface tractions of the actual state. The volume integral persisting on the right-hand side of Eq. (5) contains inertial forces over the domain O. Its transformation to a boundary integral is realized by the dual reciprocity approach.

which means that the volume integral persisting in Eq. (5) is expressible in terms of boundary integrals only. Finally, the complete boundary integral equation for elastodynamics as it results from application of the DRBEM reads: Z cij ðnÞuj ðn; tÞ þ t^ij ðxÞuj ðx; tÞ dG G

Z

u^ ij ðxÞtj ðx; tÞ dG þ r

¼ G

3. Transformation of the domain integral

Z

With a certain degree of approximation, the time- and space-dependent inertial forces ru€ j ðx; tÞ in the volume integral in Eq. (5) can be expressed by a series of products of merely time-dependent factors aðtÞ and merely spacedependent functions f ðxÞ as follows: ru€ j ðx; tÞ  r

N X

a€ n ðtÞf nj ðxÞ,

(6)

n¼1

where N is the number of terms used in the expansion. Consider now rf nj ðxÞ to be a continuous body force field defined over the domain O, acting in j-direction and centred on the nth node. As a matter of fact, these body forces induce displacements over the domain O. Denote these displacements by cnjm ðxÞ and the corresponding surface tractions by Znjm ðxÞ. Displacements cnjm ðxÞ are obtained as particular solution to the Navier’s equation3 Gcnj;mm ðxÞ þ

G cn ðxÞ þ rf nj ðxÞ ¼ 0; 1  2n m;jm

j; m ¼ 1; 2; 3, (7)

where G is shear modulus and n is Poisson’s ratio. Displacements are linked with the surface tractions via material law Zij ¼

2Gn c nj þ Gðcij;k nk þ cik;j nk Þ; 1  2n ik;k

i; j; k ¼ 1; 2; 3, (8)

with nj and nk being components of the unit outward normal. Taking into account space-independence of the aðtÞ factors, a second application of the reciprocity theorem 2 Assumption of zero body forces does not deteriorate generality of the approach, because contingent body forces can be treated in the same € tÞ. manner as inertial forces represented by the term ruðx; 3 Governing differential equation of elastostatics.

þ G

N X

 a€ n ðtÞ cij ðnÞcnjm ðnÞ

n¼1

t^ij ðxÞcnjm ðxÞ dG 

Z

G

u^ ij ðxÞZnjm ðxÞ dG

 .

ð10Þ

Discretizing the domain’s boundary G by N e boundary elements with a total of N boundary nodes, a linear, generally non-symmetric, system of algebraic equations is obtained (for details see [10]): M u€ ðtÞ þ HuðtÞ ¼ GtðtÞ,

(11)

where matrices G and H store integrated values of the fundamental solution for displacements u^ ij and surface tractions t^ij , respectively. M is a mass matrix defined by M ¼ rðGg  HwÞF 1 ,

(12)

where matrices w and g contain values of the particular solution for displacements cjm and surface tractions Zjm , respectively, due to the body force field described by the chosen inertial force approximating function f ðxÞ. Matrix F stores values of this approximating function f ðxÞ. In the outlined formulation its inverse F 1 appears. This inversion brings some additional numerical effort. However, an equivalent formulation4 without matrix inversion is also possible [11]. It is important to underline that all participating matrices are time-and frequency-independent. For a more detailed derivation of the DRBEM formulation stated only briefly in the previous text, an interested reader is referred to literature [7,12]. 4 If exclusively Neumann boundary conditions are imposed on the boundary G, the formulation without matrix inversion yields exactly identical results as the formulation with matrix inversion. If also Dirichlet or Robin boundary conditions are present, the formulation without matrix inversion is still possible, though with some additional computational effort.

ARTICLE IN PRESS 600

A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

4. Approximation of inertial forces The most vulnerable spot of the appealing DRBEM approach is its capability of a satisfactorily good approximation of the inertial forces over the domain O by the adopted expansion (Eq. (6)). In principle, improvement of the approximation can be achieved in two different ways. One way is to increase the number of terms in the expansion, i.e. to define more boundary nodes N, or, what is an issue frequently discussed in the DRBEM literature, to define also auxiliary internal nodes in order to capture the distribution of the inertial forces also in the interior of the domain. Bearing in mind, that in a general elastodynamic problem inertial forces do occur over the entire domain and not only on its boundary, the definition of internal nodes seems to be absolutely inevitable. This conclusion sounds quite discouraging especially in the context of DRBEM application to transient dynamics of infinite and semi-infinite domains. Covering the domain O by a network of internal nodes is possible in closed domains, but not in infinite and semi-infinite domains. And moreover, the necessity of internal nodes seriously degrades the attractiveness of the BEM, since the boundaryonly formulation gets lost. However, in the opinion of the authors, there is a fundamental difference between closed and infinite (also semi-infinite) domains concerning the necessity of internal nodes. In a closed domain, a source located on the domain’s boundary induces inertial forces over the whole domain and multiple wave reflections on the domain’s boundary give rise to possibly very complicated distribution patterns of inertial forces. To render these patterns by the adopted approximation (Eq. (6)), a network of internal nodes is essential, for otherwise we have absolutely no idea about what’s happening inside of the domain. In infinite and semi-infinite domains, on the other hand, provided that the source is located on the domain’s boundary, the induced waves propagate into infinity so there are no reflections. If we are interested only in the situation on the boundary, we do not need to describe waves in the domain’s interior, since we know that these are radiated into infinity, never to return and never influencing the events on the boundary once they have left it. This holds true as long as the medium remains homogeneous, without any inclusions, irregularities or inhomogeneities. Otherwise, wave reflections occur on the interfaces between two distinctive media. A second way to improve the accuracy of the approximation of the inertial forces is the proper choice of the approximating function f ðxÞ. Various functions are proposed in the literature. The most frequently used are the so-called globally supported radial basis functions based on the distance r between the so-called source point and observation point. In the further context of the present work use will be made of these functions only, therefore designation f ðrÞ will be used henceforward to indicate their dependency on the length of the position vector r.

In closed domains discretized by boundary as well as by internal nodes, the range of choice for the approximating function f ðrÞ is substantially wider than for infinite and semi-infinite domains discretized by boundary nodes only. For infinite and semi-infinite domains, the approximating function must exhibit a decay behavior of some kind, i.e. the value of the function must tend to zero as r ! 1, to reflect the effect of the ‘‘radiation damping’’ and to satisfy the regularity conditions. An example of an approximating function f ðrÞ for which the corresponding particular solution to Navier’s equation is known is f ðrÞ ¼ 1 þ r. It is a fairly simple linear function, and what is more, not possessing the demanded decay behavior. Solutions to some other functions, also not possessing the demanded decay behavior, are given in [13]. Therefore these functions are of no use in DRBEM applications to half-space dynamics. Samples of possible functions are 1 f ðrÞ ¼ ; r

f ðrÞ ¼

1 ; r2

f ðrÞ ¼

1 ; 1þr

f ðrÞ ¼ er .

(13)

The most difficult task is finding the corresponding particular solution to the Navier’s equation for body forces described by functions with a decay behavior (like the ones shown above) and acting in only one direction of the Cartesian coordinate system. Possible ways leading to the derivation of a new particular solution are presented and shortly discussed in the next section. 5. Particular solutions to the Navier’s equation Navier’s equation (7) is a system of three coupled partial differential equations of second order. Its solution yields three components of the displacement vector due to the prescribed boundary conditions and the acting body forces. There is no general method of solving this equation guaranteeing a closed-form solution. However, there are several methods possibly leading to the sought solution. These methods are listed below: (1) (2) (3) (4) (5)

Principle of superposition. Solution utilizing body force field potentials. Fourier transformation of the Navier’s equation. Papkovich–Neuber potentials. Inverse procedure.

The first four options have been tested thoroughly in the frame of this work. It turns out unfortunately that functions like (13) do exhibit unpleasant characters in conjunction with the Navier’s equation. The solution leads either to volume integrals or to a solution of some modified partial differential equation, either of these being unsolvable just the same. Let us consider now the fifth option, the inverse procedure. A very elegant and straightforward way possibly leading to a new particular solution is to assume a certain displacement field rather than body force field

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

601

Table 1 Particular solutions to the Navier’s equation discovered in the context of this work f ðrÞ

Particular solution for displacements uij and surface tractions tij

r3

1 5 uij ¼ 96ð1nÞG ðð3  32 10nÞdij  r;i r;j Þr    qr qr r4 dij þ ni r;j ð1  dij Þ tij ¼ 48ð1nÞ nj r;i ðð1  8nÞð1  dij Þ  6dij Þ þ 3r;i r;j þ ð7 þ 8nÞ qn qn

1 r6

1 75 1 uij ¼ 6ð1nÞGr 4 ðð100  2 nÞdij  r;i r;j Þ     qr qr 1 þ ð2 þ nÞ dij þ ð1  dij Þni r;j  r;i nj ð3dij þ ð1  dij Þð1 þ nÞÞ tij ¼ 3ð1nÞr 5 6r;i r;j qn qn 1 1 uij ¼ 2ð1nÞGr 2 ðð2  nÞdij þ r;i r;j Þ    qr qr 1 þ n dij þ ð1  dij Þr;j ni tij ¼ ð1nÞr 3 nj r;i ðdij þ ð1  dij Þð1  nÞÞ  4r;i r;j qn qn

1 r4

1 r

1 uij ¼ 16ð1nÞG ðð7  8nÞdij  r;i r;j Þr    qr qr 1 tij ¼ 8ð1nÞ nj r;i ð2dij þ ð1  dij Þð4n  1ÞÞ þ r;i r;j þ ð3  4nÞ dij þ ð1  dij Þr;j ni qn qn

p1ffi r

2 uij ¼ 45ð1nÞG ðð16  6nÞdij  r;i r;j Þr3=2  3   pffi qr qr 2 r tij ¼ 45ð1nÞ nj r;i ð5dij þ ð1  dij Þð9n  2ÞÞ þ r;i r;j þ ð7  9nÞ dij þ ð1  dij Þr;j ni qn qn

described by some special function f ðrÞ. Making use of the Navier’s equation again, if the displacement field is known, we can readily determine the body force field which has caused these displacements. Note that in doing so, nothing else but mere partial derivations of known functions are to by carried out. In this manner, assuming a reasonable displacement field,5 we may end up with an approximating function f ðrÞ as useful for our purposes as any other would have been. However, easier said than done. Following this approach, we will pretty soon find out that it is not that easy to come to a desirable result. In principle, the difficulty resides in two points. Firstly, even if the assumed displacement field might have been given in the form of a relatively simple function, after having performed all the partial derivations, we end up with a rather complicated expression for the approximating function, containing not only the position vector r as a variable but also Poisson’s ratio and some other undesired variables. Such functions are of not much use for us, for they do exhibit an erratic behavior. Secondly, recall that the body forces must act in only one direction of the Cartesian system, which in conjunction with the Navier’s equation (system of three coupled partial differential equations of second order) practically means that only one of the equations comprised in the Navier’s equation exhibits a non-zero right-hand side, whereas right-hand sides of the other two must vanish. The probability that our assumption for the

5 Keep in mind, that the approximating function f ðrÞ must exhibit a decay behavior of some kind. Therefore we are not utterly free in the choice of the displacement field. It must also have certain properties, reflecting the demanded properties of the approximating function describing the acting body force field.

displacement field will satisfy this demand drops almost to zero. Despite of all the involved difficulties, the perseverance finally paid off in finding of coveted particular solutions. Some very promising particular solutions found in frame of the present work are listed in Table 1. Apart from the first one listed in the Table 1, all the approximating functions f ðrÞ to which the particular solution has been found, do possess the demanded decay behavior, each of a different rate. From a point of view of numerical implementation, the last two are especially promising, because for these the corresponding particular solutions do not exhibit singular behavior, even if the approximating functions themselves do (f ðr ¼ 0Þ ! 1). Looking at the form of the particular solution of the displacements it is predictable that for certain functions f ðrÞ a general form of the solution can be written (no rigorous mathematical proof is attached): f ðrÞ ¼

uij ¼

1 , r2n

1 ððC 2 þ C 3 nÞdij  r;i r;j Þ, C 1 ð1  nÞGr2n2 n ¼ 2; 3; 4 . . . ,

ð14Þ

where C 1 ; C 2 and C 3 are constants to be determined iteratively. For a trouble-free application of the method, it would be most desirable to have particular solutions available, for which neither the approximating function f ðrÞ, nor the corresponding particular solution do exhibit any singularities. Unfortunately, no such particular solutions have been found yet. Numerical treatment of the singularities

ARTICLE IN PRESS 602

A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

z p(t)

Zoomed in

y x

z

y p(t) 3m

32

m

x 3m

32m Discretized half-space surface

Half-space

p(kN/m2) Loading time-history 1 0

t(s)

Fig. 1. Boundary element mesh of a half-space surface.

encountered by the approximating functions listed in Table 1 is addressed in Section 8. 6. Numerical example—wave amplitude attenuation Upon knowledge of applicable particular solutions, we can employ the DRBEM to study a simple case of wave propagation in an elastic half-space and assess its performance and especially the performance of various inertial forces approximating functions. The half-space model used in this study is portrayed in Fig. 1. The half-space surface measuring 32 m  32 m is covered by a regular mesh of linear quadrilateral boundary elements with size 1 m  1 m. The half-space surface is subjected to Heaviside loading with intensity pðtÞ ¼ 1 kN=m2 applied normally in the middle of the discretized region as a surface traction on an area of 3 m  3 m ¼ 9m2 . The problem is solved in the time-domain by Newmark’s b method of constant accelerations. A certain amount of numerical damping is introduced into the solution by setting the Newmark’s parameter g ¼ 0:7. The time-step length used in the analysis is Dt ¼ 0:004 s. Matrices G and H are calculated using the static Green’s functions for the considered homogeneous elastic half-space. Material properties of the half-space are summarized in Table 2. The given loading induces all three types of waves possible in a homogeneous elastic half-space: p-, s- and R-wave. Under the given loading conditions, the only well observable wave on the surface is the Rayleigh wave. The surface amplitude of a Rayleigh wave excited by an impulsive point source decays proportionally to 1=r, where r is distance from the source [14]. Decay rates obtained by four different inertial force approximating functions f ðrÞ are presented in Fig. 2. It is clearly seen that the linear approximating function f ðrÞ ¼ 1 þ r can be immediately marked as inappropriate and excluded from any further investigation. Astonishing is the comparatively better performance of the function

Table 2 Material properties of the used half-space model Shear modulus

G ¼ 17 MPa

Poisson’s ratio

n ¼ 0:4

Density

r ¼ 1700 kg=m3

Rayleigh wave velocity Shear wave velocity Dilatational wave velocity

cR ¼ 94:2 m=s cs ¼ 100:0 m=s cp ¼ 244:9 m=s

f ðrÞ ¼ r3 . Even if this function does not possess any decay behavior, on the contrary, its value grows much more rapidly than f ðrÞ ¼ 1 þ r, the amplitude decay obtained by this function complies much better with the true decay rate. A discussion on the superiority of the approximation functions f ðrÞ ¼ r3 is given in [15]. Note the very good p agreement of the amplitude decay ffiffi obtained by f ðrÞ ¼ 1= r and the excellent agreement obtained by f ðrÞ ¼ 1=r. Is the excellent agreement just accidental? One notices immediately that the correct decay rate and the approximating function are identical. One can ask if the function f ðrÞ ¼ 1=r performs so well only under the circumstance that the correct decay rate is identical with it or can it be used as a general tool, also in cases where the true decay rate is different? It is hard to give a final answer now. But to sustain the idea that the function f ðrÞ ¼ 1=r can be used as a general tool, another test has been conducted. The same half-space model is used again, but the Heaviside loading is replaced now by a harmonic loading pðtÞ ¼ p0 sin ot applied again as a normal surface traction with amplitude p0 ¼ 1 kN=m2 over an area of 3 m  3 m ¼ 9 m2 . The load’s excitation frequency is f ¼ 10 Hz. As illustrated in Fig. 3, the acting point of the loading has been shifted closer to the border of the discretized region in order to enable the observation of the propagating waves along a longer path (29 m). The surface amplitude of Rayleigh waves excitedpby ffiffi a harmonic point source decays in accordance with 1= r.

ARTICLE IN PRESS

Vertical wave crest amplitude (m)

A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

603

2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

f(r)=1+r f(r)=r3 f(r)=1/sqrt(r) f(r)=1/r Decay by 1/r

0

1

2

3

4

5

6 7 8 9 10 Distance from the source (m)

11

12

13

14

15

Fig. 2. Amplitude decay rate for various approximating functions f ðrÞ.

z y

Zoomed in

p(t) x y

p(t) 3m

z 32

m

x 3m

32m Discretized half-space surface p(kN/m2) Loading time-history

Half-space

1 t(s)

-1

Fig. 3. Half-space surface subjected to harmonically oscillating loading.

Rayleigh wave crest amplitude ( m)

1.95 Decay according to 1/sqrt(r) f(r) = 1/r f(r) = 1/sqrt(r)

1.90 1.85 1.80 1.75 1.70 0 1

2 3

4 5 6

7 8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 Distance from the source (m)

Fig. 4. Amplitude decay of Rayleigh waves caused by a harmonically oscillating source—comparison of performance of two different inertial force approximating functions f ðrÞ.

The performance of two inertial force approximating pffiffi functions has been tested, namely of f ðrÞ ¼ 1= r and 1=r. The results obtained are shown in Fig. 4. It appears that the best results are obtained by the function f ðrÞ ¼ 1=r, which is again in a very good agreement with the true decay rate. It is worth saying a few words about the inertial force approximating functions. It has been clearly demonstrated that each approximating function yields different results. Deviation of the calculated results from the true physical

character of the approximated quantity (inertial forces in our case) becomes essential only if the approximating function and the approximated quantity do exhibit fundamentally different behavior. This was the case of the approximating function f ðrÞ ¼ 1 þ r, which has been shown as completely incapable of being used for simulation of wave propagation processes in a semi-infinite media. The observations lead to the conclusion that if no domain discretization is employed in order to reconstruct the distribution of the approximated quantity, the approximating

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

604

function used must intrinsically incorporate the natural physical conditions governing the given problem. 7. Wave reflection The superiority of the BEM over FEM in modeling infinite and semi-infinite media resides in two aspects. The first one is the reduction of the problem’s dimensionality by one, which means an immense advantage for modeling of infinite and semi-infinite domains. The second one is the implicit inclusion of the regularity condition (representation of outgoing waves). In BEM one has not to introduce ‘‘artificial boundaries’’ to truncate the infinite extension of the domain which lead to incorrect wave reflections in the classical FEM application in half-space dynamics. One could say that the DRBEM is in between. It has something in common with both, with classical BEM as well as with FEM. Since no literature could be found on the performance of DRBEM with respect to wave reflections, we address this issue here. Only a limited portion of the half-space surface can be covered by boundary elements and the discretization must be truncated at some distance. In this section we shall find out whether the discretization truncation allows undisturbed wave transmission into the outer space or causes some spurious wave reflections. The half-space model portrayed in Fig. 1 is used again. The impulsive character of the applied Heaviside loading induces waves propagating into infinity in circular patterns. The loading acts as a static load afterward, so no further waves are induced. After the initially induced waves have left the discretized half-space surface, only the static displacement field remains present. This static displacement field is portrayed in Fig. 5. Any deviation of the dynamic solution from the static displacement field hints at existence of wave reflections in the dynamic solution. The dynamic solution at 10 consecutive time-instants is captured in a series of snapshots shown in Fig. 6. Looking at the series of snapshots, it is immediately clear that the wave reflection problem exists in the DRBEM solution after all. Even worse, not only mild negligible reflections are observable, but rather wide spread multiple reflections causing abnormal magnification of the displacements. Energy carried by the surface waves seems to be trapped in the discretized portion of the half-space surface

Fig. 5. Static displacement field.

and there are apparently no means how it could be radiated into the outer space. To the best knowledge of the authors, the problem of wave reflections is not addressed in the BEM literature. For the sake of brevity, we waive any lengthy treatise on all the various concepts developed and tested in frame of the present work and confine the explanation to the concept which has proved to be successful at last. As already said, under the given loading conditions (Heaviside loading) the dynamic solution must approach the static one. Let us compare the equation system to be solved by both, the static and the dynamic analysis: static analysis:

Hu ¼ Gt

dynamic analysis: ðH þ a0 MÞu1 ¼ Gt1 þ Mða0 u0 þ a2 u_ 0 þ a3 u€ 0 Þ. ð15Þ The above equation of motion for dynamic analysis arises when Newmark’s b method is employed. Upper indices 1 and 0 refer to the current and to the previous time-step, respectively. The constants a0 , a2 and a3 result from the chosen Newmark parameters b, g and the time-step length Dt. As already said, after the last wave has left the discretized region, accelerations and velocities at all nodes must be zero. Hence, the equation to be solved by dynamic analysis reduces then to ðH þ a0 MÞu1 ¼ Gt1 þ Mða0 u0 Þ.

(16) 1

Provided that the external loading t does not vary with time anymore (e.g. Heaviside load), time-variation in the system is caused only by the difference between u1 and u0 . Because velocities as well as accelerations are zero, the difference between u1 and u0 should also be zero, i.e. u1 ¼ u0 . Consequently, the equation system for dynamic analysis reduces further to Hu1 ¼ Gt1 ,

(17)

so we finally obtain a system of equations identical with the static one. How to force the velocities and accelerations tend to zero, when the wave they are associated with has left the discretized region? In the original equation system for dynamic analysis (15) all loading terms are grouped on the right-hand side. The first loading term (Gt1 ) is the external loading applied at the current time-step 1. The second loading term ðMða0 u0 þ a2 u_ 0 þ a3 u€ 0 ÞÞ expresses the state of the system at the previous time-step 0 and physically represents inertial forces over the system. In this term, velocities u_ 0 and accelerations u€ 0 from the previous timestep 0 show up. If these velocities and accelerations are associated with a wave which has already left the region, they must be zero, for this wave propagating away from the discretized region does not influence this region anymore. Since the system does not take this intrinsically into account, we must artificially alter the loading term and set all velocities and accelerations at the nodes on the border of the discretized region to zero. Note, we alter the

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

605

Fig. 6. Surface displacements due to a Heaviside loading applied in the middle of the discretized region. The solution is deteriorated by spurious wave reflections.

velocities and accelerations from the previous time-step only in the loading term, not the velocities and accelerations themselves. Generally these are non-zero and we calculate them quite normally in the current time-step. We do not include these calculated quantities in the loading term for the next time-step and set them to zero instead. Dividing nodes of our model into N regular nodes i.e. those not lying on the border of the discretized region (n 2 1; . . . ; N) and into M border nodes (see Fig. 7) we can write the velocity and acceleration vector in the loading term as u_ 0 ¼ fu_ 01x ; u_ 01y ; u_ 01z ; . . . ; u_ 0Nx ; u_ 0Ny ; u_ 0Nz ; 0; 0; 0; . . . ; 0; 0; 0gT , u€ 0 ¼ fu€ 01x ; u€ 01y ; u€ 01z ; . . . ; u€ 0Nx ; u€ 0Ny ; u€ 0Nz ; 0; 0; 0; . . . ; 0; 0; 0gT . ð18Þ

border nodes regular nodes

Fig. 7. Division of the discretized region into regular and border nodes.

The last thing to deal with are the displacements. Shall all displacements at the border nodes be set to zero, just like the velocities and accelerations were? Presumably not. If we did so, we would actually constrain the system on the border of the discretized region which would utterly violate the kinematics of the system and would lead to

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

606

wave reflections again. But since the displacements on the border of the discretized region are not calculated correctly we have to adjust them somehow. Under the condition that all incoming waves originate from somewhere inside of the discretized region and that the properties (direction and propagation velocity) of these waves at regular nodes adjacent to the border nodes are known, we can predict (or forecast in other words) the displacements at the border nodes. Then we can prescribe the forecasted values in place of the calculated values in the loading term of the equation system for dynamic analysis. This concept is illustrated in Fig. 8. The displacement vector can be partitioned and written as 0t 0t T u0 ¼ u01x ; u01y ; u01z ; . . . ; u0Nx ; u0Ny ; u0Nz ; u0t x ; uy ; uz ; . . . g . |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} regularly calculated values

forecasted values

(19)

r d2 d1

n irectio

ation d

propag

border node

regular node

X border node

Fig. 8. Forecasting of displacements on the border of the discretized region.

S1

S2

S2 S2

S4

S3 S3

S6

S6

S1

S2

S4

S5

S5 S5

S5

Fig. 9. Simplified concept of displacement forecasting.

The time-delay t is the time-duration necessary for a wave to advance from a regular node to a point on the border of the discretized region. An exact position where the wave strikes the border can be determined (point X in Fig. 8) and the corresponding displacement contributions prescribed at the next two closest border nodes (the contributions are determined by simple linear extrapolation). In this manner, we have corrected the boundary conditions on the border of the discretized region. The system ‘‘does not know’’ that it has been truncated here. The results presented below have been obtained by using a slightly simplified method of boundary condition correction. The simplicity resides in not calculating the exact position where a wave detected at a regular node strikes the border and then extrapolating the forecasted displacement into its nodal contributions at the border nodes, but rather by assuming unchangeable regular node–border node pairs. The displacement value calculated at a regular node is then prescribed at its corresponding border node, with the corresponding time-delay t. Fig. 9 clarifies the concept. S refers to a regular node and S0 refers to a border node. The inaccuracy introduced into the concept by the simplified version resides in incorrect displacement forecasting when a wave arriving at node S inclines from the straight direction SS 0 . An inclined wave does not hit exactly border node S 0 belonging to a regular node S, but rather a different location on the border of the discretized region. This is not taken into account by the simplified concept and consequently some inaccuracy deteriorates the efficiency of the wave reflection reduction. The inaccuracy is the bigger, the more the propagation direction of the wave at a node S is inclined from the direction defined by the straight line SS 0 . This problem is illustrated in Fig. 10. To avoid the fourth case, i.e. severe inaccuracies due to the considerably inclined waves, the loading should not be applied close to the border of the discretized region. The displacement values calculated at Si nodes are prescribed at the corresponding border nodes S 0i belonging to the given Si node, with the corresponding time-delay t

S

S

S

S

S

S

S

S

Fig. 10. Inaccuracy introduced by the simplified concept. (a) No inaccuracy; the wave hits exactly the node S0 . (b) Small inaccuracy; slightly inclined wave hits the border near the node S0 . (c) Tolerable inaccuracy; noticeably inclined wave hits at least the adjacent element. (d) Severe inaccuracy; the inclined wave hits the border at some remote node.

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

607

Fig. 11. Wave propagation due to a Heaviside loading applied normally to the half-space surface. Solving procedure with wave reflection removal algorithm.

depending on the distance jS i S 0i j and the wave propagation velocity. The previous example of a half-space surface subjected to a Heaviside loading has been calculated anew using the illustrated concept. A series of snapshots at 10 timeinstances is presented in Fig. 11. Solving the given problem by the enhanced computer algorithm leads to a substantial reduction of the spurious wave reflections. A more exact assessment of the effectiveness of the developed strategy offers Fig. 12 where the time-history of the vertical displacement in the middle of the loaded area is given. The dynamic solution must approach the static one. Displacements in the dynamic solution must not vary with time after the last wave has left the discretized region. Looking at Fig. 12 we observe that the spurious wave reflections still persist and some small part of the energy is

still reflected back on the border of the discretized region. A minute discrepancy between the static and the dynamic solution is observable after the last wave has left the discretized region (approximately at time t ¼ 0:22 s). These tiny wave reflections are to be accredited to the used simplified concept of wave reflection reductions. The inaccuracy due to the adopted simplification occurs chiefly in the four corners of the discretized region, where the incoming waves incline from the assumed direction. 8. Singularity of the function f ðrÞ ¼ 1=r An undesirable property of the approximating function f ðrÞ ¼ 1=r is its singularity at r ¼ 0. The solution requires evaluation of the functional values also at poles where r ¼ 0. To avoid this difficulty could be to modify the approximation function to f ðrÞ ¼ 1=r þ 1. This way has

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

608

and the inverse of this matrix enters in the calculation of the mass matrix M). (2) If the prescribed value is too small, the solution is polluted by strong artificial oscillations, which do not exist in reality. The cause of this effect is also to be sought in the mass definition. By prescribing too small values of the inertial force approximating function, we practically define concentrated mass points. Such concentrated mass points, when disturbed and set into motion, do behave like harmonic oscillators and the surrounding media does not suffice to resist their excessive vibrations.

not been used in this work. The alternative used here is the substitution of the unbounded functional value by some prescribed empirical value. In the context of the present work it has been observed that the prescribed value has a strong impact on the solution. This impact is manifested in two aspects:

0.00E+00 -5.00E-06 -1.00E-05 -1.50E-05 -2.00E-05 -2.50E-05 -3.00E-05 -3.50E-05 -4.00E-05 -4.50E-05 -5.00E-05 -5.50E-05 -6.00E-05 -6.50E-05 -7.00E-05 -7.50E-05 -8.00E-05

The criterion used to determine the optimal prescribed value of the inertial force approximating function f ðrÞ ¼ 1=r at r ¼ 0 are the wave propagation velocities. These are

Solution incorporating wave reflection reduction measures Solution with no reflection reduction measures Static value

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45 0.5 0.55 Time (s)

0.6

0.65

0.7

0.75

0.8

Fig. 12. Time-history of the vertical surface displacements in the middle of the discretized region.

p(kN / m 2)

p(kN / m 2 )

1

1

t (s)

t (s)

-1

-1

p(t)

Model 1

Model 2 3m

p(t)

Zoomed in z

3m

y

p(t)

m

3m

m

x

32

3m

17

Vertical displacement (m)

(1) The prescribed values are directly proportional to the propagation velocities of all three types of elastic waves. The higher the prescribed value, the higher propagation velocities do the waves in the numerical solution exhibit and vice versa. From a physical point of view an artificially increased wave propagation velocity hints at an underestimated mass definition (recall that functional values of the inertial force approximating function are entries of the F matrix

17m

32m

Half-space

Half-space

Fig. 13. Used half-space models.

0.85

0.9

0.95

1

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

where Dx is the element edge length (m). Upon usage of this formula, all three wave propagation velocities (cR , cs , cp ) in the numerical solution match the true propagation velocities.

of course given by the material properties of the media and must not have anything to do with the numerical DRBEM solution. On the basis of many additional numerical experiments conducted in the context of this research a simple rule has been worked out, working seamlessly for regular meshes. The suggested formula for the calculation of the functional values of f ðrÞ ¼ 1=r at r ¼ 0 is f ðrÞðr¼0Þ 

19:5 , Dx

609

9. Numerical example—harmonically oscillating source The harmonic analysis offers a good opportunity for finding important rules for spatial discretization. Elastic waves induced by a harmonically oscillating source contain only one single frequency identical with the excitation

(20)

Table 3 Exact and obtained wavelengths f (Hz)

Wavelength (m) Exact

10 15 20 25

9.42 6.28 4.71 3.77

Relative error (%)

No. of elements per exact wavelength

Calculated Model 1

Model 2

Model 1

Model 2

Model 1

Model 2

9.4 6.7 5.6 4.8

9.4 6.3 4.8 3.9

0.0 6.6 18.9 27.3

0.0 0.0 0.0 3.4

9.42 6.7 4.71 3.77

18.84 13.4 9.42 7.54

Excitation frequency = 10 Hz 1.20E-05 Vertical displacement (m)

1.00E-05

Element size 1.0m x 1.0m

8.00E-06 R

6.00E-06

= 9.4m

Element size 0.5m x 0.5m

4.00E-06 2.00E-06 0.00E+00 -2.00E-06 R

-4.00E-06

= 9.4m

-6.00E-06 -3 -2 -1

0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 Distance from the source (m)

Fig. 14. Surface displacements at excitation frequency f ¼ 10 Hz.

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

610

Excitation frequency = 15 Hz

Vertical displacement (m)

1.20E-05 Element size 1.0m x 1.0m

1.00E-05

Element size 0.5m x 0.5m

8.00E-06 R

6.00E-06

= 6.7m

4.00E-06 2.00E-06 0.00E+00 -2.00E-06 R

-4.00E-06

= 6.3m

-6.00E-06 -3 -2 -1

0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 Distance from the source (m)

Fig. 15. Surface displacements at excitation frequency f ¼ 15 Hz.

frequency. Two half-space models are used in the subsequent investigation. They are portrayed in Fig. 13. Since the developed computer code was originally intended for transient analysis, at the current state it does not enable a genuine harmonic analysis. The loading in the following examples is defined as a sine wave for tX0 (see Fig. 13). The two half-space models use different discretizations. In the coarser model, elements of size 1:0 m  1:0 m are used. In the finer model the element size is 0:5 m  0:5 m. The time-step length used is Dt ¼ 0:005 s for the coarser model and Dt ¼ 0:0025 s for the finer model. The only well observable wave type on the halfspace surface are Rayleigh waves. Their wavelength lR is given by lR ¼ cR =f ,

(21)

where cR is the Rayleigh wave velocity and f is the excitation frequency in Hz. The value of the wave length lR is used as accuracy assessment criterion for the results. Results are calculated for four frequencies and presented in Table 3. Figs. 14–17 show snapshots with surface displacements for the calculated frequencies together with vertical cuts through the middle of the discretized region.

The obtained wavelengths are depicted as distances between two subsequent wave crests (or wave troughs) in these cuts. The obtained wavelengths cannot be read off exactly from the diagrams. Therefore a tolerance limit of 0:1 m was set for the evaluation of the relative error. Consequently, read-off values within this tolerance yield also relative errors of 0.0%. Results presented in Table 3 reveal an important finding. The relative error is zero when there are more than eight elements per wavelength. This finding matches the rule of thumb known from dynamic analysis with FEM, saying that at least eight elements per wavelength are needed to guarantee satisfactory accuracy of the analysis. Apparently, this rule is applicable to DRBEM just the same. Note also that the solution is practically devoid of wave reflections. 10. DRBEM applied to simulation of moving loads Moving loads became of great importance in half-space dynamics in context with the development of high-speed trains in the last decades. A moving load traveling on a free

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

611

Excitation frequency = 20Hz

Vertical displacement (m)

1.00E-05 8.00E-06

Element size 1.0m x 1.0m

5.6m

6.00E-06

Element size 0.5m x 0.5m

4.00E-06 2.00E-06 0.00E+00 -2.00E-06

4.8m

-4.00E-06 -6.00E-06 -3 -2 -1

0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 Distance from the source (m)

Fig. 16. Surface displacements at excitation frequency f ¼ 20 Hz.

half-space surface is a problem studied thoroughly by many researchers. In the opinion of the authors, the DRBEM has very attractive features with regard to the simulation of moving loads and the design of suitable track systems for high-speed railway conveyance. Timeindependence of the produced matrices and trouble-free incorporation of arbitrary initial conditions are surely advantageous in this numerical method. In this paper, we shall confine ourselves to a performance study of the method when applied to moving loads. In such studies it is a common practice to divide the whole range of load traveling velocities v, v 2 ð0; 1Þ into four distinct zones: (1) (2) (3) (4)

Subsonic velocities: vocR . Velocities in range: cR pvpcs . Transonic velocities: cs ovpcp . Supersonic velocities: cp ov.

The half-space model used in the study of moving loads is portrayed in Fig. 18. In the numerical solution, the moving load is applied suddenly at time t ¼ 0 s and moves at constant velocity v afterward. Exactly these conditions are assumed in the

analytical solution according to Payton [16]. In the cited paper, a closed-form solution is given for vertical displacement components in the whole range of the load traveling velocities for the special case of n ¼ 0:25. The analytical solution assumes a moving point force of 1 kN, whereas in the numerical solution the load is applied as a surface traction with intensity p ¼ 1 kN=m2 over an area of 1 m2 . The resulting force is thus also 1 kN. It would be inadequate to compare these two solutions directly under the load’s path. However, at some distance from the load’s path the singular analytical solution is smoothed and a comparison, at least qualitative, is possible. The half-space material properties together with the corresponding wave propagation velocities for all three possible types of waves are given in Table 4. The boundary element mesh covers a region with measures 72 m  20 m (in x- and y-direction). The model is meshed with a regular network of square boundary elements with edge lengths of 1 m  1 m. There are totally 1531 nodes (i.e. 4593 DOFs) and 1440 elements. Newmark’s time-integration scheme is used. Time-step length is accommodated to the load’s traveling velocity. A faster load requires comparatively shorter time-steps to ensure its smooth advance in space (the load should not pass more

612

Vertical displacement (m)

Element size 0.5m x 0.5m

6.00E-06

4.8m

4.00E-06 2.00E-06 0.00E+00 -2.00E-06 -4.00E-06 -6.00E-06 -3 -2 -1

0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 Distance from the source (m)

Fig. 17. Surface displacements at excitation frequency f ¼ 25 Hz.

ARTICLE IN PRESS

Element size 1.0m x 1.0m

8.00E-06

A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

Excitation frequency = 25 Hz 1.00E-05

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

than one element during one time-step). A certain amount of numerical damping is introduced by the Newmark’s time-marching scheme (Newmark’s parameter g ¼ 0:7). Surface displacements captured at representative timeinstances are portrayed in Figs. 19–26 supplied along vertical cuts made at a distance of 3 m from the load’s path. In contrast to the majority of analytical solutions to a moving point load on elastic half-space which assume stationary load traveling from 1 to 1, Payton’s solution assumes that the load is applied suddenly at time t ¼ 0, just like the numerical solution does. Transient effects due to impulsive character of the applied load are present in both solutions at the earlier phases of the simulation. Results obtained at four different load traveling velocities are presented.

io

n

p = 1kN/m2

Zoomed in

z

72 m

m

ov

in

g

di

re

ct

p = 1kN/m2

y x 22m Half-space

Fig. 18. Half-space model used for simulation of moving loads.

10.1. Subsonic velocity; v ¼ 50 m=s Transient dynamic effects are well observable at early stages of the simulation. All initially induced waves propagate faster than the load. After the initially induced waves have left the discretized region, no dynamic effects are to be seen in the solution and the displacement field strongly resembles a displacement field caused by a stationary moving load. The vertical funnel-shaped deflection has an elliptical form. A very good agreement between the analytical solution and the solution by DRBEM is observed. 10.2. Velocity in the range cR pvpcs ; v ¼ 96 m=s In this relatively narrow range of velocities, a steady amplification of displacements is detected. This is due to the effect commonly referred to as ‘‘critical velocity’’. It occurs if the load traveling velocity matches the propagation velocity of the Rayleigh wave in the half-space. The most striking difference between the analytical solution and the solution by DRBEM is the presence of a shallow wave trough outrunning the main Rayleigh wave crest observed in the numerical solution. This suggests that the so-called head wave or von Schmidt wave is present. It is caused by the dilatational p-wave propagating along the half-space surface. Equilibrium of forces does not allow the p-wave to exist alone on the surface and therefore its passage at any place on the surface causes new shear waves. In contrast to the analytical solution which exhibits sharp peaks and bends, the numerical solution is smoother. Steady growth of displacements is captured by the numerical solution as well. 10.3. Transonic velocity; v ¼ 150 m=s

Table 4 Used material properties Shear modulus Poisson’s ratio Density

613

G ¼ 17 MPa n ¼ 0:25 r ¼ 1700 kg=m3

cR ¼ 91:94 m=s cs ¼ 100:00 m=s cp ¼ 173:21 m=s

As the load traveling velocity increases the displacements decrease gradually. A secondary surface wave lagging behind the load is observed in both solutions. The head wave trough running before the main Rayleigh wave crest

Fig. 19. Displacement field at a load traveling velocity of v ¼ 50 m=s at time: (a) t ¼ 0:205 s and (b) t ¼ 0:885 s.

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

614

Loading velocity v = 50 m/s, t = 0.205s 1.20E-06 Rayleigh wave

Vertical displacement (m)

8.00E-07

DRBEM Payton

4.00E-07 0.00E+00 -4.00E-07 -8.00E-07

shear wave

-1.20E-06 -1.60E-06 -2.00E-06 -2.40E-06 -2.80E-06 0

2

4

6

8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 Position along the loading path (m) Loading velocity v = 50 m/s, t = 0.885s

0.00E+00 -2.00E-07 Vertical displacement (m)

-4.00E-07 -6.00E-07 -8.00E-07 -1.00E-06 -1.20E-06 -1.40E-06 -1.60E-06 -1.80E-06 -2.00E-06 -2.20E-06

DRBEM Payton

-2.40E-06 -2.60E-06 0

2

4

6

8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 Position along the loading path (m)

Fig. 20. Load traveling velocity v ¼ 50 m=s. Vertical displacements at a distance of 3 m from the load’s path at time: (a) at t ¼ 0:205 s and (b) t ¼ 0:885 s.

Fig. 21. Displacement field at load traveling velocity of v ¼ 96 m=s at time: (a) t ¼ 0:150 s and (b) t ¼ 0:400 s.

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

615

Loading velocity v = 96 m/s, t = 0.150s 1.00E-05 DRBEM Payton

Vertical displacement (m)

8.00E-06 6.00E-06 4.00E-06 2.00E-06 0.00E+00 -2.00E-06 -4.00E-06 -6.00E-06 -8.00E-06 0

2

4

6

8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 Position along the loading path (m) Loading velocity v = 96 m/s, t = 0.400s

1.00E-05 DRBEM Payton

Vertical displacement (m)

8.00E-06 6.00E-06 4.00E-06 2.00E-06 0.00E+00 -2.00E-06 -4.00E-06 -6.00E-06 -8.00E-06 0

2

4

6

8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 Position along the loading path (m)

Fig. 22. Load traveling velocity v ¼ 96 m=s. Vertical displacements at a distance of 3 m from the load’s path at time: (a) t ¼ 0:150 s and (b) t ¼ 0:400 s.

Fig. 23. Displacement field at a load traveling velocity of v ¼ 150 m=s at time: (a) t ¼ 0:150 s and (b) t ¼ 0:350 s.

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

616

Loading velocity v = 150 m/s, t = 0.150s 5.00E-06 DRBEM Payton

Vertical displacement (m)

4.00E-06 3.00E-06 2.00E-06 1.00E-06 0.00E+00 -1.00E-06 -2.00E-06 -3.00E-06 -4.00E-06 0

2

4

6

8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 Position along the loading path (m) Loading velocity v = 150 m/s, t = 0.350s

5.00E-06 DRBEM Payton

Vertical displacement (m)

4.00E-06 3.00E-06 2.00E-06 1.00E-06 0.00E+00 -1.00E-06 -2.00E-06 -3.00E-06 -4.00E-06 0

2

4

6

8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 Position along the loading path (m)

Fig. 24. Load traveling velocity v ¼ 150 m=s. Vertical displacements at a distance of 3 m from the load’s path at time: (a) t ¼ 0:150 s and (b) t ¼ 0:350 s.

Fig. 25. Displacement field at a load traveling velocity of v ¼ 220 m=s at time: (a) t ¼ 0:130 s and (b) t ¼ 0:285 s.

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

617

Loading velocity v = 220 m/s, t = 0.130s 5.00E-06 DRBEM Payton

Vertical displacement (m)

4.00E-06 3.00E-06 2.00E-06 1.00E-06 0.00E+00 -1.00E-06 -2.00E-06 -3.00E-06 -4.00E-06 0

2

4

6

8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 Position along the loading path (m) Loading velocity v = 220 m/s, t = 0.285s

5.00E-06 DRBEM Payton

Vertical displacement (m)

4.00E-06 3.00E-06 2.00E-06 1.00E-06 0.00E+00 -1.00E-06 -2.00E-06 -3.00E-06 -4.00E-06 0

2

4

6

8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 Position along the loading path (m)

Fig. 26. Load traveling velocity v ¼ 220 m=s. Vertical displacements at a distance of 3 m from the load’s path at time: (a) t ¼ 0:130 s and (b) t ¼ 0:285 s.

detected in the numerical solution is still present, but it has diminished noticeably comparing to the previous case. 10.4. Supersonic velocity; v ¼ 220 m=s The load is now faster than the fastest elastic wave. No disturbance is observed prior to the load’s passage. Diminishing of displacements accompanies further increase of the load’s traveling velocity. Formation of a Mach’s cone is easily seen. The secondary surface wave lags now markedly behind the load. 10.5. Displacements vs. load traveling velocity It is interesting to demonstrate what impact does the load traveling velocity have on the magnitude of the

displacements on the surface. Maximal positive (downward) vertical displacement due to passage of a moving load are detected on the half-space surface. These are portrayed in Fig. 27 as a function of the load traveling velocity. Distinctive velocities were computed, varying from v ¼ 0 to 260 m/s. Dots on the portrayed curve mark the calculated velocities. Results at intermediate velocities are obtained by linear interpolation. Wave propagation velocities cR , cs and cp are extra marked for a better assessment of the results. It is known that maximal displacements occur when the load traveling velocity matches the Rayleigh wave propagation velocity in the half-space (critical velocity). The obtained results are in agreement with the reality. The peak value occurs when the load traveling velocity matches the cR velocity.

ARTICLE IN PRESS A. Tosecky´ et al. / Engineering Analysis with Boundary Elements 32 (2008) 597–618

Max. vertica ldisplacement (m)

618

7.00E-06 6.50E-06 6.00E-06 5.50E-06 5.00E-06 4.50E-06 4.00E-06 3.50E-06 3.00E-06 2.50E-06 2.00E-06 1.50E-06 1.00E-06 5.00E-07 0.00E+00 0

10

20

30

40

50

60

70

80

90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 Load's traveling velocity (m/s)

Fig. 27. Maximal vertical displacement on the half-space surface vs. load traveling velocity.

11. Conclusions The application of the dual reciprocity boundary element method to half-space dynamics has been studied in the present paper. The presented examples were selected such that the numerical solutions could be compared with available analytical results. Two milestones have been laid in this work. New particular solutions to the Navier’s equation have been derived. Especially one of these solutions (f ðrÞ ¼ 1=r) has proved as very successful and versatile. It has been shown that this function provides correct wave amplitude decay in total absence of internal nodes. Furthermore, a novel strategy for wave reflection removal has been developed and successfully applied. Both the particular solution as well as satisfaction of the radiation conditions are essential for DRBEM applications to half-space dynamics. The authors are fully aware that it will take more effort and numerical evidences until the DRBEM application to half-space dynamics can be accepted as a reliable and robust method in the community of researchers and engineers. The promising results presented in this paper could draw the attention of other researchers and trigger a heightened interest in the DRBEM and its application to dynamics of infinite and semi-infinite media.

[3]

[4] [5] [6]

[7]

[8]

[9]

[10]

[11]

Acknowledgments

[12]

This work was financially supported by the German Research Foundation (DFG) under the Grant no. SCHM 546/17-1. The authors wish to express their sincere gratitude for this support.

[13]

References [1] Wolf JP, Song C. Finite-element modeling of unbounded media. New York: Wiley; 1996. [2] Pflanz G. Numerische Untersuchung der elastischen Wellenausbreitung infolge bewegter Lasten mittels der Randelementmethode in

[14]

[15]

[16]

Zeitbereich. PhD thesis, Ruhr University Bochum, Bochum, Germany; 2001 [in German]. Grundmann H, Trommer E. Transform methods—what they can contribute to (computational) dynamics? Comput Struct 2001;79: 2091–102. Kausel E. An explicit solution for the dynamic loads in layered media. Technical report R81-13, MIT; 1981. Kausel E. Thin-layer method: formulation in the time domain. Int J Numer Eng Mech 1994;37:927–41. Waas G. Linear two-dimensional analysis of soil dynamics problems in semi-infinite layered media. Doctor thesis, Berkeley: University of California; 1972. Partridge PW, Brebbia CA, Wrobel LC. The dual reciprocity boundary element method. UK: Computational Mechanics Publications & Elsevier Science Publishers; 1992. Tanaka M, Chen W. Dual reciprocity BEM applied to transient elastodynamic problems with differential quadrature method in time. Comput Methods Appl Mech Eng 2001;190:2331–47. Power H, Partridge PW. Time dependent Stokes flow for interior and exterior domains via DRBEM BEM XV, vol. 1, Worcester Massachusetts, USA; 1993. p. 31–341. Tosecky´ M. Wave propagation in homogenous elastic half-space using the dual reciprocity boundary element method. Dissertation at the Faculty of Civil Engineering, Ruhr University Bochum, October 2005; hwww.ub.ruhr-uni-bochum.de/DigiBib/DissListen/ Bauingenieurwesen/T.htmli. Bialecki RA, Jurgas P, Kuhn G. Dual reciprocity BEM without matrix inversion for transient heat conduction. Eng Anal Boundary Elem 2002;26:227–36. Dominguez J. Boundary elements in dynamics. Computational Mechanics Publications & Elsevier Science Publishers; 1993. de Medeiros GC, Partridge PW, Branda˜o JO. The method of fundamental solutions with dual reciprocity for some problems in elasticity. Eng Anal Boundary Elem 2004;28:453–62. Triantafyllidis T. Arbeitsbla¨tter fu¨r Baugrunddynamik. Technical report, Lehrstuhl fu¨r Grundbau und Bodenmechanik, Ruhr University Bochum, Bochum; 2004 [in German]. Partridge PW. Towards criteria for selecting approximating functions in the dual reciprocity method. Eng Anal Boundary Elem 2000;24: 519–29. Payton RG. An application of the dynamic Betti–Rayleigh reciprocal theorem to moving-point loads in elastic media. Quart Appl Mech 1963;21(4).