A Simplified Method to Predict Damage of Axially-Loaded Circular RC Columns Under Lateral Impact Loading

Detailed finite element (FE) models are often employed to predict the impact responses of reinforced concrete (RC) columns. However, they always require substantial investments of time and effort in modeling and analysis so that they are not widely used in practice, particularly in preliminary designs. Moreover, although some simplified models have been established for beams and slabs under impact loading, few attempts have been made on modeling RC columns. For these reasons, this paper proposes a simplified modeling method to accurately capture the impactinduced response and damage of circular RC columns. In the proposed method, a two-degree-of-freedom (DOF) system was used to describe the interaction between the impactor and the impacted column. The formulas, and procedure to estimate the force–deformation relationship with strain-rate effects were presented according to the section-based analysis. The influence of the unloading stiffness on the residual deformation was addressed, and the method to determine the unloading stiffness of circular columns was proposed. Furthermore, a fiber-based beamcolumn element modeling method was developed to estimate the force–deformation relationship of the columns with strain-rate effects. The proposed simplified method was demonstrated by the drop-hammer impact tests to be capable of predicting the impact response of RC columns well. Its accuracy in the residual deformation is superior to that of the detailed FE simulation. Parametric studies were performed to investigate the damage characteristics of axially-loaded circular RC columns under various impact scenarios.


Introduction
Circular reinforced concrete (RC) columns are widely used in bridge and building structures as primary members to carry axial loads. In addition to service loads (AASHTO LRFD Bridge Design Specifications 2015; Deng et al. 2019), these circular columns should be designed to resist extreme loads such as lateral impacts from vehicle, vessel, falling rock, and so on (Davidson et al. 2012;Sha and Hao 2013;Liu et al. 2017;Do et al. 2018;Gholipour et al. 2018;Fan et al. 2015Fan et al. , 2016Fan et al. , 2018Fan et al. , 2019. A high number of column collapse accidents owing to impact loads have been documented around the world (Harik et al. 1990;Wardhana and Hadipriono 2003;Buth et al. 2010;Fan 2012). Hence, it is crucial to accurately predict the impact-induced responses of circular RC columns in the design of infrastructure.
Flexural-dominated and shear-dominated failures in circular RC columns may be caused by impact loads (Fujikake et al. 2009). Since shear-dominated failures are brittle (non-ductile), they should be avoided by some measures (e.g., increasing transverse reinforcement ratios and the size of cross-section) in the design process. In contrast to shear-dominated failures, a certain level of flexural-dominated damages caused by impact loading is acceptable considering the low frequency of occurrence

Open Access
International Journal of Concrete Structures and Materials and the high intensity of loading. Therefore, it is crucial to reasonably evaluate flexural damage of circular RC columns under impact loading during the performancebased design of infrastructure. The residual axial loadcarrying capacity has been widely used to identify the damage levels of RC columns after impact loading. The previous study (Fan et al. 2019) showed that the residual axial load-carrying capacity is primarily dependent on residual deformation rather than the maximum deformation. Hence, it is essential to accurately predict the residual deformation in impact analysis.
With the maturity of nonlinear finite element (FE) codes (e.g., ABAQUS 2009 and LS-DYNA 2014), detailed finite element (FE) modeling methods were widely employed to investigate the performance and damage of RC members under impact loading in the past two decades (Gholipour et al. 2018;Consolazio and Cowan 2003;Thilakarathna et al. 2010;Yuan and Harik 2010;Fan and Yuan 2014;Abdelkarim and ElGawady 2017;Do et al. 2018;Cao et al. 2019;Chen et al. 2016). In this approach, three-dimensional solid elements were used to model concrete, while beam or truss elements were used to simulate reinforcing steel bars. Contact algorithms were defined to capture the interaction between impactors and RC members. However, since the method requires substantial investments of time, efforts, and computational resources in modeling and analysis, it was seldom used in practice, particularly in the preliminary design iterations (Consolazio and Cowan 2005;Fan et al. 2011Fan et al. , 2019Wang and Morgenthal 2018). More importantly, because of the limitations of material models in general-purpose contact-impact nonlinear FE codes, the accuracy of this approach was shown to be barely satisfactory for some key impact-responses such as the residual deformation of the impacted RC columns. Liu et al. (2017) indicated that the typical detailed FE models had low accuracy in the prediction of the residual deformation of axially-loaded circular RC columns after impact loading. Although an improved FE modeling method was given in (Liu et al. 2017) to improve the accuracy, its coarse approximations and assumptions (e.g., specifying different RECOV parameters in the concrete material model to consider the influence of crack closure) would place some difficulties in practical applications.
In addition to the detailed FE modeling methods, some empirical formulas were derived to estimate the residual deformation of RC beams after impacts based on the experimental data from drop-hammer impact tests (Kishi and Mikami 2012;Kishi et al. 2002;Tachibana et al. 2010;Adhikary et al. 2016). These empirical formulas were based on the input energy, the static load-carrying capacity, and the maximum deformation. The influences of impact mass and velocities on the residual deformation were omitted. However, the physical experiments conducted by Zhao (2017) indicated that the impact-induced deformation of RC beams increased with the increase of the impactor's mass for the same impact energy. On the other hand, the use of these empirical formulas derived from RC beams for axially-loaded RC columns is also questionable. In the context, some simplified models with a single-or two-degree-of-freedom (DOF) system were developed to estimate the impact-induced responses of RC beams and plates with a rectangular cross-section (Fujikake et al. 2009;Fan and Yuan 2014;Zhao et al. 2018;Bertrand et al. 2015;Lee and Kwak 2018). Although similar in some aspects, the circular RC columns have some unique characteristics such as the reinforcement layout in circular cross-section and the presence of axial load. Hence, the applicability of the simplified models developed for beams and plates needs to be examined for RC columns. More importantly, few simplified models placed emphasis on the accurate prediction of the residual deformation. The residual deformation predicted by the simplified model often had lower accuracy than the maximum one (Fujikake et al. 2009). As mentioned above, the residual deformation is critical for the damage evaluation of RC columns under impact loading.
This study aims to develop a simplified modeling method to accurately predict both the peak and residual responses of axially-loaded circular RC columns under impact loading. In the developed model, a two-DOF system was employed to capture the interaction between the impactor and the impacted column. Two different methods were developed to estimate the force-deformation relationship with strain-rate effects based on the section-based analysis and the fiber-based element model, respectively. Different approaches to determine the unloading stiffness of the impacted columns were discussed for accurately predicting the residual deformation. The drop-hammer impact tests on axially-loaded circular RC columns were used to validate the proposed model. Parametric studies based on the simplified model were also performed to investigate the damage characteristics of axially-loaded circular RC columns under various impact scenarios. The iso-damage curves of the impacted columns were presented to examine the damage characteristics.

Development of the Simplified Model
The simplified model is developed in this section to predict the response of axially-loaded circular RC columns under impact loading. First, the dynamic equation of the simplified model is presented. Subsequently, the estimations of structural resistances are discussed in detail. Figure 1 presents the simplified model that is used to estimate the responses of axially-loaded RC columns under impact loading. Like (Fujikake et al. 2009), a two-degreeof-freedom (DOF) mass-spring system was employed to represent the impacting object and the impacted column. As illustrated in Fig. 1, the mass (m 1 ) of the impacting object is assumed to link the equivalent mass (m 2 ) of the impacted column by compression-only macro-elements that describe the contact stiffness (k 1 ) and damping ( c 1 ). The mass m 2 is connected to ground by discrete spring and damping elements that represent the structural resistance during impact loading. The equations governing the motions of the 2-DOF system can be written as where u 1 , u 1 and ü 1 are the displacement, velocity, and acceleration of the impactor, respectively; u 2 , u 2 , and ü 2 are the displacement, velocity, and acceleration of the equivalent mass m 2 , respectively; c 2 is the damping of the impacted column, which can be assumed to be omitted because the influence of damping on the response to a single pulse excitation is usually not significant unless it is highly damped (Chopra 2007;Fan and Yuan 2012;Fan et al. 2016); R s is the equivalent resistance of the axiallyloaded column. Equation (1) can also be expressed in the matrix form (after omitting the damping c 2 ):

Dynamic Equation of the Simplified Model
where k 2 (t) = R s u 2 , signu 2 /u 2 , which is equivalent stiffness corresponding to the force-deformation relationship of the axially-loaded column. The solution of Eq. (2) can be readily obtained by using the step-by-step numerical integration when the mass and stiffness matrices are known. The following sections will present the details on determining the parameters required for the solution of Eq. (2).

Load-Deformation Relationship Derived from Section Analysis
It is essential to reasonably estimate the load-deformation relationship of an RC column in the solution of Eq.
(2). Much research has been devoted to determining the load-deformation relationship of rectangular RC beams and plates (Fujikake et al. 2009;Carta and Stochino 2013). However, few studies were conducted for axially-loaded RC columns with a widely-used circular cross-section under impact loading. Some additional considerations are required for circular RC columns due to the circular arrangement of longitudinal reinforcement, which is different from that of rectangular members.

Material Models of Concrete and Reinforcing Steel
The concrete cover was modeled by the unconfined material model, and the core concrete was simulated by the confined material model proposed by Mander et al. (1988). The stress-strain relationship (Fig. 2a) of the confined material model can be written as (Mander et al. 1988) where σ c is the compressive stress of concrete; ε c is the compressive strain of concrete; f ′ co is the uniaxial compressive strength of the unconfined concrete; E c is the elastic modulus of concrete; ε c0 is the strain corresponding to the ultimate compressive strength (typically, ε c0 = 0.002 ); f co is the ultimate compressive strength of the concrete. f co = f ′ co for the unconfined concrete, while f co = f ′ cc for the confined concrete, as shown in Fig. 2a. The typical elastic-plastic relationship (Fig. 2b) was employed to model reinforcing steel bars, which can be expressed as Force Deformation Force Deformation Fig. 1 Simplified Two-DOF model of axially-loaded circular RC columns Page 4 of 24 Liu et al. Int J Concr Struct Mater (2020) 14:32 where σ s is the reinforcement stress; ε s is the reinforcement strain; ε y is the yield strain of reinforcing steel; f y is the yield strength of reinforcing steel; E s is the elastic modulus of reinforcing steel; E b is the secant modulus of the reinforcing bars after yielding. The dynamic behaviors of concrete and reinforcing steel bars are sensitive to strain-rate effects when RC members are subjected to impact and impulsive loads. Therefore, in addition to the constitutive model of the confined concrete given in Eqs. (3) and (4), Mander et al. (1988) presented the strain-rate effects on the dynamic behavior of confined concrete. For consistency, the following equations suggested by Mander et al. (1988) were employed in this study: where ε is the strain rate; f co,d is the compressive strength corresponding to the strain rate ε ; E c,d is the dynamic modulus corresponding to the strain rate ε . For reinforcing bars, the following equation recommended by the Japan Society of Civil Engineers (JSCE) (1993) was adopted: where f y,d are the static and dynamic yield strength of reinforcing bars, respectively. Figure 3 shows the moment-curvature relationships of a typical rectangular RC beam section and circular column sections, which were calculated with the software XTRACT (Chadwell and Imbsen 2004). The stress states of the longitudinal reinforcement and the neutral axis corresponding to the yielding moment ( M y ) in the circular column sections were largely dependent on the axial load ratio and the reinforcement layout. The circular geometry and the corresponding reinforcement layout complicates the way to determine the stress distribution in comparison with the rectangular sections. Importantly, similar to the rectangular cross-section, the moment-curvature curves of the circular cross-sections are bi-linearized as shown in Fig. 3. In static analysis, the moment-curvature curves can be readily determined by the use of the existing software. Different from the static analysis, the strain-rate effects should be taken into account by using Eqs. (6) to (8) in impact analysis. The explicit expressions for predicting the moment-curvature relations of the circular column sections are required to include the influence of strain-rate effects on the material properties readily. Generally, the preliminary yield and failure states can be estimated by the moment-curvature analysis using the software (e.g., XTRACT) when omitting the strainrate effects. Based on these states, the expressions of the moment-curvature curves can be written in an explicit form. Subsequently, the moment-curvature relationship considering the strain-rate effects can be obtained by incorporating Eqs. (6) to (8) into these explicit expressions. Assuming that the equivalent yielding state corresponds to the yielding of the longitudinal reinforcements in the j-th layer (see Fig. 4), and the neutral axis lies between the k-th layer and the k + 1th layer of the longitudinal reinforcements, the location of the neutral axis x y and the equivalent yielding moment M y can be determined based on the equilibrium equations in the cross-section. It is very complicated to determine the layer number ( j ) of the yielding reinforcements corresponding to the yielding state by iterative analysis. For simplification, the value of j can be determined from the preliminary cross-section analysis using the software XTRACT. The trial analysis showed that the strain-rate Page 5 of 24 Liu et al. Int J Concr Struct Mater (2020) 14:32 effects have a limited influence on the value of j , indicating the applicability of the simple method. In addition, computational efficiency can be improved as a result of this treatment. The axial load ( N c,y ) supported by the concrete can be written as where ε y,j is the yielding strain of the longitudinal rebars in the j-th layer; d e is the effective height of the RC column; D is the diameter of the RC column; s j is the distance from the j-th rebars to the rebars at the bottom, as illustrated in Fig. 4; y is the coordinate value along the (9)

Moment-Curvature Analysis of Circular Column Section
Dy − y 2 dy circular cross-section, as shown in Fig. 4. The axial load carried by the reinforcing bars is (10) Moment-curvature diagram of RC member: a rectangular RC beam; b circular column with axial load ratio of 7%; c circular column with an axial load ratio of 14%; d circular column with an axial load ratio of 28% Page 6 of 24 Liu et al. Int J Concr Struct Mater (2020) 14:32 where A s,i is the area of the reinforcing bar in the ith layer. When the axial load ( N ) is given, the neutral axis can be estimated from the following equilibrium equation.
Similarly, the bending moment M c,y carried by the concrete can be expressed as The bending moment M s,y resisted by the reinforcing bar is (11) N = N c,y + N s,y Based on Eqs. (12) and (13), the yielding moment M y can be determined by According to the plane cross-section assumption, the corresponding curvature θ y can be derived by The location of the neutral axis x u and the moment M u corresponding to the ultimate (or failure) state can be obtained based on the sectional equilibrium. The ultimate state is defined as the loading point with spalling of cover concrete, rupture of longitudinal rebars, and crushing of core concrete. With the rupture of the longitudinal rebars, the axial loads N c,u and N s,u as well as the bending moments M c,u and M s,u carried by the concrete and the reinforcement, respectively, can be written as  Liu et al. Int J Concr Struct Mater (2020) 14:32 Consequently, N c,u , N s,u , M c,u and M s,u corresponding to the crushing of core concrete can be expressed as (20) Correspondingly, the ultimate curvature θ u can be written as Once the key points θ y , M y and (θ u , M u ) are obtained, the equivalent bilinear moment-curvature relationship can be determined.
In order to account for the strain rate effects, the relationship between the deflection rate ( u 2 ) of the RC column and the curvature rate ( θ ) of the section is required. In the elastic range, the relationship between the midspan deflection ( u 2,el ) and the curvature ( θ e ) of the fixedend column can be given as where EI = M y /θ y and k = N θ y /M y . In the plastic range, concentrated plastic hinges would be formed at the impacted section and the ends, as illustrated in Fig. 5. In this case, the curvature ( θ p ) can be estimated by where l p is the length of the plastic hinge (usually, l p ≈ d e ). Based on Eqs. (27) and (28), the curvature rate θ at the mid-span section of the impacted RC column can be estimated by Similarly, when the axial load ( N ) is given, x u can be predicted based on the equilibrium equation as follows: The ultimate moment M u can be determined based on Eqs. (18) and (19) or Eqs. (22) and (23) as follows: By using Eq. (29), the strain rates of concrete ( ε c ) and tensile ( ε st,i ) and compressive ( ε sc,i ) rebars are calculated through the following formulae: (30) ε c = θ · x y for the yield poinṫ θ · x u for the ultimate state Page 8 of 24 Liu et al. Int J Concr Struct Mater (2020) 14:32 where d di is the distance from the i-th layer steel bars to the outermost layer bars in tension; d upi is the distance from the i-th layer steel bars to the outermost layer bars in compression. The material properties with the strainrate effects can be obtained by using Eqs. (30) to (32) and Eqs.
(3) to (8) when predicting the moment-curvature relationship of a circular column section. Similar to the reinforcing steel bars, the equation suggested by Fujikake et al. (2009) can be used to consider the influence of the location of concrete fiber on the strain-rate effect. Numerical analyses indicated that the location of concrete fiber has a limited influence on the impact response. This is primarily attributed to the fact that the dynamic increase factor of concrete is not very sensitive to the change of the strain rate in low-velocity impacts because of the functional relationships between the dynamic increase factors and the strain rates of concrete (see Eq. (6)). For this reason, Eq. (30) recommend by Carta and Stochino (2013) was used in this study to reduce the analysis time.

Load-Deformation Relationship
Since the yielding point corresponds to the translation from linear to nonlinear behaviors, the yield load P y is readily estimated from the yield bending moment M y by means of the linear elastic theory of columns: Similarly, the ultimate load P u can also be determined from the ultimate bending moment M u by means of equilibrium conditions. Since the P-delta effects would become dominant with increasing the impact-induced displacement, there are two different cases for the determination of the ultimate load P u . Generally, the ultimate load P u can be expressed as where when 2M u > Nu 2,u , P u = 8(M u − 0.5Nu 2u )/l and u 2,m = u 2,u . On the contrary, when 2M u ≤ Nu 2,u , P u = 0 . In the case, the maximum displacement u 2,m corresponding to the ultimate load P u can be determined by where During impact loading, the velocity u 2 of the impacted section varies with time. According to Eq. (29), the corresponding curvature rate θ and strainrate ε changes with time. Hence, the force-displacement relationship should be updated step-by-step to account for the influence of strain-rate effects in the solution of Eq. (1). Similar to the previous studies (Consolazio and Cowan2005;Fan et al. 2016), the central difference method was employed in this study to solve the dynamic equation of the two-DOF system for simplification. Accordingly, when the vertical displacement u 2,I and the corresponding velocity u 2,I are derived at the I-th time step, the force-displacement relationship can be updated as follows: 1. Determining θ I and θ I based on the u 2 − θ and u 2 −θ relationships; 2. Using an iterative routine to estimate x y,I and x u,I considering the influence of strain-rate effects; into Eqs. (30) to (32) to estimate the strain rates of concrete and steel reinforcement (where, x (0) y,I is the initial value estimated from the static section analysis without strain-rate effects when J = 1); b) Using Eqs. (6) to (8)  ui for the j-th iteration by solving Eqs. (11) and (24); d) Checking if x y,I and x u,I reach the results of convergence as follows: where ξ is the tolerance, which is set to be 10 −8 in this study. e) Obtaining the convergent x y,I and x u,I when Eq. (40) is met; otherwise, Step a) to d) need to repeat until getting the results of convergence.
3. Using Eqs. (14) and (25) to estimate M y and M u with the influence of strain-rate effects based on the obtained x y,I and x u,I , respectively; 4. Using Eqs. (15) and (26) to determine θ y and θ u corresponding to M y and M u , respectively; 5. Obtaining the updated force-deformation relationship from the updated moment-curvature relationship based on Eqs. (33) and (36).
During the numerical solution of the simplified two-DOF system, the column resistance can be determined for a given displacement at the I-th time step. As a part of the numerical solution, the above procedure is summarized in Fig. 6 to update the force-deformation relationship of the impacted column.

Contact Stiffness and Damping
When the ratio of the contact stiffness ( k 1 ) to the structural stiffness ( k 2 ) ranged from 0.5 to 200, the contact stiffness ( k 1 ) primarily affected the peak impact force but had a limited influence on the maximum displacement of the RC members (Fujikake et al. 2006(Fujikake et al. , 2009. The peak impact force is primarily related to the shear damage of the RC members rather than the flexural damage. Thus, considering the main objective of this study, the approximate method suggested by (Fujikake et al. 2006(Fujikake et al. , 2009 was employed to estimate the contact stiffness. According to (Fujikake et al. 2006(Fujikake et al. , 2009, the contact force-local deformation relation is assumed to be linear elastic for simplicity, and trial analyses using the simplified model can be performed with various contact stiffness for a typical impact case. The contact stiffness ( k 1 ) can be determined based on the magnitude of the impact load measured or estimated from the empirical contact model (Fujikake et al. 2006(Fujikake et al. , 2009). As mentioned above, the damping of the impacted members ( c 2 ) can be omitted for simplification. In contrast, the contact damping between the impactor and the RC members should be included to consider the energy dissipation in the local contact zone (e.g., local damage and friction between the impactor and the impacted member) (Fan and Yuan 2014). According to the recommendation given by (Fujikake et al. 2006(Fujikake et al. , 2009, the contact damping ( c 1 ) can be assumed to be one-half of critical damping, which can be written as Biggs (1964) pointed out that the kinetic energy of the simplified two-DOF system is required to equal that of the real system. Accordingly, the equivalent mass of the impacted column can be determined by where m c (z) is the distributed mass along the column length; φ(z) is the shape function of the impacted column. Biggs (1964) suggested that the deformation shape from the static application of impact loads can be used to define the shape function. Actually, the equivalent mass should be various at different impact phases (e.g., initial contact, loading phase, unloading phase Fan et al. 2019), and different deformation stages (e.g., elastic, and elastic-plastic). The equivalent mass of the impacted columns would be changed significantly at the initial contact phase, which would primarily affect the shear response. It is beyond the scope of this study and worth exploring in the future. For the flexural responses emphasized

Equivalent Mass of the Impacted Column
Page 10 of 24 Liu et al. Int J Concr Struct Mater (2020) 14:32 in this study, the equivalent mass is primarily dependent upon the behavior of the impacted column at the loading phase. In this case, the equivalent mass of the impacted column can be estimated for two typical shape functions (i.e., the simplified deformation shape shown in Fig. 5, and the static application of an impact load). Tiny differences of the equivalent mass obtained from two assumed shape functions were observed for the regular columns. Hence, the simplified deformation shape (as shown in Fig. 5) was suggested to estimate the equivalent mass of the impacted column for convenience. Figure 6 summarizes the calculation flowchart of the simplified two-DOF system for the impact responses of circular RC columns. The software MATLAB was employed in this study to implement the numerical algorithm of the two-DOF system. The proposed simplified model was validated by the drop-hammer impact tests on axially-loaded circular RC columns in this section. The impact tests on RC columns were briefly introduced. The approaches to determine the unloading stiffness were discussed in detail.

Impact Tests on Axially-Loaded Circular RC Columns
Compared with the impact tests on rectangular RC beams, few impact tests had been carried out on axiallyloaded circular RC columns. In this context, Liu et al. (2017) recently performed the drop-hammer impact tests on circular RC columns. One of the particular characteristics was the axial loads were exerted simultaneously in the impact tests, as shown in Fig. 7. Table 1 presents the information about the column specimens and the experimental results. Primary varied parameters included impact energy (low and high) with different drop heights, axial load ratio (0%, 14%, and 28%), and reinforcement ratio ( ρ L = 1.92%, ρ T = 1.3% and ρ L = 0.96%, ρ T = 0.72% ). Figure 7b presents the two different reinforcement layouts of the test columns. More details about the impact tests can be found in (Liu et al. 2017). All impact cases reported in (Liu et al. 2017) were used to validate the rationality of the developed simplified two-DOF model. The material properties of concrete and steel bars derived from the material tests were employed to define the parameters of the column in the simplified two-DOF model. The parameters (e.g., the force-deformation relationship of the column) were Initial conditions V 0 , m 1 , m 2 , k 1 , c 1 , c 2 , and material properties of RC column Determining the initial value of k 2 using Eqs. (9) to (14) and (33) to ( (11) and (24) Obtaining the convergent and Using Eqs. (14) and (25)   determined by the above-suggested approaches. In addition to the force-deformation relationship of the column (e.g., as shown in Fig. 8), the equivalent mass of the column was 36.60 kg and the value of k 1 was 2.4 × 10 8 N/m according to the methods provided in (Fujikake et al. 2006(Fujikake et al. , 2009.

Unloading Behaviors and Discussions
The previous experimental studies (Liu et al. 2017;Fan et al. 2019) indicated that the interactions between the impactor and the impacted column usually experienced four different phases, i.e., initial contact phase, loading phase having the approximate speeds of both the impactor and the RC column, unloading phase, and free vibration phase. The above procedure only determines the force-deformation relationship of the RC columns during loading phases, as shown in Fig. 8. If there is only the loading behavior of the RC column, the impact response from zero to the maximum displacement can be estimated. However, residual displacement cannot be calculated because of the absence of unloading behaviors. Actually, the previous experimental studies (Fan et al. 2019) demonstrated that the residual deformation is more critical for estimating the damage severity and the collapse risk of the impacted column. Accordingly, the unloading behavior of the RC columns was carefully examined herein. As illustrated in Fig. 8a, the complete force-deformation relationship (i.e., hysteretic model) needs to be defined to capture the entire interaction process between the impactor and the impacted column. In seismic analyses, many hysteretic models have been developed to capture the cyclic behavior of the RC beam and columns under earthquake excitations (Chopra 2007;Priestley et al. 1996). The characteristics of impact loads (e.g., unidirectional loading with high intensity) are significantly different from those of earthquake excitations (e.g., relatively long duration and cyclic loading) ( Fig. 7 Drop-hammer impact tests on axially-loaded circular RC columns: a test setup; b reinforcement layouts (Liu et al. 2017) Page 12 of 24 Liu et al. Int J Concr Struct Mater (2020) 14:32

Table 1 Column specimens and results obtained from experiments and simplified models
For the nomenclature of the test case, the letter E is used to represent the impact energy, and its subsequent numbers (i.e., "1" and "2") denote the low-energy impact (5198.5J), and the high-energy impact (13,364.9 J), respectively. The letter F means the application of the axial load ratio, and its subsequent numbers ("1", "2", and "3") denote the axial load ratios of 0%, 14%, and 28%, respectively. The letter L is employed to represent the longitudinal reinforcement, and its followed values of "12" and "06" denote the numbers of the longitudinal steel bars used in the specimens. The letter S denotes the transverse (stirrup) reinforcement, and its followed numbers ("100" and "055") means the pitch of spiral stirrups. The nomenclature is the same as the previous study (Liu et al. 2017 2014). Like blast-induced responses, the impacted column always reaches its maximum deformation at the first positive cycle, as shown in Fig. 8b. Hence, refer to the multi-linear curve model for blast analysis suggested by (Krauthammer et al. 1988), Fig. 8a presents the hysteretic model for RC columns under impact loading. In the hysteretic model, the initial loading curve ("o-r-a" curve) was obtained from the procedure mentioned in Sect. 2.2.3. The sequences of unloading and re-loading (i.e., "a-b", "c-d", "d-e", and "e-f ") were linearized for simplification.
To discuss the influence of the unloading stiffness ( k un ) on the impact responses of circular RC columns, four impact tests on axially-loaded circular RC columns were used herein, and the remaining four tests were adopted to validate the proposed approach in the following section. As a straightforward treatment, the unloading stiffness is usually assumed to equal the initial loading stiffness ( k el ) of RC members in the elastic range. For example, the treatment has been used in the impact analysis of the rectangular RC beam (Fujikake et al. 2006(Fujikake et al. , 2009Krauthammer et al. 1988). However, it was observed from Fig. 8c, d that the numericallyobtained residual deformations of the impacted columns were much larger than the experimental data when k un = k el . It means that the initial elastic stiffness is not appropriate for the circular RC column. This is due to the fact that there are differences between rectangular RC beams and circular RC columns in terms of the reinforcement layout, the boundary condition, and including axial loads or not. For example, the reinforcing steel bars in the rectangular RC beams yield together because of the identical distances to the neutral axis for these reinforcing steel bars. However, the reinforcing steel bars in the circular RC column yield gradually. From the point of view of energy, the unloading phase can be characterized by the release of elastic strain energy ( S el ) that is accumulated in the loading phase. As shown in Fig. 8a, k un can be readily estimated from the release of elastic strain energy and the column resistance P un at the loading to unloading point (a). For Fig. 8b, the point (a) means that the column displacement is maximum, and the velocity of the column and the impactor is approximately zero. For the impact tests in (Liu et al. 2017), the gravity of the hammer is less than 4% of the column resistance. Hence, the experimentally-measured impact force versus displacement data in the unloading phase ( Fig. 9) can be used to approximate the force-deformation relationship of the column during unloading and to estimate k un as follows where P un is the impact force corresponding to the critical point from loading to unloading (Fig. 9). Based on the experimental data and Eq. (43), the values of k un are presented in Table 2 for the tested columns. As shown in Fig. 8c, d, good agreements were achieved between (43) k un = 2S el P 2 un the experimental data and the numerical results when using these values of k un . This means that the approach using Eq. (43) to estimate k un is reasonable. However, the method requires the experimental data, which would block its applicability in practices because tests cannot always be performed due to high economic and time costs. Therefore, it is evident to develop a simple method independent of experimental data to reasonably estimate the unloading stiffness.
The critical implication from the above analysis is that the unloading stiffness is dependent upon the elastic strain energy. As mentioned above, it is simple for the rectangular beam to determine the status transition from elasticity to plasticity because almost all the reinforcing steel bars in the tension part began to yield at the same time. Distinct to the rectangular beam, although the steel bars at the bottom yield in the circular crosssection member, many steel bars are in the elastic range. It implies that the elastic strain energy is still accumulated after the steel bars at the bottom yield. It is one of the crucial reasons why the initial loading stiffness is not appropriate to approximate the unloading stiffness in the circular cross-section member. Similar to the rectangular beam, the nominal critical point can be defined when all the reinforcing bars below (or above) the neutral axis ultimately yield in the circular cross-section. Accordingly, the equivalent unloading stiffness (k un,eq ) can be defined as where P eq is the column resistance at the nominal critical point; u 2,eq is the column deformation at the crucial point. As shown in Fig. 10 and Table 2, the equivalent unloading stiffnesses using Eq. (44) are close to those estimated from Eq. (43) and the experimental data. The numerical results obtained using Eq. (44) also have a good agreement with the test results. In contrast to the above approach, the method using Eq. (44) is independent of the experimental data. Therefore, the method using Eq. (44) was recommended for circular columns to (44) k un,eq = P eq u 2,eq Impact force (

Results and Discussions
Figure 1l presents the impact force-time histories and the midspan displacement obtained from the simplified two-DOF models along with the experimental data. The primary peak responses obtained from the simplified models are shown in Table 1 for comparison with the experimentally-measured results. Generally, good agreements were observed between the results obtained from the two-DOF model and the experimental data for all impact cases, indicating the rationality of the developed two-DOF model. In addition to the results of the two-DOF model, the FE results obtained using the conventional modeling method are given in Fig. 11 for comparison. Apparently, compared with the FE results obtained from the detailed model, the simplified model provided much better accuracy in the prediction of the impact-induced responses, particularly for the durations of impact forces, the residual and maximum deformations, as shown in Fig. 11. Although the improved FE modeling method was provided in the previous study (Liu et al. 2017) to improve the accuracy in the prediction of the residual deformation, its coarse approximations and assumptions (e.g., specifying different RECOV parameters in the concrete material model to consider the influence of crack closure) due to the limitation of the concrete material model would place some difficulties in practical applications. Also, the computational efficiency was improved when using the simplified model instead of the detailed FE model. About 11 h per an FE run were reduced to a couple of hours using the simplified model with the force-deformation relationship derived from sectional analysis.

An Efficient Method to Obtain Structural Resistances
In the above analysis procedure, the sectional analysis was used to estimate structural resistance. Meanwhile, because the strain-rate effects were taken into account step-by-step, the iteration process (as shown in Fig. 6) was required to update the force-deformation relationship of the column. For this reason, although the computational time is reduced by the use of the simplified model, computational efficiency is expected to be further improved. On the other hand, coding to implement the section analysis would place some difficulties in the application of the simplified two-DOF model. In the context, fiber-based beam-column element models were Page 16 of 24 Liu et al. Int J Concr Struct Mater (2020) 14:32 developed to determine the force-deformation relationship of the column efficiently.

Fiber-Based Beam-Column Element Model
For the impact tests on circular RC columns in (Liu et al. 2017), Fig. 12 shows the fiber-based beam-column element model that was created by using OpenSees (Mazzoni et al. 2006). In the nonlinear beam-column model, the column was fixed at the left end, and the axial degree of the column at the right end was released. The preload in the axial direction was exerted at the right end of the column. The material model Concrete04 in OpenSees was used to simulate the concrete, and reinforcing steel bars were modeled by Steel02. Similar to the treatment in the section-based analysis, the Mander's model was employed to describe the behavior of the concrete core. Also, the material properties (e.g., concrete strength, the stress versus strain curve of steel bars) listed in (Liu et al. 2017) were used to define these material models.
The force-deformation relationship of the column was obtained by gradually increasing the lateral displacement at the impact point (i.e., displacement-based loading).

Approximation of Strain-Rate Effects
If the strain-rate effects were accurately considered, the user's development of the FE code and the iteration process still could not be avoided in spite of the use of the fiber-based beam-column element model. Figure 13 shows the changes in the velocities of m 1 and m 2 with time during the impact analyses. It was found that the speeds of the column at the impact point were usually approximate to those of the impactor except the initial contact phase, and their speeds varied on a relatively small range during the impact analyses. More importantly, due to the logarithmic relationship between the dynamic increase factors and the strain rates of concrete and reinforcing steel, the dynamic increase factors of concrete and steel bars are not sensitive to the change of velocity in low-velocity impact analyses. For convenience, the average impact velocity (i.e., V 0 2 ) was employed to estimate the influence of strain-rate effect approximately. Accordingly, the average curvature rate ( θ av ) for the test column can be determined by The average curvature rate ( θ av ) was substituted into Eqs. (30) to (32) to estimate the equivalent strain rates (45) θ av =2V 0 ll p  Liu et al. Int J Concr Struct Mater (2020) 14:32 of concrete and steel bars. Since the neutral axial depth changed slightly after the initial yielding of steel bars, the values of x y and x u in static analyses were adopted for simplification. Once the strain rates of concrete and steel bars were determined, the material properties of concrete and steel bars with strain-rate effects (Eqs. (6) to (8)) could be obtained for the material definitions in the fiber-based beam-column element models. It is worth mentioning that the approximation of strain-rate effects can also be used in the section-based resistance analysis to avoid the iteration and to improve the computational efficiency further. Figure 13 presents the force-deformation curves obtained from the fiber-based element models with the approximate strain-rate effects along with the above section-based analysis results with iteratively-updating strain-rate effects. Suitable matches were observed between the results obtained from the fiber-based element model and the section-based analysis results. It implies that the suggested approximation of strain-rate effects is reasonable and acceptable in the determination of the force-deformation relationship of the impacted column.

Structural Resistances Derived from Fiber-Based Element Models
In addition to the loading portion of the force-deformation relationship, the unloading behavior of the impacted column needs to be further investigated based on the fiberbased element model, as discussed in Sect. 3.2. According to the inspiration of Sect. 3.2, the influence of the maximum deformation in the loading phase is worth exploring in the fiber-based analyses. Figure 14 presents the unloading stiffness ( K un ) versus the maximum deformation curves. Since the nonlinear behavior (e.g., concrete cracking and yield of steel bar) occurred, the unloading stiffness of the column rapidly decreased with the increase of the maximum deformation in the initial stage. Subsequently, the unloading stiffness tended to a stable value after the maximum deformation was relatively large (e.g., higher than 0.02 m). Slightly increases in k un were observed for the column with axial loads after the minimum value of k un . It is because that the presence of an axial load plays a positive role in reverse unloading, and its lateral supporting role becomes more remarkable with the increase of the lateral displacement. For the impact cases reported in (Liu et al. 2017), the maximum displacements and the corresponding values of k un were marked out by the dashed line in Fig. 14. It was found that the values of k un were close to the results obtained from Eqs. (43) or (44). The change (gradual decrease) in k un at the initial stage was consistent with the observation from Fig. 10 when the deformation of the column is not beyond the nominal critical point of complete yield. As shown in Fig. 14, the explicit formulas were derived from fitting the numerical data to apply them in the analysis of the simplified two-DOF model.  Figure 11 presents the results of the simplified two-DOF models when using the force-deformation relationships derived from the fiber-based models instead of the section-based analysis. The simplified model using the fiber-based method is named as the proposed model in Fig. 11 to distinguish it from the simplified model with iterations. Similarly, the simplified two-DOF models were shown to have excellent accuracy in the prediction of the dynamic responses of the circular axially-loaded RC columns under impact loading. These results confirmed the applicability of the fiber-based element model to obtain the force-deformation relationship, including loading and unloading portions. Because of the approximation of strain-rate effects and the avoidance of iteration, the computational time per one run of the two-DOF model was reduced from 2 h to just ~ 30 s. Indeed, compared with the section-based analysis, an additional fiber-based beam-column element model needs to be Page 19 of 24 Liu et al. Int J Concr Struct Mater (2020) 14:32 developed. Extra energy and time are required for the development of the fiber-based beam-column element model. The proposed simplified method is expected to be used in the prediction of column responses under lateral impacts from falling rock, vehicle, and vessel. Bertrand et al. (2015) pointed out that the falling rock can be regarded as a rigid body like the drop-hammer impact tests. Hence, the proposed simplified method can be used directly in rockfall impacts except that the impact point is adjusted based on the actual impact. For vehicle and vessel impacts, the nonlinear behavior of the impactor cannot be omitted, and some particular SDOF models with the inelastic springs have been developed in current studies (e.g., Cao et al. 2019;Chen et al. 2016;Consolazio and Cowan2005;Fan et al. 2011) to describe the vehicle or vessel behavior during collisions. In this case, these particular SDOF models would be used to replace the elastic contact model of the drop hammer in this study, and the modeling method of the impacted column developed in this study is still suitable.

Impact Damage Evaluation of Axially-Loaded RC Columns
Because of the high efficiency of the simplified two-DOF model, the damage of RC columns can be quickly evaluated under different impact scenarios, particularly in the preliminary design phase. The iso-damage curves are presented in this section for the RC columns under various impacts.

Analysis Matrix
A parametric study was performed using the developed two-DOF models described previously to obtain the iso-damage curve of a column subjected to various impacts. For the two kinds of RC columns studied previously, parameters that were varied in this section included axial load ratio ( α = 7%, 14%, and 28% ), initial impact velocity ( V 0 = 2 to 18 m/s), and impact mass ( m 1 = 200to 1800kg) . A total of 6534 analysis cases were run to obtain the iso-damage curve of these two kinds of RC columns.

Damage Index
Bridge or building columns are usually designed to support substantial axial loads. For RC columns that have been subjected to lateral impact loads, it is essential to evaluate the residual axial load-carrying capacity ( P resi ) that represents the damage severity and the risk of collapse. For this reason, much research has chosen the damage index based residual strength in the damage evaluation of columns after impact or blast loading. Hence, the remaining strength of columns after impact was used to describe the impact-induced damage extent. Because of the presence of the axial preload, the damage index DI based on residual strength can be defined by where N m is the axial capacity of the undamaged column; N 0 is the preload in the axial direction. A formula was developed in the previous study based on a high number

Fig. 14 Unloading stiffness versus midspan displacement
Page 20 of 24 Liu et al. Int J Concr Struct Mater (2020) 14:32 of FE results and multivariable regression analysis to predict the residual strength of circular RC column with impact-induced flexural damage as follows (Fan et al. 2019): where u 2,resi is the residual deformation of the column after impact; ρ g is the longitudinal reinforcement ratio; ρ v is the transverse reinforcement ratio. Substituting Eqs. (47) to (49) (Fan et al. 2019) to predict the residual axial strengths of the impact-damaged RC columns. In the developed method, the residual deformation shape of the impact-damaged column was first approximated by the displacement-based push loading on the impact point before the simulation of the residual axial performance. This method, based on the approximation of the postimpact state (e.g., residual deformation), has been demonstrated by the compression after impact (CAI) tests on the impact-damaged columns (Fan et al. 2019). to predict the residual axial capacity well. Extensive parametric studies were performed based on the validated method. Based on the analysis results, Eqs. (46) to (49) were derived through multivariable regression analysis. Good agreements were achieved between the developed formula (i.e., Eqs. (47) to (49)) and the numerical and experimental results. In addition, it should be noted that the required residual deformations in Eq. (47) were estimated by the proposed 2-DOF models rather than the detailed FE models. Figure 15 plots the iso-damage curves of the axially-loaded circular RC columns subjected to various impact scenarios. For comparison, the iso-impact

Results and Discussions
energy curves are shown in Fig. 15 by the dotted lines. In Fig. 15, the color portions show different damage severity under various impact mass and velocity. The empty area in the lower-left means that the column has negligible damage, while the blank space in the upper right indicates that the column has wholly lost the axial-load carrying capacity.
Although the iso-damage curves and the iso-impact energy curves have a similar trend with impact mass and velocity, some differences between them were observed, as shown in Fig. 15. In other words, the impact-induced damages of the columns were different at the same impact energy with varying combinations of mass and velocity. It is attributed to the influence of the strain-rate effect of concrete and steel bars so that the impact resistance of the column exhibit some differences. By comparing with the results in Fig. 15a-c or d-f, the axial load ratio had an essential influence on the impact-resistant performance of the columns. Under the same energy impacts, the damage of the column decreased with the increase in axial load ratios (from 7 to 28%). Apparently, the empty area in the lower-left enlarged with increasing axial load ratios. As shown in Fig. 16, the presence of the axial load increased the bending resistance of the columns. Hence, when the impact energy was relatively low, and the impact-induced responses were also small, the axial load played a positive role in the impact-resistant performance of the column. However, the empty area in the upper right also increased with increasing axial loads. It implies that the energy dissipation capacity of the column under lateral impacts reduced with increasing axial loads. Distinct to the low-energy impacts, the presence of the axial load had a negative influence on the impact-resistant performance of circular RC columns. It is attributed to the fact that the P-Delta effect becomes significant with the increase of the impactinduced deformation. By comparing with the results in Fig. 15a, d with the same axial loads, the impactresistant performance of the RC column can be significantly improved by increasing reinforcement ratios because it improved the bending resistance of the RC column. Notably, the static axial load-carrying capacity of the column with the higher reinforcement ratio was 22% higher than that of the column with the relatively smaller reinforcement ratio, while the energy dissipation capacity of the column with higher reinforcement ratio was twice as much as that of the latter one.

Conclusions
For circular RC columns, a simplified modeling method was developed to efficiently predict the impact responses. Two different approaches were presented to obtain the Page 21 of 24 Liu et al. Int J Concr Struct Mater (2020)  force-deformation relationship of the impacted columns, and the influence of the unloading stiffness on impactinduced responses was discussed in detail. The proposed simplified method was validated by the impact tests on circular RC columns. The main conclusions can be drawn as follows: 1. The developed two-DOF model was demonstrated to be capable of predicting the responses of axiallyloaded circular RC columns under impact loading. Also, the computational efficiency using the two-DOF model would increase by three orders of magnitude in comparison with the typical detailed FE simulations. Therefore, the proposed method is especially suited for the preliminary design phase. 2. For the column with a circular cross-section, the formulas and procedure for the force-deformation relationship with strain-rate effects were presented according to the section-based analysis. Two simple methods to determine the unloading stiffness of circular columns were provided and validated for the complete definition of the hysteretic model in the impact analysis. Among them, the method based on Eq. (44) can be readily applied to impact studies independent of the experimental data. 3. The fiber-based beam-column model was developed to determine the force-deformation relationship of columns considering strain-rate effects. The proposed approximation of strain-rate effects was dem-onstrated to be reasonable and acceptable. In addition to the loading portion, the developed fiber-based beam-column model is able to reasonably capture the change of the unloading stiffness with the increase of the impact-induced maximum deformation. 4. The damage characteristics of axially-loaded circular RC columns were widely investigated by using the proposed modeling method under various impact scenarios, and the corresponding iso-damage curves were plotted for two kinds of RC columns. It was found that the impact-induced damages of the columns were different at the same impact energy with varying combinations of mass and velocities due to strain-rate effects. 5. The presence of axial loads had a significant influence on the impact-induced damage of axially-loaded circular RC columns. Under the low-energy impact, the axial load usually plays a positive role in the impact resistance of the RC column. On the contrary, increasing axial loads would reduce the energy dissipation capacity of the column because the P-delta effect becomes significant with the increase of the column's deformation.
Similar to most of the simplified models, the developed method in this paper is only applicable to evaluating the flexural responses of circular RC columns under impact loading. The shear damage (failure) should be avoided by increasing the transverse reinforcement  Liu et al. Int J Concr Struct Mater (2020) 14:32 ratio of the column. Indeed, the advanced analytical model is worth developing in the future to capture both the flexural and shear-dominated damage of columns under impact loading. In addition, the rigid impactor was assumed to strike the center zone of the specimen in the impact test used for validations in this study. However, impact positions are often uncertain in real collision events. Because of the differences between the impact test and the real collision event, further studies should be performed to experimentally investigate the influence of impact position and column aspect ratio and improve the developed simplified method for applicability in the future.

Abbreviation
Nomenclature m 1 : Mass of impactor; m 2 : Equivalent mass of RC column; u 1 ,u 1 ,ü 1 : Displacement, velocity, and acceleration of t impactor, respectively; u 2 ,u 2 ,ü 2 : Displacement, velocity, and acceleration of RC column, respectively; c 1 , k 1 : Contact damping and stiffness between impactor and RC column, respectively; c 2 , k 2 : Equivalent damping and stiffness of RC column, respectively; R s : Equivalent resistance of RC column; σ c , ε c : Compressive stress and strain of concrete, respectively; f co , ε co : Compressive strength of concrete and the corresponding strain, respectively; f ′ co , f ′ cc : Compressive strength of the unconfined and confined concrete, respectively; r, E sec , ε cc : Parameters of concrete constitutive model (see Eq. (4)); E c , E s : Elastic modulus of concrete and rebar, respectively; σ s , ε s : Stress and strain of reinforcing steel, respectively; f y : Yield strength of reinforcing steel; ε: Strain rate; E b : Hardening modulus of reinforcing steel; N, M : Axial force and bending moment of the cross-section; x y , x u : Location of the neutral axis at yield and ultimate points (see Fig. 4), respectively; D: Diameter of RC column; d e : Effective height of the cross-section; d s : Distance from the rebars at the bottom to the edge of the column (see Fig. 4); d c : Thickness of concrete cover; r s : Radius of steel rebar (see Fig. 4); s j : Distance from the rebars at the jth layer to the rebars at the bottom (see Fig. 4); y: Coordinate value along the circular cross-section (Fig. 4); A s : Area of steel rebar; A: Areas of RC column; θ, ϕ: Curvature and rotation of RC column, respectively; EI: Flexural rigidity; k: Stiffness parameter defined in Eq. (27); l : Length of RC column; l p : Length of plastic hinge; P: Lateral load of RC column; d di : Distance from the i-th layer bars to the outermost layer bars in tension; d upi : Distance from the i-th layer bars to the outermost layer bars in compression; m c (z): Distributed mass along the column length; φ(z): Shape function of RC column; S el : Elastic strain energy accumulated in the loading phase; P un : Impact force corresponding to the critical point from loading to unloading; k un : Unloading stiffness of RC column; k un,eq : Equivalent unloading stiffness of RC column; α: Axial load ratio; DI: Damage index; ρ g , ρ v : Longitudinal and transverse reinforcement ratios, respectively; β: N resi N m . Residual.
Authors' contributions BL developed the simplified modeling method of this research and concluded the proposed method in this research; WF proposed the main idea of this research, conducted the physical experiments for validations and wrote the manuscript. XH participated in the development of the simplified modeling and the manuscript preparation; SD reviewed the entire manuscript critically and participated in the discussion of numerical and experimental results; LK participated in the validation of the proposed simplified modeling method and reviewed the entire manuscript critically. All authors read and approve the final manuscript.