# A Blast-Resistant Method Based on Wave Converters with Spring Oscillator for Underground Structures

Authors: Yu Zhang; Yuanxue Liu; Runze Wu; Jichang Zhao; Ming Hu; Yizhong Tan
Journal: Shock and Vibration
 (2017)

---

## Abstract

Researches on blast-resistant measures for underground structures such as tunnels and underground shopping malls are of great importance for their significant role in economic and social development. In this paper, a new blast-resistant method based on wave converters with spring oscillator for underground structures was put forward, so as to convert the shock wave with high frequency and high peak pressure to the periodic stress wave with low frequency and low peak pressure. The conception and calculation process of this new method were introduced. The mechanical characteristics and motion evolution law of wave converters were deduced theoretically. Based on the theoretical deduction results and finite difference softwareFLAC3D, the dynamic responses of the new blast-resistant structure and the traditional one were both calculated. Results showed that, after the deployment of wave converters, the peak absolute values of the bending moment, shear force, and axial force of the structure decreased generally, which verified the good blast-resistant effect of the new blast-resistant method.


---

## Body

## 1. Introduction

Tunnels, subway stations, underground shopping malls, and so forth play an important role in the economic and social development, as well as the personal and property safety. In recent years, there have been a variety of researches on the blast-resistant measures for underground structures, mainly focused on the design and optimization on structures and materials.Blast-resistant measures based on the structure optimization mainly contain the increase of stiffness or adoption of structures good for the reflection, diffraction, and scattering of the stress wave. Usually the arch structures and structures containing holes have better performance of wave dissipation than rectangular structures or structures without holes, so these kinds of structures have attracted many attentions of scholars [1–3]. A good blast-resistant performance is also embodied in the plate-foam composite structure [4], box-shaped steel plate-reinforced concrete composite structure [5], carbon fiber reinforced composite structure [6, 7], and prestressed structure [8]. Via the explosion model tests, Yang et al. studied the antiexplosion effect of prestressed cable reinforced caverns. The results showed that prestressed structures are good for antiexplosion [9]. Kobielak et al. experimentally investigated the influence of an attenuation barrier on soil stresses and pressures acting on a buried silo, caused by underground explosions at different distances. Test results verified the blast-resistant effectiveness caused by the barrier composed of the cylindrical tubes [10]. The dynamic response of a buried silo caused by underground explosion was also studied experimentally via the measured pressures [11].Blast-resistant measures based on the material optimization are mainly concentrated on the development of porous or lightweight materials of low stiffness, and the materials are developing gradually from the traditional inorganic porous materials or lightweight materials to the polymer materials and porous metal materials currently. Materials such as the rigid polyurethane foam, polypropylene fiber concrete, rubber concrete, foam concrete, foamed aluminium, and steel fiber reinforced concrete are good choices for the blast-resistant materials. Yakushin et al. investigated the properties of low-density rigid polyurethane foams with hollow glass microspheres. The tension and compression properties in relation to the content of microspheres were determined in their work [12]. Alhozaimy et al. found that the flexural toughness and impact resistance showed an increase in the presence of polypropylene fibers [13]. The mechanical properties of concrete containing tire-rubber particles were studied by Khaloo et al. [14]. They found that, unlike plain concrete, the failure state in rubberized concrete occurred gently and uniformly and did not cause any separation in the specimen. In fact, this property is good for dynamic protection. Hernández-Olivares and Barluenga also studied the fire performance of the rubber-filled concrete [15]. Lim et al. studied the compressive, splitting tensile, and flexural strengths of lightweight foamed concrete [16]. The results have shown that the foamed concrete is a good choice for dynamic protection.In the blast-resistant methods of the traditional structure optimization, the construction process is usually complex, and sometimes the function of structures may even be affected. The shock wave mainly consists of high frequency components. In the blast-resistant methods of porous or lightweight materials, the materials are easy to get damaged unrecoverably and have large deformation under the blasting load because of the low elastic modulus and the existence of the holes. Thereby the overall stability of underground structure and surrounding rock may be affected by the large deformation.In order to improve the traditional antiknock methods, a new blast-resistant method based on wave converters with spring oscillator for underground structures is put forward in this paper. The new method mainly consists of an array of wave converters and a distribution layer. Firstly, the conception of the new method is introduced, including the formation of the wave converter and distribution layer. Secondly, the calculation process of dynamic responses for underground structures adopting the new blast-resistant method is presented. Thirdly, the mechanical characteristics and motion evolution law of the wave converter are derived, including the static constitutive relation of the wave converter, dynamic response partitioning of the wave converter, differential equation of motion for the spring oscillator, displacement transfer coefficient of the wave converter, and the stress inversion of the wave converter. A case study is also conducted to verify the applicability and rationality of the new method by comparing with the traditional structure.


## 2. Conception of New Blast-Resistant Method

The new method mainly consists of an array of wave converters and a distribution layer, shown in Figure1. The wave converter includes 2 shells containing a length adjusting rod and a spring oscillator comprised of springs and a mass block, shown in Figure 2. The initial length of the wave converter can be adjusted via the adjustment of the length adjusting rod, which can control the prestress of the compression springs. Via the adjustment of the initial converter length, the wave converter can also be easily installed under different reserve space between the rock and tunnel roof. The wave converter is compressible when the load on the top of the converter exceeds the spring prestress. The distribution layer is composed of the material with relatively low wave impedance, shown in Figure 1. Multiple reflections caused by the periodic stress wave can occur in the distribution layer, which results in the further energy dissipation.Figure 1
Schematic of the new blast-resistant method.Figure 2
Schematic of the wave converter.The new blast-resistant method combines such mechanisms as the spring deformation, inertia, and periodic vibration of the mass block to provide the resistance against the dynamic load. The self-support capacity of the surrounding rock can also be fully utilized. Via the wave converter, the shock wave with high frequency and high peak pressure can be transformed to the periodic stress wave with low frequency and low peak pressure. Thereby the shock wave is dispersed and materials under the converter can be prevented from crushing. Under the blasting load, the deformation process of the wave converter can be divided into such 3 periods as rapid loading stage, rapid unloading stage, and slow unloading stage.The above 3 stages are determined by the relative displacementΔui(t) between the top and bottom of wave converter. Under the impact loads propagating in the rock, the typical curve of Δui(t) with time for the wave converter is shown in Figure 3. As is shown in Figure 3, the rapid loading stage corresponds to the sharp increasing period, and the rapid unloading stage corresponds to the fast decreasing period, while the slow unloading stage corresponds to the slow change period.Figure 3
Typical curve ofΔui(t) with time for the wave converter.


## 3. Calculation Process of Dynamic Responses for Structures with the New Blast-Resistant Method

(1)
Establish numerical models and acquire the data needed for the calculation of the wave converter’s displacements in the rapid loading stage and rapid unloading stage. Model 1 without wave converters for the finite element analysis is set up, shown in Figure4. In model 1, the distribution layer, underground structure, and rock are established. The initial pressure caused by the wave converter on the distribution layer surface (at the wave converter bottom) is applied as the lower spring force divided by the cross-sectional area of the wave converter. Then the explosion position and blasting load should be applied. After that, under the blasting load, the finite element analysis method is used to calculate the vertical displacement-time curves u0i(t) and u2i(t) of the rock surface at the wave converter top and the distribution layer surface at the wave converter bottom, respectively. The monitoring positions for displacements are suggested to adopt the tops and bottoms of 9 wave converters along the width direction of the structure, shown in Figure 4. Then model 2 of the ground without any construction such as the structure, distribution layer, and wave converters is built up, shown in Figure 5. The numerical analysis based on model 2 under the same explosion condition as model 1 is conducted to obtain the vertical rock pressure-time curve PA0(t) for the corresponding monitoring positions of the rock surface in model 1, shown in Figure 5.
(2)
Solve the differential equation of motion for the spring oscillator and obtain the law of motion in the rapid loading stage and rapid unloading stage. In model 1, the vertical displacement-time curveu0it of the rock surface at the wave converter top multiplied by a displacement transfer coefficient Ki0<Ki<1 is considered as the vertical displacement-time curve u1it of the wave converter top in the rapid loading stage and rapid unloading stage of the true situation. The true situation refers to the real underground structure with an array of wave converters and a distribution layer. The vertical displacement-time curve u2i(t) of the distribution layer surface at the wave converter bottom in model 1 is considered directly as the vertical displacement-time curve of the wave converter bottom in the rapid loading stage and rapid unloading stage of the true situation. In the true situation, the difference between u1i(t) and u2i(t) is the relative displacement Δui(t) between the top and bottom of the wave converter. The wave converter is in the rapid loading stage before the relative displacement reaches the maximum, while it is in the rapid unloading stage during the sharp decrease period after the relative displacement reaches the maximum. Taking the displacements of the wave converter top and bottom as the boundary condition, the differential equation of motion for the spring oscillator can be solved to obtain the law of motion in the rapid loading stage and rapid unloading stage. Then the stress-time curve of the wave converter top and bottom in above 2 stages can be gotten.
(3)
Solve the differential equation of motion for the spring oscillator and obtain the law of motion in the slow unloading stage. According to step  (2), the length of the wave converter, the displacement (or position coordinate) and velocity of the oscillator at the end of the rapid unloading step can be gotten. Then the differential equation of motion for the spring oscillator in the slow unloading stage can be solved. In the slow unloading stage, the length of the wave converter can be deemed as a constant. In that case the calculation results of the dynamic response are a little dangerous which would be safer for the structure design. After that, the stress-time curve of the wave converter top and bottom in the slow unloading stage can be gotten.
(4)
Calculate the dynamic response of the underground structure. Firstly, the initial pressurePin caused by the wave converter on the distribution layer surface in model 1 should be deleted. Then the stress-time curves of the wave converter top and bottom in three stages are applied on the rock surface and distribution layer surface, instead of the wave converters. After that, the blasting load is applied and the dynamic calculation is conducted to get the dynamic response of the structure.Figure 4
Schematic of model 1.Figure 5
Schematic of model 2.


## 4. Mechanical Characteristics and Motion Evolution Law of Wave Converters

### 4.1. Static Constitutive Relation of the Wave Converter

l
0 is the initial length of the wave converter, while a is the shell thickness. So the initial clear length of the wave converter can be written as ln0=l0-2a, and ln0 is shown in Figure 6. In Figure 6, the dashed line is the position under the static equilibrium, and the solid line is the position at any time under the dynamic load. k1 is the stiffness coefficient of the upper spring, whose length in free state is l10.Δx10is the initial amount of compression for upper spring. k2is the stiffness coefficient of the lower spring, and l20 is its length in free state. Δx20is the initial amount of compression for lower spring. m and l3 are the mass and height of the mass block, respectively. PA(t) is the vertical rock pressure on the wave converter top in the true situation, while PA0(t) is the vertical rock pressure at the same location in model 2. When the wave converter is installed with the initial length l0, under static equilibrium, the geometric equation is established, yielding(1)ln0=l1+l2+l3=l10+l20+l3-Δx10-Δx20.Figure 6
Deformation process of the wave converter.As for the mass block, the balance equation is(2)mg+k1Δx10=k2Δx20.Combining (1) and (2), the solutions are(3)Δx10=k2l10+l20+l3-ln0-mgk1+k2,Δx20=mgk2+k1l10+l20+l3-ln0k1+k2-k1mgk2k1+k2.F
1 and F2 are the force of upper and lower springs, respectively. s is the cross-sectional area of the wave converter. d is the cross-sectional length. If the weight of the wave converter’s shells is ignored, F1 can be expressed as(4)F1t=k1Δx10,sPAt<k1Δx10sPAt,sPAt≥k1Δx10;F
2 can be expressed as(5)F2t=mg+k1Δx10sPAt<k1Δx10sPAt+mgsPAt≥k1Δx10.So the static constitutive relation of the wave converter yields(6)IfsPAt<k1Δx10,Δl=0(7)IfsPAt≥k1Δx10,Δl=k1+k2sPAtk1k2+mg-k2l10+l20+l3-ln0k2k1+k2.


### 4.2. Dynamic Response Partitioning of the Wave Converter

The rapid loading stage, rapid unloading stage, and slow unloading stage correspond to the rapid compression stage, rapid recovery stage, and slow recovery stage, respectively. In the slow recovery stage, the length of the wave converter can be deemed as a constant. The demarcation point of stage 1 and stage 2 is that the relative displacement reaches the maximum. The duration time of stage 2 can be determined by the relative displacement-time curve. When the rock masses above the structure are in a wide range of elasticity state, the rapid unloading stage can not be ignored, but if the rock masses are in a wide range of plastic state, the rapid unloading stage can be ignored.The computing time of stage 3 is advisable for 1 or 2 vibration periods. Via a large amount of computations, it is concluded that the computing time of stage 3 can be taken as 1 vibration period if stage 2 can not be ignored; otherwise it can be taken as 2 vibration periods.


### 4.3. Differential Equation of Motion for the Spring Oscillator

#### 4.3.1. Differential Equation of Motion in Rapid Loading Stage and Rapid Unloading Stage

N is the total number of wave converters along the width direction of the structure, and i is the serial number of the wave converter. For wave converter i, u0i(t) is the vertical displacement of the rock surface at the wave converter top in model 1, and u2iti=1,2,…,n is the vertical displacement of the wave converter bottom in rapid loading stage and rapid unloading stage of the true situation. u1i(t) is the vertical displacement of the wave converter top in rapid loading stage and rapid unloading stage of the true situation, which can be expressed as u1i(t)=Kiu0i(t). Ki is the displacement transfer coefficient of the wave converter. Δui(t) is the relative displacement between the top and bottom of the wave converter, which can be written as Δui(t)=u1i(t)-u2it. At the end of the rapid unloading stage is the initial state of the slow unloading stage, which can be deemed as the fixed-length vibration, shown in Figure 7. In Figure 7, the dashed line is the position under the static equilibrium in the fixed-length vibration, and the solid line is the position at the end of the rapid unloading stage.Figure 7
Initial state of the fixed-length vibration.As is shown in Figure6, vi(t), ai(t), and zi(t) are oscillator velocity, oscillator acceleration, and position coordinate at the time of t, respectively, yielding(8)vit=dzitdt.The length of the upper spring is(9)l1it=zit-u1it-l32.The amount of the spring compression is(10)Δx1it=l10-zit+u1it+l32.The length of the lower spring is(11)l2it=ln0+u2it-zit-l32.The amount of the spring compression is(12)Δx2it=l20-ln0-u2it+zit+l32.The oscillator acceleration yields(13)ait=k1Δx1it+mg-k2Δx2itm=k1l10-zit+u1it+l3/2+mg-k2l20-ln0-u2it+zit+l3/2m.Then the differential equation of motion for the spring oscillator is(14)d2zitdt2+k1+k2mzit=k1mu1it+k2mu2it-k2l20-ln0-mg-k1l10-k1-k2l3/2m.The above equation is an ordinary differential equation of the second order, which can be solved by the Runge-Kutta method of the fourth order. This equation can be transformed to following forms:(15)zi′′t=ft,zit,zi′t,t0≤t≤tnzit0=l1+l32,zi′t0=0.Assuming thatwit=zi′t, the above equations can be written as following ordinary differential equations of the first order(16)zi′t=wit,zit0=l1+l32,t0≤t≤tnwi′t=ft,zit,wit=k1mu1it+k2mu2it-k2l20-ln0-mg-k1l10-k1-k2l3/2m-k1+k2mzit.According to the Runge-Kutta method of the fourth order [17, 18], its numerical calculation formula can be represented as(17)zij+1t0+j+1h=zijt0+jh+h6K1+2K2+2K3+K4,wij+1t0+j+1h=wijt0+jh+h6M1+2M2+2M3+M4,where(18)K1=wijt0+jh,M1=k1mu1it0+jh+k2mu2it0+jh-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh,K2=wijt0+jh+h2M1,M2=k1mu1it0+jh+h2+k2mu2it0+jh+h2-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh+h2K1,K3=wijt0+jh+h2M2,M3=k1mu1it0+jh+h2+k2mu2it0+jh+h2-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh+h2K2,K4=wijt0+jh+hM3,M4=k1mu1it0+j+1h+k2mu2it0+j+1h-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh+hK3.Given the displacement boundary conditions of the wave converter top and bottom (u1i(t), u2i(t)), the position coordinate zi(t) and velocity vi(t) of the oscillator at the time of t can be derived. Therefore, the amount of the compression Δx1it, Δx2i(t) can be obtained.


#### 4.3.2. Differential Equation of Motion in Slow Unloading Stage

In the slow unloading stage, the vertical displacements of the wave converter top and bottom have few changes over time, so the length of the wave converter can be considered as a constant.l0′ and ln0′ are the length and clear length of the wave converter at the end of the rapid unloading stage, respectively. As is shown in Figure 7, the geometry of the wave converter at the end of the rapid unloading stage is drawn in the solid line, while the dashed line shows the static equilibrium position of the mass block. The y-axis is vertical downward and its coordinate origin is the static equilibrium position of the mass block in stage 3. In Figure 7, at the end of the rapid unloading stage, vm1 is the oscillator velocity, and ym1 is the distance from the oscillator center to that of the static equilibrium position.If the oscillator is in static equilibrium under the converter length ofl0′, ln0′ yields (19)ln0′=l10+l20+l3-Δx11-Δx21.The balance equation for the mass block is(20)mg+k1Δx11=k2Δx21,where Δx11 and Δx21 are the amount of compression for upper and lower spring in static equilibrium under the converter length of l0′, respectively. Δx11 and Δx21 are(21)Δx11=k2l10+l20+l3-ln0′-mgk1+k2,Δx21=mgk2+k1l10+l20+l3-ln0′k1+k2-k1mgk2k1+k2,F
i, fIi, and Si are the active force, inertia force, and constraint counterforce, respectively. Based on the D’Alembert principle, Fi, fIi, and Si yield(22)Fi+Si+fIi=0.The active force consists of the gravity forcemg, damping force fD, and elastic restoring force fs. Then (22) can be written as (23)my¨mt+cy˙mt+k1+k2ymt=0.ym(t) is the position coordinate of the oscillator in y coordinate. The time at the end of the rapid unloading stage is assumed as tu, and then the position coordinate and velocity of the oscillator can be expressed as zi(tu) and vi(tu) in z coordinate. Via the coordinate transform, the position coordinate and velocity of the oscillator in y coordinate can be expressed as ym1 and vm1. ym1 and vm1 are the initial conditions of (23), so the solution of (23) is(24)ymt=Acos⁡ωDt+Bsin⁡ωDt,where(25)A=ym1e-ζωnt,B=vm1+ζωnym1ωDe-ζωnt,ωD=ωn1-ζ2,ωn=k1+k2m;ζ is the damping ratio. If ζ=0, the calculation results of the dynamic response are a little dangerous which would be safer for the structure design, so ζ is considered as 0 in the analysis of stage 3. Considering that ym(t) is solved, the amount of compression for springs at any time can be obtained. After that, the stress-time curve of the wave converter top and bottom in the slow unloading stage can be gotten.


### 4.4. Displacement Transfer Coefficient of the Wave Converter

In the rapid loading stage and rapid unloading stage, based on the numerical calculations in model 1 and model 2, the following equation can be derived:(26)sPA0t/wu0it=sPA0t/w-F1Kiu0it-u2itKiu0it,where w is the area ratio of the wave converter’s cross-sectional area, and w is equal to all wave converters’ cross-sectional areas divided by the area of rock surface above the structure. According to (4) and (7), when sPA(t)≥k1Δx10, the constitutive relation for the wave converter can be modified as(27)Kiu0it-u2it=F1Kiu0it-u2itk1+k2k1k2-k1+k2l10+l20+l3-ln0-mgk1+k2+k1mgk2k1+k2.Combining (26) and (27), Ki yields(28)Ki=k1+k2sPA0t+wk1k2u2it+wk1mg-k2l10+l20+l3-ln0wk1k2u0it+k1+k2sPA0t.


### 4.5. Stress Inversion of the Wave Converter

Based on the solutions on differential equations of motion in 3 stages, if the gravity force of the wave converter shell is ignored, the stress-time curves of the wave converter top and bottom in 3 stages can be gotten.In stage 1 and stage 2, the function of the stress-time curve of the wave converter top is(29)P1it=k1Δx1its=k1l10-zit+u1it+l3/2s.The function of the stress-time curve of the wave converter bottom is(30)P2it=k2Δx2its=k2l20-ln0-u2it+zit+l3/2sIn stage 3, the functions of the stress-time curves of the wave converter top and bottom are, respectively,(31)P1it=k1Δx11-ymts=k1Δx11-k1Acos⁡ωnt+Bsin⁡ωnts,P2it=k2Δx21+ymts=k2Δx21+k2Acos⁡ωnt+Bsin⁡ωnts.Then, in model 1, the initial pressure Pin caused by the wave converter on the distribution layer surface is deleted, and the stress-time curves of the wave converter top and bottom in three stages are applied on the rock surface and distribution layer surface, respectively. After that, the blasting load is applied and the dynamic calculation is conducted to get the dynamic response of the structure (shown in calculation process (4)).


## 4.1. Static Constitutive Relation of the Wave Converter

l
0 is the initial length of the wave converter, while a is the shell thickness. So the initial clear length of the wave converter can be written as ln0=l0-2a, and ln0 is shown in Figure 6. In Figure 6, the dashed line is the position under the static equilibrium, and the solid line is the position at any time under the dynamic load. k1 is the stiffness coefficient of the upper spring, whose length in free state is l10.Δx10is the initial amount of compression for upper spring. k2is the stiffness coefficient of the lower spring, and l20 is its length in free state. Δx20is the initial amount of compression for lower spring. m and l3 are the mass and height of the mass block, respectively. PA(t) is the vertical rock pressure on the wave converter top in the true situation, while PA0(t) is the vertical rock pressure at the same location in model 2. When the wave converter is installed with the initial length l0, under static equilibrium, the geometric equation is established, yielding(1)ln0=l1+l2+l3=l10+l20+l3-Δx10-Δx20.Figure 6
Deformation process of the wave converter.As for the mass block, the balance equation is(2)mg+k1Δx10=k2Δx20.Combining (1) and (2), the solutions are(3)Δx10=k2l10+l20+l3-ln0-mgk1+k2,Δx20=mgk2+k1l10+l20+l3-ln0k1+k2-k1mgk2k1+k2.F
1 and F2 are the force of upper and lower springs, respectively. s is the cross-sectional area of the wave converter. d is the cross-sectional length. If the weight of the wave converter’s shells is ignored, F1 can be expressed as(4)F1t=k1Δx10,sPAt<k1Δx10sPAt,sPAt≥k1Δx10;F
2 can be expressed as(5)F2t=mg+k1Δx10sPAt<k1Δx10sPAt+mgsPAt≥k1Δx10.So the static constitutive relation of the wave converter yields(6)IfsPAt<k1Δx10,Δl=0(7)IfsPAt≥k1Δx10,Δl=k1+k2sPAtk1k2+mg-k2l10+l20+l3-ln0k2k1+k2.


## 4.2. Dynamic Response Partitioning of the Wave Converter

The rapid loading stage, rapid unloading stage, and slow unloading stage correspond to the rapid compression stage, rapid recovery stage, and slow recovery stage, respectively. In the slow recovery stage, the length of the wave converter can be deemed as a constant. The demarcation point of stage 1 and stage 2 is that the relative displacement reaches the maximum. The duration time of stage 2 can be determined by the relative displacement-time curve. When the rock masses above the structure are in a wide range of elasticity state, the rapid unloading stage can not be ignored, but if the rock masses are in a wide range of plastic state, the rapid unloading stage can be ignored.The computing time of stage 3 is advisable for 1 or 2 vibration periods. Via a large amount of computations, it is concluded that the computing time of stage 3 can be taken as 1 vibration period if stage 2 can not be ignored; otherwise it can be taken as 2 vibration periods.


## 4.3. Differential Equation of Motion for the Spring Oscillator

### 4.3.1. Differential Equation of Motion in Rapid Loading Stage and Rapid Unloading Stage

N is the total number of wave converters along the width direction of the structure, and i is the serial number of the wave converter. For wave converter i, u0i(t) is the vertical displacement of the rock surface at the wave converter top in model 1, and u2iti=1,2,…,n is the vertical displacement of the wave converter bottom in rapid loading stage and rapid unloading stage of the true situation. u1i(t) is the vertical displacement of the wave converter top in rapid loading stage and rapid unloading stage of the true situation, which can be expressed as u1i(t)=Kiu0i(t). Ki is the displacement transfer coefficient of the wave converter. Δui(t) is the relative displacement between the top and bottom of the wave converter, which can be written as Δui(t)=u1i(t)-u2it. At the end of the rapid unloading stage is the initial state of the slow unloading stage, which can be deemed as the fixed-length vibration, shown in Figure 7. In Figure 7, the dashed line is the position under the static equilibrium in the fixed-length vibration, and the solid line is the position at the end of the rapid unloading stage.Figure 7
Initial state of the fixed-length vibration.As is shown in Figure6, vi(t), ai(t), and zi(t) are oscillator velocity, oscillator acceleration, and position coordinate at the time of t, respectively, yielding(8)vit=dzitdt.The length of the upper spring is(9)l1it=zit-u1it-l32.The amount of the spring compression is(10)Δx1it=l10-zit+u1it+l32.The length of the lower spring is(11)l2it=ln0+u2it-zit-l32.The amount of the spring compression is(12)Δx2it=l20-ln0-u2it+zit+l32.The oscillator acceleration yields(13)ait=k1Δx1it+mg-k2Δx2itm=k1l10-zit+u1it+l3/2+mg-k2l20-ln0-u2it+zit+l3/2m.Then the differential equation of motion for the spring oscillator is(14)d2zitdt2+k1+k2mzit=k1mu1it+k2mu2it-k2l20-ln0-mg-k1l10-k1-k2l3/2m.The above equation is an ordinary differential equation of the second order, which can be solved by the Runge-Kutta method of the fourth order. This equation can be transformed to following forms:(15)zi′′t=ft,zit,zi′t,t0≤t≤tnzit0=l1+l32,zi′t0=0.Assuming thatwit=zi′t, the above equations can be written as following ordinary differential equations of the first order(16)zi′t=wit,zit0=l1+l32,t0≤t≤tnwi′t=ft,zit,wit=k1mu1it+k2mu2it-k2l20-ln0-mg-k1l10-k1-k2l3/2m-k1+k2mzit.According to the Runge-Kutta method of the fourth order [17, 18], its numerical calculation formula can be represented as(17)zij+1t0+j+1h=zijt0+jh+h6K1+2K2+2K3+K4,wij+1t0+j+1h=wijt0+jh+h6M1+2M2+2M3+M4,where(18)K1=wijt0+jh,M1=k1mu1it0+jh+k2mu2it0+jh-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh,K2=wijt0+jh+h2M1,M2=k1mu1it0+jh+h2+k2mu2it0+jh+h2-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh+h2K1,K3=wijt0+jh+h2M2,M3=k1mu1it0+jh+h2+k2mu2it0+jh+h2-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh+h2K2,K4=wijt0+jh+hM3,M4=k1mu1it0+j+1h+k2mu2it0+j+1h-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh+hK3.Given the displacement boundary conditions of the wave converter top and bottom (u1i(t), u2i(t)), the position coordinate zi(t) and velocity vi(t) of the oscillator at the time of t can be derived. Therefore, the amount of the compression Δx1it, Δx2i(t) can be obtained.


### 4.3.2. Differential Equation of Motion in Slow Unloading Stage

In the slow unloading stage, the vertical displacements of the wave converter top and bottom have few changes over time, so the length of the wave converter can be considered as a constant.l0′ and ln0′ are the length and clear length of the wave converter at the end of the rapid unloading stage, respectively. As is shown in Figure 7, the geometry of the wave converter at the end of the rapid unloading stage is drawn in the solid line, while the dashed line shows the static equilibrium position of the mass block. The y-axis is vertical downward and its coordinate origin is the static equilibrium position of the mass block in stage 3. In Figure 7, at the end of the rapid unloading stage, vm1 is the oscillator velocity, and ym1 is the distance from the oscillator center to that of the static equilibrium position.If the oscillator is in static equilibrium under the converter length ofl0′, ln0′ yields (19)ln0′=l10+l20+l3-Δx11-Δx21.The balance equation for the mass block is(20)mg+k1Δx11=k2Δx21,where Δx11 and Δx21 are the amount of compression for upper and lower spring in static equilibrium under the converter length of l0′, respectively. Δx11 and Δx21 are(21)Δx11=k2l10+l20+l3-ln0′-mgk1+k2,Δx21=mgk2+k1l10+l20+l3-ln0′k1+k2-k1mgk2k1+k2,F
i, fIi, and Si are the active force, inertia force, and constraint counterforce, respectively. Based on the D’Alembert principle, Fi, fIi, and Si yield(22)Fi+Si+fIi=0.The active force consists of the gravity forcemg, damping force fD, and elastic restoring force fs. Then (22) can be written as (23)my¨mt+cy˙mt+k1+k2ymt=0.ym(t) is the position coordinate of the oscillator in y coordinate. The time at the end of the rapid unloading stage is assumed as tu, and then the position coordinate and velocity of the oscillator can be expressed as zi(tu) and vi(tu) in z coordinate. Via the coordinate transform, the position coordinate and velocity of the oscillator in y coordinate can be expressed as ym1 and vm1. ym1 and vm1 are the initial conditions of (23), so the solution of (23) is(24)ymt=Acos⁡ωDt+Bsin⁡ωDt,where(25)A=ym1e-ζωnt,B=vm1+ζωnym1ωDe-ζωnt,ωD=ωn1-ζ2,ωn=k1+k2m;ζ is the damping ratio. If ζ=0, the calculation results of the dynamic response are a little dangerous which would be safer for the structure design, so ζ is considered as 0 in the analysis of stage 3. Considering that ym(t) is solved, the amount of compression for springs at any time can be obtained. After that, the stress-time curve of the wave converter top and bottom in the slow unloading stage can be gotten.


## 4.3.1. Differential Equation of Motion in Rapid Loading Stage and Rapid Unloading Stage

N is the total number of wave converters along the width direction of the structure, and i is the serial number of the wave converter. For wave converter i, u0i(t) is the vertical displacement of the rock surface at the wave converter top in model 1, and u2iti=1,2,…,n is the vertical displacement of the wave converter bottom in rapid loading stage and rapid unloading stage of the true situation. u1i(t) is the vertical displacement of the wave converter top in rapid loading stage and rapid unloading stage of the true situation, which can be expressed as u1i(t)=Kiu0i(t). Ki is the displacement transfer coefficient of the wave converter. Δui(t) is the relative displacement between the top and bottom of the wave converter, which can be written as Δui(t)=u1i(t)-u2it. At the end of the rapid unloading stage is the initial state of the slow unloading stage, which can be deemed as the fixed-length vibration, shown in Figure 7. In Figure 7, the dashed line is the position under the static equilibrium in the fixed-length vibration, and the solid line is the position at the end of the rapid unloading stage.Figure 7
Initial state of the fixed-length vibration.As is shown in Figure6, vi(t), ai(t), and zi(t) are oscillator velocity, oscillator acceleration, and position coordinate at the time of t, respectively, yielding(8)vit=dzitdt.The length of the upper spring is(9)l1it=zit-u1it-l32.The amount of the spring compression is(10)Δx1it=l10-zit+u1it+l32.The length of the lower spring is(11)l2it=ln0+u2it-zit-l32.The amount of the spring compression is(12)Δx2it=l20-ln0-u2it+zit+l32.The oscillator acceleration yields(13)ait=k1Δx1it+mg-k2Δx2itm=k1l10-zit+u1it+l3/2+mg-k2l20-ln0-u2it+zit+l3/2m.Then the differential equation of motion for the spring oscillator is(14)d2zitdt2+k1+k2mzit=k1mu1it+k2mu2it-k2l20-ln0-mg-k1l10-k1-k2l3/2m.The above equation is an ordinary differential equation of the second order, which can be solved by the Runge-Kutta method of the fourth order. This equation can be transformed to following forms:(15)zi′′t=ft,zit,zi′t,t0≤t≤tnzit0=l1+l32,zi′t0=0.Assuming thatwit=zi′t, the above equations can be written as following ordinary differential equations of the first order(16)zi′t=wit,zit0=l1+l32,t0≤t≤tnwi′t=ft,zit,wit=k1mu1it+k2mu2it-k2l20-ln0-mg-k1l10-k1-k2l3/2m-k1+k2mzit.According to the Runge-Kutta method of the fourth order [17, 18], its numerical calculation formula can be represented as(17)zij+1t0+j+1h=zijt0+jh+h6K1+2K2+2K3+K4,wij+1t0+j+1h=wijt0+jh+h6M1+2M2+2M3+M4,where(18)K1=wijt0+jh,M1=k1mu1it0+jh+k2mu2it0+jh-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh,K2=wijt0+jh+h2M1,M2=k1mu1it0+jh+h2+k2mu2it0+jh+h2-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh+h2K1,K3=wijt0+jh+h2M2,M3=k1mu1it0+jh+h2+k2mu2it0+jh+h2-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh+h2K2,K4=wijt0+jh+hM3,M4=k1mu1it0+j+1h+k2mu2it0+j+1h-k2l20-l0-mg-k1l10m-k1+k2mzijt0+jh+hK3.Given the displacement boundary conditions of the wave converter top and bottom (u1i(t), u2i(t)), the position coordinate zi(t) and velocity vi(t) of the oscillator at the time of t can be derived. Therefore, the amount of the compression Δx1it, Δx2i(t) can be obtained.


## 4.3.2. Differential Equation of Motion in Slow Unloading Stage

In the slow unloading stage, the vertical displacements of the wave converter top and bottom have few changes over time, so the length of the wave converter can be considered as a constant.l0′ and ln0′ are the length and clear length of the wave converter at the end of the rapid unloading stage, respectively. As is shown in Figure 7, the geometry of the wave converter at the end of the rapid unloading stage is drawn in the solid line, while the dashed line shows the static equilibrium position of the mass block. The y-axis is vertical downward and its coordinate origin is the static equilibrium position of the mass block in stage 3. In Figure 7, at the end of the rapid unloading stage, vm1 is the oscillator velocity, and ym1 is the distance from the oscillator center to that of the static equilibrium position.If the oscillator is in static equilibrium under the converter length ofl0′, ln0′ yields (19)ln0′=l10+l20+l3-Δx11-Δx21.The balance equation for the mass block is(20)mg+k1Δx11=k2Δx21,where Δx11 and Δx21 are the amount of compression for upper and lower spring in static equilibrium under the converter length of l0′, respectively. Δx11 and Δx21 are(21)Δx11=k2l10+l20+l3-ln0′-mgk1+k2,Δx21=mgk2+k1l10+l20+l3-ln0′k1+k2-k1mgk2k1+k2,F
i, fIi, and Si are the active force, inertia force, and constraint counterforce, respectively. Based on the D’Alembert principle, Fi, fIi, and Si yield(22)Fi+Si+fIi=0.The active force consists of the gravity forcemg, damping force fD, and elastic restoring force fs. Then (22) can be written as (23)my¨mt+cy˙mt+k1+k2ymt=0.ym(t) is the position coordinate of the oscillator in y coordinate. The time at the end of the rapid unloading stage is assumed as tu, and then the position coordinate and velocity of the oscillator can be expressed as zi(tu) and vi(tu) in z coordinate. Via the coordinate transform, the position coordinate and velocity of the oscillator in y coordinate can be expressed as ym1 and vm1. ym1 and vm1 are the initial conditions of (23), so the solution of (23) is(24)ymt=Acos⁡ωDt+Bsin⁡ωDt,where(25)A=ym1e-ζωnt,B=vm1+ζωnym1ωDe-ζωnt,ωD=ωn1-ζ2,ωn=k1+k2m;ζ is the damping ratio. If ζ=0, the calculation results of the dynamic response are a little dangerous which would be safer for the structure design, so ζ is considered as 0 in the analysis of stage 3. Considering that ym(t) is solved, the amount of compression for springs at any time can be obtained. After that, the stress-time curve of the wave converter top and bottom in the slow unloading stage can be gotten.


## 4.4. Displacement Transfer Coefficient of the Wave Converter

In the rapid loading stage and rapid unloading stage, based on the numerical calculations in model 1 and model 2, the following equation can be derived:(26)sPA0t/wu0it=sPA0t/w-F1Kiu0it-u2itKiu0it,where w is the area ratio of the wave converter’s cross-sectional area, and w is equal to all wave converters’ cross-sectional areas divided by the area of rock surface above the structure. According to (4) and (7), when sPA(t)≥k1Δx10, the constitutive relation for the wave converter can be modified as(27)Kiu0it-u2it=F1Kiu0it-u2itk1+k2k1k2-k1+k2l10+l20+l3-ln0-mgk1+k2+k1mgk2k1+k2.Combining (26) and (27), Ki yields(28)Ki=k1+k2sPA0t+wk1k2u2it+wk1mg-k2l10+l20+l3-ln0wk1k2u0it+k1+k2sPA0t.


## 4.5. Stress Inversion of the Wave Converter

Based on the solutions on differential equations of motion in 3 stages, if the gravity force of the wave converter shell is ignored, the stress-time curves of the wave converter top and bottom in 3 stages can be gotten.In stage 1 and stage 2, the function of the stress-time curve of the wave converter top is(29)P1it=k1Δx1its=k1l10-zit+u1it+l3/2s.The function of the stress-time curve of the wave converter bottom is(30)P2it=k2Δx2its=k2l20-ln0-u2it+zit+l3/2sIn stage 3, the functions of the stress-time curves of the wave converter top and bottom are, respectively,(31)P1it=k1Δx11-ymts=k1Δx11-k1Acos⁡ωnt+Bsin⁡ωnts,P2it=k2Δx21+ymts=k2Δx21+k2Acos⁡ωnt+Bsin⁡ωnts.Then, in model 1, the initial pressure Pin caused by the wave converter on the distribution layer surface is deleted, and the stress-time curves of the wave converter top and bottom in three stages are applied on the rock surface and distribution layer surface, respectively. After that, the blasting load is applied and the dynamic calculation is conducted to get the dynamic response of the structure (shown in calculation process (4)).


## 5. Case Study

### 5.1. Case Set-Up

Based on the software ofFLAC3D, 2 kinds of cases are set up to conduct the dynamic analysis. Case  1 is the traditional underground structure without the new blast-resistant method, while case  2 is the underground structure with the mentioned new blast-resistant method. Compared with case  2, the difference in case  1 is that the wave converter and distribution layer are not set.Figure8 shows the structure size and monitoring sections for internal forces in case  2, and the monitoring sections for internal forces are the same as case  1. The size of the structure is designed according to literature [19], and the buried depth is 10 m. Mechanical parameters for the rock and structure are listed in Tables 1 and 2. The distribution layer is made up of the foam concrete which has the density of 799 kg/m3 and thickness of 0.05 m. The mechanical parameters of the foam concrete are shown in Table 3 [20]. In case  2, the tops and bottoms of wave converters 1~5 are chosen as monitoring positions for displacements and stresses, shown in Figure 8, and the interpolation is used to get the stresses of other wave converters. The wave converter is a cubic structure with the side length of 0.2 m. l0=0.2 m, a=0.05 m, w=100%, k1=k2=100 kN/m, l10=l20=0.05 m. l3=0.14 m, m=38.9 kg.Table 1
Mechanical parameters for the rock.
Density (kg/m3)
Elastic modulus (GPa)
Poisson ratio
Cohesion (MPa)
Internal friction angle (°)
Tensile strength (MPa)
2400
13
0.28
1.1
45
0.91Table 2
Mechanical parameters for the structure.
Density (kg/m3)
Elastic modulus (GPa)
Poisson ratio
Cohesion (MPa)
Internal friction angle (°)
2400
30
0.21
3.68
58.7Table 3
Mechanical parameters for foam concrete.
Density (kg/m3)
Elastic modulus (GPa)
Poisson ratio
Cohesion (MPa)
Internal friction angle (°)
Tensile strength (MPa)
799
0.342
0.1
0.17
29
0.2Figure 8
Structure size and monitoring sections for internal forces in case  1.The width, height, and thickness of the numerical models are 35.2 m, 28.85 m, and 1 m, respectively. The blasting load is assumed as a triangle wave acting on the ground surface (in Figure9), and the loading scope is from −3 m to 3 m on x-axis. The lifting duration and drop duration of the blast pressure-time curve are set as 1 ms and 6 ms according to the literature [21], respectively. The peak of the shock wave pressure pm is 0.6 GPa. Model 1 built up based on FLAC3D for numerical calculation is shown in Figure 10.Figure 9
Curve of the blasting load.Figure 10
Model 1 built up based onFLAC3D for numerical calculation.According to the symmetry, the monitoring positions for displacements are suggested to adopt the tops and bottoms of 5 wave converters along the width direction of the structure, shown in Figure8. In order to obtain the internal forces such as the bending moment, axial force, and shear force of the structure, a total of 15 monitoring sections are set up, shown in Figure 8.


### 5.2. Analysis of the Calculation Results

The curves of the displacement transfer coefficientKi with time in the rapid loading and rapid unloading stages are shown in Figure 11, and i represents the serial number of the wave converter. The curves of u0i(t), u1i(t), u2i(t), and Δui(t) with time are shown in Figures 12–15, respectively. In all 3 stages of the wave converter deformation, the stress-time curves of the wave converter top are shown in Figure 16, while the stress-time curves of the wave converter bottom are shown in Figure 17. The internal forces such as the bending moment, axial force, and shear force are listed in Tables 4–6.Table 4
Internal forces of the structure for monitoring sections  1–5.
Monitoring section
1
2
3
4
5
Case
1
2
1
2
1
2
1
2
1
2
Bending moment (N·m)
Minimum
733
−43847
640
−45141
365
−56001
243
−92795
−290760
−129805
Maximum
221068
163402
224395
148737
222567
95049
118881
8908
348
2039
Shear force (N)
Minimum
−27317
−15193
2360
−3988
1194
−1956
1177
0
3474
23
Maximum
−426
1037
164052
64656
439202
138742
1041408
214362
2040240
282142
Axial force (N)
Minimum
−41094
−83637
−42568
−81050
−57842
−81420
−104766
−83599
−348627
−91759
Maximum
1666072
921440
1596560
921280
1320330
921600
788366
932400
409633
949376Table 5
Internal forces of the structure for monitoring sections 6–10.
Monitoring section
6
7
8
9
10
Case
1
2
1
2
1
2
1
2
1
2
Bending moment (N·m)
Minimum
−246269
−74504
−60566
−34165
−20880
−16684
−20016
−12038
−30226
−27187
Maximum
−62
59448
11794
32440
20808
19901
4810
3787
10929
8644
Shear force (N)
Minimum
−837708
−99838
−371508
−85756
−165964
−58982
−78672
−57528
−100307
−77861
Maximum
−1537
211056
−1517
90487
−953
68296
27172
105168
51451
117882
Axial force (N)
Minimum
−3218640
−814500
−2655720
−846588
−2229240
−891240
−1992840
−936840
−1573560
−834600
Maximum
−4890
−1196
−1170
−2138
−1175
−1223
−620
−1812
−5041
−983Table 6
Internal forces of the structure for monitoring sections 11–15.
Monitoring section
11
12
13
14
15
Case
1
2
1
2
1
2
1
2
1
2
Bending moment (N·m)
Minimum
230
664
261
614
295
555
97
419
−54210
−45028
Maximum
35258
23027
36594
23699
39747
24339
32743
15492
15725
14552
Shear force (N)
Minimum
−421
−92
−19776
−16431
−65196
−51572
−212922
−153744
−533254
−356511
Maximum
4128
3608
319
411
−428
450
2836
718
5536
6486
Axial force (N)
Minimum
−3124
−4749
−3185
−5153
−2953
−18050
−32879
−115840
−192165
−250929
Maximum
680160
426831
667594
419751
615300
394217
518211
354017
381111
332100Figure 11
Curves of the displacement transfer coefficient with time.Figure 12
Curves ofu0i(t) with time.Figure 13
Curves ofu1i(t) with time.Figure 14
Curves ofu2i(t) with time.Figure 15
Curves ofΔui(t) with time.Figure 16
Stress-time curves of the wave converter top.Figure 17
Stress-time curves of the wave converter bottom.In Tables4–6, the bending moment resulting in the tensile stress in the inner element of the structure is positive. The positive shear force is by counterclockwise, while the axial force to tension is positive. It can be concluded that the peak absolute value of the bending moment in case  2 is generally lower than that of case  1, and the maximum drop in the roof, side wall, and floor is 57.3%, 69.7%, and 52.7%, respectively.The peak absolute value of the shear force for monitoring sections in case  2 is also generally lower than that of case  1. The maximum drop in the roof, side wall, and floor is 86.2%, 75.6%, and 33.1%, respectively. The peak absolute value of the axial force for monitoring sections in case  2 is remarkably lower than that of case  1. The maximum drop in the roof, side wall, and floor is 44.7%, 74.7%, and 37.2%, respectively. For the roof, the decrease of the axial tensile force near the midspan is obvious, and the maximum drop occurs to the span center. The decrease of the axial tensile stress presents that the tensile failure in the roof can be alleviated via the wave converters.The curves of the vertical normal stresses for inner and outer elements in the span center with time are shown in Figure18. The positive value stands for the vertical tensile stress, while the negative value stands for the vertical compressive stress. After the adoption of wave converters, the peak vertical tensile stress for inner element in the span center drops from 0.055 MPa to 0.019 MPa, which means that the possibility of spalling damage for roof is reduced. The peak vertical compressive stress for outer element in the span center drops from 1.73 MPa to 0.153 MPa, which means that the impact load acting on the roof is reduced.Figure 18
Vertical normal stress of inner and outer elements in the span center.
(a)
Inner element
(b)
Outer elementThe peak horizontal tensile stresses of monitoring sections for roof in 2 cases are shown in Figure19. For case  2 with wave converters, the peak horizontal tensile stresses of inner elements for roof are generally lower than that of case  1, with a maximum decrease of 60%. Though the peak horizontal tensile stresses of outer elements for roof become higher than that of case  1, the peak horizontal tensile stresses after increasing are not very large on the whole. Via increasing the spring stiffness, the increase of the peak horizontal tensile stresses of outer elements for roof can be adjusted. In the design process of wave converters, the vibration period of the spring oscillator should not be next to the vibration period of the underground structure.Figure 19
Peak horizontal tensile stress of monitoring sections for roof in 2 cases.


## 5.1. Case Set-Up

Based on the software ofFLAC3D, 2 kinds of cases are set up to conduct the dynamic analysis. Case  1 is the traditional underground structure without the new blast-resistant method, while case  2 is the underground structure with the mentioned new blast-resistant method. Compared with case  2, the difference in case  1 is that the wave converter and distribution layer are not set.Figure8 shows the structure size and monitoring sections for internal forces in case  2, and the monitoring sections for internal forces are the same as case  1. The size of the structure is designed according to literature [19], and the buried depth is 10 m. Mechanical parameters for the rock and structure are listed in Tables 1 and 2. The distribution layer is made up of the foam concrete which has the density of 799 kg/m3 and thickness of 0.05 m. The mechanical parameters of the foam concrete are shown in Table 3 [20]. In case  2, the tops and bottoms of wave converters 1~5 are chosen as monitoring positions for displacements and stresses, shown in Figure 8, and the interpolation is used to get the stresses of other wave converters. The wave converter is a cubic structure with the side length of 0.2 m. l0=0.2 m, a=0.05 m, w=100%, k1=k2=100 kN/m, l10=l20=0.05 m. l3=0.14 m, m=38.9 kg.Table 1
Mechanical parameters for the rock.
Density (kg/m3)
Elastic modulus (GPa)
Poisson ratio
Cohesion (MPa)
Internal friction angle (°)
Tensile strength (MPa)
2400
13
0.28
1.1
45
0.91Table 2
Mechanical parameters for the structure.
Density (kg/m3)
Elastic modulus (GPa)
Poisson ratio
Cohesion (MPa)
Internal friction angle (°)
2400
30
0.21
3.68
58.7Table 3
Mechanical parameters for foam concrete.
Density (kg/m3)
Elastic modulus (GPa)
Poisson ratio
Cohesion (MPa)
Internal friction angle (°)
Tensile strength (MPa)
799
0.342
0.1
0.17
29
0.2Figure 8
Structure size and monitoring sections for internal forces in case  1.The width, height, and thickness of the numerical models are 35.2 m, 28.85 m, and 1 m, respectively. The blasting load is assumed as a triangle wave acting on the ground surface (in Figure9), and the loading scope is from −3 m to 3 m on x-axis. The lifting duration and drop duration of the blast pressure-time curve are set as 1 ms and 6 ms according to the literature [21], respectively. The peak of the shock wave pressure pm is 0.6 GPa. Model 1 built up based on FLAC3D for numerical calculation is shown in Figure 10.Figure 9
Curve of the blasting load.Figure 10
Model 1 built up based onFLAC3D for numerical calculation.According to the symmetry, the monitoring positions for displacements are suggested to adopt the tops and bottoms of 5 wave converters along the width direction of the structure, shown in Figure8. In order to obtain the internal forces such as the bending moment, axial force, and shear force of the structure, a total of 15 monitoring sections are set up, shown in Figure 8.


## 5.2. Analysis of the Calculation Results

The curves of the displacement transfer coefficientKi with time in the rapid loading and rapid unloading stages are shown in Figure 11, and i represents the serial number of the wave converter. The curves of u0i(t), u1i(t), u2i(t), and Δui(t) with time are shown in Figures 12–15, respectively. In all 3 stages of the wave converter deformation, the stress-time curves of the wave converter top are shown in Figure 16, while the stress-time curves of the wave converter bottom are shown in Figure 17. The internal forces such as the bending moment, axial force, and shear force are listed in Tables 4–6.Table 4
Internal forces of the structure for monitoring sections  1–5.
Monitoring section
1
2
3
4
5
Case
1
2
1
2
1
2
1
2
1
2
Bending moment (N·m)
Minimum
733
−43847
640
−45141
365
−56001
243
−92795
−290760
−129805
Maximum
221068
163402
224395
148737
222567
95049
118881
8908
348
2039
Shear force (N)
Minimum
−27317
−15193
2360
−3988
1194
−1956
1177
0
3474
23
Maximum
−426
1037
164052
64656
439202
138742
1041408
214362
2040240
282142
Axial force (N)
Minimum
−41094
−83637
−42568
−81050
−57842
−81420
−104766
−83599
−348627
−91759
Maximum
1666072
921440
1596560
921280
1320330
921600
788366
932400
409633
949376Table 5
Internal forces of the structure for monitoring sections 6–10.
Monitoring section
6
7
8
9
10
Case
1
2
1
2
1
2
1
2
1
2
Bending moment (N·m)
Minimum
−246269
−74504
−60566
−34165
−20880
−16684
−20016
−12038
−30226
−27187
Maximum
−62
59448
11794
32440
20808
19901
4810
3787
10929
8644
Shear force (N)
Minimum
−837708
−99838
−371508
−85756
−165964
−58982
−78672
−57528
−100307
−77861
Maximum
−1537
211056
−1517
90487
−953
68296
27172
105168
51451
117882
Axial force (N)
Minimum
−3218640
−814500
−2655720
−846588
−2229240
−891240
−1992840
−936840
−1573560
−834600
Maximum
−4890
−1196
−1170
−2138
−1175
−1223
−620
−1812
−5041
−983Table 6
Internal forces of the structure for monitoring sections 11–15.
Monitoring section
11
12
13
14
15
Case
1
2
1
2
1
2
1
2
1
2
Bending moment (N·m)
Minimum
230
664
261
614
295
555
97
419
−54210
−45028
Maximum
35258
23027
36594
23699
39747
24339
32743
15492
15725
14552
Shear force (N)
Minimum
−421
−92
−19776
−16431
−65196
−51572
−212922
−153744
−533254
−356511
Maximum
4128
3608
319
411
−428
450
2836
718
5536
6486
Axial force (N)
Minimum
−3124
−4749
−3185
−5153
−2953
−18050
−32879
−115840
−192165
−250929
Maximum
680160
426831
667594
419751
615300
394217
518211
354017
381111
332100Figure 11
Curves of the displacement transfer coefficient with time.Figure 12
Curves ofu0i(t) with time.Figure 13
Curves ofu1i(t) with time.Figure 14
Curves ofu2i(t) with time.Figure 15
Curves ofΔui(t) with time.Figure 16
Stress-time curves of the wave converter top.Figure 17
Stress-time curves of the wave converter bottom.In Tables4–6, the bending moment resulting in the tensile stress in the inner element of the structure is positive. The positive shear force is by counterclockwise, while the axial force to tension is positive. It can be concluded that the peak absolute value of the bending moment in case  2 is generally lower than that of case  1, and the maximum drop in the roof, side wall, and floor is 57.3%, 69.7%, and 52.7%, respectively.The peak absolute value of the shear force for monitoring sections in case  2 is also generally lower than that of case  1. The maximum drop in the roof, side wall, and floor is 86.2%, 75.6%, and 33.1%, respectively. The peak absolute value of the axial force for monitoring sections in case  2 is remarkably lower than that of case  1. The maximum drop in the roof, side wall, and floor is 44.7%, 74.7%, and 37.2%, respectively. For the roof, the decrease of the axial tensile force near the midspan is obvious, and the maximum drop occurs to the span center. The decrease of the axial tensile stress presents that the tensile failure in the roof can be alleviated via the wave converters.The curves of the vertical normal stresses for inner and outer elements in the span center with time are shown in Figure18. The positive value stands for the vertical tensile stress, while the negative value stands for the vertical compressive stress. After the adoption of wave converters, the peak vertical tensile stress for inner element in the span center drops from 0.055 MPa to 0.019 MPa, which means that the possibility of spalling damage for roof is reduced. The peak vertical compressive stress for outer element in the span center drops from 1.73 MPa to 0.153 MPa, which means that the impact load acting on the roof is reduced.Figure 18
Vertical normal stress of inner and outer elements in the span center.
(a)
Inner element
(b)
Outer elementThe peak horizontal tensile stresses of monitoring sections for roof in 2 cases are shown in Figure19. For case  2 with wave converters, the peak horizontal tensile stresses of inner elements for roof are generally lower than that of case  1, with a maximum decrease of 60%. Though the peak horizontal tensile stresses of outer elements for roof become higher than that of case  1, the peak horizontal tensile stresses after increasing are not very large on the whole. Via increasing the spring stiffness, the increase of the peak horizontal tensile stresses of outer elements for roof can be adjusted. In the design process of wave converters, the vibration period of the spring oscillator should not be next to the vibration period of the underground structure.Figure 19
Peak horizontal tensile stress of monitoring sections for roof in 2 cases.


## 6. Conclusions

In this paper, a new blast-resistant method based on wave converters with spring oscillator for underground structures is put forward. The conception and calculation process of this new method are introduced. The mechanical characteristics and motion evolution law of the wave converter are derived. The dynamic responses of the traditional underground structure and the new blast-resistant one are also calculated to verify the blast-resistant effect of the new method. The following conclusions can be drawn through the study.(1)
After the deployment of wave converters, the peak absolute values of the bending moment, shear force, and axial force decrease generally. The decrease of the peak internal forces means that smaller size and less steel are needed in the design of the structure, which could help reduce the costs.(2)
After the adoption of wave converters, the peak vertical tensile stress for inner element and the peak vertical compressive stress for outer element in the span center drop remarkably, which means that the possibility of spalling damage for roof is reduced, and the impact load acting on the roof is decreased.(3)
With wave converters, the peak horizontal tensile stresses of inner elements for roof are generally lower than that of the traditional structure, which could reduce the amount of reinforcing bars.



---

*Source: 2014726-2017-05-10.xml*