## Abstract

We study numerically the propagation of seismic waves in a three-dimensional block medium. The medium is modelled by a spatial lattice of masses connected by elastic springs and viscous dampers. We study Lamb's problem under a surface point vertical load. The cases of both step and pulse load are considered. The displacements and velocities are calculated for surface masses. The influence of the viscosity of the dampers on the attenuation of perturbations is studied. We compare our numerical results for the block medium with known analytical solutions for the elastic medium.

## 1. Introduction

Until recently, geomechanics and geophysics widely used the theory of deformation of rock masses which are treated as homogeneous media. In this theory, dynamical properties of the medium are described by the well-developed linear theory of elastic waves. On the basis of this theory, methods are developed for calculating the stress state of rocks near the mines as well as for processing and interpretation of seismological data in geophysics and mining. Serious reason to revise established views given in the papers [1–3], which show that, in the mathematical models of rock masses, it is necessary to take into account the block-hierarchical structure of the medium. Especially strong influence on the development of such an approach has been made by a fundamental concept proposed by Sadovskiy [1], according to which, rock masses are treated as systems of nested blocks of different scales, connected with each other by interlayers composed of weaker fractured rocks. Such a structure of the medium affects the process of wave propagation. The presence of the interlayers with weakened mechanical properties leads to the fact that the deformation of the block medium is mainly due to the deformation of the interlayers. As noted in [2], the structure of a block medium is the cause of various dynamic phenomena that are absent in a homogeneous medium and, therefore, cannot be described by its models. Among those dynamic phenomena we distinguish the propagation of the pendulum waves [2,3] having a low velocity of propagation (that is significantly less than the velocity of elastic waves in the blocks), long length (even under a pulse action) and weak damping.

In [4–9], we studied one-dimensional mathematical models of viscoelastic deformation of block media and have shown that the presentation of the blocks as massive rigid bodies allows us to distinguish from a complex dynamic state of the block medium that part of it, which is determined by the deformation of the interlayers between the blocks. The presentation of the blocks as massive rigid bodies made it possible to accurately describe the low-frequency pendulum waves arising under the impact load. In [4–9], it is shown that impulse load generates broadband vibrations in a block medium which, with the spread, is divided into high-frequency oscillations, corresponding to the eigen oscillations of the blocks, and low-frequency oscillations. Laboratory experiments on physical models of block media described in [3–7] have shown that high-frequency waves attenuate rapidly and the main seismic effects are caused by low-frequency waves. In contrast with the studies in [2–7], which focuses on the study of the results of laboratory and field experiments in block media, in [8,9], a theoretical analysis is carried out of the propagation and spectral characteristics of waves under the influence of the block-hierarchical structure of the rock mass.

The dynamic behaviour of a two-dimensional block medium was studied in [10–13]. In these papers, the same approach has been used that was employed in [9] for a one-dimensional block medium. In [10–13], it was assumed that the rectangular blocks interact with each other via compressible elastic interlayers. In [10], the blocks were assumed to be rigid, while in [11–13] they were assumed to be elastic. In [12,13], the equations of the orthotropic Cosserat continuum were also used for the description of the dynamic properties of a two-dimensional block medium.

A simpler model of a two-dimensional block medium can be obtained if the blocks are treated as point masses connected with each other by springs and dampers. In this case, a regular block medium can be represented as a periodic lattice of the masses. Different versions of this model were used in [14–22] in order to describe the plane and antiplane deformations of discrete media. However, in [22], the emphasis is on continuous models describing discrete microstructures.

Despite the obvious limitations of the usage of periodic lattices as models for the description of real-world block media, they have an advantage of making it possible to use analytical and numerical methods, as well as to describe qualitatively the dynamic effects inherent to these media. Initially the theory of waves in periodic structures appeared in the theory of crystal lattice dynamics [23–25], later it found applications in the mechanics of composites [26], but for a long time it was not used for description of the dynamic properties of rock masses. However, the study of wave propagation in periodic block structures is important for applications in seismology [1–3].

Below we draw the attention of the reader to some publications devoted to the propagation of seismic waves which are relevant to our study. For the first time, the problem of the dynamic impact of a vertical point force applied onto the surface of an elastic half-space was considered by Lamb [27]. In that paper he showed that, on the surface of the half-space, the Rayleigh waves [28] propagate along with the P- and S-waves. A short review of results on plane Lamb's problem can be found, e.g. in [18]. The spatial transient Lamb's problem and Rayleigh waves were studied in many monographs (e.g. [29,30]) and articles (e.g. [20,31–39]).

In [30], Cagniard obtained a solution to the three-dimensional Lamb's problem using the method of integral transformations. To calculate the inverse integral transforms, he proposed a method that now bears his name. Later that method was developed and generalized in [31]. In [32], Ogurtsov studied Lamb's problem for both vertical and horizontal impact onto the surface of an elastic half-space. In [33–35], the authors studied the propagation of seismic waves in an elastic half-space caused by a transient spherically symmetric source embedded into the half-space. In [33], the emphasis is on the theoretical calculations of the displacements on the surface of the half-space. In [35], analytical representations are obtained for the displacements and stresses in the Rayleigh wave on the surface of the homogeneous elastic half-space and P-wave inside it. In [36], the influence of inhomogeneity of the elastic half-space was studied. In [37], the analysis of the Rayleigh waves was conducted in the framework of the Cosserat continuum. In [38], the authors used the finite-element method for plane and spatial Lamb's problems for the case of harmonic loading. In [39], Kausel generalized analytical results obtained for Lamb's problem for an elastic medium for the cases of the vertical load, horizontal load and embedded source. The paper [39] contains many references related to the spatial Lamb's problem for an elastic medium. In [20], a solution is given to Lamb's problem in a discrete medium. That solution was obtained for two-dimensional and three-dimensional lattices in the form of multiple integrals, which are similar to the integral representations of the Bessel functions.

In the present paper, a block medium is modelled as a three-dimensional lattice of masses connected by springs and dampers in the axial and diagonal directions. For this model, we (i) study dispersion properties; (ii) solve Lamb's problem using a finite-difference method; (iii) study the degree of attenuation of numerical solutions to Lamb's problem on the surface of the block medium as a function of the distance from the impact point and the viscosity of the dampers; and (iv) compare that numerical solution with the results of analytical solutions for a homogeneous elastic half-space presented in [39].

## 2. Setting of the problem

We study the transient spatial Lamb's problem on the impact of a vertical point load on the surface of a half-space filled by a block medium.

The block medium is modelled by a uniform three-dimensional lattice consisting of the point masses, connected by springs and dampers in the directions of the axes *x, y, z*, and in the diagonal directions in the planes *xy, xz, yz*, as shown in figure 1. In that figure, *u*, *v* are the horizontal displacements in the directions of the axes *x*, *y*; *w* is the vertical displacement in the direction of the *z* axis; *n*, *m, k* (*x*, *y, z.* Our theoretical description of the deformation of interlayers is based on the rheological model by Kelvin–Voigt [40]. A vertical point load is applied suddenly onto the surface of the block medium (

Using the notation
*n, m, k* (

where *M* is the mass of a block; *x*, *y, z*;

In order to derive the equations of the motion at the boundary points, we consider that the stresses are equal to zero on the free surface. More precisely, using the notation
*n, m,* 0, located on the surface of the half-space in the following form:

where

The parameters *M*,

Since the planes

The initial conditions for equations (2.2) and (2.3) are supposed to be zero.

## 3. Dispersion properties

Let *l* be the length of a spring in the direction of the axis. In (2.2), replace each difference expression (2.1) by its differential approximation, e.g.

Assuming *l* small enough, equations (2.2) turn into the equations of the orthotropic elasticity theory:

If

Below we always assume

Using the above equation and (3.2), we calculate the velocity of the Rayleigh waves of an isotropic elastic medium with Poisson's ratio

In order to obtain the dispersion relations for the model described by equations (2.2), we apply to each of these equations the Laplace transform (of parameter *p*) with respect to time *t* (denoted by the superscript L), and the discrete Fourier transforms (of parameters *n*, *m*, *k*, respectively (denoted by the superscripts

The Laplace–Fourier transform of equation (2.2) is as follows:

Substituting *A* is a matrix whose entries *J* is the identity matrix, we obtain the dispersion equation for determining the frequency

The roots of the cubic equation (3.4) can be found by Cardano's formula.

Let us investigate the dispersion equation (3.4) in the special case

Let

The dispersion equation (3.5) has the following roots:

and

Since *A* is symmetric, all its eigenvalues are real. Numerical calculations show that all the principal minors of the matrix *A* are non-negative for

The vectors of the phase and group velocities and their norms are as follows:

Denote by *c*_{f,j} the vector of the phase velocity corresponding to the frequency

Note that, as might be expected, formulae (3.2) and (3.10) coincide with each other.

For arbitrary values of

Numerical calculations show that the vector of the group velocity is equal to zero at the following points:

At these points, the frequencies

It follows that there are only three different values of the eigenfrequency:

Figure 3 depicts graphs of the frequency versus wavenumbers

## 4. Stability of the difference scheme

We solve equations (2.2) and (2.3) with conditions (2.4) and zero initial conditions by a finite-difference method using an explicit scheme. For the second derivatives with respect to time, we use the second central difference approximation of the second order of accuracy:

For the first derivatives with respect to time, we use the backward difference approximation of the first order of accuracy:
*s* is the number of the time step in the finite-difference scheme. We use the same approximations for the displacements

The stability condition of the difference equations, corresponding to equation (2.3), is as follows:

Thus, the stability condition of the difference equations of Lamb's problem for the block medium is as follows:

Numerical calculations confirm the validity of this stability condition.

## 5. Results of numerical experiments: the case of a step load

In this section, we present the results of numerical calculations of disturbances in Lamb's problem for a lattice under the action of the vertical step load

In [41], an analytical solution is given to static spatial Lamb's problem for an isotropic elastic medium which is known as Boussinesq's problem. In our notation, the solution of the Boussinesq's problem, given in [41], is as follows on the surface *E* is the Young modulus and

For

Figure 5 shows the time dependence of the radial *w* displacements and their velocities *w*, shown in figure 5 (as well as in figure 7), the horizontal dashed lines correspond to the static solution (5.1). The vertical dashed lines in figure 5 (as well as in figures 7–11) correspond to the arrival time of the quasi-fronts of the longitudinal, shear, and Rayleigh waves at the point with the coordinates

As can be seen from figure 5, the displacements *w* behind the front of the Rayleigh wave oscillate, with frequency

More generally, we can say that numerical dispersion appearing in the transition from the equations of the theory of elasticity to their difference analogues is of the same nature as dispersion inherent to the discrete periodic media. However, we have to reduce parasitic effects of the discretization, while we have to take into consideration dispersion inherent to the block medium and to study its impact on the wave process.

The solid curves in figure 6 depict the plots of the maximal amplitudes of the perturbations *w*, *n* calculated for *A* and *f* at the points *w*, *f* are written next to the curve. Figure 6 shows that the maximal amplitude of the velocities *n*, tends to zero as *n* is increasing. Moreover, the rate of the descent increases with the increasing of the viscosity of the interlayers.

Thus, figures 5 and 6 show that the dissipative properties of the interlayers result in additional damping of the low-frequency pendulum waves and in more rapid decay of the high-frequency waves, which is more consistent with real seismograms.

In table 1, we present the power-law dependences of the damping of the amplitudes of the disturbances versus coordinate *n*, corresponding to

Figure 7 shows plots of the quantities *w*,

In the electronic supplementary material to [39], a program is given for calculating the vertical and radial displacements on the surface of an elastic half-space. Using that program, we calculate the vertical and radial displacements, *w* and

The results of the calculations of the vertical *w* and radial

## 6. Results of numerical experiments: the case of an impulse load

In this section, we present the results of numerical calculations of disturbances in Lamb's problem for a block medium under the action of the vertical impulse load given by the formula

Figures 8–11 show the plots of the radial *w* displacements and their velocities

We use Duhamel's integral

Figures 8–11 show that the qualitative and quantitative agreement between the results obtained for the block and elastic media are improved as the duration of the impulse (6.1) increases. The best agreement takes place for the displacements *w* (figures 8 and 9).

The results of numerical calculations for loads which have time dependencies different from those considered in §§5 and 6 can be found in the electronic supplementary material.

## 7. Conclusion

In this paper, we proposed a three-dimensional mathematical model of block rock masses. This model is based on the idea that the dynamic behaviour of a block medium can be roughly described as the motion of rigid blocks due to compressibility of interlayers between them, and that the deformation of the interlayers can be approximately described by Kelvin–Voigt's model. Using this model, we solved Lamb's problem numerically, specifically, we studied the propagation of seismic waves in a block medium under the effect of a vertical transient point load applied at the surface of the half-space. For Lamb's problem, we have shown that the main contribution to the wave process on the surface of the block medium is made by low-frequency waves in the vicinity of the front of the Rayleigh wave and that high-frequency waves are observed behind the front of the Rayleigh wave.

It is shown that the presence of the block structure in the medium leads to the following changes in its behaviour in comparison with what the model predicts of a homogeneous elastic medium, whose mechanical properties are obtained by averaging the mechanical properties of the block medium:

— In the block medium, waves propagate with dispersion, which is absent in the homogeneous elastic medium.

— On the surface of the block medium, low-frequency longitudinal and Rayleigh waves propagate with velocities which are much smaller than the corresponding velocities in the blocks.

— The velocity of propagation of low-frequency waves and the degree of attenuation is determined by the weight of the blocks, their dimensions and properties of interlayers.

— The dissipative properties of the interlayers result in additional damping of the low-frequency waves on the surface of the half-space in the vicinity of the front of the Rayleigh wave and in more rapid decay of the high-frequency waves behind the front of the Rayleigh wave.

— The difference in the behaviour of the block medium from the behaviour of the homogeneous elastic medium is especially considerable for the case of a short impulse. When the impulse duration is increased, the difference in the behaviour of the block medium from the behaviour of the homogeneous elastic medium becomes less noticeable.

The results obtained in this article show that in the study of seismic waves, it is necessary to take into account the block structure of the rocks and rheological properties of the interlayers.

## Competing interests

I declare I have no competing interests.

## Funding

This work was not supported by a grant.

- Received February 16, 2016.
- Accepted July 15, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.