Numerical Simulation of the Shed-Tunnel Structure Dynamic Response under the Rockfall Repeated Impacts with FEM-SPH Method

1. Introduction

Rockfall in the mountainous area has been the third-largest geological disaster after landslides and debris flow [1,2]. Because of its high frequency of occurrence and strong randomness of the movement characteristics, the accurate predictions of the rockfall occurrence time, movement trajectory and its damage to the protective structure are very difficult [3,4,5,6]. Therefore, the protective structures such as the rockfall walls, flexible nets and shed-tunnels are often used in the practical engineering to protect the roads and buildings [7,8]. Among them, the shed-tunnel structure has been widely used because of its good rockfall impact resistance, small damage to the natural mountains and simple maintenance [9]. Therefore, the dynamic response of the shed-tunnel structure under rockfall impact has attracted extensive attention from the scholars at

home and abroad, and has been deeply studied from three aspects, e.g. the model and field test, theoretical analysis and numerical simulation.

In the early time, the researches mainly focus on the common shed-tunnel, which consists of two parts: the reinforced concrete main structures and the cushion material on its top, and the latter plays a vital role in mitigating the rockfall impact force [10,11]. Therefore, most researches are concentrated on the effect of the cushion material property such as its type and depth on the dynamic response of the shed-tunnel. Kawahara and Muro [12] studied the effect of dry density and thickness of sand cushion on the rockfall impact pressure, and found that it increased with increasing the dry density and rapidly decreased with the increase in the soil thickness. With the test and numerical simulations, Sun et al. [13] found that the tire cushion layer could effectively reduce peak acceleration and the maximum impact force. Filling sand and gravel in the tire could improve tire stiffness and energy absorption capacity but will decrease cushion effect due to its large density. Zhao et al. [15] proposed a new cushion layer composed of a foam material (expanded polystyrene (EPS) or expandable polyethylene (EPE)) and granular material (soil or sand), and then the large-scale rockfall impact experiments are carried out on a reinforced concrete shed-tunnel with three types of cushion: sand, sand-EPS, and sand-EPE. It is found that sand-EPE is the best to resist rockfall impact and protect the shed-tunnel structure. With the model test and numerical simulation, Shen et al. [16] studied the effect of the soil particle size on the cushion buffer efficient of the rockfall impact force. Wu et al. [17] conducted the rockfall impact test to study the relationship between the impact force of the reinforced concrete slab induced by rockfall and the cushion thickness. In order to improve the ductility and the impact resistance ability of the shed-tunnel, Zhao et al. [14] proposed two types of geofoam, e.g. EPE and EPS, to replace part of sand, forming a new composite cushion. The tests showed that EPE geofoam had better resilience than EPS one, and it was more suitable for improving the capability of the composite cushion to resist multi-impacts and protect the shed-tunnel under multiple rockfall impacts. Zhao et al. [18] conducted a series of laboratory tests and numerical investigations on buffer performance of geofoam under the impact of different rock shapes, and found that the geofoam thickness and the rockfall shape both had an obvious effect on the maximum deformation and vertical stress in the geofoam. Moreover, the rockfall shape plays a crucial role in evaluating the buffer performance of the composite cushion. However, because of the construction difficulty and high cost in poor-quality bedrock, the reinforced concrete is gradually replaced by a new type of the flexible shed-tunnel, which consists of the flexible nets and vaulted steel structures. Its main advantages are simple to build and emergent construction [19]. Shi et al. [20] proposed a new concept of the flexible shed-tunnel and conducted an impact experiment on the 1:1 prototype to evaluate its performance. Wang et al. [21] conducted the full-scale tests on a flexible three-module shed-tunnel with a rockfall kinetic energy ranging from 100 to 250kJ to study its mechanical characteristics, and put forward the useful suggestions for optimal design.

Next, the theoretical studies on the dynamic response of the shed-tunnel under rockfall impact are also made. Olsson et al. [22] proposed a dynamic response model for the anisotropic composite plate under the small mass sphere impact based on Hertz contact theory. On basis of it, Wang et al. [23] studied the elastoplastic dynamic response mechanism of the reinforced concrete shed-tunnel under the rockfall. In view of the viscoelastic contact theory, Chen et al. [24] proposed a rockfall impact mechanics model, and the general solution form of the time-history law of rockfall impact load is obtained.

Finally, many researchers study the dynamic response of the shed-tunnel with the numerical simulation. Shen et al. [25] studied the effect of the rockfall shape on the response of rockfall impact against a soil buffer layer with discrete element method, and found that the penetration depth decreased linearly with increasing the sphericity degree of the rockfall. Ouyang et al. [26] proposed a new shed-tunnel absorption cushion composed of geogrid reinforced soil and EPS, and studied its dynamic mechanical behavior under rockfall impact with finite element method. With the SPH-FEM coupling method, Zhong et al. [27] studied the shed-tunnel mechanical behavior under rockfall impacts from three aspects, e.g. the rockfall impact force, the inertia effect coefficient and the structural damage evaluation. Wang et al. [28] conducted the peridynamics simulation of structural damage characteristics in the shed-tunnel under rockfall impact, and found that the non-circular rockfall would cause larger damage to it.

However, the existing studies mainly focus on the rockfall single impact on the tunnel structure, and few scholars consider the cumulative damage caused by the rockfall repeated impacts on the shed-tunnel structure. In the practical engineering, the shed-tunnel structure has a long service life, so it will inevitably encounter the rockfall repeated impacts. Therefore, with ANSYS/LS-DYAN code, the FEM-SPH coupling numerical model is established and the full restart technology of ANSYS/LS-DYAN code is adopted to conduct a more systematic and in-depth simulation study on the dynamic response of the shed-tunnel structure under the rockfall repeated impacts. The effects of four factors, e.g. rockfall mass, rockfall impact velocity, rockfall shape and rockfall impact angle, on the dynamic response of the shed-tunnel structure under the rockfall repeated impacts are studied.

2. Establishment of the numerical calculation model

2.1. The FEM-SPH coupling numerical model

The numerical simulation is conducted according to the model test conducted by Wu et al. [17]. The size of the concrete support and the reinforced concrete plate is 0.2m×1.6m×0.25m and 2.4m×1.6m×0.25m respectively. Two layers of φ14@200mm rebar mesh are laid vertically and orthogonal inside the reinforced concrete plate, whose protective layer is 20mm thick. The buffer layer is 2.4m×1.6m×0.3m. The bottom of the support is fixed in all directions and four sides of the buffer layer are fixed in normal direction. The numerical model is shown in Figure 1.

Although it has high computational efficiency in calculating the structure dynamic response, FEM will have errors because of the occurrence of the distorted meshes in the sand or soil buffer layer induced by the rockfall repeated impacts. Therefore, early in 1977, Lucy [29] and Gingold et al. [30] respectively proposed a Lagrangian meshless particle method, namely SPH, which can effectively avoid the mesh distortion caused by material local large deformation. However, its computational efficiency is much lower than that of FEM. Therefore, in order to make full use of the advantages of these two methods, many researchers proposed the FEM-SPH coupling method, which has been successfully used in the field of material failure simulation such as blasting and impact in recent years [27,31,32]. Therefore, the FEM-SPH coupling method is adopted here, in which the concrete supports and reinforced concrete plate are simulated with FEM, and the buffer layer is simulated with SPH to avoid the mesh distortion. The algorithmic principle of the FEM-SPH coupling method has been described in detail in many literatures [33] and we do not state it again.

In conclusion, the numerical model of rockfall impact the shed-tunnel structure is established with the FEM-SPH coupling method, which can not only effectively avoid the mesh distortion, but also improve the reliability and calculation efficiency of the simulation.

2.2. The calculation parameters of the materials

In the numerical model shown in Figure 1, the materials involved are rockfall, buffer layer and reinforced concrete, where the reinforced concrete is simulated with the separate model, namely, it is regarded as a combination of concrete and rebar. In view of the existing constitutive models in ANSYS/LS-DYAN code and the relevant researches [27], the constitutive models for these above four materials are as follows. Rockfall, buffer layer, concrete and rebar are assumed to be the rigid model, MAT_FHWA_SOIL model, MAT_CSCM model and PLASTIC_KINEMATIC model, respectively. Their corresponding calculation parameters are shown in Table 1.

2.3. Verification of the numerical model

In order to valid the numerical model, the rockfall impact force on the buffer layer is firstly studied. Assume the rockfall vertically impacts the buffer layer, and three impacts are carried out with the falling heights of 5, 6 and 7m respectively. Accordingly, the rockfall impact velocity when reaching the buffer layer is 9.9m/s, 10.84m/s and 11.71m/s respectively. The simulation results of the rockfall impact force on the buffer layer are compared with those tested by Wu et al. [17] and those obtained with the often-used calculation methods shown in Tab.2 to valid this numerical model. The comparison is shown in Figure 2, from which it can be seen that the simulation results agree well with the test ones. The rockfall impact force almost increases linearly with increasing the rockfall falling height. However, the impact force calculated with these methods in Tab.2 is different, whose maximum difference is about 5 times. This is because the considered factors in each method are different. For example, the Japan Road Association’s method [34] and Switzerland’s method [35] do not consider the effect of the buffer layer thickness, resulting in larger calculation values than other methods. In general, the impact force calculation results obtained with the proposed numerical model are between those of the Switzerland’s method [35] and the tunnel manual’s method [36], which indicates that the calculation results are reasonable.

Table 2. The often-used calculation methods of the rockfall impact force.

Methods	Calculation formulae
The tunnel manual’s method [36]	$P={\displaystyle \frac{Qv}{2gh}}\sqrt{{\displaystyle \frac{E\left(1-\mu \right)}{\rho \left(1+\mu \right)\left(1-2\mu \right)}}}$
the method in “Specifications for design of highway subgrades” [38]	$P=v{\gamma }^{{\displaystyle \frac{1}{2}}}{Q}^{{\displaystyle \frac{5}{6}}}{\left({\displaystyle \frac{9\pi }{16{\gamma }_r}}\right)}^{{\displaystyle \frac{1}{6}}}\sqrt{{\displaystyle \frac{4{\tan }^4\left({45}^0+\frac{\phi }{2}\right)-2}{g}}}$
Yang and Guan’ method [37]	$P={\displaystyle \frac{100\xi m\sqrt{2gH}}{0.097mg+2.21h+0.045/H+1.2}}$
Switzerland’s method [35]	$P=1.765M_E^{{\displaystyle \frac{2}{5}}}{R}^{{\displaystyle \frac{1}{2}}}{\left(QH\right)}^{{\displaystyle \frac{3}{5}}}$
Japan Road Association’s method [34]	$P=2.108{\left(mg\right)}^{{\displaystyle \frac{2}{3}}}{\lambda }^{{\displaystyle \frac{2}{5}}}{H}^{{\displaystyle \frac{3}{5}}}$

Table 2. The often-used calculation methods of the rockfall impact force.

Methods	Calculation formulae
The tunnel manual’s method [36]	$P={\displaystyle \frac{Qv}{2gh}}\sqrt{{\displaystyle \frac{E\left(1-\mu \right)}{\rho \left(1+\mu \right)\left(1-2\mu \right)}}}$
the method in “Specifications for design of highway subgrades” [38]	$P=v{\gamma }^{{\displaystyle \frac{1}{2}}}{Q}^{{\displaystyle \frac{5}{6}}}{\left({\displaystyle \frac{9\pi }{16{\gamma }_r}}\right)}^{{\displaystyle \frac{1}{6}}}\sqrt{{\displaystyle \frac{4{\tan }^4\left({45}^0+\frac{\phi }{2}\right)-2}{g}}}$
Yang and Guan’ method [37]	$P={\displaystyle \frac{100\xi m\sqrt{2gH}}{0.097mg+2.21h+0.045/H+1.2}}$
Switzerland’s method [35]	$P=1.765M_E^{{\displaystyle \frac{2}{5}}}{R}^{{\displaystyle \frac{1}{2}}}{\left(QH\right)}^{{\displaystyle \frac{3}{5}}}$
Japan Road Association’s method [34]	$P=2.108{\left(mg\right)}^{{\displaystyle \frac{2}{3}}}{\lambda }^{{\displaystyle \frac{2}{5}}}{H}^{{\displaystyle \frac{3}{5}}}$

Besides, the failure process of the reinforced concrete plate is also adopted to valid the numerical model. Assume the rockfall vertically impacts the reinforced concrete plate for several consecutive times, the impact interval is 0.05s, the impact velocity is 20m/s, and the damage of the reinforced concrete plate is observed. As shown in Figure 3, the numerical simulation results show that the main cracks in the span firstly appear at the bottom, and then continue to propagate and coalesce, which is in good agreement with the test results. It indicates that the numerical model is suitable to simulate the failure process of the reinforced concrete plate under the rockfall Impact.

3. Dynamic response simulation of the shed-tunnel structure under the rockfall repeated impacts

The numerical simulation of the dynamic response of the shed-tunnel structure under the rockfall repeated impacts adopts the full restart technology of LS-DYNA code. The calculation result of the previous rockfall impact is taken as the calculation initial condition of the next rockfall impact, in which the initialization process is achieved with the keyword STRESS_INTIALIZATION to ensure the information continuity of the numerical calculation. Then the dynamic response laws of the shed-tunnel under the rockfall repeated impacts, such as the impact force, impact depth, vertical displacement of the shed-tunnel roof and plastic strain, are studied under different conditions such as the rockfall mass, shape, impact velocity and impact angle. The numerical calculation schemes are shown in Tab.3. Among them, the impact force refers to the time-history relationship curve of the interaction between the rockfall and the buffer layer, from which the peak impact force on the buffer layer and its change law can be judged. Impact depth refers to the vertical displacement of the surface nodes of the buffer layer after the rockfall impact, which directly indicates the impact effect, and the point directly contacted with the rockfall, namely, the center point of the upper surface of the buffer layer, is set as the displacement monitoring point. The vertical displacement of the shed-tunnel roof refers to the maximum vertical displacement of the nodes on the upper surface of the reinforced concrete plate, and the center of the upper surface is selected as the monitoring point. The shed-tunnel plastic strain directly indicates the failure of the reinforced concrete plate.

Table 3. The numerical calculation schemes.

variables	Calculation scheme	The adopted calculation model
m/kg	70.3/137.31/237.26	sphere, v=20m/s, n=5, △t =0.05s
v/m/s	15/20/25	sphere, R=0.25m, n=5, △t =0.05s
rockfall shape	sphere (R=0.25m)/square 1 (0.403m×0.403m×0.403m)/ square 2 (0.570m×0.570m×0.201)	v=20m/s, n=5, △t =0.05s
θ/°	30°/45°/60°/75°/90°	sphere, R=0.25m, v=20m/s, n=5, △t =0.05s

m, v, R are the rockfall mass, impact velocity and radius respectively, θ is rockfall impact angle, which is the angle between the rockfall impact direction and the horizontal plane; n is the number of the rockfall repeated impacts and △t is time interval of the two neighboring impacts.

Table 3. The numerical calculation schemes.

variables	Calculation scheme	The adopted calculation model
m/kg	70.3/137.31/237.26	sphere, v=20m/s, n=5, △t =0.05s
v/m/s	15/20/25	sphere, R=0.25m, n=5, △t =0.05s
rockfall shape	sphere (R=0.25m)/square 1 (0.403m×0.403m×0.403m)/ square 2 (0.570m×0.570m×0.201)	v=20m/s, n=5, △t =0.05s
θ/°	30°/45°/60°/75°/90°	sphere, R=0.25m, v=20m/s, n=5, △t =0.05s

m, v, R are the rockfall mass, impact velocity and radius respectively, θ is rockfall impact angle, which is the angle between the rockfall impact direction and the horizontal plane; n is the number of the rockfall repeated impacts and △t is time interval of the two neighboring impacts.

3.1. Effect of the rockfall mass

In order to study the effect of the rockfall mass on the dynamic response of the shed-tunnel, the rockfall impact velocity and direction are fixed to be 20m/s and perpendicular to the buffer layer respectively. The rockfall radius is assumed to 0.2m (5 impacts), 0.25m (5 impacts), and 0.3m (4 impacts, because the model is completely destroyed at the fifth impact, and the calculation cannot converge), and accordingly, the rockfall mass is 70.3kg, 137.31kg and 237.26kg, respectively. The interval between two neighboring impacts is 0.05s.

3.1.1. The impact force on the buffer layer

The following results can be obtained from the calculation results shown in Figure 4 and Table 4. First of all, after each impact, the pattern of the time history curve of the impact force is basically the same. The impact force firstly increases to the maximum and then decreases to 0 with time going on. Second, when the rockfall mass is fixed, with increasing the impact number, the peak impact force increases rapidly after the second impact, and then tends to be relatively stable. When the rockfall is 70.3kg, the peak impact force increases from 221.32kN of the first impact to 259.12kN of the second impact, whose increase extent is 17.08%. And then the increase extent of subsequent impacts is not as high as that of the second impact. The peak impact force reaches the maximum 298.77kN in the fifth impact, which increased by 15.3% compared with the second impact. When the rockfall mass is 137.31kg, the peak impact force rapidly increases from 311.36kN of the first impact to 445.18kN of the second impact, whose increase extent is 42.98%. The peak impact force reaches the maximum 484.25kN at the fifth impact, which increased by 8.78% compared with the second impact. When the rockfall mass is 237.26kg, the peak impact force rapidly increases from 546.89kN of the first impact to 622.54kN of the second impact, whose increase extent is 13.83%. The peak impact force reaches the maximum 677.58kN at the third impact, which increased by 8.84% compared with the second impact. This is mainly because the buffer layer is not compressed at the first impact, and its buffer performance is the best. With increasing the impact number, its buffer capacity gradually weakens, so the peak impact force will increase at the first several times. However, with the impact continuing, the buffer layer can hardly be compressed anymore, so its buffer performance tends to be stable, resulting in the peak impact force also tends to be stable. At the same time, it can also be seen that the larger the rockfall mass is, the less the impact number required to reach the maximum is. Finally, the lager the rockfall mass is, the greater the peak impact force is.

3.1.2. The impact depth in the buffer layer

The calculation results of the time history curve of the impact depth in the buffer layer with different rockfall mass is shown in Figure 5, from which the following conclusions can be obtained. First of all, after each impact, the buffer layer will have a certain degree of rebound, because it is an elastic-plastic material, and the elastic deformation will recover after the impact force is unloaded. Second, when the rockfall mass is fixed, with increasing the impact number, the maximum impact depth gradually increases, but its increase extent is less than that of the first impact. When the rockfall mass is 70.3kg, the impact depth of the first impact is 0.099m, and that of the second impact is 0.035m, which is only 35.35% of the first one. When the rockfall mass is 137.31kg, the impact depth of the first impact is 0.115m, and that of the second impact is 0.039m, which is only 33.91% of the first one. When the rockfall mass is 237.26kg, the impact depth of the first impact is 0.119m, and that of the second impact is 0.185m, which is only 55.46% of the first one. This is mainly because the soil is compacted after the first impact, resulting in the impact depth of the first impact is much larger than that of the subsequent impacts. Meanwhile when the rockfall mass is small (70.3kg), the rebound ratio of each impact is large. However, when the rockfall mass is large (237.26kg), the rebound ratio is little. Finally, the larger the rockfall mass is, the greater the impact depth of the buffer layer is.

3.1.3. The maximum plastic strain of the rebar

The calculation results of the time history curve of the maximum plastic strain of the rebar with different rockfall mass is shown in Figure 6, from which the following conclusions can be obtained. First of all, the maximum plastic strain of the rebar shows a monotonically increasing trend on the whole. Second, the plastic strain of the rebar all increases rapidly first, and then tends to be stable. This is because the deformation in the reinforced concrete plate will drive the rebar to produce plastic strain. When the reinforced concrete plate reaches balance, the rebar no longer generates more plastic strain and achieves stability. Third, the larger the rockfall mass is, the faster the maximum plastic strain of the rebar increases. The maximum plastic strain of the rebar is sensitive to the rockfall mass, and the plastic strain produced by five rockfall impacts of 70.3kg is less than that produced by one rockfall impact of 237.26kg. When the rockfall mass is small (70.3kg), the plastic strain is not generated in the first few impacts, but when it is large (237.26kg), the plastic strain is generated in the first impact. When the rockfall mass is 237.26kg, the plastic strain produced by four impacts is 0.0673. When the rockfall mass is 70.3kg, the plastic strain produced by five impacts is 0.00467, which is only 6.94% of the rockfall mass of 237.26kg. When the rockfall mass is 137.31kg, the plastic strain produced by five impacts is 0.0464, which is 68.95% of the rockfall mass of 237.26kg.

3.1.4. The vertical displacement of the shed roof

The calculation results of the vertical displacement of the shed roof with different rockfall mass is shown in Figure 7, from which the following conclusions can be obtained. First of all, there is a rebound phenomenon in the vertical displacement of the shed roof, which is mainly speculated to be caused by the existence of the rebar. Because concrete is a typical brittle material, no rebound will occur in it in theory, so it is speculated that the existence of the rebar leads to the rebound of the shed roof. To verify this hypothesis, a comparative experiment is designed. The rockfall radius is set to 0.25m, and accordingly its mass is 137.31kg. The impact velocity is set to 20m/s, and the impact direction is vertically downward. The time history curve of vertical displacement of the shed roof is compared in one model with rebar and the other model without rebar. It can be seen from the calculation result shown in Figure 8 that when it has the rebar, the shed roof has obvious rebound, while when it has no the rebar, the shed roof vertical displacement increases monotonically, and there is no rebound phenomenon. It shows that the rebar has a great influence on the strength and deformation of the reinforced concrete plate. If the shed roof has no rebar, it is very easy to fail. Second, when the rockfall mass is large (137.31kg, 237.26kg), the vertical displacement of the shed roof gradually increases with increasing the impact number. When the rockfall mass is small (70.3kg), the rebound ratio of the vertical displacement of the shed roof after each impact is large, and the vertical displacement of the shed roof after stability has no much change compared with the initial state. However, when the rockfall mass is large (237.26kg), the rebound ratio is small. Finally, with increasing the rockfall mass, the vertical displacement of the shed roof increases significantly.

3.1.5. The maximum axial force of the rebar

The calculation results of the maximum axial force of the rebar with different rockfall mass is shown in Figure 9, from which the following conclusions can be obtained. First of all, when the rockfall mass is fixed, in the first few impacts, the axial force of the rebar basically shows a trend of increasing from zero to the maximum and then decreasing to zero. However, after several subsequent impacts, the axial force of the rebar will not recover to zero, because the plastic deformation of the rebar cannot recover, so that a certain axial force remains. Compared with the upper and lower rebars, the axial force of the lower rebars is obviously larger than that of the upper rebars. Because the lower rebars bear the larger tensile stress, it is important to strengthen the reinforcement of the lower part of the reinforced concrete plate. What is more, comparing the calculation results with three kinds of rockfall mass, it can be seen that the larger the rockfall mass is, the larger the variation of the rebar axial force is.

3.1.6. The plastic strain of the shed-tunnel

The calculation results of the plastic strain contour of the shed-tunnel with different rockfall mass is shown in Figure 10, from which the following conclusions can be obtained. After many impacts, large plastic deformation occurs in the shed-tunnel, which mainly concentrates in and around the center of the shed-tunnel. The larger the rockfall mass is, the larger the plastic zone is and the greater the damage caused by the rockfall to the shed-tunnel is. The plastic strain of the shed-tunnel is sensitive to the rockfall mass, and the damage range of the shed-tunnel produced by five impacts of the rockfall with 70.3kg is far less than that produced by four impacts of rockfall 237.26kg. The damage of the shed-tunnel begins from the bottom of the span and cracks gradually propagate towards the upper and both sides of the reinforced concrete plate.

3.2. Effect of the rockfall impact velocity

3.2.1. The impact force on the buffer layer

The calculation results of the time history curve of the impact force on the buffer layer with different rockfall impact velocity and its peak are shown in Figure 11 and Table 5 respectively. The following conclusions can be obtained. First of all, when the impact velocity is fixed, the peak impact force increases rapidly after the second impact and then tends to be stable with increasing the impact number. When the rockfall impact velocity is 15m/s， the peak impact force increases from 266.52kN of the first impact to 326.66kN of the second impact, whose increase extent is 22.56%. And then the increase extent of subsequent impacts is not as high as that of the second impact. The peak impact force reaches the maximum 363.20kN in the fifth impact, which increases by 11.19% compared with the second impact. When the rockfall impact velocity is 20m/s, the peak impact force rapidly increases from 311.36kN of the first impact to 445.18kN of the second impact, whose increase extent is 42.98%. The peak impact force reaches the maximum 484.25kN at the fifth impact, which increased by 8.78% compared with the second impact. When the rockfall impact velocity is 25m/s, the peak impact force rapidly increases from 421.33kN of the first impact to 516.98kN of the second impact, whose increase extent is 22.70%. The peak impact force reaches the maximum 569.11kN at the fourth impact, which increases by 10.08% compared with the second impact. It indicates that the larger the rockfall impact velocity is, the less the impact number required to reach the peak impact force is. What is more, comparing the time history curves of the impact force with these three velocities, it can be seen that the lager the rockfall mass is, the greater the peak impact force is.

3.2.2. The impact depth in the buffer layer

The calculation results of the time history curve of the impact depth in the buffer layer with different rockfall impact velocity is shown in Figure 12, from which the following conclusions can be obtained. First of all, when the rockfall impact velocity is fixed, with increasing the impact number, the maximum impact depth gradually increases, but its increase extent is less than that of the first impact. When the rockfall impact velocity is 15m/s, the impact depth of the first impact is 0.095m, and that of the second impact is 0.027m, which is only 28.42% of the first one. When the rockfall impact velocity is 20m/s, the impact depth of the first impact is 0.115m, and that of the second impact is 0.039m, which is only 33.91% of the first one. When the rockfall impact velocity is 25m/s, the impact depth of the first impact is 0.140m, and that of the second impact is 0.056m, which is only 40% of the first one. However, it can be seen that the increase extent increases much more after the fourth impacts. This is mainly because the shed-tunnel is seriously destroyed, which leads to much displacement of the buffer layer. Meanwhile when the rockfall impact velocity is small (15m/s), the rebound ratio of each impact is large. However, when the rockfall impact velocity is large (25m/s), the rebound ratio is little. Finally, the greater the rockfall impact velocity is, the greater the impact depth of the buffer layer is.

3.2.3. The maximum plastic strain of the rebar

The calculation results of the time history curve of the maximum plastic strain of the rebar with different rockfall impact velocity is shown in Figure 13, from which the following conclusions can be obtained. First of all, the maximum plastic strain of the rebar shows a monotonically increasing trend on the whole. Second, the plastic strain of the rebar all increases rapidly first, and then tends to be stable. Second, with increasing the rockfall impact velocity, the maximum plastic strain of the rebar increases, which is sensitive to the rockfall impact velocity, and the plastic strain produced by five rockfall impacts of 15m/s is less than that produced by one rockfall impact of 25m/s. When the rockfall impact velocity is small (15m/s), the plastic strain is not generated in the first few impacts, but when the rockfall impact velocity is large (25m/s), the plastic strain is generated at the second impact. When the rockfall impact velocity is 25m/s, the plastic strain produced by five impacts is 0.058. When the rockfall impact velocity is 15m/s, the plastic strain produced by five impacts is 0.00556, which is only 9.59% of the rockfall impact velocity is 25m/s. When the rockfall impact velocity is 20m/s, the plastic strain produced by five impacts is 0.0464, which is 79.93% of the rockfall impact velocity is 25m/s.

3.2.4. The vertical displacement of the shed roof

The calculation results of the vertical displacement of the shed roof with different rockfall impact velocity is shown in Figure 14, from which the following conclusions can be obtained. First of all, when the rockfall impact velocity is fixed, the vertical displacement of the shed roof gradually increases with increasing the impact number, and it has the rebound phenomenon. When the rockfall impact velocity is little, the displacement rebound is rather larger. While when the rockfall impact velocity is large, the displacement rebound is rather less. Meanwhile, with increasing the rockfall impact velocity, the vertical displacement of the shed roof increases significantly.

3.2.5. The maximum axial force of the rebar

The calculation results of the maximum axial force of the rebar with different rockfall impact velocity is shown in Figure 15, from which the following conclusions can be obtained. First of all, when the rockfall impact velocity is fixed, in the first few impacts, the axial force of the rebar basically shows a trend of increasing from zero to the maximum and then decreasing to zero. However, after several subsequent impacts, the axial force of the rebar will not recover to zero, because the plastic deformation of the rebar produces, so that a certain axial force remains. Comparing with the upper and lower rebars, the axial force of the lower rebars is obviously larger than that of the upper rebars, which indicates that the lower rebars bear the larger tensile stress. What is more, comparing the calculation results with three rockfall impact velocities, it can be seen that the larger the rockfall impact velocity is, the larger the variation of the rebar axial force is.

3.2.6. The plastic strain of the shed-tunnel

The calculation results of the plastic strain contour of the shed-tunnel with different rockfall impact velocity is shown in Figure 16, from which the following conclusions can be obtained. After many impacts, large plastic deformation occurs in the shed-tunnel, which mainly concentrates in and around the center of the shed-tunnel. The larger the rockfall impact velocity is, the larger the plastic zone is and the greater the damage caused by the rockfall to the shed-tunnel is. The plastic strain of the shed-tunnel is sensitive to the rockfall impact velocity, and the damage range of the shed-tunnel produced by five impacts of the rockfall with 15m/s is far less than that produced by five impacts of rockfall 25m/s. When the rockfall impact velocity is 25m/s, the shed-tunnel almost completely fails after five impacts.

3.3. Effect of the rockfall shape

Assume the rockfall mass is all 137.31kg, accordingly three kinds of rockfall shapes are sphere with radius 0.25m, cuboid 1 with 0.403m×0.403m×0.403m and cuboid 2 with 0.570m×0.570m×0.201m respectively. The contact faces of the latter two rockfalls are 0.403m×0.403m and 0.570m×0.570m respectively. The other calculation parameters are shown in Tab.3.

3.3.1. The impact force on the buffer layer

The following results can be obtained from the calculation results shown in Figure 17 and Table 6. First of all, for the spherical rockfall, with increasing the impact number, the peak impact force increases rapidly from 311.36kN of the first impact to 445.18kN of the second impact, with an increase of 42.98%. And the growth rate of the subsequent impacts is all lower than that of the second one. The peak impact force reaches 484.25kN at the fifth impact, which increases by 8.78% compared with the second impact. For the cuboid rockfall, the peak impact force of the first impact is the largest, but that of the subsequent impacts decreases rapidly. When the bottom surface is 0.403m×0.403m, the impact force generated by the first impact is the largest, reaching 1180.86kN, and then it decreases with increasing the impact number. It reaches the minimum of 518.67kN at the fifth impact, which is only 43.92% of the first one. When the bottom surface is 0.570m×0.570m, the peak impact force produced by the first impact is 1960.09kN, which is much larger than the subsequent ones. The peak impact force produced after the second impact reaches the minimum of 1090.65kN, which is only 55.64% of the first one. Meanwhile, it can also be found that the impact force generated by the cuboid rockfall is much larger than that of the spherical rockfall, and the larger the contact area of the cuboid is, the greater the impact force is.

3.3.2. The impact depth in the buffer layer

The calculation results of the time history curve of the impact depth in the buffer layer with different rockfall shape is shown in Figure 18, from which the following conclusions can be obtained. First of all, with increasing the impact number, the maximum impact depth gradually increases, but its increase extent is less than that of the first impact. When the rockfall shape is sphere, the impact depth of the first impact is 0.0115m, and that of the second impact is 0.039m, which is only 33.91% of the first one. For the cuboid rockfall, when its bottom surface is 0.403m×0.403m, the impact depth generated by the first impact is 0.081m, and that of the second impact is 0.033m, which is only 40.74% of the first one. when its bottom surface is 0.570m×0.570m, the impact depth generated by the first impact is 0.069m, and that of the second impact is 0.026m, which is only 37.68% of the first one. Meanwhile, the rebound displacement caused by the cuboid rockfall is much larger than that of the spherical one. Finally, the impact depth caused by the cuboid rockfall is less than that of the spherical one. For the cuboid rockfall, the larger the contact area is, the less the impact depth is.

3.3.3. The maximum plastic strain of the rebar

The calculation results of the time history curve of the maximum plastic strain of the rebar with different rockfall shape is shown in Figure 19, from which the following conclusions can be obtained. First of all, the maximum plastic strain of the rebar shows a monotonically increasing trend on the whole. The plastic strain of the rebar all increases rapidly first, and then tends to be stable. Second, when the rockfall mass and impact velocity is fixed, the effect of the rockfall shape has much effect on the maximum plastic strain of the rebar. Especially for cuboid rockfall, the larger the contact face is, the less the maximum plastic strain of the rebar is. When the rockfall is sphere, the plastic strain produced by five impacts is 0.0464. When the rockfall is cuboid with the bottom surface is 0.403m×0.403m, the plastic strain produced by five impacts is 0.0481, which is 103.66% of the spherical rockfall. When the rockfall is cuboid with the bottom surface is 0.570m×0.570m, the plastic strain produced by five impacts is 0.0294, which is 63.36% of the spherical rockfall.

3.3.4. The vertical displacement of the shed roof

The calculation results of the vertical displacement of the shed roof with different rockfall shape is shown in Figure 20, from which the following conclusions can be obtained. First of all, it can be seen that with increasing the impact number, the vertical displacement of the shed roof gradually increases. When the rockfall is sphere, the rebound ratio of each impact is large. However, when the rockfall is cuboid, the proportion of rebound is relatively small. Moreover, it can also be seen that when the rockfall is cuboid, the vertical displacement of the shed roof increases faster than that of spherical rockfall. Meanwhile the larger the contact area is, the slower the vertical displacement of the shed roof increases.

3.3.5. The maximum axial force of the rebar

The calculation results of the maximum axial force of the rebar with different rockfall shape is shown in Figure 21, from which the following conclusions can be obtained. First of all, when the rockfall shape is fixed, the upper and lower layers of the rebars fluctuate in two states of tension and compression. The lower rebar is mainly subjected to tension, while the upper rebar is mainly subjected to compression, which is because the lower rebar bears the tensile stress of the lower part of the shed-tunnel structure. Meanwhile, it can be seen that the rockfall shape has small effect on the axial force of the rebar.

3.3.6. The plastic strain of the shed-tunnel

The calculation results of the plastic strain contour of the shed-tunnel with different rockfall shape is shown in Figure 22, from which the following conclusions can be obtained. After many impacts, large plastic deformation occurs in the shed-tunnel, which mainly concentrates in and around the center of the shed-tunnel. The plastic strain of the shed-tunnel is sensitive to the rockfall shape, and the plastic deformation of the shed-tunnel produced by the cuboid rockfall is far less than that produced by the spherical rockfall. For the cuboid rockfall, the less that contact area is, the larger the plastic deformation and damage of the shed-tunnel is. When the impact number reaches five, the shed-tunnel almost completely fails.

3.4. Effect of the rockfall impact angle

3.4.1. The impact force on the buffer layer

The following results can be obtained from the calculation results shown in Figure 23 and Table 7. First of all, with increasing the impact number, the peak impact force increases rapidly from the second impact, and then tends to be stable. When the impact angle is 30°, the peak impact force increases from 153.25kN of the first impact to 339.96kN of the second one, an increase of 53.40%. It reaches 292.45kN at the fifth impact, which increases by 24.40% compared with the second one. When the impact angle is 45°, the peak impact force increases from 215.55kN of the first impact to 333.96kN of the second one, an increase of 54.93%. It reaches 354.67kN in the fifth impact, which increases by 6.20% compared with the second one. When the impact angle is 60°, the peak impact force increases by 45.47% from 268.76kN of the first impact to 390.97kN in the second impact. The peak impact force reaches 432.65kN at the fourth impact, which increases by 10.66% compared with the second impact. When the impact angle is 75°, the peak impact force increases from 301.69kN of the first impact to 409.63kN of the second impact, an increase of 35.78%. The peak impact force reaches 457.75kN at the fifth impact, an increase of 11.75% compared with the second one. When the impact angle is 90°, the peak impact force increases from 311.36kN of the first impact to 445.18kN of the second one, an increase of 42.98%. The peak impact force reaches 484.25kN at the fifth impact, which increases by 8.78% compared with the second one. Moreover, it can be seen that the greater the impact angle is, the larger the peak impact force is.

3.4.2. The impact depth in the buffer layer

The calculation results of the time history curve of the impact depth in the buffer layer with different rockfall impact angle is shown in Figure 24, from which the following conclusions can be obtained. First of all, with increasing the impact number, the maximum impact depth gradually increases, but its increase extent is less than that of the first impact. When the rockfall impact angle is 30°, the impact depth of the first impact is 0.056m, and that of the second impact is 0.027m, which is only 48.21% of the first one. When the rockfall impact angle is 45°, the impact depth of the first impact is 0.082m, and that of the second impact is 0.027m, which is only 32.93% of the first one. When the rockfall impact angle is 60°, the impact depth of the first impact is 0.084m, and that of the second impact is 0.044m, which is only 52.38% of the first one. When the rockfall impact angle is 75°, the impact depth of the first impact is 0.086m, and that of the second impact is 0.045m, which is only 52.33% of the first one. When the rockfall impact angle is 90°, the impact depth of the first impact is 0.115m, and that of the second impact is 0.039m, which is only 33.91% of the first one. Moreover, it can also be found that the larger the impact angle is, the greater the impact depth of the buffer layer is.

3.4.3. The maximum plastic strain of the rebar

The calculation results of the time history curve of the maximum plastic strain of the rebar with different rockfall impact angle is shown in Figure 25, from which the following conclusions can be obtained. First of all, the maximum plastic strain of the rebar shows a monotonically increasing trend on the whole. The plastic strain of the rebar all increases rapidly first, and then tends to be stable. Second, it can also be found that the larger the rockfall impact angle is, the faster the maximum plastic strain of the rebar increases. The maximum plastic strain of the rebar is sensitive to the rockfall impact angle. When the impact angle is small (30°), the plastic strain is nearly not generated, but when it is large (90°), the plastic strain is generated in the first impact. When the rockfall impact angle is 90°, the plastic strain produced by five impacts is 0.06187. When the rockfall impact angle is 30°, the plastic strain produced by five impacts is 0.0001931, which is only 0.31% of the rockfall impact angle of 90°. When the rockfall impact angle is 45°, the plastic strain produced by five impacts is 0.01445, which is 23.36% of the rockfall impact angle of 90°. When the rockfall impact angle is 75°, the plastic strain produced by five impacts is 0.04636, which is 74.93% of the rockfall impact angle of 90°.

3.4.4. The vertical displacement of the shed roof

The calculation results of the vertical displacement of the shed roof with different rockfall impact angle is shown in Figure 26, from which the following conclusions can be obtained. First of all, it can be seen that with increasing the impact number, the vertical displacement of the shed roof gradually increases. When the rockfall impact angle is small, the rebound ratio of each impact is large. However, when the rockfall impact angle is large, the proportion of rebound is relatively small. Moreover, it can also be seen that the larger the impact angle is, the faster the vertical displacement of the shed roof increases.

3.4.5. The maximum axial force of the rebar

The calculation results of the maximum axial force of the rebar with different rockfall impact angle is shown in Figure 27, from which the following conclusions can be obtained. First of all, when the rockfall impact angle is fixed, the upper and lower layers of the rebars fluctuate in two states of tension and compression. The lower rebar is mainly subjected to tension, while the upper rebar is mainly subjected to compression, which is because the lower rebar bears the tensile stress of the lower part of the shed-tunnel structure. Meanwhile, it can be seen that the impact angle has small effect on the axial force of the rebar.

3.4.6. The plastic strain of the shed-tunnel

The calculation results of the plastic strain contour of the shed-tunnel with different rockfall impact angle is shown in Figure 28, from which the following conclusions can be obtained. After many impacts, large plastic deformation occurs in the shed-tunnel, which mainly concentrates in and around the center of the shed-tunnel. The larger the rockfall impact angle is, the larger the plastic zone is and the greater the damage caused by the rockfall to the shed-tunnel is. The maximum plastic strain of the shed-tunnel is sensitive to the rockfall impact angle, and the plastic strain of the shed-tunnel produced by five impacts of the rockfall impact angle 30° is far less than that produced by five impacts of rockfall impact angle 90°.

4. Conclusions

In view of the restart function of ANSYS/LS-DYNA code, the dynamic response of the reinforced concrete shed-tunnel structure under rockfall repeated impacts is studied. Firstly, the FEM-SPH coupling numerical calculation model for the reinforced concrete shed-tunnel structure under rockfall repeated impacts is established, and then its validity is verified with the model test results. Then, the parametric sensitivity analysis is adopted to study the effect of the following four factors, e.g. rock mass, shape, impact angle and impact velocity, on the dynamic response of the shed-tunnel structure under rockfall repeated impacts.

In all, this research provides a new way to study the dynamic response of the reinforced concrete shed-tunnel structure under rockfall repeated impacts, and the research results is useful for the design of the shed-tunnel structure.

However, it is found that some simplified hypotheses are made in establishing the numerical model, which will lead to the numerical simulation results are not in good agreement with the practical ones. Therefore, the study should be made further in the future.