SciELO - Scientific Electronic Library Online

vol.13 issue6Path Planning for Unmanned Underwater Vehicle in 3D Space with Obstacles Using Spline-Imperialist Competitive Algorithm and Optimal Interval Type-2 Fuzzy Logic ControllerNonlinear Dynamic Analysis of Telescopic Mechanism for Truss Structure Bridge Inspection Vehicle Under Pedestrian Excitation author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand




Related links


Latin American Journal of Solids and Structures

Print version ISSN 1679-7817On-line version ISSN 1679-7825

Lat. Am. j. solids struct. vol.13 no.6 Rio de Janeiro June 2016 


Fracture Failure of Reinforced Concrete Slabs Subjected to Blast Loading Using the Combined Finite-Discrete Element Method

Z. M. Jainia 

Y. T. Fengb 

D. R. J. Owenb 

S. N. Mokhatarc 

aJamilus Research Center, University Tun Hussein Onn, 86400 Johor, Malaysia

bZienkiewicz Centre for Computational Engineering, Swansea University, SA18EN Wales, United Kingdom,,

cDepartment of Structure and Material Engineering, University Tun Hussein Onn, 86400 Johor, Malaysia


Numerical modeling of fracture failure is challenging due to various issues in the constitutive law and the transition of continuum to discrete bodies. Therefore, this study presents the application of the combined finite-discrete element method to investigate the fracture failure of reinforced concrete slabs subjected to blast loading. In numerical modeling, the interaction of non-uniform blast loading on the concrete slab was modeled using the incorporation of the finite element method with a crack rotating approach and the discrete element method to model crack, fracture onset and its post-failures. A time varying pressure-time history based on the mapping method was adopted to define blast loading. The Mohr-Coulomb with Rankine cut-off and von-Mises criteria were applied for concrete and steel reinforcement respectively. The results of scabbing, spalling and fracture show a reliable prediction of damage and fracture.

Keywords: Fracture; the combined finite-discrete element method; reinforced concrete slab; blast loading


Blast loading from an explosion is defined as a large-scale, rapid and sudden release of energy that able to jeopardize any targets at catastrophic damage. Therefore, blast loading is known as the most destructive dynamic load that a structure may experience (Xu and Lu, 2006). If an explosion occurs in close proximately to the structure, blast loading that travel as shock wave will setup an internal compressive stress wave and subsequently produce shattering effects on the material. This leads to spalling at the surface facing the explosion and scabbing at the rear end. Thereafter, small fragments from the failure zone propagate through the thickness of the structure. For reinforced concrete structures, this normally results in a large number of fractures over a range of size and velocity, which can cause serious consequences. The fracture initiation is fully dependent on the structure form, material properties and the condition of blast loading. The fracture can be directly observed during the separation of concrete at the scabbing zone. However, spalling also produces concrete fragments due to high-intensity of blast loading particularly that happen under contact detonation.

Fracture is well known induced by the growth and coalescence of cracks. (Wang et al., 2009) stated that there is a correlation between the fracture size distribution and the crack size function. Some fraction of energy from blast loading transferred to the concrete fragments, hence the fracture travels with high velocity. Fracture velocity of disengaged elements may be up to 60m/s (Cormie et al., 2009). Currently, there is very limited available information about the fracture and its size distribution resulting from spalling and scabbing of reinforced concrete structures loaded by blast loading. A series of blast testing were carried out by (Razaqpur, 2007), (Schenker et al., 2008), (Silva and Lu, 2009), (Kim et al., 2010), (Alonso et al., 2011) and others, but these studies only examine the dynamic responses and damage behavior of concrete structures. (Wu et al., 2009) conducted an experimental study to determine the fracture size and fragment distributions of reinforced concrete slabs. In addition, (Gronsten and Forsen, 2007) conducted a field test of an explosion using a high speed camera to capture the fragment size and ejection velocity. The high energy releases and significant short duration of blast loading make the observation of fracture very unlikely and difficult to be conducted. Until recently, there is no precise method that can be employed to estimate the ascertaining of fracture size.

The phenomenal successes of the finite element method in numerical modeling of concrete structures under dynamic loadings have led to the prestigious works in the prediction of structural responses and damage mechanisms. However, this established method may unable to produce convincing results of post-failures, particularly fracture. Previous works by (Luccioni and Luege, 2006) and (Elsanadedy et al., 2011) suggested the application of the finite element method to model the local and global damage of concrete structures subjected to high loading rate. These studies were relied on the continuum failure which not really represents the mechanism of fracture and its post-failures. On the other hand, (Rabezuk and Eibl, 2006) introduced a meshless method to model the post-cracking of concrete. Although, the meshless method has proven able to produce elegant results in multi-fracturing, but still facing major problems in the dynamic responses. Since many studies on the concrete structures under blast loading concentrate on the localized responses and post-failures, hence, the continuum failure based on the stress distribution is not longer sufficient to represent such mechanisms. The evolution of the finite element method has revealed the consideration of the discrete failure that can be achieved through a special algorithm so-called the combined finite-discrete element method.

The combined finite-discrete element method has enabled simultaneous initiation and propagation of a larger number of cracks and is indeed suitable for the modeling extensive fracture processes that result in larger number of separate interacting fragments. This hybrid method is an extended version of the finite element method and the discrete element method. The combined finite-discrete element method is often needed to model the transition of continuum to discontinuous solid bodies and hence provide accurate numerical modeling of fracturing and particulate flow. Originally, the combined finite-discrete element method has been applied by (Steadm et al., 2004), (Eberhardt et al., 2004) and (Pine et al., 2006) to solve the geo-mechanics problems. In the concrete structure field, (Munjiza et al., 2000), (May et al., 2006), (Frangin et al., 2006), (Smoljanovic et al., 2013) and (Seman et al., 2013) used this method to model reinforced concrete structures under high loading rate, subsequently represent the cutting-edge in the modeling of quasi-brittle materials. Meanwhile, this paper presents numerical modeling of fracture failure of concrete structures under blast loading using the combined finite-discrete element method. Several main challenges in numerical modeling of concrete structures under blast loading using the combined finite-discrete element method are also addressed.


The combined finite-discrete element method uses a similar algorithm as in the explicit dynamic finite element method with additional parameters for the discrete element insertion. From the basic equation that governs the linear dynamic response of any structural system of finite elements, given by:


The spatially discretised form of the dynamic equilibrium equations at time tn can be rewritten as:


where M, C and K represent the global mass, damping and stiffness matrices respectively. The fint is the vector of the internal resisting forces, and fext is the velocity applied forces. The standard notations ü, and u signify the acceleration, velocity and displacement vectors in the global domain. The subscript n implies a quantity at time tn. In the context of dynamic time integration for nonlinear problems, the governing equation generalizes the stiffness dependent elastic forces to include both geometric and constitutive nonlinearities to provide internal forces with the structure (Cottrell, 2002).

In the combined finite-discrete element method, concrete as a quasi-brittle material is represented as a continuum from which cracking virtually occurs during the deformation process. The crack then propagates to all directions when the concrete strength is degraded. Before experiencing any failure, the material remains in the homogenous elastic state. The formation of crack occurs in the direction that attempt to maximize the strain energy density when the limiting tensile stress is reached, after which the material follows a softening or damaging response that governed by appropriate relation (Klerck et al., 2004). Generally, crack formation can be modeled through the adaptation of wide range models such as the fictitious crack, the cohesive crack, the fixed crack, the non-orthogonal multiple crack and smeared crack. However, many of these models were observed to be overly constrained with induced shear stress or considered overly complex and computationally expensive for the arbitrariness of the expected solution. This recovers the so-called rotating crack model which enforces coincident rotation of the principal axes of orthotropy and the principal strain axes (Bhattacharjee and Leger, 1994; Guzina et al., 1995).

However, when dealing with concrete in the combined finite-discrete element method, the commonly accepted Rankine and rotating crack models can be used to simulate crack formation within a continuum description under tensile conditions. The advantage of applying these models is that it has been extended to include rate dependent effects that governed by tensile strength and anisotropic softening. Therefore, the modeling of crack formation and fracture become very much an engineering approach. The initial failure for the Rankine and rotating crack models can be represented by:


where σi is the principal stress invariant, and ft is the tensile strength of concrete. After the initial yield, the rotating crack formulation represents the anisotropic damage evolution by degrading the elastic modulus in the direction of the major principal stress invariant (Owen and Feng, 2001; Rockfield, 2009; De Borst et al., 2012), and the stress can be evaluated using the following equation:


where is the damage parameter and is directly dependent upon the current state of strain softening within the crack band, and εn is the local strain at the local coordinate system associated with principal stress. The damage parameter is dependent on the fracture energy denoted by:


where KIC is the critical stress intensity factor, and E is the elastic modulus of concrete. After the onset of failure, developing damage degrades the elastic modulus and the damaged elastic modulus, ED is given as:


In the combined finite-discrete element method, crack formation and fracture are obtained when the topology of the mesh is updated by insertion of the discrete fracture in the failed region. A discrete fracture is introduced when the tensile strength in the principal stress direction reaches zero and is oriented orthogonal to this direction. Cracks are presumed to rotate in order to maintain this orthogonality condition upon further loading. The discrete fracture can be inserted along the failure plane (intra-element fracture), resulting in the creation of new nodes and formation of new elements in the finite element system, or along the boundaries of the existing elements (inter-element fracture) as can be seen in Figure 1. This evolution process is continued until either the system comes to equilibrium or up to the time of interest. Details about the insertion of the discrete fracture can be referred in (Cottrel et al., 2003) and (Klerck et al., 2004).

Figure 1 Crack insertion procedure, a) initial state, b) through element, or c) along element boundary (Klerck et al., 2004). 

Fracture development and potential material degradation into discrete fracture require a special treatment of mechanical forces on the contact interfaces. In the combined finite-discrete element method, the contact interaction laws are adopted in term of a penalty method. Strength parameters on the fracture interface are defined in term of cohesion and friction angle based on a linear Mohr-Coulomb criterion (Vyazmensky et al., 2007). Following the insertion of a discrete fracture, the crack surface assumes a frictional contact response governed by contact parameters based on the theory of the contact surface conditions. The contact parameters are introduced to control the element separation from finite element to discrete element and in assuring the maintenance of energy balance during this transition period.


When an explosive detonates, the blast wave is initiated by the rapid release of a large amount of energy. The blast wave expanding spherically and propagates under compressed air within a fraction of second. The time of blast wave moves outward from the detonation source to the rigid surface is known as the time of arrival. The blast wave then loaded on structure as an incident overpressure that rises abruptly to a peak value and decays exponentially to the ambient pressure within the short duration of positive phase. The overpressure continues to decline and drop below the ambient pressure so-called the negative phase of blast wave. After acting on the surface of structure, the incident overpressure is transformed as a reflected overpressure and travel back through air that has been heated and compressed by the passage of preceding blast wave. Two distinct of the incident and reflected overpressures that act simultaneously are known as the static overpressure. Besides that, the blast wave is also accompanied by the dynamic pressure known as wind blast.

In numerical modeling, the Jones-Wilkins-Lee equation of state is widely used to represent blast wave. Blast wave is transferred to the surface of structure using coupling of Arbitrary Lagrangian-Eulerian model for multi-material fluids, and the Lagrangian model for structure (Souli et al., 2000; Bouamoul and Dang, 2008; Chafi et al., 2009; Deng and Jin, 2009). The Jones-Wilkins-Lee equation of state defines the pressure of explosion as:


where p is pressure, e is internal energy, V is the specific volume, while A, R1, B, R2 and ω are all material constants than can be referred in (Dobartz, 1981) and (Zhou et al., 2008). The Jones-Wilkins-Lee is generally considered as a pure empirical equation of state with a non-physical coefficient considered to be constant. As a result, the exponential terms can give unrealistic behaviour at high pressures. The Jones-Wilkins-Lee equation of state is normally used for modeling the pressure in application involving elastic models, isochoric plasticity and visco-plasticity models. Therefore, (Jaini and Feng, 2011) found that this equation of state cannot be used with compressible models like Rankine, rotating crack, Mohr-Coulomb and Drucker-Prager models.

Alternatively, the mapping method is the most convenient approach to model the blast loading. The determination of blast loading involves the simultaneous calculation of the incident overpressure, the reflected overpressure and wind blast. All parameters are primarily dependent on the amount of energy released by an explosion in the form of blast wave and standoff distance. The incident overpressure occurs after the blast wave hits the structure. It is basically determined as the pressure-time history using Friedlander's formula as shown in Equation (8). The reflected overpressure can be determined based on Rankine's formula as depicted in Equation (9), while (Ngo et al., 2007) suggested Equation (10) to determine wind blast.




where Pi is the incident overpressure, Pa is the ambient pressure typically taken as 101kPa, Ps is the peak overpressure that can be obtained using various graphical and empirical methods, ta is time arrival, td is the positive duration of blast loading, and b is a decay coefficient. Meanwhile, Pr and Pq refer to the reflected overpressure and wind blast respectively.


The model of reinforced concrete slab was based on an experimental study by (Zhou et al., 2008). The concrete slab has dimension of 1300mm length, 1000mm width and 100mm thickness. There is respectively 7 and 13 one-way longitudinal steel reinforcements with a diameter of 16mm located at the top and bottom of the concrete slab. Figure 2 shows the experimental setup and details of concrete slab. The concrete slab was detonated with explosive charge of the Composition B, where the weight of explosion is 0.5kg that equivalent to 0.545kg of TNT. The explosive charge was located 100mm above the top surface of concrete slab. In numerical modeling, the concrete slab was modeled in three-dimensional to ensure the accuracy of non-uniform load distribution and to avoid the effects of boundary surface especially the propagation of multi-reflected waves.

Figure 2 Experimental study of concrete slab under blast loading (Zhou et al., 2008). 

The concrete slab was divided into three volumes and discretized with unstructured 4-noded solid tetrahedral element. Figure 3 shows the finite-discrete element model which is composed of 42481 meshes. The tetrahedral element was used since it can be generated automatically using the advancing front or Delaunay mesh generation technique, geometrically versatile and fairly represent a complex three-dimensional domain. Moreover, this type of element is desirable from the consideration of the number of nodes in conjunction with the number of degree-of-freedom that have the direct impact on the solution time, especially in the modeling that involve crack formation under high loading rate and characterized by extremely high strains (Timothy et al., 1991; Owen et al., 2007).

Figure 3 Finite-discrete element model of reinforced concrete slab. 

Meanwhile, the steel reinforcements were located along the boundaries of the tetrahedral between the volumes of concrete slab. The steel reinforcements were spatially meshed as 2-noded simo beam consist of 1302 elements. The beam elements were also formulated to present bending deformation with the definition of the material stiffness. Material stiffness is represented by the flexural and shear rigidities in term of the elastic material characteristics and geometric properties such as the moment of area, the cross-sectional area and the effective shear area. On the other hand, the local effects such as gradual bond deterioration, dowel action and pullout effects were represented indirectly by the degradation of the surrounding continuum. The perfect bond was assumed between the steel reinforcements and concrete material. This assumption of compatibility of displacement and strains between steel and concrete allows the reinforcement to be treated as an integral part of the three-dimensional finite-discrete element model (Cerveta and Hinton, 1986).

4.1 Material Models

In order to achieve both continuum and fracture failures, two types of plasticity models namely the non-linear material model and brittle-fracture model were employed. The non-linear material model involves the isotropic Mohr-Coulomb with Rankine cut-off for concrete and von-Mises for steel reinforcement. It should be noted here that the limitation of the isotropic Mohr-Coulomb with Rankine cut-off is that the deformation analysis requires a flow rule to determine the orientation of the strain-increment vector with respect to the yield condition (Labuz and Zang, 2012). The compression failure is a notoriously difficult one to be obtained as the non-linear material model has no cap on its shear strength and normal stress. Meanwhile, the brittle fracture model was defined in the concrete material using the rotating crack model with rate dependency. The use of non-linear material and brittle fracture models are based on a traditional Lagrangian approach that has been enhanced by a number of numerical procedures to enable effective simulation of a combination of large strain and high pressure in quasi-brittle materials.

The material properties for concrete and steel reinforcement can be referred in Table 1. The compressive strength, tensile strength and density of concrete and steel were followed the data that provided by (Zhou et al., 2008). Meanwhile, the Young's modulus of concrete was simply calculated using the following equation:


where E is the Young's modulus, and fc is compressive strength of concrete taken as 50MPa. The fracture energy of 100N/m was defined on concrete. This value of fracture energy was calculated using an equation proposed by (Oh et al., 1999):


where Gf is the fracture energy, ft is tensile strength, and da is the size of aggregate considered as 15mm to 20mm. Previous experimental studies by (Ulfkjaer and Brincker, 1995), (Rao and Prasad, 2001), (Bazant and Becq-Giraudon, 2001), and (Amir et al., 2014) showed that the fracture energy of concrete with compression strength of 50MPa is around 90N/m to 120N/m. On the other hand, the hardening properties of concrete, as can be seen in Table 2, were envisaged to follow the non-associated flow rule. These data of cohesion stress, friction angle and dilatation angle were adopted based on the suggestion by (May et al., 2006) and (Bere, 2007).

Table 1 Material properties of concrete and steel reinforcement. 

Material Properties Concrete Steel
Young's Modulus, E (GPa) 33 210
Poisson's Ratio, ν 0.2 0.29
Density, ρ (kg/m3) 2400 7830
Ultimate Strength, ft (MPa) 4.0 560
Cohesion Stress, C (MPa) 15 -
Friction Angle, Ø (°) 45 -
Dilatation Angle, ψ (°) 15 -

Table 2 Hardening properties of concrete. 

Strain Cohesion Stress, C (MPa) Friction Angle, Ø (°) Dilatation Angle, ψ (°)
0.0 15 45 15
0.03 15 45 5
1 15 45 0

In addition to the required data in the fracture failure, the rate dependency of concrete as can be seen in Figure 4 was also introduced through the use of the dynamic tensile strength with respect to the strain rate effects. The dynamic increase factor was used to increase the tensile strength of concrete at zero strain rate until 105s-1, instead of many previous studies just adopted the strain rate effects up to 103s-1. Data of the dynamic increase factor for the tensile strength of concrete was obtained from various experimental and numerical based literatures, as example (Malvar and Ross, 1998), (Lambert and Allen, 2000), (Brara et al., 2001) and (Schuler et al., 2006). It was identified that the dynamic increase factor for strain rate of 1s-1, 5s-1, 20s-1, 50s-1, 102 s-1, 103 s-1 is 1.3, 2.5, 4.5, 8.5, 12.5 and 16.0 respectively, and remain constant at 17.0 for strain rate within the range of 104 s-1 to 105s-1. The dynamic tensile strength can be obtained by multiplying the tensile strength with dynamic increase factor at respective strain rate. For the steel reinforcement, it is acknowledged that the strain rate effects does result in an increase of the yield stress of steel but the effect is known to be relatively minor and therefore is not considered in numerical modelling. (Malvar and Crawford, 1998) described that the dynamic increase factor for high yield plastic strain at 0.005 and 0.022 is around 1.0 to 1.1 only. By assuming no strain rate effects for the steel reinforcement, the numerical modeling can be assumed to be conservative.

Figure 4 Strain rate dependency of concrete. 

4.2 Contact Surface and Discrete Element Data

The contact interaction of discrete bodies was based on the edge-to-edge algorithm. The contact resolution was provided through the linear penalty method as suggested by (Yu, 1999). An additional contact damping is also included through the use of velocity/momentum-dependent damping and global damping, set at 60% and 5% respectively. Others global discrete element contact data are listed in Table 3. These contact surface and discrete element data were initially followed the recommendation of (Bere, 2004) and (Rockfield, 2009). In general, the value of discrete element contact field is normally 10% to 20% of the minimum side length of the finite elements in the mesh system. Field sizes that are too small may result in penetration of a node through an element. The discrete contact zone is the size of buffer zone for contact detection. This is used to locate nodes close to any surface during contact detection and the ideal size should be equal to the average size of the side lengths of the finite element mesh.

Table 3 Global discrete element contact parameters for concrete. 

No. Parameter Value
1 Discrete element contact damping (%) 60.00
2 Adaptive dynamic relaxation damping (%) 5.00
3 Discrete element contact field, lf (mm) 0.02
4 Discrete element contact zone, lb (mm) 0.75
5 Smallest discrete element, ls (mm) 0.00
6 Global normal penalty coefficient, α n (GPa) 300.00
7 Global tangential penalty coefficient, α t (GPa) 30.00
8 Normal contact stress, Sn (GPa) -10
9 Initial tension cut-off contact stress, So 0
10 Coulomb friction coefficient, μ 0.2

Meanwhile, the smallest discrete element is a compulsory data for fracturing model. If the element size becomes smaller than this element, then the fracture will only occur between the element boundaries and not through the element. By setup the smallest discrete element equal to zero, the fracture can be obtained without relying on the specific size. In the dynamic-transient problems, the value of normal penalty coefficient is typically in the range of 0.5 to 2.0 times of Young's modulus. The tangential penalty is normally taken as the friction coefficient multiple with Young's modulus or usually taken as an order of magnitude approximately 0.1 of the normal penalty coefficient. It should be noted that the contact penalty coefficients applied in the concrete slab are considerably higher than what have been used for typical dynamic-transient simulations to date due to far higher pressure being generated from the blast loading. A high contact penalty stiffness was used in the concrete slab to ensure that the correct amount of energy is transmitted from the high concentration of blast pressure to the concrete material and the energy can be dissipated correctly to the steel reinforcements and support materials. This is also to ensure smooth formation of superposition of the shock wave and multi-reflected waves.

4.3 Blast Pressure-Time History

Blast loading was modeled using the multiple pressure-time history in space and time-dependent function. Therefore, the blast loading behaves non-uniformly across the surface of concrete slab. In numerical modeling, all overpressures are assumed to act simultaneously on the surface of concrete slab. As a result, the blast loading was modeled as a cumulative overpressure of the ambient pressure, incident overpressure, reflected overpressure and wind blast. This is a practical approach in the air blast simulation even though wind blast moves behind the incident overpressure. The reflected overpressure and wind blast move forward in the spherical three-direction whilst the ambient and incident overpressures are only propagate perpendicular inward the surface of concrete slab. Since the magnitude of the ambient and incident overpressures are very small compared to the reflected overpressure and wind blast, all overpressures for a single pressure-time history can be considered moving in the same direction within a similar duration.

Since the Jones-Wilkins-Lee equation of state is incompatible with the compressible material models, the mapping method was incorporated to define the blast loading on the surface of concrete slab. Although the mapping method consumes much longer preparation time, it is capable of saving computational run time and yet produces accurate results. The surface of concrete slab was mapped onto fifteen areas (area a until n) according the arrival time and the positive duration of blast loading. Beyond the area n, only ambient pressure of 101kPa was considered. For a comparison, an empirical formula proposed by (Urgessa, 2009) was also adopted to verify the area that should be imposed by blast loading. Figure 5 shows the blast pressure-time history. Meanwhile, the mapped area can be seen in Figure 6.

Figure 5 Blast pressure time-history for various mapped areas. 

Figure 6 Mapping area of blast on the surface of concrete slab. 


The numerical modeling using the combined finite-discrete element method allows comprehensive evaluations of the results in both continuum and fracture failures. When assessing the structure under high loading rates particularly blast loading, numerical results in term of scabbing, spalling and fracture behavior are of prime important. When an explosion occurs near to the surface of concrete slab, the blast pressure or shock wave is first applied very locally at the point on the concrete slab closest to the explosive charge and then vary with the distances across the surface of concrete slab. As a result, the damage of concrete slab can be separated into local and overall responses which occur at different time periods (Duranovic and Watso, 1994). The transition phase of effective strain rate from high to low at 10s-1 can be used to distinguish these two forms of predominant damage (Zhou et al., 2008). It is found from Figure 7 that the concrete slab started experiences low effective strain rate at time 0.5msec. This effective strain rate was evaluated based on the strain rate components varies with time that occurring during the crack formation and fracturing process at the center of beneath zone of concrete slab.

Figure 7 Strain rate-time history. 

5.1 Spalling and Scabbing

During the interaction of blast pressure on the concrete slab, which is happen approximately within 0.27msec at relatively high strain rate, the damage is localized as characterized by scabbing, spalling and shear plug. When the time reaches 0.5msec, the scabbing fully develops in a crater with maximum depth of 50mm and damage area approximately 480×300mm2. Meanwhile, the spalling at the top surface of concrete slab has damage area about 160×140mm2. The displacement of the concrete slab is 8.82mm. A comparison between numerical and experimental results at time 0.5msec shows a favorable agreement and gives a reliable prediction on the scabbing and spalling as depicted in Figure 8. Further damage and progressive failure occur under the multi-reflected waves and continuously exist due to the deformation owing by dynamic vibrations. This can be clearly observed by the fluctuation of low strain rate.

Spalling and scabbing damage are always formed during the duration of shock wave and reflected wave. During the high strain rate stage, the scabbed crater is formed from the scabbing damage and the propagation of cracks along the steel reinforcement bars. The diameter of scabbed crater increases under the reflected wave and remains constant during the dynamic vibration. It is found that the dynamic vibration has a little effect on the scabbing damage. However, due to the combination of reflected wave and dynamic vibration, approximately at time of 0.5msec to 0.75msec, cracks are propagate from the epicenter of concrete slab towards the free edges. Unlike scabbing damage which occurs soon after the concrete slab was loaded by blast pressure, spalling damage happens only during the transition of high to low strain rate. After that, the spalling damage propagates quickly but has no contribution to the fracture and fragmentation of the concrete slab. The diameter of the spalling and scabbing damage in the time history is shown in Figure 9. It can be clearly observed that the spalling and scabbing damage are almost constant under the dynamic vibration of more than 0.75msec.

Figure 8 Spalling and scabbing damage from experimental and numerical results. 

Figure 9 Spalling and scabbing damage. 

Consider the initial of spalling damage and scabbed crater of 0.5kg explosive at 100mm and the scaling law established that any linear dimension of the crater can be expressed as a constant multiplied by weight of the explosion (Ambrosini and Luccioni, 2012), simple formulas to evaluate the diameter of the spalling and scabbing damage can be written as:



where Sc is the apparent scabbed crater diameter in meter, Ds is the spalling damage diameter that also in meter, and W is the equivalent TNT explosive weight in kilogram. Both equations have a similar form of crater dimension defined by (Kinney and Graham, 1985) and (Luccioni et al., 2009) in the prediction of damage and crater that caused by a spherical blast wave. However, these equations are considered conservative as damage and crater behavior were defined in term of weight of explosive, whereas the damage behavior is highly likely to be influenced by the standoff distance and incidence angle in the form of scaled-distance function. However, when Equation (13) is plotted against the weight of explosion as shown in Figure 10, the results of scabbed crater diameter are in the middle of the maximum and minimum curves produced by (Xu and Lu, 2006). The experimental data from (Yuan et al., 2008) has lower scabbed crater diameter at weight of explosion of 5kg, 8kg and 10kg. This is probably because of the ratio of span to thickness (L/d) of the concrete slab used in experimental study is smaller than 10. It should be noted here that the L/d of the concrete slab in numerical modeling is 13.

Figure 10 Prediction of scabbed crater diameter of reinforced concrete slab. 

5.2 Fracture Failure

Shear plug and fracture onset are two failure mechanisms that cannot be captured directly through the experimental work. However, the numerical modeling using the combined finite-discrete element method has successfully produced these failure mechanisms. Figure 11 shows the formation of shear plug and fracture that happen during the explosion phase and the dynamic vibration at time instants ranging from 0.5msec to 5.0msec. Red color indicates the continuum body while blue color represents the fractured elements. The shear plug and fracture were observed along the mid-cross section of long span direction. At the time of 0.5msec during the transition from high to low strain rate, the initial region of potential fracturing is observed to be appeared at the mid-layer of the concrete slab and some fracture occurs at the beneath zone.

Figure 11 Shear plug and fracture formation at time 0.5 msec to 5.0 msec. 

A shear plug in the pattern of a conical frustum is formed by incline cracking through the thickness of concrete slab. This phenomenon corresponds closely to the formation of shear plug investigated by (Yankelevsky, 1994) on the concrete slab under low velocity impact. Then, fracture aggressively forms along the concrete thickness as new cracks are generated at the middle of span. It is found that the shear plug significantly contributes to the separation of elements and lead to the formation of scabbed fracture. In addition, fracture also happen due to the cracks of the bond-slip phenomenon between steel reinforcements and concrete. In this numerical modeling, there is no fracture that induces from spalling damage. However, (Yongxiang et al., 2009) argued that a concrete structure under contact explosion may experience serious spalled fracture.

Currently, the size of the fractured elements cannot be directly quantified. Unlike scabbing and spalling, the concrete fragments from numerical modeling are not based on the actual condition. The termination system designed in this numerical modeling is based on the maximum steady state energy, thus when the energy being reached, not all scabbed particles will be significantly departed from the concrete slab. In addition, there are still no relevant sources that can be taken as references in predicting the size of fracture and fragment weight of concrete structure under blast loading. Although (Wu et al., 2009) proposed an equation to estimate the size of fracture, but this will still be complex as parameter of fragment size depends on the size of mesh. In fracture failure, it is difficult to define a universally acceptable threshold since fracture depends on many factors and is subjective to the ability of elements to be separated. In this study, an effort to investigate the fragment weight of concrete slab was performed by assessing fracture together with scabbed damage. The grievance of this approach is that the fragment weight of concrete may unjustified and overestimated, but such approach is extremely valuable. Figure 12 displays the relation of the velocity with weight of fractured concrete. When the fracture ejected from the concrete slab, it is moving down with high velocity approximately 58m/sec, similar as suggested by (Cormie et al., 2009).

Figure 12 Mass and velocity of fractured concrete. 


This paper presents the fracture failure of reinforced concrete slabs subjected to blast loading. The fracture failure was achieved by employing the combined finite-discrete element method with designated material models. The material model of Mohr-Coulomb with Rankine cut-off was defined to concrete whilst steel reinforcement uses the von-Mises criterion. In addition, the rotating crack material model was used to simulate crack formation and fracture within a continuum description. The concrete slab was imposed with multiple blast pressure-time history using a mapping method. The blast loading is non-uniform across the surface of concrete slab and formed simultaneously by the combination of incident overpressure, the reflected overpressure and wind blast. It is demonstrated that the combined finite-discrete element method with accurate definition of material model, discrete element contact conditions and blast loading shows a high confidence in the numerical prediction of damage mechanisms and post-failures. The comparison of the numerical results with that obtained from experimental study shows a reliable prediction of spalling and scabbing. It is also found that fracture occurs due to rapid ejection of shear plug that contribute to the separation of elements. After experiences complete separation, the fracture travels with high velocity. However, the fracture failure in term of its size and distribution are still characterized by an enormous amount of uncertainty and the conjecture unconfirmed. Although an attempt was made to determine the weight of fractured concrete directly from the scabbed damage elements, the shortcoming of this approach is that the weight of fracture may be unjustified and overestimated. However, this approach is extremely valuable and the accuracy of results can be considered convincing.


Alonso, G.M., Cendon, D.A., Galvez, F., Erice, B. and Galvez, V.S. (2011). Blast response analysis of reinforced concrete slabs: experimental procedure and numerical simulation, Journal of Applied Mechanics 78(5): 510-522. [ Links ]

Ambrosini, R.D. and Luccioni, B.M. (2012). Craters produced by explosions on, above and under the ground, Advances in protective structures research, Taylor and Francis (London), 365-396. [ Links ]

Amir. S., van der Veen, C. and Walraven, J.C. (2014). Numerical investigation of the bearing capacity of transversely prestressed concrete deck slabs, Computational modelling of concrete structures, Taylor and Francis (London), 99-108. [ Links ]

Bazant, Z.P. and Becq-Giraudon, E. (2002). Statistical prediction of fracture parameters of concrete and implications for choice of testing standard, Cement and Concrete Research 32(4): 529-556. [ Links ]

Bere, A.T. (2004). Computational modelling of large-scale reinforced concrete structures subjected to dynamic loading, Ph.D. Thesis, Swansea University, United Kingdom. [ Links ]

Bhattacharjee, S.S. and Leger, P. (1994). Application of NLFM models predict crcaking in concrete gravity dams, Journal of Structural Engineering 120(4): 1255-1271. [ Links ]

Bouamoul, A. and Dang, T.V.N. (2008). High explosive simulation using arbitrary lagrangian-eulerian formulation, Technical Report TM 2008-254, Defence Research and Development, Canada. [ Links ]

Brara, A., Camborde, F., Klepaczko, J.R. and Mariotti, C. (2001). Experimental and numerical study of concrete at high strain rates in tension, Mechanics of Materials 33(1): 33-45. [ Links ]

Cerveta, M. and Hinton, E. (1986). Computational modelling of reinforced concrete structures, Pineridge Press (Swansea). [ Links ]

Chafi, M.S., Karami, G. and Ziejewski, M. (2009). Numerical analysis of blast-induced wave propagation using FSI and ALE multi-material formulation, International Journal of Impact Engineering 36(10-11): 1269-1275. [ Links ]

Cormie, D., Mays, G. and Smiths, P. (2009). Blast effects on buildings, Thomas Telford (London). [ Links ]

Cottrell, M.G. (2002). The development of rational computational strategies for the numerical modelling of high velocity impact, Ph.D. Thesis, Swansea University, United Kingdom. [ Links ]

Cottrell, M.G., Yu, J., Wei, Z.J. and Owen, D.R.J. (2003). The numerical modelling of ceramics subjected to impact using adaptive discrete element technique, Engineering Computations 20(1): 82-106. [ Links ]

De Borst, R., Crisfield, M.A., Remmers, J.J.C. and Verhoosel, C.V. (2012). Non-linear finite element analysis of solids and sructures, John Wiley and Sons (West Sussex) [ Links ]

Deng, R.B. and Jin, X.L. (2009). Numerical simulation of bridge damage under blast loads. WSEAS Transaction in Computers 9(8): 1564-1574. [ Links ]

Dobartz, B.M. (1981). LLNL explosives handbook: properties of chemical explosives and explosive simulation, University of California Press (California). [ Links ]

Duranovic, N. and Watso, A.J. (1994). Impulsive loading on reinforced concrete slabs-blast loading functions, Transaction on the Built Environment 8: 63-73. [ Links ]

Eberhardt, E., Steadb, D. and Coggan, J.S. (2004). Numerical analysis of initiation and progressive failure in natural rock slopes, International Journal of Rock Mechanics and Mining Sciences 41(1): 69-87. [ Links ]

Elsanadedy, H.M., Almusallam, T.H., Abbas, H., Al-Salloum, Y.A. and Alsayed, S.H. (2011). Effect of blast loading on cfrp-retrofitted rc columns-a numerical study, Latin American Journal of Solids and Structures 8(1): 55-81. [ Links ]

Frangin, E., Marin, P. and Daudeville, L. (2006). On the use of combined finite/discrete element method for impacted concrete structures, Journal de Physique IV 134(1): 461-466. [ Links ]

Gronsten, G.A. and Forsen, R. (2007). Debris launch velocity from overloaded concrete cubicles, Proceeding of the International Symposium on Interaction of the Effects of Munitions with Structures, Orlando. [ Links ]

Guzina, B.B., Rizzi, E., Willam, K. and Pak, R.Y.S. (1995). Failure prediction of smeared-crack formulations, Journal of Engineering Mechanics 121(1): 150-161. [ Links ]

Jaini, Z.M. and Feng, Y.T. (2011). Computational modelling of damage behaviour of reinforced concrete slabs subjected to blast loading, Proceeding of the 8th International Conference on Structural Dynamics, Leuven. [ Links ]

Kim, J.H.J., Yi, N.H., Oh, I.S., Lee, H.S., Choi, J.K. and Cho, Y.G. (2010). Blast loading response of ultra performance concrete and reactive powder concrete slabs, Fracture mechanics of concrete and concrete structures, Korean Concrete Institute (Seoul), 1715-1722. [ Links ]

Kinney, G.F. and Graham, K.J. (1985). Explosive shocks in air, Springer-Verlag (New York). [ Links ]

Klerck, P.A., Sellers, E.J. and Owen, D.R.J. (2004). Discrete fracture in quasi-brittle materials under compression and tensile stress states, Computer Methods in Applied Mechanics and Engineering 193(27-29): 3035-3056. [ Links ]

Labuz, J.F. and Zang, A. (2012). Mohr-Columb failure criterion, Rock Mechanics and Rock Engineering 45(6): 975-979. [ Links ]

Lambert, D.E. and Allen, R.C. (2000). Strain rate effects on dynamic fracture and strength, International Journal of Impact Engineering 24(10): 985-998. [ Links ]

Luccioni, B.M. and Luege, M. (2006). Concrete pavement slab under blast loads, International Journal of Impact Engineering 32(8): 1248-1266. [ Links ]

Luccioni, B.M., Ambrosini, D., Yuen, S.C.K. and Nurick, G.N. (2009). Evaluating the effect of large and spread explosives loads, Argentina Association for Computational Mechanics 28: 529-552. [ Links ]

Malvar, L.J. and Crawford, J.E. (1998). Dynamic increase factors for steel reinforcing bars, Proceeding of the 28th DDESB Seminar, Orlando. [ Links ]

Malvar, L.J. and Ross, C.A. (1998). Review of strain rate effects for concrete in tension, ACI Materials Journal 95(6): 735-739. [ Links ]

May, I.M., Chen, Y., Owen, D.R.J., Feng, Y.T. and Thiele, P.J. (2006). Reinforced concrete beams under drop-weight impact loads, Computers and Concrete 3(2-3): 79-90. [ Links ]

Munjiza, A., Latham, J.P. and Adrews, K.R.F. (2000). Detonation gas model for combined finite-discrete element simulation of fracture and fragmentation, International Journal of Numerical Methods in Engineering 49(12): 1495-1520. [ Links ]

Ngo, T., Mendis, P., Gupta, A. and Ramsay, J. (2007). Blast loading and blast effects on structures: an overview, Electronic Journal of Structural Engineering Special Issue: 76-91. [ Links ]

Oh, B.H., Jang, S.Y. and Byun, H.K. (1999). Prediction of fracture energy of concrete, KCI Concrete Journal 11(3): 211-221. [ Links ]

Owen, D.R.J. and Feng Y.T. (2001). Parallelised finite/discrete simulation of multi-fracturing solids and discrete systems, Engineering Computations 18(3/4): 557-576. [ Links ]

Owen, D.R.J., Feng, Y.T., Cottrell, M.G. and Yu, J. (2007). Computational issues in the simulation of blast and impact problems: an industrial perspective, Springer (Netherlands), 3-35. [ Links ]

Pine, R.J., Coggan, J.S., Flynn, Z.N. and Elmo, D. (2006). The development of a new numerical modelling approach for naturally fractured rock masses, Rock Mechanics and Rock Engineering 39(5): 395-419. [ Links ]

Rabezuk, T. and Eibl, J. (2006). Modelling dynamic failure of concrete with meshfree methods, International Journal of Impact Engineering 32(11): 1878-1897. [ Links ]

Rao, G.A. and Prasad, B.K.R. (2002). Fracture energy and softening behavior of high-strength concrete, Cement and Concrete Research 32: 247-252. [ Links ]

Razaqpur, A.G., Tolba, A. and Contestabile, E. (2007). Blast loading response of reinforced concrete panels reinforced with externally bonded gfrp laminates, Composites Part B 38: 535-546. [ Links ]

Rockfield Software Limited. (2009). Elfen explicit manual version 4.4: data file specification-material models, Rockfield Software Limited (Swansea). [ Links ]

Schenker, A., Anteby, I., Gal, E., Kivity, Y., Nizri, E., Sadot, O., Michaelis, R., Levintant, O. and Ben-Dor, G. (2008). Full-scale field tests of concrete slabs subjected to blast loads, International Journal of Impact Engineering 35: 184-198. [ Links ]

Schuler, H., Mayrhofer, C. and Thoma, K. (2006). Spall experiments for the measurement of the tensile strength and fracture energy of concrete at high strain rates, International Journal of Impact Engineering 32(10): 1635-50. [ Links ]

Seman, M.A., Nasly, M.A., Feng, Y.T. and Jaini, Z.M. (2013). Computational modelling of reinforced concrete wall subjected to transformer tank rupture, Research and applications in structural engineering, mechanics and computation, Taylor and Francis (London), 195-196. [ Links ]

Silva, P. and Lu, B. (2009). Blast resistance of reinforced concrete slabs, Journal of Structural Engineering 135(6): 708-716. [ Links ]

Smoljanovic, H., Zicaljic, N. and Nikolic, Z. (2013). Nonlinear analysis of engineering structures by combined finite-discrete element method, Gradevinar 65(4): 331-344. [ Links ]

Souli, M., Ouahsine, A. and Lewin, L. (2000). ALE formulation for FSI problems, Computer Methods in Applied Mechanics and Engineering 190(5-7): 659-675. [ Links ]

Steadm, D., Coggan, J.S. Eberhardt, E. (2004). Realistic simulation of rock slope failure mechanisms: the need to incorporate principles of fracture mechanics, International Journal of Rock Mechanics and Mining Sciences 41(1): 563-568. [ Links ]

Timothy, P.P., Shah, M.Y. and Robert D.C. (1991). Solid elements with rotational degrees of freedom: part ii-tehrahedron elements, International Journal of Numerical Methods in Engineering 31(3): 593-610. [ Links ]

Ulfkjaer, J.P. and Brincker, R. (1995). Fracture energy of normal strength concrete, high strength concrete and ultra high strength ultra ductile steel fibre reinforced concrete, Aedificatio Publishers (Freiburg), 31-44. [ Links ]

Urgessa, G.S. (2009). Finite element analysis of composite hardened walls subjected to blast loads, Journal of Engineering and Applied Sciences 2(4): 804-811. [ Links ]

Vyazmensky, A., Elmo. D. and Stead, D. (2007). Combined finite-discrete element modelling of surface subsidence, Proceeding of the 1st Rock Mechanics Symposium, Vancouver. [ Links ]

Wang, M., Hao, H, Ding, Y. and Li, Z.X. (2009). Prediction of fragment size and ejection distance of masonry wall under blast loading using homogenized masonry material properties, International Journal of Impact Engineering 36(6): 808-820. [ Links ]

Wu, C., Nurwidayati, R. and Oehlers, D.J. (2009). Fragmentation from spallation of rc slabs due to airblast loads, International Journal of Impact Engineering 36(12): 1371-1376. [ Links ]

Xu, K. and Lu, Y. (2006). Numerical simulation study of spallation in reinforced concrete plates subjected to blast loading, Computer and Structure 84(5-6): 431-438. [ Links ]

Yankelevsky, D.Z. (1994). Local response of concrete slab to low velocity missile impact, International Journal of Impact Engineering 19(4): 331-343. [ Links ]

Yongxiang, D., Shunshan, F., Changjing, X. and Lele, G. (2009). Dynamic behaviour of concrete sandwich panel under blast loading, Defence Science Journal 59(1): 22-29. [ Links ]

Yu. J. (1999). Contact interaction framework for numerical simulation of multi-body problems and aspects of damage and fracture for brittle material, Ph.D. Thesis, Swansea University, United Kingdom. [ Links ]

Yuan, L., Gong, S. and Jin, W. (2008). Spallation mechancism of rc slabs under contact detonation, Transaction Tianjin University 14(6): 464-469. [ Links ]

Zhou, X.Q., Kuznetsov, V.A., Hao, H. and Wachl, J. (2008). Numerical prediction of reinforced concrete response to blast loading, International Journal of Impact Engineering 35(10): 1186-1200. [ Links ]

Received: August 11, 2015; Revised: January 07, 2016; Accepted: February 19, 2016

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License