Fragility and Vulnerability Analysis of an RC Building with the Application of Nonlinear Analysis

Abstract

In this paper, the seismic response of a five-story reinforced concrete (RC) frame system building is analysed through fragility analysis. The structure is designed in accordance with structural Eurocodes EN1990, EN1991, EN1992 and EN1998, as a high-ductility (DCH) system. For the analysis of the response of a structural system to earthquake actions, the methods of nonlinear static (NSA) and nonlinear dynamic analyses (NDA) are applied and, based on the obtained results, fragility curves are constructed using statistical methods. A relationship between the intensity measure (IM) and engineering demand parameters (EDPs) is needed in order to estimate a fragility curve. Fragility functions represent a possibility for different states of damage to occur in a certain structural systems at the observed value of the specified IM. Ten accelerograms, used in NDA, are selected and scaled, according to EN1998 provisions, for the chosen elastic response spectrum and referent PGA. Obtained results are used for the statistical analysis and construction of fragility curves. Structural damage state threshold parameters are determined based on the Park and Ang modified damage index methodology and provisions given in FEMA, HAZUS, VISION 2000 and EN codes. Comparative analysis of the structural damage probability for the analysed RC building, calculated using different methodologies to determine damage states, is applied. The fragility analysis results showed the difference between the probabilities of the damage states to occur when different calculation methods are used. This reflects on the assessment of vulnerability curves as well. The obtained results, calculated using different methods, are analysed using comparative analysis.

1. Introduction

During high-intensity earthquake actions and/or incident actions, the response of RC structures is accompanied by the occurrence of cracks and various kinds of damage. Contemporary concepts of structural design, which employ nonlinear analysis methods, facilitate a better analysis of structural behaviour after the failure and transition into the post-elastic behaviour domain.

Designs of seismically resistant structures, designs that withstand seismic action, that use the controlled reduction of the system bearing capacity necessitate providing a nonlinear ductile response of the system. This results in the occurrence of damage to the structural and non-structural elements and installations during the nonlinear behaviour of the structure [1]. Depending on the degree of structural damage, there is a risk of considerable material losses, compromising of normal services of buildings and even the loss of human life. Even the area outside of the building, in its immediate vicinity, can be affected by the damage or collapse of the structure. For that reason, there are has been research conducted in recent years, striving to provide the probability of the occurrence and degree of structural damage and loss. To that end, fragility and vulnerability functions are used to facilitate the analysis of building damage, collapse risk and loss during seismic and other incident actions.

The design concepts of seismically resistant structures and the theoretical basis of the problem, as well as the system response analysis, have been treated in many papers. The review of the implementation of the methodology for the assessment of structure designs based on the seismic performance is the subject of the paper by Günay and Mosalam [2]. Challenges and progress in performance-based earthquake engineering—PBEE—have been presented in the paper by Krawinkler [3]. In addition, a detailed review of PBEE is the subject of papers by Moehle and Deierlein [4], Porter [5] and Pristley [6].

In the domain of seismic engineering, the fragility functions can be used to calculate the possibility for different states of damage to occur in certain structures at an observed value of a specified intensity measure. According to Porter [7], one can define the fragility function as a mathematical function that expresses the probability that some undesirable event occurs as a function of some measure of environmental excitation. A fragility functions represent the cumulative distribution function of the capacity of an asset to resist an undesirable limit state.

The structural performance and damage assessment of structural systems using fragility curves are presented in papers by Borele and Datta [8], Čokić et al. [9], Čokić et al. [10] and Folić and Čokić [11]. Vulnerability and risk assessment of the RC frame structure were calculated and analysed in a paper by Olteanu [12], in which the methodology for the damage state assessment proposed in RISK-UE, by Milutinović and Trendafiloski [13], was applied. The fragility of tall RC buildings was analysed in the doctoral dissertation of Pejović [14] and the paper by Pejović and Janković [15]. Seismic fragility estimates for reinforced-concrete-framed buildings were the topic of a doctoral dissertation by Ramamoorthy [16]. Many other authors investigated the fragility assessment of RC buildings. A review of fragility curves for seismic risk assessment of buildings in Europe was done by Maio et al. [17] In their papers, Vazurkar and Chaudhari [18], Yeh et al. [19] and Zemeeruddin and Sangle [20] focused on the development of fragility curves and damage assessment of RC buildings as well. Fragility curve assessment methods and the theoretical basis and clarification of statistical methods used in the process are presented in the papers by Porter [7], Baker [21] and Sfahani et al. [22]

The development of seismic vulnerability assessment methodologies over the past 30 years is described in a paper by Calvi et al. [23]. A review of the state of assessment of seismic vulnerability caused by natural disasters is the subject of a paper by Vamvatsikos et al. [24] Giovinazzi [25] proposed the method for the vulnerability assessment and the damage scenario in seismic risk analysis. Seismic vulnerability assessment, damage scenarios and loss estimation are investigated in the research of Vicente et al. [26] Damage assessment of current buildings on a territorial scale is the subject of the research by Giovinazzi and Lagomarsino [27] and Lagomarsino and Giovinazzi [28], and an overview of research activities and present assessments with guidelines for the future were considered in a paper by Atalić et al. [29]

Fragility and vulnerability analyses of RC-framed buildings, considering the effects of their structural characteristics, were performed by Olteanu et al. [30] Medium- and high-rise structures were considered, with four and eight stories, in order to determine the building height in the vulnerability index. Three plan configurations were compared: square, rectangular and L. It was concluded that the structures with fewer levels were more vulnerable than the structures with more levels. The structures with the L-shaped building plan configuration were more vulnerable than the ones that were regular in plan. In the research conducted by Pejović and Janković [15], 20-story, 30-story and 40-story RC high-rise buildings with core wall structural systems were chosen. The derived fragility curves, are shown separate for each considered prototype RC high-rise building (20-story, 30-story and 40-story), leading to the conclusion that fragility increases with a building’s height, but with no significant differences in fragility for certain heights of the building.

In this paper, fragility curves for different types of damage (in term of severity) are calculated based on the relationship of the intensity measure (IM), represented through values of peak ground acceleration (PGA), and the engineering demand parameters (EDPs), which are represented through the response of the analysed structure. The limit state (LS) values used in this paper as EDPs, for the determination of structural damage, are functions of inter-story drift (IDR) [31,32,33], damage index (DI) [34] and PGA [35,36]. The seismic structural response is analysed using the pushover nonlinear static analysis (NSA) and nonlinear dynamic analysis (NDA) methods. Fragility curves were constructed using mathematical statistic and probability methods, based on the obtained results and IM–EDP relationship.

Fragility curves are calculated for a five-story reinforced concrete (RC) building that exhibits the properties of the frame structural system [35]. The structural system was designed according to the set of structural Eurocodes [35,37,38,39]. Based on the analysis, obtained results are compared and final remarks and conclusions are formulated.

In this paper, a comparative analysis of the performance of the RC frame structure is performed from the perspective of the application of different methods for the assessment of the damage states and therefore their fragility and vulnerability. Since the procedure for the calculation based on the damage index assessment method is more complex in relation to the other methods applied in this paper, it was used as a reference value for comparing the obtained results. Based on the results obtained by a fragility assessment, a vulnerability calculation was performed with a comparative analysis of the results, using the results determined by the damage index calculation procedure as a reference method.

2. Methodology of the Analysis and Structural Modelling

2.1. Geometric and Material Properties of the Structure

The subject of the analysis is an office/residential building (Figure 1) with 5 levels (ground floor + 4 stories). The structural system exhibits the properties of a frame structural system [35]. The plan view and the 3D model of the structure are shown in Figure 1. The length of one span in both directions is 4.8 m, which makes the total length of the building 19.2 m in both directions. The height of the first story is 3.6 m and the heights of the other stories are 3.2 m, which makes the total height of the building 16.4 m. In order to simplify the modelling and calculation process, all vertical elements are fixed at the bottom level of the structure, i.e., soil–structure interaction is not included in the calculation and design.

The design of the structure is done according to the recommendations given in a set of structural Eurocodes [35,37,38,39] using linear elastic analysis methods. Material properties of concrete C35/45 and reinforcing steel class C (f_yk = 500 MPa and E_y = 200 GPa) [39] have been adopted for model analysis.

The structural design is done according to European building design standards, and as such the structure displays high-ductility (DCH) behaviour [35]. The calculations are performed using ETABS [40]. Geometric and reinforcement characteristics of the cross-section properties of the beams and columns are shown in Figure 2.

2.2. Loads and Actions

The loads acting on the structure are as follows: permanent loads (G_i)—self-weight of structural elements and an additional permanent load; the variable live load (Q_i); and the seismic load (S_i). The adopted value of the permanent constant load is g_pl = 3.0 kN/m² and the load intensity of the variable live load amounts to q = 2.0 kN/m² on all floors except the roof, where it amounts to q_R = 1.0 kN/m². The self-weight load of façade elements, which is imposed on all façade elements except on the roof, is equal to g_f = 10.0 kN/m. The value of the reduction factor of the live loads is ψ_2,i = 0.3 [38]. To calculate the impact of an earthquake on the structure, an elastic response spectrum, type 1, was used for the ground type C, with a PGA of a_g = 0.2∙g. The behaviour factor q of the design response spectrum for a frame structural system is equal to 5.85 [35].

2.3. Modal Analysis

The modal analysis was performed to determine the fundamental periods of vibration of the system, modes (Figure 3) and system rigidity or flexibility [35]. It was established that the system was torsionally rigid and that translational modes were dominant. Rayleigh viscous (mass–tangent stiffness) (Figure 4) proportional damping was used in THA. The damping matrix of the system is a combination of the mass and stiffness matrices, as shown in the following equation:

$$\left[C\right]={\alpha }_M·\left[M\right]+{\alpha }_K·\left[K_T\right]$$

(1)

where [C], [M] and [K_T] are the damping, mass and tangential stiffness matrices, respectively. Parameters α_M and α_K represent proportional coefficients of the damping of mass and stiffness, and they are equal to:

$${\alpha }_M=4\pi ·\frac{{\xi }_{1,i}·T_{1,i}-{\xi }_{2,i}·T_{2,i}}{T_{1,i}^2-T_{2,i}^2}, {\alpha }_K=\frac{T_1·T_2}{\pi }·\frac{{\xi }_{2,i}·T_{1,i}-{\xi }_{1,i}·T_{2,i}}{T_{1,i}^2-T_{2,i}^2}$$

(2)

where T_1,i and T_2,i are the first and last periods of vibration of interest for the analysis of structural response. ξ_1,i and ξ_2,i are the corresponding relative damping coefficients with the adopted values of 0.05. The Rayleigh damping function is calculated by applying the expression:

$${\xi }_i={\xi }_{M,i}+{\xi }_{K,i}=\frac{{\alpha }_M}{2·{\omega }_i}+\frac{{\alpha }_K·{\omega }_i}{2}=\frac{{\alpha }_M·T_i}{4·\pi }+\frac{{\alpha }_K·\pi }{T_i},\left({\omega }_i=\frac{2\pi }{T_i}\right)$$

(3)

where ω_i is the angular frequency for the corresponding inherent form of vibration. An overview of modern seismic analyses with different damping models is explained in Ćosić et al. [41]

The appropriate values of periods of vibration in both directions for T₁ correspond to the first translation period in X or Y direction (the structure is orthogonal in the plane and its behaviour will be the same in both main directions). The value of T₂ corresponds to the period of vibration in which the structure reaches at least 90% of the sum of the effective modal masses [35]. Values of the used periods are T₁ = 0.735 s (80.41% mass) and T₂ = 0.213 s (92.69% mass in sum).

2.4. Nonlinear Dynamic Analysis

THA is conducted by using ten particular accelerograms in only one direction at a time and not by combining their N-S, W-E and vertical components for 3-directional analysis, because of its simplicity and the wide application of this method of analysis. The accelerograms (Figure 5) were chosen from the European Strong-Motion Database (ESMD) [42] and ORFEUS database [43] (

Table 1. Main properties of the earthquakes that were used in NDA.

ID	Earthquake Location	Earthquake ID (Component/Orientation)	Station ID/Code	Date/Time	MW	Original PHA (cm/s2)
EQ01	Spitak, Armenia	213 (Y)	173	07/12/1988 07:41:24	6.7	179.580
EQ02	Manjil, Western Iran	230 (Y)	189	20/6/1990 21:00:08	7.4	87.045
EQ03	Umbria Marche, Central Italy	286 (Y)	221	26/9/1997 09:40:30	6.0	218.340
EQ04	Umbria Marche, Central Italy	286 (Y)	224	26/9/1997 09:40:30	6.0	106.660
EQ05	Alkion, Greece	559 (X)	214	15/6/1995 00:15:51	6.5	55.501
EQ06	Düzce, Turkey	497 (Y)	3139	12/11/1999 16:57:20	7.2	112.320
EQ07	Umbria, Central Italy	EMSC-20161030_0000029 (N-S)	CNE	30/10/2016 06:40:18	6.5	288.280
EQ08	Emilia-Romagna, Italy	IT-2012-0011 (N-S)	MOG0	29/5/2012 07:00:02	6.0	167.075
EQ09	Adana, Turkey	TK-1998-0063 (E-W)	0105	27/6/1998 13:55:53	6.2	271.955
EQ10	Emilia-Romagna, Italy	IT-2012-0011 (N-S)	MIR08	29/5/2012 07:00:02	6.0	242.970

) and scaled according to elastic response spectra for the intensity level of 10% possibility of exceedance in 50 years.

EQ01–EQ06 are selected from ESMD [42] and EQ07–EQ10 are selected from ORFEUS [43]. The criteria for time history data selection was that magnitude, M, was >5.5 Ms (type 1 RS [35]), the records correspond to soil type C and v_s,30 = 180–360 m/s [35]. In addition, records from ORFEUS [43] are obtained using the REXELite tool that allows searching for a suite of waveforms compatible with a target spectrum generated according to [35]. The main criterion for the selection of earthquake records was that the mean of the selected and then scaled records matched the elastic EN1998-1 [35] response spectrum used in the analysis. According to EN1998-1, [35] the conditions $S_{a,\mu }^{scaled}\left(T_0\right)\ge a_g·S$ and $S_{a,\mu }^{scaled}\ge 0.9·S_{a,EC8}^{Elastic}$ on the interval of $\left[0.2·T_1-2·T_1\right]$ are satisfied (Figure 5). Selected earthquake data are shown in Table 1. THA data was scaled (Figure 5) with the common scale factor F_S = 1.61, which was obtained using the least square method (LSM). With this approach, the matching of the calculated average spectrum to the targeted one is quite good, since all-time histories are amplified by the same factor [44]. The selection and scaling process of the earthquake ground motions for THA is also described in NIST GCR [45] and Bisch et al. [46] Scaled accelerograms are used for NDA, with the increment of ΔPGA = 0.1 g, in a total scaling factor range of 0.1–1.0 g.

3. Structural Model

3.1. Model for Linear Elastic Analysis

For the calculation and design of the structure in ETABS [40], a spatial (3D) model was used. The following parameters, assumptions and simplifications were adopted:

• The calculation includes the effects of second-order logic (P-Δ);
• Occurrence of cracks in structural elements was included in the calculation with the stiffness reduction of the elements according to EN1998-1 [35];
• The elastic bending stiffness and shear stiffness of columns and beams was reduced to 50%;
• Torsion stiffness of columns and beams were reduced to 10% of their elastic stiffness;
• The elastic bending stiffness of the RC plate was reduced to 50%.

3.2. Model for Nonlinear Analysis

In models for post-elastic analysis of the structural response to the removal of individual vertical elements, the following assumptions and simplifications were used:

• The calculation includes the effects of the second-order logic (P-Δ);
• To describe the nonlinear behaviour of the material, the nonlinear properties of the material were used to describe the behaviour of concrete [39,47] and reinforcement steel (Figure 6) [39,48];
• The behaviour of RC is described by the Takeda hysteresis model, while the kinematic model of hysteresis was used for reinforcement. Both models are an integral part of the software package ETABS [40];
• Parameters describing the appearance of cracks as a result of elastic bending stiffness in structural elements from the linear elastic model were not included in the nonlinear model for P1–M2–M3 interaction, because plastic hinges are modelled as fibre elements, whereas the properties of fibres are described by stress–strain relations in concrete and reinforcement steel (Figure 6);
• Columns and beams were modelled as confined RC elements with a protective layer of concrete [47];
• The beams are modelled as “L” and “T” cross-sections, with the effective width of the RC plate;
• RC plates are modelled as rigid diaphragms.

Characteristics of confined concrete in structural elements are calculated after Mander et al. [47] The value of the stress in concrete is determined by the expression:

$$f_c=\frac{f_{cm,c}·x·r}{r-1+x}$$

(4)

where:

$$x=\frac{{\epsilon }_c}{{\epsilon }_{cm,c}}$$

(5)

$f_c$ is the compressive stress in concrete. $f_{cm,c}$ is the maximum value of compressive strength of confined concrete and ${\epsilon }_{cm,c}$ is the corresponding strain value. ${\epsilon }_{cm}$ equals:

$${\epsilon }_{cm}={\epsilon }_{c2}·\left[1+5·\left(\frac{f_{cm,c}}{f_{cm}}-1\right)\right]$$

(6)

where $f_{cm}=43$ Mpa (C35/45) is the maximum value of compressive strength of unconfined concrete and ${\epsilon }_{c2}=2‰$ [39] is the corresponding strain value. $E_c$ [39] and $E_{sec}$ are tangent and secant modules of the elasticity of concrete in equations:

$$r=\frac{E_c}{E_c-E_{sec}}$$

(7)

and

$$E_{sec}=\frac{f_{cm,c}}{{\epsilon }_{cm,c}}$$

(8)

The value of the confined concrete core strain ${\epsilon }_{cu,c}$ is calculated after the expression provided in Paulay and Pristley [49]:

$${\epsilon }_{cu,c}=0.004+1.4·\frac{{\rho }_{yh}·f_{yh}·{\epsilon }_{su}}{f_{cm,c}}$$

(9)

where the volume coefficient of reinforcement steel, ${\rho }_{yh}$, is calculated according to the simple expression, which also favours safety because it adopts the minimum value of the transverse reinforcement coefficient in one of two principal directions $\left({\rho }_{w,x},{\rho }_{w,y}\right)$ [50]:

$${\rho }_{yh}=2·\min \left({\rho }_{w,x},{\rho }_{w,y}\right)$$

(10)

where $f_{yh}$ represents the yield stress of confinement reinforcement, while ${\epsilon }_{su}$ is the ultimate value of the strain of confined concrete. The calculation of ${\rho }_{yh}$ as a square root of the product of the squares of ${\rho }_{w,x}$ and ${\rho }_{w,y}$ has been considered, but Expression (10) has been adopted because it favours safety [50]. Figure 6 shows functions of the stress–strain ratio for the description of nonlinear behaviour of confined and unconfined concrete.

Characteristics of reinforcement steel for the nonlinear analysis are shown in Figure 6. Characteristic values corresponding to class C reinforcement steel are adopted. Values of stress and strain in steel at the yield point are $\left({\sigma }_{sy}=575 \mathrm{MPa}, {\epsilon }_{sy}=0.2875‰\right)$, while these values at the point of ultimate strain are $\left({\sigma }_{su}=660 \mathrm{MPa}, {\epsilon }_{su}=7.5‰\right)$. The entire procedure of the design of the reinforcement model for the nonlinear analysis is described in Annex E of EN1998-2 [48], where the stress in steel for nonlinear analysis is adopted according to the formula:

$$\frac{f_{ym}}{f_{yk}}=1.15$$

(11)

where $f_{yk}$ is the value of steel at the yield point for the linear analysis and $f_{ym}$ is the value of steel at the yield point for the nonlinear analysis. The ratio of the stress at the yield point, ${\sigma }_{sy}$, and the ultimate stress, ${\sigma }_{su}$, in steel is $k=1.15=\frac{{\sigma }_{su}}{{\sigma }_{sy}}$, while the strain at the ultimate stress is ${\epsilon }_{su}=7.5‰$ [39]. Figure 6 shows the function of the ratio of stress and strain for the nonlinear behaviour of reinforcement steel.

3.3. Properties of Plastic Hinges

Plastic hinges are modelled as fibre cross-sections. They are modelled by an automatic selection of fibre division in the cross-section of elements in ETABS [40]. Among the many expressions for the calculation of plastic hinge length [51], and because of the inconsistencies among the values obtained by different expressions, Equations (12) and (13), suggested by Park and Paulay [52] and Priestley et al. [53], are the most practical for the modelling and the analysis. A schematic presentation of plastic hinge on a column element is displayed in Figure 7.

$$L_{pl}=0.5·h$$

(12)

$$L_{pl}=0.08·L+0.022·f_y·d_{bL}\ge 0.44·f_y·d_{bL}$$

(13)

where h is the height of the cross-section of the element, $L_{pl}$ is the length of the plastic hinge, $L=M/V$ (shear length), $f_y$ is the stress at the yield point and $d_{bL}$ is the mean value of the longitudinal reinforcement profile diameter.

It is estimated that the lengths of plastic hinges, calculated according to Park and Paulay [52] and Pristley et al. [53], correspond approximately to the relative lengths of columns and beams of 0.1 L, where L is the length of the element. Therefore, the locations of the hinges are assigned as 0.05 L and 0.95 L to columns and beams in ETABS [40]. Behaviour of the fibre plastic hinges used in the calculations is a function of an auto-discretised group of fibres and their stress–strain relationship for the materials which are used in the formation of the structure’s cross-sections (Figure 8).

4. Nonlinear Analysis Results and Calculation of Fragility and Vulnerability Curves

4.1. NDA

The results of the NSA and NDA for mass proportional and modal load distributions are shown in Figure 9. The pushover curve for the modal load distribution has a better fit to NDA values, so it was chosen as a referent curve for the calculation of fragility curves, using DI_PA^MOD as an EDP, according to Ghosh et al. [54] To perform the calculation of DI_PA, it was necessary to conduct a bilinear approximation of the NSA pushover curve, using a modified equivalent energy elastic-plastic (EEEP) method (Figure 9, right), and determine yielding (d_Y,F_Y) = (8.19 cm, 8283.23 kN) and ultimate capacity (d_U,F_U) = (28.39 cm, 8283.23 kN) points.

Displacement (d), IDR and their mean values, obtained using NDA, are displayed in Figure 10.

4.2. Damage State Performance Points

Damage of a structural system may be quantified through threshold performance points, which represent the values of EDPs, which are obtained by NSA and NDA for the damage index (DI) calculations and by NDA for other described methods in this paper. There are several methods with which to define damage state performance points. The ones described in [31,32,33,34,35,36] were used to calculate fragility curves.

A modified Park and Ang damage index [34], proposed by Ghosh et al. [54], is calculated by the following equation:

$$DI=\frac{D_M-D_Y}{D_U-D_Y}+\beta ·\frac{1}{F_Y·D_U}·\underset{E\left(t=0\right)}{\overset{E\left(t=n\right)}{{\displaystyle \int }}}dE\ge 0 ; (\mathrm{if} F_U<F_Y, F_Y \mathrm{is} \mathrm{replaced} \mathrm{by} F_U)$$

(14)

where $D_M$ represents the maximum roof displacement under an earthquake in THA; $D_U$ represents the ultimate roof displacement on the bilinearized pushover curve (Figure 9, right); $F_Y$ represents the yield roof displacement on the bilinearized pushover curve (Figure 9, right); and ${{\displaystyle \int }}_{E\left(t=0\right)}^{E\left(t=n\right)}dE$ represents the absorbed hysteretic energy during the earthquake (from $t=0$ to $t=n$, where $n$ is the last integration point of the accelerogram). The value of $D_Y$ is adopted as the first yield displacement from THA, at the point of the first non-zero value of $E_h$. $\beta$ is a non-negative parameter representing the effect of cyclic loading on structural damage, usually equal to 0.15 for RC structures [55]. A more accurate and complex modified DI model was proposed and discussed in Lađinović et al. [55]

Hysteretic energy (E_h) absorbed during the earthquakes is displayed in Figure 11, left, and the first yield displacement values (d_Y) obtained using NDA are shown in Figure 11, right.

The relationship between IDR and roof displacement, d, is shown in Figure 12, left. Based on the NDA results, it was also possible to establish the relationships between DI and d (Figure 13, left), as well as between DI and IDR (Figure 13, right). These relationships were later used to calculate and compare damage state thresholds as a function of IDR in

Table 2. Structural damage state threshold values.

Damage Index	SD (DS1)	MD (DS2)	ED (DS3)	CD (DS4)
DI [54]	0.100	0.250	0.400	1.000
IDR (%)	0.355	0.679	1.008	2.365
HAZUS [32]	SD (DS1)	MD (DS2)	ED (DS3)	CD (DS4)
IDR (%)	0.333	0.600	1.533	4.000
VISION 2000 [33]	FO (DS1)	OP (DS2)	LS (DS3)	NC (DS4)
IDR (%)	0.200	0.500	1.533	4.000
FEMA [31]	IO	LS		CP
IDR (%)	1.000	2.000		4.000
EN [35]	DL	SD		NC
${\gamma }_{DS}$	0.780	1.000		1.735
$PG{A}_g^{DS}$ (g)	${\gamma }_{DL}·PG{A}_{EL}=0.156$	${\gamma }_{SD}·PG{A}_{EL}=0.2$		${\gamma }_{NC}·PG{A}_{EL}=0.347$
IDR (%)	0.380	0.508		1.050

.

In terms of the severity of the structural damage, structural damage states are defined as:

• Modified DI: SD—slight damage, MD—moderate damage, ED—extensive damage and CD—complete damage [54];
• HAZUS: SD—slight damage, MD—moderate damage, ED—extensive damage and CD—complete damage [32];
• VISION 2000: FO—fully operational, OP—operational, LS—life safe and NC—near collapse [33];
• FEMA: IO—immediate occupancy, LS—life safety and CP—collapse prevention [31];
• EN1998: DL—damage limitation, SD—significant damage and NC—near collapse [35,36].

Results displayed in Table 2 are shown in Figure 14 for a better visual comparison of damage state thresholds from the perspective of IDR.

The highest similarity between damage state thresholds is noticeable in the case of the HAZUS [32], VISION 2000 [33] and modified DI [54] models for damage assessment. The use of the FEMA 356 [31] methodology in damage state threshold calculations will give the results which are the least conservative, and the use of the methodology described in Eurocode 8 [35,36] will result in the most conservative values, which are on the side of safety.

4.3. Statistical Analysis of the Results

It is generally assumed that a fragility curve is a lognormal distribution function, which means that “If a variable is log-normally distributed, its natural logarithm is normally distributed. Which means it must take on a positive real value, and the probability of it being zero or negative is zero.” [7]. Using Kolmogorov–Smirnov and Anderson–Darling tests in EasyFit [56] on the results obtained through the NDA method, it was established that for each IM and EDP distribution, the values fit the lognormal distribution. That also means that the relationship between ln(IM) and ln(EDP) has a normal distribution for each ln(IM) in log–log space (l.l.s.). If each probability density function (PDF) in log–log space can be described as $Y~N\left(\mu ,{\sigma }^2\right)$, and each PDF in arithmetic space (a.s.) as $X~\ln N\left(\mu ,{\sigma }^2\right)$, conversion of l.l.s. values to a.s. can be performed by utilising the following expressions:

$${\mu }_{a.s.}=e^{\left({\mu }_{l.l.s.}\right)} ; {\sigma }_{a.s.}^2=\left|e^{{\sigma }_{l.l.s.}^2}-1\right|·e^{\left(2·{\mu }_{l.l.s.}+{\sigma }_{l.l.s.}^2\right)}$$

(15)

where ${\mu }_{a.s.}$ and ${\sigma }_{a.s.}$ represent the mean and the standard deviation values of the variable in arithmetic space, while ${\mu }_{l.l.s.}$ represents the mean and median values and ${\sigma }_{l.l.s.}$ represents the standard deviation value of the variable in the log–log space.

4.4. Calculation of Fragility Curves

In the case of the calculation of fragility curves, an EDP is used as a referent value in order to calculate the probability of the occurrence of a defined damage state at a particular intensity measure value $\left(P_{D{S}_i|I{M}_j}\right)$, using the expression:

$$P_{D{S}_i|I{M}_j}\left({\mu }_{LN,I{M}_j}^{ED{P}_i},{\sigma }_{LN|I{M}_j}^{ED{P}_i}\right)=1-\Phi \left(\frac{\ln EDP-{\mu }_{LN,I{M}_j}^{ED{P}_i}}{{\sigma }_{LN|I{M}_j}^{ED{P}_i}}\right)$$

(16)

where ${\mu }_{LN,I{M}_j}^{ED{P}_i}$ and ${\sigma }_{LN|I{M}_j}^{ED{P}_i}$ are the mean and standard deviation in l.l.s. of the PDF of the variable $\ln EDP$ for a particular $\ln I{M}_j$ value. $\ln ED{P}_i$ is the lognormal value of a DS threshold. These fragility curves are calculated using methods of statistical analysis and probability described in Porter [7] and Baker [21]. Probability values are calculated at each $I{M}_j$ for each $D{S}_i$. When all the probability values are calculated, the set of obtained points is fitted for each DS by using the maximum likelihood estimation (MLE) method [7,21]. After applying MLE, mean and standard deviation $\left({\mu }_{LN|D{S}_i}^{IM},{\sigma }_{LN|D{S}_i}^{IM}\right)$ as a function of the IM can be obtained, for both log–log and arithmetic space. Because the relationship between PGA and IDR or d is calculated, the fragility curves for DSs according to EN1998-1 [35] and EN1998-3 [36] can be converted to functions which have PGA as the IM, using the described procedure.

Since it is established that the values fit the lognormal distribution for each distribution of IM or EDP, the PDF will be expressed with the equation:

$$f_{LN,D{S}_i}^{IM}=\frac{1}{IM}·\frac{1}{{\sigma }_{LN|D{S}_i}^{IM}·\sqrt{2\pi }}·e^{-\frac{{\left(\ln IM-{\mu }_{LN|D{S}_i}^{IM}\right)}^2}{2·{\sigma }_{LN|D{S}_i}^{IM}{}^2}}$$

(17)

where ${\mu }_{LN|D{S}_i}^{IM}={\mu }_{l.l.s,i}^{IM}$ and ${\sigma }_{LN|D{S}_i}^{IM}={\sigma }_{l.l.s.,i}^{IM}$ represent the median and the standard deviation value of the variable $\ln IM$ for each damage state (DS).

In the case of the calculation of fragility curves, using IM as a referent value for the DS threshold, the fragility functions are fitted using an analytical cumulative distribution function (CDF) for lognormal distribution:

$$P_{D{S}_i|IM}\left({\mu }_{LN|D{S}_i}^{IM},{\sigma }_{LN|D{S}_i}^{IM}\right)=\Phi \left(\frac{\ln IM-{\mu }_{LN|D{S}_i}^{IM}}{{\sigma }_{LN|D{S}_i}^{IM}}\right)$$

(18)

where Φ is the cumulative distribution function of the standard normal distribution.

For the calculation of fragility curves according to EN1998-1 [35] and EN1998-3 [36], the IM is equal to d. Because of the performed NDA for PGA = 0.1–1.0 g and PGA = PGA_DL, PGA_SD, PGA_NC, which is 120 THA, the probabilities of each calculated d_i for each IM_j = 0.1–1.0 g can be obtained for the referent PGA = PGA_DL, PGA_SD, PGA_NC, because the relationship between d and PGA is known. Using this relationship, the fragility curves calculated according to EN1998-1 [35] and EN1998-3 [36] can be obtained as a probability—PGA functions. Upper and lower boundaries of the confidence interval (CI) are calculated using the equation:

$$CI={\mu }_{LN}\pm Z· \frac{{\sigma }_{LN}}{\sqrt{n}}$$

(19)

where ${\mu }_{LN}$ and ${\sigma }_{LN}$ are mean and standard deviation in l.l.s. of PDFs. Z is the standard score or the Z-score of the normal distribution function in the l.l.s. and n is the number of samples, or in this case, the number of EDPs per IM, which is equal to the number of used accelerograms. Lower and upper bounds of the confidence interval are adopted as the 16th and 84th percentile of the normal distribution function in the l.l.s.

When all the values of the probability of the occurrence of a certain damage state are calculated $\left(P_{D{S}_i|I{M}_j}\right)$, the set of the obtained points is fitted for each DS using the Maximum Likelihood Estimation—MLE—method described by Porter [7] and Baker [21].

$$Likelihood={\displaystyle \prod }_{j=1}^m\left(\begin{array}{c}{n}_j\\ z_j\end{array}\right)\varphi {\left(\frac{\ln \left(I{M}_j/{\theta }_{D{S}_i}^{IM}\right)}{{\beta }_{LN|D{S}_i}^{IM}}\right)}^{z_j}{\left[1-\varphi \left(\frac{\ln \left(I{M}_j/{\theta }_{D{S}_i}^{IM}\right)}{{\beta }_{LN|D{S}_i}^{IM}}\right)\right]}^{n_j-z_j}$$

(20)

where m is the total number of IM levels, ${\theta }_{D{S}_i}^{IM}$ is the mean value of distribution in arithmetic space and ${\beta }_{LN|D{S}_i}^{IM}$ is the value of the standard deviation in logarithmic form of the lognormal cumulative probability distribution function of the fitted fragility curve for the reference value of the IM for the corresponding damage state, DS. The calculation of the parameters of fragility curves is achieved by maximizing the likelihood function (20). This can be achieved, making the procedure easier in the process, by maximizing the logarithm of the maximum likelihood function, so Expression (20) has the following form:

$$\left\{{\widehat{\theta }}_{D{S}_i}^{IM},{\widehat{\beta }}_{LN|D{S}_i}^{IM}\right\}=\arg \max {\displaystyle \sum }_n^{i=1}\left\{\ln \left(\begin{array}{c}{n}_j\\ z_j\end{array}\right)+z_j\ln \varphi \left(\frac{\ln \left(I{M}_j/{\theta }_{D{S}_i}^{IM}\right)}{{\beta }_{LN|D{S}_i}^{IM}}\right)+\left(n_j-z_j\right)\ln \left(1-\varphi \left(\frac{\ln \left(I{M}_j/{\theta }_{D{S}_i}^{IM}\right)}{{\beta }_{LN|D{S}_i}^{IM}}\right)\right)\right\}$$

(21)

An example of the process of assessing the fragility curve for a particular damage state of ED according to the HAZUS methodology [32] is displayed in Figure 15, Figure 16 and Figure 17. After NDA is performed, all the results represented as black dots (lnIDR,i and IDR,i) in Figure 15 and Figure 16 are used for further analysis. By using the EasyFit software package [56], it is determined that the data fits the lognormal distribution. Because of that, the ln values (lnIDR,i) of the original data can be represented in the log–log space using the normal distribution functions (Figure 15) or in the arithmetic space using the lognormal distribution (Figure 15). The PDFs in log–log and arithmetic space are displayed for each lnIMj or IMj near the results of the NDA. The straight blue dashed line represents the boundary above which ED occurs (Figure 15 and Figure 16). The probability (P) of ED to occur is equal to the area of each PDF at a corresponding PGA, above the mentioned boundary limit of ED (Figure 15 and Figure 16). These values are calculated using Equation (16) and displayed as the blue dots in Figure 17. After the probability values for each IMj are calculated, the fragility curve of ED is fitted for Function (18), using the MLE and Expression (21). The fragility curve is fitted and mean and standard deviation values are determined using Expression (21) for CDF (18) (Figure 17). Upper and lower boundaries of the CI are calculated using Equation (19). Calculated fragility curves are shown in Figure 18, Figure 19, Figure 20, Figure 21 and Figure 22.

Probability density functions for the occurrence of different states of damage are calculated using the equations described in research by Porter [7] and Baker [21]:

$$\begin{gathered} P_{D{S}_0}=1-P_{D{S}_1}\left[I{M}_j,{\mu }_{LN|D{S}_1},{\sigma }_{LN|D{S}_1}\right] \\ {P}_{D{S}_i}=P_{D{S}_i}\left[I{M}_j,{\mu }_{LN|D{S}_i},{\sigma }_{LN|D{S}_i}\right]-P_{D{S}_{i+1}}\left[I{M}_j,{\mu }_{LN|D{S}_{i+1}},{\sigma }_{LN|D{S}_{i+1}}\right] \\ {P}_{D{S}_n}=P_{D{S}_n}\left[I{M}_j,{\mu }_{LN|D{S}_n},{\sigma }_{LN|D{S}_n}\right] \end{gathered}$$

(22)

where $P_{D{S}_0}$ is the probability of no damage to occur and $i=1,\dots ,n$ and $I{M}_j=\left(0-1.0 \mathrm{g}\right)$. $i$ is an index of a particular DS, and $j$ is an index of a particular IM (PGA). $n$ is the total number of damage states.

The results of the analysis, the calculation of fragility curves and the probability of the occurrence of different states of damage with their upper and lower CI bounds are shown in Figure 18, Figure 19, Figure 20, Figure 21 and Figure 22.

The relationship between fragility curves fits the aforementioned relationship between DS threshold values. The biggest similarity between damage state thresholds is noticeable in the cases of the HAZUS [32], VISION 2000 [33] and modified DI [34] models for damage assessment. The use of the FEMA 356 [31] methodology in damage state threshold calculations will give the results which are the least conservative, and the use of the methodology described in Eurocode 8 [35,36] will result in the most conservative values, which are on the side of safety.

The comparison of fragility models for the structure in total can also be obtained through the assessment and analysis of vulnerability curves, which represent the cumulative distribution of the total repair cost of the structure. The transformation of the fragility curves into vulnerability curves can be conducted by using the following total probability relation, according to D’Ayala et al. [57]:

$$E\left(C|IM\right)={\displaystyle \sum }_{i=0}^{n}E\left(C|D{S}_i\right)·P\left(D{S}_i|IM\right)$$

(23)

where n is the number of damage states $\left(D{S}_i\right)$ considered, $P\left(D{S}_i|IM\right)$ is the probability of a building sustaining $D{S}_i$ given intensity, IM; $E\left(C|D{S}_i\right)$ is the complementary cumulative distribution of the cost (loss) given the $D{S}_i$; and $E\left(C|IM\right)$ is the complementary cumulative distribution of cost (or loss) given a level of intensity, IM [57]. The values of $E\left(C|D{S}_i\right)$ are adopted from the tables in the research of D’Ayala et al. [57] The values given in

Table 3. Damage factor functions of building typology, for the calculation of vulnerability curves, according to D’Ayala et al. [57].

Damage Scale	Damage Factor Functions of Building Typology $\mathit{E}\mathbf{\left(}\mathit{C}\mathbf{|}\mathit{D}{\mathit{S}}_{\mathit{i}}\mathbf{\right)}$
Damage State	Slight	Moderate	Extensive	Complete	Comments
	2%	10%	50%	100%	Country: Spain, France and General Europe. DF definition: Repair/replacement cost. Comment: Part of RISK-UE project. DFs mainly based on expert judgement and national “experience” with no stated method. Values are identical to HAZUS99.
Damage State	Green Tag	Yellow Tag		Red Tag	Comments
	9.8%	25.6%		100%	Country: Greece. DF definition: Repair cost/pre-earthquake market value. Non-structural damage is assumed to be included. Comment: Based on results of post-earthquake surveys carried out after the 1999 Athens earthquake. A “representative” sample of 150 building blocks, or 983 buildings, corresponding to 10% of the total building population were surveyed.

correspond to the fragility analysis damage states with four [32,33,54] and three states of damage [31,35,36].

The results of the analysis and the calculation of vulnerability curves and the mean damage factors (MDFs) of the occurrence of different states of damage are shown in Figure 23 and Figure 24.

Vulnerability curves and MDFs for PGA are displayed in Figure 24, left. Results in Figure 24, right, show the values of MDFs with their upper and lower CI bounds for the design PGA = 0.2 g. It is evident that, as expected, the relationship between the vulnerability curves is the same as the relationship between DS thresholds for the used fragility assessment methodologies.

5. Discussion of the Results

To compare the differences between the derived fragility curves and calculated vulnerability curves, the results obtained using the modified Park and Ang method [54] were used as referent values. This method was chosen because it is the most complex one among the selected methods and because it integrates the results of ultimate displacement and yield shear force at the base of the structure, obtained using the pushover NSA, and maximum roof displacement, absorbed hysteretic energy and yield displacement values, obtained using NDA, to calculate the damage index values of the structure. In Figure 25 and Figure 26, the differences between fragility curves calculated using other methods and the ones calculated using the modified DI method are presented for different damage states.

The differences are calculated using the equation:

$$\mathsf{\Delta }{P}_{DI|I{M}_j}^m=P_{D{S}_i|I{M}_j}^m-P_{DI|I{M}_j}$$

(24)

for fragility curves, and:

$$\mathsf{\Delta }MD{F}_{DI|I{M}_j}^m=MD{F}_{D{S}_i|I{M}_j}^m-MD{F}_{DI|I{M}_j}$$

(25)

where P represents the probability of the occurrence of a particular damage state, the index $D{S}_i$ represents a particular damage state, $I{M}_j$ is the value of the intensity measure, DI represents the damage index, MDF is the mean damage factor and m is the method used for the estimation of the fragility and vulnerability of the structure.

To compare the differences between the calculated vulnerability curves, the results obtained using the modified Park and Ang method [54] were used as referent values. In Figure 27, the differences between vulnerability curves calculated using other methods and the ones calculated using the modified DI method are presented.

The differences for the value of the design PGA = 0.2 g compared to the modified damage index are shown in

Table 4. The differences for the value of the design PGA compared to the modified DI.

Difference	Fragility Analysis				Vulnerability Analysis
Δ (−)	SD	MD	ED	CD	MDF
HAZUS	0.06	0.16	−0.07	−0.03	−0.03
V2000	0.25	0.26	−0.07	−0.01	−0.01
FEMA	−0.46	−0.21	−0.11	−0.03	−0.08
EC8	0.24	0.55	0.64	0.40	0.27

.

The results of the fragility analysis show the differences when different DS thresholds are used, which reflects on the assessment of fragility curves and the calculation of vulnerability functions as well. Fragility and vulnerability curves as well as DS occurrence probabilities calculated according to FEMA 356 [31] will result in the least conservative values, as opposed to the procedure according to Eurocodes [35,36], which gives the results that are the most on the side of safety among all four procedures. Fragility and vulnerability curves as well as DS occurrence probabilities calculated according to HAZUS [32], VISION 2000 [33] and the modified Park and Ang damage index [54] give similar results, and the values of vulnerability calculations can be described as $\mathrm{E}\left(\mathrm{EC}8\right)>\mathrm{E}\left(\mathrm{DI}\right)>\mathrm{E}\left(\mathrm{VISION} 2000\right)>\mathrm{E}\left(\mathrm{HAZUS}\right)>\mathrm{E}\left(\mathrm{FEMA} 356\right)$.

6. Final Remarks and Conclusions

In this paper, a comparative analysis of the seismic fragility and vulnerability of a RC frame structure, with damage state thresholds calculated according to different methods described in [31,32,33,34,35,36,54], is performed. The structure is designed as a DCH system, with behaviour factor q = 5.85. The system’s fragility curves were derived from the results of NSA and NDA, using statistical methods. Different EDP values were selected for the assessment of damage states, while PGA was selected as the intensity measure, IM. Probabilities of the occurrence of each DS for all mentioned codes were determined and compared.

Based on the comparative analysis, it can be concluded that based on the approach for the determination of EDP values, according to a different code, results may vary significantly in the range where it can be concluded that the structure is over-strengthened [31] or designed close to its full bearing capacities [35,36] for the design PGA. More descriptive and slightly more complex procedures are described in HAZUS [32], VISION 2000 [33] and Park and Ang [34] (Ghosh et al. [54]), where the fragility and vulnerability results are related to IDR [32,33] as well as displacement and hysteretic energy [54] values, are more in-between the values obtained by using methods described in FEMA 356 [31] as well as EN1998-1 [35] and EN1998-3 [36] and give the most satisfactory results in the fragility analysis.

The methodologies described in this paper can be applied to the fragility and vulnerability assessment of different types of buildings. This comparative study can be applied to RC-framed structures that are regular in plan and similar to the analysed model. For different types of structures, including masonry-infilled frames and RC load-bearing masonry wall system cases, additional analyses should be applied before making any conclusions related to the analysis conducted in this research.