Abstract
This paper presents studies that focus on fire and explosion-induced damage of tunnel structures by employing the Discrete Element Method (DEM). By assuming a two-dimensional aggregate distribution and reconstructing the digital representation of the experimental concrete blocks, a numerical model of the tunnel lining concrete was established in the PFC2D program. The temperature distribution and the shock wave pressure at the surface of the tunnel lining were obtained by using Fluent and LS-Dyna separately; the final damage simulation of concrete section under different conditions was carried out in PFC2D. The results showed that PFC2D cooperatively provided more accurate and effective modeling and visualization of impact damage of concrete blocks. The visualizations of damage indicated the degree of damage more clearly and more intuitively. These findings also provide a potential method for further study of the damage assessment for entire tunnel lining structures.
1. Introduction
Tunnels have special characteristics such as long longitudinal depth, few entrances and exits, and a semi-closure structure. Thus, a major accident often causes fire and explosion and subsequently leads to serious spalling and secondary disasters [1]. Previous researchers usually simulated and analysed the fire and explosion processes separately when studying the accidents in tunnels [2, 3], but in actual accidents, fire and explosion often occur simultaneously, and different degrees of collaborative damage phenomenon can be observed. Especially when vehicles transporting dangerous chemical materials crash in a tunnel, there will be much higher potential threats to the concrete lining structure of the tunnel due to fire and explosion. Therefore, it is of great significance to evaluate the damage levels of tunnel lining structure caused by both fire and explosion.
Concrete is a porous material exhibiting discontinuity, anisotropy, heterogeneity, and nonlinearity at the same time [4]. The constitutive relation based on the continuum hypothesis reflects the macro response of the system through a specific mathematical model. But this model cannot represent the essential characteristics of local discontinuity of the system, which makes it difficult for the continuum mechanics model to simulate the local instability. Especially, it is hard to describe the meso-scale damage and the discontinuous behaviours of concrete.
Considering the nature of discontinuous structure of concrete, the Discrete Element Method (DEM) provides a new idea for studying the damage behaviour of concrete in meso-scale [2]. In this work, we focus on simulating the impact damage analysis of concrete lining structure at different positions under varied temperatures, including deformation and cracking. The aim of this study was to establish the concrete block model with detailed aggregate shape in PFC2D. To calibrate the different contact parameters between the clustered particles, the stress-strain curves of concrete blocks under different temperatures were evaluated in the following experiments.
2. Simulation of the Temperature Field of the Tunnel Lining Structure under Fire Condition
Simulation is the most feasible method to determine the temperature distribution of a tunnel under fire condition [5]. In this study, we established the combustion model of dangerous chemical-transporting vehicles in the computational fluid dynamics software (FLUENT) to obtain the temperature distribution at critical points of the tunnel.
2.1. The Choice of the Combustion Model and the Turbulence Model
The choice of the combustion model has a great influence on the flue gas flow field. Due to the complexity of combustion, chemical reaction, and thermal radiation, the gas turbulence will be coupled with temperature and the shape of the tunnel. For the choice of fire combustion models, previous scholars mainly adopted the PDFEBU, EDC, PPDF, SHS, and VHS models [6]. The VHS model of volume heat source was used in this study.
The most common hazardous substances in transportation are gasoline, kerosene, and other hydrocarbons. In this simulation, a methane air mixture model was selected in the combustion model, and CO2, O2, H2O, and N2 were added in the material phase. After the fluid was set as mix phase, the chemical reaction equation was set, and the full combustion was assumed in this simulation. The main flue gas products were carbon dioxide and water vapor, and the heat transfer effect of thermal radiation was considered. Since a large amount of natural gas leakage and deflagration were simulated, the turbulent visibility was not considered here.
2.2. Detailed Geometry of the Tunnel and the Boundary Conditions
The parameters of the tunnel model were referred to full-scale simulation by previous researchers [7]. The cross section size of the tunnel was approximately established as a semicircular structure with a radius of 7.5 m and a longitudinal length of 50 m. According to the fact that drivers often stop beside the wall after an accident, the whole simulation model established the fire source on the side of the tunnel for combustion (d = 1.5 m from the lining structure), as shown in Figure 1.
Because the tunnel has a closed and narrow tubular structure, the smoke and dust in the air will diffuse far away after combustion. Considering the critical wind speed, air flow, and the impact of real tunnel fire, this study simulated the longitudinal distance of 50 m.
2.3. Heat Release Rate (HRR)
Heat release rate is one of the most important parameters in both experimental and simulative studies. The heat release curve needs to be obtained according to different conditions (vehicle size, tunnel ventilation environment, and obstacles during combustion). In highway tunnels, the mixed traffic flow is composed of cars and various heavy trucks. The fire power of a car is usually between 3 MW and 5 MW. The fire power of large vehicles transporting dangerous goods, especially oil tank trucks, often exceeds 30 MW and can reach a maximum of 100 MW [6]. We obtained the standard RABT heat release curve (GA/T 714-2007) for this simulation as this curve is suitable for the combustion of hydrocarbon materials [7]. The HRR for the RABT curve was huge at the beginning of the accident in Figure 2 as the temperature approached to 1200°C. After the explosion combustion, the oxygen content in the tunnel space dropped down dramatically and the temperature gradually decreased in the following 200 minutes.
2.4. Temperature Variation in Different Sections along the Longitudinal Direction
When the concrete lining structure is exposed to the high temperature generated by the vehicle fire, the temperature field inside the lining segment will change greatly and rise to varying degrees. As a thermal inert material, concrete has a large specific heat capacity (1000 J/kg-K) and small thermal conductivity (1 W/m-K), so the heat transfer in concrete is relatively slow. Concrete has obvious thermal inertia, so the heat energy requires time to penetrate the concrete lining. As the average thickness of the concrete lining structure is 500 mm, we placed the monitoring points at the middle layer to acquire the temperature inside of the concrete lining structure. The longitudinal temperature detection section was set every 5 m, and the overall detection distance was 0–30 m. The reason why we choose 30 m is that the temperature in tunnel remains relatively stable after the first 20 m. Also, we focused on the combination of explosion and temperature damage modeling and the damage dominated by the shock wave only remained 10–15 m long.
It can be seen from Figure 3, in the longitudinal temperature distribution of the tunnel, that the section with the highest temperature shown up 5 m away from the rear of the fire vehicle. Under the initial combustion conditions, the oxygen concentration was sufficient, so there was a deflagration phenomenon, and the temperature rose sharply; when the local oxygen at the fire point in the tunnel was exhausted, the combustion intensity decreased; then, the flame developed along the downstream direction of the airflow, and the corresponding maximum temperature position gradually shifted backward. When the combustion was stable, the temperature gradually increased along the longitudinal direction. The peak and stable value of the temperature at the four points illustrated in Figure 1 are recorded in Tables 1 and 2.
(a)
(b)
(c)
Most of the high-temperature regions were at the left side (points B and C) and the right side remained at the ambient temperature. The temperature increased significantly at the downstream (around 10 –5 m) of the fire scenarios. From Table 2, we could find that the maximum temperature of concrete lining at 1 m height is relatively low.
3. The Properties of the Explosion Shock Wave
The chemical explosion is a high-speed dynamic process that takes place in a very short time, and it usually takes only tens of milliseconds or even a few milliseconds from the beginning to the end [8]. Therefore, it is difficult to observe the complete process by using the traditional experimental method, and it is hard to repeat the test due to the influence of various factors. At present, the numerical simulation of tunnel explosion is an effective way to get the pressure curve.
The shock wave is produced by compressing the surrounding air. Because the deflagration and explosion process is too complex, there is no recognized pattern to follow. It is convenient for researchers to use the TNT equivalent method to predict and assess the damage of explosion [9].
3.1. Material Property and the Equation of State (EOS) in LS-Dyna
The summary of the material parameters assigned to the simulation was presented in this section. In the LS-DYNA simulation, the Jones-Wilkins-Lee (JWL) state equation is usually used to describe the relationship between pressure P and energy V per unit volume of the TNT charge. The description equation of pressure was obtained from the isentropic curve of cylinder test and the initial density change of explosive (equation (1)):
In Table 3, ρ0 is the density of explosive, is the internal energy per unit volume of explosive, D is the velocity of explosive wave, and V is the current relative volume. A, B, R1, R2, and are the parameters of the JWL state equation, which need to be determined in experiments.
In this simulation, the air (Table 4) was simulated as a “MAT-NULL” material and was described by linear polynomial state equation.where E is the internal energy of the material, ， ρ is the current density of the air, ρ0 is the initial density of the air, and C1–C6 are the parameters of the EOS of the air. The specific parameters are shown in Table 4. In order to simplify the calculation, the surrounding rock simulated in this simulation is considered as homogeneous rock-clay material.
The state equations of explosive and air were defined in the multi-material coupling element. For solid elements, the Euler Lagrange coupling method was used to calculate; the Euler grid method was used to describe air element and the Arbitrary Lagrangian Euler (ALE) algorithm was used for element calculation. Lagrange grid and Euler grid were coupled at the interface in the same analysis model, and the boundary field of material in Euler grid was defined through the transfer of mutual force. After the explosion, the shock wave firstly propagated through the fluid medium and then interacted with the lining structure and the carriageway.
3.2. The Explosion Model and Mesh Size
The size of the element will directly affect the accuracy of the computational results. Theoretically, decreasing the grid size can improve the accuracy of numerical results, but at the same time, it can also greatly improve the scale of numerical calculation thus decreasing the computational efficiency. In addition, the affected area is not very large, so it is unnecessary to divide the full-scale model into small grids when considering the influence of grids on the calculation scale.
In most traffic accidents, the explosion source was located at the side of the tunnel carriageway, and the distance between the explosion source and the lining structure was 1.5 m as, shown in Figure 4. For the concrete lining structure around the center explosion area on the left side, the size of the element unit far from the damage area was set to 100 200 100 mm; the mesh was divided more finely, with sizes of 20 20 20 mm, which are more conducive to obtain more accurate results.
The contact explosion was calculated by adding LOADSEGMENT and DEFINESEGMENT to set the equivalent mass of TNT explosive. The TNT equivalent of explosion source was set as 70 kg according to previous research [10]. The tunnel lining structure freely contacted the surrounding rock surface, and the friction and slip were ignored. The damage of lining concrete was expressed by plastic deformation when it was subjected to explosion stress. The element structure of lining concrete was deleted directly because of the failure after exceeding the plastic deformation threshold of the material.
3.3. The Mach Reflection of the Explosive Shock Wave
In real explosion scenes, the shock wave propagates freely in the air in the form of spherical wave. Regular normal reflection will occur when the shock wave front that is perpendicular to the lining structure normally contacts the lining structure; however, with the distance between the explosion center and the normal projection getting farther and farther, the incident angle will increase. When the incident angle is equal to or greater than the critical angle, the reflected wave front and the incident wave front coincide to form the shock wave along the longitudinal direction of the tunnel, and the combined shock wave is called Mach wave.
According to Figure 5, the reflection less than the distance from the explosion center belongs to normal reflection; as it gradually exceeds the regular normal reflection area, the incident angle between the incident wave and the vertical normal increases continuously, and it belongs to oblique reflection. When the incident angle becomes larger and larger, the angle between the incident wave front and the reflected wave front is smaller, and even overlaps to form a single composite Mach wave wave.
After the explosion, the shock wave expanded in a spherical shape in the center of the explosion. The main propagation mode of air was radial compression, and the shock wave first acted on the projection position of the detonation center at the arch waist. With the increasing time, the front of the shock wave became larger and the energy per unit area decreased gradually. The pressure-time curve at different points is shown in Figure 6.
When the air shock wave encountered the vertical rigid wall, the velocity of air particles on the wall instantly decreased to zero, which made the air particles within the stagnation point rapidly accumulate, and the pressure and density also increased rapidly. As shown in Figure 6, when the air shock wave contacted the rigid wall, the equivalent stress of the wall surface at different observation points increased sharply. When the explosion time was from 0.49 to 0.99 ms, the detonation of explosive had basically completed, the air was the regular positive reflection, and the tunnel lining structure facing the detonation part of the explosive formed the crushing area. When the time was from 0.99 to 13.4 ms, the detonation shape changed from the normal reflection to the Mach wave. The Mach Wave was the compressed reflected air wave, and the propagation direction changed to the longitudinal propagation along the tunnel lining structure after the air propagated to the free surface of the tunnel lining structure. From 13.4 ms to 20 ms, the shock wave formed by the detonation had declined and the Mach wave propagating along the tunnel lining was attenuated.
4. The PFC2D Numerical Model and Damage Assessment under High-Temperature Explosion Shock Wave
4.1. Information Extraction of the Concrete Aggregate Structure
The Discrete Element Method (DEM) allows the relative motion of the internal elements and does not need to meet the continuous displacement and deformation conditions. It was especially suitable for the calculation of large deformation and large displacement. In recent years, the development direction of discrete element simulation begins to focus on the transformation from continuous media to discontinuous media. The damage and failure of brittle materials such as concrete with many pores under dynamic load has become the research direction of many scholars [11].
Among the factors that affect the strength of concrete, the shape and distribution of the aggregate are the most important factors [12]. In DEM simulation, the treatment methods of aggregate are divided into two categories: the first is to directly simplify the aggregate into spherical or ellipsoidal particles for stacking, but it cannot simulate the enhancement of concrete strength due to the mutual connection between the aggregates; the second is to construct the polyhedral shape of aggregate, and then fill in spherical or ellipsoidal particles. Concrete consists of cement mortar, coarse and fine aggregate, and their transition zone [13]. The aggregates in concrete present the characteristics of random distribution. In this study, the digital image processing technology was used to model the particle flow of the concrete test blocks.
The standard concrete test blocks were made in the laboratory and the cement mix parameters are shown in Table 5.
After being stored in a constant temperature conditioning chamber for 28 days, the sections with clear aggregate distribution on the side were selected for photography to generate RGB images (Figures 7–11).
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(d)
The weight average method was used for the gray binarization of the RGB images. Binarization was used to distinguish the pixels representing aggregate from those representing cement mortar to give them different particle flow parameters.
The purpose of binarization is to transform the structural information of concrete blocks into the modeling information of meso structure. In binary processing, we need to select the local gray threshold T as the judgment standard. After the images were divided into several independent regions, we need to find different gray T values. If the gray level of a pixel is greater than T, it is white; if it is less than T, it is black. However, although this method can be completed automatically by MATLAB, it still needs to deal with the redundant and unreasonable details manually, such as isolated points and pseudo points.
After adopting the binary process through manual intervention and adjusting the contrast ratios, the concrete images were imported into the AutoCAD software and the sizes of the images were scaled to the size required for PFC simulation (150 mm 150 mm). The reason why this study chose 2D model in PFC is that the deformation of concrete in the longitudinal direction in each section is assumed to be zero. The 2D model also can provide a clear damage assessment of concrete structure under different compressive waves and different temperatures. The closed boundary line was drawn from the beginning to the end along the concrete aggregate edge by using the polyline.
4.2. Formation and Initial Stable State of the Aggregate Particle Clusters
In the PFC2D, the basic computing element was the disk element [14]. In order to distinguish the aggregate from cement mortar, it was necessary to establish different shapes of the particle clusters. In real concrete, it was impossible to have only round particles, so the surface properties should be specified separately according to the actual situations of the aggregates. Particle clusters were independent of each other, but they influenced each other in the process of the mechanical solution. The basic disk particles in the cluster were regarded as the rigid bodies, but the disks were allowed to overlap inside the particle clusters. The irregular arrangement of particles needs to fill the gap with 2D particles in a given area and keep the overall balance of the model. In the PFC2D model, the particles can be divided into regular arrangement and irregular arrangement. Regular arrangement can be used to describe the structural part of simulation, while irregular arrangement can be used to simulate solid and internally irregular granular materials.
Particle clusters in this simulation were divided into mortar and aggregate which are both composed of irregular particles. Although the aggregate particles were randomly arranged, the structural characteristics of the whole model will still be affected by the existence of anisotropy and weak structural planes.
When producing the aggregate model, if the proportion of the aggregate particles in PFC2D is higher than the actual concrete block, the compression strength of the whole concrete model will not match with the real concrete block and raise up greatly. The aggregate characteristics are used and extracted from the real concrete sections to guarantee that the same area ration of aggregates can effectively avoid this problem [15].
There are two ways to fill particles with specific radius to achieve the required porosity. The first is to establish the boundary of the closed area, and then generate a series of noncontact particles in the closed area. Then, the required porosity can be achieved by moving the boundary. The disadvantages of this method are that the geometry of the boundary will eventually change, and that the distribution trend of the final particles will not be uniform. The second method uses the generate command and Gaussian distribution, which specify the upper and lower limits of particle radius (radius of 0.0008 and 0.001) and assign a standard deviation at the same time, and then select and generate particles. The particles were generated after the sample area size (domain extent −0.5 0.5–0.5 0.5), porosity (0.19), and particle distribution were specified.
In equation (3), where n is the number of particles, Vp is the void area, Vt is the total area, Va is the particle area, and the void ratio Rp is defined as the ratio of the void area to the total area. When we specified the radius and size of particles, we further used equation (4) to determine the number of particles:where R1 is the minimum particle radius (lower limit), R2 is the maximum particle radius (upper limit), and N is the final number of particles, is the average radius of the particles.
After the aggregate cluster and the mortar were filled, not all particles were in proper contact, which may result in errors in the calculation of the model. Therefore, it was necessary to define the contact gap as suspended particles, and the radius of these suspended particles can be expanded to meet the contact status between these particles. Finally, the combined sphere with a certain porosity was obtained, and then it was necessary to assign the contact parameters according to the concrete materials. Since the number of particles remained unchanged before and after the amplification, the amplification factor A of the particles can be obtained in equation:where R0 is the final radius of the particles.
4.3. Boundary Condition of the Particles and the Control of Wall
In PFC2D, the calculation area needs to be determined first, and the modeling calculation cannot be carried out outside the calculation area [16]. The wall needs to be built at the boundary of the calculation area. The moving speed of the wall was controlled by the servo control mechanism to update the position of the base point of the wall. The wall cannot be assigned the force values directly, and only the velocity, angular velocity, and rotation center of the wall can be specified to realize the simulation of uniaxial compression. From an energy point of view, the final external energy values exerted by constant velocity loading and constant load loading are the same, so the loading speed of the standard wall is in line with the actual experimental results. In the simulation of PFC uniaxial compression, the stress of the wall was as follows:
Fw is the force exerted on the wall by a single particle; d is the height of the model; L is the width of the model; and N is the number of particles. The relationship between the moving speed of the wall and the stress on the wall is controlled by the parameter G. is the moving speed of the wall; the contact stiffness is related to the relative displacement by the following equation: is the normal stiffness of particles, is the number of particles, and is the normal displacement.where is the tangential stiffness, and the tangential force is controlled by displacement increment.
In the process of Uniaxial Compression Test (UCT), reasonable selection of wall stiffness can increase the accuracy of simulation results [16]. If the wall stiffness is too large, the initial contact force between particles will be too large, which increases the time for the model to achieve the initial equilibrium; if the wall stiffness is too small, particles will pass through the wall and the simulation fails. It is reasonable to set the wall stiffness as 10 times of the particle stiffness according to previous studies [15].
4.4. The Calibration of the Particle Contact Parameters of the Concrete Blocks at Different Temperatures
It was necessary to calibrate the particle contact parameters of the PFC2D model before the explosion impact simulation. In our model, the aggregation of particles bonded together by the pb-bonding model. For the concrete structure, its compressive strength, shear strength, and other macro mechanical properties can only be determined by defining the interaction between circular particles [17]. PFC2D provides two methods for value assignment: the first is direct assignment and the second is indirect assignment. It includes effective modulus (e-mod), normal shear stiffness ratio (kratio), radius multiplier (Pb radius), friction coefficient (fric), and so on. In order to study the damage of concrete under different temperatures and blast waves more effectively, the second method was used in this study. The mechanical parameter between particles in the parallel bond model can only be determined by repeated calibration in the simulation process when the simulated stress-strain curve matches the experimental one.
After the concrete test blocks were prepared in the laboratory, the universal hydraulic testing machine in the laboratory was equipped with upper and lower temperature-control copper plates to adjust and control the temperature. The uniaxial compression of the same batch of concrete blocks which take the PFC2D aggregate model was tested at room temperature, 100°C, 200°C, and 300°C. The reason why we chose 300°C as the maximum test temperature is because the concrete block will be considered totally broken after this boundary, and it is unnecessary to consider the shock wave bearing capability of these types of concrete blocks. In order to keep the temperature of uniaxial compression within the specified range, it is necessary to stick high-temperature insulation cotton (black color) around the concrete blocks.
In order to ensure the accuracy of the experiment, three groups of experiments using the same batch were carried out for each test block. As a kind of thermal inert material, the thermal conductivity of concrete is poor, so it is necessary to keep the heating time at 2 hours to ensure that the whole concrete test block reaches the preset temperature before the uniaxial compression. When the concrete was damaged under compression, the side spalling was obviously observed, and the main crack spread through the test block. Shear failure was the main failure mode in numerical simulation as well as uniaxial compression of concrete. The stress-strain curve was recorded for each set of concrete test sample (Figures 12–15).
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
The UCT results simulated by PFC2D manifested that the deformation and slip of the mortar particles (ball-ball) and aggregates (clump) in the model were obvious during the damage, and the internal crack progression from the opening to the end through the damage can be clearly detected. At the initial stage of loading, due to the small stress, micro cracks appeared around some particles and they mainly existed at the interface between cement matrix and coarse aggregate. When micro cracks appeared between particles, stress redistribution occurred in the whole concrete blocks, and stress concentration was prone to form around the cracks, which leads to the continuous expansion of cracks. As the loading continued, the initial interface cracks gradually extended and widened, new cracks appeared in other aggregate interfaces, and the cracks gradually connected with each other. When the concrete specimen reached the maximal loading, the particle displacement increased obviously, and many macro cracks appeared between particles. As the macro cracks continued to expand, the bearing capacity of residual concrete gradually decreased, and oblique cracks formed between adjacent cracks (Figure 16).
It should be further noted that the meso parameter of the numerical model established by PFC cannot be directly equal to the macro parameter of the real material. The main reason is that the effective modulus and Poisson’s ratio obtained by indirect characterization methods are not equivalent to the actual elastic modulus and Poisson’s ratio, but they are indeed positively correlated. Based on the comparison of the numerical simulation and the laboratory tests, the digital model is in good agreement with the final failure mode, which indicates the rationality of the calibration parameters.
During the calibration process, the ball-ball contact parameters need to be adjusted continuously. After the stress-strain intensity curves of the test blocks were obtained from the control experiments, the calibrated stress-strain curves at different temperatures were obtained. For the three different temperatures, the key contact parameters of concrete mortar and aggregate in the ball-ball bond model were recorded, respectively, as shown in Tables 6 and 7.
Two parameters (e_mod and pb_emod) were used to characterize the elastic modulus, and k_ratio was used to characterize the Poisson’s ratio; between the particles, pb_coh is the cohesive force, pb_fa is the friction angle, and fric is the friction coefficient.
In our calibration process, the internal coefficient of aggregate was kept unchanged, and the stress-strain capacity of concrete block was changed by varying the friction coefficient and cohesion between the particles in cement mortar. The uniaxial compression model was set up in PFC2D to calibrate the concrete blocks at different temperatures and the output curve was compared with our experiment results in Figures 12–15.
4.5. Explosion Impact Simulation of Concrete Blocks in PFC2D with Variant Temperatures
The insensibility of the temperature in the concrete lining structure has been classified into three groups (100°C, 200°C, 300°C). After equal sampling of the compressive wave curve, 100 data points were generated and incorporated into PFC2D for impact loading. Combining it with the temperature distribution, at the first 5–15 m, the compressive pressure was relatively high and the temperature gradually increased from the ambient temperature to over 900°C (most of the downstream ceiling). Several critical points were studied to compare the meso-scale explosion damage with and without the variant temperatures (Figures 17–19).
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
By comparison, the damage of the concrete block was obviously enlarged due to the high temperature after being impacted by the same compressive wave.
From the two monitoring points, B0 and B4, which are slightly away from the origin of explosion shock, the shock wave showed a secondary peak, and this phenomenon is completely consistent with the propagation rule of air shock wave in the tunnel.
It can be seen from Figures 20–22 that the damage was mainly located at both sides of the surface and then developed into the depth. The damage in the central part was mainly due to the diagonal interface of the aggregates, which leads to the expansion of the cracks on both sides to the center. In the process of simulation, it can be clearly seen that the intensity of damage was obviously related to the peak value of the stress wave.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
The failure modes of experiments and PFC2D simulation models can be divided into three types. First, the aggregate is broken, which has a low probability and generally occurs on the path of the main crack, indicating that the instantaneous pressure at this stage of stress density exceeds the strength of the aggregate itself. Second, the internal failure of cement mortar usually occurs when the mortar strength is low. Third, the breakup of binding connection on the interface between aggregate and cement mortar is very common. Moreover, the crack development direction will change when encountering aggregates in the cracking path. Therefore, the section strength of the micro cement mortar and aggregate ultimately determines the final strength of the macro concrete blocks.
According to our simulation results, the tunnel concrete lining can be divided into the area of explosion damage, the area of high temperature damage, and the area of mixed damages. The boundaries of the three zones are influenced by the internal structure of the tunnel, slope, ventilation speed, explosion equivalent, and the distance between vehicles and the lining structure.
5. Conclusions
(1)The damage process of the concrete material subjected to explosive impacts at different temperatures has been simulated by the discrete element method. The aggregate generated by the clump method can truly reflect the inherent characteristics of aggregates and mortars in concrete. The damage process of the concrete blocks is a dynamic process of energy accumulation and release. Our model can intuitively visualize the fracture originated at the weak interface layer. In the process of PFC2D calibration, the parallel bond model of cement mortar directly affects the peak stress of the stress-strain curve of concrete, while the strength of aggregate itself has little effect, which is consistent with our experimental observations.(2)The damage of tunnel lining structure after explosion and combustion is mainly caused by high temperature, overpressure impact, and gases such as carbon monoxide. The increase of the temperature has amplified the compress damage significantly.(3)More attention should be paid to the improvement of vault and arch waist when repairing the middle and rear parts of the lining structure [18].(4)The simulation results of the discrete element method combined with LS-DYNA and FLUENT can more accurately determine the damage degree of the concrete blocks. These methods will be helpful for providing valuable information needed in making repair decisions and predicting the residual lifetime of a damaged tunnel.Data Availability
All the data, models, and code generated or used during the study are included in the submitted article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The research presented in this paper was funded by the Open Fund Project of State Key Laboratory of Disaster Prevention in Civil Engineering, Tongji University, China.