1 INTRODUCTION
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.
2 FRACTURE MODELING
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 t_{n} can be rewritten as:
where M, C and K represent the global mass, damping and stiffness matrices respectively. The f^{int} is the vector of the internal resisting forces, and f^{ext} 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 t_{n}. 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 f_{t} 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 K_{IC} 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, E^{D} 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}).
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.
3 BLAST LOADING
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, R_{1}, B, R_{2} 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 P_{i} is the incident overpressure, P_{a} is the ambient pressure typically taken as 101kPa, P_{s} is the peak overpressure that can be obtained using various graphical and empirical methods, t_{a} is time arrival, t_{d} is the positive duration of blast loading, and b is a decay coefficient. Meanwhile, P_{r} and P_{q} refer to the reflected overpressure and wind blast respectively.
4 MODELING OF REINFORCED CONCRETE SLAB SUBJECTED TO BLAST LOADING
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.
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}).
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 f_{c} 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 G_{f} is the fracture energy, f_{t} is tensile strength, and d_{a} 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}).
Material Properties | Concrete | Steel |
---|---|---|
Young's Modulus, E (GPa) | 33 | 210 |
Poisson's Ratio, ν | 0.2 | 0.29 |
Density, ρ (kg/m^{3}) | 2400 | 7830 |
Ultimate Strength, f_{t} (MPa) | 4.0 | 560 |
Cohesion Stress, C (MPa) | 15 | - |
Friction Angle, Ø (°) | 45 | - |
Dilatation Angle, ψ (°) | 15 | - |
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 10^{5}s^{-1}, instead of many previous studies just adopted the strain rate effects up to 10^{3}s^{-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}, 10^{2} s^{-1}, 10^{3} 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 10^{4} s^{-1} to 10^{5}s^{-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.
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.
No. | Parameter | Value |
---|---|---|
1 | Discrete element contact damping (%) | 60.00 |
2 | Adaptive dynamic relaxation damping (%) | 5.00 |
3 | Discrete element contact field, l_{f} (mm) | 0.02 |
4 | Discrete element contact zone, l_{b} (mm) | 0.75 |
5 | Smallest discrete element, l_{s} (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, S_{n} (GPa) | -10 |
9 | Initial tension cut-off contact stress, S_{o} | 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.
5 RESULTS AND DISCUSSION
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.
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×300mm^{2}. Meanwhile, the spalling at the top surface of concrete slab has damage area about 160×140mm^{2}. 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.
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 S_{c} is the apparent scabbed crater diameter in meter, D_{s} 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.
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.
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}).
6 CONCLUSIONS
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.