Earthquake Response Analysis of Tall Reinforced Concrete Chimneys considering Eccentricity

A numerical algorithm is presented to analyze earthquake response of tall reinforced concrete (RC) chimneys based on stick multidegree-of-freedom models. The algorithm considers the eccentricity phenomena between spatial discrete nodes and corresponding centroids of investigated lumps. The spatial discrete segments of the chimney are used to construct the investigated lumps. The equations of dynamic equilibrium of the investigated lumps are derived, and the numerical calculation procedure is implemented. Phenomena of eccentricity are studied for 150 m and 210 m RC chimneys. Seismic stresses and effects of vertical ground motion for the two chimneys are also studied. Numerical results show that the tensile and compressive stresses of the seismic control cross sections of the chimneys may increase under the actions of several specific earthquake waves by considering existing eccentricities. The effect of eccentricity on the earthquake responses of tall RC chimney should be considered, and stresses caused by vertical ground motion should not be neglected to obtain accurate earthquake response of chimneys.


Introduction
1.1. Literature Review. Many tall reinforced concrete (RC) chimneys have been built in places such as coal-fired power stations and refinery because of the requirements of air pollution control standards. Slender tall chimneys are sensitive to earthquake. Strong earthquake threatens the safety of the chimney and other buildings around it because of the nonuniformity of the diameter along the height of the chimney and the variation in dead weights at different heights. Research and calculations of earthquake response are necessary for the designs of tall RC chimneys under the action of earthquakes.
Destruction and collapse of tall chimneys bring huge economic losses. Research on seismic damage of chimneys indicates that the destructive effect of vertical earthquake on chimneys should be considered [1]. Liu and He [2] analyzed the seismic response of a chimney under the action of combined horizontal and vertical ground motions. ey concluded that vertical stress cannot be neglected completely although horizontal seismic stress is the controlling part. Wang and Liao [3] developed a calculation technique for analyzing the vertical earthquake response of a chimney based on multiple reflections of traveling waves in the chimney. Liu et al. [4] conducted earthquake transient response of multidegree-of-freedom (MDOF) systems by the model superposition method. Ding et al. [5] conducted a vertical impact to the bottoms of chimney models. e test results showed that only one circle crack appeared in the chimney models, and the crack had relations with the ratio of input wave frequency to the natural frequency of the chimney model. Chen et al. [6] used a shaking table to study the fracture phenomenon of a chimney under the action of vertical earthquake. ey found that vertical seismic action should be considered in the structural design of chimneys. A 115 m high RC chimney collapsed on August 17, 1999, during the Izmit earthquake in Turkey [7][8][9][10]. Huang et al. [7,8] investigated a collapsed chimney with methods of linear dynamic response spectrum and nonlinear static capacity spectrum. Akinci [9] used shell elements to analyze a collapsed chimney. e analysis results indicated that low ductility and inadequate capacity due to a large duct were the two significant factors for the collapsed chimney. Huang and Gould [10] further proposed a new three-dimensional (3-D) pushover analysis method to study the dynamic response of the stack due to an earthquake motion recorded at a nearby site. Lopes et al. [11] constructed a realistic numerical model to evaluate the seismic vulnerability of a chimney and simulated its damaged state. Minghini et al. [12] analyzed a brickwork chimney and a shear failure mechanism of the upper part of the chimney damaged by the 2012 Emilia earthquake. Zhou et al. selected an existing 240 m tall RC chimney as the research object and analyzed the seismic vulnerability curve of the chimney by considering structural damage [13], structural nonlinearity [14], and multidimensional earthquake [15] and evaluated and analyzed the vulnerability of the chimney structure. e modal frequency and mode shape are the main dynamic characteristics of chimneys. Chen et al. [16] used the Adomian decomposition method to decompose dimensionless modal shape and calculated the dynamic characteristics of a high-rise chimney. Sancibrian et al. [17] used the continuum structure system to analyze the dynamic characteristics of a chimney, including natural frequencies and mode shapes. Tatara and Ratajewicz [18] evaluated the dynamic characteristics of a typical RC industrial chimney that belongs to a heating plant and determined the lowest vibration frequency of the structure.
Soil-structure interaction (SSI) has a great influence on the seismic performance of tall chimney structures [19,20]. Halabian and Kabiri [21] studied the effect of foundation stiffness of different types of RC chimneys on the ductility reduction coefficient. A series of inelastic and elastic time history analyses up to the collapse status for a large variety of RC chimneys subjected to five scaled artificial records were conducted. Jisha et al. [22] analyzed 3-D SSI of tall RC chimneys with annular raft foundation and evaluated the effect of SSI by changing the stiffness of supporting soil and foundation [23]. Zhou et al. [24,25] studied the effect of SSI on the seismic collapse resistance of chimney structure for a 240 m high RC chimney.
Time period is an important approach of evaluating the structural response of a structure subjected to large-scale earthquake vibrations. e seismic response of chimney structures uses time history analysis considering earthquake, which is a long-duration earthquake impulse [19]. Sarkar et al. [26] studied the effects of different base conditions on the design value of shear force and bending moment of different sections of a 275 m tall RCC multiflue chimney by using a response spectrum method. e spectrum compatibility time history of specific engineering site was considered, and the time history dynamic analysis of the chimney structure was conducted in detail. Wilson [27] studied the inelastic response of 10 RC chimneys ranging from 115 m to 301 m subject to earthquake excitation. Liu et al. [28] studied the structural dynamic performance and elastic-plastic seismic response of a 240 m high special-shaped RC chimney. e weak position of the chimney was discussed by analyzing its interlayer displacement angle.

Motivation.
When the stick MDOF structural model is used to study the dynamic response of tall chimney structures, the chimney structure is usually discretized into beam elements of a certain length. in beam elements are good, but they are not sufficiently thin in practical analysis. e cross section area of a chimney structure is variable at different heights. erefore, even if the discrete length of space is equal, the vertical eccentricity of mass will occur because the mass of each element will be uneven. e vertical eccentricity of mass must be considered to improve the accuracy and validity of dynamic response analysis.
A numerical algorithm is presented in this study to analyze the earthquake response of tall RC chimneys. e chimney is spatially discretized into several segments of finite length, and the investigated lumps are then constructed by using these segments. e calculating formulae of axial force, shear force, and bending moment acting on the median cross section of each segment are derived. ese interior forces are associated with end displacement and angle. Not only the inertial forces of the investigated lumps are considered but also the rotary inertias of investigated lumps. Existing eccentricity phenomena between spatial discrete nodes and corresponding centroids of investigated lumps are considered. Dynamic equations and relations are applied alternately in the time domain to calculate the seismic response of the chimney. Two RC chimneys of 150 m and 210 m height are studied. e numerical results show that the effects of such eccentricities on earthquake response should not be neglected for tall chimneys, as well as the vertical ground motion.

Governing Equations of Investigated Lumps.
A chimney with variable cross section can be discretized into many segments of different cross section areas along the vertical direction in practical calculation based on the stick MDOF model. ese segments are used to construct the investigated lumps. Figure 1 shows a typical investigated lump i (the nodal numbering i is also used to describe the investigated lump) consisting of half the upper segment and half the lower segment around the discrete node i. Sign C is the centroid of the investigated lump. e centroid of the investigated lump is commonly not at the same point of the corresponding discrete node i for the chimney with variable cross-sectional area. Here, the distance between the discrete node and the centroid of the investigated lump within the same investigated lump is called eccentricity. e vertical equilibrium equation for the centroid of investigated lump i is where m i is the mass of investigated lump i, € w ci is the vertical acceleration of centroid C, and N l− 1 and N l are the axial forces acting on investigated lump i.

2
Shock and Vibration e transverse equilibrium equation for the centroid of investigated lump i is where € u ci is the horizontal acceleration of the centroid of investigated lump i and V l and V l− 1 are the shear forces acting on investigated lump i.
On the basis of mechanics theory, the rotational equilibrium equation should be established about the horizontal axis (y-axis) through centroid C of investigated lump i, that is, where J ci is the moment of inertia of investigated lump i about the y-axis through its centroid; € θ i and z ci are the angular acceleration and the eccentricity of investigated lump i, respectively; M l and h l are the bending moments and the length of segment l, respectively; and u i is the transverse displacement of node i. Figure 2 shows that the following relation holds: where € u i is the horizontal acceleration of node i. In Figure 2, € w i is the vertical acceleration of node i. Figure 3 shows the free-body diagrams of the upper and lower halves of segment l by an artificial cut at the central section of segment l. In Figure 3, N i ′ and N i+1 ′ are the axial compressive forces at the bottom and top ends of segment l, respectively, because of its vertical deformation; V i ′ and V i+1 ′ are the shear forces; M i ′ and M i+1 ′ are the bending moments; N ig is the axial force coming from the dead load of the structure above node i; and N i+1g is the dead load above node i + 1.

Forces and Bending Moment on the Central Cross Section of the Segment.
Following the beam theory considering shear deformation and axial force, we have is the coefficient of shear deformation in the transverse direction. (EI) l , (GA) l , and (EA) l are the bending stiffness, shear stiffness, and tensile/compressive stiffness of segment l, respectively. N l � N ig + N i ′ , where N ig is the axial force coming from the dead weight of the chimney above node i. A l is the cross-sectional area of segment l. w i and θ i are the vertical displacement and rotational angle of node i, respectively. w i+1 and θ i+1 are the vertical displacement and rotational angle of node i + 1, respectively.

Shock and Vibration 3
Given the equilibrium of the upper half of segment l, we have Given the equilibrium of the lower half of segment l, we have where N l , V l , and M l are the central axial force, central shear force, and central bending moment acting on the central section of segment l, respectively. We substitute equations (10) and (12) into equation (9), and we obtain the central bending moment of segment l, shown as follows: where K 5l � ((EI) l /h l ) − (N l h l /12).

Implementation of Algorithm.
Horizontal and vertical earthquake waves are inputted from the foundation of the tall RC chimney. e algorithm description diagram is shown in Figure 4, and the calculating procedure is as follows.
Step 1. Calculate the angular acceleration € θ i of investigated lump i at time t using equation (3), the vertical acceleration € w ci at time t using equation (1), € w i at time t by considering € w i being nearly equal to € w ci , the horizontal acceleration € u ci at time t using equation (2), and € u i using equation (4).
Step 2. Calculate angle θ i , vertical displacement w i , and horizontal displacement u i for investigated lump i at time t + Δt by integration in the time domain.
Step 3. Calculate axial force N i ′ at time t + Δt using equation (5) and then obtain central axial force N l using the formula N l � N ig + N i ′ . Central shear force V l and central bending moment M l are calculated using equations (11) and (13), respectively, at time t + Δt.
Step 4. Calculate € θ i , € w ci , and € u ci at time t + Δt as performed in Step 1. e circulation from Step 1 to Step 4 is implemented in the time domain by a computer program. e numerical algorithm is provided to calculate the earthquake response of the tall RC chimney, considering eccentricities z ci (i � 1, 2, . . . , n).
Wilson θ method is used to solve the conventional dynamic equations of equilibrium in which eccentricities are not considered. eir results are called the conventional results in this study.
Here, the rotation of the foundation is not considered.

Corresponding Relations.
For simplicity in deriving the stability condition of the algorithm, we assume that (EI) l � (EI), h l � h, and α l � α for all segments. Axial force N l is omitted in K jl (j � 1, 2, . . . , 5). Eccentric phenomena are neglected in equation (3); and m i � m and J ci � J i � J for all investigated lumps. On the basis of these assumptions, vertical governing equations are decoupled from the governing equations of transverse and rotational deformations. us, we only study the coupling equations of transverse and rotational deformations to provide the stability condition for the algorithm presented in this study.
On the basis of the above assumption, the corresponding relations between the coefficients used in Timoshenko beam theory and those used in the proposed algorithm are as follows (see Appendix A):

Stability Condition Based on Timoshenko Beam
Model. e flexural wave equation of Timoshenko beam is where . Substituting a flexural wave u � u 0 e i(ωt− kz) into the finite difference formula of Equation (15), we obtain the following equation (see equation (A.12) in Appendix A): Solving equation (16) with sin 2 (ωΔt/2) as unknown yields the following: 4 Shock and Vibration Given sin 2 (ωΔt/2) ≤ 1, we approximately have the following equation (see equation (A.16) in Appendix A):

Stability Condition for the Proposed
Algorithm. e stability condition used in the algorithm presented in this paper can be provided by using equation (18) and the corresponding relations (14) as follows:

Method Validation
We study bending wave propagation of a beam to verify the effectiveness of the new algorithm. e beam with equal cross-sectional area is designed to avoid eccentricity. e result given by the finite difference method [29] is used to test the proposed numerical algorithm. e elastic modulus of the beam is E � 209 GPa, density is ρ � 8000 kg/m 3 , Poisson's ratio is v � 0.3, and shear crosssectional correction factor is k′ � 0.886. e length of the beam is X � 6.35 m, and the radius is 0.127 m. e number of discrete segments is 500 in the length direction for the  Shock and Vibration 5 uniform beam. On the basis of the calculation parameters of the beam, the time step of the beam satisfying the stability condition of the method is Δt ≤ 1.51 μs obtained by equation (19), and the calculation time step is Δt � 1.25 μs. e length from the left end of the beam is x, and the bending moment at the section is M. After normalization, the time t is (1/X) ����� (E/ρ)t, the length is x � (x/X), and the bending moment is M � (MX/EI). e inclination moment of the left end of the beam is Figure 5 shows the comparison curves of the bending moments M at different cross sections of the beam at t � 7.0554. e calculation results of the new numerical method are nearly the same as those of the finite difference method. e comparison results show that the presented algorithm is valid for calculating the earthquake responses of tall RC chimney.

Earthquake Records Used.
Eleven earthquake acceleration records on four types of soil are used as earthquake load to calculate the structural earthquake response. Four types of soil are denoted by rock, soft rock, stiff soil, and soft soil. Soil type is characterized by shear wave velocity. e corresponding earthquake records are selected from the Pacific Earthquake Engineering Research Center website [30]. Table 1 shows the 11 earthquake records.
According to the peak value of design earthquake acceleration of 8 degrees seismic intensity in the literature [31], the peaks of the horizontal earthquake accelerations are adjusted to 0.2 g in practical calculations. e peaks of the vertical earthquake accelerations are adjusted to 0.13 g. e damping is omitted. [32] of 150 m and 210 m height are studied. Figure 6 shows the two models of RC chimneys.  Table 3 shows the parameters of the 210 m high RC chimney. According to the literature [33], the standard value of the compressive strength of C25 concrete is 16.7 N/mm 2 , the design value is 11.9 N/mm 2 , and the elastic modulus is E C � 2.8 × 10 4 N/mm 2 . e standard value of the compressive strength of C35 concrete is 23.4 N/ mm 2 , the design value is 16.7 N/mm 2 , and the elastic modulus is E C � 3.15 × 10 4 N/mm 2 . e elastic modulus of HRB335 reinforcement is E S � 2.0 × 10 5 N/mm 2 , the standard value of strength (yield strength) is 335 N/mm 2 , and the design value of strength is 300 N/mm 2 .

Two Models of RC Chimneys. Two RC chimneys
e discrete segments are all 10 m long in the vertical direction for each of the two chimneys. Each chimney will become part of such calculated model consisting of several cylinders of different cross-sectional areas because an actual chimney structure has variable cross-sectional areas. us, the phenomena of eccentricities between the nodes and the corresponding centroids of the investigated lumps appear in the calculation of the earthquake response of the actual chimney. e mass of each investigated lump consists of half the upper segment and half the lower segment around the discrete node. e moment of inertia of each investigated lump about the horizontal axis through its centroid is calculated. On the basis of the calculation parameters of the 150 m and 210 m high chimneys, the time step satisfying the stability condition of the method is obtained by equation (19), Δt ≤ 3.4 ms and Δt ≤ 3.1 ms in the 150 m and 210 m high chimneys, respectively. e time step is 1.0 ms.

Effect of Eccentricity.
Eleven earthquake waves are used. e effect of eccentricity on normal earthquake stress of the chimney are studied when horizontal and vertical components are inputted simultaneously. Figure 7 shows the spatial distribution curves of the maximum earthquake normal stresses of different cross sections for the 150 m high chimney under the actions of Landers (No. 1 in Table 1 denotes that the phenomena of eccentricities are considered. "Conventional" denotes the conventional results, in which the phenomena of eccentricities are neglected, for the earthquake responses of the real RC chimneys. Dashed lines ("tensile") denote tensile stresses, whereas solid lines ("press") denote compressive stresses. A comparison of the red solid lines and the blue solid lines in Figures 7 and 8 shows that the increased phenomena of peak values of the normal compressive stresses appear when eccentricities are considered.
A comparison of the red dashed lines and the blue dashed lines in Figures 7 and 8 shows that the increased phenomena of peak values of the normal tensile stresses appear when eccentricities are considered. Table 4 shows the details of the increase in percentage of maximum earthquake normal tensile/compressive stresses when phenomena of eccentricities are considered. Table 5 shows maximum earthquake compressive stress, mean value, and standard deviation of different sections of the 150 m high chimney under 11 earthquake      Shock and Vibration 9

Shock and Vibration
wave actions when vertical eccentricities are considered. Table 6 shows maximum earthquake compressive stress, mean value, and standard deviation when vertical eccentricities are ignored.
In Table 5, the comparison of the data shows that the maximum earthquake compressive stress obtained under different earthquake waves is significantly different, and the degree of the difference can be expressed by the standard deviation. To quantify the degree of difference, its standard deviation S is where n is the number of samples, that is, the number of selected earthquake waves; σ l is normal stress on the central section of segment l, that is, the sum of bending stress and axial stress, which is σ l � (M l /W l )+ (N l /A l ); W l is the bending-resistance modulus of the segment l; and σ is the mean normal stress of the central section.       Figure 9(a) is the comparison of the mean values of earthquake compressive stresses between Tables 5 and 6 and the corresponding standard deviation. e left side of Figure 9(a) is the comparisons of the mean values of earthquake tensile stresses and the corresponding standard deviation. e data in Tables 7 and 8 (see Appendix  B) are provided regarding maximum earthquake tensile stress, mean value, and standard deviation of the 150 m high chimney when considering vertical eccentricities and neglecting vertical eccentricities. In Figure 9(a), the red and pink columns are, respectively, the mean values of earthquake compressive and tensile stresses considering vertical eccentricities. e blue and green columns are, respectively, the mean values of earthquake compressive and tensile stress neglecting vertical eccentricities. e red and pink lines are, respectively, the standard deviations of different earthquake compressive and tensile stresses considering vertical eccentricities. e blue and green lines are the standard  deviations of different earthquake compressive and tensile stresses neglecting vertical eccentricities. Figure 9(b) shows the case of the 210 m high chimney (data in Tables 9-12). Figure 9 shows that the earthquake normal stresses increase evidently when considering vertical eccentricities. Significant differences of the maximum earthquake normal stress are observed under different seismic waves actions. Figure 10 is the boxplot of the maximum earthquake compressive stress distribution at different cross sections of 150 m and 210 m high chimneys under different seismic wave actions when considering eccentricities. ey show the mean value, minimum, maximum, and 25%-75% of earthquake compressive stresses. e distribution characteristics of the earthquake compressive stress data are also shown in the figures. e stress distributions near the bottom of the chimney are relatively intensive, whereas those at the middle and upper parts of the chimney are relatively scattered.    e results studied above show that significant differences of maximum earthquake normal stress under different earthquake wave actions. e comparisons demonstrate that the conventional results (eccentricities neglected) might be less than their actual tensile/compressive stresses in the seismic control cross section of the chimneys under the actions of several specific earthquake waves.
e calculation results show that eccentricities between the spatial discrete nodes and the corresponding centroids of investigated lumps should be considered in earthquake response of tall RC chimneys to obtain accurate calculation results.

Effects of Vertical Earthquake.
e phenomena of eccentricities are considered to study the effect of vertical earthquake on earthquake normal stresses. Two types of seismic load are applied for each of the four earthquake waves. In the first case, the combined horizontal and vertical earthquake waves (see red lines in Figures 11 and 12) are inputted. In the second case, only the horizontal earthquake wave (see blue lines in Figures 11 and 12) is inputted.
For the 150 m high chimney, Figure 11 shows that the peak normal tensile stress increases by 4% for the Landers wave (No. 1) and 10% for the Loma Prieta wave (No. 2) if vertical earthquake waves are considered. e peak normal compressive stress increases by 7% for the Loma Prieta wave (No. 2) and 10% for the Chichi wave (No. 4) if vertical earthquake waves are considered.
For the 210 m high chimney, Figure 12 shows that the peak normal tensile stress increases by 12% for the Landers wave (No. 1) and 9% for the Loma Prieta wave (No. 2) if vertical earthquake waves are considered. e peak normal compressive stress increases by 9% for Landers wave (No. 1) and 13% for Loma Prieta wave (No. 2) if vertical earthquake waves are considered.     Figure 13 shows the spatial comparison curves of the mean values of the maximum normal earthquake stresses for the 150 m and 210 m high chimneys under two types of seismic load. Figure 13(a) is obtained from the average values of Tables 5, 7, 13, and 14. Figure 13(b) is drawn from the mean values of Tables 9, 11, 15, and 16. Figure 13 shows that the normal stress of the chimney can be increased considering the vertical seismic acceleration. e results studied above show that the vertical ground motion should not be neglected, especially for tall RC chimneys built on rock sites or soft rock sites.

Conclusions
A numerical algorithm that can consider the existing phenomena of eccentricities is presented to calculate the earthquake responses of tall RC chimneys. e numerical results demonstrate that the tensile and compressive stresses may increase in the seismic control cross sections of the chimneys under the actions of several specific earthquake waves when existing eccentricities are considered. e phenomena of eccentricities should be considered to obtain accurate calculation results of earthquake responses for the seismic designs of tall RC chimneys. e maximum normal stresses of earthquake response appear within the upper half of the chimneys, and their values are considerably larger than those within the lower half of the chimneys. Vertical ground motion should not be neglected, especially for chimneys on rock sites and soft rock sites. and G are the elastic and shear moduli of the beam, respectively; ρ 0 is the mass density; A 0 is the cross-sectional area; and A s is the effective cross-sectional area for the shear deformation of the beam. e finite difference formula of equation (A.9) is (A.10)

B. Partial Calculation Results
In view vertical eccentricities, Table 13 shows the maximum earthquake compressive stresses, mean values, and standard deviations of different sections of the 150 m high chimney under 11 seismic waves when only horizontal seismic acceleration is inputted. Table 14 shows the results of maximum earthquake tensile stresses of the 150 m high chimney.
In view of vertical eccentricities, Table 15 shows the maximum earthquake compressive stresses, mean values, and standard deviations of different sections of 210 m chimney under 11 seismic waves when only horizontal seismic acceleration is inputted. Table 16 shows the results of the maximum earthquake tensile stresses of the 210 m high chimney.
Data Availability e analysis result data used to support the findings of this study are included within the article. e calculation data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest.