Analytical study of local damage on concrete slab subjected to a contact explosion

2.1. Model assumption and initial condition

According to local damage tests on concrete slabs under explosive loads, as shown in Fig. 2 [20], there are some typical phenomena: for a thick slab, shown as Fig. 2(a), there will be a compression crater on the front-face when it is subject to a contact explosion, but there will not be any irreversible deformation on the back-face. This kind of phenomenon is called explosive cratering, and it is not necessary to consider reflective wave tension effects from the back-face of the thick slab. For a slab which is not thick enough (Fig. 2(b)), when it is subjected to a contact explosion, a compression crater will form on the front-face and a spalling crater forms on its back-face, but they are isolated from each other. This is known as spalling. For a thin slab (Fig. 2(c)), a compression crater forms on the front-face and a spalling crater forms on the back-face: the two are mutually connected. This kind of phenomenon is known as explosive perforation.

A model which depicts the impulse on the surface of a medium under contact explosion load is shown in Fig. 3. To make the mathematics tractable, we will make some simplifying assumptions:

a) The chemical reaction is completed instantly, and shock waves are quasi-instantaneous jumps in all the variables, from their initial, to final, states.

b) The radius of the spherical charge is $R_w$, and the density of the charge is ${\rho }_w$. It is in contact with the intended concrete targets.

c) Under contact explosive loads the concrete will reach a plastic state instantly and dissipate explosive energy through large deformation. Compared with the plastic deformations suffered, the elastic deformation can be ignored.

d) The impulse of the explosive detonation transmitted to the target is $I$, and it is firstly transmitted as a hemispherical distributed impulse (over the FGF area), which is named as initial expansion body.

Fig. 2Photographs of local damage tests on concrete slabs under explosive loads [20]

a) Cratering on the front-face

b) Spalling on its back-face

c) Perforation

Fig. 3The impulse on the surface of medium under contact explosive loads

a)

b)

According to the law of the conversation of momentum, the vertical impulse of the detonation is equal to the vertical momentum of the initial expansion body:

where, $M$ denotes the mass of the initial expansion body, and $M=\frac{2}{3}\pi \rho a^3$; $\rho$ denotes the density of target; $v_0$ denotes the radial velocity thereof.

We assumed that the dimensions of the target are much bigger than the radius of the charge, and the kinetic energy can be written as:

Based on the conservation of energy, the relationship between the initial expansion body and the effective explosive energy can be shown to be as follows [21]:

where, $Q_v$ denotes the specific energy coefficient of the explosive, $\eta$ is the fractions of energy radiated into the target, which can be acquired as follow.

The theoretical solution of the problem of radiation of the energy of a contact explosion is an extremely difficult one. It is characterized by the fact that it takes place at the interface of two media: air and a solid half-space. However, the problem can be solved experimentally in a relatively simple fashion. V. D. Alekseenko [22] analysed lots of experimental data and summarized them in Table 1, which indicate that the material of targets and also the position of the charge relative to the free surface have a significant effect on the distribution of the energy of a contact explosion.

Based on theoretical analyses, O. V. Nagornov and V. E. Chizhov [23] found the energy fraction of a contact explosion:

where ${\rho }_{01}$, ${\gamma }_1$ are respectively the density and the adiabatic index of the upper medium (air), while ${\rho }_{02}$, ${\gamma }_2$ are respectively the density and the isentropic adiabatic index of the lower solid materials. However, for materials with $\gamma >$ 3, this formula is not applicable.

Generally, it is possible to obtain the fractions of energy radiated into (pumped into) the concrete target which is subjected to a contact explosion.

Therefore, it is possible to calculate the local damage effect of the target subjected to a contact explosion.

2.2. Basic theoretical methods

Based on the conversation of energy, internal energy, which is dissipated in the concrete matches the external work done, as is supplied by the explosion. Material yields and fractures, pieces rub and slide against each other, which dissipates energy. The deformation and displacement in the plastic stage are much greater than those in the elastic stage, therefore it is reasonable to ignore the energy dissipated in the elastic stage. In the area near the explosion source, material undergoes large irreversible deformation, which in addition to the presence of the free surface, means that the elastic displacement can be ignored [24]. Therefore, it is acceptable to use a rigid-plastic model to describe the dynamic behaviour of concrete in such cases.

According to the conservation of mass, the incompressibility of the material, and the applied boundary conditions, it is possible to build a permissible kinematic velocity field which can characterise the deformation and displacement in the plastic zone. Combined with extremum principles and energy equilibrium rules [24], it is possible to deduce an expression for material resistance to explosive load on slabs with different thicknesses. Based on the non-dimensional expressions of material resistance, curves of material resistance may be plotted to find the threshold radius of the crushed zone which leads to spalling or punching. Combined with the initial conditions, the final crushed zone radius can be calculated by integrating the equations of motion for the expansion process.

It can be supposed that the material resistance is $P$, so the equations of motion can be written as:

where $M$ is the mass of the initial expansion body, and $R$ is the radius of the crush zone. Combined with the initial conditions of the expansion body:

The final radius of crush zone can be calculated, similarly, the two kinds of threshold radius of crush zone can be obtained. It is possible to predict what kind of damage effects will occur in certain thicknesses by comparing the final radius and the two threshold radii.

The material resistance can be acquired by seeking the limiting load. When the load reaches the limit for an ideal plastic-rigid body (with a yield plateau), free plastic flow occurs. According to extremum principles and energy equilibrium requirements, the rate of work could be written as follows [24, 25]:

where $S$ and $V$ respectively are the areas and the volumes of the deformed body, $P_i^k$ are components of surface forces, $v_i$ are components of the possible kinematic velocity fields, ${\xi }_{ij}^k=\frac{1}{2}\left(\frac{\partial {\upsilon }_i^k}{\partial x_j}+\frac{\partial {\upsilon }_j^k}{\partial x_i}\right)$ are components of strain rates, ${\sigma }_{ij}^k$ are components of the actual stress, $S_{\tau }$ are surfaces of velocity discontinuity, ${\tau }_s$ is the shear yield strength of material, and $\left[v_{\tau }^k\right]$ denote the jump $v_{\tau +}^k-v_{\tau -}^k$ on $S_{\tau }$.

By virtue of the Saint Venant-von Mises equation, we know that vectors ${\sigma }_{ij}^k$ and ${\xi }_{ij}^k$ will be parallel and the quantity ${\sigma }_{ij}^k{\xi }_{ij}^k$, being the scalar product of two parallel vectors, will be equal to the product of their moduli, i.e. ${\sigma }_{ij}^k{\xi }_{ij}^k={\tau }_S{H}^k$, where ${\tau }_S$ is the yield limit for shear, $H^k$ is the intensity of shear strain-rate. Then Eq. (7) can be written as:

According to Eq. (8), the upper bound for the limit load of $P_i^k$ could be calculated from the permissible kinematic velocity field. Therefore, the problem becomes one of how to establish an accurate kinematic velocity field from the experimental data.

3.1. Material resistance of thick slabs under point-source explosion load

When the slab is sufficiently thick, compression craters are likely to be generated on the front-face while no deformation is supposed to occur on the back-face under the applied contact explosion load. This can be regarded as a point-source explosion on the surface of the half-space medium. In this paper, a rigid-plastic model is used to describe the crushing mechanism as a shear failure as dictated by Eq. (7). According to experimental observations and slip line field theory [25], a kinematically admissible velocity field (Fig. 4) can be established and the medium can be divided into three regions.

As shown in Fig. 4, the zone (FGF) nearest the explosion source is an initial expansion body with radius $a$ and average radial velocity $v_0$. The medium in region I is in a state of plastic flow. By combining the initial conditions, and according to mass conservation and the incompressibility of the material, $v_R=a^2/R^2{v}_0$ and $v_{\phi }=v_{\theta }=0$ can be acquired. According to the given velocity field, we can deduce the strain rate field and the strain rate is computed as $H^1=2\sqrt{3}{a}^2{v}_0/R^3$.

Fig. 4A kinematically admissible velocity field for thick slab

The kinematically admissible velocity fields in regions II and III can be calculated according to plasticity theory. The plasticity theory proposed by Hill [25] presents solutions to plane problems, and can describe the horizontally cylindrical charge effect on the surface of a material. According to symmetry, the velocity slip line field of the plane problem is rotated along the $Z$-axis to obtain the kinematically admissible velocity field under the action of a point-source explosion. As the rigid-plastic model complies with the associated flow rule, the strain slip line is overlapped with the velocity slip line so that the slip line in the model is not only the trajectory of the medium but also the direction of the shear stress (Fig. 5).

Fig. 5The velocity fields in regions II and III

Given that the curve BEB is the collision contact surface, and the change of its radius shows geometrical similarity with the expansion process, and as there are no shear stresses on BE, the slip line forms an angle of 45° with the surface. All $\beta$ are straight lines and the tangent to the virtual circle with centre of O and radius $R/\sqrt{2}$ are developed: $\alpha$ is the involute of this virtual circle. In region I, the expanding speed of the radius $R$ of the crushed zone is $v_R=a^2/R^2{v}_0$. Then the velocity along line $\alpha$ on BE is $v_{\alpha }=\sqrt{2}{v}_R$, while those on lines $\alpha$ and $\beta$ in region BDE are $v_{\alpha }=\sqrt{2}{v}_R$ and $v_{\beta }=0$ separately. Thus, it can be known that region BCD moves along DC at a speed of $v_{\alpha }=\sqrt{2}{v}_R$.

Exploiting symmetry, the right-hand part of the slip line field is acquired (Fig. 6). The radius of the virtual circle is $r=R/\sqrt{2}$ and it is the involute of the cluster line $\alpha$; DE and BD are the outer boundaries of region BDE, and belong to clusters $\alpha$ and $\beta$, respectively. If the expressions for curves DE and BD are known, we can calculate the volume of the region BDE after it is rotated around the $Z$-axis and the surface area corresponding to DE.

The curve ED is the involute of the virtual circle presented in Fig. 6. By randomly selecting a point $M\left(x,y\right)$ on the curve to make MF tangential to the circle at $F$, then $\angle xOF=\phi$, $OF=r$ and $KF=MF$. The parametric expression of curve ED can be acquired according to the geometric relationships:

Fig. 6The schematic boundary lines of plastic regions II and III

By drawing a tangent line $EN$ to the virtual circle at point $E$ and letting $OE=R$ and $ON=R/\sqrt{2}\text{,}$ then $EN=R/\sqrt{2}$ and $\angle EON=\pi /4$ can be acquired. Moreover, $\widehat{KN}=EN=R/\sqrt{2}$, therefore $\phi =$ 1 and $\gamma =1-\pi /4$. The deflection angle between the $Z$-$X$ plane and the $x$-$y$ plane is $\gamma =1-\pi /4$. Similarly, $BD=r\pi /2$. In the $Z$-$X$ plane, let $OM$ be the radius $\rho$ in polar coordinates and $\angle ZOM$ be the deflection angle $\phi$ measured from the $Z$-axis, then:

Given that the rate of cavity expansion along the line $\alpha$ is $u=\sqrt{2}{v}_R$ ($v_{\alpha }=\sqrt{2}{v}_R$), the velocity component in the polar coordinate can be obtained as follows:

By rotating the plastic zone around the $Z$-axis, we express the plastic zone in spherical coordinates with rotation angle $\theta$. The velocity fields in all regions, in spherical coordinates ($\rho ,\phi ,\theta$), can then be obtained.

For the kinematically admissible velocity field in region II: $v_{\rho }^{II}=v_{\alpha }/\sqrt{1+{\varphi }^2}\text{,}$$v_{\phi }^{II}=\varphi v_{\alpha }/\sqrt{1+{\varphi }^2}$ and $v_{\theta }^2=0\text{;}$ as to the strain rate field: ${\xi }_{\rho }^{II}=-v_{\alpha }/r(1+{\varphi }^2)\text{,}$${\xi }_{\theta }^{II}=v_{\alpha }/r(1+{\varphi }^2)$ and ${\xi }_{\varphi }^{II}=v_{\rho }/\rho +\partial v_{\varphi }/\rho \partial \varphi =v_{\alpha }/r(1+{\varphi }^2)+v_{\alpha }/r(1+{\varphi }^2{)}^{1/2}{\varphi }^2$; the shear strain rate is $H^{II}=v_{\alpha }\sqrt{2(3{\varphi }^4+2{\varphi }^2(1+{\varphi }^2{)}^{1/2}+1+{\varphi }^2)}/r{\varphi }^2(1+{\varphi }^2)$. The velocity difference on the surface of the tangential velocity discontinuity in region II is constant at $\left|\left[v_{\tau }^{II}\right]\right|=\sqrt{2}{v}_{\alpha }$. The surface area of the velocity discontinuity is $S^{II}=8.956\pi r^2$. The rectangular form of the velocity field in region III can be expressed as $v_Z^{III}=-v_R$, $v_X^{III}=v_R$ and $v_{\theta }^{III}=0$. The strain rate field and the shear strain rate can be represented by ${\xi }_Z^{III}={\xi }_X^{III}={\xi }_{\theta }^{III}=0$ and $H^{III}=0$. The velocity difference on the surface of the tangential velocity discontinuity in region III is constant at $\left|\left[v_{\tau }^{III}\right]\right|=\sqrt{2}{v}_{\alpha }\text{.}$ The surface area of the velocity discontinuity is $S^{III}=1.745\pi r^2$.

By substituting the physical quantities obtained above into Eq. (8), the material resistance on the surface of the half-space medium under application of a point-source explosion can be expressed in non-dimensional form as follows:

3.2. Calculation of the resistance of the limited-thickness slab subjected to point-source explosion

When a point-source explosion occurs on the surface of a limited-thickness slab, the back face of the target is affected. When the radius $R$ of the crushed zone reaches a certain value, the target material begins to be moved outwards from the back face instead of the front-face when it is subjected to contact explosion loading. Fig. 7 shows the regional division of a limited-thickness slab subjected to a point-source explosion. Region I is the crushed zone around the explosion source, while region II is the fracture zone caused by reflected tensile waves on the back face. Region III is the elastic deformation area in which deformation and displacement can be neglected.

Under such circumstances, the medium in region I flows towards region II. Considering the hindering effect of the medium in region II on the expansion body (the medium within AOA), region ACCA is considered as an expansion body with an initial velocity $v_R$ and the radius of the impact contact line AA is $R$. Here, shear yield is supposed to occur to the material so that line AB (the regional discontinuity) forms an angle of $\pi /4$ with the back face. There are two kinematically admissible velocity fields in region II: in the first, the velocity components in all directions change with position and therefore the medium is plastically deformed while in the second velocity field, the velocity components in all directions are position-independent and the medium moves as a rigid block.

Fig. 7The regional division of a limited-thickness slab subjected to point-source explosion

3.2.2. The second kinematically admissible velocity field

The second kinematically admissible velocity field only changes with time, shown in Fig. 8. The media, in different positions, have the same velocity, which suggests that rigid blocks are directly forced out from the funnel region by the expansion body. Initially, the rigid blocks have different velocities to the expansion body. While with continuing motion, the rigid blocks tend to match the velocity of the expansion body. This process is so short as to be negligible. Therefore, the initial velocity of rigid blocks is regarded as equal to that of the expansion body. The medium is mainly affected by the frictional force on the surface of each discontinuity.

Fig. 8The regional division of a thin slab subjected to point-source explosion

As the velocity component in the region II is $v_r=v_R$, $v_{\phi }=0$, $v_{\phi }=0$, then $H^k=0$ and $\left|\left[v_{\tau }^k\right]\right|=v_R$. The surface discontinuity forms a part of the conical surface: $S=\sqrt{2}\pi (L^2-R^2)$. By substituting these parameters into ${\int }_S{P}_{ni}^k{v}_{i}dS={\tau }_s{\int }_V{H}^{k}dV+\sum _{k=1}^N{\tau }_s{\int }_{S_k}\left|v_{\tau }^k\right|d{S}_k$, an equation similar to Eq. (13) is obtained as follows:

14

$\frac{P^k}{{\tau }_s}=\frac{\sqrt{2}}{2(1-\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_0)}\left({\left(\frac{L}{R}\right)}^2-1\right).$

4.1. Non-dimensional of explosion resistance

By non-dimensional Eqs. (12), (13), and (14), the resistance to a point-source explosion on the surface of the half-space is obtained, as follows, by using the given radius $a$ of the charge and $\tilde{L}=L/a$:

The resistance of a slab with limited thickness to a point-source explosion in the first velocity field is expressed as:

The resistance of a limited-thickness slab to a point-source explosion in the second velocity field is:

Fig. 9Material resistance curves

a) Perforation before spalling

b) Explosive punching

By regarding $P^k/{\tau }_s$ as a function related to $\tilde{R}$, the $P^k/{\tau }_s$-$\tilde{R}$ curves in the three velocity fields can be calculated. Suppose that the slabs are $L=6a$ and $L=1.5a$ thick, as shown in Fig. 9(a) and 9(b) separately: curves 1, 2, and 3 represent the impact resistances in the half-space velocity field as well as the first and second velocity fields of the medium with limited thickness, respectively. Curves 1 and 2 intersect at point A, showing that when the radius of the fracture zone extends to ${\tilde{R}}_A$, free surfaces are formed by the expansion body; curves 2 and 3 intersect at point B, meaning that when the radius of the fracture zone extends to ${\tilde{R}}_B$, the upper media are forced out from the ground surface by the expansion body and plugging occurs. For each $\tilde{R}$ value, three impact resistance values can be found. Thereinto, the minimum values are selected as the actual resistances, namely, the values from curves 1 and 2 to the left and right of point A, respectively.

In the explosion, the sum of the material resistance to the surrounding plastic zone subjected to impact on the contact surface is expressed as follows:

where, $S_r$ is the area of the impact contact surface. In the half-space, it can be regarded as $S_r=2\pi r^2$ without considering the influence of the free surfaces. While in a limited thickness slab, $S_r=2\pi r^2(1-\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_0)$, ${\phi }_0=\pi /4$.

4.2. Method of calculating the threshold thickness on spalling

Based on the impact resistance curves, the limiting material resistance can be calculated. By combining the initial conditions to solve the equation of motion, the radius of the crushed zone can be computed. In fact, the wave problem is converted to, and solved as, the resistance. Under threshold spalling conditions, the explosion resistance is $P_1=p^{k}2\pi a^3{\tilde{R}}^2$ and the mass $M_1=M=\frac{2}{3}\pi \rho a^3$ of the initial expansion body is a constant. By combining Eqs. (5), (15), and (18), the non-dimensional equation of motion can be obtained as follows:

The initial conditions are:

Let ${\omega }_1=6.92\pi {\tau }_{s}a/M_1$ and ${\omega }_2=11.68\pi {\tau }_{s}a/M_1$, the equation of motion is converted to:

As this differential equation, cannot be solved directly, the equation needs to be converted to $-\tilde{v}\frac{d\tilde{v}}{d\tilde{R}}={\omega }_1\mathrm{l}\mathrm{n}\left(\tilde{R}\right){\tilde{R}}^2+{\omega }_2{\tilde{R}}^2$ form and then, separating the variables and integrating gives:

Processing Eq. (22), we obtain:

Let $\xi =\frac{4.5}{{\omega }_1}{\left(\frac{v_0}{a}\right)}^2=\frac{0.65M_1{v}_0^2}{\pi {\tau }_s{a}^3}$, then solving Eq. (23) gives:

The thickness of the medium at intersection point A can be calculated by using Eqs. (15) and (16):

By substituting Eq. (24) into (25), the threshold spalling thickness can be acquired:

Eq. (26) is a function of the non-dimensional factor $\xi$. According to Eq. (2), the relationship between the non-dimensional factor $\xi$ and the mass $Q$ of the charge can be obtained:

According to the physical model of a contact explosion (Fig. 3), $R_w=a$ can be obtained so that $\xi =1.73\eta Q_v\delta {\rho }_w/{\tau }_s$. By performing dimensional analysis, $\xi =\lambda \left(\frac{{\rho }_w}{\rho }\right)\left(\frac{{\sigma }_d}{{\tau }_s}\right)\left(\frac{Q_v\delta }{C_p{v}_m}\right)$. As a non-dimensional combination of the polytropic index, it represents the internal relationship between the density of the charge, explosion heat, effective energy ratio, the mechanical equivalent of heat, the density of the material, its shear yield strength, peak pressure sustained, deformation wave velocity, and maximum particle velocity.

There are LambertW functions, recorded as $W\left(x\right)$, in Eqs. (24) and (26). When the variable $x>-e^{-1}$, $W\left(x\right)$ is an increasing function with real, unique solutions. Therefore, the solutions of all the functions in this paper are real, single-value, numbers.

By converting the threshold spalling thickness in its non-dimensional form to its scaled distance form, it is found that:

Under the condition in which the parameters of the charge and media are determined, the non-dimensional factor $\xi$ is also determined. By combining Eq. (26), this scaled distance can be calculated. When the model size is changed, it is a similar constant that serves as a basis for determining whether, or not, spalling occurs. For small-scale underground chemical explosions where the gravitational acceleration can be ignored, this result corresponds to the similarity law governing explosions.

4.3. Method of calculating the threshold thickness on perforating

To calculate the threshold perforating thickness, it is necessary to obtain the displacement and velocity of point A from the impact resistance curves. At this time, owing to the appearance of the back surfaces (free surfaces) when affected by the expansion body, the perforation body begins to be formed. The material resistance in the equation of motion is calculated using the first velocity field in the medium with its limited thickness. When the radius of the crushed zone develops to reach point B, perforation occurs. According to the relationship between the radius of the crushed zone at point B and the slab thickness, the threshold perforating thickness can be calculated.

According to the mechanical model for explosive perforation, the explosion resistance and mass of the expansion body can be obtained as $P_2=p^{k}2\pi a^2{\tilde{R}}^2(1-\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_0)$ and $M_2=2\pi \rho a^3(1-\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_0)/3$, respectively. By combining Eqs. (5), (16), and (18), the equation of motion when the perforation body begins to be formed can be acquired:

The initial condition is:

and the threshold perforating condition is:

Let ${\omega }_3=17.63{\tau }_s/\rho a^2$, the equation of motion may be converted to the following form: $\ddot{\tilde{R}}=-{\omega }_3\mathrm{l}\mathrm{n}\left(\frac{\tilde{L}}{\tilde{R}}\right){\tilde{R}}^2$. By separating the variables and then integrating $-\tilde{v}\frac{d\tilde{v}}{d\tilde{R}}={\omega }_3\mathrm{l}\mathrm{n}\left(\frac{\tilde{L}}{\tilde{R}}\right){\tilde{R}}^2$, from the initial state to the threshold perforating state, we get:

Based on Eq. (22), we can obtain the relationship between the velocity and displacement under the initial conditions at which perforation begins to occur:

By combining Eqs. (30) and (31), the following formulae can be acquired by removing the term in $v_A$:

According to the aforementioned definitions of ${\omega }_1$, ${\omega }_2$, ${\omega }_3$:

We can obtain:

Based on Eq. (25), the relationship between the radius of the crushed zone and the slab thickness at point $A$ (in its initial state) may be found to be:

In addition, according to Eqs. (16) and (17), the relationship between the radius of the crushed zone and the slab thickness at point $B$ at the threshold perforating state can be acquired:

By combining Eqs. (32) to (35), an equation relating $\tilde{L}$ with minimus quantities being ignored can be acquired by sorting the equation of motion of that perforation formed by the explosive expansion body on the back-face of the slab and the initial and boundary conditions:

In Eq. (36), $\xi$ represents the non-dimensional impact factor and has an identical meaning to that used in the calculation of the threshold spalling thickness. The root of this equation is the threshold perforating thickness that can simultaneously satisfy the equation of motion for perforation, the initial conditions, and boundary conditions. However, as analytical solutions cannot be obtained, we calculate an approximate solution using a numerical method. By combining Eqs. (28) and (36), the relationship between the threshold spalling thickness and the non-dimensional impact factor $\xi$ (in scaled distance form) can be acquired.

5.1. Theoretical analysis of two kinds of threshold thicknesses

Suppose that the shear yield strength of concrete slabs ranges from 2 MPa to 6 MPa; the density of the TNT charge is ${\rho }_w=$ 1630 kg/m³; the heat liberated in the explosion is $Q_v=$1010 kcal/kg; the mechanical equivalent of heat is $\delta =$4187 J/kcal. According to Table 1, the fraction energy radiated into the slabs is 0.1, i.e.$\eta \approx$ 0.1.

By combining Eqs. (26), (27), (28), and (29), the threshold thickness of spalling and perforating have been calculated and listed into Table 2. Besides, the plot of $\tilde{L}~\xi$ in non-dimensional form, and that of $L/Q^{1/3}~\xi$ in scaled distance form can be acquired as shown in Fig. 10(a) and 10(b), respectively.

Fig. 10The relationship between threshold thickness and non-dimensional impact fact

a) The threshold thickness in non-dimensional form

b) The threshold thickness in scaled distance form

As shown in Fig. 10, the threshold spalling thickness and the threshold perforating thickness increase with increasing non-dimensional impact factor $\xi$. The larger the value of $\xi$, the larger the threshold thicknesses, and vice versa. The non-dimensional impact factor is a comprehensive index covering the density of the charge, the heat liberated by the explosion, the effective energy ratio, and the shear yield strength of the material. It is proportional to the impact load or energy transmitted from the explosion source to the medium, but is inversely related to the strength, or resistance, of the medium. The higher the impact load applied by the explosion, the thicker the material must be to resist impact; the higher the strength of the material, the greater the explosion-resistance of the material, the thinner it need be to resist the impact.

After calculation, it can be found that the threshold spalling thickness $L_z$ is about twice the threshold perforating thickness $L_P$, which agrees with the relationship $L_z=2L_P$ as shown in the BRL formula and the revised version thereof.

5.2. Numerical analysis of two kinds of threshold thicknesses

5.2.1. Numerical model

In order to verify the accuracy of theoretical results, the numerical analysis was performed with the hydrocode AUTODYN v12.0 [27]available from ANSYS. In the numerical simulation, the concrete slabs subjected to contact explosions are modelled with an axial symmetric two-dimensional model with an element size of 5 mm [28], as shown in Fig. 11. The slab is modelled by a Lagrange subgrid, in which the coordinates move with the material. In order to describe the infinite boundary of slab, the lateral Lagrange subgrid is set as a non-reflecting boundary, which eliminates the reflections of wave from the boundary. On the other hand, the air and high explosive are modelled by the Euler subgrid, wherein the grid is fixed, with material allowed to flow through [29, 30]. The Euler-Lagrange interface interaction is considered, which means that the Lagrange subgrid imposes a geometric constraint on the Euler subgrid whereas the Euler subgrid provides a pressure boundary to the Lagrange subgrid [31]. The boundary condition of the Euler subgrid is set as an outflow boundary. The erosion technique is adopted to simulate the damage process.

Fig. 11Numerical model: a) different materials in the model; and b) numerical mesh.

Combined with Table 2, the theoretical results of threshold thickness of slab under contact explosions have been calculated and shown in Table 3, in which the shear yield strength of material is 6 MPa and the mass of charge is 0.853 kg.

Table 3Threshold thickness by theoretical calculation

Shear yield strength of material (MPa)	6
Radius of charge (mm)	50
Mass of charge (kg)	0.853
Threshold spalling thickness (m)	0.416
Threshold perforating thickness (m)	0.226

In numerical models, the charge of explosives and material properties of slab are set according to Table 3, which are both constant for all models. But the thickness of slab is different in each model, which varies from 120 mm to 520 mm. More detailed properties for numerical model have been listed in Table 4, 5, 6, and 7.

Table 4Numerical model

Slabs

Material: concrete Shear yield strength: 6 MPa Thickness: 120 mm, 220 mm, 320 mm, 420 mm, 520 mm Subgrid: Lagrange

Charge of explosives

High explosive: TNT Density: 1630 kg/m3 Radius: 50 mm Mass: 0.853 kg Subgrid: Euler

5.2.5. Results and discussion

There are a series of numerical models have been built, and the numerical results of damage contour are shown in Table 8.

The numerical results (in Table 8) show that as the thickness increases, the damage of the slab decreases, which is reasonable and consistent with common sense. Due to the thicker slab can stand higher loads, under the same blast loads, the lighter damage will occurs. In order to discuss the degree of damage, it is essential to define some categories of damage. Even though there are lots of experimental and numerical study on damage of slabs subjected to blast loads, the available literature does not provide a definition for the damage criterion under which the spallation or perforation failure occurs. Besides, most of the literatures [35, 36] are concerned on the spall damage and fracture of ductile and brittle metals, which are less about the construction materials such as concrete and rock [37, 38].

Lu Fangyun et al. [34] divided the damage of one-way square slabs into the following categories: (a) low damage, with only small cracks in both surface of the slab; (b) moderate damage, with spallation occurring on bottom of the slab; (c) high damage, with perforation from the upper to the bottom surface; and (d) collapse damage, which is the punching failure of the slab. Due to their works focus on the structural holistic performance of slabs subjected to standoff blast loads, they defined damage criteria correlated to the displacement or the support rotation angle. However, these categories cannot be applied to judge the local damage effects of slabs under contact explosions, which are much different from the former.

According to McVay M. K. [39], the spall damage of concrete slabs can be divided into three categories: (I) Light damage: from initial state to a few barely visible cracks; (II) Threshold for spall: from a few cracks and a hollow sound to a large bulge in the concrete with a few small pieces of spall on the bottom surface, which is equivalent to threshold spallation in this paper; and (III) Medium spall: from a very shallow spall to spall penetration up to one third of the slab thickness. But he did not define the more serious damage, which was perforation or punching. Therefore, it is essential to add two damage categories: (IV) Threshold for perforation: from a spall to perforate the slab, which means the spall crater connects to the blast crater; (V) Perforation or punching: blast crater on the top of slab extends to the bottom, which means a big hole occurs on the slab.

Therefore, combined with the damage contour of numerical results, it is possible to classify the damage degree of each target under the same blast load (0.853 kg TNT) according to categories of spall damage of concrete slabs, which has been listed in Table 8.

As shown in Table 8, the colour scale describes the degree of damage, where red represents the most serious damage. The thickness of No. 2 and No. 4 targets are respectively threshold perforating thickness and threshold spalling thickness of slab under contact explosion (0.853 kg TNT), which are very close to the theoretical calculation by this paper, shown in Table 3. To an extent, it can proved that the method of this paper to calculate the local damage of slabs subjected to a contact explosion is reasonable.

Table 8Damage of numerical results

Targets	Thickness (mm)	Damage contour		Local damage on the bottom	Categories
No. 1	120			Blast crater extends to the bottom of the target	Perforation
No. 2	220			Damaged area just connected to the blast crater	Threshold perforation
No. 3	320			Damage or failure up to one third of the slab thickness	Spallation
No. 4	420			A few small pieces of spall	Threshold spallation
No. 5	520			A few barely visible cracks	Light damage