Numerical Analysis of a High-Velocity Projectile’s Impact on Shallow Steel Tunnels in Soft Sandstone

Abstract

Tunnels are underground infrastructures intended for diverse community applications as well as military applications. During impact loading due to high-velocity projectiles such as ballistic missiles, materials experience a high strain rate. Moreover, there is a superficial augmentation of the dynamic strength when geomaterials such as rock are subjected to a high strain rate. Despite this strength enhancement, tunnels can get damaged by the impact load of a projectile hitting at a high velocity if they are present at a shallow depth. The present study is an effort to comprehend the response of a shallow tunnel in soft sandstone due to the impact load by a ballistic projectile using the FEM-based software ABAQUS/CAE 6.11. The Drucker–Prager damage model and the Johnson–Cook damage model were used to define the properties of the rock mass and steel tunnel lining, respectively. The crown of the 3 m diameter tunnel was kept at different depths from 1 m to 5 m from the surface. A striking velocity of 1000 m/s at a normal position to the target was given to the projectile. The projectile caused noticeable damage to the tunnel lining up to 3 m crown depth. Increasing the crown depth had a positive effect on the maximum depth of the projectile penetration up to 4 m tunnel crown depth, after which the effect reversed, making the tunnel safer. The maximum von Mises stress on the tunnel lining reduced in a logarithmic trend with an increase in the crown depth, gradually lowering to an impact load lesser than the yield stress of the tunnel lining material after a crown depth of 4.5 m.

1. Introduction

Tunnels are underground infrastructures or geostructures intended for diverse community applications such as transportation, conveyance of cables (electricity, telephone), sewers, recreation, etc., as well as military applications such as strategic defense, storage, etc. Tunnels are considered shallow when the tunnel crown is close to the ground surface [1]. Ballistic threats are those that intend to damage the target by penetration (e.g., missile), striking at a high velocity. The possibility of damage to a buried structure by such a round depends on its ability to penetrate overlying soil or rock [2]. Zukas [3] categorized ballistic threats into kinetic energy threats (small armor projectiles, high caliber rounds) and chemical energy threats (shaped charges, explosively formed projectiles). Shallow tunnels at risk of kinetic energy threats such as high-velocity ballistic missiles, during emergency situations such as war, terrorist attack, etc., will be subjected to impact load during the strike, affecting the stability of the tunnel. It is also possible to blast inside a tunnel after penetrating it during such an event by applying chemical energy threats [4].

Impact load is a class of non-periodic dynamic loads applied for a short duration of time [5]. Dynamic events such as impact and blast are characterized by transient responses in terms of stress and strain states. The duration of the event plays a significant role in these processes. As the duration of the event gets smaller (i.e., increasing strain rate), the response of the material diverges from the quasi-static behavior to dynamic equilibrium [3].

Materials experiencing impact load undergo a high strain rate, and geomaterials such as rock and concrete subjected to a high strain rate show a superficial augmentation in dynamic strength [6]. Therefore, the dynamic strength increase in geomaterials has drawn enormous concern in structural design under impact loading conditions. Taking into account the dynamic strength variations of geomaterials when designing structures subjected to impact loading conditions can enhance safety and cost-effectiveness. However, despite this strength improvement, shallow underground structures may still be vulnerable to damage from the impact force of a high-speed ballistic projectile. Therefore, it is essential to evaluate the optimum tunnel crown depth for the safety of a shallow tunnel subjected to impact load.

Examples of experimental, analytical, and numerical studies related to impact loading on different materials (metals, concrete, soil, rock, and rubber) individually can be found in the literature, although these studies conducted in the past related to impact loading on metals and geomaterials separately [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]. At present, there is a rise in experimental and numerical studies associated with tunnels subjected to impact loading [4,24,25,26,27,28,29,30,31,32,33]. However, these studies are mostly related to the study of crack propagation, deformation on the tunnel lining, and rock mass properties. Therefore, literature related to the influence of crown depth on the impact load on such geo-structures is not yet extensive.

The availability of sophisticated software and high-speed computers has made it easier to interpret complex rock mass behavior with the help of numerical modeling. The numerical methods in rock engineering are classified either as continuum or as discontinuum. A continuum model, as the name suggests, can be used for numerical analysis when there is enough evidence that the actual physical system of the rock mass can be represented in a continuum. Continuum models will be representative only if the rock mass in the field conditions acts like a continuum, i.e., only a few fractures or joints are present in the system. FEM is a numerical modeling technique based on continuum models.

The FEM approach involves dividing a complex region into discrete elements that are represented by simple geometric shapes. Mathematical equations express the properties and governing relationships of the elements, with unknown values at specific points known as nodes. These individual elements are then linked together using an assembly process to form the complete system. Load and boundary conditions are taken into account to obtain a set of linear or non-linear algebraic equations. The solution of these equations yields an approximate representation of the system’s behavior. The stress–strain behavior of the elements is approximated by constitutive relationships. FEM is the first numerical method with enough flexibility for the treatment of material heterogeneity and non-linear deformability, mainly plasticity. Among the various software packages available for numerical modeling, ABAQUS/CAE does not have predefined material properties, and thus, it can be defined independently. For this reason, ABAQUS/CAE was selected for use in this study. Preprocessing, simulation, and postprocessing are the three major components of a complete ABAQUS analysis (Figure 1). Preprocessing involves defining the model of the physical problem and creating the input file. The simulation stage runs in the background, solving the problem; postprocessing involves the display and evaluation of the results in a suitable viewer.

The present study is an effort to comprehend the response of shallow steel tunnels in soft sandstone due to the impact load by a ballistic projectile using FEM-based software ABAQUS/CAE 6.11. A 3 m diameter tunnel with the crown placed at different depths (from 1 m up to 5 m) from the rock surface. A striking velocity of 1000 m/s at a normal position to the target was given to the non-deformable projectile. The depth of penetration, velocity of the projectile, and impact load/von Mises stress on the tunnel lining corresponding to different crown depths were evaluated.

2. Methodology

2.1. Geometric Properties

The geometry of the projectile, tunnel lining, and surrounding rock were created in their local coordinates in the part module. The part module allowed us to create, edit, manage the parts, define features (e.g., solids, shells, rounds, etc.), and add reference points to a rigid part.

2.1.1. Projectile Geometry

A solid, non-deformable steel projectile, as shown in Figure 2, was modelled with a single node reference point to assign mass and initial velocity for the analysis. A striking velocity of 1000 m/s at a normal position to the target was assigned to the projectile.

2.1.2. Tunnel and Surrounding Rock Geometry

Steel tunnels in soft sandstone were considered with different crown depths from 1 m to 5 m below the overburden rock. The tunnel lining material was considered to be 50 mm thick steel. The tunnel section was considered to be 3 m in diameter. The vertical side boundaries of the soft sandstone were kept at a 2 m distance from the tunnel center. The upper half part of the tunnel and surrounding rock assembly were considered for the analysis. A schematic of the tunnel at 1 m crown depth with the surrounding rock mass is shown in Figure 3.

2.2. Material Properties

2.2.1. Constitutive Damage Models

The Drucker–Prager hardening plasticity model was used to define the damage model for the soft sandstone and the Johnson–Cook damage model was used to define the damage model for the steel in their respective material properties in the library.

Johnson–Cook Damage Model

The Johnson–Cook damage model is an elasto-viscoplastic material model proposed by Johnson and Cook [35] to predict the flow and fracture behavior of different metals under a wide range of strain, high strain rate, high temperature, and high pressures. The expression for the damage model is given by

$$\sigma =\left[A+B{\epsilon }^n\right]\left[1+Cln{\dot{\epsilon }}^{*}\right]\left[1-T^{*m}\right]$$

(1)

where σ is von Mises stress, ε is equivalent plastic strain, $\dot{\epsilon }$* is the dimensionless plastic strain rate, T* is homologous or current temperature, A is the yield stress constant, B and n are strain hardening constants, C is the viscous effect constant, and m is the thermal softening constant.

The general expression for the strain rate at fracture is given by

$${ϵ}^f=\left[D_1+D_{2}exp{D}_3{\sigma }^{*}\right]\left[1+D_{4}ln{\dot{ϵ}}^{*}\right]\left[1+D_5{T}^{*}\right]$$

(2)

where D₁, D₂, D₃, D₄, and D₅ are fracture strain constants; σ* is dimensionless pressure-stress ratio.

Drucker–Prager Yield Criterion

The Drucker–Prager yield criterion (Figure 4) is an isotropic smooth approximation of the Mohr–Columb yield criterion proposed by Drucker and Prager used for geomaterials [36]. Mathematically,

$$\left(\sigma \right)=\alpha I_1+ \sqrt{J_2}-k=0$$

(3)

where I₁ = first invariant of stress tensor = ${\sigma }_1^{\prime }+{\sigma }_2^{\prime }+{\sigma }_3^{\prime }$; J₂ = second invariant of stress tensor = $\frac{1}{6}\left[{\left({\sigma }_1^{\prime }-{\sigma }_2^{\prime }\right)}^2+{\left({\sigma }_2^{\prime }-{\sigma }_3^{\prime }\right)}^2+{\left({\sigma }_3^{\prime }-{\sigma }_1^{\prime }\right)}^2\right]$; ${\sigma }_1^{\prime }$, ${\sigma }_2^{\prime },$ and ${\sigma }_3^{\prime }$ are the principal effective stresses; α and k are material constants dependent upon the internal friction angle, $\varphi$ and cohesion, c defined as

$$\alpha =\frac{2sin\varphi }{\surd 3(3-sin\varphi )}$$

(4)

$$k=\frac{6ccos\varphi }{\surd 3(3-sin\varphi )}$$

(5)

This elasto-plastic failure criterion becomes a modification of the von Mises failure criterion when α = 0 [37]. The Drucker–Prager yield criterion tends to overestimate the strength of intact rock, which makes it suitable for defining the strength augmentation of the soft sandstone in the event of impact loading. For this reason, the present study opted for the utilization of this criterion to define the rock mass properties [6,36,37].

2.2.2. Material Properties

The material properties were created for soft sandstone and steel lining in the material library. For soft sandstone, the material properties given in

Table 1. Material parameters for soft sandstone. Reproduced from [10] published by DTIC (USA), 1992.

Density, ρ	Young’s Modulus, E	Poisson’s Ratio, υ	Angle of Friction, ϕ	Flow Stress Ratio	Dilation Angle, ψ
2291 kg/m3	13.6 × 109 N/m2	0.35	30°	0.778	3°

were used. The properties used for steel are given in

Table 2. Material parameters of steel tunnel lining. Reproduced with permission from [39] published by Elsevier, 2002.

Parameter		Value
Modulus of elasticity, E		200 × 109 N/m2
Poisson’s ratio, υ		0.3
Density, ρ		7850 kg/m3
Yield stress constant, A		490 × 106 N/m2
Strain hardening constants	B	807 × 106 N/m2
	n	0.73
Viscous effect constant, C		0.001
Thermal softening constant, m		0.94
Reference strain rate, ε0		1.000 s−1
Melting temperature, θmelt		1800 K
Transition temperature, θtransition		293 K
Fracture strain constants,	D1	0.0705
	D2	1.732
	D3	0.54
	D4	−0.015

. The stress–strain parameters and dynamic increase factor (DIF) values for soft sandstone were adopted from Chakraborty [38]. The projectile was assumed to be a rigid non-deformable steel projectile with a weight of 33.665 kg. The material properties from the library were assigned to different parts after creating sections in the property module.

2.3. Boundary Conditions and Loading

After defining the material properties, all the parts created were assembled to fix their position in the global coordinate system in the assembly module. Subsequently, an initial step and an explicit step with a 0.01 time period (10 ms) were created in the step module, and boundary conditions were applied to all the parts. Explicit dynamic analysis is a superior and efficient computation technique for evaluating the dynamic responses of simulations of considerable-sized models within a short time span. This procedure performs a large number of small time increments efficiently using an explicit central-difference time integration rule [40]. The top horizontal boundary was free to displace in all directions, while the bottom horizontal boundary was restrained against both vertical and horizontal displacements. All the side boundaries were allowed for vertical displacements but restrained against horizontal displacements. Gravity loading was applied to the rock and tunnel, and the projectile was given a velocity of 1000 m/s in the predefined fields of the load module. General contact was defined between the projectile and target in the interaction section. Hard contact was employed in the interaction module to define the interaction between the rock/steel interface during the assembly of the different parts in the model. This allowed the transfer of stress between the two surfaces in contact using a kinematic contact algorithm that conserved the momentum between two contacting bodies. Separation was allowed after contact. The assembly at 1 m depth with boundary conditions and loading is shown in Figure 5a.

2.4. Mesh

Adaptive meshing (ALE) was used to define the mesh region of the parts in the step module. This meshing technique is a combination of the features of pure Lagrangian analysis and pure Eulerian analysis that aids in maintaining a superior grade mesh during the entire analysis by allowing the mesh to alter irrespective of the material, even in the case of excessive strain/deformation or material depletion. Arbitrary mesh sizes were assigned to the different parts in the mesh module, and a task was created in the optimization module. Thereafter, the model was analyzed in the job module to observe the results. After the mesh convergence study, C3D8R type elements of size 0.05 for rock, R3D3 type elements of size 0.025 for steel lining, and R3D3 type elements of size 0.01 for the projectile were defined in the mesh module for all the analyses. The meshed assembly at 1 m depth is shown in Figure 5b.

3. Results and Discussion

3.1. Damage Modes of Rock and Tunnel Lining

3.1.1. Tunnel Crown Located at 1 m Depth

The deformation, as well as the stress distribution after 10 ms of the impact on the overburden rock and steel tunnel lining with crown at 1 m depth, are shown in Figure 6. A crater region was visible on the rock after impact, indicating perforation of the rock, as shown in Figure 6a. The projectile was observed to be coming out of the lining after perforating the latter at the crown, as clearly seen in Figure 6b. The steel lining was damaged significantly with a hole in the crown. The tunnel lining was distorted downwards along with the projectile indicating that the stress generated in the tunnel lining was due to the momentum of the projectile.

3.1.2. Tunnel Crown Located at 2 m Depth

The deformation of the rock and steel tunnel lining with a crown at 2 m depth, along with the respective stress distribution 10 ms after the impact, are illustrated in Figure 7. A crater region with circular stress distribution after 10 ms of the strike was displayed on the rock part after impact, as shown in Figure 7a. The projectile could not perforate the steel lining but damaged it extensively, as seen in Figure 7b. Due to the propagation of shock waves, the tunnel lining was subjected to a high stress, distributed over a combination of two elliptical areas.

3.1.3. Tunnel Crown Located at 3 m Depth

Figure 8 displays the deformation with a stress distribution 10 ms after striking the rock and steel tunnel with 3 m crown depth. The impact created a crater with circular stress distribution on the rock part as shown in Figure 8a. The projectile was not able to perforate the lining at the crown after striking. Only a marginal deformation without any significant damage was observed at the crown from Figure 8b due to the propagation of shock waves. The stress distribution pattern on the tunnel lining was observed to be a combination of two elliptical areas.

3.1.4. Tunnel Crown Located at 4 m Depth

Figure 9 displays the deformation as well as the stress distribution of the rock and steel tunnel lining at 4 m crown depth. A crater with a circular stress distribution pattern was created on the overburden rock by the impact of the projectile, as visible in Figure 9a. The steel lining was not damaged as the projectile could not penetrate the crown, as seen in Figure 9b, although it was subjected to a circular stress distribution due to the formation and propagation of shock waves.

3.1.5. Tunnel Crown Located at 4.5 m Depth

The deformation, along with the stress distribution on the rock and steel tunnel with the crown at 4.5 m depth 10 ms after the impact, are shown in Figure 10. The impact resulted in a crater region with circular stress distribution on the overburden rock, as depicted in Figure 10a. As the projectile could not penetrate the lining, no permanent deformation was observed from Figure 10b at the crown with minor stress distributed over the whole lining due to shock wave propagation.

3.1.6. Tunnel Crown Located at 5 m Depth

The deformation of the rock and steel lining at 5 m crown depth, along with the respective stress distribution, is displayed in Figure 11. A crater region with circular stress distribution similar to the previous cases was visible on the rock part after the impact, as seen from Figure 11a. The projectile could not penetrate the lining. No permanent deformation was observed from Figure 11b at the crown, and a scanty of stress distribution regions due to shock wave propagation were observed scattered over the tunnel lining.

3.2. Stresses on the Tunnel Lining

The variations of von Mises stress/impact load on the tunnel crown with time (stress–time history curve) for different crown depths are represented in Figure 12. The impact load or von Mises stress generated on the tunnel lining 5 ms after the impact at the rock surface for all the tunnel depths due to the propagation of shock waves. The linings at 1 m, 2 m, and 3 m crown depths were subjected to a high stress/impact load, and the stress-time history plot showed typical impact load behavior with propagation of the projectile through the tunnel lining. The stress–time history plots at 1 m and 2 m crown depths displayed two von Mises stress peaks with maximum values of 800.7 MPa and 525.1 MPa at 5 ms and 4 ms after the impact, respectively. The maximum von Mises stress values at 3 m and 4 m crown depths dropped to 486.8 MPa and 220.7 MPa occurring 5.5 ms and 10 ms after the impact, respectively. The impact load/von Mises stress at 4.5 m and 5 m depths were measured as 67.7 MPa and 45.1 MPa at 9 ms and 9.5 ms after the impact, respectively; however, the projectiles were not able to penetrate the tunnel crown. This indicates the propagation of shock waves after the impact, which were not able to influence the tunnel lining at such depth significantly and were non-destructive to the structure.

The maximum von Mises stress/impact load generated on the tunnel lining due to the projectile impact is plotted in Figure 13. It was observed that the maximum stress on the tunnel lining started to reduce in a logarithmic trend with an increase in the tunnel crown depth. The regression coefficient (coefficient of determination, R²) for the logarithmic equation of the trend line was calculated as 0.922, which is higher than the acceptable range (>0.7). The impact load became asymptotic with a value less than the yield stress of the tunnel lining material after a tunnel depth of 4.5 m. Therefore, it could be inferred that to minimize the effect of the ballistic load due to a 33.655 kg projectile at a velocity of 1000 m/s, a 3 m diameter tunnel with 50 mm thick steel lining should be located at a crown depth below 4.5 m.

3.3. Displacement of the Projectile Nose

The displacements of the projectile nose for different tunnel depths are plotted in Figure 14. At 1 m tunnel depth, the projectile nose displacement reached a peak of 1.250 m at 2.2 ms, after which the rate of change of displacement became negligible, indicating that the projectile was attached to the lining after perforation. The maximum depth of penetration in this condition was 1.268 m at 10 ms. At 2 m tunnel crown depth, the projectile nose was observed to be displaced to a maximum depth of 1.627 m after 4 ms, which was greater than that of 1 m crown depth. In this case, the rate of change of displacement decreased after the peak, indicating a reversal of projectile direction. At 3 m tunnel depth, the projectile was observed to be displaced to a maximum of 2.066 m after 6.5 ms, showing similar behavior as the crown depth at 2 m. At 4 mand 4.5 m tunnel crown depths, the maximum displacements were 2.698 m, and 2.168 m, respectively. At these crown depths, the projectile nose displacements continued to increase, indicating the advancement of the projectile through the rock mass. At 5 m tunnel depth, the projectile nose was observed to be displaced to a maximum of 1.837 m after 6.5 ms and followed by a reversal of the direction of displacement.

The maximum displacements of the projectile for the tunnel crown at different depths are presented in Figure 15. The maximum depth of penetration showed an increasing trend with increasing crown depth up to 4 m, after which it decreased. At a very shallow crown depth (<4 m), the strength of the steel lining dominated in penetration resistance despite the dynamic strength augmentation in the soft sandstone. The maximum projectile penetration depth was observed as 2.698 m when the tunnel was at a depth of 4 m. The dynamic strength augmentation of soft sandstone was employed to resist the projectile penetration after the depth of 4 m, thus reducing the projectile penetration depth for the crown depth after 4 m.

3.4. Velocity of the Projectile

The variations of the velocity of the projectile with time for different tunnel depths are plotted in Figure 16. At 1 m, 4 m, and 5 m tunnel depths, the velocity–time history curves of the projectile were found to be irregular or wavy in nature. At 2 m and 3 m crown depths, the velocity of the projectile became negative at 4 ms and 7 ms indicating a reversal of direction which was also inferred from the projectile nose displacement. The resistance from steel was contributing to this reversal of direction of displacement. At 4 m and 4.5 m crown depths, the velocity of the projectile continued to decrease but did not become negative. As a result, the displacement direction was not reversed. The resistance from rock and its strength augmentation under the dynamic load were contributing to this. At a 5 m crown depth, the rock reduced the velocity of the projectile to a negative value at 6.5 ms, due to which the direction of displacement reversed.

4. Conclusions

In this study, the response of a 3 m diameter shallow tunnel in soft sandstone rock subjected to impact load by a projectile was observed with the help of numerical analysis in ABAQUS/CAE 6.11. The tunnel crowns were placed at different depths from 1 m to 5 m. The Drucker–Prager damage criteria and the Johnson–Cook damage model were used to define the properties of the rock mass and steel tunnel lining, respectively. The non-deformable projectile strikes normal to the rock surface with a velocity of 1000 m/s.

The projectile caused noticeable deformation or damage to the tunnel lining up to a crown depth of 3 m. Moreover, the shape of the stress distribution pattern due to the impact load at the crown appeared unique depending upon the tunnel crown depth from the overburden rock and propagation of shock waves.

At depths of 1 m, 4 m, and 5 m, the velocity–time history curves of the projectile were irregular or wavy in nature. The reversal of the projectile direction was observed at crown depths of 2 m, 3 m, and 5 m. Increasing crown depth had a positive effect on the maximum depth of the projectile penetration up to 4 m tunnel crown depth, after which the effect reversed, making the tunnel safer. The projectile penetrated to a maximum depth of 2.698 m when the tunnel crown was at a depth of 4 m.

The maximum stress at the tunnel lining reduced in a logarithmic trend with an increase in the crown depth, gradually lowering to an impact load less than the yield stress of the tunnel lining material after a crown depth of 4 m. At a very shallow crown depth (< 4 m), the strength of the steel lining dominated the penetration resistance despite the dynamic strength augmentation in the soft sandstone. The influence of dynamic strength augmentation of soft sandstone instigated after a crown depth of 4 m, hindering the projectile penetration into the rock.

From this study, to mitigate the impact of a ballistic load caused by a 33.655 kg projectile traveling at a velocity of 1000 m/sec, it is recommended that a 3 m diameter tunnel, with a 50 mm thick steel lining, should be positioned at or below 4.5 m crown depth.