Nonlocal Formulation for Numerical Analysis of Post-Blast Behavior of RC Columns

Residual axial capacity from numerical analysis was widely used as a critical indicator for damage assessment of reinforced concrete (RC) columns subjected to blast loads. However, the convergence of the numerical result was generally based on the displacement response, which might not necessarily generate the correct post-blast results in case that the strain softening behavior of concrete was considered. In this paper, two widely used concrete models are adopted for post-blast analysis of a RC column under blast loading, while the calculated results show a pathological mesh size dependence even though the displacement response is converged. As a consequence, a nonlocal integral formulation is implemented in a concrete damage model to ensure mesh size independent objectivity of the local and global responses. Two numerical examples, one to a RC column with strain softening response and the other one to a RC column with post-blast response, are conducted by the nonlocal damage model, and the results indicate that both the two cases obtain objective response in the post-peak stage.


Introduction
Reinforced concrete (RC) columns are important loadbearing elements to prevent progressive collapse of the reinforced concrete structures. However, in the design of the reinforced concrete structures, the extreme loads generated by terrorist attacks or accidental explosions were usually ignored that made the structures particularly vulnerable when they are subjected to these loading conditions. When a specific blast load is applied, RC columns usually suffer localized failure (spalling) and exhibit global deformation due to shear and flexural bending that depend on load characteristics and structural parameters (Crawford and Magallanes 2011). In addition, certain residual axial load carrying capacity may still remain, which provides the main basis for assessing of damage level and collapse risk for structures under blast loading (Brunesi et al. 2015;Parisi 2015;Petrone et al. 2016;Russo and Parisi 2016). Therefore, it is quite essential to develop a more precise numerical model for post-blast response of RC columns.
Currently, only a few experiments for post-blast residual axial capacity of RC columns have been reported. Li et al. (2012) obtained the residual axial capacity of limited seismic (LS) and non-seismic (NS) RC columns under quasi-static loads, while the blast damage is applied through the predicted residual lateral deflection calculated by a numerical model. Roller et al. (2013) took a deep investigation on residual axial capacity of conventional and hardened RC columns subjected to contact or close-range detonation loading, and the blast damage is mainly due to material losses. With great progress on numerical method in recent years, a lot of research work has been carried out to study the dynamic response of RC components and structures under blast loading (Shi et al. 2008;Bao and Li 2010;Jayasooriya et al. 2011;Park et al. 2014;Brunesi et al. 2015;Li et al. 2016;Lim et al. 2016). In most of these numerical studies, a numerical model was firstly developed and validated by the experimental results. Then, successive parametric studies were carried out with this validated numerical model. However, the convergence test of post-blast response is not included in these papers, as a consequence, once the element size changes, the accuracy of postblast residual capacity may be compromised.
As it is well known, the stress strain relation of concrete materials contains a descending branch, which is called strain softening, such that the finite element results always become mesh-dependent (Jirásek and Bažant 2002). In classical continuum mechanics, the introduction of strain softening generally leads to erroneous results. From a mathematical point of view, the governing partial differential equations lose hyperbolicity in dynamics and ellipticity in quasi-statics. Due to the lack of a length scale in the constitutive model, the localization zone approaches to zero volume where energy dissipation vanishes in the fracture process since a finite amount of dissipated energy per unit volume is defined through the strain softening law, which is unacceptable from a physical point of view. In the numerical sense, results show an extreme sensitivity to the spatial discretization, and deformation localizes into an infinitely narrow band upon mesh refinement.
This sensitivity of the mesh size on the numerical results has attracted enormous attention of researchers in the field of solid mechanics, and different solutions have been proposed to deal with the ill-posed problem (Bažant and Oh 1983;Pijaudier-Cabot and Bažant 1987;Bažant and Jirásek 2002;Armero and Ehrlich 2006). Among which, the crack band model is practically the simplest but a crude method to avoid the pathological sensitivity to mesh refinement by scaling of the constitutive law. And it has been incorporated into the constitutive model of concrete in some commercial finite element software, e.g., LS-DYNA (Hallquist 2007). Although this model ensures correct energy dissipation in the localized damage band, it is only appropriate for mode-I fracture. However, as pointed by Ožbolt et al. (2011), the failure mode tends to change from mode-I to mixed modes as the loading rate increases. In order to achieve a thoroughly regularization, the nonlocal theory is widely adopted in the analyses related to mesh sensitivity strain softening problems since the pioneering work by Bažant (Pijaudier-Cabot and Bažant 1987). The essential characteristic of the nonlocal model is the incorporation of a length scale, by which the constitutive law at a material point of a continuum depends on not only this point but also its certain neighborhood. In fact, concrete has a complicated internal microstructure, the change of which is crack bridging, and can be captured by the enriched continuum model through the characteristic length. Moreover, the embedded strong discontinuity approach developed in recent years successfully overcame the mesh-dependent problem by incorporating a localized dissipative mechanism through introduction of strong discontinuity in kinematics and defining a softening cohesive law (Armero and Ehrlich 2006;Jukić et al. 2014). However, it is not straightforward to implement into the general finite element codes such as LS-DYNA.
This paper aims to eliminate the pathological mesh size dependent problem of post-blast residual axial capacity for RC columns, which will result in a more precise dynamic simulation of RC structures subjected to blast loads. Firstly, a comprehensive analysis of the post-blast residual axial capacity of a RC column is conducted using the currently available concrete models in LS-DYNA that would demonstrate the deficiencies of these models. Then, a nonlocal integral formulation is incorporated into the Mazars damage model (Mazars and Pijaudier-Cabot 1989) to maintain the objectivity of post-peak responses for RC columns. Furthermore, the applications of two numerical examples demonstrate the feasibility of the nonlocal damage model in dealing with the pathological mesh size dependent responses involving concrete strain softening. It is worth noting that the dynamic response of structural components to blast loading strongly change with the type of explosion, and only the type of detonation loading is considered in this study.

Assessment of Post-Blast Responses in Considering Concrete Strain Softening
RC columns typically suffer certain damages when the RC structures are subjected to blast loading. Its subsequent response of the structure mainly depends on the post-blast residual axial capacity of critical RC columns. Therefore, it is of importance to get an appropriate post-blast residual axial capacity such that correct dynamic response of RC structures can be generated. Generally, a high-fidelity physics-based (HFPB) model was developed in LS-DYNA to do some research work related to blast analysis of RC columns as several concrete constitutive models in the material library are available and proved to be suitable for high strain rate and intensely nonlinear behavior of concrete (Shi et al. 2008;Bao and Li 2010;Park et al. 2014;Li et al. 2016), especially the Karagozian & Case concrete (KCC) model (*MAT_072R3) and the continuous smooth cap (CSC) model (*MAT_159) (Hallquist 2007). In this section, two constitutive models are firstly introduced. Then, comparative studies are carried out to investigate the influence of the different mesh sizes on the post-blast residual axial capacity of a RC column. Corresponding discussions are given based on these comparative studies.

Karagozian & Case Model
As known, the volumetric and deviatoric responses of concrete are decoupled in the KCC model, in which a user defined equation of state (EOS) defines volumetric strain versus pressure relationship and a movable failure surface limits the deviatoric stress. In order to govern the evolution of the failure surface, three fixed independent failure surfaces termed as yield, maximum and residual are defined. For hardening behavior, the failure surface is interpolated between the yield and maximum failure surface based on a damage variable, and for softening behavior, a similar interpolation is performed between the maximum and residual failure surfaces. The strain rate dependent behavior is considered through a user defined dynamic increase factor (DIF) curve (Hallquist 2007).
Three parameters are used in the KCC model to control the damage evolution, they are b 1 , b 2 and x for compression damage, tensile damage and volume expansion, respectively. It is a common sense that the numerical results will not be objective upon mesh refinement if the strain softening behavior of the constitutive model is not associated with a localization limiter or characteristic length. To eliminate the pathological mesh size dependent phenomenon, the crack band method is incorporated into this model by forcing the area under the stress-strain curve to be a constant value G f /h, where G f is the fracture energy and h the element characteristic size. G f can be obtained through experiment (Lee and Lopez 2014). However, this method is only valid for tensile softening and becomes invalid when the element size exceeds 250 mm or the element size is smaller than the localization width (Magallanes et al. 2010).

Continuous Smooth Cap Model
The CSC model is originally developed for use in roadside safety simulations, and is also frequently used to model concrete materials subjected to impact and blast loads. Its yield surface combines a shear surface with a hardening compaction surface smoothly and continuously through the multiplicative formulation, which avoids the numerical complexity of dealing with ''corner'' regions. As softening is a remarkable behavior of concrete in tension and low confinement compression, the strain-based damage formulations are incorporated into this model, and brittle damage and ductile damage are distinguished. In order to consider the strain rate effect for impact and blast analysis, the viscoplastic formulation is adopted. Details of the description of the model can be referred to the reports Murray 2007).

Validation of the Numerical Model
In all the HFPB simulations, concrete is modeled by 8-noded solid element with single integration point. The reinforcing steel is modeled with 2-noded beam element. And perfect bond between concrete and reinforcement is assumed by sharing nodes. The KCC or CSC model is used for concrete material and an elasto-plastic material (*MAT_PICEWISE_LINEAR_PLASTICITY) is used to model steel. The erosion technique is implemented to consider the severely damaged concrete elements.
In order to validate the reliability of the numerical model, two numerical analyses of the CTEST20 column (Crawford et al. 2012) and B40 beam (Magnusson et al. 2010a) under blast loading are conducted and the numerical results are compared with experimental data.

CTEST20 Column
The geometry and reinforcement details of the column are shown in Fig. 1. The unconfined compressive strength of concrete is 41.4 MPa. The yield stress and hardening modulus of reinforcement are 475 and 751 MPa, respectively. The convergence tests with 12.7, 25.4 and 50.8 mm mesh sizes show that the mesh size of 25.4 mm is appropriate for the numerical analysis.
The mid-height deflection time histories of the column are shown in Fig. 2a, as can be observed in this figure, the numerical model with both KCC and CSC models present reliable peak deflection as compared with the measured value. However, the predicted residual deflection by the numerical model is higher than the experimental result. This might be due to the erosion of concrete material at the two ends and the simplified modeling of blast load by adopting LOAD_-BLAST_ENHANCED in LS-DYNA (Hallquist 2007). Figure 2b shows the calculated final deformation state of the column by KCC and CSC models, there is a slight difference between the results although the same erosion criterion is used. The most possible reasons are that the automatic generated material parameters are adopted for CSC model such that the dynamic increase factor (DIF) might be different from that of KCC model. Moreover, different damage definitions are adopted in the two models. Since strong nonlinear behavior existed in the tested RC column under blast loading, the numerical predictions conducted in this section could be considered within an acceptable level.

B40 Beam
In the authors' opinion, the deviations between the simulated and tested results of CTEST20 column mainly generate from simplified modeling of blast load. As a result, a RC beam (type B40) tested by Magnusson et al. (2010a) that also presented the accurate blast load-time history, is simulated in this section. Some primary geometry parameters of the beam are shown in Fig. 3, and more detailed information can be found in the related references (Magnusson et al. 2010a, b). In the finite element model, the concrete compressive strength and reinforcement yield strength are defined as 43 and 604 MPa, respectively. And the blast load with reflecting pressure of 1250 kPa and impulse density of 6.38 kPa s is adopted in this study as reported by Magnusson et al. (2010a).
The dynamic response of the RC beam is simulated by both KCC and CSC models. As can be observed from Fig. 4 displacement is about 16.5 mm in both cases due to shear failure of the RC beam, which is obviously shown by the calculated damage distribution and final deformation mode in Fig. 5. And the shear failure mode is just the same as that reported by Magnusson et al. (2010a).

Post-Blast Response Analysis of A RC Column
The KCC and CSC models have been frequently used in dynamic response analysis of RC columns under blast loading, and it is generally considered to be converged when the displacement response does not vary a lot upon mesh refinement. However, there is limited information on the studies that considered the mesh size influence on the postblast residual axial capacity to date.
A specific RC column is taken as an example here to calculate the post-blast residual axial capacity by varying the element sizes. Figure 6 shows the scheme of the RC column, the clear height of which is 2750 mm. And the section width and height are both 300 mm with a concrete cover of 40 mm. A total of four U20 mm longitudinal rebar are uniformly distributed along the perimeter of the section, and the stirrup reinforcement with a diameter of 10 mm is spaced 200 mm along the length. The uniaxial compressive strength for the concrete is 30 MPa. The yield and ultimate strength of steel reinforcement are 335 and 480 MPa, respectively. The loading scheme consists of three stages, as shown in Fig. 6. In stage 1, an initial axial load of 821 kN is applied on the top end to produce an initial stress state in the column. After a steady state is reached, an idealized blast load with peak overpressure of 10 MPa and duration time of 1 ms is applied to the column along the full height and certain damages will be generated till the dynamic response is damped out in stage 2. Finally, the displacement control loading scheme is used to calculate the post-blast residual axial capacity in stage 3. Four different element sizes, i.e., 50, 25, 12.5 and 6.25 mm, are used in this study. Both KCC and CSC models are adopted for the concrete material. And the piecewise linear isotropic plasticity model is employed to model the steel reinforcement. The mid-span transverse displacement time history curves are described in Fig. 7 for different mesh sizes. As can be observed, for the results calculated by both KCC and CSC models, the 50 mm element size seems a little large for the results to achieve convergence, while the displacement time history curves are almost identical for cases with element sizes of 25, 12.5 and 6.25 mm before loading stage 3. In stage 3, the displacement curves diverge a lot, which indicates different post-blast response of the RC column for analysis results with different mesh sizes. The post-blast residual axial capacity can be easily obtained from the axial load time history curves as shown in Fig. 8. It can be observed that the residual axial capacity has a pathological dependence on the mesh size. This is mainly because that the concrete material models are not equipped with a complete regularization formulation, hence the damages generated in the local areas still have a mesh dependent character even if the constant fracture energy method is introduced. Figure 9 shows the damage distribution of the RC column at the end of stage 2, the local damages generated during blast loading are obviously very different, and leading to different residual axial capacity. It is no doubt that the different residual axial capacity will affect the dynamic response of RC structures under blast loading. Therefore, it is essential to take measures eliminating the mesh dependent issue in order to ensure reliable analysis results for RC structures under blast loading.

Proposed Numerical Model
As described in Sect. 2, the KCC and CSC models with the constant fracture energy method still fail to obtain objective results of post-blast residual axial capacity. To solve this issue, an enrichment to the Mazars damage model is proposed by means of introduction of a length scale, in which the local equivalent strain is replaced by the nonlocal equivalent strain through a nonlocal formulation. And the enriched Mazars damage model is incorporated into the general finite element code LS-DYNA so that the numerical study of post-blast response for RC columns and even RC structures subjected to blast loads can be quite convenient. Because the three dimensional concrete constitutive model with nonlocal formulation is computational demanding even for RC components, the uniaxial form of the Mazars damage model is adopted for the computations of the mesh-dependent problem in this study.

Description of Mazars Damage Model
Damage mechanics has been proved to be efficient for modelling concrete behavior in recent years (Heo and Kunnath 2013; Ren et al. 2015). Mazars damage model belongs to one of the first models based on this framework and has been widely used in simulating the damage behavior of concrete (Mazars and Pijaudier-Cabot 1989). It is described in Fig. 10 by the relation where r is the stress tensor, e is the strain tensor, C is the stiffness matrix of material and D is the scalar damage which varies from 0 to 1 with 0 represents the intact material and 1 for fully failure of the material. And the damage criterion inspired from the St. Venant criterion of maximum principal strain is adopted, its equation is: where j represents the current damage value which is the maximum value of equivalent strain ever reached in time history, and the initial value of j is j 0 , which indicates the initiation of damage.
As extensions play a significant role in the damage behavior of concrete, Mazars (1986) suggested the equivalent strain evaluating the local intensity of extensions, which is defined as: where e i are the principal strains and Á h i is the Macaulay bracket.
In order to describe the different behaviors of tension and compression, two separated damage variables denoted as D t and D c are defined, and the total scalar damage D is calculated from the weighted combination of them, namely: With where, A i and B i are characteristic parameters of the concrete material identified from test data. In Eq. (4), the weighting coefficients a t and a c are functions of the strain state (a t ¼ 0 for pure compression and a c ¼ 0 for pure tension, a t and a c vary from 0 to 1 in other cases).

Nonlocal Formulation
The post-blast residual axial capacity of RC columns is closely related to the blast damages of concrete material. As stated in Sect. 2, the blast damages show pathological sensitivity to the mesh size, which directly results in the mesh-dependent residual axial capacity. The integral type nonlocal model provides a fully regularization for the strain softening problem through introducing a characteristic length in the local constitutive laws (Pijaudier-Cabot and Bažant 1987), which is achieved by weighted spatial averaging of a suitable state variable. Generally, state variables such as the damage energy release rate, damage value, total strain and equivalent strain are chosen for applying the nonlocal formulation. In this paper, the equivalent strain in Mazars damage model is taken as the nonlocal state variable, and with its nonlocal counterpart replacing e eq in Eq. (5), which is defined as where e andẽ are the nonlocal and local equivalent strains, respectively. V is the spatial domain of interest, which is decided by the characteristic length. And aðx; nÞ is the nonlocal weight function that defines the interaction coefficients between point x and n. The weight function is typically normalized by R V aðx; nÞ dn ¼ 1 for all x, then the normalized weight function can be obtained: where a 0 ðsÞ is a function which depends on the distance s between point x and n. In this study, the Bell-type weight function is adopted (Jirásek and Bazant 2002): where, R is the characteristic length that limits the interaction zone, which makes the function exactly zero for s ! R.

Numerical Implementation
LS-DYNA is a general finite element code which provides abundant interfaces for users developing their own subroutines to achieve research objectives (Hallquist 2007). In this study, before invoking the local constitutive model, the local equivalent strain must be replaced by the nonlocal counterpart of the current time step. The scheme of the beam element formulation is shown in Fig. 11, only one integration point is used for one element to save computational time in explicit analysis. The values of the nonlocal equivalent strain must be traced at every individual integration point, whose coordinates are denoted as x k (k ¼ 1; 2; . . .; N; N represents the total number of integration points). The numerical algorithm can be determined as with where w l is the integration weight of the integration point number l, and the value of which is 1 as single integration point is used for all the beam elements. It is worth noting that a kl would be zero when the distance between k and l is larger than the interaction radius R (characteristic length). It could be calculated only once before the first step and then stored as it does not vary during the simulation. Figure 12 shows the numerical procedure of the nonlocal formulation. The transformation matrix from local to nonlocal equivalent strain is defined as b, in which b kl ¼ w l a kl . Then, before invoking the local constitutive model at every new time step, the local equivalent strain is transformed to the nonlocal one, and stresses of each integration point are updated, the corresponding history variables are stored for next step.

Numerical Examples
In this section, the numerical model of RC columns is developed using Hughes-Liu beam element (Hallquist 2007). The cross section is discretized into concrete and reinforcing steel fibers, and perfect bond between concrete and steel is assumed. The severe damage of concrete or steel can be considerer by the definition of ''failure strain'' in the constitutive models.

Softening Response of A Tested RC Column
One of the tests on RC columns conducted by Tanaka and Park (Tanaka 1990), labeled specimen 7, is analyzed in this example due to the considerable axial load applied, which makes the overall behavior entering post-peak region. The geometry and section details are shown in Fig. 13. The unconfined and confined compressive strength of the concrete are f 0 c ¼ 32 MPa and f c ¼ 39 MPa, respectively. A bilinear stress-strain relationship is used for the reinforcement with yield strength of 510 MPa and 1% strain hardening ratio.
In order to have a deep insight into the influence of concrete strain softening behavior on post-peak response of the RC column with mesh refinement, the analysis is carried out with 10, 20 and 40 elements and both the local and nonlocal damage models are adopted for concrete. The main material parameters for concrete are A c ¼ 0:501, B c ¼ 707:1, A t ¼ 0:0005 and R = 1100. The displacement control scheme is applied on the top node of the column. The computed forcedisplacement curves are shown in Fig. 14. As can be seen in the figure, the post-peak stage of the curves become more and more brittle as the number of elements increases for the local damage model. This is mainly because that the inelastic deformation localizes in one element, resulting in the increasing strains in the extreme element as the column enters the post-peak stage. However, the post-peak results computed by the nonlocal damage model are almost the same with mesh refinement, and match well with the envelope of the cyclic response from the experiment.  Fig. 13 Configuration of the tested RC column.

Recalculation of the RC Column in Sect. 2.4
In Sect. 2.4, the post-blast residual axial capacity of the RC column is pathologically sensitive to the mesh size as the constant fracture energy method used in the concrete constitutive models fails to capture objective damage values under blast loading. In this section, the RC column is used as an example to demonstrate the feasibility of the nonlocal damage model in calculating the post-blast residual axial capacity. Finite element models with 20, 40 and 80 elements are analyzed with the same load in Sect. 2.4. The main material parameters for concrete are A c ¼ 1:5154, B c ¼ 1767:77, A t ¼ 0:1 and R = 600. And the strain rate effect of concrete is accounted by the dynamic increase factor (DIF) calculated from Malvar (Malvar and Ross 1998).
The computed axial force versus time history curves are presented in Fig. 15. It is obvious that the curve calculated by the finite element model with 20 elements shows a minor difference with the rest two curves, which means that the finite element model with 20 elements does not guarantee the convergence of the numerical results. This is consistent with the consensus that a sufficiently fine mesh is generally required for nonlocal model. The results computed from the finite element models with 40 and 80 elements are almost the same, indicating that the proposed numerical model can not only capture identical damage values during blast loading but also ensure objective post-blast residual axial capacity of the RC column.
Furthermore, the finite element model with 40 elements is herein taken for a series of calculations by varying the characteristic length R in the nonlocal model. According to Bažant (Pijaudier-Cabot and Bažant 1987), the material characteristic length cannot be less than the beam depth h, which is 300 mm in this example. As can be observed in Fig. 16, when the characteristic length R equals to 200 mm, the RC column loses its loading capacity under blast loading. And it is obvious that adopting a smaller characteristic length will lead to underestimate the dissipated energy within the damage zones as well as the post-blast residual axial capacity. The proper calibration of the nonlocal internal length is not a trivial task. Although this quantity is usually related to the theoretical width of the fracture process zone (material characteristic length), its actual definition is not completely clear. In this study, the simulation results seem to match the experiment well when twice of the column height is used for R. However, some specifically-designated experiments are still needed to calibrate this important parameter in the future.

Conclusions
Although the KCC and CSC models are equipped with the constant fracture energy method to consider strain softening behavior of concrete, the calculated post-blast residual axial capacity of a RC column still shows strong sensitivity to spatial discretization because the constant fracture energy method in the concrete constitutive models fails to capture objective damage values, which may significantly affect the post-blast residual axial capacity of RC columns. As a consequence, the nonlocal formulation is incorporated into Mazars damage model to overcome the problem. The validity of the nonlocal model is shown through two numerical examples, one to a tested RC column with strain softening response and the other to a RC column with postblast response. And both the two cases obtain objective results with mesh refinement. It is worth mentioning that the characteristic length R in the nonlocal model plays a dominant role in accounting for the microstructure of the concrete material, which could significantly affect the numerical results and needs to be calibrated by experiment in the future. And the presented numerical model should be further validated against other experimental results, exploring other types and magnitudes of blast loading. Nevertheless, it is still potential for analysis of RC structures when softening or post-blast response is expected.