Numerical Investigation of the Dynamic Response of a Sand Cushion with Multiple Rockfall Impacts

: A shed cave structure with a sand cushion is often used as a protective structure for rockfall disasters. Because of the randomness and unpredictability of rockfall disasters, the cushions of shed caves often suffer multiple impacts from rockfalls. Aiming at the problem of multiple impacts of rockfall, this paper uses the three-dimensional discrete element method to study the dynamic response of multiple rockfall impacts on sand cushions from different heights. Before conducting large-scale simulation studies, the input parameters in the numerical model are veriﬁed with data from laboratory experiments. Analyzing the simulation results shows that when the same point is impacted multiple times, the maximum impact force and the maximum penetration depth will increase with the number of impacts. According to the numerical results, a calculation formula of the maximum impact force that considers the number of impacts is ﬁtted. At the same time, considering the impact response when the rockfall impacts different positions multiple times, the distance range that the subsequent impact is not affected by the previous impact is given. The signiﬁcance of studying the multiple impacts of rockfalls is shown by a numerical study of rockfalls impacting a sand cushion multiple times from different heights, and it provides a reference for the design of rockfall disaster-protection structures in practical engineering.


Introduction
Rockfall disasters are one of the most common geological disasters near mountain roads. The occurrence of rockfall disasters has had a great impact on the economic development of mountainous areas. In China, with the advancement of the western development strategy, it is even more necessary to reasonably and effectively protect against rockfall disasters. In addition, casualties and damage to roads and railways due to rockfall disasters should be avoided. The shed cave structure with a buffer layer is a typical structure for rockfall-protection in mountainous areas and it is widely used in rockfall disaster-protection in mountainous areas [1][2][3].
Sand is a cushion material commonly used in rockfall-protection engineering [4]. It is a conveniently obtained material that has good durability, high economic benefit, and good buffering performance. Many scholars have studied the buffer performance of sand buffer layers. A series of laboratory experiments were conducted to investigate the effect of the dry density and thickness of the sand pad on the impact pressure applied to the soil surface, the earth pressure at the bottom of the mold, and the transferability of the impact pressure [5]. In recent years, the discrete element method (DEM) has emerged as a suitable numerical tool for analyzing the impact of rockfalls from micro-to macroscales [6]. The calibration 2 of 21 of validated DEM models with corresponding experimental data provides researchers with data that are nearly impossible to obtain experimentally. Using the commercial DEM software PFC3D 5.0, the initial kinetic energy of a rockfall can be set to 5000 kJ, which is difficult to achieve experimentally, to study the impact of the rockfall on bunkers covered by a soil buffer layer [7]. It is also possible to study the energy propagation and block bouncing in the process of a rock block impacting a granular medium, as well as the evolution of the shock-induced force chain and its relationship with the global mechanical response of the granular buffer layer through the 3D discrete element model [8,9]. At the same time, the effect of the particle size on the buffering efficiency of the soil buffer layer can be studied by combining the experiments [10]. The influence of the sand buffer layer thickness and porosity on the buffer performance can also be analyzed by the discrete element numerical method [11]. In these studies, the DEM was found to be an effective method for studying the impact response of sand buffers.
Because of the high frequency and randomness of rockfall disasters, it is difficult to conduct monitoring and provide feedback in time. Often, the shed-hole cushion cannot be repaired and replaced in time, and it will be impacted again. In addition, rockfall impact is often accompanied by multiple rocks. Therefore, the shed cavity cushion is often repeatedly impacted. Such an occurrence will compact the sand buffer layer and reduce the ability to disperse the impact force of falling rocks, which will lead to the destruction of the shed cave structure in the long run. At present, few reports are available on the multiple impacts of rockfall on shed cave cushions. A new type of energy dissipator is designed for the place where the rockfall occurs so that the structural energy-dissipating rock shed can withstand the multiple impacts of the rockfall [12]. Tests of spherical rockfall impacting a sand buffer layer and foam composite cushion layer were performed. In these multiple impact laboratory tests, after the first impact, with an increasing number of impacts of the buffer layer, the acceleration of the falling weight is greater [13,14]. The impact force is a key parameter in the structural design of shed tunnels. In the case of multiple impacts, few studies have reported on the impact force calculation after each impact. Simultaneously, under multiple impacts, the impact position is not fixed. If the distance between two adjacent impact positions is close, the effect of the second impact will be affected by the first impact. However, at present, research is lacking on the influence of the distance between different impact locations on the impact effect for multiple impacts.
This paper studies and analyzes the phenomenon of multiple shocks in rockfall disasters. A small-scale laboratory multiple shock test was carried out to verify the discrete element model. The dynamic process of a large-scale rockfall impacting a sand buffer layer multiple times is simulated using a calibrated numerical model. The difference between the dynamic response of multiple rockfall impacts and a single impact is explored, and a calculation method for the impact force during multiple impacts is given in combination with the existing rockfall impact force calculation formula. Furthermore, considering the influence of the distance between impact points on the impact effect during two adjacent impacts, the distance effect on the impact effect at different impact positions is given. The experimental and numerical study on the multiple impacts of rockfall on the sand buffer layer provides a basis for considering multiple impact problems in practical engineering.

Reduced-Scale Impact Test
Firstly, the dynamic response of the sand cushion under multiple impacts of rockfall is studied by a laboratory impact test.

Overview of the Test
In this paper, the impact test of rockfall on a sand buffer layer is performed in a self-designed drop weight impact test device. In this laboratory test, spherical rockfalls made four consecutive impacts on the center point of the sand buffer layer from a height of 1 m. The rockfall has a radius of 0.057 m, a mass of 1.7 kg, and a density of 2187 kg/m 3 . The

Overview of the Test
In this paper, the impact test of rockfall on a sand buffer layer is performed in a selfdesigned drop weight impact test device. In this laboratory test, spherical rockfalls made four consecutive impacts on the center point of the sand buffer layer from a height of 1 m. The rockfall has a radius of 0.057 m, a mass of 1.7 kg, and a density of 2187 kg/m 3

Test Method
In the laboratory test, the acceleration time-history curve of the falling weight is collected by the acceleration sensor installed on the upper surface of the drop hammer, and the impact force time-history curve is obtained according to Newton's second law. In order to avoid errors in the test as much as possible, three groups of tests with multiple impacts under the same conditions were carried out. Judging from the signals collected by the accelerometer, the acceleration trends in the three groups of tests are the same, and the peak accelerations are slightly different, but the difference is small. The impact force results are shown in Section 3.2.

Numerical Model Establishment and Validation
As a common and reliable scientific research method, numerical simulations are often used in civil engineering research [15][16][17]. As discrete element software, PFC3D is mainly used to study the micromechanical behavior between particles [18,19]. The ball and the wall are the basic components of the software, and different numerical models are established by assigning them different parameters. The choice of the particle radius in the model is crucial and determines the speed of the numerical simulations. After the model is built, mesoscopic parameters must be assigned to the particles to achieve the macroscopic properties of the material, that is, to match the macroscopic parameters by selecting the appropriate mesoscopic parameters. No clear conversion relationship holds between macro parameters and meso parameters. In this paper, the model parameters in the numerical model are corrected by performing laboratory experiments to determine the particle size, stiffness, damping coefficient, and other simulation parameters.
Simulating the rockfall impact buffer layer mainly involves three types of rockfall, buffer layers, and protective structures. Considering that the rockfall is less likely to break after impacting the sand, the rockfall is simulated by a single rigid sphere (ball) [20], which is colored in red. The sand buffer layer is generated by a Gaussian distribution of spheres with different radii, and the thickness of the buffer layer is set by controlling the number and position of the spheres [21]. The size of the spherical particles is a key parameter. The

Test Method
In the laboratory test, the acceleration time-history curve of the falling weight is collected by the acceleration sensor installed on the upper surface of the drop hammer, and the impact force time-history curve is obtained according to Newton's second law. In order to avoid errors in the test as much as possible, three groups of tests with multiple impacts under the same conditions were carried out. Judging from the signals collected by the accelerometer, the acceleration trends in the three groups of tests are the same, and the peak accelerations are slightly different, but the difference is small. The impact force results are shown in Section 3.2.

Numerical Model Establishment and Validation
As a common and reliable scientific research method, numerical simulations are often used in civil engineering research [15][16][17]. As discrete element software, PFC3D is mainly used to study the micromechanical behavior between particles [18,19]. The ball and the wall are the basic components of the software, and different numerical models are established by assigning them different parameters. The choice of the particle radius in the model is crucial and determines the speed of the numerical simulations. After the model is built, mesoscopic parameters must be assigned to the particles to achieve the macroscopic properties of the material, that is, to match the macroscopic parameters by selecting the appropriate mesoscopic parameters. No clear conversion relationship holds between macro parameters and meso parameters. In this paper, the model parameters in the numerical model are corrected by performing laboratory experiments to determine the particle size, stiffness, damping coefficient, and other simulation parameters.
Simulating the rockfall impact buffer layer mainly involves three types of rockfall, buffer layers, and protective structures. Considering that the rockfall is less likely to break after impacting the sand, the rockfall is simulated by a single rigid sphere (ball) [20], which is colored in red. The sand buffer layer is generated by a Gaussian distribution of spheres with different radii, and the thickness of the buffer layer is set by controlling the number and position of the spheres [21]. The size of the spherical particles is a key parameter. The actual size of the sand particles is approximately 0.1 mm. However, simulating a sand pad with a particle size of 0.1 mm is impractical, especially on a practical engineering scale. Considered comprehensively, the particle size range of the sand buffer layer in this model will be larger than the actual particle size range [4,22]. The main analysis in the simulation is the impact of rockfall on the buffer layer. The protective structure is directly treated as a wall element. That is, the surroundings of the sand buffer layer are constrained by rigid wall elements. In this study, a linear model was used to simulate the sand buffer layer, the contact between the sand buffer layer and the rockfall, and the contact between the sand buffer layer and the protective structure. When applying PFC3D software to simulate the dynamic impact process of rockfall under vertical falling conditions, factors such as gravity, viscous resistance, and viscous damping need to be considered [23]. The acceleration of gravity in the model is 9.8 m/s 2 , the damping is set to 0.01, and the friction coefficient is set to 0.5. The normal contact force of the rigid sphere in the simulation is the rockfall impact force. Since this paper considers multiple impacts of falling rocks, the rigid sphere is deleted after each impact, and a new rigid sphere is established to continue the impact.

Numerical Model Establishment
Before conducting the analysis and research of a rockfall impacting a sand buffer layer multiple times, a numerical model of the same size as the laboratory test was established to verify the feasibility of using PFD3D software to calculate the sand buffer layer of the rockfall impact.
According to the description of the test, rockfalls and cushions of the same size were established in PFC3D software, and the dynamic impact process was simulated for a height of 1 m. The numerical model is shown in Figure 2. The density of rockfall shall refer to the value in the test. The size of the falling rock shall be replaced by a sphere with the same volume as the falling hammer in the test. Only one particle is used to represent the falling stone, which is because the falling rock is regarded as a rigid body. The size of the sand buffer layer is 1 m × 1 m × 0.3 m. The material parameters of the falling rock and the sand cushion are the same as those in the laboratory experiment, and the other input parameters in the model are shown in Table 1. The bottom and four sides of the cushion is a fixed constraint.
the contact between the sand buffer layer and the rockfall, and the contact be sand buffer layer and the protective structure. When applying PFC3D softwa late the dynamic impact process of rockfall under vertical falling conditions, fa as gravity, viscous resistance, and viscous damping need to be considered [23]. eration of gravity in the model is 9.8 m/s 2 , the damping is set to 0.01, and the fr ficient is set to 0.5. The normal contact force of the rigid sphere in the simula rockfall impact force. Since this paper considers multiple impacts of falling rock sphere is deleted after each impact, and a new rigid sphere is established to co impact.

Numerical Model Establishment
Before conducting the analysis and research of a rockfall impacting a s layer multiple times, a numerical model of the same size as the laboratory test lished to verify the feasibility of using PFD3D software to calculate the sand b of the rockfall impact.
According to the description of the test, rockfalls and cushions of the sam established in PFC3D software, and the dynamic impact process was simu height of 1 m. The numerical model is shown in Figure 2. The density of rockfal to the value in the test. The size of the falling rock shall be replaced by a sphe same volume as the falling hammer in the test. Only one particle is used to rep falling stone, which is because the falling rock is regarded as a rigid body. The sand buffer layer is 1 m × 1 m × 0.3 m. The material parameters of the falling ro sand cushion are the same as those in the laboratory experiment, and the othe rameters in the model are shown in Table 1. The bottom and four sides of the c fixed constraint.

Numerical Model Validation
The sphere rockfall with a mass of 1.7 kg impacts the sand cushion multiple times at the falling height of 1 m. Figure 3 compares the time-history curve of the impact force obtained by the laboratory test under multiple impacts with the curve that was extracted by the numerical model.

Numerical Model Validation
The sphere rockfall with a mass of 1.7 kg impacts the sand cushion multiple times at the falling height of 1m. Figure 3 compares the time-history curve of the impact force obtained by the laboratory test under multiple impacts with the curve that was extracted by the numerical model. The instant when the rockfall starts to contact the cushion is set to zero. The impact force increases rapidly when the rockfall contacts the sand buffer layer and quickly decreases to zero after reaching the maximum impact force. The entire shock process is very brief. Under the four impacts, the numerical simulation results obtained by using discrete element software are close to the experimental results, and the error is small. In general, the average value of the maximum impact force results of the three groups of laboratory tests was compared with the numerical results. The two results differ in maximum impact force by 10.76%, 7.04%, 2.79%, and 5.59%, respectively. In terms of the duration of the impact process and the change trend of the impact force, the simulation results agree well with the test results.   The instant when the rockfall starts to contact the cushion is set to zero. The impact force increases rapidly when the rockfall contacts the sand buffer layer and quickly decreases to zero after reaching the maximum impact force. The entire shock process is very brief. Under the four impacts, the numerical simulation results obtained by using discrete element software are close to the experimental results, and the error is small. In general, the average value of the maximum impact force results of the three groups of laboratory tests was compared with the numerical results. The two results differ in maximum impact force by 10.76%, 7.04%, 2.79%, and 5.59%, respectively. In terms of the duration of the impact process and the change trend of the impact force, the simulation results agree well with the test results.
Therefore, the analytical model and method proposed above can more reliably reflect the dynamic response of rockfalls impacting a sand buffer layer. This is applied to the study in the next section considering the multiple impacts of rockfall on a sand buffer layer at the engineering scale.

Results and Discussion
Now, consider the dynamic process of rockfall impacting a sand buffer layer multiple times at the engineering scale. The impact energy of rockfall disasters along Japanese railways was calculated. Rockfall energies are approximately normally distributed. The proportion of rockfall events with an impact energy of less than 100 kJ was 68%, and the proportion of rockfall events with an impact energy of less than 1000 kJ reached 90% [24]. A study conducted off the east coast of New Wales, Australia, found that the average diameter of rockfalls in the sandstone and granite areas of the basin was 0.45 m, and 95% of the rockfalls had an impact energy of less than 1340 kJ [25]. Therefore, in the research of this paper, the diameter of the rockfall is 0.5 m, the density is consistent with that in the laboratory experiment, and the mass of the rockfall is 1.144 t. The drop heights are 10 m, 20 m, 30 m, 40 m, and 50 m, respectively. The impact energy corresponds to 112 kJ, 224 kJ, 336 kJ, 448 kJ, and 560 kJ, respectively.

Large-Scale Impact Numerical Model
In the research for this paper, the influence of the constraints around the sand buffer layer on the impact effect is ignored, and the size of the buffer layer is guaranteed to be large enough. When the ratio of the size of the buffer layer to the diameter of the rockfall exceeds five, the effect of the lateral boundary constraints of the buffer layer can be ignored [26]. Therefore, the size of the cushion is selected as 5 m × 5 m × 2 m. Since the size of the buffer layer is much larger than that of the laboratory test in Section 2, the soil radius has also been adjusted accordingly, ranging from 0.03 m to 0.06 m. Except for the adjustment of the soil radius in the cushion layer, the parameters of the cushion layer do not change. The number of spherical particles in the buffer layer in the numerical model is 76801. The model diagram is shown in Figure 4. the dynamic response of rockfalls impacting a sand buffer layer. This is app study in the next section considering the multiple impacts of rockfall on a s layer at the engineering scale.

Results and Discussion
Now, consider the dynamic process of rockfall impacting a sand buffer lay times at the engineering scale. The impact energy of rockfall disasters along Jap ways was calculated. Rockfall energies are approximately normally distribute portion of rockfall events with an impact energy of less than 100 kJ was 68 proportion of rockfall events with an impact energy of less than 1000 kJ reache A study conducted off the east coast of New Wales, Australia, found that t diameter of rockfalls in the sandstone and granite areas of the basin was 0.45 m of the rockfalls had an impact energy of less than 1340 kJ [25]. Therefore, in th of this paper, the diameter of the rockfall is 0.5 m, the density is consistent with laboratory experiment, and the mass of the rockfall is 1.144 t. The drop height 20 m, 30 m, 40 m, and 50 m, respectively. The impact energy corresponds to 11 336 kJ, 448 kJ, and 560 kJ, respectively.

Large-Scale Impact Numerical Model
In the research for this paper, the influence of the constraints around the s layer on the impact effect is ignored, and the size of the buffer layer is guara large enough. When the ratio of the size of the buffer layer to the diameter of exceeds five, the effect of the lateral boundary constraints of the buffer layer nored [26]. Therefore, the size of the cushion is selected as 5 m × 5 m × 2 m. Sin of the buffer layer is much larger than that of the laboratory test in Section 2, the has also been adjusted accordingly, ranging from 0.03 m to 0.06 m. Except for ment of the soil radius in the cushion layer, the parameters of the cushion la change. The number of spherical particles in the buffer layer in the numeric 76801. The model diagram is shown in Figure 4. Next, we discuss the dynamic response of rockfall on the sand cushion f impacts. The analysis focuses on the impact force and penetration depth on the ion. The purpose is to reveal the importance of considering the multiple impa falls.  Next, we discuss the dynamic response of rockfall on the sand cushion for multiple impacts. The analysis focuses on the impact force and penetration depth on the sand cushion. The purpose is to reveal the importance of considering the multiple impacts of rockfalls.  Figure 5 shows the cloud map of the change in the position of the rockfall and the cushion after each impact when the rockfall falls from the same height and impacts the sand cushion several times, for a total of four impacts. To observe the changes more intuitively in the position of the rockfall and the impact point on the top surface of the cushion after the impact of the rockfall, the model was cut along the xOz plane to observe the cross-section. The upper black line is the starting surface of the sand cushion, and the lower black line is the plane of the maximum penetration depth after impact. The vertical distance between the two black lines is the maximum penetration depth. The area enclosed by the white line and the upper black line is the range of particles affected by the impact. The impact of the rockfall caused a certain penetration into the top surface of the cushion, and some particles near the impact point were lifted. With the increase in the number of impacts, the height of the rockfall relative to the top surface of the cushion gradually decreased, and the affected particles in the lower part of the impact point in the cushion and near the impact point on the top surface of the cushion tended to increase.  Figure 6 shows the time-history curve of the impact force of the spherical rockfall impacting the sand cushion multiple times from different heights. The instant when the falling rock and the cushion first contacted was taken as time zero. The change rule of the rockfall impact force is that at the moment of contact with the cushion, the impact force rapidly increases to the peak impact force and then decreases relatively slowly to zero. Since each impact will cause a certain penetration into the cushion, the time of contact with the cushion will be delayed for the next impact compared with the previous impact. With the increase in the number of impacts, the time lag of the peak impact force and the magnitude of the impact force increase, but the duration of the entire impact process is gradually shortened. The maximum impact force and impact duration are important factors to be considered in actual rockfall-protection projects. It can be seen that the impact number will affect the impact response of rockfall. Therefore, it is obvious that when discussing the dynamic response to rockfall impact, it is not comprehensive enough to consider the consequences of only one rockfall impact.  Figure 6 shows the time-history curve of the impact force of the spherical rockfall impacting the sand cushion multiple times from different heights. The instant when the falling rock and the cushion first contacted was taken as time zero. The change rule of the rockfall impact force is that at the moment of contact with the cushion, the impact force rapidly increases to the peak impact force and then decreases relatively slowly to zero. Since each impact will cause a certain penetration into the cushion, the time of contact with the cushion will be delayed for the next impact compared with the previous impact. With the increase in the number of impacts, the time lag of the peak impact force and the magnitude of the impact force increase, but the duration of the entire impact process is gradually shortened. The maximum impact force and impact duration are important factors to be considered in actual rockfall-protection projects. It can be seen that the impact number will affect the impact response of rockfall. Therefore, it is obvious that when discussing the dynamic response to rockfall impact, it is not comprehensive enough to consider the consequences of only one rockfall impact.   Figure 7 shows the results of the maximum impact force of the spherical rock falling on the sand cushion multiple times from five falling heights. The impact force increases with the drop height. Regardless of the drop height, the impact force results of the four impacts show a gradually increasing trend, but the increasing range gradually decreases. For example, the impact force results for the drop height of 50 m are 2299 kN, 2645 kN, 2807 kN, and 2819 kN for the four impacts. The results of the second, third, and fourth impacts were increased by 15%, 22.1%, and 22.6%, respectively, compared with the first impact. The fourth impact force is slightly larger than the third impact force. Therefore, only four shocks are considered in this paper to study multiple shocks. For different drop heights, the total increase in the impact force of the four impacts is different. For heights of 10 m, 20 m, 30 m, 40 m, and 50 m, the fourth impact force increased by 33.9%, 32.7%, 27.4%, 23.1%, and 22.6%, respectively, compared with the first impact force. Because of the different falling heights, the impact energy of the falling rock varies greatly. With the increase in the number of impacts, the main reason for the gradual increase in the impact force is the gradual increase in the compactness of the cushion. However, the increase in the compactness of the cushion has a limit. After the first impact, the compactness of the cushion increases more at higher drop heights. During the last few impacts, the change in cushion disorientation decreases with an increasing drop height. Therefore, a phenomenon is observed in which an increasing falling height results in an impact force that is smaller under multiple impacts than under a single impact. However, since the maximum impact force increases with the increasing impact height, a higher impact height must still be considered.

Impact Force at the Center of Multiple Impacts
stainability 2023, 15, x FOR PEER REVIEW of 10 m, 20 m, 30 m, 40 m, and 50 m, the fourth impact force increased by 3 27.4%, 23.1%, and 22.6%, respectively, compared with the first impact forc the different falling heights, the impact energy of the falling rock varies grea increase in the number of impacts, the main reason for the gradual increase force is the gradual increase in the compactness of the cushion. However, th the compactness of the cushion has a limit. After the first impact, the compa cushion increases more at higher drop heights. During the last few impacts, t cushion disorientation decreases with an increasing drop height. Therefore non is observed in which an increasing falling height results in an impact smaller under multiple impacts than under a single impact. However, since t impact force increases with the increasing impact height, a higher impact hei be considered.  The impact force of rockfall on the shed cavity cushion is a crucial inspection index in the design of shed cavity structures. No existing impact force algorithm considers the calculation of the impact force during multiple impacts. The maximum impact force algorithm of rockfalls proposed by Labiouse et al. [27] based on Hertzian contact theory is modified by introducing a coefficient, α, related to the number of impacts, t, through the numerical results of the maximum impact force under multiple impacts.
The maximum impact force equation can be expressed as the following: In the formula, P max is the maximum impact force (kN); M E is the modulus of the subgrade reaction obtained from a standardized plate bearing test on the soil cushion (kN/m 2 ); R is the radius of the falling block in contact with the cushion (m); W is the weight of the falling rock (kN); and H is the falling height (m).  Comparing the two curves in Figure 8 shows that the maximum impac lated by the formula proposed by Labiouse et al. agrees well with the impac of the numerical model in this paper. The error of the calculated result can within 10%. Therefore, this formula can be combined with the numerical resu the calculation of the maximum impact force when the rockfall hits the buffe times. The numerical results for the impact force of the falling rock impacting t ion multiple times from heights of 10 m, 20 m, 30 m, 40 m, and 50 m are a discussed. Consider the ratio of the impact force after the second, third, an pacts to the impact force of the first impact for a given height. Considering heights comprehensively, the expression for α corresponding to the numbe is obtained by performing a fit. Figure 9 compares the numerical results o force with the results of the fitting formula for different drop heights. Comparing the two curves in Figure 8 shows that the maximum impact force calculated by the formula proposed by Labiouse et al. agrees well with the impact force result of the numerical model in this paper. The error of the calculated result can be controlled within 10%. Therefore, this formula can be combined with the numerical results to explore the calculation of the maximum impact force when the rockfall hits the buffer layer many times.
The numerical results for the impact force of the falling rock impacting the sand cushion multiple times from heights of 10 m, 20 m, 30 m, 40 m, and 50 m are analyzed and discussed. Consider the ratio of the impact force after the second, third, and fourth impacts to the impact force of the first impact for a given height. Considering the five drop heights comprehensively, the expression for α corresponding to the number of impacts t is obtained by performing a fit. Figure 9 compares the numerical results of the impact force with the results of the fitting formula for different drop heights.
where t is the number of impacts and α is an expression that only relates to the number of impacts.
ainability 2023, 15, x FOR PEER REVIEW Figure 9. Comparison of impact force results between numerical and fitted for Figure 10 compares the calculated results of the improved maxi algorithm with the numerical simulation results. Impact force algorithm for multiple impacts: The numerical results for the impact force of the falling rock impacting the sand cushion multiple times from heights of 10 m, 20 m, 30 m, 40 m, and 50 m are analyzed and discussed. Consider the ratio of the impact force after the second, third, and fourth impacts to the impact force of the first impact for a given height. Considering the five drop heights comprehensively, the expression for α corresponding to the number of impacts t is obtained by performing a fit. Figure 9 compares the numerical results of the impact force with the results of the fitting formula for different drop heights. Figure 10 compares the calculated results of the improved maximum impact force algorithm with the numerical simulation results.   The impact forces of the impacts in the figure correspond to drop heights of 10 m, 30 m, 40 m, and 50 m, respectively, from left to right. Solid color, unpatterne represent the numerical results, and diagonal striped fills represent the calculation of formulas. The results calculated by the improved impact force algorithm agre with the numerical simulation results as a whole, and the maximum error of the calc results is 5.87%. For the more dangerous, high-impact energy cases, the error is con within 5%. Therefore, the maximum impact force algorithm modified by fitting t merical simulation results in this paper can be considered reasonable.

Penetration Depth at the Center of Multiple Impacts
In addition to the maximum impact force of the rockfall, the maximum pene depth of the rockfall into the buffer layer is also an important parameter to be cons in engineering design [20,28]. During the experiment, because of the special back phenomenon of the sand cushion, the depth of the cushion after the rockfall im difficult to obtain. However, the change in the penetration depth can be ob The impact forces of the impacts in the figure correspond to drop heights of 10 m, 20 m, 30 m, 40 m, and 50 m, respectively, from left to right. Solid color, unpatterned fills represent the numerical results, and diagonal striped fills represent the calculation results of formulas. The results calculated by the improved impact force algorithm agree well with the numerical simulation results as a whole, and the maximum error of the calculated results is 5.87%. For the more dangerous, high-impact energy cases, the error is controlled within 5%. Therefore, the maximum impact force algorithm modified by fitting the numerical simulation results in this paper can be considered reasonable.

Penetration Depth at the Center of Multiple Impacts
In addition to the maximum impact force of the rockfall, the maximum penetration depth of the rockfall into the buffer layer is also an important parameter to be considered in engineering design [20,28]. During the experiment, because of the special backfilling phenomenon of the sand cushion, the depth of the cushion after the rockfall impact is difficult to obtain. However, the change in the penetration depth can be observed intuitively in the numerical simulation. Figure 11 shows the time-history curve of the penetration depth of the spherical rockfall impacting the sand cushion multiple times from different heights. After the rockfall is in contact with the cushion, the penetration depth gradually increases. After reaching the peak penetration depth, the penetration depth recovers to a certain extent because of the rebound phenomenon of the rockfall. Similarly, with an increasing number of impacts, the maximum penetration depth appears slightly later and gradually increases. This result confirms that the size of the maximum penetration depth can provide a basis for selecting the thickness of the buffer layer of the shed cavity structure. It also confirms the necessity of considering multiple impacts of rockfalls. depth recovers to a certain extent because of the rebound phenomenon of the rockfall. Similarly, with an increasing number of impacts, the maximum penetration depth appears slightly later and gradually increases. This result confirms that the size of the maximum penetration depth can provide a basis for selecting the thickness of the buffer layer of the shed cavity structure. It also confirms the necessity of considering multiple impacts of rockfalls.  Figure 12 shows the results of the maximum penetration depth of the spherical rockfall with multiple impacts on the sand cushion for five drop heights. As with the maximum impact force results, the maximum penetration depth increases with the number of impacts regardless of the drop height. However, unlike the results for the maximum impact force, the maximum penetration depth increases with the same magnitude as the number of impacts increases. For example, when the drop height is 10 m, the maximum penetration depths are 0.181 m, 0.236 m, 0.293 m, and 0.343 m, respectively. Compared with the first impact, the maximum penetration depth of the latter three increased by 30.4%, 61.8%, and 89%, respectively. At different drop heights, the total increase in the penetration depth under the four impacts is roughly the same. However, the penetration depth increases slightly with an increasing drop height. At heights of 10 m, 20 m, 30 m, 40 m, and 50 m, the ratios of the fourth penetration depth to the first penetration depth were 1.89, 1.97, 2.02, 2.06, and 2.08, respectively. Figure 12 shows the results of the maximum penetration depth of the sph fall with multiple impacts on the sand cushion for five drop heights. As wit mum impact force results, the maximum penetration depth increases with the impacts regardless of the drop height. However, unlike the results for the ma pact force, the maximum penetration depth increases with the same magni number of impacts increases. For example, when the drop height is 10 m, the penetration depths are 0.181 m, 0.236 m, 0.293 m, and 0.343 m, respectively. with the first impact, the maximum penetration depth of the latter three in 30.4%, 61.8%, and 89%, respectively. At different drop heights, the total inc penetration depth under the four impacts is roughly the same. However, the depth increases slightly with an increasing drop height. At heights of 10 m, 20 m, and 50 m, the ratios of the fourth penetration depth to the first penetration 1.89, 1.97, 2.02, 2.06, and 2.08, respectively.

Distance between Different Impact Positions
At present, the only research on multiple impacts of rockfalls is limited to multiple impacts of rockfalls on the same location. Notably, when a rockfall curs, the impact position of the multiple impacts of the rockfall on the cu change. After the impact position changes, it needs to be considered whether t response between two adjacent shocks is affected by the distance between the tions. To address this problem, this paper studies the effect of the distance b impact points on the impact response when the rockfall impacts different i tions. The approximate distance that the second impact is not affected by the when two consecutive impacts are at different positions is given.
Considering the rockfall radius R and the size of the cushion, in addition impact at the center of the cushion, six impact positions are also considered. Th 1R, 1.5R, 2R, 3R, and 4R from the center point O, corresponding to six points A

Distance between Different Impact Positions
At present, the only research on multiple impacts of rockfalls is limited to consecutive multiple impacts of rockfalls on the same location. Notably, when a rockfall disaster occurs, the impact position of the multiple impacts of the rockfall on the cushion may change. After the impact position changes, it needs to be considered whether the dynamic response between two adjacent shocks is affected by the distance between the shock positions. To address this problem, this paper studies the effect of the distance between the impact points on the impact response when the rockfall impacts different impact locations. The approximate distance that the second impact is not affected by the first impact when two consecutive impacts are at different positions is given.
Considering the rockfall radius R and the size of the cushion, in addition to the first impact at the center of the cushion, six impact positions are also considered. They are 0.5R, 1R, 1.5R, 2R, 3R, and 4R from the center point O, corresponding to six points A, B, C, D, E, and F, respectively, as shown in Figure 13. To compare the effect of the distance between shock locations on the shock response, two shock simulations were performed. Taking point A as an example, the other five points are the same. For the convenience of comparison, the first one impacts point A only once. The second carries out the simulation of rockfall impacting the sand cushion twice: the first impact is at point O, the center point of the cushion, and the second impact is at point A. The numerical results of the two shocks were compared and analyzed, and the distance between the two shock locations that did not affect each other was explored.

Impact Force at Different Positions
Taking the falling height of 50 m as an example, Figure 14 shows the time-history curve of the impact force of the rockfall at different positions on the top surface of the sand cushion. The left picture is the curve of impacting the selected impact position only once, and the right picture is the time-history curve of impacting the center point of the cushion and then impacting the selected position. The time-history curve of the impact force at each point in the left figure shows the same change trend. The instantaneous impact force of the rockfall contacting the cushion rapidly increases to the maximum impact force and then gradually decreases to zero. The time-history curves of different shock points in the right figure are quite different, especially point A. The main reason for this is that after the first impact on point O, the particles of the cushion layer are partially arched, as shown in Figure 5, the sand particles here are relatively loose, and point A is located in the arched part of the particles. Therefore, the maximum impact force is delayed. The impact force time-history curves of other impact locations are consistently in trend with those with only one impact, and the magnitude of the peak impact force is slightly different. To compare the effect of the distance between shock locations on the shock response, two shock simulations were performed. Taking point A as an example, the other five points are the same. For the convenience of comparison, the first one impacts point A only once. The second carries out the simulation of rockfall impacting the sand cushion twice: the first impact is at point O, the center point of the cushion, and the second impact is at point A. The numerical results of the two shocks were compared and analyzed, and the distance between the two shock locations that did not affect each other was explored.

Impact Force at Different Positions
Taking the falling height of 50 m as an example, Figure 14 shows the time-history curve of the impact force of the rockfall at different positions on the top surface of the sand cushion. The left picture is the curve of impacting the selected impact position only once, and the right picture is the time-history curve of impacting the center point of the cushion and then impacting the selected position. The time-history curve of the impact force at each point in the left figure shows the same change trend. The instantaneous impact force of the rockfall contacting the cushion rapidly increases to the maximum impact force and then gradually decreases to zero. The time-history curves of different shock points in the right figure are quite different, especially point A. The main reason for this is that after the first impact on point O, the particles of the cushion layer are partially arched, as shown in Figure 5, the sand particles here are relatively loose, and point A is located in the arched part of the particles. Therefore, the maximum impact force is delayed. The impact force time-history curves of other impact locations are consistently in trend with those with only one impact, and the magnitude of the peak impact force is slightly different. To compare the effect of the distance between shock locations on the shock response, two shock simulations were performed. Taking point A as an example, the other five points are the same. For the convenience of comparison, the first one impacts point A only once. The second carries out the simulation of rockfall impacting the sand cushion twice: the first impact is at point O, the center point of the cushion, and the second impact is at point A. The numerical results of the two shocks were compared and analyzed, and the distance between the two shock locations that did not affect each other was explored.

Impact Force at Different Positions
Taking the falling height of 50 m as an example, Figure 14 shows the time-history curve of the impact force of the rockfall at different positions on the top surface of the sand cushion. The left picture is the curve of impacting the selected impact position only once, and the right picture is the time-history curve of impacting the center point of the cushion and then impacting the selected position. The time-history curve of the impact force at each point in the left figure shows the same change trend. The instantaneous impact force of the rockfall contacting the cushion rapidly increases to the maximum impact force and then gradually decreases to zero. The time-history curves of different shock points in the right figure are quite different, especially point A. The main reason for this is that after the first impact on point O, the particles of the cushion layer are partially arched, as shown in Figure 5, the sand particles here are relatively loose, and point A is located in the arched part of the particles. Therefore, the maximum impact force is delayed. The impact force time-history curves of other impact locations are consistently in trend with those with only one impact, and the magnitude of the peak impact force is slightly different.   to right. A solid, unpatterned fill represents the first impact simulation, and a horizontal stripe fill represents the second impact simulation. Since the particles in the mat are spheres, this results in slightly different levels of compaction at each location in the mat. At the same time, as the impact position approaches the cushion boundary, boundary constraints will also have a certain impact. Therefore, in the first impact simulation, the impact force slightly differs between impact locations. Taking the drop height of 50 m as an example, the impact force at point A to point F is 2306 kN, 2327 kN, 2297 kN, 2356 kN,  2287 kN, and 2238 kN, respectively. The maximum error is 5%. This error can be cancelled out when comparing two shock results. This method can still be considered reasonably reliable for investigating the effect of the distance between different impact locations. and 50 m from left to right. A solid, unpatterned fill represents the first impact simulat and a horizontal stripe fill represents the second impact simulation. Since the particle the mat are spheres, this results in slightly different levels of compaction at each locat in the mat. At the same time, as the impact position approaches the cushion bounda boundary constraints will also have a certain impact. Therefore, in the first impact sim lation, the impact force slightly differs between impact locations. Taking the drop hei of 50 m as an example, the impact force at point A to point F is 2306 kN, 2327 kN, 2  kN, 2356 kN, 2287 kN, and 2238 kN, respectively. The maximum error is 5%. This er can be cancelled out when comparing two shock results. This method can still be con ered reasonably reliable for investigating the effect of the distance between different pact locations. The abscissa in Figure 15 is from point A to point F from left to right, and the dista from point O increases gradually. That is, the distance between the impact positions creases gradually during multiple impacts. First, it can be observed that at different d heights and different impact positions, an impact is smaller if the center point of the cu ion is impacted first. This is mainly because after the first impact on the center point, particles near the impact point are smashed, while the particles at other positions beco sparser than the initial state. Therefore, the impact force is slightly larger when the imp occurs only once. Taking the drop height of 50 m as an example, the difference betw the impact forces without and with a prior impact on the center point is 3.45%, 8.4 7.65%, 1.53%, 0.04%, and 0.08%. Judging from the difference in the impact force, poin which is a radius from point O, is the most affected point, followed by point C, poin and point D. Starting from point E at 1.5 m from point O, that is, when the distance fr the impact position is 3R, the impact force results in consistency between the two imp modes. This finding means that when the distance between the impact positions of t adjacent impacts is greater than 3R, the maximum impact force during the second imp is not affected by the first impact. The abscissa in Figure 15 is from point A to point F from left to right, and the distance from point O increases gradually. That is, the distance between the impact positions increases gradually during multiple impacts. First, it can be observed that at different drop heights and different impact positions, an impact is smaller if the center point of the cushion is impacted first. This is mainly because after the first impact on the center point, the particles near the impact point are smashed, while the particles at other positions become sparser than the initial state. Therefore, the impact force is slightly larger when the impact occurs only once. Taking the drop height of 50 m as an example, the difference between the impact forces without and with a prior impact on the center point is 3.45%, 8.46%, 7.65%, 1.53%, 0.04%, and 0.08%. Judging from the difference in the impact force, point B, which is a radius from point O, is the most affected point, followed by point C, point A, and point D. Starting from point E at 1.5 m from point O, that is, when the distance from the impact position is 3R, the impact force results in consistency between the two impact modes. This finding means that when the distance between the impact positions of two adjacent impacts is greater than 3R, the maximum impact force during the second impact is not affected by the first impact. The same conclusion is obtained for drop heights of 10 m, 20 m, 30 m, and 40 m.

Penetration Depth at Different Positions
Taking the falling height of 50 m as an example, Figure 16 shows the time-history curve of the penetration depth of the rockfall at different positions on the top surface of the sand cushion. This curve is the same as the time-history curve of the impact force. The left picture shows the curve of only one impact on the selected impact position, and the right picture shows the time-history curve of the impact on the center point of the cushion and then the selected position. The time-history curve of the penetration depth of each point in the left figure shows the same trend of change. After the rockfall contacts the cushion, the penetration depth gradually increases to the maximum penetration depth and then recovers to a certain extent. The time-history curves of different shock points in the right figure are quite different, particularly point A. The penetration depth is substantially different from other impact locations at point A. The maximum penetration depth is quite different from the result when only point A is impacted, as shown in the left picture, and is much larger than the penetration depth of other impact points, as shown in the right picture. At other impact locations, the time-history curves of the penetration depth are consistently in trend with those with only one impact, and the magnitude of the peak penetration depth is slightly different. In order to achieve more intuitive understanding, the vertical displacement nephogram of particles at different impact positions is shown in Figure 17. and is much larger than the penetration depth of other impact points, as shown in the right picture. At other impact locations, the time-history curves of the penetration depth are consistently in trend with those with only one impact, and the magnitude of the peak penetration depth is slightly different. In order to achieve more intuitive understanding, the vertical displacement nephogram of particles at different impact positions is shown in Figure 17.    The abscissa in Figure 18 is from point A to point F from left to right, and the distance from point O increases gradually; that is, the distance between the impact positions increases gradually during multiple impacts. Taking the drop height of 50 m as an example, the difference between the penetration depth results without and with a prior impact at the center point is 22.86%, 0.40%, 7.44%, 4.53%, 0.60%, and 0.15%, respectively. From the difference in the penetration depths, point A at the shortest distance from point O is the most affected point, followed by point C. After the first impact on point O, the sand particles near the impact point are relatively loose. Point A is located in the loose part of the particles, so the phenomenon of the maximum depression depth at point A changes greatly. Likewise, starting from point E at 1.5 m from point O, that is, at a distance of 3R from the impact location, the difference between the penetration depth results for the two impact modes is very small. This observation shows that when the distance between the impact positions of two adjacent impacts is greater than 3R, the maximum penetration depth of the second impact is slightly affected by the first impact. The same conclusion is obtained for drop heights of 10 m, 20 m, 30 m, and 40 m. tion depth results when the drop height is 50 m as an example, the penetration de from point A to point F are 0.2903 m, 0.2925 m, 0.2918 m, 0.2863 m, 0.2875 m, and 0. m, respectively. The maximum error is 1.9%. The distance between the impact loca where the impact effect is affected between two consecutive impacts can still be expl by comparing the maximum penetration depths of the two impact modes. The abscissa in Figure 18 is from point A to point F from left to right, and the dist from point O increases gradually; that is, the distance between the impact posi

Conclusions
This paper focuses on the dynamic response of spherical rockfalls' multiple impacts on a sand cushion from different falling heights. First of all, carry out the reduced-scale impact test. Secondly, a numerical model is established and compared with the test results to verify the reliability of the model. On this basis, an engineering scale numerical model is established to discuss the effect of multiple impacts on the impact response, including multiple impacts at the same location and multiple impacts at different locations. In addition, based on the numerical simulation results of the maximum impact force, the existing impact force algorithm is modified to provide a reference basis for the design of rockfall prevention engineering. The specific conclusions are as follows: (1) When rockfall disasters occur, there are often multiple impacts. Therefore, it is necessary to consider the multiple impacts of rockfall. Both the test and numerical simulation show that the maximum impact force and the maximum penetration depth increase to varying degrees with the increase in the number of impacts. For a given drop height, the increase in the maximum impact force decreases gradually, and the increase in the maximum penetration depth is almost constant. With the increase in the impact height, the increase in the impact force decreases and the increase in the penetration depth slightly increases. (2) Based on the numerical simulation results, the maximum impact force equation proposed by Labiouse et al. was revised, and a coefficient related to the number of impacts was introduced to obtain the impact force after each impact under multiple impacts. It can provide a more direct reference for the design of shed tunnel structures. (3) In the rockfall disaster, there is the possibility that the two adjacent impact positions are different during multiple impacts. Therefore, with the help of numerical simulation, this paper explores the impact of the distance between impact positions on the response to the impact. The numerical results of the maximum impact force and the maximum penetration depth show that the impact response of the second impact is not affected by the first impact when the distance between the impact positions of two adjacent impacts is three times the radius of the falling rock. Under this condition, the impact response of rockfall can be ignored for the first impact and only the second impact is considered. Data Availability Statement: This study did not report any data.