Abstract
In the rockfall prevention and control project, the reinforced concrete (RC) slab and sand (gravel soil) soil cushion layer are commonly used to form the protection structure, thereby resisting the rockfall impact. Considering that the oversized deformation of the cushion layer under impact load using the finite element simulation cannot converge, this article establishes a numerical calculation model using smoothed particle hydrodynamics–finite-element method coupling (SPH–FEM). First, the standard Lagrange finite-element mesh is established for the whole model using ABAQUS, and then the finite-element mesh of the soil cushion layer is converted to SPH particle at the initial moment of the calculation, and finally the calculation results are solved and outputted. The results indicate that, compared with the results of the outdoor rockfall impact test, the relative errors of the rockfall impact force and the displacement of the RC slab are within 10%, which proves the rationality of the coupling algorithm; moreover, in terms of the numerical simulation, the SPH–FEM coupling algorithm is more practical than the finite element for reproducing the mobility of the rockfall impacting the sand and soil particles. In addition, at an impact speed of less than 12 m·s−1, the cushion layer is able to absorb more than 85% of the impact energy, which effectively ensures that the RC slab is in an elastic working state under small impact energy and does not undergo destructive damage under large impact energy; the peak impact force of the rockfall is approximately linear with the velocity, and the simulated value of the peak impact force is basically the same as that of the theoretical value of Hertz theory; the numerical simulation is good for reproducing the damage process of the RC slab in accordance with the actual situation. The SPH–FEM coupling algorithm is more justified than the FEM in simulating the large deformation problem, and it can provide a new calculation method for the design and calculation of the rockfall protection structure.
1 Introduction
Rockfall is one of the common geological hazards in high mountain and canyon areas, which poses a huge threat to life and property in the affected area due to its high energy and abrupt dynamic characteristics [1,2,3]. Generally, rockfall prevention is divided into active and passive protection, but most hazardous rock formations develop on steep slopes, making active protection very difficult. Reinforced concrete (RC) shed tunnels and pile plate retaining walls are widely used passive protection structures due to their high stiffness and good protective performance. Both structures consist of a composite cushion layer in contact with the rockfalls and RC slabs. The cushion layer of the shed tunnel structure is placed on top of the RC slab, while the cushion layer of the pile plate retaining wall is placed on the front side of the concrete slab. Because the cushion layer is in direct contact with the rockfalls, it needs to be composed of materials that have excellent energy dissipation performance, such as sand and gravel soil particles and flexible EPS and other energy-consuming substances. In the actual project, how to carry out suitable protection structure design and calculation is the focus and difficulty of the current scholars [4].
The existing research on the dynamic response of protection structures under rockfall impact loads mainly includes theoretical, experimental, and numerical methods. Lei et al. [5] studied the shear resistance of concrete slabs under rockfall loads based on theoretical methods. Wang et al. [6] investigated the dynamic response of the rockfall impacting a concrete slab using Olsson’s fluctuation control equation and verified it with numerical simulations. Wang et al. [7] and Su et al. [8] analyzed the damage modes of RC slabs under low-velocity impacts. However, these studies only explored the impact resistance of a single RC structure and did not consider the comprehensive dynamic response of structures equipped with an energy dissipation layer (i.e., cushion layer). In addition, most of the theoretical methods currently used for analysis and calculation are based on Hertz contact theory, which cannot consider the elastic–plastic deformation of the structure and have limitations. In terms of physical model experiments, Mougin et al. [9] established a large-scale RC slab model and studied the dynamic response of the slab under 3,360 kJ impact energy. Kishi et al. [10] conducted experimental research on the ultimate impact resistance of prestressed RC protection structures. Bhatti et al. [11] carried out an experimental study on the impact effect of shed slabs under different impact energies. Zhang et al. [12] conducted multi-working condition rockfall impact tests based on a large-scale rockfall impact platform and systematically analyzed the change in impact force and deformation mode of RC slabs. Based on 14 square RC slabs, Gu et al. [13] summarized the deformation damage patterns of the slabs under different impact conditions. Zineddin and Krauthammer [14] illustrated the effect of reinforcement rate on the damage pattern of RC slabs based on modeling experiments. Although model experiments can be used to effectively analyze the dynamic response of rockfall impact protection structures, the limited experimental conditions cannot fully represent the complex actual impact situations, and the experiments are costly and time-consuming.
With the rapid development of computers, numerical calculation has gradually become an important research tool for solving impact problems due to their economy and efficiency. The finite-element method (FEM) of the shed-cushion layer composite structure has received the most attention in related studies [15,16,17,18,19,20,21]. However, these numerical simulations are carried out using the FEM, focusing on the energy dissipation effect of the cushioning layer, weakening the failure modes of the RC slab, and the results are completely dependent on numerical simulation, lacking effective experimental verification. Relevant studies have shown that the FEM has a good accuracy for small deformation problems and can obtain better results. However, for large deformation instability flow and free surface problems, the convergence and accuracy of numerical calculations lack effective guarantees. Taking the sandy cushion layer of the fluid–solid body as an example, under impact vibration loads, it exhibits obvious flow and large deformation characteristics, and the use of FEM may lead to grid distortion and calculation failure [22]. The discrete element method has complex parameter calibration restrictions [23]. Smoothed particle hydrodynamics (SPH) is a meshless continuum mechanics calculation method. Due to the absence of grid relationships between particles, it can overcome the difficulties of traditional FEM calculations for large deformations and can use the FEM in small deformation areas to improve calculation speed. Therefore, the SPH–FEM coupling algorithm can combine the advantages of both methods [24].
In the design of protective structures, it is commonly necessary to scientifically determine the key characteristic parameters such as the impact energy, peak impact force, and the penetration depth of rockfalls. However, the study of rockfall impacting protection structure is a typical impact dynamic problem, which involves complex large deformation and elastoplastic mechanics. Theoretical formulas and numerical simulations can provide calculation methods. The commonly used theoretical formulas are empirical formulas based on limited test data, which cannot describe the impact effect under all complex impact conditions, and the obtained impact characteristic values are very discrete. Among them, the conservative calculation model will lead to the waste of materials, while the progressive calculation method may lead to the insufficient bearing capacity of the structure. How to choose a reasonable calculation mode is a great challenge for engineers and technicians. Numerical simulation provides a more reliable solution for the design parameters of protective structures. Among them, the FEM is the most widely studied, but for static and quasi-static problems, it can obtain better results. For large deformation problems involving impact and explosion, FEM may face the disadvantages of poor accuracy and convergence.
In this article, the soil cushion layer-RC slab is taken as the research object, and the numerical calculation model is established based on the SPH–FEM coupling algorithm. In the computational model, FEM elements are used for reinforced steel and concrete materials, while SPH particles are used for the sandy soil cushion layer involving large deformation. In order to verify the reliability of the simulation, this article compared and analyzed the simulation results with the test results of the outdoor large-scale rockfall platform model, which verifies the scientific validity of SPH–FEM. In addition, the energy dissipation characteristics of the cushion layer and the deformation mode of the RC slab are analyzed in detail. The research in this article can provide new calculation ideas for the design theory and application of the related rockfall protection structure.
2 Computational theory
2.1 SPH theory
SPH is a pure Lagrangian meshless particle method [25]. Although the traditional Lagrangian FEM has a high computational efficiency, it is prone to grid distortion. SPH can effectively avoid grid distortion and naturally simulate some phenomena, such as material fracture and splashing, but its computational efficiency is lower than that of FEM [26]. In this method, SPH allows the given equation group to be discretized directly through interpolation without defining a grid. The calculation domain is discretized into a series of particles, and the field variables carried by the particles are ultimately transformed into particle motion equations through partial differential equations [27]. The field variables of the particles are obtained through function integration:
where W is a smooth function, and h is a smoothing length used to define the influence area of the smooth function and to determine how many particles affect the interpolation of a specific point.
In order to finally obtain the particle control equations, the continuous integrals in Eq. (1) are replaced by a discretized summation of the particles in the computational domain:
where j represents the number of particles existing in the coordinate x, m j and ρ j are the mass and density of the particles, respectively (Figure 1).
Kernel function.
2.2 SPH–FEM coupling algorithm
In the SPH–FEM coupling algorithm, the contact between particles and finite elements is a key issue in computation. Typically, the interaction force between SPH and FEM is coupled by using the exchange velocity field and displacements field on the interface. In some adaptive algorithms, particles can be defined as nodes and FEM elements as the primary surface. During the coupling computation process, the physical quantities of the finite element and smoothed particle flow are calculated according to all the known velocities and displacements of the nodes. Based on the material model, the stress and strain rate of the material are determined, and the forces on the element nodes and the forces on the particles and their surrounding adjacent nodes are determined after repeated iterations to determine the solution for the entire calculation time. The coupled calculation process is shown in Figure 2.
Flowchart of the SPH–FEM coupling algorithm.
In this work, when considering the impact of rockfalls, the rockfalls are divided into finite elements, and the contact between the cushion layer particles and the rockfalls, as well as the contact with the RC slab elements, is based on the penalty function contact algorithm. When the particles contact with the FEM elements, force contact and transmission are achieved through variables such as particle and element displacement and velocity.
3 Rockfall impact model test
3.1 Test materials
An outdoor large-scale rockfall impact test platform is constructed (Figure 3). Based on the shed structure and pile–slab rock fence structure, two RC slabs with the same configuration are made, one of which is used to study the single impact response, and the other is used to study the damage mode of the RC slab under repeated impacts. The RC slabs' dimensions are 2.4 m × 1.6 m × 0.25 m. The concrete used is P.O 42.5 Portland cement with a strength grade of C30 and a mix ratio of cement:water:sand:stone = 1:0.5:1.5:2.8. The coarse aggregate used is crushed stone with a continuous grain size distribution of 5–15 mm, and the fine aggregate is natural river sand. Two layers of d 14 mm @ 200 mm steel reinforcement meshes are laid orthogonally inside the plate, and a 20 mm thick concrete protective layer is added (Figure 4). For the simulated rockfalls in the field collapse, a concrete cube is poured inside a steel mold, and the numbered rockfall block used in the simulation is labeled C, with a side length of 0.5 m and a mass of 290.8 kg.
Test platform diagram: (a) outdoor impact platform and (b) schematic diagram of impact platform.
RC slab model: (a) schematic diagram of reinforcement for a concrete slab and (b) dimensions of RC slab.
3.2 Test methods
The test data acquisition system mainly includes (1) a 1,000 g piezoelectric sensor placed on the top of the rockfall hammer with a precision of 17.37 Pc·g−1. The falling rock impact force-time curve can be obtained based on the falling rock acceleration; (2) a self-recovering displacement sensor with a range of 100 mm and a repeatability accuracy of 0.01 mm installed at the center point of the lower surface of the concrete slab. The acceleration sampling frequency is 10 kHz, and the displacement sampling frequency is 50 kHz. To ensure the uniformity of the sample and protect the strain gauges and force sensors attached to the surface of the slab, larger pieces of crushed stones are appropriately removed. After each impact test, the soil within a larger area affected by the impact is excavated and then refilled and compacted with the same pressure every 10 cm. In addition, considering the impact caused by the refilling of soil under different impact conditions, the data acquisition instrument needs to be reset before each test. Therefore, the obtained results reflect the load effect of a single impact.
4 Rockfall impact model test
4.1 Numerical model
According to the experiment, the maximum impact height of the rockfall relative to the surface of the cushion layer is about 8.0 m. In the numerical model, the impact height is converted to a maximum impact velocity of 12 m·s−1, which is applied to the surface of the 0.2 m thick cushion layer. The impact point of the rockfall hammer is located at the center point of the upper surface of the cushion layer. Due to a large number of working conditions and the fact that the vertical collision of the structural surface is the most unfavorable condition under the same conditions, oblique collision tests are not conducted in this study. The final established SPH–FEM model is shown in Figure 5.
SPH–FEM numerical model.
In the simulation calculation, except for the cushion layer soil, which is modeled using the particles, all other materials are modeled using the FEM elements. A penalty function is used to simulate the relevant contact between the falling rock and the cushion layer, the cushion layer and the RC slab, and the friction coefficient is 0.4 for the falling rock and the cushion layer and 0.48 for the cushion layer and the RC slab [28].
4.2 Calculation parameters and boundary condition
The previous studies studies have shown that soil is a typical elastic–plastic material with complex deformation characteristics [29] and a very short elastic phase. In order to describe the deformation characteristics of soil, dozens of constitutive models have been proposed, such as the classical Mohr–Coulomb (M–C) model, the Drucker–Prager (D–P) model, the Cam–Clay model [30], the Duncan–Chang model [31], and the like. With the deepening of the research, more scientific soil constitutive models have been proposed, including the 3D fractional elastoplastic constitutive model, which can describe soil under true three-dimensional stress conditions [32], and the fractional order elastoplastic model, which can consider different granular soils under drained and undrained conditions [33]. In this work, to simplify the computation, the cushion layer soil is modeled using the Drucker–Prager constitutive model, which follows the non-associated flow rule. In addition, existing studies indicate that rebar is a typical rate-dependent material. At different loading rates, the mechanical properties of rebar are different. Especially at high strain rates, such as explosion and impact loads, the rebars will show different strength and deformation characteristics than those under static loads. At present, the constitutive models used to simulate reinforcement materials mainly include ideal elastic–plastic model, linear strengthening model, three-line model, and power strengthening model. Among them, the three-line model is more suitable for reinforcing bars with obvious flow amplitude [34]. The three-line model is divided into elastic segment, yield plateau segment, and strengthening segment. The specific expression is as follows:
where E is the elastic modulus of the rebar,
The natural rockfall in the movement will undergo slope friction, air resistance, and mutual collision, causing the speed of the rockfall to be generally about 20 m·s−1. This means that the loading speed of the rockfall impacting protection structure will not reach the strain rate of high-speed impact. Therefore, according to the research of Lin et al. [35], the corresponding parameters of the rebar at a strain rate of 8.4 s−1 are selected in the model of this work (Table 1).
Material parameters of concrete and cushion layer soil [28]
| Material type | Density (kg·m−3) | Modulus of elasticity (MPa) | Poisson’s ratio | Shear expansion angle (°) | Eccentricity | f b0/f co | K | Stickiness parameters |
|---|---|---|---|---|---|---|---|---|
| Concrete | 2,500 | 30,000 | 0.20 | 38 | 0.1 | 1.16 | 0.667 | 0.00001 |
| Cushion layer soil | 1,540 | 35 | 0.28 | Friction angle (°) | Yield strength (kPa) | Shear expansion angle (°) | — | — |
| 28 | 200 | 0 | — | — | ||||
| Rebar | 7,800 | 2.0 × 105 | 0.27 | Yield strength (MPa) | Ultimate strength (MPa) | Extreme tensile strain | — | — |
| 503 | 662 | 0.130 | — | — |
To analyze the damage and failure state of concrete during the impact process, the concrete concrete damage plasticity constitutive is first proposed by Lubliner et al. [36]. This model mainly describes the degradation degree of stiffness of concrete under dynamic loading by defining a damage factor. The relationship between the damage factor and the initial undamaged elastic modulus is given by
where E 0 is the initial elastic modulus of the concrete, and the damage factor d is a function of the uniaxial tensile and compressive damage variables d t and d c.
The specific calculation method of the damage factor can be obtained by using the Simpson integration method proposed by Wang and Yu [37] based on damage mechanics. The relevant material parameters of the numerical model in this research are shown in Table 1.
In addition, due to the significant stiffness difference between the rockfall and the cushion layer, the rockfall hardly deforms during the impact when it contacts with the cushion layer. Therefore, in the simulation, the rockfall is modeled as a rigid body.
The boundary conditions of the model are the same as the actual situation and restrained the three degrees of freedom at the support bottom. To save computation costs, velocity loading is used in the numerical calculation, and the velocity is applied to the center of gravity of the rockfall. The dynamic analysis step is set according to the actual impact time, with a duration of 0.03 s.
4.3 Stability and computational settings of SPH
Although the SPH method is widely used in solving large deformation problems, studies have shown that when it is used to solve the problem of material strength, the unstable motion of particles under stretching or compression will occur, resulting in abnormal aggregation of particles or even complete computational collapse, which is called stress instability [38,39]. According to the stability analysis of Swegle et al. [40], the stress instability of the SPH method depends on the plus or minus of the stress state and the product of the second derivative of the kernel function. Sufficient conditions for stress instability of the SPH method are
where W″ is the second derivative of the kernel function, and σ is the stress state at a point.
The kernel function takes the most common form as the cubic B-spline function (Figure 6):
where s = |x − x′|/h, h is the smooth length, and αd is the normalization coefficient of the function in space.
Stability of cubic spline function.
The stress instability of the cubic spline function is shown in Figure 6. As can be seen from the figure, for cubic splines, since the value of smooth length is generally equal to or slightly greater than the particle spacing, the particles are usually in the tension instability region, which reveals why the tension state is more prone to stress instability than the compression state. In addition, for the calculation of this article, the smooth particle flow area in the soil can only withstand the compressive stress. Therefore, there is no tensile instability in the calculation.
In order to eliminate the compression instability in the calculation process of this article, artificial viscosity can be introduced to solve the problem. The artificial viscosity is mainly composed of linear viscosity Π = −αρlc∇·v , which eliminates the post-wave oscillation, and secondary viscosity Π = βρl 2 c(∇·v)2, which inhibits the numerical oscillation of the strong discontinuity surface of the shock wave, where α and β are dimensionless coefficients, l is the wave propagation length, and c is the sound velocity.
The above two kinds of viscosity are further combined to form a mixed artificial viscosity suitable for the SPH method [41]:
where α Π and β Π are constants, ϕ ij is used to reduce the numerical divergence generated when the particles approach, ρ̅ ij is the average density, and c̅ ij is the average sound velocity.
The mixed artificial viscosity can be divided into two parts, where α Π can be considered as shear or volumetric viscosity, and β Π can be used to prevent particles from penetrating each other. In this work, α Π = 1.0 and β Π = 10 are used.
In general, there are three ways for meshes to be transformed into smooth particles: (1) Time standard. The method can specify that the affected units are transformed according to the given time, regardless of deformation. (2) Strain standard. It specifies the maximum absolute value of the principal strain corresponding to the mesh in the transition area. With the deformation of the element, if the absolute value of the maximum principal strain is greater than the specified threshold, the mesh elements are converted to particles. (3) Stress standard. It specifies the absolute value of the maximum principal stress of the mesh element in the conversion area. As the element deforms, if the absolute value of the maximum principal stress is greater than the specified threshold, the mesh elements is converted to particles.
Corresponding to the actual situation, this article determines the way of converting mesh elements into particles according to the time conversion criterion, that is, from time t = 0, the C3D8R elements in the SPH region are discretized into PC3D type elements. The final single mesh is a 0.02 m hexahedron, the total number is 96,000, the number of nodes is 107,811, and the number of PC3D elements converted into each mesh is 2 × 2 × 2. Therefore, the total number of particles is 7.68 × 105. The kernel function of the particle is the cubic B-spline function, and the mass, density, and smooth length of the particle are calculated automatically by the program according to the mesh element [22].
5 Results and analysis
5.1 Numerical simulation and experimental comparison analysis
Under the same cushion layer, the outdoor test carried out a total of six impact conditions at different speeds (Table 2). The experimental and simulated values of peak impact force and central point displacement of the RC slab under different working conditions are compared and analyzed. The results show that the peak value of rockfall impact force and the displacement at the midpoint of the RC slab increase with the increase in impact velocity. Under the impact condition of 2 m·s−1, the peak impact force is 131.5 kN, and the force amplification effect is up to 45 times compared with the static weight of 2.9 kN. Taking the impact velocity of 12 m·s−1 and the thickness of the cushion layer of 0.2 m as an example, Figure 7 shows the comparison results of the test and numerical simulation of the impact force of rockfall and the duration curve of the vertical displacement of the RC slab. The results indicate that the numerical results are in good agreement with the experimental results. The duration curve of rockfall impact force quickly increases to its peak and then rapidly decreases. As for the displacement of the center point of the RC slab, it continues to increase as the rock is pushed into the cushion layer. Once it reaches its peak, the slab rebounds, causing the displacement to decrease. Finally, the concrete slab undergoes plastic permanent deformation and becomes stable. The peak rockfall impact force obtained from numerical simulation is 471.8 kN, and the measured result is 499.5 kN, with both lasting around 11 ms. The peak displacement of the center point of the RC slab obtained from numerical simulation is 17.8 mm, and the measured value is 18.5 mm.
Comparison between experimental and numerical simulation results
| Cushion layer thickness (m) | Velocity (m·s−1) | Peak force (kN) | Deflection of slab (mm) | ||
|---|---|---|---|---|---|
| Test results | Simulation results | Test results | Simulation results | ||
| 0.2 | 2.0 | 131.5 | 120.6 | 4.8 | 4.2 |
| 4.0 | 206.4 | 198.7 | 6.1 | 5.7 | |
| 6.0 | 267.1 | 279.3 | 8.9 | 9.1 | |
| 8.0 | 352.9 | 343.1 | 12.6 | 12.0 | |
| 10.0 | 408.3 | 386.7 | 15.6 | 15.1 | |
| 12.0 | 499.5 | 471.8 | 18.5 | 17.8 | |
Comparison of numerical simulation results and experimental values: (a) duration curves of impact force and (b) displacement of concrete slab center point.
In general, according to the comparison results of the impact force and the peak displacement of the center point of the concrete slab under different working conditions summarized in Table 2, it can be seen that the relative error between the numerical simulation results and the test results is within 10%, which indicates that the SPH–FEM coupling algorithm adopted in this study is reasonable.
Figure 8 summarizes the duration curve of energy with a cushion layer thickness of 0.2 m and an impact velocity of 12 m·s−1. ALLAE in the figure represents hourglass energy, which is the pseudo-strain energy introduced to prevent the singularity of the global stiffness matrix. Generally, when ALLAE accounts for less than 15% of the total energy of the system, the results are considered reliable; otherwise, the grid needs to be refined, or the hourglass energy should be controlled. ALLFD represents the change in energy due to friction during impact. ALLKE represents the change process of kinetic energy of the system, and for this research, it can specifically describe the change of kinetic energy in the impact process of rockfall. ALLPD represents the change process of the plastic strain energy of the system. In this model, the main plastic strain energy is mainly due to the large deformation of the cushion layer. ALLSE represents the elastic recoverable strain energy of the system. In this article, ALLSE mainly refers to the elastic strain energy generated by RC slabs.
Variation of system energy in numerical simulations (h = 0.2 m, v = 12 m·s−1).
The results show that the kinetic energy of rockfalls is mainly transformed into the plastic deformation energy, frictional energy, and internal energy of the system. According to the results in Figure 8, the total energy of the system under the 12 m·s−1 simulation condition is 20.9 kJ, and the maximum hourglass energy generated during the impact process is 0.4 kJ, accounting for only 1.9% of the total energy. Previous studies have shown that the pseudo-strain energy accounts for less than 15% of the total energy, and the numerical simulation results are accurate. In addition, the kinetic energy of the rockfall is mainly converted into the plastic strain energy of the system. Under the 20.9 kJ impact energy, the plastic strain energy of the system is 16.1 kJ, accounting for more than 77.0% of the total energy.
5.2 Energy dissipation analysis
In the impact process, the total energy of the system is mainly converted into plastic deformation energy. Further analysis is conducted on the energy dissipation law of the system at different velocities. Figure 9(a) shows the plastic energy consumption characteristics of the cushion layer at different velocities. The results show that when the impact velocity is less than 6 m·s−1, all the plastic deformation energy produced by the system is absorbed by the cushion layer. Once the velocity is greater than 6 m·s−1, the plastic deformation energy produced by the cushion layer increases with the increase in velocity. In addition to being converted into the plastic deformation energy of the cushion layer, a small portion of the system energy is also converted into the plastic deformation energy of the RC slab. According to the relationship between the plastic deformation energy of the cushion layer and the plastic deformation energy of the system shown in Figure 9(b), when the impact velocity is 8 m·s−1, the plastic deformation energy produced by the system is 6.8 kJ, and the plastic deformation energy absorbed by the cushion layer is 6.3 kJ, accounting for 92.3%; when the impact velocity is 10 m·s−1, the plastic deformation energy produced by the system is 10.9 kJ, and the plastic deformation energy absorbed by the cushion layer is 9.8 kJ, accounting for 89.9%; when the impact velocity is 12 m·s−1, the plastic deformation energy produced by the system is 16.1 kJ, and the plastic deformation energy absorbed by the cushion layer is 14.2 kJ, accounting for 88.2%. With the increase in impact energy, the soil particle cushion layer is compressed and hardened, resulting in a decrease in the proportion of plastic deformation energy absorbed by the cushion layer, and a portion of the plastic deformation energy is mainly absorbed by the RC slab. Overall, within the scope of this study, due to the existence of the cushion layer, the sand cushion layer can absorb more than 85% of the energy, effectively protecting the RC slab from undergoing destructive damage and achieving graded energy dissipation against impacts, thereby ensuring that the RC slab is in an elastic working state under small impact energy; under large impact energy, the RC slab does not undergo destructive damage.
Energy consumption characteristics of cushion layers in numerical models (h = 0.2 m, v = 12 m·s−1): (a) duration curves of structural plastic energy dissipation at different velocities and (b) percentage of structural plastic energy consumption.
5.3 Changes in impact force
Figure 10 summarizes the pattern of change in impact force at different impact velocities. The results showed that the peak impact force during the rockfall impact increased linearly with velocity (Figure 10(a)). In addition, the duration of the impact force at different velocities varied, showing the characteristics that the greater the rockfall impact velocity, the shorter the impact time and the greater the impact force. This is consistent with the results of classical impact dynamics theory. Furthermore, the maximum impact force of the rockfall at different impact velocities calculated using the Hertz theory is compared with the numerical simulation results (Figure 10(b)). When the velocity is low, the two results are in good agreement, but as the velocity increased, the relative error between the results obtained from the Hertz contact theory and the simulation results increased, which is consistent with the results of Wang et al. [42].
Trends of impact force: (a) simulation results of impact force duration curve and (b) comparison of experimental and theoretical values of impact force.
5.4 Analysis of rockfall impact process
The soil cushion layer is composed of granular particles with typical fluid–solid properties, and traditional FEM cannot intuitively simulate the flow characteristics of sand under impact. In addition, previous studies have shown that the use of FEM to simulate large deformations caused by impact and explosion often leads to convergence problems due to grid distortion. Figure 11 summarizes the results of using the SPH–FEM coupling algorithm to simulate the process and characteristics of a rockfall pressed into a 0.2 m thick cushion layer at an impact velocity of 12 m·s−1. When the rockfall contacts the cushion layer, the sand flows outward due to compression and shear, especially at the edges and corners of the rockfall where there is greater stress concentration. As the rockfall is further pressed into the cushion layer, the stress on the sand continues to increase until the rockfall reaches a depth of 3.22 cm after 5 ms of impact. When the impact time reaches 8 ms, the rockfall reaches the maximum depth of 4.82 cm, and the maximum stress on the sand is 2.5 MPa. Later, when the rockfall begins to bounce back, the impact depth and stress on the cushion layer decrease. At the end of the impact, the depth of the impact crater on the cushion layer is 4.64 cm, and the corresponding Mises stress of the sand is 0.9 MPa, while the measured depth of the cushion layer in the field is 5.2 cm (Figure 11(b)). In addition, Figure 11(a) shows the maximum equivalent plastic strain of the cushion layer corresponding to an impact time of 8 ms. One can see that the maximum plastic strain of the cushion layer is 1.59, indicating that the soil material undergoes a large deformation during the impact. Figure 11(c) shows a magnified contour map of the deformation of the cushion layer, which indicates that after the impact, the shape of the impact crater on the cushion layer is the same as that on the bottom surface of the falling rock. Since the particles can slip and slide in the flow process, the total deformation of the soil in Figure 11(c) is slightly larger than the vertical displacement. Furthermore, the shape of the impact crater exhibits raised sand at both ends and slight bulging in the central area due to rebound after compaction. Overall, the numerical simulation using the SPH–FEM coupling algorithm is more realistic in reproducing the flowability of sand particles during the rockfall impact than using the FEM.
Process and characteristics of rockfall pressed into the cushion layer. (a) Equivalent plastic strain, (b) crater morphology of test results, and (c) deformation of the cushion layer (magnification 5×).
Figure 12 shows the duration curve of the motion trajectory of the center of mass of the rockfall. It should be noted that, in order to increase calculation accuracy, the model sets the initial distance between the lower surface of the rockfall and the upper surface of the cushion layer as 2 cm. The results indicate that the motion trajectory of the center of mass of the rockfall is basically consistent with the variation in the impact depth of the cushion layer. During the impact, the rockfall first continuously presses into the cushion layer, and their velocity decreases with time. When the velocity of the rockfall decreases to 0 m·s−1, the deformation of the cushion layer reaches the maximum, and the displacement of the center of mass of the rockfall also reaches the maximum. After that, the cushion layer rebounds, and the block center-of-mass displacement direction reverses and eventually comes to rest.
Trajectory of block center of mass.
5.5 Damage characteristics of RC slabs
In order to study the failure mode of concrete slabs under the low-speed impact of falling rocks, seven repeated impact tests are conducted while keeping the thickness of the cushion layer at 0.2 m and the velocity at 12 m·s−1. This study uses a restart function to simulate the cumulative damage process of the structure under multiple impact loadings [43]. Figure 13 shows the evolution of the failure process of the RC slab under low-speed impact loading of rockfalls. The results show that under the second impact, tensile cracking occurs at the centerline of the RC slab (Figure 13a). When the number of impacts increases to 4, the longitudinal centerline crack further extends to the surface of the concrete slab, and more closely spaced tensile cracks occur between the first cracks (Figure 13b), with tensile damage factor near the centerline reaching above 0.8. At the end of the impact, the centerline crack extends completely to the top of the slab, and the vertical cracks on both sides further expand, forming significant shear cracks (Figure 13c). The maximum deflection of the slab during the test is 9.7 cm, and the maximum deflection obtained from the numerical simulation is 8.3 cm. Overall, the deformation mode of the slab has both bending and shearing characteristics, and the numerical simulation results are basically consistent with the experimental observations.
Failure characteristics of RC slab: (a) initiation of cracking in RC slab, (b) crack propagation in the concrete slab, and (c) the final failure mode of the concrete slab.
6 Discussion
In this work, the SPH–FEM coupling algorithm is used to study the dynamic response of cushion layer-RC slab under the rockfall impact. The results of the numerical simulation agreed with those of the physical model, which proved that the SPH–FEM coupling algorithm combines the advantages of FEM and SPH. For impact and explosion models, the SPH method can reflect the large deformation of the structure more accurately and avoid the non-convergence problem caused by mesh distortion. SPH particles and FEM elements are used to simulate the process of rockfall pressing into the cushion layer, as shown in Figure 14.
In the initial propagation stage, stress waves transfer forces to particles/FEM elements through contact with each other, thus causing particles to move or elements to deform.
In the compression stage, violent collisions occur between particles, resulting in dislocation and irregular movement between particles due to the different strength of particles at different locations. Unlike particle motion, the FEM elements do not intrude into each other and only exhibit significant deformation under the action of stress waves. Obviously, given the continuity of FEM, the SPH method is more reasonable for the actual compression process.
In the stable stage, the final displacements corresponding to both methods are consistent. The random motion of SPH particles may be an important reason that the final plastic deformation of the cushion layer is slightly larger than that of the FEM method.
The process of rockfall impacting protection structure. (a) Original state. (b) Compressed state. (c) Stable state.
This study mainly considered the contact between the rockfall plane and the cushion layer plane, no matter in the experiments or numerical simulations. However, in fact, the shape of the natural rockfall is irregular. For example, when a sharp rockfall impacts the protective structure, a puncture collision is more likely to occur. Compared with the contact between two planes, puncture collision is more likely to cause large mesh deformation and even computational failure [44]. In order to draw a comparison with the test, the puncture damage is not reflected in this research. In addition, oblique impact is more realistic than frontal impact. However, the SPH–FEM coupling algorithm adopted in this work still has important reference value for the engineering case of low-speed rockfall collision. Finally, the impact velocity has great influence on the convergence of numerical simulation results. In this study, a low-speed impact condition with a maximum impact velocity of 12 m·s−1 is considered. Therefore, the dynamic response of the material under high strain rate conditions still needs to be further investigated.
7 Conclusion
The SPH–FEM coupling method can fully utilize the advantages of FEM and SPH to solve the convergence problem of large deformation. In this study, the SPH–FEM coupling algorithm is used to investigate the dynamic response of the cushion layer-RC slab subjected to the impact of rockfalls, and the numerical simulation results are compared and verified with the results of a large-scale rockfall impact test, and the following conclusions are drawn:
When the impact velocity is 12 m·s−1, the numerical simulation results and the measured results of the rockfall impact force and the displacement of the RC slab have good agreement. The peak impact force obtained from numerical simulation is 471.8 kN, and the measured value is 499.5 kN; the impact duration is approximately 11 ms. The peak displacement of the center point of the RC slab obtained from the numerical simulation is 17.8 mm, while the measured result is 18.5 mm.
The soil cushion layer can absorb more than 85% of the impact energy, effectively protecting the RC slab from damage. When subjected to small impact energy, the RC slab remains in an elastic state; when subjected to large impact energy, the RC slab does not suffer catastrophic damage.
During the impact, the maximum plastic strain of the cushion layer is 1.59, indicating that the sand body undergoes large deformation during the impact. The impact crater formed by the cushion layer had the same shape as the bottom surface of the rockfall. The crater exhibited raised sand at both ends and slight bulging in the central area due to rebound after compaction. The numerical simulation using the SPH–FEM coupling algorithm is more practical in reproducing the flowability of sand particles during the rockfall impact than using the FEM.
The cumulative damage process of the RC slab is studied. At the end of the impact, the centerline crack completely extended to the top of the slab, the vertical cracks on both sides further expanded, and significant shear cracks are formed. Therefore, the deformation mode of the slab has both bending and shearing characteristics.
Acknowledgments
The authors thank the faculty of Geosciences and Environmental Engineering Laboratory of Southwest Jiaotong University for the support of experiments in this research.
-
Funding information: This work is supported by the Shandong Province Natural Science Foundation Youth Branch (Grant No. ZR2021QE209), Foundation of China Postdoctoral Science Foundation (Grant No. 2022M723687), and Doctoral Science and Technology Startup Foundation of Shandong University of Technology (Grant No. 420048).
-
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
-
Conflict of interest: The authors state no conflict of interest.
References
[1] Wu, C., J. Ding, and G. Huang. Dynamic response analysis of a new composite shed-tunnel structure under the rockfall impact. Building Structure, Vol. 48, No. S2, 2018, pp. 546–549.Search in Google Scholar
[2] Matsimbe, J. Comparative application of photogrammetry, handmapping and android smartphone for geotechnical mapping and slope stability analysis. Open Geosciences, Vol. 13, No. 1, 2021, pp. 148–165.10.1515/geo-2020-0213Search in Google Scholar
[3] Haberler, M., M. Huber, T. Wunderlich, and C. Kandler. Fuzzy system based analysis of deformation monitoring data at Eiblschrofen rockfall area. Journal of Applied Geodesy, Vol. 1, No. 1, 2007, pp. 16–25.10.1515/jag.2007.003Search in Google Scholar
[4] Pichler, B., C. Hellmich, and H. A. Mang. Protection of infrastructure subjected to rockfall. Journal of the Mechanical Behavior of Materials, Vol. 16, No. 1–2, 2005, pp. 103–112.10.1515/JMBM.2005.16.1-2.103Search in Google Scholar
[5] Lei, P., T. Zhou, X. Wang, Z. Zhen, and Y. Kou. Research on punching failure of shed under impact load. Chinese Journal of Underground Space & Engineering, Vol. 1504, 2019, pp. 1091–1097.Search in Google Scholar
[6] Wang, D., Y. Liu, X. Pei, and S. Shi. Elasto-plastic dynamic responses of reinforced concrete slab sunder rockfall impact. Journal of Southwest Jiaotong University, Vol. 5106, 2016, pp. 1147–1153.Search in Google Scholar
[7] Wang, M., C. M. Song, D. R. Wang, and Y. F. Pan. Calculational methods of reinforced concrete slabs to low velocity impact. Acta Armamentarii, Vol. 2004, No. 06, 2004, pp. 672–678.Search in Google Scholar
[8] Su, H., M. Wang, and C. Song. A unified method for critical thickness of reinforced concrete plates under low velocity impact. Protective Engineering, Vol. 3805, 2016, pp. 50–56.Search in Google Scholar
[9] Mougin, J., P. Perrotin, M. Mommessin, J. Tonnelo, and A. Agbossou. Rock fall impact on reinforced concrete slab: An experimental approach. International Journal of Impact Engineering, Vol. 31, No. 2, 2005, pp. 169–183.10.1016/j.ijimpeng.2003.11.005Search in Google Scholar
[10] Kishi, N., H. Konno, K. Ikeda, and K. G. Matsuoka. Prototype impact tests on ultimate impact resistance of PC rock-sheds. International Journal of Impact Engineering, Vol. 27, No. 9, 2002, pp. 969–985.10.1016/S0734-743X(02)00019-2Search in Google Scholar
[11] Bhatti, A. Q., S. Khatoon, A. Mehmood, A. Dastgir, and N. Kishi. Numerical study for impact resistant design of full scale arch type reinforced concrete structures under falling weight impact test. Journal of Vibration & Control, Vol. 18, No. 9, 2011, pp. 1275–1283.10.1177/1077546311419176Search in Google Scholar
[12] Zhang, S., J. Zhang, L. Tao, H. Liu, and L. Gao. Tests for cushion performance of a protective layer with S-shaped steel joist and sandwich slab under rockfall impact. Journal of Vibration and Shock, Vol. 3624, 2017, pp. 148–155.Search in Google Scholar
[13] Gu, S., F. Peng, Z. Yu, and J. Li. An experimental study on the damage effects of the concrete slabs under low-velocity impact. Journal of Vibration & Shock, Vol. 3824, 2019, pp. 107–114.Search in Google Scholar
[14] Zineddin, M. and T. Krauthammer. Dynamic response and behavior of reinforced concrete slabs under impact loading. International Journal of Impact Engineering, Vol. 34, No. 9, 2007, pp. 1517–1534.10.1016/j.ijimpeng.2006.10.012Search in Google Scholar
[15] Wang, X., B. Ren, Q. Wang, and W. Huang. Refined dynamic response simulation of sand-EPE-tunnel shed roof with rockfall impact. Journal of Beijing Jiaotong University, Vol. 4603, 2022, pp. 118–127.Search in Google Scholar
[16] Wang, J., P. Zhao, S. Yuan, L. Li, and L. Xie. Numerical study on impact resistance of steel shed gallery with composite cushion. China Civil Engineering Journal, Vol. 51, 2018, pp. 7–13.Search in Google Scholar
[17] Wang, M., S. Shi, and Y. K. Yang. Numerical simulation of a flexible rock shed under the impact of a rockfall. Engineering Mechanics, Vol. 31, No. 5, 2014, pp. 151–157.Search in Google Scholar
[18] Hu, X., X. Mei, Y. Yang, G. Luo, and J. Wu. Dynamic response of pile-plate rock retaining wall under impact of rockfall. Journal of Engineering Geology, Vol. 2701, 2019, pp. 123–133.Search in Google Scholar
[19] Li, L., S. Yuan, L. Xie, and P. Zhao. Study on the cushioning properties of shed with EPE cushion under rock-fall load. Sichuan Building Science, Vol. 4203, 2016, pp. 46–49.Search in Google Scholar
[20] Wang, D., S. He, Y. Wu, and X. Li. Cushioning effect of rock sheds with EPS cushion on rock-falls action. Journal of Vibration & Shock, Vol. 33, No. 4, 2014, pp. 199–203 + 214.Search in Google Scholar
[21] Yan, P., J. Zhang, and Q. Fang. Numerical simulation of the effects of falling rock’s shape and impact pose on impact force and response of RC slabs. Construction & Building Materials, Vol. 160, No. Jan. 30, 2018, pp. 497–504.10.1016/j.conbuildmat.2017.11.087Search in Google Scholar
[22] Luo, J., J. Xiao, K. Ma, and J. Mao. SPH-FEM coupled method for analyzing a hemispherical shell impact soil. Journal of Vibration & Shock, Vol. 36, No. 17, 2017, pp. 195–199 + 230.Search in Google Scholar
[23] Wang, Y., J. Li, Z. Li, G. Feng, H. Wu, and J. He. Assessment of rockfall impact force by particle flow code numerical simulation based on discrete element model. Journal of Southwest Jiaotong University, Vol. 51, 2016, pp. 22–29.Search in Google Scholar
[24] Johnson, G. R. and R. A. Stryk. Conversion of 3D distorted elements into meshless particles during dynamic deformation. International Journal of Impact Engineering, Vol. 28, No. 9, 2003, pp. 947–966.10.1016/S0734-743X(03)00012-5Search in Google Scholar
[25] Lucy, L. B. A numerical approach to the testing of the fission hypothesis. Astrophys Journal, Vol. 812, 1977, pp. 1013–1024.10.1086/112164Search in Google Scholar
[26] Huang, Y. and L. Hao. The state of the art of SPH method applied in geotechnical engineering. Chinese Journal of Geotechnical Engineering, Vol. 30, No. 2, 2008, pp. 256–262.Search in Google Scholar
[27] Liu, G. R. and M. B. Liu. Smoothed particle hydrodynamics-a mesh free particle method, World Scientific Publishing Company, New Jersey, 2003.10.1142/9789812564405Search in Google Scholar
[28] Luo, J. Impact resistance of composite beam and steel vierendeel plate. Doctoral dissertation, Gui zhou University, Gui zhou, 2018.Search in Google Scholar
[29] Dafalias, Y. F. Bounding surface plasticity. I: Mathematical foundation and hypoplasticity. Journal of Engineering Mechanics, Vol. 112, No. 9, 1986, pp. 966–987.10.1061/(ASCE)0733-9399(1986)112:9(966)Search in Google Scholar
[30] Roscoe, K. H., A. N. Schofield, and A. Thurairajah. Yielding of clays in states wetter than critical. Geotechnique, Vol. 13, No. 3, 1963, pp. 211–240.10.1680/geot.1963.13.3.211Search in Google Scholar
[31] Duncan, J. and C. Y. Chang. Nonlinear analysis of stress and strain in soils. Proc. Paper: J Soil Mech & Found Division ASCE, Vol. 96, 1970, pp. 1629–1653.10.1061/JSFEAQ.0001458Search in Google Scholar
[32] Lu, D., J. Liang, X. Du, C. Ma, and Z. Gao. Fractional elastoplastic constitutive model for soils based on a novel 3D fractional plastic flow rule. Computers & Geotechnics, Vol. 105, 2019, pp. 277–290.10.1016/j.compgeo.2018.10.004Search in Google Scholar
[33] Sun, Y. and Y. Xiao. Fractional order plasticity model for granular soils subjected to monotonic triaxial compression. International Journal of Solids & Structures, Vol. 118–119, 2017, pp. 224–234.10.1016/j.ijsolstr.2017.03.005Search in Google Scholar
[34] Zhou, H., M. Li, J. Duan, and C. Song. Research on dynamic constitutive model and model parameters of steel bars. Journal of Ordnance Equipment Engineering, Vol. 4308, 2022, pp. 193–202.Search in Google Scholar
[35] Lin, F., X. Gu, X. Kuang, and X. Yin. Constitutive models for reinforcing steel bars under high strain rates. Journal of Building Materials, Vol. 11, No. 1, 2008, pp. 14–20.Search in Google Scholar
[36] Lubliner, J., J. Oliver, S. Oller, and E. Oñate. A plastic-damage model for concrete. International Journal of Solids & Structures, Vol. 25, No. 3, 1989, pp. 299–326.10.1016/0020-7683(89)90050-4Search in Google Scholar
[37] Wang, Z. and Z. Yu. Concrete damage model based on energy loss. Journal of Building Materials, Vol. 7, No. 4, 2004, pp. 365–369.Search in Google Scholar
[38] Dyka, C. T. and R. P. Ingel. An approach for tension instability in smoothed particle hydrodynamics (SPH). Computers & Structures, Vol. 57, No. 4, 1995, pp. 573–580.10.1016/0045-7949(95)00059-PSearch in Google Scholar
[39] Chen, J. K., J. E. Beraun, and C. J. Jih. Improvement for tensile instability in smoothed particle hydrodynamics. Computational Mechanics, Vol. 23, No. 4, 1999, pp. 279–287.10.1007/s004660050409Search in Google Scholar
[40] Swegle, J. W., D. L. Hicks, and S. W. Attaway. Smoothed particle hydrodynamics stability analysis. Journal of Computational Physics, Vol. 116, No. 1, 1995, pp. 123–134.10.1006/jcph.1995.1010Search in Google Scholar
[41] Monaghan, J. J. and R. A. Gingold. Shock simulation by the particle method SPH. Journal of Computational Physics, Vol. 52, No. 2, 1983, pp. 374–389.10.1016/0021-9991(83)90036-0Search in Google Scholar
[42] Wang, X., T. ZHou, J. Shi, and Y. Xia. Theoretical and LS-DYNA simulation study of based on the theory of free-fall rockfall’s impact on soil layer. Journal of Beijing Jiaotong University, Vol. 4304, 2019, pp. 9–17.Search in Google Scholar
[43] Meng, Y., W. Yi, and F. Huang. Experimental study on impact behavior of hinged RC beams. Journal of Building Structures, Vol. 39, No. 4, 2018, pp. 82–90.Search in Google Scholar
[44] Chuzel, Y., R. Ortiz, and A. Combescure. Three dimensional SPH–FEM gluing for simulation of fast impacts on concrete slabs. Computers & Structures, Vol. 89, No. 23, 2011, pp. 2484–2494.10.1016/j.compstruc.2011.06.002Search in Google Scholar
© 2024 the author(s), published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.
