Assessment of Soil–Structure Interaction Effects on the Beirut Port Silos Due to the 4 August 2020 Explosion: A Coupled Eulerian–Lagrangian Approach

Abstract

Blast loadings have become the subject of research in recent decades due to the threats they pose to the surrounding medium. On 4 August 2020, a huge explosion occurred in the Port of Beirut that led to massive damages in the medium surrounding it. Researchers have conducted studies in order to estimate the equivalent explosive mass as well as the damage extent left on structures; however, the studies considered the soil–structure interaction by simple methods. For that, this paper aims to understand the effect of explosion on the grain silo structure present at the port with an emphasis on the soil–structure interaction effects. The structure consists of a group of silos resting on a raft footing that is supported by group of driven piles. A soil–structure model analysis is performed in order to investigate the soil behavior, the damage extent in piles, and the soil–structure interaction due to the Beirut explosion using the CEL (Coupled Eulerian–Lagrangian) approach that suits events involving large deformation. The analysis is performed using the ABAQUS/Explicit FEM software (version 6.14) taking into account the properties of soil medium, the contact algorithm at the soil–structure interface, and the boundary conditions in order to better simulate the real field conditions and ensure accurate results. The work is primarily validated through site data such as the crater size and silo damage.

1. Introduction

In recent years, explosions have drawn great attention from governments and people due to the increase in their occurrence and the massive threats they pose to lives and properties. Furthermore, the substantial risk that blasts might impose on highly representative momentums that can become assault targets has made it vital to study the mechanical effects of a blast on such structures [1,2,3]. An explosion is primarily defined as an intensively rapid and exothermal reaction that persists within milliseconds. It is associated with a tremendously vigorous energy release, typically with high temperature generation and high-pressure gas release. The hot gases expand rapidly, and in hot temperatures cause the air around to expand as well, forming a thin layer of air propagating spherically. This layer is known as blast shock wave and has a speed greater than that of sound waves through which the pressure becomes discontinuous [4]. Compared to static property, dynamic loads imposed from explosions cause a short time dynamic behavior of materials due to the induction of strain rates in the order of 10⁻¹ to 10³ s⁻¹. Such dynamic behavior is represented through materials’ significant strength and escalated stiffness [5].

Two main categories can be distinguished when it comes to investigating blasts: experimental methods and numerical methods. Experimental methods have been widely used to study explosion characteristics and effects. However, with the advanced development of computers and numerical methods, simulations of explosion have become more accurate and cost-efficient where the development time is reduced and protection is optimized. Several research works have been carried out in order to characterize blasts and validate numerical solutions by experimental procedures, and it has been found that numerical mockups of blast loading in the near-field are of adequate accuracy.

Temsah et al. [6] studied the dynamic response of reinforced concrete beams subjected to blast loading using the single degree of freedom method (SDOF) and later using a finite element method tool such as ABAQUS where the nonlinearity of materials, the stiffness degradation, and strain rate effects were implemented. Experimental tests were also conducted to verify the models. It was concluded that the CONWEP model built-in ABAQUS can successfully model the true explosion wave without reflection phenomenon [7,8].

Previous research and many military services have focused on the effects of blast shocks on structural components and on developing blast-resistant designs [9]. However, less coverage has been applied to the effects of blast loadings on substructures. Foundations serve as the most vital part of a structural system for they represent the base of a structure that transmit the loads from upper structures down to the soil. For that, prior to any restoration or recovery operations, it is important to investigate the effects of explosions on the behavior of foundations. One of the important factors to be considered in the behavior of super- and substructures subjected to blasts is the soil–structure interaction system which represents how the soil affects the motion of structures and vice versa.

Soil–structure interaction (SSI) is associated with the wave propagation in a coupled system that consists of buildings erected on the soil surface. Its beginnings can be traced back to the late 1800s, and it has steadily developed and matured over the decades and through the first half of the twentieth century. The SSI phenomenon then exhibited escalated progress in the second half of the 20th century due to the needs of the nuclear power and offshore industries, the introduction of powerful computers and modeling techniques such as finite elements, and the need for seismic safety improvements [10].

Pile–soil–structure interaction (PSSI) is a combination of soil dynamics, structural dynamics, nonlinear behaviors, computational mechanics, computer technology, and many other disciplines. The nature of pile–soil–structure interaction is three-dimensional, and the soil is a nonlinear and anisotropic medium, further complicating the issue. Since the mid-1980s, this issue has been the focus of research in the field of civil engineering and many experts have done a lot of work on this [10]. Different studies have been performed on the soil–pile interaction under seismic loading. However, referring back to the fact that blast loading is very far from seismic loading, only literature based on blast loading will be involved.

Jayasinghe et al. [11] studied the blast response of a pile base in saturated sand using explicit nonlinear finite element analysis. The results indicated that the upper part of the pile appears to be vulnerable, and the pile response decreases significantly as the distance from the explosion increases. In 2014, Jayasinghe et al. [12] conducted a numerical study to comprehend the response of pile foundations subjected to surface explosion using the Arbitrary Lagrange–Euler method. The study showed that the pile’s behavior is affected by the number of piles, the spacing between piles, and the boundary conditions of the pile’s head [11,12]. Bo Huang [13] used the ANSYS/LS-DYNA to investigate the complex response of the pile–soil–structure interaction (PSSI) system under blasting load. Results showed that the peak-shear stress at the top of the pile is greater than the tip, and the distribution laws are inverted triangles.

Ambrosini and Luccioni [14] analyzed the effect of underground explosion on soils and structures through hydrocode and simple soil models that can simulate such phenomena. The soil ejected and crater formation were some of the main outputs of the simulation that validated the numerical work through acquiring close dimensions from experimental results. Souli and Shahrour [15] analyzed the ALE (Arbitrary Lagrangian–Eulerian) method by first simulating the propagation of shock waves caused by an explosion and further by considering a soil–structure coupling problem.

De [16] modeled the effect of explosion on soil and the pressures produced by it using a fully coupled Euler–Lagrange interaction. The model was used to investigate a 2D axisymmetric crater-forming event, and a 3D surface explosion above an underground tunnel case. The results of numerical analyses matched those of other analyses, field experiments, and centrifuge model tests very closely. Hao et al. [17] performed a computational method to explore the pile–soil system elastic and inelastic responses under blasts in addition to the pile’s failure. The pile–soil system is modeled as a beam–column element on Winkler’s foundation. Results included the moment–axial force interaction curve that represented the pile’s yield value and the pile node plastic hinge rotations, which refer to the pile’s failure. Nagy et al. [18] used a coupled numerical model to simulate the behavior of a reinforced concrete buried structure under the effect of blast loads. The findings showed that with the same conditions, buried explosions result in major effects on the buried structure rather than surface explosions.

On 4 August 2020, a dramatic explosion obliterated part of Beirut’s port causing widespread damages to buildings, warehouses, and grain silos. It caused more than 204 deaths, 6500 injuries, and property damage of about USD 15 billion. The blast was initiated as a fire in warehouse number 12 in Beirut Port and shortly after the fire, a colossal explosion occurred forming a mushroom cloud in the air and a supersonic blast wave radiating through the city. The blast was considered one of the most powerful ammonium nitrate explosions in history and it tore the city into shreds [19]. The port of Beirut, being an important and critical point in the country where 48 connected cylindrical silos exist which are based in the ground by a piled-raft foundation, exhibited serious damage which called the structures’ stability into question. Through extensive structural analysis, the true mass of equivalent TNT was found to be 220 tons which corresponds to 564 tons of ammonium nitrate. This amount represents up to 20.5% of the total ammonium nitrate stored. Moreover, it was proved that the silos present in the port didn’t diffract the wave away from the buildings behind where the total force degenerated by them was found to be less than 0.5% of the total energy [20].

Different researchers have worked on the outgoing effects of the explosion on the structures and the surrounding residential buildings through numerical and non-numerical approaches. However, most ignored the destructive effects of the blast on the soil system and consequently on the foundation. Thus, the study of the dynamic soil–structure interaction system is crucial in order to understand the behavior of the silos, the soil, the pile foundation, and their interactions. Better understanding of soil–structure interaction behavior under blast loading will lead to a better comprehension of the consequences of a blast and the stability of the superstructures.

2. Materials and Methods

2.1. Site Characterization

2.1.1. Silos

The Beirut Port has been reconstructed many times to suit the increasing scale of modern container ships. The need for expansion was critical to ensure that the port could handle the largest container ships. The Beirut port arena consists of grain silos that were constructed in the 1960s–1970s to bear a grain storage capacity of 105,000 tons. It involved a total of 42 cylindrical tubes that were constructed in two stages, each acquiring a diameter of 8.5 m, a wall thickness of 0.17 m, and a height of 48 m. Each cylinder possessed a bearing capacity of 2500 tons of grains. The walls of the cylinders converge in three rows that are axially offset by 2.16 m. The first phase involved constructing three rows of silos, each with eight silos, for a total of 24 silos. In the second phase, three additional rows of silos were added as an extension to the existing three rows, each with six silos, for a total of 18 silos. The second-phase silos were 1.2 m apart from the first-phase silos. Later, in 1997, rehabilitation works were carried out due to concrete carbonization and the capacity of the silos increased to 120,000 tons of grains through the addition of six new silos [21] (Figure 1).

The internal arrangement of the silos is depicted in Figure 2. Stiffeners were used to link the silos in both longitudinal and transversal directions. The Longitudinal Stiffeners (L.S.) were 4.4 m apart and linked the silos in a transversal direction. From the ground level (0 m) to level +3.45 m, LS were provided by 65 × 17 cm columns attached to the silos’ walls. For levels above +3.45 m up to the roof level, L.S. were established using 17 cm thick walls connecting the silo walls, as shown in Figure 3. Moreover, 3.0 m long Transversal Stiffeners (T.S.) were used to link the silos in the longitudinal direction with thickness ranging from 17 cm in the middle zone to 52.5 cm at the edge zones. These stiffeners ran all the way from the ground to the roof slab. Except for the last two silos of each row, all silos were filled with grains based on data acquired from various sources and field visits.

Later, between the years 2000 and 2002, the silos experienced rehabilitation works due to concrete deterioration in the two external long rows of silos. Using scabblers and bush hammers, the internal surfaces of the external rows were sandblasted and roughened. The following figures depict the reinforcement details (Figure 4 and Figure 5).

Referring to the study done by Temsah et al. [20], the average compressive strength of concrete specimen extracted from the remaining silos was found to be 15.8 MPa and the average tensile strength was 3 MPa. Furthermore, the density for each specimen was also obtained from its mass and volume, and the average density was found to be 2350 kg/m³. Consequently, the concrete modulus of elasticity was determined to be 19,467 MPa from the preceding parameters and according to ACI 318-14 equation [22]. Moreover, the steel properties were determined using the uniaxial tensile test according to ASTM A615 requirements [23] and they are summarized in

Table 1. Steel rebar properties.

Rebar Diameter in mm	Average Yield Strength (MPa) of Steel Rebars from Silos Prior to Rehabilitation Works (2000–2002)	Average Yield Strength (MPa) of Steel Rebars from Silos after Rehabilitation Works (2000–2002)
10	407	522
12	505	447
14	455	464

:

2.1.2. Foundation

The arena’s foundation is designed using a piled-raft system, comprising 3287 piles with depths varying from 14 to 16 m. The raft component is constructed from type B reinforced concrete, with depths spanning from 400 mm to 3200 mm. In areas where the raft’s thickness is 400 mm, it has been reinforced with 12 mm rebars, positioned at 125 mm intervals in both axes. Meanwhile, in sections where the thickness extends to 3200 mm, 25 mm rebars, also spaced at 125 mm in both directions, have been employed. The piles, as depicted in Figure 6 by a cross symbol, have a square cross-sectional dimension of 300 mm × 300 mm. Each pile is reinforced with four 18 mm rebars and is supplemented with shear reinforcement using 7 mm rebars spaced at intervals of 180 mm.

2.1.3. Soil

Based on the borehole drawings provided, the soil underneath the foundation consists mainly of 2 m of backfill and the rest 15 m of sand and silt.

Table 2. Summary of soil properties.

Properties	First Layer (0–2 m)	Second Deep Layer (2–15 m)
USCS Classification	SP—Poorly graded sand with gravel	SP—Poorly graded sand
Specific Gravity	2.75	2.55
Wet Unit Weight (kN/m3)	17.3	18.85
Angle of friction Φ (°)	37	36
Dilation Angle ψ (°)	7	6
Degree of Saturation (%)	50.4	100

provides a summary of the soil properties.

2.2. Numerical Analysis Method

In the realm of explosion process numerical simulations, various methods such as CEL (Eulerian–Lagrangian Coupled), ALE (Arbitrary Lagrangian–Eulerian), SPH (Smoothed Particle Hydrodynamics), and others have been widely employed. Each of these approaches possesses its own set of merits and drawbacks. The choice of the CEL method for this paper was deliberate and is underpinned by several key considerations. Firstly, CEL is well regarded for its versatility, as it allows for the efficient modeling of complex fluid–structure interactions and shockwave propagation. Secondly, it excels in capturing the dynamic behavior of materials and the dispersion of blast-related effects, making it particularly suitable for the study of explosion processes. Additionally, the CEL method has demonstrated a robust track record in accurately simulating various explosion scenarios, lending credibility to its application in this research [20,24,25]. While other methods such as ALE and SPH are valuable in their own right, CEL was chosen due to its demonstrated suitability for the specific objectives of this study, where precision in modeling explosion phenomena and their effects is paramount.

2.2.1. Geometric Model

The whole geometric model is classified into five distinct features resembling the soil, the air, the TNT, the raft–pile foundation, and the superstructure. The TNT is placed at a distance of 75 m from the front face of the silos based on a previous study done by Temsah et al. [20] (Figure 7). A Eulerian medium is defined to represent the air, soil, and explosive material. The volume fraction tool available in ABAQUS is used herein to create a scalar discrete field by carrying out a Boolean comparison between the Eulerian part instance (medium) and other reference part instances (soil, air, TNT) that cross the Eulerian instance.

The Lagrangian parts are used to define the rest of the model including the piles, raft, and the superstructure. The eight-node brick element C3D8 is used to model the soil and the grains in the 3D explicit model. Four-node shell element S4R is used to model the silo walls with its reinforcement built in and defined by layers, material of each layer, cross-section as per rebar, spacing, and orientation. The lateral and transversal stiffeners were detailed as shell elements. The raft was also modeled using the shell elements with its reinforcement defined as rebar layers inside the shell. Owing to the substantial number of piles present beneath the raft and in order to optimize computational processing time, these piles were depicted as frame elements with two specific profiles. The first profile was rectangular, symbolizing the concrete section, while the second profile was a box-shaped representation of the reinforcement integrated within the piles, as referenced in source [26]. The dimensions of this box profile were formulated based on the area of the reinforcement within the pile’s perimeter.

In order to select the adequate depth of the Eulerian medium, different sizes of Eulerian instances were tried in order to reach the optimum size of the medium at which stresses become zero in both directions. The dimensions of the Eulerian medium were found to be 450 m in the direction parallel to the x-direction, 300 m in the direction parallel to the z-direction, and 170 m in depth (y-direction).

2.2.2. Materials Model

Concrete Model

The mechanical behavior of concrete is modeled using the Concrete Damaged Plasticity (CDP) model which represents the inelastic behavior of concrete by combining isotropic scalar damaged elasticity with isotropic tensile and compressive plasticity. The CDP model is particularly intended for applications involving monotonic, cyclic, and/or dynamic stress applied to plain or reinforced concrete at low confining pressures [27]. In addition, the effect of strain rate on concrete mechanical characteristics was taken into account using CEB-FIP Model relations [28].

Table 3. Concrete properties defined in ABAQUS.

Factors	Symbol	Concrete of Original Silos	Concrete of Rehabilitated Silos	Concrete of Raft Foundation and Piles
Elastic Modulus (MPa)	E	19,467	22,800	21,538
Poisson’s ratio	υ	0.2	0.2	0.2
Density (kg/m3)	ρ	2400	2400	2400
Compressive strength (MPa)	F′c	15.8	23.7	21
Peak Compressive strain (mm/m)	ac	1.06	1.13	1.1
Tensile Strength (MPa)	ft	3	3	3
Strain rate (s−1)	$\dot{\in }$	0.24–4	0.24–4	0.24–4
Dilation angle (°)	ψ	36	36	36
Eccentricity	ɛ	0.1	0.1	0.1
Bi-axial to Uni-axial strength ratio	fb0/ft0	1.16	1.16	1.16
Second stress invariant ratio	K	0.67	0.67	0.67
Viscosity parameter	μ	0	0	0

provides the properties of the materials used.

Steel Model

Both longitudinal and transverse reinforcement were modeled with Von Mises type elasto–plastic material models with user-defined hardening [29]. The model is based on the concept that the steel material will behave elastically until reaching its yielding point, and then it will acquire a perfectly plastic behavior beyond that point until reaching a failure point.

TNT Model

In explosives modeling, the Jones–Wilkins–Lee (JWL) equation of state (EoS) is often used to describe the explosive properties. The JWL expression mimics the pressure (P) caused by the expansion of the detonation and can be expressed in terms of the initial energy (Emo) per unit mass as follows [30]:

$$\mathrm{P}=\mathrm{A}\left(1-\frac{\mathrm{w}\rho }{R_1{\rho }_0}\right)\exp \left(-R_1\frac{{\rho }_0}{\rho }\right)+\mathrm{B}\left(1-\frac{\mathsf{\omega }\rho }{{R_2\rho }_0}\right)\exp \left(-R_2\frac{{\rho }_0}{\rho }\right)+\frac{\mathrm{w}{\rho }^2}{{\rho }_0}{\mathrm{E}}_{\mathrm{mo}}$$

(1)

where A, B, R₁, and R₂ are material constants that have been discovered through dynamic testing for numerous common explosives, ρ is the explosive density, and ρ₀ is the density of the detonation products.

Because the initial relative density (ρ/ρ₀) in the JWL equation is expected to be one, nonzero starting specific energy E_m0 values must be supplied. The TNT charge in this analysis has been modeled by JWL equation of state with its properties attained from [31] and summarized in

Table 4. JWL parameters for TNT explosive.

Parameter	Value	Units
Mass Density	1630	kg/m3
Detonation Wave Speed	6930	m/s
A	3.738 × 1011	Pa
B	3.747 × 109	Pa
R1	4.15	-
R2	0.9	-
ω	0.35	-
Specific Energy	3.68 × 106	J/kg
Ambient Pressure	101,325	N/m2
Specific Gas Constant	-	J/kg·K
Specific Heat	-	J/kg·K
Viscosity	-	kg/m·s

.

Air Model

The air is modeled as an ideal gas with a linear equation of state as a model. The parameters required for an ideal gas model in Abaqus/Explicit are presented in

Table 5. Ideal gas EOS parameters.

Parameters	Values	Units
Mass Density	1.293	kg/m3
Ambient Pressure	101,325	N/m2
Specific Gas Constant	286.9	J/kg·K
Specific Heat	717	J/kg·K
Viscosity	1.82 × 105	kg/m·s

[31].

Soil Model

In this research, the soil is modeled based on Mohr–Coulomb which is widely used for applications in the geotechnical engineering area to simulate material response under essentially monotonic loading. Since the first layer of soil (2 m) is mainly backfill, it acquires properties close to those of the deep layer, and since its total depth of 2 m is negligible compared to the 70 m of soil provided in the model, the properties of the deep layer will be considered in ABAQUS.

2.2.3. Interaction

The “General Contact” feature in ABAQUS is adopted in this study since it allows the definition of contact between all or various areas of the model. The general contact in ABAQUS/Explicit enables the formation of interactions for all exterior faces, analytical rigid surfaces, shell perimeter edges, edges based on beams and trusses, and Eulerian material boundaries present in the model [32]. Due to the large number of piles modeled as frame elements which complicates the definition of contact, the “All with Self” strategy available in ABAQUS is utilized and the “Excluded Pairs” criteria is enabled where all contact pairs other than the Piles’ skin–Eulerian medium pair are omitted.

The raft–pile and the raft–silo connections are resembled through the tie constraint available in ABAQUS which enables a bond between two separate surfaces so that there is no relative motion among them despite the fact that the meshes formed on the surfaces of the regions may be dissimilar. The tie constraint is also based on the concept of master–slave surface assignment. Since the raft and piles are of the same material, the piles’ surface is considered the master surface due to the fact that it acquires a coarser mesh. In raft–silo connections, the silo is the harder material (higher f′_c and modulus of elasticity); thus, it is considered the master surface.

2.2.4. Boundary Conditions

In dynamic analysis of soil–structure interaction, the boundary conditions on the sides of the soil should be defined in a way that prohibits the reflection of waves in the medium. Since the soil model is generated in this study as a fraction of the Eulerian instance, a large size of the Eulerian medium was considered such that it adequately permits the wave propagation without reflecting it back to the medium. A Eulerian size of 450 × 300 × 170 m was found to be adequate for the analysis model. The non-reflective boundaries were used at the surrounding of the Eulerian instances and the bottom part was considered fixed in all directions.

2.2.5. Mesh Size

As mentioned before, the overall geometric model is divided into different regions representing the soil, the air, the pile–raft foundation, the superstructure, and the explosive material. The size of element mesh was chosen in mesh sensibility analysis to reach the optimal mesh size that leads to high accuracy in the analysis yet a reasonable computational run time. The Eulerian medium attained a mesh size of 3 m to provide adequate accuracy, i.e., to achieve damage levels in silos similar to those in the real case scenario. The Lagrangian mesh was used to model the piles, raft, and the superstructure. The silos’ structural elements adopted a global uniform mesh size of 1.5 m, and the raft foundation attained a 1 m mesh size. Piles were meshed with 1.4 m for 14 m piles and 1.6 m for 16 m piles (

Table 6. Mesh sensitivity analysis for Lagrangian parts.

Raft		Piles (14 m)			Silos
Mesh Size (m)	Maximum Vertical Displacement (cm)	Mesh Size (m)	Maximum Stress (N/m2)	Maximum Lateral Displacement (cm)	Mesh Size (m)	Damage Level
3	−1.75 × 10−2	3 m	1.19 × 106	32	3	Part of the 1st row + Extension. The last silo in all rows is destroyed
1	−1.42	1.4	2.09 × 108	38	1.5	First two rows of silos + last two silos in 3rd row. The third row exhibited damages in its internal face
0.8	−1.425	0.4	2.10 × 108	38.1	1	First two rows of silos + last two silos in 3rd row. The third row exhibited damages in its internal face

).

3. Results and Discussion

3.1. Model Validation

One of the crucial aspects of numerical modeling is to validate the model through real-life simulation. In the case of the Beirut explosion, the dimensions of the crater formed at the seaside of the port were used to validate the model. A crater is a bowl-shaped depression that forms as a result of an explosion and is connected to the ejection and removal of soil material from the ground’s surface. The complex properties of air, soil, and the soil–air interface make the crater formation a complicated mechanism. According to field investigations, the crater formed at the port acquired an elliptic shape with dimensions 106 × 90 m as shown in Figure 8.

The crater generated by the numerical model attained dimensions of 105 × 87 m as shown in Figure 9 and Figure 10. The plot of the equivalent plastic strain (vertical displacement) of soil against the true distance of the paths highlighted in red (Figure 11) was used to indicate the crater’s size. The model attained a difference of 0.9% from the real dimensions of the crater in x-direction and a difference of 3% in z-direction. It is worth mentioning that the crater’s dimensions in this study reached values closer to the real dimensions compared to when fixed boundaries were used at the silos’ bases [20] where a difference of 13% in x-direction and 4% in z-direction between the model and real dimensions were achieved. This also signifies the importance of modeling the substructures and mainly considering the soil–structure interaction in the model where more accurate results are attained.

The damage level of the silos serves as a crucial validator for the accuracy and reliability of the numerical model used in assessing the structural integrity of the port. Site visits to the port have played a pivotal role in this validation process, providing tangible evidence of the extent of damage sustained by the silos in different rows.

First and foremost, the observation that the first and second rows of silos are completely demolished is a significant finding (Figure 12). This suggests that the model’s predictions regarding these specific silos were highly accurate, as it aligns perfectly with the real-world conditions observed during the site visits. Such a precise match between the model’s predictions and actual damage levels underscores the model’s ability to simulate and predict structural outcomes in a complex and dynamic environment such as a port.

Additionally, the third row of silos exhibited a particular pattern of damage, with full destruction of the last two silos and varying degrees of destruction in the others (Figure 13). Again, the numerical model’s ability to replicate this damage pattern confirms its effectiveness in capturing the intricacies of real-world structural behavior (Figure 14). This observation further strengthens confidence in the model’s capabilities, as it successfully mimics the specific damage distribution within the third row of silos.

3.2. Superstructure’s Response

3.2.1. Damage

In order to understand the behavior of the silos along with the extension building, damages in concrete were visualized using the output of the CDP model. The tensile damage in the overall structures is depicted in Figure 12 and Figure 13 in order to precisely show the destruction level happening in the concrete structures. As can be seen in Figure 12, the silos along with the extension building exhibited severe damage in their structure.

According to the results demonstrated in Figure 12 and Figure 13, the first row of silos is totally destroyed and the second row of silos is significantly damaged. This is compatible with the real case that occurred at Beirut Port (Figure 14). Figure 13 illustrates that the third row exhibits partial structural degradation, with complete collapse observed in the last two silos of the row, while the remaining silos in this row sustained partial internal damage, primarily along their surfaces facing the point of detonation. These results were in accordance with the real case scenario and the ones concluded from the 3D scanning of the site [33].

3.2.2. Displacement

The vertical displacements were studied for the rest of the remaining silos in the third row (silos 1 to 14). Figure 15 shows the silo numberings in sideview. The vertical displacement at each of the remaining silos was obtained and plotted as shown in the graph below (Figure 16). According to the plotted graph, it can be deducted that the silos encountered settlement of maximum value −3 cm at silo number 9.

In addition to the vertical displacements determined for the silos, the side sway of each silo was also determined through extracting the lateral displacements in the blast direction (in z-direction) at the top of each silo at different time instances (Figure 17). At time t = 0.25 s, the maximum lateral displacement was localized in silo number 8 with a value of 22.7 cm and then in silos number 9, 10, and 11 with a value of 20 cm. Later, at time t = 0.4 s, the remaining silos of the third row exhibited maximum lateral displacements in the direction of the blast of values ranging between 20 and 30 cm, and the maximum lateral displacement was attained in the first and the eighth silos with values of 30 cm and 27.2 cm, respectively. The location of the maximum lateral displacements at the maximum and end times in silos number 1, 8, and 14 in a sequential manner may refer to the fact that they resemble ends of the connected structures. In other words, when the explosion occurred, silo number 8, which is located at the separation joint between phase 1 and phase 2 and 3 silos, experienced the highest lateral displacement. Later, after the destruction of the extension building and silos number 15 and 16, the first, thirteenth, and fourteenth silos resembled the ends of the overall structure and thus experienced the highest lateral displacements. These results differ from those obtained when the silos’ bases were considered fixed [20] where the maximum displacement occurred in silo number 9 and the lateral displacements of all silos exhibited lower values.

The maximum lateral displacement proves that the silos are not stable since the lateral sway reached a value greater than those acceptable in the codes. Moreover, due to the scarce availability of lateral displacement limits in the case of non-building structures such as silos, a research paper written by Kunpeng Guo [34] is taken as a reference to compare drift ratios. The limit states in the latter are determined through different peak ground accelerations and IDA (incremental dynamic analysis) curves.

Table 7. Definition of limit states for selected RC silo [34].

Structural Damage State	Description of Damage	Drift Ratio
None or slight damage	No or slight damage to the RC silo.	0.001074
Moderate damage	Cracks of RC silo shell can be seen remarkably; the stored materials do not contribute to developments of existing cracks.	0.003265
Severe damage	The stored materials contribute to the development of cracks of the silo shell; a sudden collapse of silo structure may occur.	0.005571
Total damage	Partial or complete collapse of the silo occurs.	0.008000

provides the limit states provided for the RC silo based on drift ratios. According to the limit states provided by [34] or an RC silo, the silo is said to be totally damaged for a drift ratio of 0.008. The maximum drift ratio occurring in the silo of the third row is plotted against the determined limit states in Figure 18. It was found that the drift ratio provided in the silo is 0.03, which indicates the total damage of the silo.

3.3. Substructure’s Response

3.3.1. Behavior of the Raft Foundation

Damage

The average tensile damage under each silo row has been assessed and is presented in Figure 19. The raft supporting the silos exhibits damage at multiple locations across its cross-section. In close proximity to the explosion site, linked to the first row of silos, substantial tensile damage has been recorded, totaling 70%. In contrast, the last row of silos demonstrates damage of approximately 60%. Remarkably, the middle row of silos has incurred a lesser degree of damage in comparison to the outer rows, with an estimated damage level of 40%.

Displacement

Figure 20 provides valuable insights into the vertical displacements of the silo’s raft across three distinct rows. Notably, these displacement patterns can be attributed to the lateral load effect of the blast loading, which exerts varying forces on different rows of silos. In the first row, a consistent uplifting trend is observed, with displacements ranging from 1.11 cm to 1.74 cm. This upward movement is most pronounced from Silo 5 to Silo 16, and it stems from the lateral forces imposed by the explosion, causing the silos in this row to be lifted upwards. Conversely, in the second row, a clear settlement pattern emerges, as evidenced by negative displacements intensifying from −0.93 cm (Silo 5) to −1.45 cm (Silo 16). Similarly, the third row exhibits more substantial settlement, with values ranging from −1.85 cm for Silo 5 to a significant −3.63 cm for Silo 16. These settlement patterns result from the lateral load effect of the blast, which induces a downward settling motion in the second and third rows. Given these disparities and the lateral load influence, particularly at the extreme ends of each row, it is imperative to assess the structural implications associated with these vertical displacements.

3.3.2. Behavior of Piles

In saturated soil environments, when a reinforced concrete pile is exposed to blast loading, it is susceptible to various forms of damage, including bending failure, compressive damage of the concrete, shear failure, and significant settlement. Due to the transient nature of blast loads—possessing high intensity but short duration—shear stresses develop more rapidly than flexural deflections. This implies a higher likelihood of shear failure occurring at pile heads, given their fixity. To better understand the extent of pile damage, stress thresholds have been established based on the concrete pile’s compressive and tensile limiting stresses, referenced from ACI code [22]. Specifically, the compressive strength is denoted as 0.33 f′_c = 6.93 × 10⁶ N/m², while the tensile strength is 0.67 fy = 8 × 10⁷ N/m². Figure 21 illustrates that numerous piles have stress values surpassing these strengths in both compression and tension, indicative of profound damage. Notably, piles beneath the first row of silos exhibited severe tensile damage, aligning with the raft foundation’s observations of uplifting in this region. This uplifting effect has seemingly exerted an upward pull on the underlying piles, causing this tensile distress. Conversely, the piles underneath the third row of silos predominantly faced compressive damage. This is consistent with the observed settlement of the raft in this region, which presumably exerted downward pressure on the piles, leading to pronounced compressive impairment.

Figure 22 shows the average lateral displacement of piles beneath three silo rows following exposure to blast loading. For the piles under the first silo’s row, a displacement of 10.4 cm at the head and 1.20 cm at the tip is evident, suggesting that these piles experienced the highest lateral force from the blast. Piles beneath the second row exhibited a slightly reduced displacement, with values of 8.3 cm at the head and 0.96 cm at the tip, while those under the third row displayed the least displacement, with values of 7.2 cm and 0.84 cm, respectively. The consistent trend of diminishing displacement from head to tip across all rows highlights the attenuating impact of soil depth on the blast’s energy. Moreover, the descending pattern of displacement across the silo rows might imply a differential distribution of blast energy, possibly due to varying distances from the blast epicenter or the influence of intervening structural or geotechnical barriers. This underscores the importance of understanding the nuanced interactions between blast forces and infrastructure components to ensure structural resilience.

3.3.3. Behavior of Soil

Stress and Strain in Soil

Figure 23 and Figure 24 visually depict the evolution of plastic zones in the soil, manifested as plastic strains. Specifically, Figure 23 showcases these strains in the vicinity of the crater while Figure 24 highlights the corresponding observations around the silos. For reference, PE11, PE22, and PE33 represent the plastic strains in the x, y, and z orientations, respectively.

An analysis of Figure 23a–c reveals the greatest plastic strain at the crater’s epicenter, with a pronounced maximum of 1.762 in the x-direction. Conversely, Figure 24a shows the peak plastic strain of 0.0712 in the x-direction, notably at the front and western boundaries of the silos. Furthermore, the y- and z-directions manifest their respective maximum plastic strains at the silo’s frontage, with magnitudes of 0.09 and 0.0257, as indicated in Figure 24b,c. The presence of these nonzero plastic strain measurements suggests that the soil has undergone irreversible deformations.

To further assess the condition of the soil, a Mohr’s circle was constructed using the maximum and minimum principal stresses. As illustrated in Figure 25, the soil beneath the silos approaches the failure envelope but does not surpass it. This suggests that the soil under the silos has not undergone failure but was severely damaged.

Soil–Pile Interface

The data presented in Figure 26 highlights the average contact shear stress between the piles and the surrounding soil at different depths, both at the initial time (t = 0 s) and 0.8 s post-explosion. At the outset (t = 0 s), the shear stress exhibits a mild fluctuation with depth. Starting at 38,900 kPa at the surface, it peaks at 45,000 kPa at a depth of 2 m, before slightly decreasing to 43,900 kPa by 14 m. However, in the aftermath of the explosion, there is a marked reduction in shear stress across all depths, which is indicative of the profound impacts of the explosive event on soil–pile interactions. The first silo row exhibits shear stresses ranging from 6900 kPa (at the surface) to 15,700 kPa at a depth of 14 m. Remarkably, the second row demonstrates more amplified contact stresses, culminating at 21,400 kPa at 14 m depth.

In Figure 27, the relationship between the initial and subsequent contact stress post 0.8 s of the explosion at varying pile depths is distinctly illustrated. For the piles situated beneath the first row, there is a discernible ratio ranging from 0.18 at the pile head to 0.36 at its tip. When considering the piles under the second row of silos, the ratio fluctuates between 0.36 and 0.49. Concurrently, for the third row, the ratio varies from 0.17 to 0.39. A critical analysis of these findings strongly suggests that the reduction in contact between the pile and the soil is most pronounced under the first row of silos, which is in immediate proximity to the explosion source. A further implication from the data underscores that, across all three scenarios, the pile head experiences a more significant loss of soil contact compared to the pile tip.

4. Conclusions

In this research, an advanced nonlinear numerical analysis was conducted to assess the response of silos, raft, piles, and soil to significant blast loading, with the Beirut explosion serving as the primary case study. The study involved comprehensive site investigations, detailed geometry and material parameter determination, and numerical assessments using ABAQUS software (version 6.14) with the finite element method and the CEL approach. Key findings from the study highlight the importance of including the substructure model, emphasizing the significance of soil–structure interaction (SSI), and revealing structural instability indicators.

The study revealed that including the substructure model is crucial for accurate results, as demonstrated by precise crater dimensions and reduced variances compared to models without a foundation. It also emphasized the importance of SSI through observations of vertical displacements in the silos and extreme sway values at their summit. Furthermore, the research highlighted the intrinsic link between a structure’s stability and the integrity of its underlying substructure, showcasing severe damage to the piles and raft–pile foundation.

Additionally, the study identified stress levels exceeding the failure threshold in the immediate vicinity of the explosion, indicating soil degradation. However, soil near the silos remained within the failure envelope. There was also a notable loss of contact between the piles and the surrounding soil.

In light of these findings, the research underscores the urgency of prioritizing post-blast rehabilitation efforts and raises concerns about the suitability of specific damaged components for future use. It calls for interdisciplinary collaboration among geotechnical engineers, structural engineers, and blast experts to develop robust post-blast restoration techniques and strategies. Ultimately, this collaborative approach may lead to the design of structures with inherent resilience to blasts and similar high-impact events, contributing to enhanced structural safety and functionality.