Introduction. In the 70s of the last century, the finite element method (FEM) began to be extensively used in the analysis of the static stability of slopes and flanks [1–3]. According to the practice of geotechnical calculations, FEM, in contrast to simplified methods, provided taking into account such important factors as the real geometry and layered structure of the objects under study, as well as the presence of anti-landslide structures and the physical nonlinearity of the material with a previously unknown area of plasticity [4]. Currently, many researchers start using FEM to analyze the dynamic stability of actual slopes, as well as forecasting the risks of rock falls in the underground workings and from slopes located along highways and railways. Also, one of the topical areas of mountain dynamics is finite element modeling of earthquake consequences, taking into account the kinematic instability of conglomerations [5]. Today, there are the following methods of dynamic calculations in the field of geotechnics:

1. SRM (strength reduction method) is a method of reducing strength [6, 7]. It is designed to calculate the safety margin of a mountain range in a physically nonlinear formulation. In SRM, the actual strength parameters of the soil are iteratively reduced through dividing by some factor greater than 1:

where c, ф — actual values of the adhesion and the angle of internal friction of the soil, respectively; — the adhesion and the angle of internal friction of the soil, respectively, after their reduction relative to the actual values; kSRM — the shear strength reduction factor. Value kSRM, corresponding to the limiting state of the system, determines the lower limit of the strength parameters of the material.

2. LEM (limit equilibrium method) is a method of limiting kinetics based on D'Alembert's principle [8][9]. It is focused on the analysis of dynamic stability of articulated rock massifs.

3. TLEM (thin layer element method of FEM) is a thin-layer finite element method [10] in which elastic-plastic elements of a thin layer are used to model the behavior of kinematically unstable structures.

Analysis of the results obtained using the SRM, LEM and TLEM methods showed that there is currently no unified vision of mathematical modeling of the behavior of structurally unstable geotechnical systems under nonstationary external influence. This determines the topicality of developing a methodology for dynamic analysis of systems of the “base — weakened layer — block” type in a finite element formulation through a new approach to simulating sliding planes.

Materials and Methods. The equation of motion of a mechanical system in a finite element formulation can be given as follows [11]:

(1)

where [M], [C], [K] — matrices of masses, damping and stiffness of the ensemble of finite elements, respectively;

— vectors-columns, respectively, of nodal accelerations, velocities, displacements;  — vectors-columns of specified static and dynamic loads, respectively, at a time t. In the future, we assume that the matrices [M] and [K] are consistent.

For the numerical integration of equation (1), we use Newmark method [12]. Here, we assign the integration step along the time axis t so that the contributions of physically significant proper pairs are considered with sufficient accuracy. In the future, we will consider the kinematic methods of excitation of vibrations, set using either model seismogram , or model accelerogram . With this method of setting the dynamic effect, the second term of the right side of equation (1) will be zero:

Consider the method of excitation of mechanical vibrations from a model seismogram. Function can be written as [13]:

, (2)

where A — initial amplitude;  — attenuation coefficient; — angular frequency of external influence. Figure 1 shows a graph of function for the values: A = 0,01553 m; = 0,7143; = 5 s–1.

Research Results. As the first model example, consider the problem on forced vibrations of the slope located on the base (Fig. 2). The boundary conditions of the problem are shown in Figure 3, where letter S denotes the point under study.

Relationships between the geometric parameters of the slope and the base (Fig. 2) are presented in Table 1.

Table 1

Slope – Base Geometrics Relationship

Mechanical characteristics of the slope and base material are ae follows: deformation modulus Е = 21 MPa; Poisson's ratio v = 0.3; specific gravity y = 1702 кг/м3; adhesion с = 45 кПа; internal friction angle ф = 15°.

To model the slope and the base, we use SOLID45 volumetric finite elements of the ANSYS Mechanical software package. The finite element model for the variant with the parameters: ls = 2Hs, Hs = 10 m, assigned to the global Cartesian coordinate system, is shown in Figure 4.

The finite element grid is constructed so that on the contact surface, adjacent nodes of the base and slope have the same coordinates, but different numbers. This is done in order to arrange a weakened layer in this place. The kinematic effect in the form of a model seismogram (2) is set at each integration step ti in the form of nodal displacements on the end surfaces of the model with parameters: X = 0 and X=l1 + ls + l2.

Modeling of the weakened layer (Fig. 4) is performed using elastic-viscous combined finite elements COMBIN14 [14]. The two-node element COMBIN14, consisting of a spring with stiffness k and a liquid friction damper with a damping coefficient сv, is shown in Figure 5. In the case under consideration, this element works only for tensioncompression.

In each node of the contact surface (Fig. 4), along the global X, Y, Z axes, we introduce elements of COMBIN14. Parameters of combined elements are:

k = 30 кН/м; ky = kz = 9,44∙10 7 кN/m; сv = 0,5.

In this example, we further introduce the assumption of the natural undeformed state of the “base — weakened layer — slope” system. For calculations, we use the nonlinear solver of the ANSYS Mechanical complex.

The finite element modeling results in the form of visualization of the deformed state of the “base — slope” system with the account of the maximum horizontal displacement and the distribution of amplitude horizontal displacements Wx(t) are shown in Figures 6 and 7. The integration step of equation (1) t = 0,01 s. As can be seen, the introduction of 3D elastic-viscous elements makes it possible to simulate the effect of kinematic instability of the “base — weakened layer — slope” mechanical system with the kinematic method of excitation of vibrations.

The amplitude value of the displacement at point S was WxS max  = 1,7 cm. For the slope option ls = Hs (fig. 2) WxS max = 1,1 cm. The graphs of base and slope vibrations at the studied point S (Fig. 3) in the direction of X-axis are shown in Figure 8.

Based on the above graphs, it can be seen that starting from the moment of time t > 1.5 s, there is a mismatch of base and slope vibrations.

Let us consider the behavior of the “base — weakened layer — slope” system (Fig. 3) when vibrations are excited using a model accelerogram . To this end, we differentiate expression (2) twice. As a result, we get:

(3)

Figure 9 shows the graph of function (3) for parameters: A = 0,01553 m; X = 0,7143; = 5 s–1.

The kinematic effect in the form of model accelerogram (3) , by analogy with seismogram (2), is set at each integration step ti in the form of nodal accelerations on the end surfaces of the model with parameters: X = 0 and X=l1 + ls + l2. Figure 10 shows the graphs of vibrations at the studied point S (Fig. 3) under kinematic action in the form of a model accelerogram. Comparing the vibration graphs shown in Figures 8 and 10, we establish that they almost coincide. This indicates the correctness of the developed finite element model, which allows describing the behavior of the “base — weakened layer — slope” system with various methods of unsteady kinematic action.

As a second model example, let us consider the problem on forced slope vibrations with a kinematically unstable wedge-shaped inclusion (Fig. 11). Due to the symmetry of the configuration, only 1/2 part of the slope and inclusions are taken into account in the computational scheme. The boundary conditions for the accepted design scheme are shown in Figure 12. Here, letter S denotes the point under study, which belongs simultaneously to the slope base and the wedgeshaped inclusion.

The finite element model of slope and wedge-shaped inclusion is shown in Figure 13. As in the previous example, in this case, we use SOLID45 and COMBIN14 elements with the same material characteristics.

The function describing the model accelerogram has the form:

where A — acceleration amplitude; — frequency of external influence. Figure 14 shows graph  at A = 2,5 m/s2, = 2 Hz. The values of acceleration at the i-th step of integrating the equation of motion (1) are applied to the nodes of the model surface with coordinate X = 0 (Fig. 12).

The simulation result in the form of a distribution of the amplitude values of displacements Wx (t) is shown in Figure 15. Integration step t = 0,01 s. Vibration graphs of the slope base and wedge-shaped inclusion at the studied point S (Fig. 12) in the direction of X-axis are shown in Figure 16.

As can be seen from Figure 15, with a given kinematic effect, the continuity of the slope array breaks along the weakened layer, and the wedge-shaped inclusion shifts relative to the slope base along X-axis.

“Drift” in Figure 16 is due to the fact that this finite element model has no connections that prevent displacements along X-axis. As shown in [11], it is possible to solve the “drift” problem by subtracting the displacement of the slope base, which represents the displacement “as a rigid whole”, from the displacement values  and . It should be noted that the obtained amplitude values of the displacements provide evaluating the dynamic parameters of the “slope — weakened layer — wedge-shaped inclusion” system.

Conclusion. A finite element model has been developed and validated to study the dynamic behavior of kinematically unstable slopes in a three-dimensional formulation, taking into account the physical nonlinearity of the material.