Blast Load Analysis and Simulation of Unreinforced Concrete Masonry

This study presents an innovative blast analysis tool that can be used in assessing blast performance of concrete masonry walls. The analysis tool proposes a structural model of unreinforced masonry (URM) walls, and a simplified procedure to generate load histories on structures subjected to blast loading. This study focuses on hollow concrete masonry units (CMU) walls. Three-dimensional masonry walls are modeled using discrete rigid bodies connected by nonlinear springs. Connections between two masonry blocks attached by mortar joint are modeled by treating adjoining springs as springs in series. The blast load histories are calculated using UFC 3-340-02. The proposed tool generates individual load-time histories for each discrete elements of a wall. The individual load-time histories, and the peak pressure distribution on a wall surface, are compared to experimental measurements to validate the blast-load-generating algorithm. The two modules (blast-load, and numerical model generator) are then used to simulate an actual experiment on URM wall for validation. For this study, experiments on URM wall arching is simulated. The proposed tool is able to approximate the wall response, particularly, the maximum positive and negative displacements. The good agreement proved the applicability of the proposed analysis tool to the analysis of URM walls.


Introduction
Blast events have become a common occurrence in today's world. Explosions that of interest to designers can either be intentional or accidental. Examples of intentional explosions are bomb blasts, which have become a common method for terrorists in spreading panic. Examples of accidental explosions are gas explosions in residential homes with gas leakage.
Since blast events have become more common, there is a growing interest in the analysis of structures under blast loading. Current industry practices in blast analysis typically use static loading. Most analysis approach only considers the structural frames of the structure and masonry walls are often disregarded in the analysis. Studies show, however, that even unreinforced masonry walls offer blast protections that are greater than expected [1]. Another approach in blast analysis is with the use of more advanced computational fluid dynamics (CFD) software packages coupled with dynamic structural response. This approach is less conservative, however, most commercial packages used to implement such analyses are costly. Also, there are limitations when it comes to conducting blast experimentations due to high cost. It is much more economical to perform numerical studies in place of actual blast experiments.
This study aims to simulate unreinforced masonry walls subjected to blast loading, for the purpose of analyzing the performance of said walls. Particularly, the study aims to (1) develop an algorithm that automatically generates time-varying forces acting on individual elements of a discretized wall, (2) propose a numerical model of a wall made of hollow concrete masonry units (CMU) that is suitable for simulating response of wall to blast loading, and (3) 3 properties of individual discrete elements will result in a set of macroscale mechanical properties different from the actual test results.
Gu, et al. [9], proposed formulas to derive mesoscale mechanical properties from values derived from tests by numerical regression analysis based on size effects laws in 2D models. The proposed numerical model discrete element mechanical properties, specifically, compressive strength (݂ , in ‫,)ܽܲܯ‬ tensile strength (݂ ௧ ), Poisson's ratio (ߥ ), and elastic modulus ‫ܧ(‬ ), with size effects adjustments are as shown below.

Arching response of URM walls subjected to blast loading
In a series of published studies by Abou-Zeid, et al. [1] [10] [11], the effects of enforcing composite action between the infill panel and the surrounding frame as a simple and economical approach for hardening existing buildings are investigated. The schematic drawing of the specimen is shown in Figure  1 (a). Each specimen is made of CMU, ͳͻͲ ݉݉ ‫ݔ‬ ͳͻͲ ݉݉ ‫ݔ‬ ͵ͻͲ ݉݉ in size, and 11 courses high and 2.5 courses wide, which are face-shell mortar bedded with ͳͲ ݉݉ mortar joints. Abou-Zeid [10] gave the mechanical properties as per the material tests conducted for the study (see Figure 1 (b)).
The hardening effect that is described above is due to the concept of a wall acting as an arch. This concept assumes that after the development of tension cracks at the ends and the center of the wall, the wall acts as two rigid segments, with each segment rotating about its end until either the masonry crushes or the two segments snap through [10]. Several specimen walls are subjected to a number of different weights of ammonium nitrate/fuel oil (ANFO), with TNT equivalency factor of 0.80 [12]. Figure 1 (c) shows the matrix of the explosive weights that are ignited. For this study, W2 to W5 walls are simulated.

Blast load computation
This study makes use of an automation algorithm to solve for the blast load-time histories of an explosive event. The load-time histories are built using the method prescribed by UFC 3-340-02 [2]. The full outline of the procedure for solving for the individual points in the load-time history are laid out by Tan [13], and is implemented using xPLoDE (Explosive Pressure Load on Discrete Elements), which is specifically developed for this study.
The blast loading on a wall depends on the relative position of the blast source. If the external surface of a wall is acted upon by direct blast pressure -if the blast source is positioned such that the incident shock wave has the external surface of the wall in consideration as the first surface that it impinges (other than the ground) -it experiences higher pressures compared to walls that are shielded from the direct pressure.
For walls experiencing direct pressures, UFC 3-340-02 outlines the procedure in Section 2-15.3.2 of the code, where this case is described as front wall loading. This load-time history prescribed by UFC 3-340-02 is a simplification of the actual blast loading history. For roofs and walls next to the front walls -called side walls by UFC 3-340-02 -the pressures are much lower.

Numerical model
The proposed method represents a structure as an assembly of rigid elements that are connected by multiple springs that models the response of a structure from the elastic range up to fracture. For this study, an algorithm is built to automate the assembly of the numerical model with springs. The algorithm is called MWBLOCk (Masonry Wall Block Layer with On-Surface Spring Constant k).
This method is based on the rigid body-spring method (RBSM) in that it models the connections between elements as multiple springs, and it is also based on the discrete element method (DEM) in that it models the contact and re-contact between elements. The proposed rigid body-spring -discrete element method (RBS-DEM) is implemented using the virtual environment of wxAERO. The spring model used in this study is based on the proposal by Furukawa, et al. [14], which was developed for unreinforced solid masonry blocks. For this study, however, since hollow masonry blocks are investigated, a modified spring constant derivation is proposed.

Individual spring constants for solid elements
If two adjacent elements are continuous, a restoring force is generated. The force is calculated by multiplying the spring constant between elements with the relative displacement in the normal ‫ݑ(‬ ) and shear ‫ݑ(‬ ௦ ) direction of the surfaces. Normal force ݂ and shear force ݂ ௦ (see Figure 2 (b)) due to the displacement, are assumed to act at the center of the surface segment ‫.ܣ݀‬ Since the element is rigid, it is assumed that the relative displacement is zero inside the pyramids and only has the values of ‫ݑ‬ and ‫ݑ‬ ௦ at the surface segment [14].

Springs between elements
CMU walls are built by joining units with mortar that acts as binder between units. When a connection between two elements is continuous, spring constants between them are equivalent to their individual spring constants in series. For two CMU sub-elements in contact, and bound with mortar, the equivalent spring constant per area ݇ ത is: ͳ where ݇ ǡ ݇ ǡ ݇ ெ -spring constants per area of elements A, B, and mortar. Substituting equations (3) to (4), the equivalent normal spring constant per area ݇ ത and equivalent shear spring constant per area ݇ ത ௦ , between element A, mortar element M, and element B, are defined by Furukawa, et al. [14], as: where ݈ ǡ ݈ -element centroid distance to surface, for elements A and B, ‫ݐ‬ ெ -mortar thickness, ‫ܧ‬ ǡ ‫ܧ‬ ǡ ‫ܧ‬ ெ -Young's modulus for elements A, B, and mortar, and ߥ ǡ ߥ ǡ ߥ ெ -Poisson's ratio for elements A, B, and mortar.
Note that in the proposed method, a single CMU is considered as rigid element, therefore, the mortar properties are incorporated in the model by introducing a separate spring connecting the springs of adjacent CMU elements. In equation (5), the middle term of the denominator corresponds to the spring constant of the mortar.

Spring constants for hollow concrete masonry walls
Furukawa, et al. [14], derived spring parameters by summing forces acting on a volume pyramid formed by a segment area and the centroid of the block in consideration. At the end of the derivation, Furukawa, et al. [14], was able to express the spring parameters in terms of the surface segment area, and the restoring springs between continuous elements are derived by considering an equivalent constant of springs in series. The derivation, however, only considered solid blocks. For CMUs, voids or grout infills affect the spring development if the method of Furukawa, et al. [14], is used directly. This study proposes an alternative approach to solving the spring constants for hollow concrete masonry units.
To solve for the spring constants for a hollow CMU, this study proposes to subdivide a CMU further into smaller elements and develop the springs for each sub-elements. Figure 3 (b) shows the proposed subdivision, where face shells are treated as single sub-elements, and where the webs and grout in-fills are also individually treated as rigid sub-elements. Applying the development of surface segment areas to the smaller CMU sub-elements, we can arrive at spring constants on the surfaces of each subelements. Figure 3 (c) shows the springs of one face shell and one void/grout in-fill sub-element.
The black springs shown in Figure 3 (c) is considered as internal springs for this study. Since one CMU is considered as rigid element, these black springs do not deform, and are not relevant for this study. In the following discussions, the vertical springs, which connects the rigid block to the block above or below it, and the horizontal springs, which connects the rigid block to the block to its left or right, for the web sub-elements are derived. Equation (5) is rewritten to express the equivalent spring constants per area in terms of transformed dimensions instead of the pyramid height. Specifically, for vertical (red) springs: Similarly, the equivalent spring constants per area of the horizontal (green) springs are as follows: where ‫ܧ‬ ǡ ‫ܧ‬ ǡ ‫ܧ‬ ெ -Young's modulus for elements A, B, and mortar, ߥ ǡ ߥ ǡ ߥ ெ -Poisson's ratio for elements A, B, and mortar, ‫ܪ‬ ǡ ‫ܪ‬ -height of homogenized elements A, and B, ܹ ǡ ܹ -width of homogenized elements A, and B, and ‫ݐ‬ ெ -mortar thickness.

Springs at web.
For the vertical (magenta) springs, it needs to be ensured that in each segment area, the material in one element (one constituent spring) is uniform to easily derive the equivalent spring using equation (5). This means that, for each of the block overlap, the surface sub-elements are made small enough to ensure material uniformity on each constituent spring of the springs in series (see Figure  5 (a)). Once material uniformity is ensured, the equivalent spring constants per area for vertical (magenta) springs for web elements take the form of equation (6).
For the horizontal (blue) springs, the equivalent spring constant per area is actually for more than three materials in series. Referring to Figure 5 (b), the web spring is actually a series of web springs and void/grout in-fill springs. For this study, it is assumed that the exterior voids ‫ݐ(‬ ଵ ǡ ‫ݐ‬ ଶ , see Figure 5 (b)) of the CMU are always filled with grout. Also, for this study, the contribution of the spring connection between the webs and face shells is ignored. The equivalent spring constant per area for the horizontal (blue) springs for web elements can be expressed as: where ‫ݐ‬ ௩ ൌ ͲǤͷͲ ݈ ௩ଶ ݈ ௩ଷ -total void length for element A, ‫ݐ‬ ௩ ൌ ݈ ௩ଵ ͲǤͷͲ ݈ ௩ଶ -total void length for element B, ‫ݐ‬ ௪ ൌ ݈ ௪ଷ ݈ ௪ସ -total web thickness for element A, and ‫ݐ‬ ௪ ൌ ݈ ௪ଵ ݈ ௪ଶ -total web thickness for element B.

Constitutive model for the equivalent springs.
The constitutive model of springs used for this study is based on the model proposed by Gu, et al. [9]. The normal spring is assumed to be elastic from compressive strength ߪ ௫ (at strain ߝ ), up to tensile strength ߪ ௧௫ (at strain ߝ ௧ ). When the compressive strength is reached, there is increase in strain without any increase in the load carried. On the other hand, when the tensile strength is reached, the line slope reverses, and the strain continues to increase until fracture is reached at ߝ ௧௨ . For the shear springs, for both direction (negative or positive), when the shear strength ߪ ௦௫ (at strain ߝ ௦ ) is reached, the strain continues to deform even without For this proposed method, wxAERO is used. wxAERO defines the springs in terms of forces and deformations. A force-deformation constitutive model based on the stress-strain model by Gu, et al. [9], is used for this proposed method (see Figure 6). The force-deformation model for CMU spring (red line) and mortar spring (green line) are shown in Figure 6 (a), and the equivalent spring model (blue line) is shown in Figure 6 (b). For the individual models shown in Figure 6 (a), the elastic range has slope equal to the total spring constants for elements A, B, and mortar, respectively: where, ‫ܭ‬ Ǥ ǡ ‫ܭ‬ Ǥ ǡ ‫ܭ‬ Ǥெ -total normal spring constant for elements A, B, and mortar, ‫ܭ‬ ௦Ǥ ǡ ‫ܭ‬ ௦Ǥ ǡ ‫ܭ‬ ௦Ǥெ -total shear spring constant for elements A, B, and mortar, ݇ Ǥ ǡ ݇ Ǥ ǡ ݇ Ǥெ -normal spring constant per area for elements A, B, and mortar, ݇ ௦Ǥ ǡ ݇ ௦Ǥ ǡ ݇ ௦Ǥெ -shear spring constant per area for elements A, B, and mortar, and ‫ܣ݀‬ -surface area for element segment in consideration.
For the equivalent spring, the slope of the elastic range is the equivalent total spring constant: Alternatively, equation (10) can be expressed in terms of the individual spring constant per area: where ݇ ത ǡ ݇ ത ௦ -equivalent spring constant per area for normal and shear spring. For the equivalent spring shown in Figure 6 (b), the elastic range limit of the force-deformation curve is set at the instant that the minimum of the strengths of the constituent springs is reached. The maximum strength of the equivalent spring is the minimum strength of the constituent springs: where ߪ Ǥ௫ ǡ ߪ ௧Ǥ௫ ǡ ߬ ௦Ǥ௫ -equivalent spring compressive, tensile, and shear strength, ߪ Ǥ ǡ ߪ Ǥ ǡ ߪ Ǥெ -compressive strength for elements A, B, and mortar, ߪ ௧Ǥ ǡ ߪ ௧Ǥ ǡ ߪ ௧Ǥெ -tensile strength for elements A, B, and mortar, and ߬ ௦Ǥ ǡ ߬ ௦Ǥெ ǡ ߬ ௦Ǥ -shear strength for elements A, B, and mortar.
From the maximum strengths, the maximum forces is solved: where, ݂ Ǥ௫ ǡ ݂ ௧Ǥ௫ ǡ ݂ ௦Ǥ௫ -maximum spring compressive, tensile, and shear force. Lastly, the deformation at the instance of maximum force is solved. At the said instance, the deformation is the sum of the deformation of individual springs. Alternatively, the total equivalent spring constants can be used in solving the total deformation. This is done by rearranging equation (1) into the following equations: where ‫ݑ‬ Ǥ௧௧ ǡ ‫ݑ‬ ௧Ǥ௧௧ ǡ ‫ݑ‬ ௦Ǥ௧௧ -total deformations at maximum compressive, tensile, and shear force, ݂ Ǥ௫ ǡ ݂ ௧Ǥ௫ ǡ ݂ ௦Ǥ௫ -maximum spring compressive, tensile, and shear force, and ‫ܭ‬ ഥ ǡ ‫ܭ‬ ഥ ௦ -equivalent total normal and shear spring constant.
When springs between rigid elements fail, the continuous material is converted into discontinuous one. The contact and friction forces need to be defined. For this study, these forces are generated by the virtual environment of wxAERO.

Blast load validation
Normal blast load, or the peak pressure load on the point of a wall nearest to the blast source, is checked by comparing the output load-time history from past experiments. In a work by Rigby & Sielicki [15], the TNT equivalence of a PE4 charge is investigated. In a series of experimental set-ups, the pressure histories of a normally reflected hemispherical PE4 explosion are measured by a pressure gauge installed perpendicular to the blast source.
Included in the Rigby & Sielicki's study [15] is the pressure-time history for a hemispherical ʹͷͲ ݃ PE4 explosion at a standoff distance of Ͷ ݉. Given the blast loading parameters, and considering the TNT equivalence of 1.20 for PE4 explosive, a load-time history is generated by xPLoDE (Explosive Pressure Load on Discrete Elements. The load histories from the experiments and xPLoDE are overlaid as shown in Figure 7, which shows close correspondence between the two curves. The study by Rigby, et al. [16], measured the blast pressure at a point in the wall directly normal to the blast source (G1), at points ʹ ݉ to the left (G2) and above (G3) of point G1, and at a point on the ground ͵ ݉ to the left of point G1 (G4). The peak reflected pressure, impulse, and duration of both the positive and negative phases of each test setup are measured, and tabulated. To validate the xPLoDEgenerated load-time history for shockwaves approaching a surface at an angle of incidence, the test cases are individually run using xPLoDE. Tan [13] tabulated the xPLoDE-computed values for peak pressures, impulse, and phase duration, and compared them to the test results. The blast pressure distribution generated by xPLoDE is plotted and compared to the test results (see Figure 8).

Response validation
The validity of the proposed tool, consisting of xPLoDE and MWBLOCk, for blast simulations is tested by replicating an experiment by Abou-Zeid, et al. [10].
The typical numerical model of the walls is shown in Figure 9 (a). Included in the numerical model is the wall made of unreinforced CMU (dark gray blocks), eleven courses high and two and half courses The frame, and the top and bottom support blocks are fixed in place, while the blocks for the wall are connected together with springs, with values as shown in Figure 9 (b) and (c). The blocks interact with the frames and the supports only through contact. The supports also prevent the wall from falling since the wall is subjected to gravity loads.
Due to some material parameters not given in the reference study [10], and due to some model parameters that has to be determined, a parametric study is performed for the numerical model to     11 replicate the actual test results as closely as possible. The full parametric study is presented in Tan [13]. As a result of the study, the recommended material properties are presented in Table 1, and the recommended model parameters are presented in Table 2. The resulting response of the numerical model is compared to the actual test measurements as show in Figure 10.

Conclusions
The blast analysis tool proposed in this study, is composed of two modules: xPLoDE (Explosive Pressure Load on Discrete Elements) and MWBLOCk (Masonry Wall Block Layer with On-Surface Spring Constant k). xPLoDE computes the blast loading on individual elements of a discretized structure using UFC 3-340-02. MWBLOCk generates a numerical model of a CMU wall as an assembly of rigid blocks and computes the normal and shear constants of the springs that connect individual rigid blocks.  Notes: (a) Abou-Zeil, et al. [10], (b) Seyedrezai [18], (c) Gu, et al. [9]. The blast load-time histories and peak pressure output by xPLoDE are validated by comparing them to measurements from experiments done by other studies. And the validity of the proposed tool is investigated by comparing the response computed by the proposed tool to the response measured from experiments.
It is demonstrated that xPLoDE gives a very close estimate of the blast loading on walls exposed to direct pressures (front wall pressures). The peak pressure distribution on a front wall surface, and the load-time history for a point in a wall is shown to be in close correspondence to experimental values.
It is also demonstrated that the maximum wall response can be adequately approximated by the proposed model. It is demonstrated that the proposed tool gives adequate approximation of the maximum wall response displacements, especially for one-way walls with arching mechanism.