An effective model for analysis of reinforced concrete members and structures under blast loading

Abstract

The traditional fiber beam model has been widely used in the seismic analysis of reinforced concrete members and structures. However, the inability to capture shear failure restricts its application to blast loadings. In this article, a numerical model that considers both rate-dependent shear behavior and damage effect is proposed based on the traditional fiber beam element. This is achieved using the modified compression-field theory with a concrete damage model and bilinear steel model in the principal directions. Meanwhile, a condensed three-dimensional stress–strain relation from the isotropic hardening plasticity model is implemented to simulate longitudinal reinforcement bars, as large shear strain would be produced under severe blast loads. The proposed model is validated by comparing the numerical and test results. The high-fidelity physics-based finite element model, validated by the same experiment, is also used in the study to prove the efficiency of the proposed model. Case studies of a reinforced concrete beam and a six-story reinforced concrete frame structure subjected to blast loads are then carried out. The results indicate that the proposed model is reliable compared with the high-fidelity physics-based model. In addition to the accuracy, comparisons of the computational time show an excellent performance with respect to the efficiency of the proposed model.

Keywords

blast loading diagonal shear failure fiber beam model numerical simulation reinforced concrete

Introduction

With the increasing terrorist attacks and accidental explosions, which may lead to a large number of casualties and great economic losses, an increasing number of researchers have turned their attention to the response analysis of reinforced concrete (RC) members and structures under blast loading to provide guidance for protective designs (Crawford and Magallanes, 2011; Krauthammer, 2008, 2014; Li and Shi, 2015; Williams and Williamson, 2011). When an RC member is subjected to an intense transverse dynamic load, for example, blast or impact, the propagation of elastic and plastic waves along the transverse direction may lead to a local failure by spalling, as illustrated in Figure 1. For beam/column members, their transverse dimension (thickness) is usually much smaller than the main dimension (length) such that it takes much less time for a stress wave to propagate through the former dimension of the member, typically a few microseconds. The transverse stress wave propagates back and forth several times and dramatically reduces to zero, followed by the movement of each cross section as an entirety, generating dynamic deformation along the element. This relatively longer-term dynamic behavior is called global structural response and usually takes several milliseconds or much longer depending on the structural properties. Under blast loading, the failure modes of columns generally include breach, spalling, direct shear, diagonal shear, flexure, or combined patterns (Crawford and Magallanes, 2011). Among these failure modes, diagonal shear failure is very common, as can be observed in some blast experiments by Karagozian & Case (K&C) (Crawford et al., 2012). The brittle failure of RC columns can cause catastrophic collapse of structures, because they are critical load-bearing components, resulting in a serious threat to the safety of structures, personnel, and properties (Krauthammer, 2014).

Figure 1.

Responses of RC columns under blast loading.

Currently available blast analysis methods for RC components are typically based on the single-degree-of-freedom (SDOF) system (Biggs, 1964), Euler–Bernoulli beam theory (EBT) (Carta and Stochino, 2013), Timoshenko beam theory (TBT) (Krauthammer et al., 1993a, 1993b, 1994), and high-fidelity physics-based (HFPB) finite element (FE) methods (Bao and Li, 2010; Crawford et al., 2012; Crawford and Magallanes, 2011; Hao and Hao, 2014; Magnusson et al., 2010; Shi et al., 2008). Although the SDOF method has the capability to account for the combined effects of flexure, shear, and axial force using a specifically derived resistance function, its overly simplified assumption generally leads to unreliable and inaccurate results (El-Dakhakhni et al., 2009). The EBT has been used for the study of RC beams that failed in flexure under blast loads (Carta and Stochino, 2013), while this theory is inadequate for impact type loadings because its governing equations are dispersive: both the phase velocity and the group velocity are proportional to the wave number, which means that the velocities become infinite as the wavelength approaches zero (Flügge and Zajac, 1959). However, experimental observations indicate that the transverse motions of a beam have a finite maximum wave speed irrespective of the wavelength (Jones, 2011). As for the TBT, it is able to account for not only rotatory inertia but also shear deformation (Timoshenko, 1921); its governing equations can be solved using the finite difference (Krauthammer et al., 1993a, 1993b, 1994) or FE (Dragos and Wu, 2014; Wu and Sheikh, 2013) formulation, but the latter is more straightforward for a subsequent whole structural analysis when it is implemented in a general FE code, such as LS-DYNA (Hallquist, 2007). Due to the high cost of blast experimental programs and some limitations from security concerns, the HFPB method has been widely adopted by researchers for the analysis of RC columns (Bao and Li, 2010; Shi et al., 2008; Williams and Williamson, 2011) and structures (Jayasooriya et al., 2011; Luccioni et al., 2004; Shi et al., 2010) against blast loads. Despite the popularity of the HFPB method, it must be admitted that this method is complex and requires excessive computational effort. Meanwhile, considering various constitutive models of concrete and steel materials, users may obtain very different analysis results for the same problem if the chosen material models differ. Moreover, being limited by current computer power, it is not practical to adopt the HFPB method when analyzing the responses of large frame structures subjected to blast loadings.

The fiber beam element model is a good compromise between accuracy and efficiency for the analysis of RC structures. It has been widely used in seismic analysis (Neuenhofer and Filippou, 1997; Taucer et al., 1991), and a state-of-the-art flexure–shear fiber beam element has been reviewed by Ceresa et al. (2007). However, only a few researchers have attempted to analyze RC members under blast loads using the fiber beam element model. Valipour et al. (2009) proposed an improved fiber beam element model, which introduced a shear cap at the section level to account for possible shear failure, but it was unable to calculate the correct shear demand excited by a blast load. Guner and Vecchio (2012) developed a nonlinear dynamic analysis method for shear-critical frames under impact, blast, and seismic loads, which had the primary advantage of accounting for shear effects, coupled with axial and flexural behaviors. Nevertheless, the applications for RC components and structures under blast loading are still absent and the damage effect has not been well considered. Therefore, it is necessary to develop an appropriate numerical model based on the traditional fiber beam (TFB) model in the general FE code that can consider shear and damage effects and be suitable for practical applications ranging from nonlinear dynamic response to the progressive collapse analysis of RC structures.

This study aims to improve the TFB model by considering shear and damage effects for the nonlinear dynamic analysis of RC members and RC structures under blast loading and to capture the global response (diagonal shear, flexure–shear, and flexure) of the members rather than local effects such as spalling. In section “Proposed numerical model,” the proposed numerical model is described in detail. The proposed and HFPB models are calibrated in section “Model validations,” using an RC column in a field test. Then, the reliability and efficiency of the proposed model are demonstrated thoroughly via analysis of an RC beam and an RC frame structure in sections “Numerical study of a blast-loaded RC beam,”“Response of a six-story RC frame structure under blast loads,” and “Time consumption.” Discussions and conclusions from the numerical results are summed up in the end.

Proposed numerical model

The proposed numerical model is based on the TFB model available in the general FE code LS-DYNA (Hallquist, 2006), which is incrementally objective (rigid body rotations do not generate strains, which allows for the treatment of finite strains that occur in many practical applications) and computationally efficient. Furthermore, it also includes finite transverse shear strains that take the shear effect of components into consideration if a proper shear model is adopted. Moreover, applications to structural analysis (even progressive collapse) are straightforward in the framework of LS-DYNA when the proposed model is implemented.

RC components generally suffer severe shear deformations and material damage due to intense blast loads, possibly leading to the progressive collapse of RC structures. As a consequence, the primary function of the proposed model is to properly account for the shear and damage effects of RC composite material. Moreover, correctly considering the shear behavior of longitudinal reinforcement bars is also essential under large shear strains. These can be achieved using the modified compression-field theory (MCFT) (Vecchio and Collins, 1986) with a rate-dependent concrete damage model and bilinear steel model in the principal directions. Moreover, the constitutive behaviors for longitudinal reinforcement bars are condensed from the 3D stress–strain of the isotropic hardening plasticity model to account for correct shear demands under large shear strains.

Constitutive models

To take the damage effects of RC components into account in the proposed model for blast analysis, the concrete damage model proposed by Mazars and Pijaudier-Cabot (1989) is adopted. The model is described by the stress–strain relation

σ = (1 - D) C : ε

(1)

damage law

D = g (κ)

(2)

and loading–unloading conditions

f (ε, κ) = ε_{eq} (ε) - κ ⩽ 0, \overset{\cdot}{κ} ⩾ 0, \overset{\cdot}{κ} f (ε, κ) = 0

(3)

with

ε_{eq} = \sqrt{\sum_{i = 1}^{3} {(〈 ε_{i} 〉)}^{2}}

(4)

where $σ$ and $ε$ are the stress and strain tensors, respectively; $C$ is the elastic stiffness tensor, D is the scalar damage variable, which varies from 0 to 1 (0 for intact material and 1 for full failure of the material), $ε_{eq}$ is the scalar measure called equivalent strain, $ε_{i} (i = 1, 2, 3)$ are the principal strains, and $κ$ is an internal variable that equals the maximum equivalent strain ever reached in the history of the material.

The shape of the stress–strain curve is controlled by the specific expressions of the damage function in Figure 2(a). In this study, except for the damage law originating from Mazars’ model, a simple damage law in exponential form for tension is also incorporated in the code, as shown in Figure 2(b). The strain rate effect is considered using the empirical relationship proposed by Malvar and Ross (1998) through the adjustment of the parameters in the constitutive model.

Figure 2.

Schematic of concrete damage model: (a) Mazars model and (b) exponential form for tension damage.

When RC members suffer significant shear deformations near the supports, the dowel action of the longitudinal reinforcement bars may contribute much to the shear resistance. Therefore, the isotropic hardening plasticity model (Chen and Han, 2012) is implemented in this study to model the longitudinal reinforcement bars, because the shear stress is considered in the yield function. The yield condition of this model is as follows

ϕ = σ_{i}^{2} - σ_{sy}^{2} = \frac{3}{2} s_{ij} s_{ij} - σ_{sy}^{2} = 0

(5)

where

σ_{sy} = σ_{s 0} + β E_{p} ε_{eff}^{p}

(6)

and the effective plastic strain increment can be derived as follows

Δ ε_{eff}^{p} = \frac{σ_{i} - σ_{sy}}{3 G + E_{p}}

(7)

Scaling back the stress deviators

s_{ij} = s_{ij}^{*} - \frac{3 G Δ ε_{eff}^{p}}{σ_{i}} s_{ij}^{*}

(8)

where $σ_{i}$ is the stress intensity, $σ_{s 0}$ and $σ_{sy}$ are the initial and current yield limit, respectively, $s_{ij}$ is the deviatoric stress tensor, $β$ is the hardening parameter, $E_{p}$ is the plastic hardening modulus, $ε_{eff}^{p}$ is the effective plastic strain, and $s_{ij}^{*}$ is the trial deviatoric stress tensor.

A condensation has to be made to the steel 3D stress–strain relation for adoption in beam elements. The first step is to compute the trial stress, assuming that the incremental strains are elastic. When the trial stress is outside the yield surface, a secant iteration is used to solve equation (8) for the strain increments $Δ ε_{2}$ and $Δ ε_{3}$ to produce zero normal stresses $σ_{2}$ and $σ_{3}$ . The secant iteration scheme is formulated in equation (9), and the normal strains in the y- and z-directions are calculated until the energy norm estimated by equation (10) is kept within the designated tolerance. Moreover, the dynamic increase factor (DIF) of the yield strength of steel is obtained from Malvar (1998), depending on the strain rate

Δ ε_{k}^{n + 1} = Δ ε_{k}^{n - 1} - \frac{Δ ε_{k}^{n} - Δ ε_{k}^{n - 1}}{σ_{k}^{n} - σ_{k}^{n - 1}} σ_{k}^{n - 1} (k = 2, 3)

(9)

\sqrt{σ_{2}^{2} + σ_{3}^{2}} < tolerance (e . g . 1.0 e - 8)

(10)

Stress updating formulation

A fixed shear strain pattern is assumed in the section level, as it is impossible to determine the actual distribution of shear strains and stresses in advance. The fixed strain pattern is inspired by the plane section hypothesis, which tries to establish relationships between the section’s shear deformation and shear strain distribution, sectional internal forces, and shear stress distribution. As mentioned in Garcia and Bernat (2007), although the fixed strain approach does not guarantee compatibility between the fibers, it still gives satisfactory results in the analysis of the shear resistant mechanism of RC sections; in this way, the computational intensity and time decrease dramatically compared with the inter-fiber equilibrium approach. Moreover, this method is stable and applicable for the post-peak stage, which is very necessary for the possible failure of RC components and collapse of RC structures. Therefore, the strain field of a cross section is a superposition of the classical plane section hypothesis for the normal axial strain and fixed constant shear strain, as shown in Figure 3.

Figure 3.

Schematic of Hughes-Liu beam element model.

In this study, the secant stiffness formulation is adopted, because it has the outstanding characteristics of excellent convergence and numerical stability. The stress updating process stated by Guner and Vecchio (2010) is adopted, which is simply described as follows.

The total strains are assumed to consist of concrete net strain and concrete plastic offset strain, and a perfect bond between the concrete and reinforcement bar is assumed, resulting in the following expression

[ε] = [ε_{c}] + [ε_{c}^{p}] = [ε_{s}] + [ε_{s}^{p}] = [\begin{matrix} ε_{x} \\ ε_{y} \\ γ_{xy} \end{matrix}]

(11)

where $ε$ is the total strain; $ε_{c}$ and $ε_{c}^{p}$ are concrete net strain and concrete plastic strain, respectively; $ε_{s}$ and $ε_{s}^{p}$ are steel net strain and steel plastic strain, respectively; $ε_{x}$ , $ε_{y}$ , and $γ_{xy}$ are the components of the total strain in the global system.

$ε_{x}$ and $γ_{xy}$ are known from the main program. To calculate the principal strains in the fiber, $ε_{y}$ should be assumed in advance; then, the principal strains, $ε_{c 1}$ and $ε_{c 2}$ , can be calculated easily. The corresponding principal stresses, $f_{c 1}$ and $f_{c 2}$ , are calculated using the uniaxial constitutive model as described in section “HFPB model.” The concrete and reinforcement material secant moduli are obtained based on Figure 4 as follows

{\bar{E}}_{c 1} = \frac{f_{c 1}}{ε_{c 1}}

(12)

{\bar{E}}_{c 2} = \frac{f_{c 2}}{ε_{c 2}}

(13)

{\bar{G}}_{c} = \frac{{\bar{E}}_{c 1} \times {\bar{E}}_{c 2}}{{\bar{E}}_{c 1} + {\bar{E}}_{c 2}}

(14)

{\bar{E}}_{sx} = \frac{f_{sx}}{ε_{sx}}

(15)

{\bar{E}}_{sy} = \frac{f_{sy}}{ε_{sy}}

(16)

where ${\bar{E}}_{c 1}$ and ${\bar{E}}_{c 2}$ are the secant moduli of concrete in the principal directions, respectively; ${\bar{G}}_{c}$ is the shear modulus; ${\bar{E}}_{sx}$ and ${\bar{E}}_{sy}$ are the secant modulus of steel in the x- and y-directions, respectively.

Figure 4.

Determination of secant moduli for concrete and steel reinforcement.

Then, the concrete material stiffness matrix in the principal directions can be written as follows

[D_{c}]' = [\begin{matrix} {\bar{E}}_{c 1} & 0 & 0 \\ 0 & {\bar{E}}_{c 2} & 0 \\ 0 & 0 & {\bar{G}}_{c} \end{matrix}]

(17)

and can be transformed to the global axes as follows

[D_{c}] = [T_{c}]^{T} \times [D_{c}]' \times [T_{c}]

(18)

where $T_{c}$ is the transformation matrix, given by

[T_{c}] = [\begin{matrix} \cos^{2} θ & \sin^{2} θ & \sin θ \cos θ \\ \sin^{2} θ & \cos^{2} θ & - \sin θ \cos θ \\ - 2 \sin θ \cos θ & 2 \sin θ \cos θ & \cos^{2} θ - \sin^{2} θ \end{matrix}]

(19)

The steel reinforcement material stiffness matrix for the global axes can be assembled as follows

[D_{s}] = [\begin{matrix} ρ_{x} \times {\bar{E}}_{sx} & 0 & 0 \\ 0 & ρ_{y} \times {\bar{E}}_{sy} & 0 \\ 0 & 0 & 0 \end{matrix}]

(20)

where $ρ_{x}$ and $ρ_{y}$ are the smeared reinforcement ratios in the x- and y-directions, respectively.

The resulting composite material stiffness matrix is assembled as follows

[D] = [D_{c}] + [D_{s}]

(21)

The stress matrix required is calculated using the following expression

[σ] = [D] \times [ε] - [σ^{0}]

(22)

where

[σ^{0}] = [\begin{matrix} S_{01} \\ S_{02} \\ S_{03} \end{matrix}]

is the stress matrix derived from the plastic offset strain as shown in Figure 4.

The assumption that there are no clamping stresses in the transverse direction is considered here, as demonstrated in section “Introduction,” such that the stress along the thickness dramatically decreases to zero after the wave reflects and unloads several times. According to equation (22), the expression of $σ_{y}$ can be obtained as follows

σ_{y} = D_{21} \times ε_{x} + D_{22} \times ε_{y} + D_{23} \times γ_{xy} - S_{02}

(23)

After forcing $σ_{y} = 0$ , a new value of $ε_{y}$ is calculated using equation (24)

ε_{y} = \frac{- D_{21} \times ε_{x} - D_{23} \times γ_{xy} + S_{02}}{D_{22}}

(24)

Using the newly calculated $ε_{y}$ , new principal stresses are derived and the above calculations are repeated until a specific tolerance of $σ_{y}$ is achieved. Then, normal stress $σ_{x}$ and shear stress $τ_{xy}$ at the fiber state are updated using the following expressions

σ_{x} = D_{11} \times ε_{x} + D_{12} \times ε_{y} + D_{13} \times γ_{xy} - S_{01}

(25)

τ_{xy} = D_{31} \times ε_{x} + D_{32} \times ε_{y} + D_{33} \times γ_{xy} - S_{03}

(26)

Implementation

The proposed numerical model is implemented in the general FE code LS-DYNA through user-defined interfaces. The flow chart of the numerical solution procedure is presented in Figure 5. Because an explicit formulation is used, there is no need to generate the tangent stiffness matrix. Hence, it is possible to capture the strong nonlinear behavior of the RC members when the material undergoes severe damage.

Figure 5.

Flow chart of the proposed numerical model in LS-DYNA.

Model validations

To validate the proposed model, test results reported in the reference (Crawford et al., 2012) are used. Simulations using the TFB model are also conducted to demonstrate the improvement and accuracy of the proposed model. Moreover, to prove the efficiency of the proposed model, the HFPB method, commonly used by a number of research groups (Bao and Li, 2010; Jayasooriya et al., 2011; Luccioni et al., 2004; Shi et al., 2008, 2010; Williams and Williamson, 2011), is also adopted as the benchmark simulation so that the computational times of different cases, as well as the computational accuracies, can be compared. Therefore, the experimental data from the test are also used to check the reliability of the HFPB model. The properties of the specimen are briefly introduced below.

The CTEST20 column is taken from a series of standalone column tests (Crawford et al., 2012) and has a clear height of 3277 mm. The column has a squared cross section with a size of 355.6 mm. Eight #8 (25.4 mm (1 in.)) longitudinal reinforcements are uniformly spaced along the perimeter of the section, and #3 (9.5 mm (0.375 in)) stirrup reinforcements are distributed along the length with a spacing of 324 mm. An axial load of 444.8 kN is applied to the column first, and then the top end of the column is fixed; finally, the explosives are detonated to generate the nonlinear dynamic response. As the mass of an explosive charge is not released for security concerns, a blast load with a peak overpressure of 49.1 MPa and duration of 0.62 ms, which were computed in the reference, is used in this study.

The unconfined compressive strength of concrete measured in the test is $f'_{c} = 41.4 MPa$ . A yield stress and tensile stress of $f_{y} = 475 MPa$ and $f_{u} = 751 MPa$ , respectively, are used for reinforcement bars in the numerical analysis, as recommended by K&C.

HFPB model

For all the HFPB simulations in this article, concrete is modeled with solid elements and reinforcement bars are modeled with beam elements, while the interaction between concrete and rebar is implicitly captured by the behavior of the concrete. As is well known, the numerical results are mesh-dependent if the descending branch of the constitutive model is considered. On one hand, a very fine mesh is needed for shock wave propagation to transmit the blast energy to the structure. On the other hand, concrete materials soften after the peak stress, which causes localization of the strain to a narrow band if the constitutive model is not equipped with a characteristic length, leading to FE solutions exhibiting spurious mesh sensitivity, losing objectivity with regard to the choice of the mesh (Pijaudier-Cabot and Bažant, 1987).

The third release of the K&C model, generally used in blast analysis for concrete material and available in LS-DYNA as *MAT_CONCRETE_DAMAGE_REL3 (Hallquist, 2007), is chosen for the simulation of concrete material because it has been proven reliable in many previous studies. Shear strength enhancement factor versus effective strain rate is defined by a curve, which is obtained from Malvar and Ross (1998). There are a total of three parameters in this concrete model that controls the damage evolution. They are b₁, b₂, and $ω$ , that is, compression damage, tensile damage, and volume expansion, respectively.

As for the compression softening behavior, the default value of compression softening parameter b₁ is correct only for a 4-in (101.6 mm) element, which means it must be adjusted by the users to reflect the element sizes. An empirical expression for b₁ was developed by K&C (Crawford et al., 2012), which is as follows

b_{1} = 0.34 h + 0.79

(27)

where h is the element size (in).

The tensile damage exhibited by the model is considered using the fracture energy and tension softening parameter b₂ and is affected by the element and aggregate sizes. As is stated by Magallanes et al. (2010), the regularization method ceases to be valid when the elements are much smaller than the localization width $w_{L}$ (the parameter “LOCWID” on the input cards), and as is well known, quite small element sizes are especially required in the blast analysis of structural responses. Therefore, an expression based on Malvar’s work (Wu et al., 2014) is used in this study

\begin{matrix} b_{2} = (0.09 w_{L}^{2} - 0.98 w_{L} + 3.06) \\ (1 - 0.004 f'_{c}^{2} + 0.097 f'_{c} - 0.484) \end{matrix}

(28)

where $w_{L}$ is in the units of “inch” and $f'_{c}$ is the unconfined compressive strength of concrete in units of “ksi.”

$ω$ is the very important parameter that controls the volume expansion of concrete to affect the shear-dilatancy behaviors, which has been discussed deeply by Crawford et al. (2013). The suggested value is 0.75 without confinement and 0.90 in a confined case.

Reinforcement bars play an important role in the analysis of RC components and structures under blast loading. A true stress and strain relationship is generally needed in simulations. After the engineering stress–strain data are transformed into a true form, it is implemented in LS-DYNA using the material model *MAT_PIECEWISE_LINEAR_PLASTICITY (Hallquist, 2007), which allows an arbitrary stress versus strain curve to be defined and a plastic failure strain to be specified to mimic the fracture. The strain rate effect is accounted for using the Cowper and Symonds model (Hallquist, 2006).

TFB model

To illustrate the advantages of the proposed numerical model, the TFB model available in LS-DYNA is also used for all numerical examples. In this study, the material model *MAT_CONCRETE_BEAM is used for concrete fibers. An effective stress versus effective plastic strain curve is defined to model the concrete nonlinear behavior, and the strain rate effect is also considered using a load curve obtained from Malvar and Ross (1998). The longitudinal reinforcement fibers are modeled using *MAT_PIECEWISE_LINEAR_PLASTICITY, as stated in section “HFPB model.” Note that the Hughes-Liu beam element formulation requires that different material models used in the same section have the same number of history variables; here, both have 5 (Hallquist, 2007).

Validation of models

Mesh size convergence tests are carried out using mesh sizes of 12.7, 25.4, and 50.8 mm for both the concrete and reinforcements in the HFPB model, and it is found that the mesh size of 25.4 mm is appropriate for the simulations. And the mesh size effect on the material softening is also considered by adjusting the damage parameters b₁, b₂, and $ω$ , according to the corresponding equations in section “HFPB model” (listed in Table 1), while the other concrete parameters are automatically generated by LS-DYNA. To capture the possible spalling of the concrete material, the erosion technique is implemented for the concrete cover. A maximum principle stress of 10 MPa is defined as the erosion criterion in this study. The top end is free to move along the z-axis to allow the initial axial load to be applied in the first load stage, and then it is fixed at the time when the blast load is applied.

Table 1.

Values for key parameters in HFPB model.

b ₁	b ₂	b ₃	$ω$	Loc width (mm)
1.13	2.07	1.15	0.75	25.4

A beam element size of approximately 83 mm (a total of 40 elements) is adopted in both the TFB model and the proposed numerical model. The parameter that controls the strain rate effect in the proposed numerical model is set to 1 (ignore strain rate effect for 0). The key parameters for concrete materials are listed in Tables 2 and 3.

Table 2.

Values for key parameters in TFB model.

$ε_{p}$	$σ$ (MPa)
0.0	12.4
9.3e−4	41.4
0.1611	0.0
0.2	0.0

Table 3.

Values for key parameters in proposed model.

$ε_{0}$	$A_{c}$	$B_{c}$	$ε_{f}$
0.0015	1.1054	1178.51	0.0058

The transverse displacements at the mid-height of the column from the test and simulations are compared in Figure 7. As can be observed in Figure 6, both the HFPB model and the proposed model predict the maximum displacement at mid-span excellently with the test result. However, the maximum displacement calculated by the TFB model is only 64.8 mm, which is much smaller than the test result of 175.3 mm. The numerically simulated residual displacements are 113, 98, and 55 mm for the HFPB model, the proposed model, and the TFB model, respectively, and the tested residual displacement is 109 mm. Obviously, the HFPB model predicts the best result for both the maximum mid-span displacement and the residual displacement, as the concrete material model captures both the shear effect and damage level very well. The proposed model also satisfactorily predicts the maximum mid-span displacement, but the residual displacement is 10% less compared with the test results. The most possible reason is that the Mazars model used in the principal directions of the RC material does not account for any plastic strain in concrete, which leads to a smaller residual value. A more precise concrete plastic damage model may improve the predictions, although it may also increase both the complexity and computational cost.

Figure 6.

Comparison of mid-height displacement time history.

Figure 7 shows the base shear force histories obtained from the HFPB model and the proposed model. The peak values predicted by the two models are 2137.1 and 2077.4 kN, respectively, and the difference is only 2.8%, indicating a good performance of the developed model. Note that the shear force calculated by the TFB model is not included in this figure because the peak shear is 6800 kN, which is a considerably larger value compared with 2137.1 kN.

Figure 7.

Comparison of shear force histories.

A very important aspect of analyzing RC columns under blast loads involves compression-membrane behavior. It is a resistance mechanism attributed to the fixity of a column’s supports, also referred to as arching, that generates from column growth. As stated earlier, the top end of the column is fixed when the lateral blast load is applied; thus, the compression-membrane behavior can be captured by all three models. Figure 8 illustrates that a large compressive axial load exists in the column, causing an overestimation of the compression-membrane mechanism when the TFB model is used. The axial force history curves of both the HFPB model and the proposed model show that the axial force first increases in compression and then changes to tension, which implies the transformation of the compression-membrane mechanism to a tension-membrane mechanism, indicating the occurrence of a shear failure under blast loading. Whereas the TFB model fails to capture such shear failure, as shown in Figure 8, the proposed model successfully predicts the transformation, with, however, significantly lower tensile force (1092.6 kN) compared with 2456.5 kN of the HFPB model. This can possibly be attributed to the differences in considering the concrete damage effect between the two models. Figure 9 shows the final failure mode of the column from the proposed model, HFPB model, and TFB model. Apparently, the HFPB model well predicts a diagonal shear failure mode, and even the spalling of the concrete material is well captured. The failure mode obtained from the proposed model matches very well with the HFPB model, whereas only a flexural deformation mode is observed from the TFB model, demonstrating a significant improvement and the accuracy of the proposed model when it is used in the blast analysis of RC columns under severe blast loading.

Figure 8.

Comparison of axial force histories.

Figure 9.

Final deformation state of the column, obtained from the (a) proposed model, (b) HFPB model, and (c) TFB model.

Numerical study of a blast-loaded RC beam

In section “Validation of models,” the HFPB model and the proposed model were validated through the analysis of a RC column subjected to a specific blast loading scenario. To comprehensively demonstrate the feasibility of the proposed model, a 4-m long RC beam with two fixed ends (Figure 10) subjected to blast loading is used as an example. The beam element has a rectangular cross section of 200 mm × 400 m, with four 20 mm longitudinal reinforcement bars at each corner. Stirrup reinforcements with a diameter of 10 mm are placed along the beam with a constant spacing of 200 mm. The load is assumed to be uniformly distributed over the whole length of the beam. The natural vibration period of the beam is 10.8 ms, which is obtained from modal analysis in LS-DYNA. Four load cases are considered in this study, as outlined in Table 4, of which case 1 corresponds to an impulsive range with a relatively small peak overpressure, case 2 follows in the impulsive range but with a large peak overpressure, case 3 doubles the peak overpressure of case 2, and case 4 belongs to the dynamic range, with the same impulse of case 2 but with a much smaller overpressure and much longer duration (Table 4).

Figure 10.

Geometry and section details of the RC beam.

Table 4.

Blast load parameters for different cases.

Case	Peak overpressure (MPa)	Duration (ms)
1	1	1
2	20	1
3	40	1
4	1	20

A 20 mm element size is chosen for the HFPB model based on a convergence study, and the damage parameters for the concrete constitutive model are determined according to section “HFPB model.” For TFB and the proposed models, an element length of 50 mm is adopted, as it was found that using smaller element size had a negligible effect on the results, and the softening parameters (plastic strain in the descending branch for the TFB model and $ε_{f}$ for the proposed model) of the constitutive models are adopted according to the crack band theory (Bažant and Oh, 1983).

Figure 11 shows the time histories of mid-span displacement, shear force, and axial force from the simulations using different models under loading case 1. As observed, the results match very well for all the models. This is because the impulsive load case has relatively small peak overpressure and duration, which means a small impulse. The peak displacement is only approximately 1.2 mm and the shear stress in the outer concrete fiber is below 1.5 MPa; thus, the beam is still in the linear elastic stage.

Figure 11.

Displacement, shear force and axial force time history for case 1: (a) mid-span displacement, (b) shear force, and (c) axial force.

For load case 2, the peak overpressure is increased to 20 MPa, which is intense enough to induce strong nonlinear behavior of the beam, especially severe shear damage near the support. Figure 12 shows the comparison of the results of the beam subjected to load case 2 using different models. It can be observed from Figure 12(a) that the peak displacement at mid-span predicted by the HFPB, TFB, and proposed models are 57.1, 46.6, and 58.6, respectively. Compared to the HFPB model, the TFB model underestimates the response by 18.4%. The proposed model predicts a very close value compared with the HFPB model. The comparison of shear force histories near the support in Figure 12(b) shows that the TFB model generates larger maximum shear forces in the loading phase than the HFPB and proposed models, while the proposed model matches the HFPB model quite well. The axial force time histories shown in Figure 12(c) indicate that the shear failure of the RC beam is captured by both the proposed and HFPB models because the axial force decreases dramatically, while the TFB model cannot. The reason for the huge difference is obvious. The RC beam undergoes severe shear damage near the support, which can be captured easily by the HFPB model, because the constitutive model of concrete used has three shear failure surfaces. Meanwhile, because the shear and damage effects are also incorporated in the proposed model, both the maximum mid-span displacement and shear demand show a good fit, except for the residual displacement due to the ignorance of plastic deformation in Mazars’ model, as mentioned before. Nevertheless, the concrete constitutive model used in the TFB model is elasto-plastic, which cannot account for the nonlinear shear and damage effects of concrete material, resulting in underestimation of the deformation and overestimation of the shear demand near the beam support.

Figure 12.

Displacement, shear force and axial force time history for case 2: (a) mid-span displacement, (b) shear force, and (c) axial force.

A more intense impulsive load is considered in case 3 to demonstrate the capability of the proposed model in capturing the shear and damage effect under severe blast loading. As observed in Figure 13(a), the beam exhibits a shear-flexural failure mode. The maximum mid-span displacements are 278.4 and 275.7 mm for the HFPB model and the proposed model, respectively, with a discrepancy less than 1%. Moreover, the shear demands also match very well for these two models, with values of 1720.1 and 1817.2 kN, respectively, a difference of approximately 5%. The axial loading capacity of the RC beam calculated by both the proposed and HFPB models is almost lost at 10 ms, as indicated in Figure 13(c). Because the overpressure and impulse are both twice those in load case 2, the shear damage is much more severe.

Figure 13.

Displacement, shear force and axial force time history for case 3: (a) mid-span displacement, (b) shear force, and (c) axial force.

Load case 4, with a duration time of 20 ms, belongs to the dynamic range, as the natural vibration period of the RC beam is 10.8 ms. The ratio of loading duration and natural period is <2. Figure 14 shows that the first peak value of both the displacement at mid-span and shear force near the support fit quite well for all three models, while the first peak value of the axial force shows some difference between the proposed and HFPB models. And the values begin to diverge markedly from each other after the first peak, resulting in different residual displacements and free vibrations. This is caused by certain flexural damages of the RC beam. As shown in Figure 14(a), the free vibration periods of the beam are 14.22 and 15.11 ms for the HFPB model and the proposed model, respectively, while that for the TFB model still remains at approximately 10.8 ms, which indicates that no damage can be captured because the constitutive models do not account for any damage effect.

Figure 14.

Displacement, shear force and axial force time history for case 4: (a) mid-span displacement, (b) shear force, and (c) axial force.

Figure 15 shows the shear stress time history of an outer concrete fiber near the beam support, obtained from the proposed and TFB models. Because the response of the beam is in the elastic stage for case 1, its shear stress history is not plotted in the figures. The maximum shear stresses of case 2 and case 3, calculated by the TFB model, in Figure 15(b) reach 20 and 40 MPa, respectively, which are overly large values for concrete because the nonlinear shear behavior is not considered, leading to underestimated deformations and overestimated shear demand, as stated before. In contrast, the values obtained by the proposed model in Figure 15(a) are 13.0 and 15.2 MPa. For load case 4, because only a weak nonlinear behavior is expected near the support, the shear stress calculated by the proposed and TFB models matches relatively well.

Figure 15.

Shear stress time history of outer concrete fiber near column support: (a) proposed model and (b) TFB model.

Response of a six-story RC frame structure under blast loads

To further demonstrate the accuracy of the proposed model in the analysis of structural responses under explosive loadings, a six-story RC frame structure is modeled in this section. The six-story RC frame structure has three bays with a span of 4 m along the x-axis and two bays with a span of 4 m along the y-axis, as shown in Figure 16. The dimensions of all beams and columns are 300 mm × 300 mm. Beams and columns consist of four longitudinal reinforcements placed at the four corners, with diameters of 16 and 20 mm, respectively, and stirrup reinforcements with a diameter of 10 mm are spaced 200 mm apart for both beams and columns. The floor slab is 100 mm thick, with a steel layer in the middle. The compressive strength of concrete is assumed to be 30 MPa, and a bilinear stress–strain relationship is used for reinforcements, with a yield strength of 335 MPa and ultimate strength of 480 MPa. The static load applied to the slab is 2.0 kN/m² and that to the top face of the beams is 60 kN/m², considering the weight of the infill walls. The blast load is assumed to have a peak overpressure of 20 MPa and duration time of 1 ms. Only one column on the ground floor (column C1 in Figure 16) has the blast load applied to it for simplicity.

Figure 16.

Schematic of the three FE models for the RC frame structure.

Three FE models are developed to study the response characteristics of the RC frame under blast loading. As outlined in Table 5, M1 is regarded as an efficient model, with a fiber beam element representing both beams and columns and layered-shell element representing floor slabs, and the newly proposed model is used for column C1, which is directly subjected to the blast load; M2 is considered as a benchmark model, with solid elements representing concrete and beam elements representing reinforcement bars. To demonstrate the superiority of the proposed model to the TFB model, the efficient model, with column C1 simulated by the TFB model, is also analyzed, denoted as M3.

Table 5.

Descriptions of the three FE models.

FE model	Description
M1	TFB model for beams and columns, layered-shell element for slabs, and the proposed model for column C1 (14,400 shell elements and 3120 beam elements)
M2	Solid element for concrete, beam element for longitudinal and transverse reinforcements (822,816 solid elements, 128,000 shell elements, and 226,848 beam elements)
M3	TFB model for beams and columns, layered-shell element for slabs (14,400 shell elements and 3120 beam elements)

An overview of M1, M2, and M3 used in the analysis is shown in Figure 16. Label A is the node at the middle of the loaded column, while label B is the node at the middle of the joint above the loaded column. To reduce computational effort, from 0 to 50 ms, a combination of dead loads and live loads is applied to the frame; then, the blast load acting on the column is applied at t = 200 ms, and the analysis is terminated at t = 1000 ms.

Figure 17 gives the time history of the transverse displacement of node A along the y-direction; as shown in this figure, the first peak values are +88.1, +76.1, and +43.6 mm for M1, M2, and M3, respectively. Finally, residual displacements of +64.7, +65.1, and +33.7 mm are captured by the models. When node A undergoes transverse displacement, the vertical displacement of node B is increased along the z-direction; as can be observed in Figure 18, the values are approximately +11.3 and +8.0 mm for M1 and M2, respectively, while it is approximately −4.1 mm for M3. This indicates that an overly compression-membrane mechanism is formed through beam growth in the TFB model, because it cannot capture any shear damage near the supports. The axial load initially supported by column C1 has to seek an alternative load path when the column loses its load capacity. As presented in Figure 19, column C1 still has certain residual axial capacity for M1 and M2, which means that only part of the axial load sustained by column C1 needs to be redistributed by the adjacent columns. Note that the axial load sustained by column C1 in M3 is increased due to the added resistance from the compression-membrane mechanism, as can be observed in Figure 19. The most possible reason for the differences between M1 and M2 can be attributed to the difference of blast damage captured by them, which results in a different post-blast behavior. Meanwhile, simplifications of the joint regions and floor systems in the proposed model may also contribute to the divergent results, and a more reliable model that accounts for these critical aspects is currently under research. However, the improvements of the proposed model are apparent for the analysis of large-scale RC frames under blast loading.

Figure 17.

Transverse displacement time history of node A.

Figure 18.

Vertical displacement time history of node B.

Figure 19.

Axial force time history of column C1.

Time consumption

The reliability of the proposed model has been demonstrated in the previous numerical analyses. The above numerical simulations were all carried out on a desktop computer with 3.4 GHz CPU frequency and 8G memory. As shown in Table 6, the total elapsed time of the HFPB model and the proposed model of the column calculated in section “Validation of models” are 39 min 1 s and 51 s, respectively. For the RC beam in section “Numerical study of a blast-loaded RC beam,” only the computational time of case 2 is compared, because the four cases are almost identical with respect to time consumption. They are 13 min 52 s and 5 s for the HFPB model and the proposed model, respectively, proving a very significant advancement of the proposed model in saving computational time while assuring accuracy. The simulation time for the benchmark model and the proposed model of the six-story RC frame structure in section “Response of a six-story RC frame structure under blast loads” is 15 h 29 min and 3 h 7 min, respectively. Obviously, the proposed model greatly reduces the amount of computational time for the analysis, making it a promising model for the possible progressive collapse analysis of large-scale RC structures under blast loading.

Table 6.

Comparisons of computational time between the proposed and HFPB models.

	RC column	RC beam (case 2)	RC structure
Proposed model	51 s	5 s	3 h 7 min
HFPB model	39 min 1 s	13 min 52 s	15 h 29 min

Conclusion

This article proposes a numerical model to account for both diagonal shear behavior and damage effects (excluding material losses from spalling) based on the TFB model. The proposed model is implemented in the general FE code LS-DYNA, making it convenient for the dynamic analysis of both RC components and structures under blast loading. The diagonal shear failure of RC components and the transformation of resistance mechanisms from compression-membrane to tension-membrane are captured efficiently. The computational time is substantially reduced in all the analyses compared with the HFPB model. The reliability and efficiency of the proposed model provide great potential for progressive collapse analysis and the design of RC structures under blast loads.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors gratefully acknowledge the financial support of the National Natural Science Foundation of China under grant number 51238007 and the Special Research Fund for Doctoral Disciplines of Colleges of Ministry of Education of China under grant number 20120032110049 for carrying out this research.

References

Bao

(2010) Residual strength of blast damaged reinforced concrete columns. International Journal of Impact Engineering 37(3): 295–308.

Bažant

(1983) Crack band theory for fracture of concrete. Matériaux et construction 16(3): 155–177.

Biggs

(1964) Introduction to Structural Dynamics. New York: McGraw-Hill.

Carta

Stochino

(2013) Theoretical models to predict the flexural failure of reinforced concrete beams under blast loads. Engineering Structures 49: 306–315.

Ceresa

Petrini

Pinho

(2007) Flexure-shear fiber beam-column elements for modeling frame structures under seismic loading—state of the art. Journal of Earthquake Engineering 11(S1): 46–88.

Chen

Han

(2012) Plasticity for Structural Engineers. New York: Springer Science & Business Media.

Crawford

Magallanes

(2011) The effects of modeling choices on the response of structural components to blast effects. International Journal of Protective Structures 2(2): 231–266.

Crawford

Magallanes

. (2012) Use and Validation of the Release II K&C Concrete Material Model in LS-DYNA. Glendale, CA: Karagozian & Case.

Crawford

Magallanes

. (2013) The importance of shear-dilatancy behaviors in RC columns. International Journal of Protective Structures 4(3): 341–378.

10.

Dragos

(2014) Interaction between direct shear and flexural responses for blast loaded one-way reinforced concrete slabs using a finite element model. Engineering Structures 72: 193–202.

11.

El-Dakhakhni

Changiz Rezaei

Mekky

. (2009) Response sensitivity of blast-loaded reinforced concrete structures to the number of degrees of freedom. Canadian Journal of Civil Engineering 36(8): 1305–1320.

12.

Flügge

Zajac

(1959) Bending impact waves in beams. Archive of Applied Mechanics 28(1): 59–70.

13.

Garcia

Bernat

(2007) Shear-bending-torsion interaction in structural concrete members: a nonlinear coupled sectional approach. Archives of Computational Methods in Engineering 14(3): 249–278.

14.

Guner

Vecchio

(2010) Pushover analysis of shear-critical frames: formulation. ACI Structural Journal 107(1): 63–71.

15.

Guner

Vecchio

(2012) Simplified method for nonlinear dynamic analysis of shear-critical frames. ACI Structural Journal 109(5): 727–737.

16.

Hao

(2014) Influence of the concrete DIF model on the numerical predictions of RC wall responses to blast loadings. Engineering Structures 73: 24–38.

17.

Hallquist

(2006) LS-DYNA Theory Manual. Livermore, CA: Livermore Software Technology Corporation.

18.

Hallquist

(2007) LS-DYNA Keyword User’s Manual. Livermore, CA: Livermore Software Technology Corporation.

19.

Jayasooriya

Thambiratnam

Perera

. (2011) Blast and residual capacity analysis of reinforced concrete framed buildings. Engineering Structures 33(12): 3483–3495.

20.

Jones

(2011) Structural Impact. Cambridge: Cambridge University Press.

21.

Krauthammer

(2008) Modern Protective Structures. New York, USA: CRC Press.

22.

Krauthammer

(2014) Foreword: challenges for advancing the performance of military concrete. Magazine of Concrete Research 66(1–2): 3–5.

23.

Krauthammer

Assadi-Lamouki

Shanaa

(1993a) Analysis of impulsively loaded reinforced concrete structural elements—I. Theory. Computers & Structures 48(5): 851–860.

24.

Krauthammer

Assadi-Lamouki

Shanaa

(1993b) Analysis of impulsively loaded reinforced concrete structural elements—II. Implementation. Computers & Structures 48(5): 861–871.

25.

Krauthammer

Shanaa

Assadi

(1994) Response of structural concrete elements to severe impulsive loads. Computers & Structures 53(1): 119–130.

26.

Shi

(2015) Blast Analysis of Building Structures. Beijing, China: Science Press.

27.

Luccioni

Ambrosini

Danesi

(2004) Analysis of building collapse under blast loads. Engineering Structures 26(1): 63–71.

28.

Magallanes

Malvar

. (2010) Recent improvements to release III of the K&C concrete model. In: 11th international LS-DYNA users conference, Michigan, USA, 6–8 June 2010, pp. 3–37.

29.

Magnusson

Ansell

Hansson

(2010) Air-blast-loaded, high-strength concrete beams. Part II: numerical non-linear analysis. Magazine of Concrete Research 62(4): 235–242.

30.

Malvar

(1998) Review of static and dynamic properties of steel reinforcing bars. ACI Materials Journal 95(5): 609–614.

31.

Malvar

Ross

(1998) Review of strain rate effects for concrete in tension. ACI Materials Journal 95(6): 735–739.

32.

Mazars

Pijaudier-Cabot

(1989) Continuum damage theory-application to concrete. Journal of Engineering Mechanics 115(2): 345–365.

33.

Neuenhofer

Filippou

(1997) Evaluation of nonlinear frame finite-element models. Journal of Structural Engineering: ASCE 123(7): 958–966.

34.

Pijaudier-Cabot

Bažant

(1987) Nonlocal damage theory. Journal of Engineering Mechanics 113(10): 1512–1533.

35.

Shi

Hao

(2008) Numerical derivation of pressure–impulse diagrams for prediction of RC column damage to blast loads. International Journal of Impact Engineering 35(11): 1213–1227.

36.

Shi

Hao

(2010) A new method for progressive collapse analysis of RC frames under blast loading. Engineering Structures 32(6): 1691–1703.

37.

Taucer

Spacone

Filippou

(1991) A Fiber Beam-Column Element for Seismic Response Analysis of Reinforced Concrete Structures. Berkeley, CA: Earthquake Engineering Research Center, College of Engineering, University of California, Berkeley.

38.

Timoshenko

(1921) LXVI. On the correction for shear of the differential equation for transverse vibrations of prismatic bars. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 41(245): 744–746.

39.

Valipour

Huynh

Foster

(2009) Analysis of RC beams subjected to shock loading using a modified fibre element formulation. Computers and Concrete 6(5): 377–390.

40.

Vecchio

Collins

(1986) The modified compression-field theory for reinforced concrete elements subjected to shear. ACI Journal Proceedings 83(22): 219–231.

41.

Williams

Williamson

(2011) Response of reinforced concrete bridge columns subjected to blast loads. Journal of Structural Engineering: ASCE 137(9): 903–913.

42.

Sheikh

(2013) A finite element modelling to investigate the mitigation of blast effects on reinforced concrete panel using foam cladding. International Journal of Impact Engineering 55: 24–33.

43.

Crawford

Lan

. (2014) Validation studies for concrete constitutive models with blast test data. In: 13th international LS-DYNA users conference, Michigan, USA, 8–10 June 2014, pp. 1–12.