CN112069719A - Three-dimensional microscopic discrete element analysis method for rockfall impact cushion - Google Patents

Abstract

The invention discloses a three-dimensional microscopic discrete element analysis method of a rockfall impact cushion, which comprises the following steps: the sand particles in the cushion layer are assumed to be spherical particles; putting a plurality of spherical particles one by one in a designated space to generate an initial cushion layer model; compacting the initial cushion layer model under the conditions of gravity, vibration and boundary constraint; solving the motion equations of all spherical particles in the compacting process; judging whether compaction is finished according to the solution result of the motion equation to obtain a compact cushion layer model; and carrying out a falling rock impact cushion test on the compact cushion model to obtain a three-dimensional microscopic discrete element analysis result of the falling rock impact cushion. According to the method, on the basis of the initial cushion layer model, the initial cushion layer model is compacted under the action of gravity and vibration to obtain the compacted cushion layer model, and the sand cushion layer mesoscopic model with controllable density is generated.

Description

Three-dimensional microscopic discrete element analysis method for rockfall impact cushion

Technical Field

The invention belongs to the field of rockfall protection engineering, and particularly relates to a three-dimensional microscopic discrete element analysis method for a rockfall impact cushion layer.

Background

The discrete cell method can simulate the processes of extrusion, collision and friction among cushion material particles from a microscopic scale, and observe more detailed simulation results. The currently mainstream commercial discrete element analysis software includes PFC 2D/3D, EDEM, etc., and open source codes YADE, Esys-Particle, LIGGGHTS, etc. The current numerical simulation research of rockfall impact cushions also has a plurality of defects: most numerical simulation researches on falling rock impact cushions adopt a homogeneous finite element model to simulate the stress deformation of the cushions, but few researchers pay attention to the resistance change rule on the contact surfaces of the falling rocks and the cushions; at present, commercial software developed based on a discrete element method is abundant, but a model which can be used for simulating and analyzing a rockfall impact cushion layer is few; because the difference between the particle size of the cushion particles and the diameter magnitude of the rockfall is large, a cushion discrete element model with a large number of particles needs to be established, and the problem of contact judgment among a large number of particles is difficult to solve by the current contact search algorithm; at present, the analysis on the mesoscopic mechanism of the rockfall impact cushion layer is few, and the influence of factors such as strength, viscosity and density of the cushion layer material on the impact effect cannot be explained according to mesoscopic mechanical characteristics such as extrusion, collision and friction among particles in the cushion layer.

Disclosure of Invention

Aiming at the defects of the prior art, the invention aims to provide a three-dimensional microscopic discrete element analysis method of a rockfall impact cushion layer, so as to solve the problem that the microscopic mechanical analysis among particles in the cushion layer is less in the prior art.

In order to solve the technical problems, the technical scheme adopted by the invention is as follows:

a three-dimensional microscopic discrete element analysis method for a rockfall impact cushion comprises the following steps: the sand particles in the cushion layer are assumed to be spherical particles; putting a plurality of spherical particles one by one in a designated space to generate an initial cushion layer model; compacting the initial cushion layer model under the conditions of gravity, vibration and boundary constraint; solving the motion equations of all spherical particles in the compacting process; judging whether compaction is finished according to the solution result of the motion equation to obtain a compact cushion layer model; and carrying out a falling rock impact cushion test on the compact cushion model to obtain a three-dimensional microscopic discrete element analysis result of the falling rock impact cushion.

Further, the information of the spherical particles comprises the radius, the mass and the moment of inertia of the particles, the spatial position, the translation speed and the rotation speed of the particles.

Further, the designated space is a cylindrical space; the boundary of the designated space is a rigid boundary.

Further, the number and the particle size distribution of the spherical particles are obtained through assignment of a probability density distribution function.

Further, a plurality of spherical particles form a particle system in a designated space; the state change of the particle system is respectively described by a global coordinate system and a local coordinate system; describing contact search and collision among all spherical particles in the particle system through a global coordinate system; the translation and rotation of each spherical particle in the particle system is described by a local coordinate system.

Further, the method for establishing the motion equation comprises the following steps: detecting contacted spherical particles by a quick contact search algorithm; calculating a contact force between the spherical particles which are contacted according to the contact collision model; and establishing a corresponding motion equation according to the resultant force and resultant moment of the contact force applied to each spherical particle. Further, the quick contact search algorithm includes: in the global coordinate system, equally dividing the global space where the spherical particles are located into a plurality of small spaces; expanding the small space to obtain an expanded space; calculating the overlapping distance between two spherical particles according to the sphere center distance and the radius between the spherical particles; respectively searching the overlapping distances between all particles in the small space and the expansion space; judging whether all spherical particles are contacted or not according to the overlapping distance; the expanded space comprises spherical particles in a local space and spherical particles on two sides of a space boundary adjacent to the spherical particles.

Further, the calculation formula of the overlap distance is as follows:

wherein the content of the first and second substances,

is an overlap distance, r_i、r_jRadius of particles i and j, respectively, l_ijIs the interparticle distance.

Further, the contact collision dieThe type is as follows: f_ij＝F_nij+F_tij，

Wherein, F_ijAs contact force between particles, F_nijFor normal contact force, F_tijIs the tangential contact force.

Further, the formula for calculating the normal contact force is as follows:

wherein the content of the first and second substances,

as a local coordinate system L_iThe contact force in the lower normal direction is,

in order to be the elastic force,

for viscous damping, k_nThe stiffness of the hook spring is shown,

to a depth of mutual invasion, c_nThe damping coefficient is represented by a coefficient of damping,

expressed in a local coordinate system L_iNormal relative movement velocity, E, of lower particle j at contact point into particle i_i、E_jThe modulus of elasticity, r, of particles i and j, respectively_i、r_jRadius of particles i and j, respectively, l_ijThe distance between the particles is the distance between the particles,

is a unit direction vector of the centroid of particle j pointing to the centroid of particle i, (r)_i+r_j-l_ij) Indicating the overlap distance, ξ_nIs a dimensionless coefficient, m^*Represents equivalent mass, m_iAnd m_jMass of particles i and j, respectively, e_nIs a coefficient of normal velocity recovery,

is the velocity of the movement of particle j relative to particle i,

is the absolute value of the particle jIn the case of the speed,

is the absolute velocity of particle i.

Further, the calculation formula of the tangential contact force is as follows;

wherein the content of the first and second substances,

representing a local coordinate system L_iContact force in the down-cut direction, mu_sDenotes the maximum coefficient of static friction, μ_dThe coefficient of dynamic friction is expressed as,

is the tangential velocity of particle j relative to particle i at the point of contact,

and

respectively represent local coordinate coefficients L_iAnd L_jThe angular velocities of the particles i and j below each around their centroid,

and

respectively representing the vector radii of the centroids of particle i and particle j to the point of contact,

is the tangential velocity of particle i relative to particle j at the point of contact.

Further, the resultant force of all contact forces experienced by any spherical particle is:

the resultant moment experienced by any spherical particle is:

wherein the content of the first and second substances,

representing a local coordinate system L_iContact force of the lower j-th particle on the i-th particle, N_iIndicates the number of particles coming into contact with the ith particle,

indicating the resultant of all contact forces to which the particle i is subjected,

the contact moment of the jth particle on the ith particle is shown,

representing the resultant moment of all forces to which particle i is subjected.

Further, the formula of the contact moment between the spherical particles is as follows:

wherein the content of the first and second substances,

as a local coordinate system L_iThe contact force in the downward tangential direction is,

is the vector radius.

Further, in the local coordinate system, the motion equation of each spherical particle is

Wherein m is_iWhich is indicative of the mass of the particles,

representing the translational acceleration of the center of mass of the particle,

indicating that particle i is subjected to the combined forces of all interparticle contact forces,

representing the external force action of gravity, vibration load and the like,

representing the rotational acceleration of the particle i around the center of mass,

representing the resultant moment of particle i subject to all interparticle contact forces,

showing the action of external moment such as vibration load,

representing the moment of inertia of the particle.

Further, the calculation formula for the inertia moment of the spherical particles is as follows:

wherein the content of the first and second substances,

denotes the moment of inertia of the particle, m_iMass of spherical particles, r_iThe radius of a spherical particle.

Further, the solving process of the motion equation comprises:

acquiring initial conditions of cushion layer particles, wherein the initial conditions comprise mass, moment of inertia, centroid position, translation speed, rotation angle and rotation speed; searching particles which are contacted with each other in each time step according to a quick contact search algorithm; calculating contact force resultant force and resultant moment between each particle and particles in mutual contact according to the initial conditions; calculating the translation increment, the translation speed increment and the rotation speed increment of the mass center of each particle in the time step; under the local coordinate system, the time step is continuously updated, and the steps are repeated for all the particles until the movement reaches the termination condition.

Further, the calculation formula of the centroid shift increment is as follows:

the calculation formula of the translational velocity increment is as follows:

the calculation formula of the rotation speed increment is as follows:

wherein the content of the first and second substances,

in order to translate the increment of the center of mass,

representing the translational acceleration of the centroid of the particle, at is the time step,

as local coordinate coefficient L at time t_iThe lower particle i is subjected to the resultant of all contact forces,

m represents the external force action such as gravity, vibration load and the like at the moment t_iIs the mass of the spherical particles,

representing the resultant moment of contact force experienced by particle i at time t,

showing the action of external moment such as vibration load at the moment t,

representing the moment of inertia of the particle.

Further, the contact searching method of the falling rock particles and the spherical particles in the compact cushion layer in the process of impacting the cushion layer by the falling rock comprises the following steps: establishing a spherical domain space with a certain radius in a local coordinate system of the rockfall particles; calculating the distances between the falling rock mass center and all the particles in the compact cushion layer; judging whether the particles are in the spherical space or not according to the distance; and carrying out contact search on the particles in the spherical domain space and the falling rocks, and judging whether the falling rocks are contacted with the spherical particles.

Further, the radius of the spherical domain space is:

R_s＝R₀+α·r_m (2)

wherein R is₀Denotes the falling rock radius, r_mThe maximum radius of all cushion layer particles is shown, and the suggested value range of alpha is 1-2.

Compared with the prior art, the invention has the following beneficial effects: on the basis of the initial cushion layer model, the initial cushion layer model is compacted under the action of gravity and vibration to obtain a compacted cushion layer model, and a sand cushion layer mesoscopic model with controllable density is generated; an equation of motion of the spherical particles is established according to the rapid contact search calculation and the contact collision model, a global contact search algorithm among the particles is improved, and the contact judgment efficiency is greatly improved; an inter-particle normal and tangential contact force model based on linear spring damping is introduced, and the resistance on a contact surface and the extrusion and shearing action among particles are analyzed in detail on a microscopic level.

Drawings

FIG. 1 is a diagram of spherical particle models and model information in a system coordinate system;

FIG. 2 is a diagram of a model of an initial mat layer randomly placed in a designated space;

FIG. 3 is a diagram of a method for determining contact between spherical particles;

FIG. 4 is a schematic diagram of a planned space method local contact search;

FIG. 5 is a schematic diagram of a local contact search by a sphere space method;

FIG. 6 is a schematic diagram of a spherical inter-particle normal contact force model;

FIG. 7 is a model diagram of normal contact force;

FIG. 8 is a schematic diagram of a spherical inter-particle tangential contact force model;

FIG. 9 is a model diagram of tangential contact force;

FIG. 10 is a comparison of discrete element models before and after compaction of a sandy soil mat;

FIG. 11 is a schematic diagram of a discrete element model particle system cycle calculation;

FIG. 12 is a schematic view of a spherical rockfall impact cushion test device such as Calvetti;

FIG. 13 is a gravel gradation curve for test mats such as Calvetti;

FIG. 14 is a discrete element model of a rockfall impact pad;

FIG. 15 is a comparison and verification of falling rock impact force time course curve test data and DEM numerical value results;

FIG. 16 illustrates a DEM model contact force and resistance calculation method;

FIG. 17 shows compressive and shear stresses (v) of the pad layers in different regions of the DEM model₀＝14.00m/s)；

FIG. 18 is a comparison of the cushion cell compressive and shear stress calculations (v)₀＝14.00m/s)；

FIG. 19 is a comparison of the viscous force, dynamic resistance and total resistance of the mat;

FIG. 20 is a DEM model shim interparticle coefficient of friction sensitivity analysis.

Detailed Description

The present invention is further illustrated by the following figures and specific examples, which are to be understood as illustrative only and not as limiting the scope of the invention, which is to be given the full breadth of the appended claims and any and all equivalent modifications thereof which may occur to those skilled in the art upon reading the present specification.

The sandy soil cushion layer is essentially composed of discrete particles, so that the stress deformation characteristic of the sandy soil cushion layer is determined by the extrusion deformation of the particles and the contact force among the particles, and the relation between the compressive resistance stress on a contact surface and the cushion layer strength, viscosity and density impact factors can be analyzed on a microscopic level by introducing a normal contact force calculation model and a tangential contact force calculation model into an inter-particle contact model. The method combines the characteristics of the rock fall impact cushion layer model, improves a contact search algorithm and a particle compaction method, and analyzes the resistance change on the rock fall contact surface at a repeated point.

A three-dimensional microscopic discrete element analysis method for a rockfall impact cushion comprises the following steps: the sand particles in the cushion layer are assumed to be spherical particles; putting a plurality of spherical particles one by one in a designated space to generate an initial cushion layer model; compacting the initial cushion layer model under the conditions of gravity, vibration and boundary constraint; solving the motion equations of all spherical particles in the compacting process; judging whether compaction is finished according to the solution result of the motion equation to obtain a compact cushion layer model; and carrying out a falling rock impact cushion test on the compact cushion model to obtain a three-dimensional microscopic discrete element analysis result of the falling rock impact cushion.

The method specifically comprises the following steps:

(1) generating a large number of discrete sand cushion particle models according to the characteristics of the quantity, the particle size and the like of the sand cushion particles; sand in beddingThe shape of the soil particles is typically oval particles of random shape, but considering that the size of the rockfall and bedding is very different from the particle size, this results in a bedding model that will contain a large number of sand particles, and therefore too fine a simulation of the particle shape would greatly increase the computational load, and the prior art shows that spherical particles can reflect the dynamics of a sand bedding. The spherical particle model assumes that each ovoid particle i can be reduced to have a certain radius R_iAnd the model only needs to consider mass m_iAnd moment of inertia I_iTwo kinematic parameters, considering that the contact deformation between the particles is very small, the shape of the spherical particles after being deformed by force cannot be described additionally, so the spherical particles are assumed to be rigid bodies. Therefore, compared with other particles in the shape of polyhedron, the calculation of contact search and contact collision between spherical particles is simpler, the memory is less, and the cushion material particles of the invention are simulated by rigid spheres.

Since the sand cushion particles in the test of the present invention are spherical in shape, the rockfall model is also assumed to be spherical in the discrete element analysis model and can be regarded as spherical particles having a large particle size and different material characteristics from the cushion, and therefore, the inter-particle contact search and contact collision model is also used to solve the problem of contact between the rockfall model and the cushion particles.

In the sand cushion layer, the sand particles are different in size, generally meet a certain distribution probability, and can be obtained through a grading curve of particle sizes. The invention uses a spherical particle model to simulate sand particles, so that the diameter d of the spherical particles_iCan be regarded as the equivalent particle size of the sand particles and then the diameter d by a probability density distribution function according to the same grading curve or other forms (e.g. uniform distribution, normal distribution, etc.)_iAnd assigning values to generate a cushion layer particle model meeting the requirements of particle quantity and particle size distribution. As shown in fig. 1, the change of state of the particle system can be described by a global coordinate system G and a local coordinate system L, respectively, since the global coordinate system is more suitable for processing contact search and collision between all particles in the particle system and the local coordinate system is more suitable for processing translation and rotation of each particleAnd (6) moving.

The invention assumes that the spherical particle model i has a three-dimensional spatial position vector

Translation velocity vector

And rotational velocity vector

Wherein L is_iThe local coordinate system of the ith particle is shown. Therefore, each three-dimensional spherical particle has 3 translational degrees of freedom and 3 rotational degrees of freedom simultaneously, and the total number of degrees of freedom is 6.

In summary, all spherical particle model information includes the radius { R ] of the particle_i}, mass { m_iAnd moment of inertia

And spatial position

Translation speed

And speed of rotation

Wherein i is 0, 1, 2 … N _c0 denotes falling rock particles, 1, 2 … N_cShowing the mat pellets.

(2) Putting the cushion layers one by one in a designated space to generate an initial cushion layer model; in order to obtain a compact cushion discrete meta-model, particles are firstly randomly put in a specified space to generate an initial cushion model. Generally, the cushion structure is composed of walls and bottom supports surrounded by barriers, so that the discrete element model of the cushion needs to be provided with corresponding boundary conditions before being put into use. The boundaries of the discrete element model are generally divided into flexible boundaries and rigid boundaries: the flexible boundary is formed by flexibly connecting flaky particles and can be deformed; a rigid boundary is a plane that can collide with particle contact, but is not itself deformable.

The present invention employs a rigid boundary assumption, since a rigid boundary can be considered as a radius R_iThe model of contact collision between the rigid boundary and the particles can also be regarded as a model of contact collision between particles, and therefore, the inter-particle contact search and contact collision model is also used for solving the problem of contact between the boundary and the particles.

After the particle placement region is determined and the boundary conditions are set, the generated cushion particle models can be randomly placed into the designated space one by one. In particular, assume that a given cylindrical pad dimension is known, with a radius R_cThickness of h_cFirst, N is randomly generated in a designated space_cIndividual position coordinate

(i＝1、2…N_c) Wherein the designated space is a radius R_cHeight of α h_cThe cylindrical space of (2) can be slightly larger than the space of the cushion layer, and the value range of alpha is suggested to be 1-3. And then, under the global coordinate system G, shifting all the particle centroids one by one to randomly generated spatial positions to obtain an initial cushion layer model, as shown in FIG. 2.

Information for all grains of the initial pad layer model can be noted as m_i，I_i，R_i，u_i(t＝0)，v_i(t＝0)，ω_i(t ═ 0) }, obviously, the particle distribution in the initial cushion model is relatively loose, the volume rate of the cushion particles is generally not achieved, mutual invasion exists among the particles, if a compact cushion model is to be obtained, the spatial positions of all the particles need to be redistributed under the action of gravity and the constraint of a boundary, the contact among all the particles needs to be searched, the contact force among the particles is calculated, then the motion equation of all the particles is solved, and the spatial distribution { u ═ of all the particles of the compact cushion model can be obtained_i(t_m)，v_i(t_m)，ω_i(t_m)}。

(3) The loose mat particles are compacted under the action of gravity and vibration to obtain a compacted mat model. Because the particle compaction process relates to the change of the spatial position of the particles, the particle models which are contacted need to be detected through a quick contact search algorithm, and the contact force between the particles can be calculated by defining a contact collision model, so that a motion equation can be established and solved for each particle model to realize the simulation of sand cushion compaction.

Inter-particle fast contact search algorithm: when the sand cushion particles are compacted, the sand particles are contacted and separated from each other every moment, and finally, all the particles form a stable contact state through extrusion and shearing action. Therefore, if the discrete element method is adopted to simulate the sand particle compaction process, all N of the bedding layer is solved firstly_cContact between individual particles determines a problem. As shown in FIG. 3, in the method for determining contact between spherical particles, the distance between the center of sphere and the sum of the radii of the two particles are generally compared, r_i、r_jRadius of particles i and j, respectively, l_ijIs the interparticle distance, the interparticle overlap distance

The calculation formula of (2) is as follows:

if it is

It indicates that particle i is in contact with particle j, if

It indicates that no contact between the particles has occurred.

However, if a global contact search algorithm is used, it is necessary to search every particle for all other N_cDetermination of contact between 1 particle, total N_c(N_c-1) touch determination.For example, for 10⁴About 10 is required for each particle to perform a global contact search⁸The secondary contact determination results in excessive calculation cost. Therefore, the invention improves the contact search algorithm and provides a new quick contact search algorithm, which comprises the following specific processes: the intersection detection is performed only between a small number of particles in the vicinity of any particle, which may be in contact, at a time, which can greatly reduce the number of contact determinations.

Aiming at the characteristics of the sand cushion layer compacting process, the invention provides a local search algorithm which is called as a space planning method. Firstly, in a global coordinate system G, equally dividing the global space where cushion layer particles are located into N_sAnd the small spaces are equivalent to rectangular boxes, so that the contact search range between each particle and other particles is reduced from the whole space to the local space, and the contact search between the particles can be performed among a small number of particles contained in each box. If each box contains too many particles, it can be further divided into smaller boxes until each box contains the appropriate number of particles, as shown in FIG. 4 a. This approach has the advantage of being able to lock all N's quickly_cThe local space to which each particle belongs greatly reduces the contact search range and the contact judgment times among the particles; a disadvantage is that a particle at any local spatial boundary may contact a particle in its neighboring local space in addition to a particle in that space, as shown in fig. 4b, thus also requiring a contact determination across neighboring spatial particles. In order to solve the problem, when the particle i in any local space is searched in contact with other particles, the range to be searched is slightly enlarged on the basis of the original space, and the new space can be called as an expanded space. The expanded space contains both the particles inside the original local space and the particles on both sides of the boundary of the space adjacent to the expanded space, so that the contact search for each particle is not missed any more. Therefore, each box corresponds to an expansion space to be searched, and N is assumed to be contained in each expansion space on average_kParticles (N)_k≈N_c/N_s) Then the number of contact determination times is reduced to N after the space planning method is adopted_c ²/N_sThen, withN to which each particle i may be exposed_kThe-1 other granule number may be recorded in the touch index array index (i).

Aiming at the condition that the falling rocks impact the sand cushion, the space planning method is not suitable for solving the contact search problem between the falling rocks and the particles because the sizes of the falling rocks particles and the cushion particles are not in the same order of magnitude. Therefore, the invention also provides a sphere space method for searching cushion layer particles which are possibly contacted with falling rocks, as shown in fig. 5, and the basic idea of the method is as follows: in the local coordinate system L of the falling rock particles (i ═ 0)₀In (1), a certain radius R is established_sThen calculating the distance L between the falling rock mass center and all other particles_0jIf L is_0j<R_sIt means that the bedding grain is within the ball domain space for which the falling rocks are to be searched and recorded in the contact Index (0) with which the falling rocks may come into contact. Generally, a sphere radius R is suggested_sThe value taking method comprises the following steps: r_s＝R₀+α·r_m (2)

Wherein R is₀Denotes the falling rock radius, r_mThe maximum radius of all cushion layer particles is shown, and the suggested value range of alpha is 1-2. The method has the advantages of simple algorithm and capability of quickly reducing the range of contact search.

According to the improved inter-particle rapid contact search algorithm, all adjacent particle contact indexes index (i) which are possibly contacted with the particle i are obtained, then accurate intersection detection is carried out according to the formula (1) in the particle range recorded by the index (i), the times of accurate detection can be reduced, and finally the number of the particle which is actually contacted with the particle i is recorded in the contact (i) array, so that a foundation is laid for further calculating the inter-particle contact force.

Interparticle contact collision model: in order to further calculate the contact force among the particles, the invention provides a contact collision model among the particles. First, the contact force F between particles_ijCan be decomposed into normal contact force F_nijAnd tangential contact force F_tijI.e. F_ij＝F_nij+F_tijThe direction of the normal contact force is parallel to the direction of the connecting line of the mass centers of the particles and is tangentialThe direction of the contact force is parallel to the common tangential direction of particles i and j at the point of contact. As shown in figure 6 of the drawings,

expressed in the local coordinate system L of the particle i_iIn which the contact force of a particle j on a particle i is indicated by the first subscript mark representing the stressed particle and the second subscript mark representing the stressed particle, the interaction of the force according to Newton's third law has

The invention is based on contact forces acting on the particles i

The contact collision model is described as an example.

(1) Normal contact force model

The normal contact force model comprises a linear elastic model, a nonlinear elastic model, an elastic-plastic model, a viscoelastic model, a viscoplastic model and the like, and different energy dissipation mechanisms are considered.

In order to simulate the mesoscopic mechanical characteristics and the energy dissipation mechanism of the rock fall impact sand cushion and analyze the influence of the cushion material strength, viscosity and density impact factors, the normal contact force adopts a linear spring damping model (LSD for short), and the normal contact force adopts a linear spring damping model

Can be divided into two parts, one part is the depth of mutual invasion with the normal direction

Spring force of interest

Another part is the velocity of invasion from the normal direction

Related viscous damping

The general form is:

first, with respect to the elastic force

Calculation formula, k_nRepresenting the Hooke spring rate, and the elastic moduli E of particles i and j_i、E_jAnd radius r_i、r_jRelated, k_nThe calculation formula is as follows,

while

The amount of deflection of the hooke spring is shown, since the present invention assumes the particle model to be a rigid sphere, the amount of spring deflection is characterized by the overlap distance between two intersecting spherical particles, as shown in fig. 6, where the overlap is exaggerated for ease of illustration,

the calculation formula is as follows,

wherein l_ijIs the inter-particle distance (r)_i+r_j-l_ij) The distance of the overlap is indicated,

is the unit direction vector with the centroid of particle j pointing to the centroid of particle i,

representing the depth vector of the intrusion of particle j into particle i.

Second, for viscous damping

C is a calculation formula of_nRepresenting the damping coefficient by calculating^[174]，

Wherein m is^*Represents equivalent mass, m_iAnd m_jThe mass of particles i and j, respectively. Xi_nIs a dimensionless coefficient, and a normal speed recovery coefficient e_nIn which e_n＝v_r/v_i，v_rIndicating the normal rebound velocity, v_iIndicating the normal incidence velocity.

While

Expressed in a local coordinate system L_iNext, the normal relative movement speed of particle j intruding into particle i at the contact point:

wherein the content of the first and second substances,

is the velocity of the movement of particle j relative to particle i,

is the absolute velocity of the particle j,

is the absolute velocity of particle i.

Normal contact force obtained from the LSD model

And interpenetration depth vector

The relationship of (2) is shown in FIG. 7. In the figure, the dotted straight line is AND

Elastic contact force in linear relation

In superposition of viscous damping

The normal contact force shown by the solid line is obtained

The variation of the contact force can generally be divided into two phases, according to the normal speed increment: when in use

When the particles move oppositely, the compression stage is adopted; when in use

When the particles move in opposite directions, the recovery phase is achieved. At the start of the compression phase,

but because of the initial impact velocity

Therefore, there is a viscous damping force as shown in equation 5

At the end of the recovery phase, if the normal contact force continues to decrease, this results in a contact force

Less than 0 occurs, however, in the sand bedding model of the present invention it is assumed that there is no gravitational interaction between the sand particles, and therefore, when

When it is forced to make

At this time, the depth of penetration between particles

This can be seen as residual deformation of the particles. The LSD model has the advantages of simple model and high calculation efficiency, and can consider the influence of elastic force and viscous damping.

(2) Tangential contact force model

The tangential contact force model of the particles generally adopts a molar-coulomb friction criterion, and the calculation formula of the tangential contact force is as follows:

wherein the content of the first and second substances,

indicating the magnitude of the tangential contact force or friction, mu_sDenotes the maximum coefficient of static friction, μ_dRepresenting the coefficient of dynamic friction, the maximum static friction coefficient generally being slightly greater than the coefficient of dynamic friction, i.e. mu_s>μ_d。

The tangential velocity of particle j relative to particle i at the point of contact can be generally decomposed into the relative tangential velocity of the particle's centroid and the relative linear velocity of the point of contact caused by the rotation of the particle about the centroid, as shown in fig. 8. The calculation formula of the tangential relative movement speed is as follows:

wherein the content of the first and second substances,

and

respectively represent local coordinate coefficients L_iAnd L_jThe angular velocities of the particles i and j below each around their centroid,

and

respectively representing the vector radii of the centroids of particle i and particle j to the point of contact,

is the tangential velocity of particle i relative to particle j at the point of contact.

Suppose the relative tangential velocities of particle i and particle j initially

Then the sliding gradually along the tangential direction from the relative static state under the action of external force, but the static frictionForce of

Is uncertain and, usually in a static state,

is balanced and equal to the external force until the external force increases and exceeds the maximum static friction force

Then, the kinetic friction force is converted into the kinetic friction force

Because of the coefficient of dynamic friction mu_dIs a constant related to the particle material, so the kinetic friction force is easier to calculate, and the magnitude of the static friction force can be calculated by definition

In relation to the relative speed of movement

To approximately solve the problem,

wherein the content of the first and second substances,>>1 and ζ>0, as shown in FIG. 9, frictional force

Always in the direction of the tangential velocity relative to the contact point

In the opposite direction. The static friction force increases with the speed from zero and rises to the maximum static friction force

Then decreases again until gradually approachingDynamic friction force

The parameters and ζ control the slope of the rising and falling segments. The tangential contact force model has the advantage of facilitating the numerical solution of the tangential contact force by an explicit time step method.

As shown in fig. 6, the moment that particle i is subjected to particle j is equal to the cross product of the tangential contact force and the vector radius,

since any particle i may come into contact with a plurality of particles, the resultant force and resultant moment of all contact forces to which the particle i is subjected,

wherein the content of the first and second substances,

indicating the contact force of the jth particle on the ith particle, N_iIndicates the number of particles coming into contact with the ith particle,

representing the resultant of all forces to which particle i is subjected.

The contact moment of the jth particle on the ith particle is shown,

representing the resultant moment of all forces to which particle i is subjected.

Particle compaction algorithm: the inter-particle contact search algorithm solves the problem of mutual intrusion judgment among all particles, and the inter-particle contact collision model solves the problem of inter-particle contact force calculation. The essence of the process of compacting the cushion particles is that under the action of external force, the particles are mutually engaged to form a stable structure through the shearing action among the particles, so that the redistribution and the compaction of the spatial positions of the particles are realized. In order to improve the randomness of the distribution of the cushion layer particles after being densified and the densification efficiency, the invention provides a particle densification algorithm under the action of gravity and vibration.

Firstly, generating an initial cushion layer model { u } in loose distribution according to a particle generation algorithm and a space delivery algorithm_i(t＝0)，v_i(t＝0)，ω_i(t ═ 0) }, as shown in fig. 10(a), the particles of the initial cushion model are randomly distributed in the designated space, and there may have been mutual invasion between the particles, and assuming that the spherical particles are rigid bodies, the contact collision model according to the present invention defines the mechanical parameters of the particles, including the contact stiffness k_nDamping coefficient c_nMaximum coefficient of static friction mu_sCoefficient of kinetic friction mu_d(ii) a Then, the gravity g (t) is applied to all the particles, and at the same time, the vibration loads u (t) in the X-axis direction and the Y-axis direction are simultaneously applied to the lateral rigid boundary of the model, and the vibration loads can be simple harmonics U (t) u (t) controlled by displacement₀sin (t) or other form of loading; finally, according to the gravity G (t) and the vibration load U (t) borne by the particle system and the contact force between the particles, establishing a motion equation of the particle system and solving, simulating the redistribution process of the spatial position and the motion state of all the particles under the action of external force, stopping vibrating after the volume rate of the cushion particles reaches the target volume rate, enabling the cushion to be dense continuously under the action of gravity until all the particles stop moving and the spatial position does not change any more, wherein the cushion model obtained by the gravity and vibration compaction method is very close to the particle distribution and the stress state of the actual cushion, and the dense cushion model is marked as { u } after being dense_i(t_m)，v_i(t_m)，ω_i(t_m) As shown in fig. 10 (b).

The motion equation and the solving method of the particle system are as follows: root of herbaceous plantAccording to the contact search algorithm and the contact collision model, the resultant force of the contact force applied to each particle i is calculated, and in order to further calculate the spatial position change and redistribution of all the particles under the external force compaction action (the external force may include gravity, vibration load and the like, and can be uniformly recorded as F_ui(t)), equations of motion need to be established and solved for the particle system.

The solution of the motion equation of the particle system needs to disperse time into small time steps, then substitutes the initial conditions of the particles, and gradually iterates the solution of the motion equation.

It is assumed that the motion state of each particle i can pass through the spatial position

Translation speed

And angular velocity

Described, as shown in fig. 11, the basic idea of the loop calculation for the particle system is: firstly inputting all N of the initial cushion layer models_cPhysical property, spatial position and motion state information { m ] of each particle at time t-0_i，I_i，R_i，u_i(t＝0)，v_i(t＝0)，ω_i(t-0) and then into the first cycle (time cycle), analysing the state of motion { u } of all particles at time t_i(t)，v_i(t)，ω_i(t) }; then, the method enters a second cycle (particle cycle), the motion state of the particles i at the time t is known, all particle indexes contact (i) in contact with the particles i are established according to a contact search algorithm, and the contact force resultant force on the particles i is calculated according to a contact collision model

Sum and resultant moment

And gravity and boundary applied external force effects F_uiThen, according to the motion equation, the change { delta u) of the motion state of the particle i in the delta t time step is calculated_i(t)，Δv_i(t)，Δω_i(t) } (Δ t is the time step of the explicit solver), until all N_cThe analysis of each particle is finished; finally, the motion states { u ] of all the particles at the t + delta t moment are updated simultaneously according to the motion state change in the delta t time step_i(t+Δt)，v_i(t+Δt)，ω_i(t + Δ t) }, and enter the next time step cycle to repeat the above calculation until the compaction process reaches the termination time t_mOutput t_mSpatial position and motion state of all particles at a time { u }_i(t_m)，v_i(t_m)，ω_i(t_m) And fourthly, obtaining the compact cushion layer model, and introducing the specific method as follows.

(1) Equation of motion

According to Newton's second law of motion, the equation of motion for each particle in the local coordinate system is

Wherein m is_iWhich is indicative of the mass of the particles,

representing the translational acceleration of the center of mass of the particle,

indicating that particle i is subjected to the combined forces of all interparticle contact forces,

indicating the external force action such as gravity, vibration load, etc.

Representing the rotational acceleration of the particle i around the center of mass,

representing the resultant moment of particle i subject to all interparticle contact forces,

representing the action of external moment such as vibration load.

The moment of inertia of a particle is expressed, and the calculation formula for the magnitude of the moment of inertia of a spherical particle is:

(2) solving method

In discrete element simulation, the motion process of each particle in each time step delta t needs to be solved through a time step method, and the delta t of the cushion layer particles_critDepending on the minimum value of the critical time step for the contact pairs between all particles, the formula is calculated as:

Δt_critwith mass m of the particles_iAnd contact stiffness k_ijIn this connection, it can be seen from the formula that the smaller the particle mass, the greater the contact stiffness, Δ t_critThe smaller the theoretical calculation of (c). However, Δ t_critThe theoretical calculation result is only used as a reference for the value of delta t, and the specific value size can be determined only by analyzing the numerical calculation convergence.

After the time step Δ t is determined, the motion equation of each particle i can be solved by using a time step method, which includes the following detailed steps:

step 1: given cushion particlesi initial conditions. These conditions include the mass m of the particles_iMoment of inertia about the centre of mass

Initial centroid position

Initial centroid translation velocity

Initial centroid rotation angle

Initial centroid rotational velocity

Step 2: establish the index (i) of the particles with which particle i may come into contact. According to a rapid contact search algorithm, searching particles which are possibly contacted with the particles i every n multiplied by delta t time steps, and establishing a particle index (i) to facilitate the next step of accurately judging the contact of the particles in a small range. Value of n and the planning space size L_cellAnd spherical domain space radius R_sAnd the distance traveled by the particle at in the time step

When the spatial position of the particle is changed greatly and exceeds the original local search space range, local contact search and index establishment need to be carried out again. Repeated trial calculation shows that in the cushion granule compacting calculation example of the invention, when the time step length delta t is 0.1ms, the suggested value range of n is 100-1000, and in the falling rock impact cushion calculation example, the suggested value range of n is 1000-10000. Then, a few particles in the range of the particle index (i) are accurately intersected and detected every deltat time step, and a contact index contact (i) which is contacted with the particle i is obtained.

And step 3: calculating the sum of resultant forces borne by the particles i at the moment tAnd (4) resultant moment. According to the contact collision model, firstly, the contact force of the particles is calculated according to the formulas (17) and (18)

And contact torque

Then calculating the contact force resultant of the particle i

Sum and resultant moment

A general expression of a calculation formula can be written as,

and 4, step 4: calculating the grain i centroid translation increment within the time increment Δ t

Increment of translational velocity

And increment of rotational speed

And 5: in a local coordinate system L_iThen, for all N_cThe solving processes from step 2 to step 4 are repeated for each particle at the same time, and then the position coordinates of all the particles at the t + delta t moment are updated at the same time

Translation speed

And speed of rotation

Step 6: and (4) judging whether the movement reaches a termination condition, for example, the cushion granules gradually tend to be stable and stop moving under the gravity compaction, or the falling rocks impact the cushion granules.

And 7: and selectively outputting the motion result of each time step according to the simulation purpose. For example, the spatial position of the rockfall particles

Translation speed

And speed of rotation

And contact force

And contact torque

And the like.

Examples

According to the discrete element analysis model, the invention is used for carrying out numerical simulation on the spherical falling rock impact sand cushion test and verifying the accuracy and reliability of the analysis model by comparing with the test result.

Ball-shaped rockfall impact cushion test

The falling rock impact coarse soil (Granular soil) cushion impact effect under different working conditions has four conditions, namely a falling rock inclined impact loose cushion, a falling rock positive impact loose cushion, a falling rock repeated impact compact cushion and a falling rock impact lightweight ceramic aggregate (LECA), and the embodiment of the invention selects a verification test of macroscopic and microscopic analysis models of the falling rock positive impact loose cushion.

The falling rock impact test device consists of a lifter tower, a coarse-grained soil cushion layer and spherical falling rocks. As shown in fig. 12, coarse-grained soil is filled into a Reinforced Concrete (RC) annular foundation pit under the action of gravity to form a cushion layer, and the thickness of an RC support at the bottom of the foundation pit is 0.5m, so that the deformation of the pit wall and the support can be ignored and can be regarded as a rigid boundary. The obtained sand cushion layer has the diameter of 10.70m, the thickness of 2.00m and the relative density D_rThe friction angle is 30 degrees, the grading curve of the particles is shown in figure 13, and the particle size distribution range is 0.08-34.29 mm. During the test, the upper layer (about 1m thick part) of the cushion layer is removed after each impact and then backfilled to ensure that the compactness of the cushion layer is unchanged.

In the test, the radius R of the reinforced concrete falling rock is 0.45m, the mass m of the reinforced concrete falling rock is 850kg, and the triaxial accelerometer is embedded in the position of the center of the sphere and used for recording the acceleration time-course curve of the falling rock and multiplying the acceleration time-course curve by the falling rock mass m to obtain the impact force time-course curve of the falling rock acting on the cushion layer. Testing the rockfall speed and penetration time course curve by integrating acceleration with timeAnd (4) obtaining. The falling height range of the working condition design of the falling rock positive impact loose cushion layer is 5.00-18.45 m, and the corresponding initial impact velocity v₀Is 9.90 to 19.02 m/s.

Influence of the cushion layer lateral boundary conditions on the test results in the test: because the lateral pit wall of the cushion layer plays a role in rigid restraint on the cushion layer, when falling rocks impact the cushion layer, the compression wave can be transmitted to the rigid lateral boundary, then a reflected wave is formed and then is transmitted to the position of an impact point, so that the principle of judging whether the lateral boundary influences the falling rocks impact force is that the transmission round-trip time t of the compression wave is observed_cWhether the total impact action time t is longer than the total impact action time t of falling rocks_m. If t is_c>t_mIf the side boundary does not influence the impact force, the action of the impact force is finished before the reflected wave returns; if t is_c≤t_mIt is then demonstrated that the compression waves reflected by the lateral boundary have an effect on the impact force. In the test, according to the falling rock impact force time course curve and the time difference of the starting point of the bottom plate compressive stress time course curve, the one-way propagation time of the compression wave in a cushion layer with the thickness of 2m is obtained to be 0.015s, so that the propagation speed of the compression wave in the cushion layer is measured to be v_c130 m/s. Distance L of impact point from side boundary in test_c5.35m, so the round trip time t for the compression wave to propagate to the side boundary_c＝2L_c/v_c0.082s, and the average value of actually measured impact time of falling rocks is about t_m0.090s, although t_c<t_mHowever, since the influence time is close to the end time of the rock fall impact action, the influence of the lateral boundary under the test condition can still be ignored, so that the influence of the lateral boundary condition on the impact result is also eliminated when a macroscopic and microscopic model of the cushion layer is established later.

Experimental verification of mesoscopic discrete element model

The discrete element model of the sandy soil bedding course is established by a discrete element method, the problem of selecting the size of the bedding course model and the grain diameter of the bedding course grains is firstly solved, the selection principle is to reproduce the size of the model and the grain diameter grading in the test as much as possible, but the overlarge size of the model and the undersize grain diameter can generate the huge number of discrete element grains, so the calculation cost is exponentially increased, and therefore, on the premise of not influencing the overall mechanical property and the dynamic response of the bedding course, the proper size of the model and the grain diameter are usually selected, and the convergence analysis of the grain diameter is carried out.

(1) Geometric dimension of cushion layer model

For the discrete element model, the number of particles required to be established by the full-size cushion model under the simulation test condition is too large, so that the invention applies the non-reflection boundary condition on the boundary of the cushion model with the limited size to eliminate the influence of the stress wave reflected by the boundary. According to the boundary processing method, after the stress wave propagates to the artificial boundary, the boundary should completely absorb the energy of the incident wave. The condition for the medium not to generate reflected waves at the boundary is that the velocity v of the mass point propagating to the medium boundary is set_iEqual to 0, the discrete element model is processed by making the velocity of the particles in contact with the boundary in the contact collision model restore the coefficient e_nAnd (0) to make the particle speed contact with the boundary and then attenuate to zero, thereby achieving the purpose of absorbing the incident wave energy.

According to the research of the lateral boundary constraint influence of the granular medium, the ratio L of the size of the cushion model to the rockfall radius is shown_cWhen the/R is more than or equal to 5.0, the influence of the lateral boundary constraint on low-speed impact disappears rapidly, so that the invention takes the lateral dimension of the discrete element cushion layer model as L_c＝2.4m，L_c/R＝5.3>5.0, and then applying a non-reflective boundary condition on the lateral boundary. The thickness h of the bedding model was taken to be equal to the thickness 2m of the test, and a non-reflecting boundary condition was also imposed on the bottom boundary in order to eliminate the influence of the bottom boundary condition on the falling rock impact.

(2) Particle size of cushion layer particles

According to the grading curve of coarse soil particles measured by tests, the number of the estimated bedding particles reaches tens of millions, and the calculation is difficult in discrete element numerical simulation due to too much time consumption. However, current studies on the size effect of the particulate material show that when the ratio D/D of the specimen size D to the particle diameter D of the particles is sufficiently large, the static and dynamic mechanical properties of the mat are substantially no longer influenced by D/D. For example, through the research of biaxial compression test of two-dimensional particle media, the influence of particle size effect on the shear modulus of the media is reduced to the minimum when the D/D is increased to 30-40. Triaxial compression tests simulating a non-sticky granular material by the PFC3D program show that the D/D is 10, 20, 30, 40 and 50 respectively, the obtained peak intensity has convergence, and the D/D has little influence on the result after exceeding 50. The mechanical response of dynamic loads is independent of particle size, as long as the particle size is sufficiently small.

The invention selects a spherical rockfall impact soil cushion test, and the grain size range measured by the test is 0.08-34.29 mm. In the prior art, two kinds of spherical particles with the variation ranges of 100-300 mm and 50-150 mm are adopted to simulate the particles, the particle sizes of the particles are uniformly distributed in the range, the number of the particles is respectively about 10000 and 80000, numerical simulation results show that model impact force calculation results of the two kinds of particles are similar, and the accuracy of simulation results is verified by comparing with test results, which shows that although the particle size range of a discrete element model exceeds the particle size range measured by tests, the particle size range is small enough, the size effect of the particle size can be ignored, the significance of continuously reducing the particle size is not large, and finally, a particle model with the particle size distribution range of 100-300 mm with a small number is selected. In addition, the distribution form of the particle size of the discrete element particles is different from the grading curve of the particles in the test, but the analysis of the microscopic mechanism of the discrete element calculation result is not influenced. Therefore, the calculation accuracy and efficiency of discrete element analysis are comprehensively considered, spherical particles with the particle size range of 33-99 mm are selected by the cushion layer model, about 12000 spherical particles are selected in total, and the distribution form of the particle sizes is uniform.

(3) Establishment of rock fall impact cushion layer model

Based on the above test conditions, the dimensions of the discrete shim model, as well as the particle size and number, were determined. Next, a rock fall impact cushion model is generated according to the space putting and compacting algorithm proposed by the present invention, as shown in fig. 14. Firstly, randomly putting 12000 generated spherical particles in a space range of a cushion layer model; then, loosely distributed cushion particles are arranged inCompacting under the action of gravity and vibration to obtain a cushion layer with the thickness of 2.0m and the particle volume rate of about 50 percent; finally, spherical rockfall particles are generated at the central position of the surface of the compacted cushion layer, and the rockfall particles have the diameter of 0.90m and the density of 2227kg/m according to the test conditions³The mass was 850 kg.

(4) Calculation parameter determination and comparison analysis with experiment

After the discrete element model of the rockfall impact cushion is established, parameters of the discrete element model are determined according to test conditions, and the accuracy of a model prediction result is verified. Firstly, the test result with the impact speed of 14.00m/s in the test is selected to calibrate the parameters of the DEM model, and the calibration method is as follows.

First, for the physical parameters of the falling rock model, since the falling rock is assumed to be a rigid body in the discrete element model, the density, elastic modulus, poisson's ratio and static and dynamic friction coefficient with the cushion particles of the falling rock spherical particles are shown in table 5-3 with reference to the test conditions and the values of the finite element model.

The density of the individual sand particles can then be determined first for the sand mat particles and the inter-particle physical parameters. Knowing the mat density from the test conditions ρ_g＝1500kg/m³And a porosity of about 50%, from which condition the density of the mat particles can be determined to be p_i＝2200kg/m³。

Finally, according to the existing DEM model parameter calibration research, the falling rock impact force and the impact action time are mainly influenced by the Young modulus E of the particles_iDensity rho_iAnd normal coefficient of restitution e_n(ii) an effect; the interparticle friction coefficient mu has a small influence on the falling rock impact force, and has a large influence on the impact force on the bottom plate. Therefore, from the test results of the falling rock impact force time course curve (mainly the peak value and the peak reaching time of the falling rock impact force), the Young's modulus E of the cushion granules was determined_iAnd normal coefficient of restitution e_nAs shown in tables 5-3. And the interparticle friction coefficient mu is mainly related to the macroscopic expression friction angle psi of the cushion layer_gIn connection with the method for taking the value of mu, the friction of DEM cushion layers for DEM granular media with given boundary conditions and porosityThe angle and the angle of friction between spherical particles satisfy the following relationship,

tan(ψ_g)＝α+βtan(ψ_p) (31)

wherein psi_pDenotes an inter-particle friction angle, μ ═ tan (ψ)_p) And alpha and beta are empirical coefficients related to porosity, according to the experimental value psi_gThe value range of the friction coefficient mu between particles calculated at 30 degrees is 0.25-0.35, the specific value needs to be further subjected to parameter sensitivity analysis, and the numerical research result shows that the influence of mu on the falling rock impact force within the range of 0-0.5 is very small, so the value of the invention is shown in the table.

The comparison between the numerical simulation result of the falling rock impact force time-course curve obtained from the above-described determined model parameters and the test result is shown in fig. 15 (a). The peak value of the rock falling impact force obtained by the test is about 216kN, the DEM simulation result is about 226kN, and because the bottom of the DEM model is provided with the non-reflection boundary, the simulation result does not have a second peak value like the test, so that the total impact action time is about 120 ms. And (3) continuing to simulate the impact test of the falling rocks at other impact speeds according to the set of calibration parameters, wherein the result is shown in FIG. 15(b), and the result shows that the falling rock impact force time-course curve obtained by the DEM model is still very close to the test result, so that a foundation is laid for analyzing the collision mechanism between the falling rocks and the cushion layer by using a DEM method in the next step.

TABLE 5-3 calibration parameters of DEM model for rockfall impact pad

Mesoscopic discrete element calculation results and analysis

Because the macroscopic finite element model is difficult to explain the mechanism of viscous force and dynamic resistance and plays a role in total resistance, the microscopic mechanism of the rock fall impact cushion layer is further analyzed according to the microscopic discrete element calculation result.

Method for calculating viscous force and dynamic resistance in mesoscopic model

Because the sand cushion belongs to a granular medium and has solid-liquid two-phase property, the falling rocks impact the cushion in a process similar to the flow of liquid around the surface of an object, and the falling rocks are subjected to resistance similar to the flow resistance of the object in flowing viscous fluid. The flow resistance is divided into viscous resistance and differential pressure resistance: the viscous resistance is the resistance generated by the shearing force between the surface of the object and the fluid, is related to the first power of the relative speed, and the direction of the viscous resistance is parallel to the tangential direction of the surface of the object; the pressure difference resistance is the resistance generated by the pressure stress between the surface of the object and the fluid in the vertical direction, is related to the quadratic power of the relative speed, and the direction of the resistance is consistent with the normal direction of the surface of the object. The magnitude of the two resistances is determined by the shape of the object, the surface roughness and the adhesion of the cushion material. Most importantly, the resistance to a streamlined body is mainly viscous shear force (viscous force) on the surface of the body, while the resistance to a non-streamlined body (blunt body) is mainly the resultant force (dynamic resistance) of surface pressure stress related to the shape of the body, and the pressure difference resistance of the blunt body is often much larger than the friction resistance. In general, the naturally occurring falling rocks are rarely streamlined, most falling rocks are blunt bodies and comprise spherical falling rocks, flat falling rocks and the like, so that the falling rocks are subjected to resistance mainly from dynamic resistance and have relatively small influence on viscous force when impacting a cushion layer. The theoretical assumption accords with the stress deformation characteristic of a sand cushion, the cushion in the central area of the spherical rockfall contact surface is mainly in a three-way compression state, the cushion near the edge of the contact surface is mainly in a shear deformation state, and in order to verify the theoretical assumption, the compressive stress and the shear stress on the cushion and the rockfall contact surface are compared and analyzed from a microscopic scale according to the calculation result of a discrete element model.

According to the discrete element model, the stress analysis method on the contact surface of the cushion layer comprises the following steps: providing cushion material on contact surfaceIn a layer of imaginary units, the unit stress sigma can be decomposed into compressive stress sigma along the normal and tangential directions of the contact surface_nAnd shear stress sigma_tAs shown in fig. 16 (a). According to the DEM calculation principle, the stress sigma on the rockfall contact surface is equal to the contact force F of all particles in the unit_0jThe ratio of the resultant force to the contact area dS, as shown in FIG. 16(b), the compressive stress σ_nAnd shear stress sigma_tThe formula for calculating (a) is as follows,

dS＝Rdθ(2πR sinθ) (34)

wherein N is_dSIndicates the number of particles in the cell that are in contact, F_n0jAnd F_t0jThe normal contact force and the tangential contact force of the falling rocks and the particles in the unit are respectively shown, and dS represents the unit stress area at the position which forms an included angle theta with the negative direction of the z axis. According to the theoretical assumption, the dynamic resistance to the falling rocks should be equal to the compressive stress sigma_nThe viscous force on the falling rocks, integrated over the contact surface S, should be equal to the shear stress sigma_tIntegral on the contact surface S, then according to DEM model calculation principle, the dynamic resistance F can be obtained_nAnd viscous force F_tThe formula for the calculation of (a) is,

wherein N is_cIndicating the number of bedding particles that come into contact with the falling rocks. Resistance F to falling rocks_RShould equal the resultant of the contact forces for all particles coming into contact with the falling rocks:

calculation analysis of viscous force and dynamic resistance

As shown in FIG. 17, the spherical contact surface is divided into three areas (Zonei, Zone II, Zonei II) equally, the stress state of the cushion unit at the center of each area is taken as a representative, the comparison result of the compressive stress and the shear stress is shown in FIG. 18(a-c), and the shear stress peak value sigma of the sand cushion unit of Zonei is shown in FIG. 18_tAbout the peak value σ of compressive stress_n0.07 times of the ZoneI cushion unit, proves that the Zonei cushion unit mainly receives three-dimensional compression, has small shearing deformation, and the rockfall surface mainly receives the action of pressure difference resistance (dynamic resistance) and has small friction resistance (viscous force); zoneii's sand cushion unit's shear stress peak value sigma_tAbout compressive stress σ_nThe peak value of 0.93 times proves that the falling rocks are subjected to larger viscous shearing force at the edge of the contact surface, and the magnitude of the viscous force increases with the increase of theta.

And then, further quantitatively analyzing the stress distribution characteristics of different positions on the rockfall contact surface. As shown in fig. 18(d), regions I, II and III correspond to compressive stress peaks of 0.48MPa, 0.33MPa and 0.09MPa, respectively, and shear stress peaks of 0.03MPa, 0.08MPa and 0.08MPa, respectively, and the results indicate that at the same time, the compressive stress does decrease with increasing theta, and the shear stress increases with increasing theta, which is generally much greater than the shear stress.

Finally, the effect of the dynamic resistance and the viscous force, respectively, in the total resistance is comparatively analyzed. Knowing that the dynamic and viscous forces are equal to the integral of the compressive and shear stresses, respectively, on the contact surface, fig. 19 compares the dynamic resistance F_nViscous force F_tAnd the magnitude of the total resistance F, the results show the dynamic resistance F_nThe peak value of (A) is about 85% of the total resistance F, and the viscous force F_tThe peak value of (a) is about 15% of the total resistance F, which proves that under the condition of falling rock impact, the dynamic resistance plays a dominant role in the resistance of falling rocks, and the influence of viscous force is small.

Finally, the adhesion of the mat is mainly due to the different adhesion of the mat materials, which may affect the adhesion strengthThe inter-shearing effect is related, and the viscosity of the cushion layer can be characterized by the friction coefficient among particles (making mu equal to mu)_d＝μ_s). Fig. 20 compares the changes of the interparticle friction coefficient within the range of 0-0.5, and the rockfall impact force and the invasion depth time course curve, and the result shows that the calculation result is not strong in the sensitivity to the friction coefficient for the spherical rockfall impact cushion problem, and the result proves that the viscous force plays a very small role in rockfall resistance and can be ignored again.

The invention can simulate the sandy soil cushion layer consisting of discrete particles, determines the stress deformation characteristic of the sandy soil cushion layer through the extrusion deformation of the particles and the contact force among the particles, and can analyze the relation between the compressive stress on the contact surface and the strength, viscosity and density of the cushion layer on a microscopic level by introducing a normal and tangential contact force calculation model into the contact model among the particles.

The sand cushion mesoscopic model with controllable density is generated by regarding the cushion as a mesoscopic discrete system consisting of three-dimensional spherical particles with 6 degrees of freedom, providing a random putting of falling rock particles and a gravity and vibration compaction algorithm; the global contact search algorithm among the particles is improved, a rapid contact search algorithm combining a planning space method and a spherical domain space method is provided, and the contact judgment efficiency is greatly improved; an inter-particle normal and tangential contact force model based on Linear Spring Damping (LSD) is introduced, and the resistance on a contact surface and the extrusion and shearing action among particles are analyzed in detail on a microscopic level. The microscopic discrete element calculation result shows that the impact near region cushion layer unit compression stress is far greater than the shear stress, the dynamic resistance plays a leading role in the resistance of falling rocks, and the influence of viscous force is small.

The present invention is not limited to the above embodiments, and any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention are included in the scope of the claims of the present invention which are filed as the application.

Claims (10)

1. A three-dimensional microscopic discrete element analysis method for a rockfall impact cushion is characterized by comprising the following steps:

the sand particles in the cushion layer are assumed to be spherical particles;

putting a plurality of spherical particles one by one in a designated space to generate an initial cushion layer model;

compacting the initial cushion layer model under the conditions of gravity, vibration and boundary constraint;

solving the motion equations of all spherical particles in the compacting process;

judging whether compaction is finished according to the solution result of the motion equation to obtain a compact cushion layer model;

and carrying out a falling rock impact cushion test on the compact cushion model to obtain a three-dimensional microscopic discrete element analysis result of the falling rock impact cushion.

2. The method for analyzing three-dimensional microscopic discrete elements of a rockfall impact pad according to claim 1, wherein the designated space is a cylindrical space; the boundary of the designated space is a rigid boundary; the number and the particle size distribution of the spherical particles are obtained through probability density distribution function assignment.

3. The method for three-dimensional microscopic discrete element analysis of a rockfall impact pad according to claim 1, wherein a plurality of spherical particles in a designated space constitute a particle system; the state change of the particle system is respectively described by a global coordinate system and a local coordinate system; describing contact search and collision among all spherical particles in the particle system through a global coordinate system; the translation and rotation of each spherical particle in the particle system is described by a local coordinate system.

4. The method for analyzing the three-dimensional microscopic discrete elements of the rockfall impact pad according to claim 3, wherein the method for establishing the motion equation comprises:

detecting contacted spherical particles by a quick contact search algorithm;

calculating a contact force between the spherical particles which are contacted according to the contact collision model;

and establishing a corresponding motion equation according to the resultant force and resultant moment of the contact force applied to each spherical particle.

5. The method of claim 4, wherein the equation of motion of each spherical particle in the local coordinate system is

Wherein m is_iWhich is indicative of the mass of the particles,

representing the translational acceleration of the center of mass of the particle,

indicating that particle i is subjected to the combined forces of all interparticle contact forces,

representing the external force action of gravity, vibration load and the like,

representing the rotational acceleration of the particle i around the center of mass,

representing the resultant moment of particle i subject to all interparticle contact forces,

showing the action of external moment such as vibration load,

representing the moment of inertia of the particle.

6. The method for analyzing the three-dimensional microscopic discrete elements of the rockfall impact pad according to claim 5, wherein the formula for calculating the magnitude of the moment of inertia of the spherical particles is:

wherein the content of the first and second substances,

denotes the moment of inertia of the particle, m_iMass of spherical particles, r_iThe radius of a spherical particle.

7. The method for analyzing the three-dimensional microscopic discrete elements of the rockfall impact pad according to claim 1, wherein the solution process of the motion equation comprises:

obtaining initial conditions of cushion layer particles, wherein the initial conditions comprise the radius, the mass, the rotational inertia, the centroid position, the translation speed, the rotation angle and the rotation speed of the particles;

searching particles which are contacted with each other in each time step according to a quick contact search algorithm;

calculating contact force resultant force and resultant moment between each particle and particles in mutual contact according to the initial conditions;

calculating the translation increment, the translation speed increment and the rotation speed increment of the mass center of each particle in the time step;

under the local coordinate system, the time step is continuously updated, and the steps are repeated for all the particles until the movement reaches the termination condition.

8. The method for analyzing the three-dimensional microscopic discrete elements of the rockfall impact pad according to claim 7, wherein the calculation formula of the centroid translation increment is as follows:

the calculation formula of the translational velocity increment is as follows:

the calculation formula of the rotation speed increment is as follows:

wherein the content of the first and second substances,

in order to translate the increment of the center of mass,

representing the translational acceleration of the centroid of the particle, at is the time step,

as local coordinate coefficient L at time t_iThe lower particle i is subjected to the resultant of all contact forces,

m represents the external force action such as gravity, vibration load and the like at the moment t_iIs the mass of the spherical particles,

representing the resultant moment of contact force experienced by particle i at time t,

showing the action of external moment such as vibration load at the moment t,

representing the moment of inertia of the particle.

9. The three-dimensional microscopic discrete element analysis method of the rockfall impact pad layer according to claim 1, wherein the contact search method of rockfall particles and spherical particles in the compact pad layer during rockfall impact pad layer comprises:

establishing a spherical domain space with a certain radius in a local coordinate system of the rockfall particles;

calculating the distances between the falling rock mass center and all the particles in the compact cushion layer;

judging whether the particles are in the spherical space or not according to the distance;

and carrying out contact search on the particles in the spherical domain space and the falling rocks, and judging whether the falling rocks are contacted with the spherical particles.

10. The method of claim 9, wherein the radius of the sphere domain space is:

R_s＝R₀+α·r_m (2)

wherein R is₀Denotes the falling rock radius, r_mThe maximum radius of all cushion layer particles is shown, and the suggested value range of alpha is 1-2.

Images