Open Journal of Civil Engineering
Vol. 3  No. 3 (2013) , Article ID: 36198 , 7 pages DOI:10.4236/ojce.2013.33016

Seismic Damage Estimation of an Actual Reinforced Concrete Structure Using Subset MCMC

Shigeru Kushiyama

Department of Architecture, Hokkai Gakuen University, Sapporo, Japan


Copyright © 2013 Shigeru Kushiyama. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received April 27, 2013; revised May 27, 2013; accepted June 5, 2013

Keywords: Failure Probability; Nonlinear Dynamic Response Analysis; Subset Simulation; Markov Chain Monte Carlo


To estimate seismic damage of structures under strong motions is very important to know true safety of structures. However, we have to deal with very small failure probability issue to investigate quantitatively. If we use standard MCS (Monte Calro Simulation) to discuss failure probability, e.g., 1 × 10−6 order, since we have to execute nonlinear dynamic response analyses of approximate 107 times, it is not realistic. Recently, a subset simulation, which reduces the computation time by replacing small failure probability into the product of conditional failure probabilities, was proposed. In this study, the subset simulation is applied to estimate failure probability of an actual reinforced concrete building with 11 stories, and discuss the safety of the structure by checking with design criteria.

1. Introduction

The final goal of structural design is quantitative estimation of structural safety. However, quantitative damage estimation, i.e., failure probability, subject to an individual strong motion is not obtained from the procedure based on a structural design code. Because a specific ground motion is ordinary derived from a acceleration spectrum, which reflects properties of many ground motions predicted at each site, given in the design code. Therefore, we know only whether a structure is safe for the above specific ground motion.

Recently, the subset simulation for reliability estimation of structures subject to strong motions was proposed [1]. A basic idea of the subset simulation is to reduce the computation time by replacing small failure probability into the product of conditional failure probabilities.

In this paper, the above simulation is adopted for reliability estimation of an actual reinforced building with 11 stories under the conditions of uncertainty on material properties such as story stiffness and story yield strength.

2. Failure Probability Using Subset MCMC

The concept of subset simulation proposed by Siu-Kui Au are as follows [2]. Given a failure event, let be a decreasing sequence of failure events so that,.

By the definition of conditional probability, we have

. (1)

The above equation means that a failure probability is replaced by the production of sequence conditional probabilities. There, P(F1) is obtained from standard MCS (Monte Calro Simulation), and P(F2) ~ P(Fm) are obtained from subset MCMC (Markov Chain Monte Calro). Under the condition that the limit state function Z is negative, failure probability using subset methods is expressed by the following equation.

. (2)

where: number of samples generated in subset : number of samples defined in subset

, for example.

: number of samples satisfing in the last chain level,.

: level number of Markov chain.

On the other hand, a conditional probability density function is expressed by the following equation.

. (3)

where: index function, if is included in, others.

3. Ground Motions and Analytical Model

Table 1 shows input ground motions. Ground motions of No.1 - No.7 were observed at HKD180 point of K-net. No.8 - No.11 are design ground motions which are often used for dynamic analysis in Japan. In this study, these ground motions are normalized so that PGV (Peak Ground Velocity) leads to 50 kine, then we call them the level 2 ground motions in Japan. Figure 1 shows an example of normalized ground motions, and Figure 2 shows the response spectra of No.8 - No.11 ground motions.

Analytical model is an actual reinforced concrete residential building with 11 stories and two spans in X and Y direction. Table 2 shows the structural heights and the story weights of this model.

In nonlinear dynamic response analysis, it is assumed that a skeleton curve is tri-linear, and a hysteresis rule is masing type. In Rayleigh damping ratio, it is assumed that first and second damping of natural vibration mode is 5%. Time-stepping solver for the equation of motion adopts Wilson’s θ method (θ = 1.42).

As structural failure index, two indexes of MCDR (Maximum Column Drift Ratio) and MDD (Maximum Ductility Demand) are considered. In addition, it is assumed that strong sway observed on tall buildings after ground motion shake down does not occur about 11 stories building, and the following cutoff time is used.

. (4)

Figure 1. JMA Kobe 1995 NS normalized earthquake.

Figure 2. Response spectrum of No.8 - No.11. Damping ratio: h = 5%.

Table 1. Input ground motions.

where peak-time is the latest time among occurrence times of maximum absolute amplitude on PGA (Peak Ground Acceleration), PGV (Peak Ground Velocity) and PGD (Peak Ground Displacement).

In the subset simulation, it is assumed that material properties are uncertain and input ground motions are deterministic. Concretely, first and second story stiffness, story yield strength, i.e., SK1, SK2, SSy, at each floor are treated as uncertain variables. Furthermore, it is assumed that only two variables in (SK1, SK2) or (SK2, SSy) vary at a simulation in order to simplify. Table 3 shows the mean values of these three variables specified by transforming the relationship between interstory drift and story shear force, which is derived from pushover analy-

Table 2. Structural height and story weight of analytical model.

Table 3. Story yield strength and 1st and 2nd story stiffness (mean value).

sis, into a tri-linear skeleton curve. Each parameter is the set of i.i.d. (independent identically distributed) according with the following Gaussian distribution:

. (5)

where: coefficient of variation, it is assumed that is equal to 5% for the combination (SK1, SK2) and 2% for the combination (SK2, SSy).

Calculation flows are as follows. Firstly, standard MCS is executed. Response analysis of times per a simulation is repeated in MCS, and then a weakest floor is specified. In the weakest floor, frequencies that interstory drift or ductility demand shows maximum values is the biggest among all floors.

Secondly, the subset simulation is executed. Samples of the weakest floor are generated using MCMC.: the number of samples generated about each subset level, is set 600 based on convergence diagnostics by Raftery and Lewis [3]. And to the samples of other floors, the values generated in MCS are applied in order to reduce computation time for the parallel simulation mentioned later.

4. Basic Response Properties

Natural elastic period of analytical model is 0.764 (sec) in X direction, and 0.8204 (sec) in Y direction. Firstly, we describe response results subject to the 50 kine normalized 1995 Kobe NS shown in Figure 1.

Figures 3(a) and (b) show the restoring force and the ductility demand in X and Y direction frame. According to the restoring force, it is clarified that drift occurs and its degree in X direction frame is slightly larger than Y


Figure 3. Restoring force and ductility demand. (a) X direction frame; (b) Y direction frame.

direction. Under the comparison of ductility demand in X and Y direction, the damage of X direction frame is obviously larger than Y direction frame. The tendency was same about results subject to other earthquake motions. Therefore, in estimation of failure probability mentioned later, only X direction frame will be described. Incidentally, the ductility demand in X direction frame is over 1 at other floors except the top floor.

Figure 4 shows time variation of energy dissipated in X direction frame. According to this figure, yielding does not yet occur although cracking occurs in the top floor, however yielding occurs in first floor. Figure 5 shows the ductility demands of No.1 - No.11 input motions. According to this figure, the maximum values of, i.e., are over 1 for other ground motions except No.1, 4, 5, and appear not only at lower story but also at upper story or middle story.

5. Sensitivity of Parameter and Quantitative Estimation of Failure Probability

5.1. Sensitivity of Parameter

The subset simulation can automatically find out a parameter to reduce the value of limit state function, in other words to make larger failure index, without a special knowledge. However, we know only the relative intensity on sensitivity of two parameters. In sensitivity analysis, parallel simulation is executed one time since the effect by initial value’s discrepancy on parameter is


Figure 4. Time variation of energy dissipated in X direction frame. (note) Ed: damping energy, Ec: cracking energy, Ey: yielding energy, Ek: kinetic energy, Es: strain energy. (a) 11F; (b) 1F.

Figure 5. Ductility demand (X direction frame).

small. Then, response analyses are executed 3000 times (=600 in MCS + 600 × 4 in four levels of subset) per a parallel simulation.

According to the results for all input motions when the failure index is MDD, the occurrence frequency of the weakest floor is large at 1st floor: its ratio is 46.5%. The following ratios are 13% (9th floor), 12% (7th floor) and 11% (3rd floor) based on calculations of 600 times in MCS. While, when the failure index is MCDR, there was no occurrence of the weakest floor at 1st floor. The reason is that MCDR leads to relatively small since the structural height at 1st story is remarkably larger than other floor as shown in Table 2. Therefore, the results on failure index MDD will be hereinafter discussed mainly.

Table 4 shows that the combination of parameters is 1st and 2nd story stiffness for label “a”, and is 2nd story stiffness and story yield strength for label “b”.

Figure 6 shows sample parameters generated in each subset level. Here, the center of ellipse is the mean value μ of parameter, and the concentric ellipses are μ ± 1σ ~ 5σ. Subset-0 means MCS. Analytical cases name, e.g., Xcase02-b means that X is analytical frame, case02 is No.2 input motion, and “b” is the label of combination. Green


Figure 6. Scatter diagram of sample parameters on each subset level. (a) X-case02-b; (b) X-case06-b.

Table 4. Combination of parameters.

marker samples are within 10% of smaller values of limit state function Z. Paying attention to green makers of this figure, values of vertical axis SSy are smoothly small as subset is proceeding, however values of horizontal axis SK2 does not almost change. Thus, the effect for sensitivity of two parameters is that SSy is greater than SK2. This tendency is same for other input motions as shown in Table 5. In the case of failure index MCDR, the results of SSy > SK2 are obtained for other ground motions except No.1, 4, 5, 6. While, sensitivity results in combination of parameter “a” are almost SK2 > SK1 as shown in Table 5. Thus, it is seemed that relative intensity of sensitivity is SSy > SK2 > SK1.

5.2. Quantitative Estimation of Failure Probability

Firstly, let’s explain design criteria assuming in this study. It is general that design criteria are presented by interstory drift ratio. However, results of failure index MCDR are not valid because there is no occurrence of the weakest floor at 1st floor as mentioned earlier for this analytical model. So, it is assumed that design criteria presented by MCDR can be replaced by the criteria of MDD. Table 6 shows the relationship between MDD and MCDR at termination of subset-4 picked up the analyticcal cases with the common collapsed story. The mean value of MDD/MCDR (%) is about 2.2 as shown in this table.

Table 7 is an example of design criteria. After given criteria on MCDR as shown in the left column, MDD criteria result in 2.2 (the above mean value) × MCDR criteria. According to this table, the values of MDD, i.e., the maximum ductility demand, are over the definition of limit state shown in the right column. In the case of not high rise buildings like this analytical model, it is known that the response analysis gives fairly severer results than actual damages. It is seemed that its sign also appeared in the analytical results of this study.

Figures 7 (a)-(c) show relationship between ductility demand and failure probability and its ensemble mean value (black line) derived from parallel simulations. Each line shows the result of MCS (blue line), subset-1 (green line), subset-2 (red line), subset-3 (cyan line) and subset-4 (magenta line) in order from left side. Iteration of parallel simulation is 50 times, i.e., response analyses of 150,000 times (=3000 × 50), per a ground motion.

Figure 7(d) shows the convergence index by Gelman and Lubin to check whether the iteration is enough. value is pointed out at least within 1.1 ~ 1.2 according to the convergence judgment [4]. Iteration of 50 times satisfies this judgment as shown in Figure 7(d) or Table 8.

MDD of No.2, i.e., Kushiro offshore 2004 EW, is the largest in all ground motions, and No.8, i.e., JMA Kobe 1995 NS, follows next. However, note that this is the results subject to level 2 ground motion. According to Figure 7, the ensemble mean of MDD at termination of subset-4, which corresponds to failure probability = (1/601) × (0.1)4 = 1.664 × 10−7, is over 4.4 of the life safe. Now, if the annual allowable risk for residential buildings is 1 × 10−6, and level 2 ground motion has 10% exceedence in 50 years, which means that occurrence probability is 2.1 × 10−3 (=1/475, i.e., return period 475 years), the value of vertical axis of Figure 7 corresponding to the above annual risk is 4.76 × 10−4 (=1 × 10−6/2.1 × 10−3). The ductility demand corresponding to 4.76 × 10−4 is 4.3 for No.2 ground motion, and is 3.9 for No.8, these values are slightly smaller than 4.4 of life safe. However, taking into account that the response analysis gives fairly severer results than actual damages as mentioned earlier, it seemed that there is still some margin for the life safe.

Table 5. Relative intensity of sensitivity (failure index: MDD).

Table 6. Relationship between MDD and MCDR.

Table 7. Design criteria.


Figure 7. Failure probability and convergence index by Gelman-Lubin. (a) No.2; (b) No.7; (c) No.8; (d) R2: convergence index.

Table 8. Convergence index and computation time of parallel simulation.

As shown in Table 8, the computation time using note book computer: Dell Vostro 3500 with 32 bit operating system is very long for ground motion with a long cutoff time. In No.2 with cutoff time of 66.44 (sec), the computation time per a response analysis was 8.5267 (sec).

6. Conclusions

Using subset MCMC, the quantitative damage estimation of an actual reinforced concrete building with 11 stories subject to level 2 ground motion was investigated. Results are as follows:

1) It is seemed that the relative intensity of sensitivity is SSy (story yield strength) > SK2 (2nd story stiffness) > SK1 (1st story stiffness).

2) If an annual risk for residential building is 1 × 10−6, the ductility demand for allowable risk of No.2 ground motion, Kushiro offshore 2004 EW, which gave the largest damage to the analytical building, does not yet reach to the limit state of life safe.


  1. S. K. Au and J. L. Beck, “Subset Simulation and Its Application to Seismic Risk Based on Dynamics Analysis,” Journal of Engineering Mechanics, Vol. 129, No. 8, 2003, pp. 901-917. doi:10.1061/(ASCE)0733-9399(2003)129:8(901)
  2. W. R. Gilks, S. Richardson and D. J. Spiegelhalter, “Introducing Markov chain Monte Carlo,” In: W. R. Gilks, S. Richardson and D. J. Spiegelhalter, Eds., Markov Chain Monte Carlo in Practice, CRC Press, 1996, pp. 1-19.
  3. A. E. Raftery and S. M. Lewis, “The Number of Iterations, Convergence Diagnostics and Generic Metropolis Algorithms,” 1999.
  4. A. Gelman, “Inference and Monitoring Convergence,” In: W. R. Gilks, S. Richardson and D. J. Spiegelhalter, Eds., Markov Chain Monte Carlo in Practice, CRC Press, 1996, pp. 131-143.