Dynamic failure analysis of underground reinforced concrete dome bunkers subjected to external blast loading

Abstract

This study investigates the dynamic failure mechanisms of underground reinforced concrete (RC) dome bunkers subjected to both surface and subsurface blast loading, addressing a gap in literature regarding 3D monolithic geometries. A high-fidelity numerical model was developed using the Coupled Eulerian-Lagrangian (CEL) method in Abaqus/Explicit, incorporating the Johnson-Holmquist II (JH-2) constitutive model for concrete and the Mohr-Coulomb model for soil. The framework was validated against numerical and experimental data, achieving less than 5% deviation in near-field peak pressures. Parametric investigations assessed the influence of burial depth, charge orientation, and shell thickness. Results reveal a critical “geometric vulnerability” to lateral blast vectors, where horizontal detonations induced 75% more concrete volume loss than equivalent overhead blasts due to the bypass of the dome’s compressive arching action. Furthermore, a counter-intuitive “Stiffness-Damage Paradox” was identified: increasing shell thickness from 0.5 m to 1.25 m exacerbated damage under specific impulsive loads. Analysis of internal energy histories indicates that rigid, thicker shells trap elastic strain energy, leading to brittle comminution, whereas compliant shells facilitate soil-structure interaction and plastic dissipation. These findings suggest that ductility and geometric efficiency, rather than pure mass, govern survivability in underground protective design. However, readers are cautioned that the manifestation of this stiffness-damage paradox is subjected to the specific structural compliance, soil acoustic impedance, and explosive conditions investigated herein, highlighting the critical need for coupled SSI analysis in protective design.

Keywords

underground RC dome blast loading soil-structure interaction (SSI)Coupled Eulerian-Lagrangian (CEL)failure analysis stiffness-damage paradox geometric vulnerability

Introduction

The escalating prevalence of high-yield explosive threats has necessitated a paradigm shift in the protective design of critical infrastructure. Underground reinforced concrete (RC) facilities are increasingly prioritized for command centers and ammunition storage due to the natural attenuation provided by the soil overburden (Mussa et al., 2017; Zhou et al., 2021). However, the survivability of these buried structures is not governed by structural hardening alone; it is fundamentally dictated by the complex, non-linear Soil-Structure Interaction (SSI) during the blast event (Weidlinger and Hinman, 1988). The propagation of shock waves through the geomaterial creates a transient loading environment where the acoustic impedance mismatch between the soil and the structure governs the magnitude of energy coupling (Mandal et al., 2022; Wang et al., 2005). While deep burial provides protection against airblasts, modern earth-penetrating weapons (EPWs) and accidental subsurface explosions pose severe threats that can bypass the protective soil layer, subjecting the lining to extreme localized pressures (Fan et al., 2023; Yu et al., 2018).

Historically, the design of hardened facilities relied on semi-empirical methods and Single-Degree-of-Freedom (SDOF) idealizations, such as those outlined in technical manuals like TM 5-855-1 (U.S. Department of the Army, 1986). While computationally efficient, these uncoupled approaches often homogenize the geologic medium and fail to capture the spatial and temporal evolution of damage under close-in detonations (Nagy et al., 2010). Consequently, contemporary research has transitioned toward high-fidelity numerical simulations capable of capturing large deformations and fluid-structure interaction. Advanced hydrocode formulations, particularly the Coupled Eulerian-Lagrangian (CEL) approach, have proven superior to traditional Lagrangian Finite Element Methods (FEM) in handling the mesh distortion associated with the extreme discontinuities of an explosion (Keshavarz Mirza Mohammadi et al., 2023; Yang et al., 2018). Recent studies utilizing these methods have successfully characterized failure modes for rectilinear tunnels and cylindrical silos (Huang et al., 2024; Mobaraki and Vaghefi, 2015).

Recent advancements in protective structures have extensively investigated the dynamic response of underground facilities subjected to blast loading. However, the prevailing literature remains heavily skewed towards linearly extended structures, such as tunnels analyzed under plane-strain conditions (Zaid et al., 2022), or vertical cylindrical shafts and silos (Huang et al., 2024). Unlike these cylindrical or rectilinear box structures, dome geometries rely on complex, three-dimensional multi-axial shell action and compressive membrane forces to resist severe external blast pressures. While arches have been studied for symmetric loading (Chen et al., 2013a, 2013b; Chen et al., 2014), the dynamic response of 3D hemispherical domes to asymmetric, close-proximity buried blast vectors remains underexplored. Existing literature rarely quantifies how the geometric efficiency of a dome might become a liability under lateral blast loading.

Furthermore, while the general attenuating role of soil overburden is widely acknowledged, the critical influence of dynamic Soil-Structure Interaction (SSI) driven by acoustic impedance mismatch (Mandal et al., 2022; Wang et al., 2005; Weidlinger and Hinman, 1988) remains inadequately quantified for three-dimensional shell structures (Chen et al., 2013a). Previous parametric studies (Gui and Chien, 2006) often treat structural survivability as a direct monotonic function of the surrounding soil’s stiffness and the structure’s thickness, neglecting the stiffness-damage paradox (Krauthammer and Astarlioglu, 2017). Standard static design principles suggest that increasing thickness enhances resistance; however, under impulsive shock loading, the complex energy coupling between diverse geomaterials and the rigid concrete shell introduces a counter-intuitive phenomenon—a “stiffness-damage paradox” (Krauthammer and Astarlioglu, 2017). In this paradox, stiffer confining media or increased structural rigidity may inadvertently attract higher stresses, amplifying the transmitted shock wave and exacerbating brittle structural damage rather than mitigating it. Currently, there is a critical research gap in systematically quantifying how varying geomechanical properties and threat vectors dictate the transition from structural stability to catastrophic failure in buried dome bunkers.

To address these deficiencies, this study presents a rigorous numerical investigation into the failure mechanisms of underground RC dome bunkers subjected to external blast loading. A 3D non-linear dynamic finite element model is developed using the CEL approach in Abaqus/Explicit. The Johnson-Holmquist II (JH-2) constitutive model is employed to simulate the brittle behavior of concrete under high pressures, while the surrounding soil is modeled using the Mohr-Coulomb criterion to capture shear failure. This research systematically examines the influence of charge weight, burial depth, and threat orientation on structural integrity. Crucially, this study challenges conventional design wisdom by quantifying the inverse correlation between dome thickness and blast resistance for specific threat vectors, identifying critical thresholds where increased rigidity exacerbates concrete comminution rather than mitigating it. While this counterintuitive ‘stiffness-damage paradox’ qualitatively highlights the dangers of energy trapping in rigid structures, it must be emphasized that its exact occurrence and magnitude are heavily subjected to the specific combinations of soil acoustic impedance, structural compliance, and explosive conditions investigated herein.

Constitutive models of materials

The fidelity of the Coupled Eulerian-Lagrangian (CEL) simulation relies on the accurate characterization of material behavior under high-strain-rate loading. Five distinct material phases were modeled: the concrete dome, steel reinforcement, surrounding soil, TNT explosive, and air.

Concrete: Johnson-Holmquist II (JH-2) model

To capture the brittle failure mechanisms of the dome, specifically tensile spalling and comminution under blast waves, the Johnson-Holmquist II (JH-2) constitutive model was employed (Johnson and Holmquist, 1994). Unlike standard plasticity models, JH-2 accounts for the progressive accumulation of damage (D), transitioning the material from an intact state to a fully fractured state based on a polynomial equation of state (EOS). The model incorporates strain-rate effects and pressure-dependent strength, which are critical for predicting the breaching of the dome crown described in this study. The material parameters were calibrated for normal strength concrete (30 MPa) based on triaxial compression data validation by (Xu and Zuo, 2021) and are summarized in Table 1. To simulate the physical disintegration and spalling of the concrete, an element erosion algorithm was implemented within the JH-2 framework. A strict erosion criterion was applied where elements are automatically deleted (removed from the computational domain) when the scalar damage parameter (D) reaches a value of 1.0. This distinguishes between material that has failed but remains in place (cracked/pulverized) and material that has been completely ejected or crushed.

Table 1.

Johnson-Holmquist Concrete (JH-2) Model Parameters of normal concrete. Uniaxial Compressive Strength: 30 MPa (Values, including the dynamic HEL magnitude, are directly adopted from the rigorously calibrated parameters published by Xu and Zuo, 2021).

Parameter	Symbol	Value	Description
Density (kg/m³)	$ρ_{0}$	2400	Initial material density
Shear modulus (GPa)	$G$	12.5	Elastic shear stiffness
Intact strength constant	$A$	0.6304	Normalized intact strength coefficient
Intact strength exponent	$N$	0.8437	Pressure-dependent intact strength exponent
Fractured strength constant	$B$	0.2101	Normalized fractured strength coefficient
Fractured strength exponent	$M$	0.8437	Pressure-dependent fractured strength exponent
Strain rate coefficient	$C$	0.006	Strain rate sensitivity factor
Reference strain rate	${\dot{ε}}_{0}$	1.0	Quasi-static strain rate reference
Tensile strength (MPa)	$T$	3.4	Hydrodynamic tensile failure limit
Compressive strength (MPa)	$f c ’$	30	Unconfined compressive strength
Hugoniot elastic limit (GPa)	$H E L$	1.481	Onset of plastic deformation under shock
Pressure at HEL (GPa)	$P H E L$	0.811	Hydrostatic pressure at HEL
Damage constant 1	$D_{1}$	0.6	Plastic strain contribution to damage
Damage constant 2	$D_{2}$	0.1	Pressure contribution to damage
Bulk modulus term 1 (GPa)	$K_{1}$	16.667	Linear EOS coefficient (low pressure)
Bulk modulus term 2 (GPa)	$K_{2}$	73.19	Quadratic EOS coefficient
Bulk modulus term 3 (GPa)	$K_{3}$	−236.2	Cubic EOS coefficient (high pressure)

Soil: Mohr-Coulomb plasticity

The geological medium surrounding the bunker was simulated using the Mohr-Coulomb elastoplastic model (Seyedan and Seyedi Hosseininia, 2015). explicitly state that simpler models like Mohr-Coulomb are preferable when advanced model parameters cannot be accurately measured. The soil domain parameters, adopted from (Huang et al., 2022) to represent a dense, compacted medium capable of high shock coupling, are listed in Table2.

Table 2.

Geotechnical properties of modeled soils (Huang et al., 2022).

Elastic			Mohr-coulomb
Density (kg/m³)	Modulus of elasticity (MPa)	Poisson’s ratio	Angle of friction (°)	Angle of dilation (°)	Cohesion (kPa)	Abs plastic strain
2200	50	0.3	24	0.0	100	0.0

Steel reinforcement: Johnson-Cook model

The ductility and post-yield behavior of the rebar network were simulated using the Johnson-Cook (J-C) plasticity and damage model (Johnson and Cook, 1985). This model is essential for capturing the rate-dependent hardening and thermal softening of steel during the rapid transient loading of the blast. The rebars were modeled with a yield strength of 792 MPa and were embedded within the concrete host elements. The specific constants used for the reinforcement are presented in Table 3.

Table 3.

Johnson-cook model parameters for steel (plasticity and damage) (Johnson and Cook, 1985).

Elastic properties
Parameter	Symbol	Value
Young’s modulus (GPa)	$E$	200
Poisson’s ratio	$ν$	0.3
Density (kg/m³)	$ρ$	7850

Johnson-cook plasticity model
Parameter	Symbol	Value
Initial yield stress (MPa)	$A$	792
Strain hardening constant (MPa)	$B$	510
Strain hardening exponent	$n$	0.26
Strain rate coefficient	$C$	0.014
Thermal softening exponent	$m$	1.03
Melting temperature (K)	$T_{m}$	1793
Reference temperature (K)	$T r$	293.2

Johnson-cook damage model
Parameter	Symbol	Value
Failure constant 1	$D_{1}$	0.05
Failure constant 2	$D_{2}$	3.44
Failure constant 3	$D_{3}$	−2.12
Failure constant 4	$D_{4}$	0.002
Failure constant 5	$D_{5}$	0.61

Explosive and air (Eulerian Domain)

The expansion of the detonation products was governed by the Jones-Wilkins-Lee (JWL) Equation of State (Lee and Tarver, 1980), which accurately describes the pressure-volume-energy relationship of high-energy explosives. The ambient air, acting as the transmission medium for the blast wave in the Eulerian domain, was modeled using the Ideal Gas EOS. The parameters for the TNT charge and air domain are provided in Table 4 and Table 5, respectively.

Table 4.

Core Parameters values of Jones-Wilkins-Lee (JWL) Equation of state for TNT explosive (Luccioni et al., 2009).

Parameter	Symbol	Value	Description
Detonation velocity (m/s)	$C d$	6930	Speed of detonation wave propagation
Initial density (kg/m³)	$ρ_{0}$	1650	Undeformed explosive density
Initial energy (J/m³)	$E_{0}$	6.0 × 10⁹	Detonation energy per unit volume
JWL coefficient 1 (GPa)	$A$	373	High-pressure response coefficient
JWL coefficient 2 (GPa)	$B$	3.74	Low-pressure response coefficient
JWL exponent 1	$R_{1}$	4.15	High-pressure decay exponent
JWL exponent 2	$R_{2}$	0.95	Low-pressure decay exponent
Grüneisen coefficient	$ω$	0.3	Thermodynamic coupling parameter

Table 5.

Parameters values of ideal gas equation of state for air (Luccioni et al., 2009).

Parameter	Symbol	Value	Description
Gas constant $(J / (kg \cdot K))$	$R$	288.2	Specific gas constant $(R = C p - C v)$
Ambient pressure $(Pa)$	$P_{0}$	1.01 × 10⁵	Atmospheric pressure at sea level
Specific heat $(Cv) (J / (kg \cdot K))$	$C p$	717.3	Heat capacity at constant pressure
Dynamic viscosity $(P a \cdot s)$	$μ$	1.85 × 10^-5	Shear resistance
Initial density $(k g / m^{3})$	$ρ_{0}$	1.225	Density at 20°C and 1 atm

Geometry and finite element methodology

Geometric configuration and reinforcement

The baseline conceptual geometry and scale of the hemispherical dome were adapted from the protective design principles outlined by (Merritt and Newmark, 1958) for structures resisting severe blast overpressures. To rigorously isolate the effects of structural stiffness during the dynamic Soil-Structure Interaction (SSI) parametric study, the outer radius of the dome was fixed at 7.5 m. This ensured that as the shell thickness was parametrically varied from 0.5 m to 1.25 m, the external surface area interacting with the Eulerian blast wave remained completely constant. Additionally, the classical configuration was specifically modified for this study by integrating a continuous monolithic raft foundation to enhance global stability under asymmetric lateral loads. As illustrated in Figure 1, this geometry provides a realistic functional footprint typical of underground command centers. The structure was buried at depths of 5 m, 10 m, and 15 m, measured from the ground surface to the dome crown.

Figure 1.

Geometry of the reinforced concrete domed bunker and TNT charge locations.

The reinforcement rebars network was explicitly modeled to capture the confinement effects and tensile capacity during potential breaching. As detailed in Figure 2, the rebar configuration consists of a meridional and circumferential grid for the dome and a standard orthogonal grid for the foundation slab. A nominal reinforcement ratio (ρ) ranging from 0.8% to 1.2% was established at the dome base (springing line) where the circumference is largest. Due to the geometric convergence of the meridional rebar grid, the effective reinforcement ratio increases progressively towards the dome apex. This densification was explicitly captured in the discrete rebar model (T3D2 elements), providing enhanced confinement in the crown region. Because the primary focus of this study is the global dynamic response and macroscopic concrete failure rather than localized bond-slip, the reinforcement was idealized as continuous truss lines without the explicit geometric modeling of lap splices or mechanical joints. The bond between the continuous rebar and the concrete was simulated assuming perfect kinematic strain compatibility using the *EMBEDDED ELEMENT technique within the Lagrangian mesh (Wang et al., 2005).

Figure 2.

3D model depicting concrete domed bunker with the embedded steel reinforcement, (a) anatomy section of the domed bunker and rebar reinforcement, (b) reinforcement rebars of the dome shell and foundation.

Computational domain and discretization

Numerical simulations were performed using the Coupled Eulerian-Lagrangian (CEL) method in Abaqus/Explicit. This formulation mitigates mesh distortion issues common in large-deformation blast problems by allowing Eulerian materials (explosive products and soil) to flow through a fixed mesh while interacting with the Lagrangian structure (Huang et al., 2024; Keshavarz Mirza Mohammadi et al., 2023). The computational domain, shown in Figure 3, spans 60 × 60 × 80 m, comprising a 40 m air layer and a 40 m soil depth. This domain size was selected based on a convergence study to ensure that stress wave reflections at the boundaries did not interfere with the structural response during the analysis window. The model discretization utilized three distinct element types optimized for their respective material behaviors:

⁃ Eulerian Domain (Soil, Air, TNT): Discretized using 8-node linear Eulerian brick elements (EC3D8R) with a characteristic mesh size of 1000 mm.

⁃ Lagrangian Domain (Concrete Dome): Modeled using 8-node reduced-integration brick elements (C3D8R) with hourglass control with an element mesh size 125 mm.

⁃ Reinforcement: Modeled using 2-node linear 3D truss elements (T3D2) with an element size of 200 mm.

Figure 3.

Detailed modeling views of the computational domain: (a) Overall 3D model of the Eulerian domain containing the soil, air, TNT, and structure; (b) Cross-sectional cut through the soil medium illustrating the burial depth and TNT placement; (c) Isolated view of the RC dome bunker and adjacent TNT charge, clearly illustrating the true hemispherical geometry.

The concrete domain was discretized using a 125 mm mesh to accurately resolve localized damage inherent to the JH-2 continuum damage model, while the reinforcement was modeled using a nominal 200 mm mesh to capture 1D axial yielding. Strain compatibility between these differing mesh densities was automatically resolved using the Embedded Region constraint algorithm (Gaur et al., 2019). This formulation interpolates the translational degrees of freedom from the host concrete elements to fully constrain the embedded rebar, ensuring perfect kinematic coupling without requiring coincident nodes.

Mesh sensitivity and convergence analysis

To ensure the numerical results were independent of spatial discretization, a mesh sensitivity analysis was conducted on the Lagrangian domain. This convergence analysis was performed under the severe loading scenario of a 1650 kg surface detonation at a 5 m burial depth (0.75 m dome thickness bunker). Crucially, the full non-linear material constitutive models were activated during this check to ensure the mesh was evaluated against both global elastic deformation and localized plastic damage regimes, specifically tracking the damage accumulation variable (SDV_DAMAGE). Four mesh densities were evaluated with element sizes of 200 mm, 175 mm, 150 mm, and 125 mm. As illustrated in the vertical displacement history at the dome apex (Figure 4(b)), the structural response exhibited significant divergence at coarser meshes (200 mm and 175 mm). However, the deviation in peak displacement minimized significantly between the 150 mm and 125 mm models, indicating that the solution approached an asymptotic value with negligible deviation in the failure topology. Consequently, a uniform mesh size of 125 mm was adopted for the concrete structure. While it is acknowledged that this element size effectively homogenizes high-frequency spalling cracks that occur at the sub-centimeter scale, it provides a robust convergence for the global deformation response and macroscopic structural breaching (volume loss) which are the primary metrics of this study. This resolution offers an optimal balance between computational efficiency and the accurate capture of localized shear failure, aligning with mesh density recommendations used in similar underground blast studies (Koneshwaran et al., 2015; Mandal et al., 2021). Furthermore, to ensure the accurate propagation of the blast shock wave through the geomaterial and to prevent the artificial underestimation of stress transfer, a rigorous supplementary mesh sensitivity analysis was conducted for the Eulerian domain (Figure 4(c) and (d)). Characteristic Eulerian element sizes of 1250 mm, 1000 mm, and 750 mm were evaluated under the maximum threat scenario. The convergence criteria were based on both the peak incident overpressure reaching the soil-structure interface and the resulting global vertical displacement of the dome. The analysis revealed that refining the mesh from 1000 mm to 750 mm yielded only a marginal difference in peak pressure (∼4%) and structural displacement (∼5%), indicating practical convergence. While finer sub-meter Eulerian grids are standard for unconfined free-air detonations, the immense computational cost of applying such discretization across a massive 60 × 60 × 80 m fully coupled 3D SSI domain renders it prohibitive. Therefore, the 1000 mm Eulerian mesh was adopted for all parametric cases, as it provides an optimal and fully converged balance between capturing the high-fidelity physical kinematics of the blast wave and maintaining computational efficiency.

Figure 4.

Mesh sensitivity analysis for the Lagrangian and Eulerian domains: (a) location of the monitored point at the dome apex; (b) vertical displacement history of the Lagrangian domain mesh sensitivity; (c) vertical displacement history of the Eulerian domain mesh sensitivity; and (d) incident pressure history of the Eulerian domain mesh sensitivity.

Boundary conditions and modeling of soil-structure interaction

To simulate a semi-infinite geologic medium, non-reflecting flow boundaries were applied to the lateral and top Eulerian surfaces (Mussa et al., 2017), with the base constrained as a rigid bedrock stratum. Because structural survivability under buried blast loading is heavily governed by acoustic impedance mismatch, the dynamic Soil-Structure Interaction (SSI) was modeled using a Coupled Eulerian-Lagrangian (CEL) formulation. In Abaqus/Explicit, the interaction between the continuous Eulerian domain (highly deformable soil and expanding detonation products) and the Lagrangian RC dome was enforced using a general contact algorithm based on the penalty method (Hibbitt et al., 2014). This algorithm tracks the dome’s exterior surface and computes Eulerian volume fractions, applying restoring contact forces to prevent penetration and robustly transfer blast momentum.

To prevent non-physical numerical noise and hourglass energy instabilities associated with manually over-stiffened interfaces, the explicit solver’s unmodified default contact parameters were utilized. For normal behavior, the default penalty constraint (stiffness scale factor of 1.0) was maintained, allowing the Abaqus solver to autonomously calculate the mathematically optimal penalty stiffness internally based on local element stiffness and dimensions (Hibbitt et al., 2014). Tangential behavior was explicitly defined as frictionless (μ = 0) to provide a conservative, upper-bound estimate of structural damage. By deliberately neglecting beneficial interface shear dissipation, the full blast wave impulse is transferred as normal stress to the concrete dome, consistent with established numerical practices for underground blast vulnerability assessments (Gaur et al., 2019; Tiwari et al., 2020). This approach accurately resolves complex SSI phenomena—such as localized soil yielding and interface separation—while avoiding the mesh entanglement inherent to purely Lagrangian models.

Parametric investigation matrix

A comprehensive parametric study was conducted to evaluate the bunker’s survivability under varying threat scenarios. As detailed in Table 6, the study considered surface and buried explosions with charge weights of 825 kg and 1650 kg. The investigation matrix covers critical variations in charge orientation (vertical, 45∘, and horizontal offsets) to identify the most damaging threat vectors relative to the dome geometry. The selected explosive masses (825 kg and 1650 kg TNT equivalent) were specifically chosen to represent the severe yield range of massive Vehicle-Borne Improvised Explosive Devices (VBIEDs) detonated at the surface, as well as the explosive payloads of large-scale military Earth Penetrating Weapons (EPWs) detonated in the subsurface medium.

Table 6.

Parametric investigation matrix.

Surface explosions
Dome burial depth from the dome crown to ground surface (m)	Dome thickness (m)	Charge weight (kg)	Charge location description (Refer to Figure 1 for TNT locations)	Orientation angle	Number of charges
5	0.75	1650	On the ground surface, directly above the bunker’s apex (Y-axis).	N/A	1
		1650	Two explosion cases located on the ground surface, at horizontal distances of 5.0 m and 10.0 m from Y-axis.		2
		825	On the ground surface, directly above the bunker’s apex (Y-axis).		1
		825	Two explosion cases located on the ground surface, at horizontal distances of 5.0 m and 10.0 m from Y-axis.		2
10	0.75	1650	On the ground surface, directly above the bunker’s apex (Y-axis).	N/A	1
10	0.75	1650	Two explosion cases located on the ground surface, at horizontal distances of 5.0 m and 10.0 m from Y-axis.	N/A	2
15	0.75	1650	On the ground surface, directly above the bunker’s apex (Y-axis).	N/A	1
15	0.75	1650	Two explosion cases located on the ground surface, at horizontal distances of 5.0 m and 10.0 m from Y-axis.	N/A	2

Buried explosions
Dome burial depth from the dome crown to ground surface (m)	Charge weight (kg)	Dome thickness (m)	Charge location description (Refer to Figure 1 for TNT locations)	Orientation angle	Number of charges
5	1650	0.5	One explosion case buried horizontally at 1.5 m to the side of the bunker’s concrete base.	90° (horizontal)	1
		0.75	Four explosion cases buried in a line directly between the ground surface and the bunker’s apex at a distance of 1.5 m, 2.5 m, 3.5 m, 4.5 m above the bunker’s apex.	0° (vertical)	4
			Three explosion cases buried along a 45-degree line at distance of 1.5 m, 2.5 m, 3.5 m from the bunker’s outer surface.	45° (angled)	3
			Three explosion cases buried horizontally at 1.5 m, 2.5 m, 3.5 m to the side of the bunker’s concrete base.	90° (horizontal)	3
		1	One explosion case buried horizontally at 1.5 m to the side of the bunker’s concrete base.	90° (horizontal)	1
		1.25	One explosion case buried horizontally at 1.5 m to the side of the bunker’s concrete base.	90° (horizontal)	1

Numerical validation

To validate the accuracy and reliability of the fully coupled 3D Soil-Structure Interaction (SSI) numerical framework, the present model was benchmarked against the experimental field tests conducted by Tan et al. (2014). Specifically, the validation replicates “Case 2” from their study, which involved the subsurface detonation of a 102.3 g cylindrical TNT charge buried at a depth of 350 mm beneath a reinforced concrete slab (2000 × 2000 × 150 mm), as illustrated in Figure 5. This specific burial depth was selected because it represents a complex interaction zone where the failure mechanism transitions from direct air-blast cratering to confined shock-structure interaction. It is imperative to note that the primary validation of the current CEL methodology is evaluated directly against the physical experimental measurements reported in their field tests, rather than their numerical approximations. The numerical data provided by Tan et al. (2014) were utilized strictly as a supplementary verification for free-field soil stress propagation. To accurately capture the complex fluid-solid coupling and the steep pressure gradients inherent to this small-scale physical experiment, a finely resolved Eulerian mesh with a characteristic element size of 25 mm was adopted for the validation model. This high-resolution discretization ensured that the simulated blast stress propagation and the subsequent structural kinematics accurately reflected the empirical records.

Figure 5.

Schematic graph of the case study of TNT at 350 mm depth.

Verification of blast wave propagation

To verify the soil constitutive model and the CEL coupling algorithm, the free-field pressure histories were analyzed and benchmarked against the reference 2D axisymmetric numerical data provided by (Tan et al., 2014). Figure 6 compares the pressure-time histories at three-gauge locations (P1, P2, and P3) in the soil. As summarized in Table 7, The numerical model demonstrates high accuracy in the near-to-mid field regions. At Gauge P3 (nearest to the charge), the predicted peak pressure of 86 MPa shows excellent agreement with the reference value of ∼82 MPa (a deviation of less than 5%). Similarly, the mid-field gauge (P2) accurately captured the pressure decay, predicting 33 MPa against the reference 36 MPa. The arrival times of the shock front for the high-pressure zones (P3 and P2) matched the experimental records with a discrepancy of 1.1% and 17.6% respectively. This precise synchronization confirms that the soil’s elastic stiffness and density parameters were correctly calibrated to match the physical acoustic wave speed of the test medium. At the furthest gauge (P1), the model predicted a peak pressure of 4 MPa (compared to ∼3 MPa reference) and a delayed arrival time (220 µs vs 160 µs). This deviation in the far field is attributed to the energy dissipation characteristics of the linear-elastic Mohr-Coulomb model used. While the model slightly over-predicts the residual pressure and lags in arrival time at large distances, this does not compromise the validity of the structural analysis, as the primary structural damage is governed by the high-intensity impulse captured accurately in the near-field (P3 and P2) regions.

Figure 6.

Pressure-time curves of the stress gauges at the points shown in Figure 5.

Table 7.

Comparison of Peak Overpressures (Present 3D Model vs. Reference 2D Model).

Gauge location	Parameter	Reference (Tan et al., 2D)	Present model (3D CEL)	Agreement/Error
Gauge pressure-P3	Peak pressure	∼82 MPa (Est.)	86 MPa	4.8% (Good)
Gauge pressure-P3	Arrival time	∼182 μs	∼180 μs	∼1.1% (accurate)
Gauge pressure-P2	Peak pressure	∼36 MPa (Est.)	33 MPa	8.3% (Good)
Gauge pressure-P2	Arrival time	∼153 μs	∼180 μs	∼17.6% (mismatch)
Gauge pressure-P1	Peak pressure	∼3 MPa (Est.)	4 MPa	Acceptable
Gauge pressure-P1	Arrival time	∼160 μs	∼220 μs	Delayed (far field)

Structural failure mode validation

It is important to emphasize that this slab validation serves strictly to verify the CEL fluid-structure coupling and the JH-2 material erosion thresholds under confined shock loading, rather than providing a geometric comparison to the dynamic behavior of the dome. The ultimate validation criterion was the model’s ability to predict the distinct failure mode observed in the deep-burial experiment. Consistent with the experimental observations for Case 2, the numerical simulation (Figure 7) predicted the formation of radial tensile cracks initiating from the center of the slab and propagating outward. Crucially, the model correctly predicted the absence of a conoid crater on the slab surface, successfully distinguishing the deep-buried response from the spalling failure observed in shallower detonations. Quantitatively, the simulation predicted a maximum vertical displacement (bulge) of 20 mm, which is lower than the ∼130 mm bulge reported in the experiment. This discrepancy is attributed to the distinct physical phases of the explosion. As noted by Tan et al., the large experimental bulge in Case 2 was driven primarily by the late-stage quasi-static expansion of blast products (gas pressure) after the shock wave had passed. The current CEL formulation is optimized to capture the high-speed shock front and stress wave transmission, which are the governing factors for the spalling and shear failure modes investigated in this study. Consequently, the model provides a valid and conservative prediction of shock-induced damage limits.

Figure 7.

Failure patterns of the concrete slab subjected to the soil explosion with underground explosion depth of 350 mm, (a) failure mode of (Tan et al., 2014) experimental work, (b) failure mode the validation model.

Validation scaling and precautions

Validating the framework against a small-scale experiment (Tan et al., 2014) primarily verifies the CEL contact algorithm and JH-2 constitutive model. Extrapolating these mechanics to full-scale scenarios relies on the Hopkinson-Cranz scaling law (Z = R/W1/3), a standard practice for evaluating protective structures (Huang et al., 2024). However, strict precautions are required because gravity-induced geostatic confining pressure does not scale linearly in 1-g small-scale tests. Because geomaterial strength is highly pressure-dependent, a rigorous geostatic initialization step was executed in the full-scale models prior to detonation. This ensures the soil and RC dome are subjected to true, unscaled in-situ gravitational stresses. Furthermore, while micro-scale material heterogeneities do not scale perfectly, employing the homogenized JH-2 continuum damage model ensures reliable prediction of global failure topologies at the full scale.

While validation against dimensionally matched, full-scale experimental data remains the gold standard, accessing open-source physical field tests of deeply buried, large-scale RC domes subjected to high-yield detonations is practically unfeasible due to military classification and financial constraints. Consequently, this study relies on the established state-of-the-art methodology of validating the underlying complex mechanics at a scaled level before extrapolating to the full-scale domain via Hopkinson-Cranz scaling laws.

Numerical stability

Finally, the numerical stability of the explicit dynamic solution was verified by monitoring the global energy history of the Lagrangian domain (concrete slab) throughout the detonation process, as shown in Figure 8. To ensure the solution is free from non-physical modes (hourglassing), the ratio of artificial strain energy (ALLAE) to total internal energy (ALLIE) was analyzed. The results indicate a peak hourglass energy of 0.429 KJ units against a peak internal energy of 10.03 KJ. The resulting ratio of 4.3% remains consistently below the 5% threshold recommended in blast engineering literature (Abir et al., 2017; Mandal et al., 2022), confirming that the mesh discretization is physically sound and the hourglass control algorithm is effective.

Figure 8.

Energy time history of concrete slab of the validation model.

Results and discussion

Surface explosion scenarios

This section evaluates the dynamic response of the underground dome bunker subjected to surface detonations. Unlike buried explosions, surface blasts involve a complex energy transfer where the hemispherical air-blast wave couples with the ground to induce a compressive stress wave. The analysis indicates that for the specific surface detonation scenarios investigated herein (825 kg and 1650 kg TNT at burial depths of 5 m to 15 m), structural survivability was heavily reliant on the attenuation capacity of the soil overburden. However, it must be emphasized that the direct structural resistance of the RC dome remains a fundamental and equally important line of defense, particularly in scenarios involving higher-yield threats or shallower burial depths where the transmitted blast load remains strong upon reaching the structure. The peak response parameters for these scenarios are summarized in Table 8.

Table 8.

Summary of peak response parameters for surface blast scenarios, comparing peak overpressure and vertical displacement at the dome crown for varying charge weights (825 kg and 1650 kg) and burial depths.

Case ID	Charge weight (kg)	Burial depth (m)	Lateral offset (m)	Peak pressure at crown (MPa)	Peak vertical Disp. (mm)	Load attenuation (%)	Damage mode of the dome
S1	1650	5	0	10.0	−5.10	—	Tensile spalling
S2	1650	5	5	8.7	−2.50	13.9%	Minor spalling
S3	1650	5	10	8.1	−2.40	19.8%	Elastic response
S4	1650	10	0	8.6	1.62	14.9%	Elastic response
S5	1650	10	5	7.9	1.63	21.8%	Elastic response
S6	1650	10	10	6.9	1.69	31.7%	Elastic response
S7	1650	15	0	7.2	1.32	28.7%	Elastic response
S8	1650	15	5	6.5	1.39	35.6%	Elastic response
S9	1650	15	10	6.3	1.41	37.6%	Elastic response
S10	825	5	0	6.4	−3.28	36.6%	Minor spalling
S11	825	5	5	5.9	−2.50	41.6%	Elastic response
S12	825	5	10	5.8	−2.48	42.6%	Elastic response

Blast wave propagation and ground shock transmission

The transmission of explosive energy from the surface to the buried structure is characterized by a transition from high-frequency impulsive loading to quasi-static ground shock. As illustrated in the time-lapse pressure contours in Figure 9, the soil medium functions effectively as a mechanical low-pass filter.

• Initial Coupling (t < 5 ms): The detonation creates a high-intensity, hemispherical air-blast. A portion of this energy couples into the soil, generating a compressive stress wave with a sharp rise time.

• Propagation and Filtering (5 < t < 20 ms): As the wave propagates vertically, its high-frequency components are selectively damped by the soil’s hysteretic behavior and geometric spreading (attenuation). By the time the wave reaches the 10 m and 15 m depths (Figure 9), the impulsive shock has dispersed into a diffuse stress field.

Figure 9.

Attenuation of ground shock with increasing depth from 5m at top to 15m at bottom. Pressure wave propagation contours at 5 ms, 12.5 ms, and 20 ms from left to right respectively for a 1650 kg surface charge, illustrating the dissipation of energy density as the wave front penetrates the soil overburden.

This filtering effect is critical; it alters the load profile applied to the dome from a shock-dominated impulse (capable of causing local shear failure) to a distributed pressure pulse, thereby engaging the dome’s global compressive arching capacity rather than local shear resistance.

Parametric investigation of dynamic structural response

To quantify the protective efficiency of the soil-structure system, the influence of burial depth, lateral offset, and charge weight were isolated.

Effect of burial depth on load attenuation

The relationship between burial depth and structural survivability is highly non-linear, exhibiting a trend of diminishing returns. As shown in Figure 10, increasing the depth from 5 m to 10 m reduces the peak incident pressure at the crown from 10.1 MPa to 8.2 MPa (14% reduction). However, the displacement response in Figure 11 reveals the true efficacy of the overburden. The peak vertical displacement decreases from −5.1 mm at 5 m depth to −0.91 mm at 10 m depth. This suggests that the first 5 meters of soil overburden are critical for mitigating the destructive impulsive phase of the blast. Beyond 10 m depth, the load is sufficiently attenuated that the structure responds elastically, rendering deeper excavation economically inefficient for the specific threat levels investigated here.

Figure 10.

Influence of burial depth on structural response. Pressure time-histories at the dome crown for 5 m, 10 m, and 15 m depths (1650 kg charge at center), demonstrating the non-linear attenuation of shock intensity.

Figure 11.

Influence of burial depth on structural response. Vertical displacement time-histories at the dome crown for 5 m, 10 m, and 15 m depths (1650 kg charge), demonstrating the non-linear attenuation of shock intensity.

Effect of charge lateral offset

Lateral offset significantly reduces the destructive potential of surface blasts by altering the angle of incidence. As the charge is moved from the centerline (0 m offset) to 10 m offset, the stress wave travel path lengthens, and the wavefront strikes the dome obliquely. The simulation results (Figure 12) confirm that the centered charge (0 m offset) represents the worst-case scenario for surface blasts, generating symmetric compressive loading that maximizes vertical displacement. Offset blasts induce lower peak pressures (8.1 MPa at 10 m offset vs 10.1 MPa at center) and reduced vertical deflection.

Figure 12.

Influence of lateral distance on structural response. Pressure time-histories at the dome crown for 5 m depth (1650 kg charge at center, 5 m side, and 10 m side), demonstrating the non-linear attenuation of shock intensity.

Effect of explosive charge weight

The structural response scales proportionally with charge weight, consistent with Hopkinson-Cranz scaling laws. Reducing the charge weight by 50% (from 1650 kg to 825 kg) resulted in a 36% reduction in peak pressure (10.0 MPa to 6.4 MPa) and a corresponding 36% reduction in vertical displacement (Figure 13). This linear scaling indicates that for surface bursts, where soil attenuation dominates, the structural response remains in the linear-elastic or early plastic range, avoiding the severe non-linear failure modes observed in buried blast scenarios.

Figure 13.

Effect of explosive yield on structural loading. Comparison of pressure-time histories at the dome crown for 1650 kg and 825 kg charges at 5 m burial depth.

Structural integrity and damage assessment

The failure mechanism for surface explosions is distinct from buried explosions. While the overall dome structure remained largely intact (Von Mises stresses remained below the compressive limit of 30 MPa), severe localized damage was observed at the foundation. The observed failure topology reveals a predominant global structural response dictated by the filtering characteristics of the soil medium. In these scenarios, the soil overburden acts as a mechanical low-pass filter; it significantly attenuates the high-frequency, high-amplitude peak overpressure of the initial blast wave (e.g., capping the transmitted stress to 10.0 MPa at a 5 m burial depth) while largely preserving the lower-frequency global impulse. Consequently, the loading profile reaching the soil-structure interface is no longer a localized, high-intensity shock, but rather a broadly distributed, longer-duration compressive stress wave. This sustained global impulse effectively engages the entire monolithic dome, pressing the structure downwards into the underlying media. The dome safely transfers this load via compressive membrane action down to the foundation. As a result, the foundation slab is forced downward against the resistance of the underlying soil bed, inducing massive flexural stresses and a predominant global bending response across the base plate.

This global bending mechanism couples with the stress wave propagation through the structure to produce the primary failure mode: tensile spalling at the inner surface of the foundation slab. As illustrated in the damage contours in Figure 14 (Left), the incident compressive wave propagates through the concrete and reflects off the internal free surface (the bunker floor). Upon reflection, the wave inverts into a high-amplitude tensile wave. Since the concrete’s tensile strength is low (T ≈ 3.4 MPa), this reflected wave—exacerbated by the tension induced by the global downward bending of the slab—exceeds the material limit, causing immediate brittle spalling (indicated by the localized SDV_DAMAGE accumulation). Notably, the dome haunches, typically vulnerable to shear, remained elastic (stress <8.9 MPa), as confirmed by the Von Mises stress history in Figure 14 (Right). This confirms that surface blasts primarily threaten internal equipment via spalling debris at the floor rather than through global structural collapse. This behavior contrasts sharply with the buried blast scenarios discussed in the section titled 'Buried explosion scenarios', where fully coupled ground shock and soil confinement drastically alter the failure mechanics.

Figure 14.

Mechanism of tensile spalling. Left: Damage contours (SDV_DAMAGE) showing localized failure at the foundation inner surface for bunker at 5 m depth (TNT at center with 1650 kg). Right: Von Mises stress history showing peak stresses remaining within the elastic range at the haunches, confirming failure is driven by wave reflection rather than structural overload.

Buried explosion scenarios

Unlike surface detonations where a portion of the blast energy vents into the atmosphere, buried explosions result in fully coupled energy transfer (Ambrosini and Luccioni, 2020), generating high-intensity ground shocks. This section evaluates the bunker’s response to fully confined 1650 kg TNT detonations. The numerical analysis reveals that structural survivability is governed by the geometric efficiency of the dome against symmetric loading and the ductility demand imposed by asymmetric shear loading.

Geometric vulnerability: Vertical vs. horizontal load paths

A critical finding of this study is the identification of the dome’s geometric vulnerability to lateral threat vectors. Standard protective design often prioritizes overhead cover; however, the parametric results summarized in Table 9 demonstrate that the Horizontal (Side-Buried) Blast is significantly more destructive than the Vertical (Overhead) Blast. When evaluating these results, it is imperative to distinguish between the visual representation of damage in the contours (e.g., Figures 15 and 17) and the quantitative ‘Percentage Loss’ calculated in Table 9. In the damage contours, elements rendered in red indicate a severely damaged state (D≈1.0) where the concrete has lost cohesive strength but remains physically confined within the soil-structure matrix, allowing it to still transmit compressive shock waves (simulating pulverized concrete). Conversely, the ‘Percentage Loss’ reported in Table 9 is highly sensitive to the numerical erosion limits applied in the explicit solver.

Table 9.

Summarize the analysis of damage results of each case.

Charge position (refer to Figure 1)	Deleted volume (m3)	Percentage loss (%)
Case-4 (vertical, 1.5 m)	35.01	5.05%
Case-5 (vertical, 2.5 m)	4.20	0.61%
Case-6 (vertical, 3.5 m)	0.01	0.001%
Case-7 (vertical, 4.5 m)	0	0.00%
Case-8 (angled, 1.5 m)	51.65	7.45%
Case-9 (angled, 2.5 m)	12.56	1.81%
Case-10 (angled, 3.5 m)	0.06	0.01%
Case-11 (horizontal, 1.5 m)	61.38	8.85%
Case-12 (horizontal, 2.5 m)	29.98	4.32%
Case-13 (horizontal, 3.5 m)	3.51	0.51%

Figure 15.

Geometric Vulnerability Analysis. (a) Vertical blast engages compressive arching capacity. (b) Horizontal blast generates transverse shear at the haunch, causing catastrophic breaching.

To prevent non-physical mass loss, element deletion in this study was strictly governed by the JH-2 damage parameter reaching its absolute limit (D = 1.0). This strict threshold ensures that elements are only removed from the computational domain upon experiencing a complete loss of cohesive strength, representing extensive comminution or spalling. Consequently, this formulation prevents the artificial inflation of structural volume loss that would otherwise occur if those severely cracked—yet physically present and confined—elements were prematurely deleted. It is acknowledged that numerical element deletion serves as a computational proxy for severe structural breaching and comminution, rather than an exact physical measure of ejected debris mass. However, applying the strict (D = 1.0) limit ensures a conservative and standardized metric for comparing the relative severity of failure topologies across varying threat vectors.

• Vertical Blast (Symmetric Compression): As shown in Figure 15(a), when the charge is detonated directly above the apex (Case 4, 1.5 m standoff), the shock wave propagates symmetrically. The structural resistance in this scenario is governed by a dual mechanism. First, the dome geometry effectively converts this vertical load into membrane compressive stresses, utilizing the concrete’s inherent high compressive strength. Second, the geometric convergence of the meridional reinforcement toward the apex results in a localized increase in the effective reinforcement ratio well above the nominal 0.8% design value. This creates a highly confined “composite cap” at the crown that offers superior resistance to punching shear. Consequently, despite high peak reflected pressures (>500 MPa as shown in Figure 16), the total concrete volume loss was limited to 5.05%. This aligns with findings by (Chen et al., 2013b), who noted that arches subjected to uniform or symmetric loads primarily deform in a compression-dominated mode.

• Horizontal Blast (Asymmetric Shear): In contrast, the horizontal detonation at the haunch level (Case 11, 1.5 m standoff) targets the structure’s geometric weak point: the wall-foundation junction. As illustrated in Figure 15(b), this load path bypasses the compressive arching mechanism. This induces high out-of-plane shear and global bending moments, resulting in a severe global breaching and a concrete volume loss of 8.85%. This failure mode is consistent with the “bending-induced tensile failure” observed by Huang et al. (2024) in underground cylindrical structures subjected to lateral blast loading. This geometric sensitivity is further corroborated by the inclined (45°) blast scenario (Case 8), which resulted in a 7.45% concrete volume loss—an intermediate value between the vertical (5.05%) and horizontal (8.85%) extremes (Table 9). This result confirms a direct correlation between the magnitude of the lateral load component and the severity of the structural breach, demonstrating that the dome’s protective efficiency progressively degrades as the threat vector deviates from the axis of symmetry.

Critical standoff and Failure Mode Transition

The structural response exhibits a distinct transition in failure mechanisms governed by the scaled standoff distance (Z).

• Local Compression-Induced Tensile Failure ( $Z < 0.25 m / k g 1 / 3$ ): At close-in ranges (Cases 4, 5, 8, 9, 11, 12), the response is governed by shock hydrodynamics rather than structural dynamics. The peak incident pressure exceeds the dynamic compressive strength of concrete, resulting in immediate local pulverization (comminution) and full-thickness breaching before global deformation occurs (see Figure 17(a)). This aligns with the “punching shear” failure modes described by Wang et al. (2005) for close-in subsurface blasts.

• Bending-Induced Tensile Failure ( $Z > 0.28 {m / k g}^{1 / 3}$ ): As the standoff increases, the soil attenuates the shock wave intensity. At R = 3.5 m (Case 6), the failure mode shifts to tensile spalling on the interior face (scabbing), driven by the reflection of the compressive wave at the concrete-air interface (see Figure 17(b)). This transition from direct shear to flexural response is consistent with the damage evolution observed in underground structures by Krauthammer (2008).

Figure 16.

Pressure time-histories at the dome crown, inclined, right haunches, and left haunches for case-4.

However, it must be emphasized that these specific scaled standoff thresholds ( $Z < 0.25 a n d Z > 0.28 {m / k g}^{1 / 3}$ ) are highly dependent on the structural geometry, the 30 MPa concrete strength, and the specific acoustic impedance of the highly compacted soil medium modeled herein. While these exact quantitative threshold values cannot be universally generalized, the phenomenological trend—specifically, the transition from localized hydrodynamic punching shear to global flexural spalling as the scaled distance increases—represents a robust, generalizable behavioral framework for the assessment of buried concrete structures. (Figure 18).

Figure 17.

Transition of failure modes. (a) Catastrophic shear punching at 2.5 m standoff (Case-5). (b) Localized tensile spalling at 3.5 m standoff (Case-6).

Figure 18.

The Stiffness-Damage Paradox Analysis. (a) Damage topology of the flexible dome (thickness = 0.5 m), exhibiting ductile membrane deformation. (b) Damage topology of the rigid dome (thickness = 1.25 m), exhibiting extensive brittle fracturing and comminution. (c) Internal energy time-histories revealing a significant “energy gap” of trapped elastic strain energy in the 1.25 m dome, contrasting with the high plastic dissipation efficiency (∼90%) of the 0.5 m dome.

The stiffness-damage paradox: Influence of structural thickness

A counter-intuitive inverse relationship was observed between the structural thickness of the dome and its survivability against close-in buried explosions. This phenomenon was evaluated under the horizontal detonation vector (Case 11), which was identified in the section titled 'Geometric vulnerability: Vertical vs. horizontal load paths' as the governing “worst-case” scenario causing the greatest structural damage. As visually demonstrated in Figures 18(a) and 18(b), the flexible 0.5 m dome maintained global integrity despite large deformations, whereas the rigid 1.25 m dome suffered extensive comminution and material loss. Quantitatively, the concrete volume loss increased from 7.52% to 14.62% as the thickness increased from 0.5 m to 1.0 m (Figure 19). This phenomenon, termed here as the “Stiffness-Damage Paradox,” is governed by two fundamental mechanisms of Soil-Structure Interaction (SSI): Stiffness-Induced Load Attraction and Energy Dissipation Efficiency.

Figure 19.

Concrete volume loss increases with dome thickness due to stiffness-induced load attraction.

1. Stiffness-Induced Load Attraction: The response of a buried structure is fundamentally coupled with the surrounding soil medium. As noted by (Chee et al., 2020), a stiffer structure buried in soil tends to attract higher loads, whereas stress is diverted around more compliant structures. This phenomenon is mathematically governed by the decoupling theory established by (Weidlinger and Hinman, 1988; Wong and Weidlinger, 1983;). According to their formulation, the interface pressure (Pi) is a function of the free-field stress (σff) and a radiation damping term (ρcu˙) proportional to the structure’s velocity: $P i = 2 σ f f - ρ c u^{˙}$

• Compliant Response (0.5 m Dome): The thinner dome possesses lower flexural stiffness, allowing it to accelerate and deform rapidly (u˙>0) in phase with the soil medium. This relative motion activates the radiation damping term (ρcu˙), effectively relieving the interface pressure through soil arching (Olarewaju et al., 2012).

• Rigid Response (1.25 m Dome): Conversely, the thicker domes act as rigid inclusions within the soil. Their high stiffness resists acceleration (u˙≈0), preventing the relief term from activating. Consequently, the structure attracts the full magnitude of the reflected blast wave (2σff), resulting in significantly higher stress concentrations at the fluid-structure boundary.

2. Trapped Elastic Energy vs. Plastic Dissipation The internal energy history plots (Figure 18c) reveal a critical divergence in energy management strategies. The flexible 0.5 m dome exhibits a ductile response, where the Plastic Dissipation Energy curve closely tracks the Total Internal Energy, indicating that approximately 90% of the input shock energy is consumed through benign global deformation and membrane action. In contrast, the 1.25 m dome exhibits a massive “energy gap” (∼550 MJ) between Total Energy and Plastic Dissipation. Because the rigid dome cannot deform plastically to dissipate the load, this energy is trapped as Elastic Strain Energy. This trapped energy manifests as high-frequency structural vibrations and shear waves that exceed the material’s brittle fracture limit, leading to the extensive comminution (pulverization) observed in the damage contours. Consequently, design optimization for close-in buried threats should prioritize ductility and compliance over pure rigidity. A thinner shell allows for beneficial soil-structure interaction (load shedding) and energy dissipation, whereas a rigid shell attracts critical dynamic loads that lead to brittle failure.

These findings have significant implications for the protective design of underground facilities. They suggest that merely increasing the thickness of monolithic concrete shells may yield diminishing returns or even detrimental effects against close-in impulsive loading due to the trapping of elastic energy. A more effective design strategy shifts the philosophy from ‘rigid resistance’ to ‘controlled compliance,’ utilizing composite architectures to decouple the structural response from the ground shock. Promising avenues for future implementation include the use of ductile inner liners (e.g., steel or UHPFRC) to contain spalling debris, or the adoption of multi-layered double-shell systems that incorporate sacrificial interlayers to dissipate energy through local crushing before it reaches the critical sanctuary. This approach allows the structure to manage high-intensity impulses through global deformation and damping rather than relying solely on material hardness.

It is crucial to caution readers that this ‘stiffness-damage paradox’ is subjected to the specific structural and explosive conditions modeled in this study. Although qualitatively such a paradox may always exist due to the fundamentals of acoustic impedance mismatch, its quantitative effects and failure thresholds will vary significantly from case to case. The intention here is not to advocate against structural thickness, but rather to highlight the necessity of fully coupled SSI analysis in thick-walled protective design. For practical structural performance assessment, this paradox dictates a paradigm shift in survivability criteria. Traditional assessments that equate monolithic wall thickness with a higher safety factor may be fundamentally flawed for close-in buried detonations. Future vulnerability evaluations of existing underground infrastructure must transcend static capacity checks and incorporate fully coupled SSI metrics—specifically analyzing the ratio of plastic dissipation to trapped elastic strain energy. Consequently, accurate performance assessment requires engineers to quantify not just the geometric reinforcement, but the structure’s capacity for global ductile yielding relative to the surrounding soil’s acoustic impedance.

Conclusions

This research elucidated the non-linear response of underground RC dome bunkers to high-yield explosive threats. The following conclusions are drawn:

(1) Geometric Efficiency and Reinforcement: The dome geometry exhibits superior resistance to vertical loads via membrane compression. This resistance is further enhanced at the apex by the geometric convergence of meridional reinforcement, which creates a confined composite cap capable of withstanding high impulsive pressures. However, the structure is highly vulnerable to lateral shear; horizontal blasts targeted the haunch-foundation interface, resulting in severe global breaching and an 8.85% concrete volume loss, compared to only 5.05% for vertical blasts.

(2) Failure Mode Transition: A distinct transition in failure mechanisms was observed between threat vectors. Surface detonations were governed by wave reflections, causing tensile spalling at the foundation’s inner surface while leaving the dome haunches elastic. In contrast, close-in buried detonations resulted in direct compressive punching shear and global bending, driven by the fully coupled ground shock.

(3) Burial Depth Threshold: A burial depth of 10 m was identified as the critical efficiency limit. Beyond this depth, soil attenuation reduces the peak incident pressure and displacement by over 80% relative to shallow burial (5 m), rendering the structural response largely elastic.

(4) The Stiffness-Damage Paradox: Contrary to conventional hardening principles, increasing dome thickness reduced survivability against close-in buried threats. Thicker shells (1.25 m) acted as rigid inclusions, attracting higher interface pressures and trapping elastic energy, which led to brittle comminution. Conversely, thinner shells (0.5 m) utilized soil-structure decoupling and plastic deformation to dissipate 90% of the shock energy, preventing extensive comminution failure. However, it must be strongly emphasized that this paradox is highly dependent on the specific acoustic impedance of the soil and the structural configurations tested, and its exact manifestation will vary under different explosive conditions.

Footnotes

Acknowledgements

The authors express their sincere gratitude to the Department of Civil Engineering, Faculty of Engineering and Technology, Jamia Millia Islamia, New Delhi, India.

ORCID iD

Mahmood Mohammed Taib Maaroof

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Abir

Arumugam

Dhana Sekaran

, et al. (2017) Numerical simulation of blast wave propagation in layered soil featuring soil-structure interaction. In: Proceedings of the 6th International Conference on Computational Methods in Structural Dynamics and Earthquake Engineering (COMPDYN 2015). Rhodes Island, pp. 4752–4765.

Ambrosini

Luccioni

(2020) Effects of underground explosions on soil and structures. Underground Space 5(4): 324–338. https://doi.org/10.1016/j.undsp.2019.09.002

Chee

Krauthammer

Astarlioglu

(2020) Shallow buried RC structures behavior under airblast in the time and P-I domains. Engineering Structures 203: 109885. https://doi.org/10.1016/j.engstruct.2019.109885

Chen

Jin

Fan

(2013a) Elastic responses of underground circular arches considering dynamic soil-structure interaction: a theoretical analysis. Acta Mechanica Sinica 29(1): 110–122. https://doi.org/10.1007/s10409-013-0012-7

Chen

Xia

Zhou

, et al. (2013b) Dynamic responses of underground arch structures subjected to conventional blast loads: curvature effects. Archives of Civil and Mechanical Engineering 13(3): 322–333. https://doi.org/10.1016/j.acme.2013.04.004

Chen

Zhou

Fan

, et al. (2014) Dynamic responses of buried arch structure subjected to subsurface localized impulsive loading: experimental study. International Journal of Impact Engineering 65: 89–101. https://doi.org/10.1016/j.ijimpeng.2013.11.004

Fan

Chen

Hong

, et al. (2023) Considering explosive charge shape and embedded depth in the design of concrete shelter thickness. Defence Technology 20: 44–57. https://doi.org/10.1016/j.dt.2022.09.015

Gaur

Chakraborty

Matsagar

, et al. (2019) Response of box-type structure and effect of structural parameters on reflected pressure under blast load. International Journal of Protective Structures 10(3): 359–379. https://doi.org/10.1177/2041419619854769

Gui

Chien

(2006) Blast-resistant analysis for a tunnel passing beneath Taipei Shongsan airport—a parametric study. Geotechnical and Geological Engineering 24(2): 227–248. https://doi.org/10.1007/s10706-004-5723-x

10.

Hibbitt

Karlsson

Sorensen

(2014) ABAQUS User-Manual Release 6.14. Dassault Systèmes Simulia Corp.

11.

Huang

Zhang

, et al. (2022) Deformation characteristics and damage assessment of prefabricated frame tunnels after central and external explosions. Sustainability 14(16): 9942. https://doi.org/10.3390/su14169942

12.

Huang

Yue

Zhu

, et al. (2024) Failure analysis of underground concrete silo under near-field soil explosion. Tunnelling and Underground Space Technology 147: 105696. https://doi.org/10.1016/j.tust.2024.105696

13.

Johnson

Cook

(1985) Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures. Engineering Fracture Mechanics 21(1): 31–48. https://doi.org/10.1016/0013-7944(85)90052-9

14.

Johnson

Holmquist

(1994) An improved computational constitutive model for brittle materials. In: High-Pressure Science and Technology: Proceedings of the Joint International Association for Research on High Pressure. Colorado Springs, CO, pp. 981–984.

15.

Keshavarz Mirza Mohammadi

Khalilpour

Sareh

(2023) Simulating the response of buried structures to external blast loads: methods, challenges, and advances. Engineering Reports 5(6): e12607. https://doi.org/10.1002/eng2.12607

16.

Koneshwaran

Thambiratnam

Gallage

(2015) Blast response of segmented bored tunnel using coupled SPH-FE method. Structures 2: 58–71. https://doi.org/10.1016/j.istruc.2015.02.001

17.

Krauthammer

(2008) Modern Protective Structures. CRC Press (Taylor & Francis Group).

18.

Krauthammer

Astarlioglu

(2017) Direct shear resistance models for simulating buried RC roof slabs under airblast-induced ground shock. Engineering Structures 140: 308–316. https://doi.org/10.1016/j.engstruct.2017.02.056

19.

Lee

Tarver

(1980) Phenomenological model of shock initiation in heterogeneous explosives. Physics of Fluids 23(12): 2362–2372. https://doi.org/10.1063/1.862940

20.

Luccioni

Ambrosini

Nurick

, et al. (2009) Craters produced by underground explosions. Computers & Structures 87(21–22): 1366–1373. https://doi.org/10.1016/j.compstruc.2009.06.002

21.

Mandal

Goel

Agarwal

(2021) Surface and buried explosions: an explorative review with recent advances. Archives of Computational Methods in Engineering 28: 4193–4208. https://doi.org/10.1007/s11831-021-09553-2

22.

Mandal

Goel

Agarwal

(2022) Dynamic response of underground tunnel in soft soil under surface and subsurface explosion. Practice Periodical on Structural Design and Construction 27(2): 04021081. https://doi.org/10.1061/(asce)sc.1943-5576.0000663

23.

Merritt

Newmark

(1958) Design of Underground Structures to Resist Nuclear Blast. Vol. 2. Structural Research Series No. 149. University of Illinois at Urbana-Champaign.

24.

Mobaraki

Vaghefi

(2015) Numerical study of the depth and cross-sectional shape of tunnel under surface explosion. Tunnelling and Underground Space Technology 47: 114–122. https://doi.org/10.1016/j.tust.2015.01.003

25.

Mussa

Mutalib

Hamid

, et al. (2017) Assessment of damage to an underground box tunnel by a surface explosion. Tunnelling and Underground Space Technology 66: 64–76. https://doi.org/10.1016/j.tust.2017.04.001

26.

Nagy

Mohamed

Boot

(2010) Nonlinear numerical modelling for the effects of surface explosions on buried reinforced concrete structures. Geomechanics and Engineering 2(1): 1–18. https://doi.org/10.12989/gae.2010.2.1.001

27.

Olarewaju

Kameswara Rao

NSV

Mannan

(2012) Response of underground pipes to blast loads. In: Moustafa

(ed) Earthquake-Resistant Structures - Design, Assessment and Rehabilitation. InTech, pp. 507–524.

28.

Seyedan

Seyedi Hosseininia

(2015) Significance of soil compaction on blast resistant behavior of underground structures: a parametric study. Civil Engineering Infrastructures Journal 48(2): 373–386.

29.

Tan

Zhang

Cho

, et al. (2014) Failure mechanisms of concrete slab–soil double-layer structure subjected to underground explosion. Shock Waves 24(5): 545–551. https://doi.org/10.1007/s00193-014-0518-3

30.

Tiwari

Chakraborty

Matsagar

(2020) Analysis of curved tunnels in soil subjected to internal blast loading. Acta Geotechnica 15(2): 509–528. https://doi.org/10.1007/s11440-018-0694-x

31.

U.S. Department of the Army (1986) Fundamentals of Protective Design for Conventional Weapons. Technical Manual TM 5-855-1. Washington, DC.

32.

Wang

Hao

, et al. (2005) A full coupled numerical analysis approach for buried structures subjected to subsurface blast. Computers & Structures 83(4–5): 339–356. https://doi.org/10.1016/j.compstruc.2004.08.014

33.

Weidlinger

Hinman

(1988) Analysis of underground protective structures. Journal of Structural Engineering 114(7): 1658–1673. https://doi.org/10.1061/(asce)0733-9445(1988)114:7(1658)

34.

Wong

Weidlinger

(1983) Design of underground protective structures. Journal of Structural Engineering 109(8): 1972–1979. https://doi.org/10.1061/(asce)0733-9445(1983)109:8(1972)

35.

Zuo

(2021) Study on the JH-2 model parameters for metro shield cutting reinforced concrete pile. Geotechnical and Geological Engineering 39: 1–12. https://doi.org/10.1007/s10706-021-01830-y

36.

Yang

Wang

, et al. (2018) Numerical modeling of surface explosion effects on shallow-buried box culvert behavior during the water diversion. Thin-Walled Structures 133: 153–168. https://doi.org/10.1016/j.tws.2018.09.039

37.

Chen

Fang

, et al. (2018) Blast mitigation effect of the layered concrete structure with an air gap: a numerical approach. International Journal of Protective Structures 9(4): 432–460. https://doi.org/10.1177/2041419618766951

38.

Zaid

Sadique

Alam

(2022) Blast resistant analysis of rock tunnel using abaqus: effect of weathering. Geotechnical and Geological Engineering 40(2): 809–832. https://doi.org/10.1007/s10706-021-01927-4

39.

Zhou

Liu

, et al. (2021) Evaluation of blast-resistant ability of shallow-buried reinforced concrete urban utility tunnel. Engineering Failure Analysis 119: 105003. https://doi.org/10.1016/j.engfailanal.2020.105003