# INNOVATIVE METHODOLOGIES FOR RESILIENT BUILDINGS AND CITIES

EDITED BY : Izuru Takewaki, Masayuki Kohiyama, Tomaso Trombetti, Solomon Tesfamariam and Xinzheng Lu PUBLISHED IN : Frontiers in Built Environment

#### Frontiers Copyright Statement

© Copyright 2007-2019 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use. ISSN 1664-8714 ISBN 978-2-88963-072-1 DOI 10.3389/978-2-88963-072-1

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# INNOVATIVE METHODOLOGIES FOR RESILIENT BUILDINGS AND CITIES

Topic Editors:

Izuru Takewaki, Kyoto University, Japan Masayuki Kohiyama, Keio University, Japan Tomaso Trombetti, University of Bologna, Italy Solomon Tesfamariam, University of British Columbia, Canada Xinzheng Lu, Tsinghua University, China

Resilient buildings and cities are in the center of common interests in modern academic communities for science and engineering related to built environment. Resilience of buildings and cities against multidisciplinary risks, e.g. earthquakes, strong winds, floods, etc., is strongly related to the sustainability of buildings and cities in which reduction of damage during a disaster and fast recovery from the damage are key issues. The reduction of damage is related to the level of resistance of buildings and the time of recovery is affected by the amount of supply of damaged members, assurance of restoration work, etc. Robustness, redundancy, resourcefulness, and rapidity are four key factors for supporting the full realization of design and construction of resilient buildings and cities. This research topic gathers cutting-edge and innovative research from various aspects, e.g. robustness of buildings and cities against earthquake risk, structural control and base-isolation for controlling damage risks, quantification of resilience measures, structural health monitoring, innovative structural engineering techniques for higher safety of buildings, resilience actions and tools at the urban scale, etc.

Citation: Takewaki, I., Kohiyama, M., Trombetti, T., Tesfamariam, S., Lu, X., eds. (2019). Innovative Methodologies for Resilient Buildings and Cities. Lausanne: Frontiers Media. doi: 10.3389/978-2-88963-072-1

# Table of Contents


Nurbaiah Mohammad Noh and Solomon Tesfamariam


Koki Makita, Kyoichiro Kondo and Izuru Takewaki

*73 Optimal Viscous Damper Placement for Elastic-Plastic MDOF Structures Under Critical Double Impulse*

Hiroki Akehashi and Izuru Takewaki *90 Influence of the Brace Configurations on the Seismic Performance of* 

*Steel Concentrically Braced Frames* T. Y. Yang, Hediyeh Sheikh and Lisa Tobber

*103 Experiments and Fragility Analyses of Piping Systems Connected by Grooved Fit Joints With Large Deformability*

Tao Wang, Qingxue Shang, Xi Chen and Jichao Li

*117 Application of Association Rules to Determine Building Typological Classes for Seismic Damage Predictions at Regional Scale: The Case Study of Basel*

Lorenzo Diana, Julien Thiriot, Yves Reuland and Pierino Lestuzzi

# Editorial: Innovative Methodologies for Resilient Buildings and Cities

Izuru Takewaki <sup>1</sup> \*, Masayuki Kohiyama<sup>2</sup> , Tomaso Trombetti <sup>3</sup> , Solomon Tesfamariam<sup>4</sup> and Xinzheng Lu<sup>5</sup>

*<sup>1</sup> Department of Architecture and Architectural Engineering, Kyoto University, Kyoto, Japan, <sup>2</sup> Department of System Design Engineering, Keio University, Yokohama, Japan, <sup>3</sup> Department of Civil, Chemical, Environmental and Materials Engineering, University of Bologna, Bologna, Italy, <sup>4</sup> Faculty of Applied Science, University of British Columbia, Okanagan, BC, Canada, <sup>5</sup> Department of Civil Engineering, Tsinghua University, Beijing, China*

Keywords: resilience, robustness, earthquake risk, structural control, monitoring

#### **Editorial on the Research Topic**

#### **Innovative Methodologies for Resilient Buildings and Cities**

Resilient buildings and cities are in the center of common interests in modern academic communities for science and engineering related to built environment. Resilience of buildings and cities against multidisciplinary risks, e.g., earthquakes, strong winds, floods, etc., is strongly related to the sustainability of buildings and cities in which reduction of damage during a disaster and fast recovery from the damage are key issues. The reduction of damage is related to the level of resistance of buildings and the time of recovery is affected by the amount of supply of damaged members, assurance of restoration work, etc. Robustness, redundancy, resourcefulness, and rapidity are four key factors for supporting the full realization of design and construction of resilient buildings and cities. This Research Topic gathers cutting-edge and innovative research from various aspects, e.g., robustness of buildings and cities against earthquake risk, structural control and base-isolation for controlling damage risks, quantification of resilience measures, structural health monitoring, innovative structural engineering techniques for higher safety of buildings, resilience actions and tools at the urban scale, etc. After a detailed review, the following eight papers have been published in this Research Topic, i.e., collapse risk assessment of reinforced-concrete buildings, robustness evaluation of building structures considering the whole generation process of earthquake ground motions, multi-hazard prevention and mitigation in building structures, robustness evaluation of building structures with long natural period, optimal viscous damper placement for elastic-plastic structures subjected to the critical double impulse, proposal of configurations in concentrically steel braced frame, deformability of pipes in buildings as non-structural members, assessment of seismic vulnerability of buildings at regional scale.

In the first paper, Noh and Tesfamariam present the risk assessment for collapse of reinforced concrete moment resisting frame (RC MRF) buildings located in Vancouver, Canada which were designed in compliance with the building code. They investigate three- and six-story RC MRF building systems that include unreinforced masonry infill walls by comparing with those without such infill walls. These building systems represent low- to mid-rise structures. They conduct incremental dynamic analysis and develop the seismic fragility curves. Then mean annual frequency of collapse is evaluated via a sophisticated combination of fragility curves and hazard curves. It is shown that the bare RC buildings without infill wall are not sufficiently resistant to collapse to earthquakes in the case that the number of stories increases. Furthermore, it is demonstrated that the URM infill walls are apt to influence the behavior of the framed structure to collapse. It may be important to investigate not only the collapse risk of building structures but also the detailed behavior to collapse when developing earthquake-resisting systems for upgrading the earthquake resilience.

#### Edited by:

*Andrea Belleri, University of Bergamo, Italy*

Reviewed by: *Andre R. Barbosa, Oregon State University, United States*

> \*Correspondence: *Izuru Takewaki takewaki@archi.kyoto-u.ac.jp*

#### Specialty section:

*This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment*

Received: *06 May 2019* Accepted: *28 June 2019* Published: *12 July 2019*

#### Citation:

*Takewaki I, Kohiyama M, Trombetti T, Tesfamariam S and Lu X (2019) Editorial: Innovative Methodologies for Resilient Buildings and Cities. Front. Built Environ. 5:94. doi: 10.3389/fbuil.2019.00094*

**4**

In the second paper, Makita et al. focus on the robustness evaluation under the detailed consideration of the process of theoretical ground motion generation. In the process of theoretical ground motion generation, three stages are taken into account, i.e., (i) the process of fault rupture, (ii) the characteristics of wave propagation from the fault to the earthquake bedrock, (iii) the site amplification above the earthquake bedrock. They investigate the uncertainty in the fault rupture slip (slip distribution and rupture front). The process of wave propagation from the fault to the earthquake bedrock is expressed by the stochastic Green's function method in which the Fourier amplitude of the ground motion at the earthquake bedrock from a fault element is represented by Boore's model and the phase angle is modeled by the phase difference method. The robustness for the uncertain fault rupture slip and the uncertain fault rupture front is evaluated seen to play a key role for resilient building design. Since the critical ground motion produces the most detrimental building response among possible scenarios, it is expected that the proposed method can be a reliable tool for use in resilient building design.

In the third paper, Lin et al. present a research on multihazard prevention and mitigation in building structures which is of much interest in the civil engineering field. They propose an analytical model that facilitates to calculate the structural resistance of a type of multi-hazard resilient prefabricated concrete (MHRPC) frame which is subjected to simulated earthquake ground motions and subject to scenarios of column removal. The MHRPC frame consists of unbonded posttensioning (PT) tendons, energy-dissipating steel angles, and large rotational shear plates in addition to prefabricated RC beams and columns. Through a set of laboratory experiments, it is shown that the MHRPC frame exhibits low damage features and possesses self-centering properties of under seismic loading. On the other hand, when subjected to scenarios of column removal, the MHRPC frame is proven to possess high resistance against progressive collapse. It is concluded that this paper provides useful data for the design of MHRPC frame structures which are primarily designed for earthquake and progressive collapse. It seems important to take into account the concept of multi-hazard prevention and mitigation to enable the development of structural design methods for more robust and resilient reinforced concrete building structures.

In the fourth paper, Makita et al. focus on the use of the threedimensional finite difference method (3DFDM) for generation of simulated ground motions with rather long periods. The 3DFDM is taken advantage for capturing the critical ground motion which should be considered in the design stage of structures with a rather long natural period. Since it is understood that the 3DFDM is usually time-consuming, its use in a simple sensitivity algorithm is not preferred. This is because an independent response sensitivity is to be calculated for many design parameters. To remedy this problem, the authors introduce the method of bi-cubic spline interpolation for uncertain parameter distributions (seismic moment distribution). Then the response surface method is used effectively for predicting the response surface approximately from some sampling points. The authors treat the fault rupture slip distribution described in terms of seismic moments as an uncertainty parameter. It is shown that the most critical ground motion for structures with a rather long natural period can be captured by the proposed method in an effective manner. It seems desirable to take into account the whole system of ground motion generation with a certain level of accuracy in order to develop a resilient structural design method of buildings for unpredictable earthquake occurrence.

In the fifth paper, Akehashi and Takewaki propose a new method for optimal viscous damper placement for elastic-plastic multi-degree-of-freedom (MDOF) structures subjected to the critical double impulse as a representative of near-fault ground motions. The quantities in terms of the maximum interstory drift along the building height and the sum of the maximum interstory drifts in all stories are selected as the objective function and constraint. The corresponding optimization algorithm is investigated in which time-history response analyses and sensitivity analysis are effectively used. A novel concept of double impulse pushover (DIP) is proposed for determining the input velocity level of the critical double impulse. It is demonstrated that the combination of two algorithms, one for effective reduction of the overflowed maximum interstory drift via the concentrated allocation of dampers and the other for effective allocation of dampers via the use of a stable objective function, is effective for finding a stable optimal damper placement. It seems that the optimal damper placement, i.e., effective and efficient use of damper materials, leads to the design of controlled building structures with increased robustness, redundancy, and resilience.

In the sixth paper, Yang et al. investigate the concentrically braced frame (CBF) which is a prevalent earthquake resistant system. While concentric braces have large stiffness, the check of buckling is important in its use. Several configurations of these bracing systems are proposed for different building codes available worldwide. These codes require the satisfaction of some details in the structural members and connections. However, it is true that there is no guidance in the selection of the best bracing configuration for design. The authors systematically examine the effect of the bracing configuration on the seismic response of a five-story prototype office building located in Vancouver, Canada. The authors design five different bracing configurations and investigate the candidates which follow the National Building Code of Canada and CSA S16 standard. They conduct the structural analysis in detail and systematically compare structural responses, initial costs, and life-cycle costs of the prototype building with five different bracing configurations. They derive the results that the different bracing configurations lead to different results in sizing the structural members. It seems that this result affects the initial material usage and the overall lifecycle cost of the building. It appears that the appropriate selection of bracing configuration leads to a reliable structural design of buildings with larger resilience.

In the seventh paper, Wang et al. study the pipes of a diameter of 150 mm (called DN150) which are often connected by grooved fit joints and employed as the stem pipelines. They enable to convey water vertically to different building stories and distribute them horizontally to different rooms. They are usually called the non-structural members which do not contribute to the structural safety for disturbances, but may influence significantly the functionality of buildings. Large deformability is needed in the grooved fit joints to allow the deformation between adjacent stories during an earthquake. The grooved fit joint can be improved by introducing a wider groove which enables the achievement of such large deformability. However, the authors claim that its seismic performance has never been thoroughly studied yet. This study conducted quasi-static tests on 12 DN150 grooved fit joints, including four elbow joints and eight DN150- DN80 Tee joints. The mechanical behavior, rotational capacity and failure mode were examined and discussed. Finally, seismic fragility analysis of DN150 stem pipeline system in a 10-story building is conducted. It is demonstrated that the improved joints are sufficiently resistant for the maximum plausible earthquake ground motion and the leakage can be prevented in a reliable manner. It seems that the non-structural elements are important to assure the resilience of buildings in view of the business continuity.

In the eighth paper, Diana et al. explain that assessing seismic vulnerability at large scales requires accurate attribution of individual buildings to more general typological classes that are representative of the seismic behavior of the buildings sharing the same attributes. One-by-one evaluation of all buildings is a process of requiring long time and a lot of cost. It is shown that detailed individual evaluations are only suitable for important buildings, such as hospitals and other buildings, which play a central and emergent role in the post-earthquake phase. For other ordinary buildings, more simplified approaches seem sufficient. The central issues for reliable seismic assessment at large scales are to define a taxonomy that contains the most widespread typological classes as well as performing the attribution of the appropriate class to each building. The authors explain that a fast and accurate survey process is needed to attribute a correct class to each building which composes the urban system. The authors use the association-rule learning (ARL) to find links between building attributes and typological classes. They evaluate the accuracy of wide spreading these links learned on <250 buildings of a specific district in terms of class attribution and seismic vulnerability prediction. They make it clear that the time required to provide seismic vulnerability scenarios at city scale is significantly reduced, while accuracy is reduced by <5%. The authors enable this by considering only three attributes available on public databases (i.e., period of construction, number of floors, and shape of the roof). It seems important to assess the seismic vulnerability of buildings at regional scales by conducting the definition of the accurate attribution of individual buildings to more general typological classes that are representative of the seismic behavior of the buildings with the same attributes. This enables a reliable achievement of the construction of resilient buildings with less vulnerability.

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Takewaki, Kohiyama, Trombetti, Tesfamariam and Lu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Seismic Collapse Risk Assessment of Code-Conforming RC Moment Resisting Frame Buildings Designed With 2014 Canadian Standard Association Standard A23.3

#### Nurbaiah Mohammad Noh and Solomon Tesfamariam\*

School of Engineering, The University of British Columbia, Okanagan Campus, Kelowna, BC, Canada

#### Edited by:

Luigi Di Sarno, University of Sannio, Italy

#### Reviewed by:

Emanuele Brunesi, Fondazione Eucentre, Italy Fabrizio Mollaioli, Università degli Studi di Roma La Sapienza, Italy

\*Correspondence: Solomon Tesfamariam solomon.tesfamariam@ubc.ca

#### Specialty section:

This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment

Received: 10 August 2018 Accepted: 21 September 2018 Published: 18 October 2018

#### Citation:

Noh NM and Tesfamariam S (2018) Seismic Collapse Risk Assessment of Code-Conforming RC Moment Resisting Frame Buildings Designed With 2014 Canadian Standard Association Standard A23.3. Front. Built Environ. 4:53. doi: 10.3389/fbuil.2018.00053 This paper presents collapse risk assessment of code-conforming reinforced concrete moment resisting frame (RC MRF) buildings located in Vancouver, Canada. This assessment investigates the three- and six-story of regular RC MRF building systems, with and without unreinforced masonry infill wall, representing low- to mid-rise structures. These buildings are designed according to the current National Building Code of Canada and detailed based on the 2014 Canadian Standards Association A23.3 standard provision for high seismic regions. Two different ductility classes of seismic building design, namely ductile and moderately ductile, are considered to identify the capability, equality, and/or difference of the seismic performance of these designed buildings. Nonlinear dynamic analysis is applied in the performance-based seismic assessment procedures to assess the collapse response of structural for the set of 50 pair ground motion records. Next, the seismic fragility curves are developed through incremental dynamic analysis. Finally, mean annual frequency of collapse is calculated through combination of fragility curve and hazard curve. The results indicate that the bare RC buildings are vulnerable to earthquake-induced collapse when the number of the story increased. The presence of the URM infill walls significantly influence the collapse behavior of the frame structure. Compared to moderately ductile MRF buildings, ductile MRF buildings show a better collapse performance, are strongly influenced by the capacity of the building system.

Keywords: moment resisting frame, reinforced concrete buildings, collapse risk, National Building Code of Canada, code-conforming building, masonry infill walls, fragility curve, nonlinear analysis

# INTRODUCTION

Collapse of buildings is identified as main contribution to injuries, casualties, and economic loss from past earthquakes (Porter, 2016). Current seismic design philosophy of Canadian and other countries, under rare earthquake events, focuses on buildings collapse prevention (DeVall, 2003; Haselton et al., 2011; Iervolino et al., 2017). In the current National Building Code of Canada (NBCC), reinforced concrete (RC) frame structures can be designed for two different ductility levels: ductile and moderately ductile (Mitchell et al., 2003; NBCC, 2010). The ductile and moderately ductile designs should have similar seismic performance under the site-specific design earthquake. However, the presence of unreinforced masonry infill (URM) infill walls, which is normally not considered in analysis and design process, may affect the performance and risk of collapse of infilled-RC frames in a severe earthquake to be greater or lesser than bare RC frames depending on the characteristics of URM infill walls (Mehrabi et al., 1996; Asteris et al., 2016; Sattar and Liel, 2016b). Thus, designing the buildings, by following a prescriptive design code, the collapse safety is uncertain (e.g., Burton and Deierlein, 2013; Jeon et al., 2015). In this paper, collapse risk of Canadian buildings is quantified, with seismicity of Vancouver, British Columbia, for the two ductility classes (ductile and moderately ductile) and with and without consideration of URM infill wall.

Various studies are reported on code-conforming RC moment-resisting frames (MRF). Haselton et al. (2011) implemented performance-based earthquake engineering (PBEE) approach to assess the collapse risk of ductile RC-MRF buildings designed according to ASCE 7-02 (ASCE 2002, 2005) and ACI 318-05 (2002). All models were developed using lumped-plasticity element to capture strength and stiffness degradation. Nonlinear dynamic analysis through incremental dynamic analysis (IDA) method was conducted to quantify collapse safety of 30 RC MRF varying from 1- to 20-stories. The collapse risk accounted for the effects of both record-torecord (RTR) variability and modeling uncertainty. The seismic assessment results show that the collapse probability of the buildings for the 2% in a 50-years hazard level, P(C|Sa2/50), ranges from 3 to 20%. The mean annual frequency of collapse (λc) values range in 0.7 × 10−<sup>4</sup> to 7.0 × 10−<sup>4</sup> (i.e., probability of collapse in 50 years P<sup>c</sup> = 0.4–3.4%). Koopaee et al. (2017) investigated the effect of ground motion selection method on seismic collapse fragility of 10-story RC MRF, representing mid-rise building, designed to New Zealand standard (NZS 1170.5, 2004; NZS 3101, 2006). In this study, fiber-based model was used to account the loss in vertical load carrying capacity of columns. Nonlinear dynamic analysis was carried out using multiple stripe analysis, using ground motion records selected at various hazard levels. The results indicate that the case study building has a 23% of P(C|Sa2/50), which is higher than value reported in Haselton et al. (2011), particularly when comparing with average value estimated from 8- to 12-story RC frame buildings [i.e., P(C|Sa2/50) = 17.5%]. This study found that considering loss of vertical load carrying capacity in structural model contributed to the higher value of collapse rate. Iervolino et al. (2017) quantified the seismic risk of code-conforming RC MRF buildings designed with the Italian NTC 2008 (Italian building code, 2008) design code. All the structures are designed for five sites with increasing seismic hazard. The results show that the RC MRF buildings considered in this study have λ<sup>c</sup> between 1.0 × 10−<sup>5</sup> and 8.47 × 10−<sup>5</sup> and on average λ<sup>c</sup> = 2.4 × 10−<sup>5</sup> (i.e., P<sup>c</sup> = 0.12%). The studies highlighted that as the site hazard level increases, the probability of collapse increases for all building systems.

Several studies focused on the seismic performance assessment of the RC frame with URM infill wall, which includes non-collapse or collapse, compared with bare RC frame buildings. Past earthquakes have indicated that most of the RC frame with URM infill wall buildings experiencing poor seismic performance (e.g., Zhao et al., 2009; Kam et al., 2010). Presence of URM infill wall can increase strength, stiffness and energy dissipation of the RC frame that lead to better seismic performance of buildings, but may introduced brittle or shear failure in the column due to interaction between infill and surrounding frame (e.g., Klingner and Bertero, 1978; Mehrabi et al., 1996; Stavridis et al., 2012; Stylianidis, 2012). The column shear failure is caused by the force from compressive diagonal strut in the URM infill wall. The columns shear failure can lead to a loss of the axial gravity load carrying capacity of the columns and typically occurs prior to flexural yielding of the columns at low drift levels (Burton and Deierlein, 2013). Neglecting columns shear failure due to the interaction between infill and surrounding frame in the seismic collapse assessment may lead to increase the uncertainty in the system response (Burton and Deierlein, 2013; Jeon et al., 2015; Sattar and Liel, 2016a). Therefore, a numerical modeling of infilled-RC frame that capable to simulate interaction between infill and frame and column shear failure is required to explicitly assess the collapse risk of RC frame with URM infill walls.

In this paper, the collapse risk assessment of code-conforming RC MRF buildings designed according to 2014 Canadian building code is carried out. The collapse assessments are conducted for ductile and moderately ductile RC MRF buildings. **Figure 1** shows four set of building typologies considered in this study, which focusing on two different ductility class of RC MRF buildings with and without URM infill walls. The first step of this study is characterizing a set of typical structure by the archetypical buildings concepts. The selection of the archetypical structures is based on the FEMA P695 guidelines (FEMA-P695, 2009) to quantify the building's response and performance factor in the NBCC 2010 seismic design code. The details of the archetypical RC MRFs buildings with and without the masonry infill are outlined in sections Selection of Buildings: Characteristics and Configuration and Seismic Design Procedure. Next, the nonlinear structural analysis model is developed based on the selected building designs by considering the typical and critical failure modes of the structures as described in section Nonlinear Simulation Models: RC Frames and URM Infill Walls. Then, the ground motion is selected and scaled corresponded to a target seismic hazard level as discussed in section Ground Motion Selection. Next, the nonlinear model is used to simulate the response of the building until structural collapse through incremental dynamic analysis (IDA) [section Incremental Dynamic Analysis (IDA)]. Finally, as described in section Collapse Fragility Assessment, a collapse fragility curve is developed using IDA results. For this study, the seismic collapse risk assessment involves; (1) the evaluation of the median collapse capacity [g]; (2) the collapse margin ratio (CMR), (3) the mean

annual frequency of collapse, λ<sup>c</sup> , and (4) probability of collapse in 50 years.

the equivalent static force procedure, as (Mitchell et al., 2003):

# CODE-CONFORMING RC FRAME STRUCTURES AND GROUND MOTION SELECTION

## Selection of Buildings: Characteristics and Configuration

For seismicity of Vancouver, BC, two sets of office building occupancy code-conforming RC MRF buildings (ductile and moderately ductile detailing) and building heights (3- and 6 stories) are designed. These buildings are founded on soft rock, which is the reference ground condition (class C) in the NBCC 2010. The Canadian code specifies that contribution URM infill wall is not considered in the design (DeVall, 2003). However, in this study, collapse risk assessment is carried for the RC building with and without consideration URM infill walls (**Figure 2B**). The building is regular in plan and elevation as shown in **Figure 2**. The RC office building considered has 7–6 m bays in N-S direction and 3 bays in the E-W direction consisting of 2– 9 m and a bay (**Figure 2**). The story height for each level is 3.65 m. The following section explains the seismic design procedure and detailing according to the current NBCC and CSA A23.3-14 (CSA, 2014) standard provision, respectively.

## Seismic Design Procedure

Each MRF RC building is designed based on the 2010 NBCC (NBCC, 2010) seismic design provision. The minimum design base shear force, of each studied building, was computed using

$$\mathbf{V} = \frac{\mathbf{S}(\mathbf{T\_a})\mathbf{M\_v}\mathbf{I\_E}\mathbf{W}}{\mathbf{R\_d}\mathbf{R\_o}}\tag{1}$$

where, S (Ta) is design spectral response acceleration expressed as a ratio of gravitational acceleration, g, at the design period Ta, M<sup>v</sup> is the factor to account for higher mode effects, I<sup>E</sup> is the earthquake importance factor, W is the building seismic weight including 25% of snow load, and R<sup>d</sup> and R<sup>o</sup> are the ductilityrelated and overstrength-related force modification factors. The empirical fundamental lateral period, T<sup>a</sup> is given by 0.075h3/<sup>4</sup> n , where h<sup>n</sup> is height of building. However, since the dynamic fundamental periods usually greater than T<sup>a</sup> value (more than 150%), therefore, the fundamental lateral period value used to select the design spectral response acceleration is taken as 1.5Ta.

The buildings are designed with importance factor, I<sup>E</sup> = 1.0, on very dense soil and soft rock (Soil class C) and assumed fixed at ground level. Two different levels of ductility factors, which are ductile (R<sup>d</sup> = 4.0 and R<sup>0</sup> = 1.7) and moderately ductile (R<sup>d</sup> = 2.5 and R<sup>0</sup> = 1.4), are also considered (Mitchell et al., 2003). Since there is no eccentricity in the buildings, the effects of accidental torsional eccentricity is calculated by applying the lateral forces F<sup>x</sup> at an accidental eccentricity of ±0.1Dnx =4.245 m, where Dnx is a plan dimension of the building in the computed eccentricity direction that is perpendicular to the direction of seismic loading. The same floor plan is used at all levels as illustrated in **Figure 2**. Typical super imposed load values of 1.0, 0.5, 0.5 kN/m<sup>2</sup> are used as partition, mechanical services, and roofing material, respectively. The office floor live load of 2.4 kN/m<sup>2</sup> and snow load of

2.2 kN/m<sup>2</sup> are considered. Concrete compressive strength, f<sup>c</sup> ′ and concrete modulus of elasticity, E<sup>c</sup> was taken as 30 and 24,700 MPa, respectively. High yielding strength deformed rebar with yield strength, f<sup>y</sup> = 400 MPa is used. The E-W moment resisting frame was chosen for the current study, as shown in **Figure 2B**. ETABS 2015, which has an option for current Canadian code CSA A23.14, is utilized to analyze and design three-dimensional structure system. The modal response spectrum analysis procedure including second order P-Delta effects is considered in the design.

**Table 1** provides a summary of the design characteristics for beams and columns. The interior and exterior column is 500 × 500 and 450 × 450 mm, respectively. The columns are assumed fixed at the base with ignoring soil-structure interaction. The beams of both the N-S and E-W frames are 400 mm wide × 600 mm deep for first three stories and 400 × 550 mm for other stories. The two-way slab floor system consists of a 210 mm thick slab. For the infilled-RC frame buildings, the URM infill walls are considered as a non-structural element and neglected in the structural design process. Thickness and compressive strength of infill panels are 150 mm and 6.9 MPa, respectively.

#### Nonlinear Simulation Models: RC Frames and URM Infill Walls

A two-dimensional (2-D) models of code-conforming RC frame buildings with and without URM infill walls (**Figure 3A**) are developed using the Open System for Earthquake Engineering Simulation (OpenSees) software package (OpenSees, 2009). The 2-D model employs a leaning column to account for P-1 effects. For the RC members, the model used lumped- plasticity element models and rigid beam-column joints zones. The lumpedplasticity model can be introduced using an elastic beam-column element with two zero-length elements at both element ends, as shown in **Figure 3B**. The zero length-elements are related to a rotational hinge model with a hysteretic response to represent the flexural behavior for the elements. The rotational hinge behavior is attached by the multiple uniaxial materials to display the moment-rotation or force-deformation relationship. This simulation approach is used for the collapse prediction to capture deterioration of the steel reinforcing bars due to the rebar buckling and low cyclic-fatigue, as well as to record the strength and stiffness deterioration in assessing the global collapse. In this study, Rayleigh damping model is used, which is a damping ratio of 5% is assigned to the first and third modes of the structure. For nonlinear dynamic analysis, the Newton algorithm is used to solve the system equations.

To describe the flexural behavior of the RC beam-column element, the most common OpenSees implementation of the peak oriented hysteretic model developed by Ibarra et al. (2005) (modIMKmodel), is adapted. This model captures four modes of cyclic deterioration, which includes; (1) the basic strength deterioration, (2) the post-capping strength deterioration, (3)



<sup>a</sup>Beam size for level 1–3, <sup>b</sup>Beam size for level -6, <sup>c</sup> Internal column size, <sup>d</sup>External column size.

(E) hysteretic curve (cyclic behavior).

the unloading stiffness deterioration, and (4) the accelerated reloading stiffness deterioration. The hysteretic model requires six key parameters: elastic stiffness (Ke), effective yield moment (My), strain hardening ratio (Mc/My), pre-capping rotation (θcap), post-capping rotation (θpc), and cyclic deterioration parameter (λ). The moment-rotation law (backbone) according to Ibarra et al. (2005) is illustrated in **Figure 3C**. All the parameters are obtained using the predictive equations and tools developed by Haselton et al. (2007). The M<sup>y</sup> was calculated for each element from the sectional analysis.

Analytical method to model infilled-frame buildings with consideration of shear failure are limited in the literature. Crisafulli and Carr (2007) developed a new macro model for the seismic response of infilled RC frames, which is can be modeled by four-node element. This model developed based on a multiple-strut formulation that includes two parallel off-diagonal struts and a special shear spring in each direction to account the diagonal tension failure and shear mechanism in the infill walls. Although the model capable to represent different type of failure modes in shear for URM infill wall, but the bending moment and shear forces of the surrounding frame not able to predict appropriately due to the simplicity of the model. Jeon et al. (2015) adopted a three-strut model developed by Chrysostomou et al. (2002) to assess seismic fragility of lightly RC frames with URM infills. Fiber-type displacement based model was used to capture the flexural response of the surrounding frame. Two off-diagonal struts and one central strut at each direction were employed to simulate the column shear failure. A contact length between the column and infill panel is computed by using Smith (1967) approach. Two types of column shear failure were considered: shear failure and flexure-shear failure. The column shear failure was simulated using two zero-length springs at the face of the beam and the flexure-shear mode the column shear strength model adopted from Sezen and Moehle (2004). For the flexure-shear failure has been modeled using a zero-length spring at center of the column incorporating with limit state model proposed by Elwood (2004). Burton and Deierlein (2013) applied a pair of central and off-diagonal compression-only strut at both direction with zero contact length to assess seismic performance of the non-ductile RC frame with URM infill walls. Lumpedplasticity model that can be introduced using an elastic beamcolumn with two zero-length elements at both ends is used to idealize beams and columns response. The flexure hinge and shear spring in series was introduced at top and bottom column to capture the column shear failure, and has been modeled based on the flexure-shear failure type that occurred after infill failure. Sattar and Liel (2016b) assessed the collapse risk of masonryinfilled RC frame buildings using the strut modeling enhanced by finite element (FE) analysis. Two struts in each direction were employed to simulate the column shear failure, which the contact length were determined through FE model response. This study also implemented lumped-plasticity model with spring in series at the top of column only. Burton and Deierlein (2013) and Sattar and Liel (2016a) also used similar approach developed by Sezen and Moehle (2004) to calculate the column shear strength. Although all studies applied similar approach to capture the response of strut and column shear failure, but the infilled-RC frame model developed by Burton and Deierlein (2013) is less sophisticated and able to simulate most of the seismic collapse mode of buildings. Therefore, in this study, the model and tools developed by Burton and Deierlein (2013) for the infilled-frame structures is adopted (**Figure 3B**).

A pair of off-diagonal strut, as presented in the **Figure 3B**, is used to capture shear failure due to the interaction between the column and the infill walls. As illustrated in **Figure 4**, zero-length shear springs model with a rigid softening material model is attached in series with the flexural hinges at the face of the columns. The column shear strength is calculated by using a model developed by Sezen and Moehle (2004). Whereas, the column shear deformation are



obtained from the modeling parameter and acceptance criteria offered in Tables 6–8 of ASCE/SEI 41 (2007), which a proposed supplement to ASCE/SEI 41 developed by Elwood et al. (2007). Furthermore, this study assumed that collapse generated by axial column failure occur at ultimate shear deformation.

For the URM infill wall, an equivalent diagonal struts model is implemented to simulate the behavior of the masonry infill in terms of strength and stiffness. The performance of the masonry infill wall subjected to the in-plane effect is simulated by using two opposing pairs of the diagonal compression-only strut in each direction (**Figure 3B**). The monotonic backbone curve, shown in **Figure 3D**, proposed by Decanini et al. (2014) is applied to define the overall in-plane behavior of URM infill walls. The strength of the infill model is calculated using strength model proposed by Priestley and Calvi (1991) and Paulay and Priestley (1992). The distribution amount of force and stiffness of the strut are adopted from Chrysostomou (1991) who has investigated the force and stiffness distribution of central and off-diagonal strut of infilled-frame buildings using the principle of virtual displacements. The author founds that the maximum force distributed to the off-diagonal strut is ∼25% of the total strut force. Thus, in this study, the stiffness proportion assigned to the central and off-diagonal strut are 75 and 25% of the total strut stiffness, respectively (**Figure 4**). The available elements and materials in the OpenSees are used to define the model. The selection of the most appropriate hysteretic model used to calculate and capture the strength, stiffness, and behavior of the URM infill walls is presented in Noh et al. (2017). **Figure 3E** shows the cyclic behavior of single-story infilled-frame by executing monotonic curves of RC frame member and infill wall, as illustrate in **Figures 3C,D**.

**Table 2** presented the first three mode period of vibration obtaining through modal analysis conducted in OpenSees model. The dynamic fundamental period that calculated for the codeconforming RC frame models is based on the effective cracked stiffness of the structural members (**Table 2**). As shown in **Table 2**, the spectral acceleration at 1.5 s is selected and the time range used to select ground motion is between 0.2 and 2.0 for the 3-story building without infill. For the 3-story building with infill, the spectral acceleration is selected at 0.4 s. While for 6-story

building without infill the spectral acceleration is selected at 2.0 and 2.5 s for ductile and moderately ductile, respectively, ranges from 0.2 to 3.0 s. The spectral acceleration at 0.6 s is selected for a 6-story building with infill, with ranges from 0.2 to 1.0 s. As reported in **Table 2**, the dynamic fundamental period obtained for same building height is not constant as suggested in the code provision. It is also observed that the fundamental period of the infilled-RC frame decreases due to the presence of URM infill wall. According to Banik and Halder (2016) the fundamental period of the RF frame buildings not only depends on the building height, but the fundamental period also influenced by effective stiffness properties of structural members, bay width, number of bays and support conditions.

#### Ground Motion Selection

The RC buildings located in Vancouver, BC, Canada was considered due to its location and population and significant seismic activities (Onur et al., 2005; Atkinson and Goda, 2011). The selected building site consisting of very dense soil and soft rock site condition (Site class C) with the average shear-wave velocity in the upper 30 m (Vs30) between 360 and 760 m/s. Sa(T1) of the buildings with 5% damping.

In order to perform the nonlinear dynamic analysis, a set of 50 ground motion records is selected based on the seismic hazard of interest. In this study, ground motion records were selected using the multiple conditions mean spectrum (CMS) method developed by Goda and Atkinson (2011). The selected ground motion records that follow same procedure reported in Tesfamariam et al. (2015) and Tesfamariam and Goda (2015a,b), that considered three dominant sources of earthquake: shallow crustal, deep in-slab, and off-shore megathrust interface earthquakes type from the Cascadia subduction zone. **Figure 5** shows an example of the ground motion records selected for 3S\_D&MD\_NI buildings with T<sup>1</sup> = 1.5 s.

#### INCREMENTAL DYNAMIC ANALYSIS (IDA)

A series of nonlinear dynamic analyses is conducted through IDA to evaluate the seismic response and collapse fragility (Vamvatsikos and Cornell, 2002). The IDA was carried out based on multiple records of ground motions by increasing scale factor of each ground motion amplitude until cause sideways collapse. In this study, the IDA is implemented on eight RC MRF building models using 50 pairs of ground motion records selected based on the 2% PE in 50 years. Collapse is defined as the dynamic instability point when the maximum inter-story drift ratio of the buildings exceeds 0.10 (Haselton et al., 2007) and IDA curve transforms to horizontal line for each ground motion. The IDA curves provide a relationship between intensity measure level and maximum inter-story drift ratio.

In this assessment, collapse mechanisms due to column shear failure are not directly simulated in the analysis model. However, the columns are expected to experience critical shear failure when interacting with URM infill walls. Therefore, the shear critical failure modes, so-called non-simulated collapse mode, is evaluated by post-processing IDA result. In this process, forcebased limit state check is applied to predict the median collapse drift ratio (CDR) at which critical shear failure will occur. The [ CDR is defined using Equation (2) developed by [ Aslani and Miranda (2005),

$$\widehat{\text{CDR}}\_{\text{shear}} = \frac{1}{0.26 \left( \frac{\text{p}}{\text{A}\_{\text{g}} \text{f}\_{\text{c}} / \rho\_{\text{sh}}} + 25.4 \right)} \ge \frac{1}{100} \tag{2}$$

where P is the axial load on the column, A<sup>g</sup> is the gross cross-sectional area of the column, and ρsh is the transverse reinforcement ratio. The median CDR is calculated based on the properties of the column. The capacity of the deformation increases with decreasing axial load or more shear reinforcement.

**Figures 6**, **7** show IDA results for the 3- and 6-story RC MRF buildings. The circle marks on each IDA curve, in **Figures 6C,D**, **7C,D**, shows the shear failure occurs predicted by Equation (2). For columns with ductile shear design and detailing, the predicted shear failure occurs at a CDR between 0.0295 and [ 0.0339 rad, while for moderately ductile the predicted shear failure occurs at a CDR between 0.0206 and 0.0275 rad. From the [ IDA results, the collapse fragility relationships, as described in the following section, are then developed to compute the median collapse capacity.

#### COLLAPSE FRAGILITY ASSESSMENT

Seismic collapse fragility curve shows the probabilistic relationship between the frequency of failure (collapse) of buildings and peak ground motion acceleration in an earthquake (Porter, 2017). In this study, a damage measure, which represents actual response of the building under seismic excitation, is based on maximum interstory drift (i.e., maximum interstory drift ratio is 0.1 for collapse). The collapse fragility curve that constructed based on the IDA results can be represented by a lognormal cumulative distribution function (Baker, 2015), as shown in Equation (3).

$$\mathbf{P}\left(\mathbf{C}|\mathbf{IM}=\mathbf{x}\right) = \Phi\left(\frac{\ln\left(\frac{\mathbf{x}}{\mu}\right)}{\beta}\right) \tag{3}$$

where P (C|IM = x) is the probability that a ground motion with IM = x will cause the structure to collapse, 8 () is the standard normal cumulative distribution function, µ is the median of the fragility function and β is the standard deviation of lnIM (sometimes referred to as the dispersion of IM).

**Figures 8A,B** show the collapse fragility curves for 3- and 6 story RC MRF buildings derived from IDA and Equation (2). In this study, the collapse fragility assessment is also accounted for both record-to-record (RTR) variability and structural modeling (M) uncertainties for a comparison purposes. A modeling uncertainty of βln, M = 0.50 is assumed as suggested by Haselton and Deierlein (2007). The total uncertainty, βln, Total that account for both RTR and M is computed by using the square root of the sum of the squares (Equation 4),

$$
\beta\_{\text{ln,Total}} = \sqrt{\beta\_{\text{ln,RTR}}^2 + \beta\_{\text{ln,M}}^2} \tag{4}
$$

The dotted line is used to represent the RTR+M fragility curve. **Figure 8** also illustrates the collapse fragility curve for RC MRF with and without consideration of shear failure, as shown by the solid and dashed curve, respectively.

Based on these fragility curves, the collapse performance parameters, in terms of the median collapse capacity, can be obtained at the ground motion capacity with 50% of probability of reaching the collapse targets. The results in **Figure 8** shows that the median collapse capacity for ductile bare RC MRF buildings ranging from 0.15 to 0.53 g, which is greater than the moderately ductile buildings varying from 0.08 to 0.41 g. Presence of URM infill walls on the RC MRF frame considered in this study increase the median collapse capacity by 66–90% and 74 and

94% for ductile and moderately ductile building, respectively, as compared to the bare RC MRF buildings. As shown in **Figure 8**, the collapse fragility capacity with consideration of shear failure significantly reduced by 8–30% compared to the sidesway collapse mechanism. Overall, the median collapse capacity decreases with increasing building height.

Median collapse capacity results are then expended to the evaluation of collapse median ratio (CMR). The CMR is evaluated by dividing the median collapse capacity of the structure by the ground motion intensity, SMT. As recorded in **Table 2**, the ductile RC MRF buildings with and without URM infill walls induce higher CMR compared to moderately ductile RC MRF buildings, ranging from 0.65 to 2.25, and 0.45 to 1.92, respectively. Any CMR value that <1.0 indicates that the median collapse capacities are lower the SMT and vice-versa.

#### MEAN ANNUAL FREQUENCY OF COLLAPSE

The mean annual frequency of collapse in seismic risk assessment (Porter, 2017 and Ellingwood and Kinali, 2009) is given by:

$$
\lambda\_{\rm c} = \int\_0^\infty \text{F}\_{\rm R}(\mathbf{x}) \frac{\text{dH}(\mathbf{x})}{\text{dx}} d\mathbf{x} \tag{5}
$$

where FR(x) is the cumulative distribution function of seismic capacity of structural (collapse fragility curve), and dH(x) is the differential of the seismic hazard curve. Equation (5) can be expressed as

$$
\lambda\_{\mathbf{c}} = \sum \mathbf{P} \left[ \mathbf{C} | \mathbf{Q} = \mathbf{x} \right] \mathbf{P} [\mathbf{Q} = \mathbf{x}] \tag{6}
$$

where P[Q = x] denote the seismic hazard and P[C|Q = x] represents the seismic fragility of a structure with collapse damage.

The ground motion hazard curves were selected at Vancouver, BC, Canada for site classification C and a period according to building model. **Figure 9** shows the hazard curve for Vancouver at different periods.

The probability of collapse that will occur at least once in time, t (t = 50 years) is calculated using the following equation:

$$\mathbb{P}[\mathcal{C}\text{ in }t\text{ years}] = 1 - \exp(-\lambda\_c.t) \tag{7}$$

where λ<sup>c</sup> is the expected value of collapse derived from Equation (5).

Seismic collapse assessment results for 3-and 6-story RC MRF buildings with and without URM infill walls, in terms of mean annual frequency of collapse and probability of collapse at 50 years, are reported in **Table 3**. **Figure 10** shows the mean annual frequency of collapse of this study including results from different studies that used different code during seismic design process. Data in **Table 3** also reported for both RTR and RTR+M uncertainties.

As reported in **Table 3** and illustrated in **Figure 10**, the mean annual frequency of collapse for ductile RC MRF buildings is between 16 and 39 × 10−<sup>5</sup> with probability of collapse in 50 years ranging from 0.8 to 1.92%. For the moderately ductile RC MRF buildings, the calculated mean annual frequency of collapse is

design codes.


TABLE 3 | Result of seismic collapse assessment for ductile and moderately ductile RC MRF with and without URM infill walls.

between 13 and 165 × 10−<sup>5</sup> with probability of collapse in 50 years ranging from 1.05 to 7.94%. These results reveal that the seismic safety level is higher for the ductile RC MRF buildings than moderately ductile RC MRF buildings, due to the special detailing of reinforcement during design in increasing resistance for structural collapse.

Furthermore, this study also found that the 6-story bare RC MRF buildings have slightly higher mean annual frequency of collapse than the 3-story bare RC MRF building, which have λ<sup>c</sup> = 39 to 165 × 10−<sup>5</sup> , corresponding to 1.92 to 7.94% probability of collapse in 50 years. This study observes that the increase in building height will increase the overall mass of the building that becomes dominant over the increase in the lateral stiffness. Moreover, this is due to lower spectral acceleration based on the hazard site for longer periods during building design stage, which in turn may affect the collapse capacity of the buildings.

Comparison of the bare RC MRF and RC MRF with URM infill walls buildings illustrates inconsistent results among infilled-RC MRF buildings. As observed, the 6-story RC MRF building with URM infill walls show lower collapse rate compared to the 6-story bare RC MRF buildings, ranging from 16 to 21 × 10−<sup>5</sup> . Whereas, the collapse rate of 3-story RC MRF with URM infill walls buildings are 20–80 × 10−<sup>5</sup> illustrates higher collapse rate compared to the 3-story bare RC MRF buildings. The results of 6-story bare RC MRF buildings and 3-story RC MRF building with URM infill walls buildings are at least 2.5 and 1.2 times, respectively, to collapse compared to the 6-story RC MRF building with URM infill walls. Collapse performance of the RC MRF buildings with URM infill walls shows not uniform results may be effects of the geometric and material properties of the URM infill walls and surrounded frame. Besides, as the number of the stories increases, the inter-story drift ratios will be increased due to the increase in the building stiffness offered by the URM infill walls. Overall, this observation denotes that the RC MRF buildings with ductile design provide better collapse resistance and enhanced the safety of the buildings.

The results of mean annual frequency of collapse of ductile RC MRF buildings are further compared with studies conducted by Haselton et al. (2011) and Iervolino et al. (2017), whereby each building is designed based on different code provision that includes American code, and Italian code, respectively. For comparison purposes, the collapse assessment are evaluated by considering both record-to-record and modeling uncertainties. It can be observed that there is significant variation in the mean annual frequency of collapse for different designed buildings, which may due to the differences in design spectral, response reduction factors, site location, ground motion selection method, and structural configurations. As shown in **Figure 10**, RC frame building designed based on Canadian code have calculated mean annual rate of collapse between 1.6 and 3.9 × 10−<sup>4</sup> whereas American code (Haselton et al., 2011) have calculated mean annual rate of collapse between 1 to 4.7 × 10−<sup>4</sup> . The mean annual rate of collapse for buildings designed based on Italian code ranging from 0.13 to 0.85 × 10−<sup>4</sup> . The range value of mean annual frequency of collapse for Canadian buildings is similar to the American buildings, which can considerably perceived to be within a reasonable range as reported in Haselton et al. (2011). In comparison, the buildings designed for Italian codes has significantly lower mean annual rate of collapse than the others buildings designed for Canadian and American codes. The decrease in collapse probability for Italian buildings is due to different response reduction factors during design process, site location and response history analysis. In addition, different approach to nonlinear dynamic analysis is adopted, where Canadian and American buildings assessed using IDA method whereas Italian buildings assessed using multiple stripe analysis method.

## DISCUSSION AND CONCLUSION

This study evaluates the seismic collapse risk of two sets of code-conforming RC MRF buildings, considering ductile and moderately ductile design details, and with and without URM infill walls. The code-conforming RC MRF buildings were designed according to current Canadian seismic design code, NBCC 2010, and detailing, CSA A23.3-14. IDA method is used to perform nonlinear dynamic analysis. Collapse fragility relationships are developed from IDA results to represent the probability of building collapse as a function of ground motion intensity. The collapse mean annual frequency is determined by integrating seismic hazard and fragility function.

The results from this assessment identified that the mean annual frequency of collapse varied between 16 and 165 × 10−<sup>5</sup> correspond to the collapse probabilities in 50 years from 1 to 7.9%. The ductile RC MRF buildings present better seismic collapse performance, for both 3- and 6- story buildings when compared to moderately ductile RC MRF buildings. These results are expected since the ductile RC MRF buildings are designed based on special reinforcement detailing with higher overstrength that importance to against structural collapse. It was observed that structural configuration may lead to different result of collapse probability at 50 years.

Seismic collapse assessments are also conducted to show the seismic performance over building height. Comparative studies between 3- and 6-story RC MRF buildings without infill demonstrate that as the building height increases, the calculated mean annual frequency of collapse is increases. This study would be useful to extend by investigating the influence of soil-structure interaction (SSI) effect, which this effect also may contribute to collapse risk. The collapse response is expected to be different for RC frame with SSI effect between building height and structural systems because the flexibility of the soil foundation reduces the overall stiffness of the structures that lead to expanding in the fundamental period of the system.

This study also investigated the influence of URM infill walls on the seismic collapse risk of buildings. A comparative study of this building system shows that the RC MRF buildings with URM infill walls are capable to resist higher ground motion of intensities. However, there is an inconsistency in term of collapse rate, whereby the 6-story RC MRF buildings have higher collapse rate when compared among the buildings considered. The presence of infill walls, in terms of material and geometric characteristic, and number of building story are considered as the main factors associated with these outcomes. A more comprehensive study of URM infill wall is recommended to enhance understanding on the overall contribution of infill wall to collapse safety by considering both in-plane and out-of-plane URM infill walls effect.

This study further compares the mean annual frequency of collapse of ductile RC MRF buildings for different codes. It has been observed that the RC MRF buildings designed by Canadian code is in reasonable range of collapse risk for buildings. Further, there is significant variation in the mean annual frequency of collapse of the building designed for Italian codes.

The results presented in this paper only focused on the lowto mid-rise buildings. However, the results obtained from this study are helpful to comprehend seismic safety level of current Canadian buildings. This study can be extended by considering the SSI effect and both in-plane and out-of-plane URM infill walls effect to obtain better understanding on the overall seismic collapse performance of buildings in Canada.

# AUTHOR CONTRIBUTIONS

NMN carried out this research under the supervision of ST.

# ACKNOWLEDGMENTS

The first author wishes to thank the financial support by Ministry of Higher Education (MoHE) Malaysia and study leave provided by Universiti Teknologi MARA (UiTM), Malaysia to undertake this research work at The University of British Columbia. The authors acknowledge the financial support through the Natural Sciences an Engineering Research Council of Canada (RGPIN-2014-05013) under the Discovery Grant programs.

Frontiers in Built Environment | www.frontiersin.org

# REFERENCES


Zhao, B., Taucer, F., and Rossetto, T. (2009). Field investigation on the performance of building structures during the 12 May 2008 Wenchuan earthquake in China. Eng. Struct. 31, 1707–1723. doi: 10.1016/j.engstruct.2009.02.039

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Noh and Tesfamariam. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Critical Ground Motion for Resilient Building Design Considering Uncertainty of Fault Rupture Slip

Koki Makita, Kyoichiro Kondo and Izuru Takewaki\*

Department of Architecture and Architectural Engineering, Graduate School of Engineering, Kyoto University, Kyoto, Japan

The process of theoretical ground motion generation consists of (i) the fault rupture process, (ii) the wave propagation from the fault to the earthquake bedrock, (iii) the site amplification. The uncertainty in the site amplification was taken into account in the previous research (Makita et al., 2018). On the other hand, the uncertainty in the fault rupture slip (slip distribution and rupture front) is dealt with in the present paper. The wave propagation from the fault to the earthquake bedrock is expressed here by the stochastic Green's function method in which the Fourier amplitude of the ground motion at the earthquake bedrock from a fault element is represented by the Boore's model and the phase angle is modeled by the phase difference method. The validity of the proposed method is investigated through the comparison with the existing simulation result by other methods. By using the proposed method for ground motion generation and for optimization under uncertainty in the fault rupture slip, a methodology is presented for deriving the critical ground motion imposing the maximum response of an elastic SDOF model at the earthquake bedrock or at the free ground surface. It is shown that the critical response exhibits the SDOF response several times larger than that due to the average fault rupture slip model. Furthermore, the robustness evaluation with respect to the uncertain fault rupture slip and the uncertain fault rupture front is presented for resilient building design. Since the critical ground motion produces the most detrimental building response among possible scenarios, the proposed method can be a reliable tool for resilient building design.

Keywords: critical ground motion, worst input, stochastic Green's function method, fault rupture, wave propagation, phase difference, site amplification, resilience

# INTRODUCTION

Many peculiar earthquake ground motions have been observed in the world, e.g., Mexico (1985), Northridge (1994), Kobe (1995), Chi-chi (1999), Tohoku (2011), Kumamoto (2016). To model these ground motions from their occurrence mechanisms, several models have been proposed. The whole process of ground motion generation consists of (i) the fault rupture process, (ii) the wave propagation from fault to the earthquake bedrock, (iii) the site amplification. These models can be classified generally into the theoretical approach, the numerical analysis approach, the semi-empirical approach and the hybrid approach. In the theoretical approach and the numerical analysis approach, the wavenumber integration method and the finite difference method are the representatives and are suitable for the generation of directivity pulses and surface waves with the

Edited by:

Ioannis Anastasopoulos, ETH Zürich, Switzerland

#### Reviewed by:

Emmanouil Rovithis, Institute of Engineering Seismology and Earthquake Engineering (ITSAK), Greece Raffaele De Risi, University of Bristol, United Kingdom

> \*Correspondence: Izuru Takewaki takewaki@archi.kyoto-u.ac.jp

#### Specialty section:

This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment

Received: 13 July 2018 Accepted: 12 October 2018 Published: 07 November 2018

#### Citation:

Makita K, Kondo K and Takewaki I (2018) Critical Ground Motion for Resilient Building Design Considering Uncertainty of Fault Rupture Slip. Front. Built Environ. 4:64. doi: 10.3389/fbuil.2018.00064

predominant period longer than 1–2 s (Bouchon, 1981; Hisada and Bielak, 2003; Yoshimura et al., 2003; Nickman et al., 2013). On the other hand, the semi-empirical approach is suitable for the generation of random ground motions with the predominant period shorter than 1–2 s and can generate a large ground motion in terms of small-size ground motions using the scaling law of fault parameters. The empirical Green's function (Wennerberg, 1990) and the stochastic Green's function (Hisada, 2008) are often used in the semi-empirical approach. The hybrid approach is the method which combines the random ground motions of shorter predominant period with the waves of longer predominant period by using a matching filter.

Although most of the previous approaches of ground motion generation were aimed at generating ground motions for a fixed set of parameters, several parameters should be treated as uncertain numbers (aleatory or epistemic) to make the approach more reliable (Abrahamson et al., 1998; Lawrence Livermore National Laboratory, 2002; Morikawa et al., 2008; Cotton et al., 2013).

As for researches on the effect of uncertainty of parameters on response variability, Taniguchi and Takewaki (2015) derived the bound of earthquake input energy to building structures by considering shallow and deep ground uncertainties and soilstructure interaction. Okada et al. (2016) proposed a new interval analysis technique for a soil-pile-structure interaction model by taking into account the uncertainty in soil properties. Makita et al. (2018) considered a base-isolation, building-connection hybrid structural system (Murase et al., 2013, 2014; Kasagi et al., 2016; Fukumoto and Takewaki, 2017) and took into account the uncertainty in the site amplification. They treated the fault as a point source. On the other hand, the uncertainty in the fault rupture slip (slip distribution and rupture front) is dealt with in the present paper. The wave propagation from the fault to the earthquake bedrock is expressed by the stochastic Green's function method (Irikura, 1986; Yokoi and Irikura, 1991) in which the Fourier amplitude at the earthquake bedrock from a fault element is represented by the Boore's model (Boore, 1983) and the phase angle is modeled by the phase difference method (Yamane and Nagahashi, 2008). The validity of the proposed method is investigated through the comparison with the existing simulation result by other methods.

By using the proposed method for ground motion generation and for optimization under uncertainty in the fault rupture slip, a methodology is presented for deriving the critical ground motion causing the maximum response of an elastic SDOF model at the earthquake bedrock or at the free ground surface (Drenick, 1970; Takewaki, 2007). The uncertainty in the fault rupture slip is treated by using an interval analysis in which the slip distribution and rupture front are modeled as interval parameters, i.e., the parameters in the certain prescribed range can take any value in a non-probabilistic sense (Ben-Haim, 2006). It is shown that the critical ground motion imposes the maximum SDOF response which may be several times larger than that computed under the average fault rupture slip model. Furthermore, the robustness evaluation with respect to the uncertain fault rupture slip and the uncertain fault rupture front is presented for resilient building design. Since the critical ground motion produces the worst building response among possible scenarios, the proposed method can be a reliable tool for resilient building design.

## STOCHASTIC GREEN'S FUNCTION METHOD FOR GROUND MOTION GENERATION

In a previous research by the authors (Makita et al., 2018), a point-source model of the fault rupture was assumed and the fault rupture process could not be taken into account. In this paper, the stochastic Green's function method based on a planesource model of the fault rupture is introduced. The method will be explained in the following section.

# Ground Motion Generation Using Scaling Law

The generation of ground motions using the plane-source model of the fault rupture is conducted by dividing the fault plane into many fault elements and considering the delay of the fault element rupture initiation in the fault rupture process. The stochastic Green's function method is used for generating a small ground motion resulting from the rupture of a fault element.

Assume that the fault plane is divided into N<sup>L</sup> × N<sup>W</sup> fault elements (NL: number of divisions in the longitudinal direction, NW: number of divisions in the width direction) and the slip in one fault element is divided into N<sup>D</sup> slips. It was made clear by Irikura (1983) that NL, NW, and N<sup>D</sup> can be regarded almost equal to the cubic root of the product of the ratio of the seismic moment M<sup>0</sup> <sup>L</sup> of the whole fault to the seismic moment M<sup>0</sup> <sup>S</sup> of the fault element and the ratio of the stress drop 1σ<sup>S</sup> of the fault element to the stress drop 1σ<sup>L</sup> of the whole fault. In the stochastic Green's function method, this scaling law is usually used and is expressed by

$$\frac{L\_L}{L\_S} = \frac{W\_L}{W\_S} = \frac{D\_L}{D\_S} = \frac{\mathbf{r}\_L}{\mathbf{r}\_S} = \left(\frac{M\_{0L}}{\left(\Delta\sigma\_L/\Delta\sigma\_S\right)M\_{0S}}\right)^{1/3} \approx N\_L, N\_W, N\_D \tag{1}$$

where L, W, D, τ denote the fault length, fault width, fault slip, rise time (fault slip time), respectively, and ( )<sup>L</sup> , ( )<sup>S</sup> indicate the quantity related to the whole fault and that to the fault element.

When 1σL/1σ<sup>S</sup> = 1, the ground motion displacement Uij (t) due to one fault element is produced by N<sup>D</sup> slips uij (t) and is expressed by

$$\begin{split} U\_{\vec{l}\vec{j}}\left(t\right) &= f\left(t\right) \* u\_{\vec{l}\vec{j}}\left(t\right) \\ &= \sum\_{k=1}^{N\_D} u\_{\vec{l}\vec{j}}\left(t - \left(k - 1\right)\frac{\tau\_{\vec{l}\vec{j}}}{N\_D}\right) \end{split} \tag{2}$$

where ij indicates the ij sub-element in one fault element and τij is the rise time of the ij sub-element. Furthermore f (t) is the slip correction function specifying the initiation of slips in one fault element and is expressed by

$$f\left(t\right) = \sum\_{k=1}^{N\_D} \delta\left(t - \left(k - 1\right)\frac{\tau\_{\vec{ij}}}{N\_D}\right) \tag{3}$$

where δ(t) is the Dirac delta function. **Figure 1A** shows the conventional slip function D(t) and slip correction function f(t).

From these preparations, the ground motion displacement U (t) due to the whole fault may be expressed by

$$\begin{aligned} U(t) &= \sum\_{i=1}^{N\_W} \sum\_{j=1}^{N\_L} f\left(t - t\_{\vec{ij}}\right) \* u\_{\vec{ij}}(t) \\ &= \sum\_{i=1}^{N\_W} \sum\_{j=1}^{N\_L} \sum\_{k=1}^{N\_D} u\_{\vec{ij}}\left(t - t\_{\vec{ij}} - \left(k - 1\right)\frac{\tau\_{\vec{ij}}}{N\_D}\right) \end{aligned} \tag{4}$$

In this paper, the slip correction function f(t) revised by Irikura (1986) and Yokoi and Irikura (1991) is used and is described by

$$f\left(t\right) = \delta\left(t\right) + \frac{1}{n'} \sum\_{k=1}^{(N\_D-1)n'} u\_{ij} \left(t - \frac{\left(k-1\right)\pi}{\left(N\_D-1\right)n'}\right) \tag{5}$$

Irikura (1986) introduced the number n ′ of re-division to remove the effect of artificial periodicity due to the equal-size element division. **Figure 1B** presents the slip function D(t) and slip correction function f(t) revised by Irikura (1986) and Yokoi and Irikura (1991). Irikura (1994) introduced the following constraint in the setting of n ′ .

$$\frac{n'N\_D}{\pi} > 2f\_H \tag{6}$$

where f<sup>H</sup> is the upper bound of the effective frequency.

Finally the ground motion displacement U (t) due to the whole fault can be expressed in terms of the ground motion displacements uij (t) due to the fault elements.

$$\begin{aligned} U(t) &= \sum\_{i}^{N\_L} \sum\_{j}^{N\_W} u\_{ij} \left( t - t\_{ij} \right) \\ &+ \sum\_{i}^{N\_L} \sum\_{j}^{N\_W} \sum\_{k}^{(N\_D - 1)n'} \frac{1}{n'} u\_{ij} \left( t - t\_{ij} - \frac{\left( k - 1 \right) \pi}{\left( N\_D - 1 \right) n'} \right) \end{aligned} \tag{7}$$

The concept of the stochastic Green's function method used in the present study is illustrated in **Figure 2**.

Assuming that the fault rupture develops in a concentrically, tij can be expressed by

$$\begin{split} t\_{ij} &= t\_{p\ ij} + t\_{r\ ij} \\ &= \frac{r\_{ij}}{\beta} + \frac{\eta\_{ij}}{V\_r} \end{split} \tag{8}$$

where tpij: the propagation time from the fault element to the recording point at the earthquake bedrock, trij: the slip initiation time in the fault element (the slip initiation time of the initiating point = 0), rij: the distance from the fault element to the recording point at the earthquake bedrock, ηij: the distance from the slip initiation point in the whole fault to the fault element, β : the shear wave velocity of the ground, V<sup>r</sup> : the slip propagation speed in the fault.

#### Small Ground Motion From Element Fault

A small ground motion (acceleration) at the earthquake bedrock due to the slip of a fault element can be derived by setting a point source at the center of the fault element (Boore, 1983). The Fourier amplitude spectrum of the ground motion acceleration at the earthquake bedrock can be expressed by

$$\left|A\_{Sij}\left(\boldsymbol{\omega}\right)\right| = \text{Source}\left(\boldsymbol{\omega}\right) \cdot \text{Pass}\left(\boldsymbol{\omega}\right) \tag{9}$$

where Source (ω) is the term related to the source (fault) and Pass(ω) is the term related to the wave attenuation in the pass from the fault element to the earthquake bedrock.

In Boore (1983), Source (ω) and Pass(ω) are set as

$$\text{Source}\left(\omega\right) = \frac{R\_{\theta\phi}\text{FS}\cdot\text{PRTITN}}{4\pi\rho\beta^3} \left|\omega^2\dot{M}\left(\omega\right)\right|\tag{10}$$

$$\text{Pass}\left(\omega\right) = \frac{1}{r\_{\vec{ij}}} \exp\left[-\frac{\pi f r\_{\vec{ij}}}{Q\left(f\right)\beta}\right] \tag{11}$$

where Rθφ: radiation pattern coefficient, FS: amplification due to the free surface (=2), PRTITN: reduction factor that accounts for the partitioning of energy into two horizontal components, ρ : mass density of earthquake bedrock, β : shear wave velocity of earthquake bedrock, Q f : Q-value. Furthermore, M˙ (ω) is the source spectrum and is expressed by

$$\left| \dot{M} \left( \omega \right) \left( 2 \pi f \right)^{2} \right| = \frac{M\_{0\dot{\imath}\jmath} \left( 2 \pi f \right)^{2}}{1 + \left( f / f\_{c\ \dot{\jmath}} \right)^{2}} P \left( f, f\_{\max} \right) \tag{12}$$

where M0ij is the seismic moment of the fault element ij and fc ij is the corner frequency of the fault element ij. P f , fmax is a filter for reducing the higher frequency components and is expressed by

$$P\left(f, f\_{\text{max}}\right) = \frac{1}{\sqrt{1 + \left(f/f\_{\text{max}}\right)^m}}\tag{13}$$

where fmax is the cut-off frequency for higher frequency components and m = 4 is assumed according to Boore (1983).

In this paper, the phase difference method due to Yamane and Nagahashi (2008) is used for expressing the phase of ground motion. The standard deviation of the phase difference due to the fault element ijcan be expressed by

$$
\sigma\_{\rm ij}/\pi = 0.06 + 0.0003r\_{\rm ij} \tag{14}
$$

This relation refers to inland earthquakes (Makita et al., 2018). In this paper, the near-fault ground motion is assumed in which the effect of the rupture directivity is small. With this standard deviation of the phase difference, the phase spectrum is described by

$$\begin{aligned} \phi\_{k+1\ i\ j} &= \phi\_{k\ i\ j} + \Delta\phi\_{k\ i\ j} \quad \text{( $k = 1, 2, \ldots, N/2 - 1$ )}\\ \Delta\phi\_{i\ j} &= -\left(\mu + s \cdot \sigma\_{i\ j}\right) \end{aligned} \tag{15}$$

where φk ij is the k-th phase spectrum of the fault element ij and 1φk ij is the k-th phase difference spectrum of the fault element ij. N is the number of adopted frequencies. Furthermore, µ is the mean of the phase difference and s is the Gaussian random number with 0 mean and unit standard deviation. In this paper, a constant value of µ is assumed in all the fault elements.

The Fourier transform ASij (ω) of the acceleration aSij (t) at the earthquake bedrock due to the fault element ij can be expressed by

$$A\_{\mathbb{S}\circ j}(\omega) = \left| A\_{\mathbb{S}\circ j}(\omega) \right| \cdot e^{i\phi\_{\mathbb{S}\circ}(\omega)} \tag{16}$$

TABLE 1 | Similarities and differences between this research and benchmark test.


The inverse Fourier transform of AS ij (ω) leads to the acceleration aS ij (t) at the earthquake bedrock due to the fault element ij. Finally, the substitution of this aS ij (t) into Equation (7) (difference of displacement and acceleration does not matter) provides the acceleration a (t) at the earthquake bedrock due to the whole fault.

#### VERIFICATION OF THE METHOD OF GROUND MOTION GENERATION USING PHASE DIFFERENCE METHOD

The benchmark test was conducted by Kato et al. (2011) and the model S21 is used for comparison. The benchmark test uses an empirical envelope function of acceleration time histories. On the other hand, the method proposed in this paper employs the phase difference method for representing the phase. Therefore, the proposed method was validated against the above benchmark test. **Table 1** shows the similarities and differences between the proposed method and the method used in the benchmark test.

The fault plane and three recording points for the model S21 used in the benchmark test Kato et al. (2011) are shown in **Figure 3**. The recording points are three points (a), (b), (c). The fault plane is assumed to be vertical and the fault type is the rightlateral strike-slip fault. The fault length=8000m, fault width = 4,000 m, fault slip quantity = 1 m, the seismic moment =M<sup>0</sup> = 1.04 × 1018Nm, strike angleθ , dip angle δ and rake angle λ are (90◦ , 90◦ , 180◦ ). The hypocenter is located at (0, 1,000, 4,000 m) and the fault rupture is propagated concentrically with rupture velocity V<sup>r</sup> = 3000(m/s). The hypocenter of each sub-fault is assumed to be located at the center.

In Somerville et al. (1999), Eshelby (1957), and Brune (1970), the area S(km<sup>2</sup> ) of the fault, the stress drop 1σ of large earthquakes and the corner frequency f<sup>c</sup> are described by the following equations:

$$S = 2.23 \times \left(M\_0 \times 10^7\right)^{3/2} \times 10^{-15} \tag{17}$$

$$
\Delta \sigma = \frac{7}{16} \frac{M\_0}{R^3} \times 10^{-14} \tag{18}
$$

$$f\_c = 4.9 \cdot 10^6 V\_s \left(\frac{\Delta \sigma}{M\_0}\right)^{1/3} \tag{19}$$

where R(km) isthe effective radius (S = πR 2 ). In these equations, the unit of V<sup>s</sup> is km/s, that of 1σ is bar and that of M<sup>0</sup> is dyne-cm.

From Equations (17–19), 1σ = 13.95(Mpa) and f<sup>c</sup> = 0.404(Hz) are calculated, then τ = 2/f<sup>c</sup> ≈ 5.0(s) is set from Boore (1983). The soil conditions are summarized in **Table 2** and the amplification of the ground motion is evaluated by one-dimensional wave propagation theory.

The fault plane is divided into N<sup>W</sup> × N<sup>L</sup> elements. N<sup>W</sup> = 4 is set in the fault width direction and N<sup>L</sup> = 8 is set in the fault length direction. The area of sub-fault is S <sup>S</sup> = 1 (km<sup>2</sup> ). The seismic moment in each fault element (M0S) is 5.40 × 1015Nm and the stress drop (1σS) is assumed to be 13.95(Mpa). The slip D<sup>S</sup> of each sub-fault is 0.167(m) from M0<sup>S</sup> = µSSD<sup>S</sup> and N<sup>D</sup> = 6 from the ratio of fault plane to sub-fault (1/0.167 m). Thus the seismic moment after superimposing the small earthquakes (M<sup>0</sup> ′ ) is calculated as M<sup>0</sup> ′ = N<sup>W</sup> · N<sup>L</sup> · N<sup>D</sup> · M0<sup>S</sup> ≈ 1.04 × 1018(Nm), which is the same as M0. The corner frequency (fcS) is 2.33 Hz from Equation (19) and the radiation pattern (Rθφ) is set to 0.63, which is a uniform value in the frequency domain. As for the phase angle, the standard deviation of phase differences (σij/π ) are calculated from Equation (14) and its mean µ/π in each point is set to −0.140 at Point (a), −0.125 at Point (b) and −0.130 at Point (c). Regarding the horizontal component of superimposing wave, only the SH wave is generated by setting PRRITEN to 1 for simplification. Each small earthquake is generated by disassembling into the NS direction component and the EW direction component. **Table 3** summarizes the source parameters of the fault plane and sub-faults.

As described above, the amplification of ground motion at the above control points of the soil surface is evaluated by onedimensional wave propagation theory. In this model, the number of layer is one and the transfer function for describing ground amplification is defined by the following equation:

$$H\_G\left(\omega\right) = \frac{1}{\cos kH + i\alpha \sin kH} \tag{20}$$

where k, H, and α are the complex wave number, the thickness of layer and the complex impedance, defined by the following equations:

$$k = \omega \sqrt{\rho/G^\*}, \ G^\* = (1 + 2\xi i)\mathcal{G}, \ \alpha = \sqrt{\rho\_1 \mathcal{G}\_1^{\*\ast}}/\sqrt{\rho\_2 \mathcal{G}\_2^{\*\ast}} \tag{21}$$

G and G ∗ are the shear modulus and the complex shear modulus. Furthermore ξ is the hysteretic damping ratio of soil. In the benchmark test, ξ =0 is given. It is noted that, since the radiation damping is taken into account at the earthquake bedrock, the amplification divergence does not occur.

**Figure 4** shows the comparison between the proposed method and the abovementioned benchmark test (Kato et al., 2011). The upper one in **Figure 4** presents the acceleration at the free ground surface at three points. The lower one in **Figure 4** illustrates the pseudo velocity response spectrum. The numbers 1, 2, 3 in figure legend indicate the difference of uniform random numbers for phase angles in Hisada (Kato et al., 2011) and the difference of Gaussian random numbers for phase difference (Equation 15) in Makita et al. (2018). It can be observed from these figures that, while the acceleration time histories exhibit somewhat different properties, especially in its envelope, the pseudo velocity response spectra of both approaches correspond fairly well. This result supports the validity of the method used in this paper. A lessdamped response (acceleration time history) by the proposed

TABLE 2 | Soil conditions.


method with respect to the benchmark test case may result from the fact that, while an envelope function is used in the benchmark test, such function is not used in the proposed method. However, such difference does not cause serious difference in the structural response because the envelope function influences only the initial and ending parts of acceleration time histories with slight effect on the maximum structural response.

# CRITICAL FAULT RUPTURE SLIP MAXIMIZING THE STRUCTURAL RESPONSE

## Concept of Critical Setting of Fault Rupture Slip Distribution and Fault Rupture Front Maximizing the Structural Response at the Earthquake Bedrock and Free-Ground Surface

The soil model and the fault model treated in section Verification of the Method of Ground Motion Generation Using Phase Difference Method (**Tables 2**, **3**) are used again in this section. Although the soil model used in this section (the same as in the benchmark test) seems rather simple, it is noted that the principal objective of this paper is to pay attention to the influence of the fault rupture process on the response of structures on the surface ground. More detailed examination of the effect of the soil properties above the earthquake bedrock will be made in the future as discussed in the previous paper (Makita et al., 2018).

**Figure 5** presentsthe conceptual diagram of the critical setting of the fault rupture slip distribution Dij and the fault rupture front


maximizing the structural response at the earthquake bedrock and free-ground surface (Case A: elastic SDOF model at the earthquake bedrock, Case B: elastic SDOF model at free-ground surface). T is the natural period of the SDOF model and h is the damping ratio. The fault rupture front includes the fault rupture initiation time trij (related to the rupture propagation velocity in the fault) and the rise time τij of the slip in each fault element. More specifically, trij and τij are treated as independent uncertain parameters in the latter uncertainty modeling. The uncertainty in the fault rupture slip distribution (quantity of slip) for the fixed fault rupture front (concentrically) is treated in section Critical Setting of Fault Rupture Slip Distribution and the uncertainty in the fault rupture front (trij and τij) for the fixed fault rupture slip distribution is dealt with in section Critical Setting of Fault Rupture Front.

A genetic algorithm (GA) has been used for optimization (Goldberg, 1989), i.e., the maximization of the response for uncertain parameters. In this paper, a candidate model of the fault rupture slip distribution or the fault rupture front are treated as chromosomes, and the parameters of each fault element

are treated as genes. First, we generate a number of candidate models (first generation), in which fault parameters are changed randomly. Then, we evaluate these models and generate the next generation by selecting elite individuals, do mutation and conduct crossover. In this GA, the Elitist expected value model is used in which the population size is 200, the number of elite individuals is 2 and the probability of crossover is 0.8. It is noted that global and local search of the optimal solution is possible via GA.

# Critical Setting of Fault Rupture Slip Distribution

The quantities of the fault rupture slip in fault elements are selected as uncertain parameters. The fault rupture initiation time trij and the rise time τij of the slip in each fault element are fixed to the nominal values in this section, i.e., the fault rupture develops from the initiation point concentrically. The ( )<sup>C</sup> denotes the nominal value and α is the uncertain parameter. The

TABLE 4 | Parameters of interval analysis.


interval parameters of the fault rupture slip can be expressed by

$$\mathbf{D}^{I} = \left\{ \left[ D\_{\vec{ij}}{}^{C} - \alpha \Delta D\_{\vec{ij}}{}^{}{D} D\_{\vec{ij}}{}^{C} + \alpha \Delta \vec{D}\_{\vec{ij}} \right] \right\} \left( i = 1, \cdots, N\_{W}, j = 1, \cdots, N\_{L} \right) \tag{22}$$

The over-bar indicates the upper-side value and the under-bar does the lower-side value. The parameters of the interval analysis are shown in **Table 4**.

In this section, the quantities of the fault rupture slip in fault elements are varying in accordance with the following condition.

$$\begin{aligned} \mathbf{D}^l &= [D^\ell - \alpha \Delta \underline{Q} \; \; D^\ell + \alpha \Delta \bar{D}] \\ &= [0.167(\text{m}) - 0.3 \times 0.167(\text{m}) \; 0.167(\text{m}) + 0.3 \times 0.167(\text{m})] \\ &= [0.117(\text{m}) \; 0.217(\text{m})] \end{aligned}$$

Furthermore, the rise time is set to τ = W/(2Vr) = 0.67s from Day (1982) because τ = 5.0s is too long in this model as shown in Kato et al. (2011).

#### Time–History of Wave

**Figure 6** shows the critical ground surface acceleration and deformation of the SDOF model (T = 0.5 s) at three points for Case A (earthquake bedrock motion) and Case B (free-ground surface motion) with respect to uncertain fault rupture slip distribution. Furthermore, **Figure 7** presents the critical ground surface acceleration and deformation of SDOF model (T=1.0, 2.0 s) at three points for Case A (earthquake bedrock motion) and Case B (free-ground surface motion).

It can be observed from **Figures 6**, **7** that the amplification of the ground motion acceleration and the deformation response of the SDOF model of T = 0.5, 1.0 s is larger than those of T = 2.0 s. This means that the effect of the criticality in the uncertainty of the fault rupture slip is larger in the model of shorter natural periods T = 0.5, 1.0 s. In other words, the deviation of structural response between the critical and the nominal case is larger for the natural periods of the SDOF model at T = 0.5 s and T = 1.0 s compared to the case where T = 2.0 s.

#### Fourier Amplitude Spectrum

**Figure 8** shows the Fourier amplitude of critical ground-surface acceleration at three points for three SDOF models (T = 0.5, 1.0, 2.0 s) for Case A (earthquake bedrock motion) and Case B (freeground surface motion). The broken line indicates the natural frequency of the SDOF model. It can be observed that the Fourier amplitude of critical ground-surface acceleration is amplified much around the natural frequency of the SDOF model. This

phenomenon is remarkable in the SDOF model of T = 0.5 s at Point (c). It may be concluded that the critical setting of the fault rupture slip quantity makes the SDOF model resonant to the input.

#### Phase Difference Distribution

**Figure 9** presents the phase difference distribution of critical ground-surface acceleration at three points for three SDOF models (T = 0.5, 1.0, 2.0 s) for Case A (earthquake bedrock motion) and Case B (free-ground surface motion). It can be seen that the standard deviation σ/π at Point (c) is the smallest and that at Point (a) is the largest. This may be related to the forward directivity effect. Furthermore σ/π of the critical model is larger than that of the nominal model at Point (b) and that is smaller than that of the nominal model at Point (c) except (i) of Case B. In addition, µ/π of the critical model becomes larger than that of the nominal model. This means that the critical setting of the fault rupture slip quantity makes the acceleration time history delayed.

#### Fault Rupture Slip Distribution

**Figure 10** illustrates the fault rupture slip distribution maximizing the response of the three SDOF models (T = 0.5, 1.0, 2.0 s) at three points for Case A (earthquake bedrock motion) and Case B (free-ground surface motion). The seismic moment is indicated at the top of the figures. It can be observed that the seismic moment of the critical fault rupture distribution exhibits a value close to the value for the nominal model M<sup>0</sup> = 1.04 × 1018Nm. It can also be seen that the critical fault rupture distributions are different in Case A and Case B. This may be because the site amplification is included in Case B in the process of criticality.

**Figure 11A** shows wave superimposing time tij (from the fault rupture initiation in the fault to the arrival at the earthquake bedrock) for each fault element at three points and **Figure 11B** presents the grouping of fault elements with large slip maximizing the response of two SDOF models (T = 1.0, 2.0 s) at three points for Case B (free-ground surface motion). The triangle in **Figure 11A** indicates the rupture initiation point and the numbers above **Figure 11B** indicate the mean of wave superimposing time tij at grouping fault elements. It can be observed from **Figure 11B** that the mean of wave superimposing time tij at grouping fault elements is slightly shorter than the natural period of the SDOF model. This

may lead to the fact that the acceleration input reflecting the critical fault rupture slip distribution contains the component resonant to the natural period of the SDOF model and amplifies the structural response. In other words, the Fourier amplitude spectrum at the ground surface is amplified in the frequency range of the critical mean of wave superimposing time at grouping fault elements resonant to the natural period of the SDOF model. In addition, the slip of the fault element with the wave superimposing time of 3.71 s becomes large and this may induce a pulse-type wave. Furthermore, the slip distribution in **Figure 11B** may be regarded as an asperity distribution and this distribution can be used as a tool for setting an asperity in the characteristic model of the fault rupture.

## Critical Setting of Fault Rupture Front

The fault rupture initiation time trij and the rise time τij of the slip in each fault element are selected as uncertain parameters. The quantities of the fault rupture slip in fault elements are fixed to the nominal values in this section. ( )<sup>C</sup> denotes the nominal value and α is the uncertain parameter. The interval parameters of the fault rupture initiation time and the rise time of the slip can be expressed by

$$\mathbf{t\_r}^I = \left\{ \left[ t\_{r\,j}{}^C - \alpha \Delta \underline{t}\_{\,\,r\,\dot{j}}, \ t\_{r\,\dot{j}}{}^C + \alpha \Delta \overline{t}\_{r\,\dot{j}} \right] \right\} \,\, \left( i = 1, \cdots, N\_W, \, j = 1, \cdots, N\_L \right) \tag{24}$$

$$\begin{aligned} \pi\_{\mathbf{r}}^{\,I} = \left\{ \left[ \begin{aligned} \pi\_{ij}{}^{\,C} - \alpha \,\Delta \mathbb{1}\_{ij}, & \pi\_{ij}{}^{\,C} + \alpha \,\Delta \bar{\pi}\_{ij} \right] \right\} & \left( i = 1, \cdots, N\_W, \, j = 1, \cdots, N\_L \right) \end{aligned} \tag{25}$$

FIGURE 11 | Wave superimposing time t ij for each fault element and grouping of fault elements with large slip, (A) Wave superimposing time for each fault element at three points for Case A and Case B, (B) Grouping of fault elements with large slip maximizing the response of two SDOF models (T = 1.0, 2.0 s) at three points for Case B (free-ground surface motion).

where

$$
\Delta \underline{t}\_{r \, \dot{\underline{ij}}} = t\_{r \, \dot{\underline{ij}}} ^C , \ \Delta \bar{t}\_{r \, \dot{\underline{ij}}} = t\_{r \, \dot{\underline{ij}}} ^C
$$
 
$$
\Delta \underline{\underline{x}}\_{\, \dot{\underline{ij}}} = \Delta \underline{x} = \tau\_{\dot{\underline{ij}}} ^C , \ \Delta \bar{\underline{x}}\_{\, \dot{\underline{ij}}} = \Delta \bar{\underline{x}} = \tau\_{\dot{\underline{ij}}} ^C
$$

As explained before, the over-bar indicates the upper-side value and the under-bar does the lower-side value. Furthermore 1( ) indicates the parameter for normalization of variation. The parameters for interval analysis is shown in **Table 5**.

#### Time–History

**Figure 12** presents the critical ground surface acceleration and deformation of three SDOF models (T = 0.5, 1.0, 2.0 s) at three points for Case B (free-ground surface motion) with respect to uncertain fault rupture front. It can be observed that the critical TABLE 5 | Parameters of interval analysis.


acceleration record at Point (c) is amplified much as a pulse-type one and those at Point (a) and (b) are also amplified largely. It can also be seen that, while the uncertainty in the quantity of the fault rupture slip increases the response for the nominal

parameters up to about two times as observed in the previous figures, the uncertainty in the rupture propagation velocity in the fault and the rise time of the slip amplifies the response several times. This phenomenon is remarkable in the model of T = 0.5, 1.0 s. However, this is less obvious for T = 2.0 s at Point (c).

# Fourier Amplitude of Ground-Surface Acceleration

**Figure 13A** presents the Fourier amplitude of ground-surface acceleration at three points maximized for fault rupture front for three SDOF models (T = 0.5, 1.0, 2.0 s) and Case B

<sup>(</sup>free-ground surface motion), (A) Fourier amplitude, (B) Phase difference distribution.

(maximization for free-ground surface motion). The broken line indicates the natural frequency of the SDOF model. It can be observed that, as seen in the case of the uncertainty in the fault rupture slip, the Fourier amplitude of ground-surface acceleration is amplified much around the natural frequency of the SDOF model. This phenomenon is remarkable at Point (c) and this may be related to the fact that the critical acceleration record at Point (c) is amplified much as a pulse-type one in **Figure 12**. Furthermore, the amplification is larger than that considering the uncertainty in the quantity of the fault rupture slip. It should be remarked that the natural period of the SDOF model does not coincide with the amplified range of the Fourier amplitude for the model T = 2.0s at Point (c). This may cause the lower response amplification of the model T = 2.0s at Point (c) in **Figure 12**(iii) than other models.

#### Phase Difference Distribution

**Figure 13B** shows the phase difference distribution of groundsurface acceleration at three points maximized for fault rupture front for three SDOF models (T = 0.5, 1.0, 2.0 s) and Case B (free-ground surface motion). It can be seen that the standard deviation σ/π of the critical one is larger than that of the nominal one. On the other hand, the standard deviation σ/π of the critical one becomes smaller than that of the nominal one at Point (c). It may result from the fact that, while the duration time at Point (a) and (b) becomes longer, the duration at Point (c) becomes shorter as a result of pulsetype ground motion. In addition, µ/π becomes larger as a result of criticalization of the rupture front (the uncertainty in the rupture propagation velocity in the fault and the rise time of the slip). This phenomenon was also observed in considering the uncertainty in the quantity of the fault rupture slip.

#### Critical Fault Rupture Front

**Figure 14** presents the wave superimposing time tij (from the fault rupture initiation in the fault to the arrival at the earthquake bedrock) at each fault element and the rise time τij at each fault element characterizing the critical fault rupture front maximized for three points and three SDOF models (T = 0.5, 1.0, 2.0 s) corresponding to Case B (free-ground surface motion). The triangle in **Figure 14** indicates the rupture initiation point. It can be observed from **Figure 14A** that the fault element rupture occurs in a concentrated manner corresponding to the time interval of the natural period of the SDOF model. In addition, the rupture directivity effect can be observed in the model of T = 1.0, 2.0 s at Point (c) and it may be related to the fact that the pulse-type input induces larger responses. Furthermore, it can be observed from **Figure 14B** that, while the rate of τij around the nominal value is large in the model of T = 0.5 s, τij moves to the lower limit around τij = 0.47–0.49 s in the model of T = 1.0, 2.0 s.

#### Robustness Evaluation for Uncertain Fault Rupture Slip Distribution and Uncertain Fault Rupture Front

**Figure 15A** shows the robustness function αˆ, proposed by Ben-Haim (2006), with respect to the deformation of the SDOF model for uncertain parameters of quantity of fault rupture slip for Case (B). Once the value αˆ in the vertical axis is fixed, the corresponding deformation of the SDOF model in the horizontal axis indicates the maximum value for varied uncertain parameters (quantity of fault rupture slip) prescribed by αˆ. In particular, the deformation of the SDOF model for αˆ = 0 indicates the maximum response for the nominal parameters. It can be observed that the robustness becomes the smallest for the model at Point (b). This is because the response of the SDOF model is the largest at Point (b). The slope of the robustness function indicates the degree of the robustness. As the slope becomes steeper, the model becomes more robust.

**Figure 15B** presents the robustness function αˆ with respect to the deformation of the SDOF model for uncertain parameters of fault rupture front (slip initiation time and rise time) for Case (B). Once the value αˆ in the vertical axis is fixed, the corresponding deformation of the SDOF model in the horizontal axis indicates the maximum value for varied uncertain parameters (slip initiation time and rise time) prescribed by αˆ. It can be observed that the robustness becomes the smallest for the model at Point (b) as in **Figure 15A**. Compared to the case in **Figure 15A**, the slope of the robustness function becomes small.

Since the robustness is closely related to the resilience, the presented method using the robustness function seems useful for the evaluation of resilience of buildings against uncertain fault rupture slip distribution and uncertain fault rupture front.

## CONCLUSIONS

To promote a new methodology for resilient building design, a critical excitation method has been proposed in which the whole process of theoretical ground motion generation is treated. The process consists of (i) the fault rupture process, (ii) the wave propagation from the fault to the earthquake bedrock, (iii) the site amplification. The uncertainty in the fault rupture slip has been dealt with in the present paper, i.e., the quantity of the fault rupture slip, the rupture propagation velocity in the fault and the rise time of the slip. The wave propagation from the fault to the earthquake bedrock has been expressed by the stochastic Green's function method in which the Fourier amplitude at the earthquake bedrock from a fault element has been represented by the Boore's model and the phase angle has been modeled by the phase difference method. The validity of the proposed method has been investigated through the comparison with the existing simulation result by other methods incorporating the empirical envelope function of acceleration time histories. By using the proposed method for ground motion generation and for optimization under uncertainty in the fault rupture slip, a methodology has been presented for deriving the critical ground motion causing the maximum response of an elastic SDOF model at the earthquake bedrock or at the free ground surface. A genetic algorithm (GA) has been used for optimization, i.e., the maximization of the response for uncertain parameters. The following conclusions have been derived.


Since the critical ground motion produces the worst building response among possible scenarios, the proposed method can be a reliable tool for resilient building design.

## AUTHOR CONTRIBUTIONS

KM formulated the problem, conducted the computation, and wrote the paper. KK conducted the computation and discussed the results. IT supervised the research and wrote the paper.

# ACKNOWLEDGMENTS

Part of the present work is supported by KAKENHI of Japan Society for the Promotion of Science (No. 17K18922, 18H01584). This support is greatly appreciated.

#### Makita et al. Uncertainty in Fault Rupture Slip

## REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Makita, Kondo and Takewaki. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Analytical Model for Multi-Hazard Resilient Prefabricated Concrete Frame Considering Earthquake and Column Removal Scenarios

Kaiqi Lin<sup>1</sup> , Xinzheng Lu<sup>2</sup> \*, Yi Li <sup>3</sup> , Weidong Zhuo<sup>1</sup> and Lieping Ye<sup>2</sup>

<sup>1</sup> College of civil engineering, Fuzhou University, Fuzhou, China, <sup>2</sup> Beijing Engineering Research Center of Steel and Concrete Composite Structures, Tsinghua University, Beijing, China, <sup>3</sup> Key Laboratory of Urban Security and Disaster Engineering of Ministry of Education, Beijing University of Technology, Beijing, China

#### Edited by:

Ehsan Noroozinejad Farsangi, Graduate University of Advanced Technology, Iran

#### Reviewed by:

Emanuele Brunesi, Fondazione Eucentre, Italy Abbas Sivandi-Pour, Graduate University of Advanced Technology, Iran Christian Málaga-Chuquitaype, Imperial College London, United Kingdom

> \*Correspondence: Xinzheng Lu luxz@tsinghua.edu.cn

#### Specialty section:

This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment

Received: 29 September 2018 Accepted: 21 November 2018 Published: 04 December 2018

#### Citation:

Lin K, Lu X, Li Y, Zhuo W and Ye L (2018) Analytical Model for Multi-Hazard Resilient Prefabricated Concrete Frame Considering Earthquake and Column Removal Scenarios. Front. Built Environ. 4:73. doi: 10.3389/fbuil.2018.00073

Research on multi-hazard prevention and mitigation in building structures is the most recent developing trend in civil engineering. In this study, an analytical model is proposed to calculate the structural resistance of a type of multi-hazard resilient prefabricated concrete (MHRPC) frame under earthquake and column removal scenarios. The MHRPC frame is assembled using prefabricated RC beams and columns, unbonded post-tensioning (PT) tendons, energy-dissipating steel angles, and large rotational shear plates. According to the experimental results, the MHRPC frame exhibits the features of low damage and self-centering under seismic loading. Meanwhile, when subjected to column removal scenarios, the MHRPC frame is proven to demonstrate a high progressive collapse resistance. In order to calculate the seismic and progressive collapse resistance of the MHRPC frame, analytical models for the critical components in the MHRPC frame (PT tendons and steel angles) are compared and selected based on the experimental results and numerical simulations. Furthermore, calculation methods for the seismic and progressive collapse resistance of the MHRPC frame specimens are proposed. The calculation results are validated using the experimental results. This study could provide a reference for the design of MHRPC frame structures, considering both earthquake and progressive collapse.

Keywords: calculation model, multi-hazard, prefabricated concrete frame, earthquake, progressive collapse

# INTRODUCTION

Over the past several years, increased attention has been paid to multi-hazard mitigation and prevention of building structures in the engineering community. Li et al. (2011) reviewed the state-of-art research on the multi-hazard from the perspective of (1) damages and loses, (2) assessment of effects, and (3) design and mitigation strategies. The importance of life-cycle and multi-hazard design are highlighted. Gidaris et al. (2017) reviewed the multi-hazard fragility and restoration models of highway bridges for the risk and resilience assessment of regional portfolios and transportation networks. Kamath et al. (2015) and Shah et al. (2016) evaluated the residual lateral resistance of a two-story RC frame after the post-earthquake fire. The effectiveness of ductile detailing in improving the structural residual strength is validated. A multi-hazard resistant bridge pier considering earthquake and explosion load is proposed and experimentally tested by Fujikura et al. (2008). ElSayed et al. (2015) suggested the layout of seismically detailed reinforcement in concrete-block shear walls for resisting the blast load. For the most commonly constructed multi-story RC frame structures, numerous existing studies have proven that earthquake-induced collapse and progressive collapse starting from a local failure are the major failure modes of multi-story RC frames (Sozen et al., 1998; Lu et al., 2012). Hence, increasing the resistance capacity of seismic and progressive collapse is critical for improving the safety margin and collapse resistance of RC frames.

Progressive collapse of a building structure refers to the disproportionate chain collapse action of a structure, initiated by a small and localized failure that may be caused by fire, explosion or overloading (Ellingwood, 2006). A typical progressive collapse example of an RC frame is the 1995 bombing of the Murrah Federal Building in Oklahoma City (Sozen et al., 1998). Thereafter, numerous countries published progressive collapse design requirements for RC frames.

In fact, the seismic and progressive collapse designs for RC frames vary significantly in terms of the design methodology. Seismic design aims to resist a system-level lateral load and realize a "strong-column-weak-beam" failure mode under an earthquake. Hence, the seismic resistance of frame columns is critical to the seismic performance of an RC frame. In contrast, progressive collapse design needs to bridge local unbalanced vertical loads by enhancing lateral components to redistribute the unbalanced gravity load, and avoiding the initial failure propagation. According to Lin et al. (2017), an RC frame with a relatively low seismic design intensity can hardly meet the requirements of progressive collapse design, and its beams should be strengthened to prevent progressive collapse. However, after the progressive collapse design, it is found that newly added progressive collapse reinforcements in the beams may lead to an unfavorable "strong-beam-weak-column" failure mode, which will, in turn, weaken the structural seismic performance.

In order to improve the seismic and progressive collapse resistance of newly-designed RC frame structures simultaneously, a novel multi-hazard resilient prefabricated concrete (MHRPC) frame, incorporating a series of highperformance components, namely post-tensioning (PT) tendons, energy-dissipating steel angles, and shear plates, is proposed, as illustrated in **Figures 1A,B** (Lin et al., in press). Furthermore, the seismic and progressive collapse performance of the newly proposed MHRPC frame was experimentally compared with a conventional RC frame by means of seismic cyclic tests of beam-column joint specimens and progressive collapse tests of two-span substructures, as illustrated in **Figures 1C,D**. The test setups are shown in **Figure 1E**. The results indicate that, compared to the conventional RC frame, the MHRPC frame specimen exhibits substantially smaller residual deformations and less component damage following the seismic cyclic test. During the progressive collapse test, the MHRPC specimen exhibits a significantly higher progressive collapse resistance than the conventional RC specimen, and meets the chord rotational capacity requirement as stipulated in Department of Defense (2016), which demonstrates superior progressive collapse resistance. It is concluded that the MHRPC frame system provides a satisfying solution for improving the seismic and progressive collapse resistance of RC frame.

Although many experimental, numerical, and analytical studies have been conducted on the earthquake resilient RC structures (Priestley and Tao, 1993; Wolski et al., 2009; Gerami et al., 2013; Fakharifar et al., 2014; Gerami and Sivand-Pour, 2014; Song et al., 2014, 2015; Lu et al., 2015). Limited work has been reported on the progressive collapse performance of this type of structures. Moreover, existing analytical models for the progressive collapse resistance calculation are only suitable for conventional RC frames (Yi et al., 2008; Lu et al., 2018). Note that the proposed MHPRC frame is comprised of a series of high-performance components, the resistance contribution of different components needs to be quantified under both earthquake and column removal scenarios. Hence, it is necessary to propose an accurate and easy-to-implement analytical model for the seismic and progressive collapse designs of the MHRPC frame system.

Based on the experimental results, an analytical model is proposed in this study to calculate the structural resistance of the MHRPC frame, considering earthquake and column removal scenarios analytically. The results are compared to the experimental ones and exhibit strong agreement. This study could provide a reference for the design of multi-hazard resilient RC frame structures.

# RESISTANCE CONTRIBUTIONS OF KEY COMPONENTS IN MHRPC FRAME

# Seismic Cyclic Test

According to Lin et al. (in press), the measurement results of the PT tendons and steel angles in the seismic cyclic tests are illustrated in **Figure 2**. The results indicate that the tendon force increased as the joint rotation amplitude increased (**Figure 2A**). Meanwhile, the steel angles were found to yield, and dissipated energy once the joint rotation reached 1.6% (**Figure 2B**). Compared to the steel angles, the shear plates remained elastic during the seismic cyclic tests (Lin et al., in press).

# Progressive Collapse Test

The measurement results of the PT tendons, steel angles, and shear plates in the progressive collapse test are illustrated in **Figure 3**. According to the experimental observations of Lin et al. (in press), the loading process of the MHRPC specimen can be divided into the beam mechanism and catenary mechanism stages. During the beam mechanism stage, progressive collapse resistance is provided by the compressive arch action (CAA) and flexural capacities of the beams. The displacement corresponding to the peak resistance at this stage is defined as D<sup>b</sup> as shown in **Figure 3B**. During the catenary mechanism stage, the progressive collapse resistance originates from the catenary force of the PT tendons and steel angles. The displacement corresponding to code required chord rotation (i.e., 0.20 rad) is defined as D0.20 as shown in **Figure 3B** (Department of Defense, 2016). According to the strain gauge reading on the shear plate, the shear plate acts

FIGURE 1 | Schematic and test results of the MHRPC frame. (A) Deformation of MHRPC frame under column removal scenario. (B) Details of beam-column joint. (C) Comparison of the seismic cyclic performance between MHRPC and conventional RC frames. (D) Comparison of the progressive collapse performance between MHRPC and conventional RC frames. (E) Test setups (Lin et al., in press).

as a cantilever beam and resists the shear force transferred from the bolt during the test (Lin et al., in press).

It can be concluded that, for the proposed MHRPC frame, the energy-dissipating steel angles and PT tendons served as the key load-resisting components in both the seismic cyclic and progressive collapse tests. The shear plate was responsible for transferring part of the shear force at the beam ends and ensuring a large chord rotational capacity under the column removal scenario. Therefore, this study focuses on analyzing the resistance contributions of the PT tendons and steel angles. Based on the component-level analyses, an analytical model is developed to calculate the seismic cyclic and progressive collapse resistance of the tested specimens.

#### ANALYTICAL MODEL FOR PT TENDON

In this study, referring to the Design Specification for Unbonded Post-Tensioned Precast Concrete Special Moment Frames Satisfying ACI 374.1 and Commentary (ACI 550.3- 13) (American Concrete Institute, 2013), the PT tendon model proposed by Mattock (1979) is adopted. According to Mattock (1979), the stress-strain relationship of a PT tendon can be expressed by Equation (1):

$$f\_s = E\varepsilon \left\{ 0.020 + \overline{\left[ 1 + \left( \frac{E\varepsilon}{1.04\text{\textpi}} \right)^{8.36} \right]^{\frac{1}{8.36}}} \right\} \tag{1}$$

where E is the elastic modulus of the PT tendon; f py is the yield strength of the PT tendon, which can be approximated as 90% of the ultimate strength; ε is the strain of the PT tendon, which can be calculated by ε = ε<sup>0</sup> + 1ε ; ε<sup>0</sup> is the initial strain. According to the material test carried out by Lin et al. (in press), the tensile strength of the PT tendon is 1,993 MPa.

In the tests on the MHRPC frame specimens, as the prefabricated beams and columns are covered by steel sleeves at the ends, the deformation of the PT tendon is assumed to arise primarily from the relatively rigid body rotations between the prefabricated components. Hence, the incremental strain (1ε ) can be calculated by analyzing the specimen deformation modes. Thereafter, the resistance contribution of the PT tendons can be determined using Equation (1).

#### ANALYTICAL MODEL FOR STEEL ANGLE

In the MHRPC frame, the steel angle connections aid in: (1) providing flexural strength and dissipating seismic energy in the seismic cyclic tests; and (2) providing a catenary force during the catenary mechanism stage of the progressive collapse tests. It should be noted that, in the catenary stage, the steel angles will experience large deformations. Hence, the analytical model of the steel angle should be able to calculate the load-carrying capacity in both the small and large deformation stages.

The connection between the steel angles and prefabricated components is identical to the widely used semi-rigid connection with top and seat angles. The steel angles are arranged around the beam-column joint region, and connected to the prefabricated beams and columns by means of high-strength bolts. As semirigid connections with top and seat angles are used extensively in steel structures, numerous experimental, numerical and analytical studies exist (Kishi and Chen, 1990; Calado and Ferreira, 1994; Mander et al., 1994; Bernuzzi et al., 1996; Ahmed et al., 2001; Kishi et al., 2001; Garlock et al., 2003; Komuro et al., 2004; Li, 2007; Yuan, 2007; Yang and Jeon, 2009; Ahmed and Hasan, 2015; Hasan et al., 2017; Kong and Kim, 2017). Among these researches, Kishi and Chen (1990) proposed calculation methods for the initial stiffness and ultimate moment capacity of different types of semi-rigid connections, including the top and seat angle connection. Mander et al. (1994) calculated the plastic moment capacity of the connection by means of virtual work principles. For the initial stiffness, the Kishi and Chen model 1990 was adopted by Mander et al. (1994) in their research. Another widely used steel angle connection model was proposed by Garlock et al. (2003), based on cyclic load tests on different steel angles.

## Key Parameters of Steel Angle Connections

When the beam-column joint is subjected to a rotation of θ, the prefabricated beam rotates along point O and the deformation of top and seat angles are shown in **Figure 4A**. Under such a load, the top angle is subjected to a horizontal force V as shown in **Figure 4B** and the plastic hinges locate at the column side of the steel angle. Furthermore, the beam side of the seat angle rotates along the plastic hinge on the steel angle as shown in **Figure 4C**.

According to Kishi and Chen (1990), Mander et al. (1994), and Garlock et al. (2003), when the top angle is subjected to a horizontal force V at the beam side, its strength is affected by the following parameters: plastic hinge distance g, angle thickness t, and angle width b, among others. The distance between two plastic hinges is determined by the bolt size and location and the steel angle filet length. In order to provide an improved description of the analytical model for the steel angle, a series of parameters for calculating the steel angle resistance (g1, g1', g2, g2', r, L, t, and L') are defined, as illustrated in **Figure 4D**. Moreover, the geometric parameters of a steel angle connection are defined, including the steel angle leg length (l<sup>1</sup> and l2), bolt diameter (d), and bolt head diameter (D). Based on the deformation mode in **Figure 4C**, the seat angle strength depends primarily on the flexural resistance at the plastic hinge location.

#### Analytical Model for Steel Angles

The analytical model for the steel angle connections includes calculation methods for the yield moment, initial stiffness, and post-yield resistance. Hence, the proposed model is introduced in terms of these three aspects.

#### Yield Moment

(1) The Kishi and Chen (1990) and Garlock et al. (2003) method (namely, Method Y1)

When the beam-column joint reaches the yield state, according to the moment equilibrium between the external force (My) and internal force, the yield moment can be expressed as follows (Kishi and Chen, 1990; Garlock et al., 2003):

$$M\_{\circ} = M\_{\text{set}} + M\_{\text{\textit{p}}} + V\_{\text{\textit{p}}}h \tag{2}$$

where Mseat is the flexural strength at plastic hinge 3 (**Figure 4C**), which can be calculated by Equation (3); M<sup>p</sup> is the flexural strength at plastic hinge 2 (**Figure 4B**); V<sup>p</sup> is the shear force transferred from the frame beam; and h is the beam height.

Assuming that both angle legs have the same thickness t, the flexural strength of a steel angle with a width b and material strength f <sup>y</sup> can be calculated according to Equation (3).

$$M\_{\mathbb{P}} = \frac{f\_{\mathbb{Y}} b t^2}{4} \tag{3}$$

Moreover, the shear force V<sup>p</sup> can be calculated using Equation (4), as suggested by Garlock et al. (2003):

$$V\_{\mathcal{P}} = \frac{2M\_{\mathcal{P}}}{\mathcal{g}}\tag{4}$$

(2) The Mander et al. (1994) method (namely, Method Y2)

In contrast to Kishi and Chen (1990) and Garlock et al. (2003), Mander et al. (1994) deduced the expression for the yield moment by means of the virtual work principle, as indicated in Equation (5):

$$M\_{\rm Y} = m\_3 + m\_2 + (m\_1 + m\_2)\frac{h'}{\mathcal{g}}\tag{5}$$

where m1, m2, and m<sup>3</sup> are the flexural strength at plastic hinges 1, 2, and 3, respectively, which can also be calculated by Equation (3); h ′ is the distance between plastic hinge 2 and the rotation center O; and g is the distance between plastic hinges 1 and 2.

#### Initial Stiffness

The initial stiffness of the top and seat angle connection can be derived by calculating the initial stiffness of the top and seat angles, respectively. The initial stiffness contribution of the seat angle can be calculated by Equation (6) (Mander et al., 1994).

$$K\_{\text{set}} = \frac{4EI}{l\_{\text{so}}} \tag{6}$$

where E is the elastic modulus of the steel; I is the sectional moment of inertia; and lsois the distance between plastic hinge 3 and the angle leg tip, as indicated in **Figure 4A**.

Calculation methods for the initial stiffness contribution of the top angle include the Kishi and Chen (1990) and the Garlock et al. (2003) methods.

(1) The Kishi and Chen (1990) method (namely, Method S1)

According to Kishi and Chen (1990), the top angle acts as a cantilever beam and its initial stiffness contribution can be derived according to Equation (7):

$$K\_{\text{top}} = \frac{\Im EI}{1 + \frac{0.78t^2}{\mathcal{g}'^2}} \frac{h\_0^2}{\mathcal{g}'^3} \tag{7}$$

where E is the elastic modulus of the steel; I is the sectional moment of inertia; t is the angle thickness; g ′ is the distance from plastic hinge 1 to the mid-thickness of the angle leg on the beam side; and h<sup>0</sup> is the distance between the mid-thickness of the top and seat angles on the beam side. When the top and seat angles have the same thickness, h<sup>0</sup> = h+t.

(2) The Garlock et al. (2003) method (namely, Method S2)

Garlock et al. (2003) assumed that the top angle is fixed at the bolt positions on both the beam and column sides, as illustrated in **Figure 5**. The initial stiffness contribution of the top angle can be calculated by Equations (8–10). The initial stiffness can be expressed in terms of the bending and shear stiffness (Kbend and Kshear):

$$\frac{1}{K\_{\text{top}}} = \frac{1}{K\_{\text{shear}}} + \frac{1}{K\_{\text{bend}}} \tag{8}$$

$$K\_{\text{shear}} = \frac{12EI}{0.26gt^2} \tag{9}$$

$$K\_{\text{bend}} = \frac{12EI}{\mathcal{g}^3} - \frac{6EIC\_{\theta}}{\mathcal{g}^2} \tag{10}$$

$$\mathcal{C}\_{\theta} = \begin{bmatrix} \frac{\frac{3}{\xi} \left( 1 + 2 \frac{\epsilon}{\xi} \right)}{\frac{2}{\xi} \left( 1 + \frac{3 \epsilon}{2 \xi} \right) + \frac{2}{L} \left( 1 + \frac{3 \epsilon}{2L} \right)} \end{bmatrix} \tag{11}$$

where Cθ is the rotation angle corresponding to a unit movement of the steel angle heel, which can be derived by Equation (11); e is the half-length of the square rigid zone, which is equal to t/2; and g and L are two steel angle connection parameters, as defined in **Figure 5**. As suggested by Garlock et al. (2003), when calculating the initial stiffness contribution of the top angle, g<sup>1</sup> ′ and L ′ as depicted in **Figure 4D** are assigned to g and L, respectively. Taking the beam height into consideration, the flexural stiffness contribution of the top angle can be determined.

#### Post-yield Resistance

The shear force of the top angle following yielding can be calculated by Equation (12), according to Garlock et al. (2003):

$$V = \left(\frac{2M\_{\mathcal{P}}}{g - \Delta}\right)\alpha\tag{12}$$

where ∆ is the movement of the top angle heel, as illustrated in **Figure 5**; and α is the material hardening parameter of the steel.

It should be noted that Kishi and Chen (1990) also proposed a power-law-based model to simulate the post-yield behavior of the steel angle connection. However, the power index in the Kishi and Chen method 1990 should be calibrated by means of a series of experimental data, which is not practical in engineering design; therefore, it is not adopted in this study.

#### Model Selection for Steel Angle Steel Angle Database

In order to the compare the accuracies of the different methods mentioned above, a total of 45 steel angle connection specimens were collected from the literature (Calado and Ferreira, 1994; Mander et al., 1994; Bernuzzi et al., 1996; Garlock et al., 2003; Komuro et al., 2004; Li, 2007; Yuan, 2007; Yang and Jeon, 2009; Ahmed and Hasan, 2015). Detailed information regarding these specimens can be found in the database in **Appendix Table 1** in Supplementary Material.

Moreover, a benchmark FE model for the steel angle connection in the experiments of Lin et al. (in press) is constructed using MSC.Marc. Based on the benchmark model, 24 FE models with various design parameters (bolt positions, angle thicknesses, and material strengths) are built to establish the numerical database in **Appendix Table 2** in Supplementary Material.

The cyclic tests on steel angles conducted by Garlock et al. (2003) are used to validate the accuracy and feasibility of the FE models. Taking specimen L6-516-9 as an example, solid elements are used to construct the models. The steel angle model is divided into four layers along the thickness. The material parameters are assigned according to the tension coupon results provided by Garlock et al. (2003). The simulated cyclic response of specimen L6-516-9 is compared to the experimental result in **Figure 6**, which indicate that the above FE model can accurately reproduce the cyclic behaviors of the steel angle in the tests.

Based on the above experimental and numerical database, different analytical models for the steel angle connection are compared and selected. Subsequently, the selected models are used to analyze the experimental results of the MHRPC specimens.

#### Model Selection

#### **Yield strength**

The ratios between the calculated yield strengths and experimental results for the 45 specimens in the experimental database are compared in **Figure 7A**. The blue hollow triangle data points indicate the results from Method Y1 (Kishi and Chen, 1990; Garlock et al., 2003), while the red solid diamond data points represent the results from Method Y2 (Mander et al., 1994). The mean absolute errors (MAEs) of Methods Y1 and Y2 are represented by the blue and red solid lines, respectively. Note that Method Y2 is only suitable for beam-column joint specimens with top and seat angle connections, and the corresponding data points for specimens 1–7 (Garlock et al., 2003) are absent as these tests are cyclic tests of the steel angles.

The results indicate that both methods can provide an accurate estimation of the yield strength of the steel angle connections. The MAE of Method Y1 is 16.42% with a standard deviation of 0.163, while the MAE of Method Y2 is 30.91% with a standard deviation of 0.197. Moreover, additional data points from Method Y1 are within the range of ±20%, as illustrated in **Figure 7A**.

The calculated yield strengths using Method Y1 are compared to the numerical database results in **Figure 7B**. The results demonstrate that Method Y1, validated by the experimental database, can accurately calculate the yield strength of the steel angle connection model. The MAE is 7.44% with a standard deviation of 0.064. Hence, Method Y1 is adopted in this study to calculate the yield strengths of the steel angle connections.

#### **Initial stiffness**

Similarly, the ratios between the calculated initial stiffness and experimental results for the specimens in the experimental database are illustrated in **Figure 8A**. The blue hollow triangle data points denote the results from Method S1 (Kishi and Chen, 1990), while the red solid diamond data points represent Method S2 (Garlock et al., 2003). The MAEs of Methods S1 and S2 are indicated by the blue and red solid lines, respectively. The MAE of Method S1 is 89.61% with a standard deviation of 1.356, while that of Method S2 is 133% with a standard deviation of 2.251. Both Methods S1 and S2 can provide accurate predictions for certain specimens. In contrast, neither of the two models can provide an effective prediction for the remaining specimens. It should be noted that such errors may arise from either the analytical models or experimental measurement errors. The errors of the analytical models may come from the idealization of the connection, measured material strength and so on. The experimental measurement errors may be caused by the installation of the tested specimens (e.g., relative slide, content of stiffening, machining error, etc.), the accuracy of the apparatus and so on. Actually, errors larger than 500% are also found in existing literature when predicting the initial stiffness of the top and seat angle connections (Kong and Kim, 2017). Hence a numerical database is necessary for selecting a more accurate initial stiffness calculation method. Similar model selection methodology is also adopted in many existing studies (Hasan et al., 2017; Kong and Kim, 2017).

Compared to the experimental results, the FE model results can avoid errors from the experimental measurement. For all of the models in the numerical database, the calculated initial stiffness values are compared to the numerical results in **Figure 8B**. It is demonstrated that the calculated initial stiffness of Method S1 is slightly smaller than the numerical result. The MAE of Method S1 is 12.80% with a standard deviation of 0.145. In contrast, the result of Method S2 is higher than the corresponding numerical results. The MAE is 61.59% with a standard deviation of 0.247. Based on the comparison of the two methods, Method S1 is suggested for calculating the initial stiffness of the top angle in this study.

To conclude, the proposed calculation procedure for the top and seat angle connection is summarized in **Figure 9A**. Taking model T6 in **Appendix Table 2** in Supplementary Material as an example, the calculated load-displacement curve using the proposed method is compared to the numerical results, with strong agreement, as indicated in **Figure 9B**.

# ANALYTICAL MODEL FOR MHRPC FRAME Seismic Resistance

According to El-Sheikh et al. (2000), during the initial loading stage of the seismic cyclic test, the deformation of the beamcolumn joint arises from the flexural deformation of the prefabricated concrete beam until reaching the linear limit state. The initial stiffness is defined as the secant stiffness when the component reaches the linear limit state. El-Sheikh et al. (2000) suggested that the linear limit state moment could be determined according to Equations (13–15).

$$M\_{\rm ll} = \min\left(M\_{\rm ll1}, M\_{\rm ll2}\right) \tag{13}$$

$$M\_{\rm ll1} = 5f\_{\rm pi}A\_{\rm p}h/12\tag{14}$$

$$M\_{\rm ll2} = \bar{m} \left( 1 - \frac{f\_{\rm ci}}{0.85 f'\_{\rm c}} \right) = 0.5 f\_{\rm pi} A\_{\rm p} h \left( 1 - \frac{f\_{\rm ci}}{0.85 f'\_{\rm c}} \right) \tag{15}$$

where Mllis the linear limit state moment corresponding to the frame beam on one side of the beam-column joint; T is the pretension force of a PT tendon; f pi is the initial stress of a PT tendon; A<sup>p</sup> is the total area of the PT tendons; h is the beam

height; f ci is the initial stress of the concrete under pre-tension; and f <sup>c</sup> ′ is the concrete cylinder strength. Note that, in the seismic cyclic tests of Lin et al. (in press), as the specimen is loaded at the beam ends on both sides of the joint, the calculated linear limit state should be 2Mll.

Before the specimen reaches the linear limit state, the initial stiffness of the moment-rotation relationship for the MHRPC joint specimens can be approximated by Equation (16).

$$R = 2\frac{EI\_\text{g}}{L} \tag{16}$$

where E is the material elastic modulus, and for simplicity, the elastic modulus of concrete is used in this study; I<sup>g</sup> is the sectional moment inertia; and L is the effective length of the frame beam on one side of the joint, which can be taken as the distance between the edges of two steel jackets. In the seismic cyclic tests of Lin et al. (in press), L = 1.2 m. According to the drawings provided in Lin et al. (in press), the calculated initial stiffness is 8,138 kN·m for the MHRPC joint specimen. Moreover, in the first test of this specimen (initial pre-stressing level: 42%), the calculated linear limit state moment is 34.7 kN·m, corresponding to a joint rotation of 0.426%. In the second test (initial pre-stressing level: 20%), the linear limit state moment is 16.7 kN·m, corresponding to a joint rotation of 0.205%.

According to Lin et al. (in press), in the seismic cyclic test of the MHRPC frame joint specimen, the prefabricated beams rotate along the corner points at the beam-column interfaces, as illustrated in **Figure 4A**. The steel angles begin to deform and contribute to the flexural resistance once the specimen reaches the linear limit state. The deformation of the steel angle heel (∆) can be calculated according to Equation (17):

$$
\Delta = \left(\theta - \theta\_{\text{ll}}\right) h \tag{17}
$$

where θ is the joint rotation and θll is the joint rotation corresponding to the linear limit state. Note that, when θ is smaller than θll, ∆ = 0. h is the beam height. The resistance from the top and seat angle connection can be calculated according the procedure illustrated in **Figure 9A**.

Moreover, the elongation of the top PT tendon (∆tendon) can be calculated by Equation (18), according to the deformation mode of the beam-column joint, as illustrated in **Figure 10**:

$$
\Delta\_{\text{tendon}} = \theta \left( d\_1 + d\_2 \right) = \theta h \tag{18}
$$

where d<sup>1</sup> and d<sup>2</sup> are the distances from the PT tendon centers to the bottom of the prefabricated beam, respectively. Note that, in the MHPRC frame, the PT tendons are arranged symmetrically along the sectional height, and the elongation of the bottom PT tendon can also be calculated by Equation (18). When

considering the total length of the PT tendon (ltendon), the strain increment in Equation (1) can be calculated following 1ε = ∆tendon/ltendon, based on which the internal force of the PT tendon can be determined. Furthermore, the flexural resistance contribution of the PT tendons can be obtained.

The moment-rotation relationship of the MHRPC joint specimen can be approximated by summing the flexural resistance contributions from the PT tendons and steel angle connections. Note that, in the seismic cyclic tests of Lin et al. (in press), the MHRPC joint specimen is tested twice in order to verify the reparability and self-centering capacity. During the first test, the initial pre-stressing level is set as 42%, while in the second test, it is 20%. The analytical backbone curves are compared to the experimental results in **Figure 11**. It should be noted that the calculated tendon force-rotation relationships in the two tests are also compared to the experimental measurements in **Figure 11**. The results indicate that the calculated cyclic response of the MHRPC joint specimen, including the tendon force-rotation relationship, fits strongly with the test results. The proposed method can be used in the seismic resistance design and analysis of such MHRPC frame structures.

#### Progressive Collapse Test

As discussed previously, according to the load-displacement curve in **Figure 1D**, the loading process of the progressive collapse test of the MHRPC substructure can be divided into the beam mechanism and catenary mechanism stages. Hence, the analytical model also includes the beam mechanism and catenary mechanism parts.

#### Beam Mechanism

The progressive collapse resistance of the MHRPC substructure mainly originates from the CAA and flexural resistance of the frame beam. The calculation method proposed by Lu et al. (2018) is adopted in this study. The peak displacement corresponding to the peak CAA resistance can be estimated by Equation (19).

$$
\delta = 0.00050l^2/h \tag{19}
$$

where l is the total specimen length and h is the beam height. Furthermore, the progressive collapse resistance can be calculated following the Park and Gamble (2000) method (Equations 20–25).

$$P = \frac{2\left(M\_1 + M\_2 - N\delta\right)}{\beta l} \tag{20}$$

$$N = C\_\text{c} + C\_\text{s} - T = 0.85f'\_\text{c}\beta\_1 cb + C\_\text{s} - T \tag{21}$$

$$M\_1 = 0.85f'.\beta\_1 \varepsilon' b \left(0.5h - 0.5\beta\_1 \varepsilon'\right) + C\_\text{s}'(0.5h - d')$$

$$\begin{split} M\_1 &= 0.85f'\_{\text{c}}\beta\_1 \mathbf{c}' b \left( 0.5h - 0.5\beta\_1 \mathbf{c}' \right) + C'\_s (0.5h - d') \\ &+ T' (0.5h - d') \end{split} \tag{22}$$

$$M\_2 = 0.85f'\_{\text{c}}\beta\_1 c b \left(0.5h - 0.5\beta\_1 c\right) + C\_\text{s} (0.5h - d')$$

FIGURE 11 | Comparison of seismic cyclic test results. (A) First test of MHRPC joint specimen (initial pre-stressing level: 42%). (B) Second test of MHRPC joint specimen (initial pre-stressing level: 20%).

$$+\ \text{T(0.5h - d')}\tag{23}$$

$$\varepsilon' = \frac{h}{2} - \frac{\delta}{4} - \frac{\beta l^2}{4\delta} (\varepsilon + \frac{2t}{l}) + \frac{T' - T - C\_s' + C\_s}{1.7 f'\_{\varsigma} \beta\_1 b} \tag{24}$$

$$\mathcal{L} = \frac{h}{2} - \frac{\delta}{4} - \frac{\beta l^2}{4\delta} (\varepsilon + \frac{2t}{l}) - \frac{T' - T - C\_s' + C\_s}{1.7 f'\_{\text{c}} \beta\_1 b} \tag{25}$$

where the parameter definitions are illustrated in **Figure 12**. A more detailed CAA calculation procedure in RC frame beams can be found in Lu et al. (2018).

The calculated CAA resistance of the progressive collapse test specimen in Lin et al. (in press) is 58.53 kN. Combined with the previously estimated peak displacement, a typical loading point can be denoted as illustrated in **Figure 15A**.

#### Catenary Mechanism

During the catenary mechanism stage, the MHRPC substructure deforms, as illustrated in **Figure 13**. The elongation of a PT tendon can be calculated by Equation (26), according to the

geometry.

$$\Delta\_{\text{tendon}} = 2\left(\sqrt{\delta^2 + l\_{\text{n}}^2} - l\_{\text{n}}\right) \tag{26}$$

where δ is the column stub displacement and l<sup>n</sup> is the net span of the prefabricated beam. Based on the tendon elongation, the strain increment and tendon force can be calculated further. The progressive collapse resistance contribution of the PT

tendons is composed of the vertical components of the tendon forces.

Moreover, according to Yang and Tan (2013), the catenary force contribution of the steel angle connections is dependent on the ultimate tensile force of the steel angles, which can be calculated by Equation (27).

$$N\_{\mathbf{a}} = (b\_{\epsilon\overline{\mathcal{f}},a} - nd\_{b,\text{hole}})t\_{\mathbf{a}}f\_{\mathbf{y}} \tag{27}$$

where N<sup>a</sup> is the ultimate tensile force of the steel angle and beff ,<sup>a</sup> is the effective width of the steel angle, which is determined by Equations (28, 29).

$$b\_{\rm eff,a} = \min\left(d\_{\rm h} + 2m\_{\rm a}, \frac{d\_{\rm h}}{2} + m\_{\rm a} + \frac{w}{2}, \frac{b\_{\rm h}}{2}, \frac{d\_{\rm h}}{2} + m\_{\rm a} + e\_{\rm x}\right) \tag{28}$$

$$m\_{\rm a} = m - t\_{\rm a} - 0.8r\_{\rm a} \tag{29}$$

where the parameter definitions are provided in **Figure 14**.

For the steel angle used in the progressive collapse test of the MHRPC substructure, the calculated beff ,<sup>a</sup> and N<sup>a</sup> values are 62.5 mm and 130.54 kN, respectively. Furthermore, the progressive collapse resistance contribution is composed of the vertical components of tensile forces in the steel angles.

The calculated load-displacement curve of the progressive collapse test is compared to the experimental result, as illustrated in **Figure 15A**. It should be noted that the tendon force-displacement relationship is also compared in **Figure 15B**. The results indicate the proposed method can provide reasonable predictions of the progressive collapse resistances for MHRPC substructures in the beam mechanism and catenary mechanism stages.

#### CONCLUSION

A novel MHRPC frame system has been proposed to improve the seismic and progressive collapse performance of commonly used RC frames by Lin et al. (in press). In this study, based on the experimental tests of the MHRPC frame specimens, an analytical model for the design of the MHRPC frame is established and validated by means of experimental results. The main contributions of this study include the following.


This work could provide a reference for the seismic and progressive collapse design of such MHRPC frame structures.

#### AUTHOR CONTRIBUTIONS

KL and XL contributed development of the proposed model and performed the analysis. KL wrote the first draft of the manuscript; KL, XL, YL, WZ, and LY wrote sections of the manuscript. All authors read and approved the submitted manuscript version.

## ACKNOWLEDGMENTS

The authors are grateful for the financial support received from the National Natural Science Foundation of China (No. 51778341).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbuil. 2018.00073/full#supplementary-material

# REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Lin, Lu, Li, Zhuo and Ye. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Finite Difference Method-Based Critical Ground Motion and Robustness Evaluation for Long-Period Building Structures Under Uncertainty in Fault Rupture

#### Koki Makita, Kyoichiro Kondo and Izuru Takewaki\*

*Department of Architecture and Architectural Engineering, Graduate School of Engineering, Kyoto University, Kyoto, Japan*

#### Edited by:

*Ehsan Noroozinejad Farsangi, Graduate University of Advanced Technology, Iran*

#### Reviewed by:

*Aleksandra Bogdanovic, Institute of Earthquake Engineering and Engineering Seismology, Macedonia Rita Bento, Universidade de Lisboa, Portugal*

\*Correspondence:

*Izuru Takewaki takewaki@archi.kyoto-u.ac.jp*

#### Specialty section:

*This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment*

Received: *26 November 2018* Accepted: *03 January 2019* Published: *22 January 2019*

#### Citation:

*Makita K, Kondo K and Takewaki I (2019) Finite Difference Method-Based Critical Ground Motion and Robustness Evaluation for Long-Period Building Structures Under Uncertainty in Fault Rupture. Front. Built Environ. 5:2. doi: 10.3389/fbuil.2019.00002*

It is known that, while the stochastic Green's function method is suitable for generating ground motions with short periods, the three-dimensional finite difference method (FDM) is appropriate for ground motions with rather long periods. In the previous research, the stochastic Green's function method was used for finding the critical earthquake ground motion for variable fault rupture slip (slip distribution and rupture front). However, it cannot be used for ground with irregularities and for wave component with rather long periods. In responding to this request, the FDM is used in this paper for finding the critical ground motion for structures with rather long natural period. Since the FDM is time-consuming, it seems unfavorable to use it in a simple sensitivity algorithm where an independent response sensitivity is calculated for many design parameters. To overcome this difficulty, the bi-cubic spline interpolation of uncertain parameter distributions (seismic moment distribution in this paper) and the response surface method for predicting the response surface from some sampling points are used effectively in this paper. The uncertainty parameter is the fault rupture slip distribution described in terms of seismic moments. It is found that the critical ground motion for structures with rather long natural period can be found effectively by the proposed method.

Keywords: critical ground motion, worst input, fault rupture, finite difference method (FDM), response surface method, spline interpolation, long-period structure, robustness

# INTRODUCTION

There was common understanding that earthquake ground motions can be classified into a few types (Abrahamson et al., 1998). However, many earthquake ground motions of peculiar characteristics have been observed recently in the world (for example, Mexico, 1985; Northridge, 1994; Kobe, 1995; Chi-chi, 1999; Tohoku, 2011). After such earthquakes occur, we feel that unpredictable ground motions can occur and powerful theoretical approaches are inevitable to respond to those ground motions. One of the approaches may be the critical excitation method (see Drenick, 1970; Takewaki, 2007) which enables the search of the worst earthquake ground motion among possible ones. To tackle the worst ground motion under the consideration of fault rupture models, tools for producing ground motions in terms of fault rupture models may be necessary (Makita et al., 2018b).

It is well recognized that, while the empirical Green's function method (Irikura, 1986; Yokoi and Irikura, 1991) or the stochastic Green's function method (Wennerberg, 1990; Hisada, 2008) is suitable for generating ground motions with short periods, the three-dimensional finite difference method (FDM) is appropriate for ground motions with rather long periods (Bouchon, 1981; Hisada and Bielak, 2003; Yoshimura et al., 2003; Nickman et al., 2013). To enhance the usability of the FDM, an open software (GMS: Ground Motion Simulator) is available (Aoi and Fujiwara, 1999; Aoi et al., 2008; Maeda et al., 2012, 2016; Tanaka et al., 2014). The combination of these two-type motions with the use of a matching filter is acknowledged as the most powerful and reliable method for generating earthquake ground motions under the consideration of fault ruptures and surface waves. Since the parameters used in these methods for ground motion generation contain various uncertainties, i.e., aleatory uncertainty and epistemic uncertainty (Taniguchi and Takewaki, 2015; Okada et al., 2016), the treatment of such uncertainties are essential for the reliable estimation of ground motions (Abrahamson et al., 1998; Lawrence Livermore National Laboratory, 2002; Morikawa et al., 2008; Cotton et al., 2013).

In the previous research (Makita et al., 2018a), the effect of the fault rupture was taken into account simply by introducing the phase difference method. The robustness of a new building structural system consisting of base-isolation and building connection (Murase et al., 2013) was investigated for uncertain ground models. However, the fault rupture mechanisms cannot be considered in detail. In another previous research (Makita et al., 2018b), the stochastic Green's function method was used for finding the critical earthquake ground motion for variable fault rupture slip (slip distribution and rupture front). However, it cannot be used for ground with irregularities and for wave component with rather long periods.

Since the FDM is time-consuming, it seems unfavorable to use it in a simple sensitivity algorithm where an independent response sensitivity is calculated for each design parameter. To overcome this difficulty, the bi-cubic spline interpolation of uncertain parameter distributions (seismic moment distribution in this paper) through the control points and the response surface method for predicting the response surface from some sampling points are used effectively in this paper. The uncertain parameter is the fault rupture slip distribution described in terms of seismic moments at the control points. It is found that the critical ground motion for building structures with rather long natural period can be found effectively by the proposed method.

To investigate the effect of uncertainty level in the fault rupture on the robustness of building structures using the robustness function (Ben-Haim, 2006), several uncertainty levels are set and the critical fault rupture model is sought. Then the maximum story ductility is obtained for each uncertainty level.

## OUTLINE OF PROPOSED METHOD

In this paper, the uncertain parameter is the fault rupture slip distribution described in terms of seismic moments M0. First of all, the nominal distribution of fault rupture slip is given. The rupture front is assumed to develop concentrically as shown in **Figure 1A**. To reduce the degree of freedom in the setting of variable fault rupture slip distribution, some control points are selected in the fault. Sampling points in the uncertain parameter range are planned by the experimental design method at all control points. Then the uncertain parameters at all points are interpolated from the values at the control points by introducing bi-cubic spline interpolation of uncertain parameter distributions (seismic moment distribution in this paper). In the next, ground motions are generated by using the FDM. The response of a structure under the generated ground motion is computed. Then the response surface is obtained by the least-squares method. The maximum value of the response surface is determined by using the Sequential Quadratic Programming (SQP) method. Finally, the critical fault model is found and the corresponding ground motion is generated. The earthquake response analysis is conducted under the critical ground motion. The flowchart for finding the critical earthquake ground motion is shown in **Figure 1B**. To investigate the robustness of a building structure with rather long natural period with respect to the variability of the fault rupture distribution, the robustness function due to Ben-Haim (2006) is introduced. The relation of the maximum response of the structure with the uncertainty level of variable fault rupture slip distribution provides the quantitative evaluation of the robustness of the structure against the uncertain environment.

As mentioned above, in the generation of ground motions, FDM is used for producing ground motions with rather long periods (usually longer than about 1–2s). For ground motion components with shorter periods, the stochastic Green's function method is often used. Then both ground motion components are combined with the matching filter. Since a building structure with a rather long fundamental natural period is treated here, FDM is employed for generating ground motions.

## FINITE DIFFERENCE METHOD (FDM)

The three-dimensional finite difference method (FDM) is often used as a useful numerical method for generating earthquake ground motions on the ground with irregularities. It can also take the fault rupture mechanism into consideration. In the research group of ground motion generation, an open source

bi-cubic spline function. (E) Set the asperity from the obtained seismic moment distribution. (F) Termination of assignment of seismic moment.

(GMS: Ground Motion Simulator: http://www.gms.bosai.go.jp/ GMS/) can be used (Aoi and Fujiwara, 1999; Aoi et al., 2008; Maeda et al., 2012, 2016; Tanaka et al., 2014). In this paper, such open source software is used. The reliability and accuracy of this software will be investigated through the comparison with actual earthquake events and the benchmark tests (see **Appendix**).

# OPTIMIZATION IN PROPOSED METHOD

# Response Surface Method

The response surface method (RSM) is often used as an efficient and reliable method for prediction of responses of structures with many parameters (Khuri and Cornell, 1996). The procedure of the RSM can be summarized as follows: (i) Select the control

points, (ii) Plan sampling points by the experimental design method for the control points, (iii) Interpolate the uncertain parameters at all points from the values at the control points, (iv) Generate ground motion using the FDM, (v) Calculate the response of a structure under the generated ground motion (vi) Estimate the response surface by the least-squares method, (vii) Search the maximum value of the response surface using the SQP method.

While the earthquake ground motions have to be generated repeatedly in the design procedure based on the conventional critical excitation method after the change of design conditions (the change of uncertainty level in the fault as treated in this paper or the change of superstructures etc.), those do not need to be generated in the design procedure based on the proposed critical excitation method using the RSM. This is because the earthquake ground motions for given sampling points have been generated and those can be used repeatedly.

Let x<sup>i</sup> denote the seismic moment at the control point i. The response surface in terms of quadratic functions can be expressed as

$$\mathbf{y} = \sum\_{i} c'\_{i} \mathbf{x}\_{i} + \sum\_{i} c\_{i} \mathbf{x}\_{i}^{2} + \sum\_{i} \sum\_{j \neq i} c\_{ij} \mathbf{x}\_{i} \mathbf{x}\_{j} + c\_{0} + \varepsilon \tag{1}$$

where the first term is a linear term, the second term is a quadratic term, the third term is a cross term, the fourth term c<sup>0</sup> is a constant, and the fifth term ε is an error term.c<sup>i</sup> , c ′ i , cij are their coefficients.

Although the second-order approximation is said to be inferior to the third-order approximation, it has some merits, (i) the required number of sampling points for a given accuracy is small, (ii) the solution is stable, (iii) the computational load for the increasing number of input factors is within a reasonable range. The coefficientsc<sup>i</sup> , c ′ i , cij, c<sup>0</sup> are determined by using the well-known least-squares method.

In the next, let us explain the sampling method for the second-order approximation. The representative methods are (i) Latin Hypercube Sampling (LHS), (ii) Central Composite Design (CCD), (iii) Box-Behnken Design (BBD) (see Box and Behnken, 1960). Among these sampling methods, CCD, and BBD are designed for the evaluation of the second-order approximation. **Figure 2A** shows the CCD sampling method for three parameters. Each axis indicates the variation of the corresponding uncertain parameter. In the CCD method, the necessary number of sampling points for n uncertain parameters is (n + 2)(n + 1)/2 (Ohbuchi et al., 2011). In this paper, CCD method is employed.

In CCD, three types of sampling points exist, (i) Central point, (ii) Epaxial point, (iii) Factorial design. Let M0ij and M¯ <sup>0</sup>ij denote the seismic moment and its nominal value. The central point indicates a nominal value. The Epaxial point, (M0ij/M¯ <sup>0</sup>ij) − 1, is on an axis and it varies the range (−1, 1) as shown in **Figure 2B**. **Figure 2B** shows an example such that the value only at the point (1, 5) varies. Since the Factorial design is intended to interact with each other, it makes each parameter vary ±1/ √ n as shown in **Figure 2C**. The objective of the Factorial design is to investigate the interaction between parameters.

#### Seismic Moment Distribution Using Spline Interpolation

If the number of divisions in the fault plane is small, the FDM cannot simulate the smooth fault rupture and keep the computational accuracy in a wide frequency range. When we consider many uncertain parameters in a fault plane, the robustness evaluation needs formidable amount of computational load. Therefore, some techniques are needed to reduce the computational load.

The seismic moment distribution at all points is obtained by using the bi-cubic spline interpolation for the seismic moments at the control points. Consider the fault element in the region [pi , pi+1] × [q<sup>j</sup> , qj+1]. The seismic moment in this region can be expressed by

$$f\_{\vec{l}\vec{j}}(p,q) = \sum\_{k,l=0}^{3} a\_{\vec{ij}}^{kl} (p - p\_i)^k (q - q\_j)^l \tag{2}$$

where a kl ij is the coefficient. The method of setting of the bi-cubic spline functions is explained in **Figures 3A–F**. The respective procedures (a)-(f) can be summarized as follows.

(a) Set the rectangular fault model (nominal model). All the fault elements in this nominal model have a constant seismic moment. Select some control points.

#### TABLE 1 | Soil conditions.


The constraint at the stage (f) is introduced following the research by Ishii et al. (2000) and Somerville et al. (1999). Ishii et al. (2000) defined 70% from the largest of the fault elements as the principal rupture region and Somerville et al. (1999) reported that the mean area of the asperities in inland earthquakes is 22%.

The rise time at the above-mentioned stage (f) is set following the research by Day (1982).

$$
\pi = W / (2V\_r) \tag{3}
$$

The width of the asperity is substituted into Equation (3). In this paper, the area S<sup>a</sup> of the asperity is calculated at the stage (f) and a square fault is assumed. Then the equivalent width W<sup>a</sup> is substituted into Equation (3).

Compared to the previous works, the proposed method enables the reduction of the number of uncertain parameters and the smooth setting of the parameters on the fault.

#### Constraint on Parameter Variation Using nth-Order Hypercube and Hypersphere

Let **x** = {x1, · · · , xn} T and **x** denote a set of uncertain parameters x<sup>i</sup> and their nominal set. n is the number of uncertain parameters (the number of control points in this paper). In the uncertainty

analysis, the simplest constraint on parameter variation is a box type described by

$$R\_1(\mathbf{x}, \bar{\mathbf{x}}, \alpha) = \left\{ \mathbf{x} \, \middle| \, \left| (\mathbf{x}/\overline{\mathbf{x}}) - 1 \right| \le \alpha \right\} \tag{4}$$

where α is a given value representing the uncertainty level. R1is a hyper cube of order n. Although R<sup>1</sup> is suitable for problems with a relatively larger uncertainty level, it is often the case that the final solution goes to its end of the range. Furthermore, the setting of a larger uncertainty level is apt to cause a heavy computational load. To remedy this, a hyper sphere constraint is often used which can be defined by

$$R\_2(\mathbf{x}, \bar{\mathbf{x}}, \Sigma, \mathcal{c}) = \{ \mathbf{x} | (\mathbf{x} - \bar{\mathbf{x}})^T \Sigma^{-1} (\mathbf{x} - \bar{\mathbf{x}}) \le \mathcal{c} \}\tag{5}$$

where 6 is the covariance matrix. Since (**x** − ¯**x**) <sup>T</sup>6−<sup>1</sup> (**x** − ¯**x**) follows the χ <sup>2</sup> distribution of order n, the probability of **x** ∈ R2(**x**¯, 6, β) becomes Fn(c). Fn(−) is the probability distribution function of the χ <sup>2</sup> distribution of order n. If we define the confidence region by β = Fn(c), Equation (5) can be re-expressed by

$$R\_2(\mathbf{x}, \bar{\mathbf{x}}, \Sigma, F\_n^{-1}(\beta)) = \langle \mathbf{x} | (\mathbf{x} - \bar{\mathbf{x}})^T \Sigma^{-1} (\mathbf{x} - \bar{\mathbf{x}}) \le F\_n^{-1}(\beta) \rangle \tag{6}$$

When **x** follows the normal distribution, Equation (6) indicates that the probability of (**x** − µ) <sup>T</sup>6−<sup>1</sup> (**x** − µ) ≤ F<sup>n</sup> −1 (β) is β . The setting of variable regions in the fault may be possible by substituting the parameters defining the inhomogeneity of the fault parameters into the covariance matrix 6.

There are some researches on the variability of fault parameters (Somerville et al., 1999; Ishii et al., 2000). In this paper, the result by Ishii et al. (2000) on "inland faults" is used. Ishii et al. (2000) investigated the inversion of 15 inland earthquakes (seismic moment, rupture velocity etc.) and derived the mean and the coefficient of variation of the ratio (the principal rupture region to the overall fault) of the seismic moment. In this paper, the mean is µ = 2.1 and the coefficient of variation is CV = 0.36. The obtained variance σ <sup>2</sup> = (2.1 × 0.36)<sup>2</sup> = 0.57 is substituted into the covariance 6 in Equation (6). In this paper, the covariance between different faults is treated as 0.

#### NUMERICAL SIMULATION

Consider a numerical example using the FDM for generating the earthquake ground motions. The critical fault model is investigated by the proposed method. In addition, to evaluate the validity and the degree of criticality of the obtained critical ground motion, two models (Recipe 1 and Recipe 2) are considered based on the strong ground motion prediction recipe (Earthquake Research Committee, 2017). The search of the critical excitation for several levels of uncertainty in the fault is conducted using the proposed method and the corresponding structural responses are clarified.

#### Modeling of Ground and Fault

Consider a quarter grid model of FDM as shown in **Figure 4A**. The fault length and width are 30 km and 18 km. The other parameters are strike =90◦ , dip =90◦ , rake =180◦ , and the original base point is located at (0, −15, 2)(km). The seismic size is assumed to be M<sup>w</sup> = 6.8 and M<sup>0</sup> = 1.8 × 1019(Nm). It is assumed that the rupture propagates concentrically from the rupture initiation point H(Hx, Hy, Hz) = (0, −10.2, 15.8)(km) with the rupture propagation velocity V<sup>r</sup> = 2800(m/s). The time shift tshift in each fault element is given by

$$\begin{aligned} t\_{\text{shift}} &= t\_{\text{start}} + \xi / V\_r \\ \xi &= \sqrt{(X - H\_x)^2 + (Y - H\_y)^2 + (Z - H\_z)^2} \end{aligned} \tag{7}$$

where tstart is the source rupture initiation time.

In the setting of area source in the three-dimensional FDM, it is necessary to approximate this by multiple point sources placed at the difference grid points. In this paper, multiple point sources are placed at the difference grid points on the source layer and the seismic moment is released by considering the time delay due to the rupture propagation. For this purpose, the fault plane is divided into the small fault size dx = 0.6(km) and small faults of 50 × 30 = 1500 are placed on the source layer. To express the sequential rupture of the divided area sources, it is necessary that dx is sufficiently smaller than the wave length λ(km) of the rupture front. The wave length of the rupture front can be computed by λ = Vr/f = 2.8/1 = 2.8(km) with the effective frequency f = 1(Hz). It can be understood that dx ≪ λ is satisfied and the sequential rupture can be expressed in a sufficient manner.

The following triangle function is employed as the source time function in the fault element.

$$f(t) = \begin{cases} 0 & \text{( $t \le -1/2f\_c$ )}\\ 4f\_c^2 t + 2f\_c & \text{( $-1/2f\_c \le t \le 0$ )}\\ -4f\_c^2 t + 2f\_c & \text{( $0 \le t \le 1/2f\_c$ )}\\ 0 & \text{( $1/2f\_c \le t$ )} \end{cases} \tag{8}$$

where f<sup>c</sup> is the inverse of the rise time τ and τ indicates the width of the bottom of the triangle in f(t). The source time function is shown in **Figure 4B**.

The three-dimensional difference grid is set as 120km × 150km × 60km (−60km ≤ X ≤ 60km−75km ≤ Y ≤ 75km0km ≤ Z ≤ 60km) as shown in **Figure 4A**. **Figure 4A** shows 1/4 of the total region. The material properties of the soil layer and the source layer are shown in **Table 1**. In the software "GMS," the inhomogeneous grid is used as the difference grid. The grid interval in the source layer is triple of that in the soil layer. In this paper, the grid interval in the source layer is set so as to satisfy the condition on the effective frequency 0–1 Hz. The fourth-order accurate scheme is used in the difference operator. In the fourth-order accurate scheme, 5– 6 grids are required in one wave length (shear wave). In the soil layer, the one grid length 200 (m) leads to the effective

frequency f ≤ Vs/(5H) = 2000/(5 × 200) =2(Hz). On the other hand, in the source layer, the grid length 600(m) leads to the effective frequency f ≤ Vs/(5H)= 3400/(5 × 600) = 1.13(Hz). These parameters satisfy the condition on the effective frequency. The absorbing zone of 12 km is placed at the side and bottom of the object region to damp the reflected wave as shown in **Figure 4A**. The time duration is 30 (s) and the time increment is 0.015 (s). The number of time steps is 2,000.

#### Modeling of Superstructure

Consider a 20-story steel building frame. The simplest and most efficient model for vibration analysis of building frames is a shear building model. Usually the shear building modeling is completed by doing a static lateral force analysis for obtaining the story shear-drift relation. However, a shear building model with a predetermined stiffness and strength distribution is assumed here for simple presentation of the proposed critical excitation method. The shear building modeling is conducted and the elastic-plastic response is assumed here. The floor mass is 3.0 × 10<sup>6</sup> (kg) and the fundamental natural period is 2.4 (s). The story stiffness distribution is assumed to be trapezoidal. The 2% stiffness-proportional structural damping is assumed. The story shear-drift relation is assumed to be bilinear and the yield drift is assumed to be 0.02 (m). The post-yield stiffness ratio to the initial stiffness is 0.05. The target story ductility is 2 and the maximum ductility along height is the objective function.

Although a 20-story steel building frame is treated here for a simple presentation of the proposed excitation method, other types of building structures (reinforcedconcrete building structures, taller building structures, etc.) can be dealt with in a similar manner so long as they possess a positive post-yield stiffness. This is because, if building structures with negative post-yield stiffness are treated, the earthquake response may be unstable and the analysis of response sensitivity may cause some difficulties.

story ductility factors for Recipe model 1, 2, Nominal model, and Critical model.

#### Problem Formulation

In the experimental design planning, CCD, the base point on the axis (Epaxial point) is set to 1.0. The factorial design is planned so that the distances of all the sampling points from the base point are equal. 512 points are prepared in the fractional factorial design. The points of the same number 512 were also added as the sampling points from the uniform random numbers in the hyper sphere.

Consider the following optimization problem by the SQP.

$$\begin{array}{ll}\text{Find} & \mathbf{x} \\ \text{which maximizes} & h(\mathbf{x}) \\ \text{subject to} & \mathbf{x} \in R\_1(\mathbf{x}, \overline{\mathbf{x}}, \alpha) \\ & \mathbf{x} \in R\_2(\mathbf{x}, \overline{\mathbf{x}}, \Sigma, F\_n^{-1}(\beta)) \\ & \sum \mathbf{x} = 0 \end{array} \tag{9}$$

In this paper, the maximum story ductility along height is employed as the objective function h(**x**). Following Ishii et al. (2000), α = 1.1 and β = 0.95 are assumed. In case ofR2(**x**, **x**, 6, F<sup>n</sup> −1 (β)), the maximum value of **x** is 1.1 and the maximum norm of **x** is 4.56. The fault model is shown in **Figure 5A**. Twenty-four control points are chosen and the seismic moment ratios M<sup>0</sup> ij/M¯ <sup>0</sup> ij (M<sup>0</sup> ij, seismic moment at the fault element, M¯ <sup>0</sup> ij, seismic moment of the nominal model at the fault element) at the control points are selected as **x**. The seismic moments at the control points are varied in this paper.

To evaluate the validity of the critical ground motion, two models (Recipe 1 and Recipe 2) are considered based on the strong ground motion prediction recipe (Earthquake Research Committee, 2017). The corresponding fault models are shown in **Figures 5B,C**. The areas of Asperity 1, Asperity 2, and back ground of the recipe model are Sa<sup>1</sup> = 90(km<sup>2</sup> ), Sa<sup>2</sup> = 34(km<sup>2</sup> ), S<sup>b</sup> =416(km<sup>2</sup> ), respectively. The seismic moments of those are M0a<sup>1</sup> = 6.73 × 1018(Nm), M0a<sup>2</sup> = 1.55 × 1018(Nm), M0<sup>b</sup> = 9.73 × 1018(Nm). The rise times of those are τa<sup>1</sup> =1.61(sec), τa<sup>2</sup> = 1.07(sec),τ<sup>b</sup> = 3.21(sec).

#### Result

Using the proposed critical excitation method, the critical seismic moment ratios of the critical fault model are obtained. The flowchart of the proposed method is explained in **Figure 1** and the detailed explanation is made in section Optimization in Proposed Method.

**Figure 6A** shows the distribution of critical seismic moment ratios of the critical fault model which produces the maximum story ductility factor. The area ratio of the asperity to the total fault is 22% and the rise time in the asperity is assumed to be τ = 1.95(s). This rise time is slightly longer than that in the recipe model. On the other hand, **Figure 6B** presents the distribution of asperity. It can be observed from **Figure 6B** that the asperity with a large area is produced near the observed site and the seismic moment is concentrated at the bottom and left side in the fault. This phenomenon of multiple asperities was seen in the previous work (Makita et al., 2018b).

**Figure 7** presents the time histories of ground velocity and acceleration (component of transverse) for each fault model.

**Figure 7A** shows the velocities for four models (Recipe model 1, 2, Nominal model, and Critical model). On the other hand, **Figure 7B** presents the accelerations for such four models. A pulse-type wave can be observed in all fault models. In addition, the maximum velocity and the maximum acceleration were observed in the recipe model 1. Furthermore, the amplitude of the afterward wave is large in the acceleration of Recipe model 2 [(b) in (ii)].

**Figure 8** shows the time histories of inter-story drift and story shear for each fault model (component of transverse). The quantities in 2, 4, 6, 8th stories are presented here. It can be observed that, while the response to the nominal model remains elastic, the responses to Recipe model 1, and the critical model become large and go into the plastic range mainly in the lower stories. The fact whether the building model remains elastic or goes into the plastic range can be found from **Figure 9**.

**Figure 9** presents the response for four ground motions. **Figure 9A** shows the hysteretic loops in 2, 4, 6, 8th stories for Recipe model 1, 2, the Nominal model and the Critical model. It can be observed that the size of the hysteretic loop is large in the order of the Critical model, Recipe model 2, Recipe model 1. **Figure 9B** illustrates the maximum story ductility factors for Recipe model 1, 2, the Nominal model, and the Critical model. It can be observed that, while the story ductility factors to Recipe model 1, Recipe model 2, and the Nominal model are within the design limit 2, the story ductility factors to the Critical model are beyond the design limit in 2–5th stories. The story ductility factor to the Critical model is 3.2 times larger than that to the Nominal model and 1.1 times larger than that to Recipe model 2.

#### Robustness Evaluation by Changing Uncertainty Level

In this section, the uncertain parameter α in R<sup>1</sup> is changed to investigate the influence of the uncertainty level in the fault on the robustness in the response of the superstructure. The concept of the robustness function due to Ben-Haim (2006) is used here.

**Figure 10** shows the distribution of seismic moment ratio at each uncertainty level α (1.0–1.5). In the prediction of the fault

model, the response surface obtained in the previous analysis is used. The critical fault model for each uncertainty level α (1.0– 1.5) is searched via the SQP. This enabled efficient analysis of the robustness function. Focusing on the maximum asperity, the maximum value of M0ij/M¯ <sup>0</sup>ij is approximately proportional to α up to α =1.2. On the other hand, when α is larger than 1.3, the maximum value does not change much. This indicates that the proposed method avoids the production of the excessive asperity. In addition, it can be observed that some asperities are located at the edge of the fault except the location near to the observation site.

**Figure 11A** shows the maximum value and norm of the control point parameter for various uncertainty levels α . The maximum value is influenced primarily by R<sup>1</sup> and the maximum norm is influenced primarily by R2. It can be observed that the maximum value attains the corresponding value α and the maximum norm attain the upper limit of R<sup>2</sup> except α = 1.0, 1.4.

**Figure 11B** shows the comparison of ductilities under several ground motions produced for various uncertainty levels. It can be observed that the uncertainty level influences much the responses in the lower stories in this case. **Figure 11C** presents the robustness function αˆ with respect to the ductility (Ben-Haim, 2006). It can be seen that the ductility is over the design criterion (story ductility factor = 2) in the models with the uncertainty level larger than or equal to α = 1.1. This means that the setting of α influences greatly the robustness of the building structure.

# VERIFICATION OF PROPOSED METHOD USING ACTUAL EARTHQUAKE

In this section, the validity and applicability limit of the proposed method are investigated by comparing with the actual earthquake case. The object earthquake is the Osaka North earthquake in 2018 (M<sup>W</sup> = 5.67) (Earthquake Research Committee, 2018).

**Figure 12A** shows the source, observation stations and region of the finite-difference model (GMT). The grid interval is 140 (m). The effective frequency is 0–0.5Hz. The difference grids lower than 30 grids are treated as inhomogeneous ones and the

grid interval is made triple. The energy absorbing zone was set at the bottom and the side (60 grids). **Figure 12B** presents the fault model of this earthquake [(i) Control point for critical model, (ii) Recipe model]. The fault length is 4 km and the fault width is 6 km. The other parameters are strike = 351◦ , dip = 50◦ , rake = 52◦ . The seismic size is M<sup>0</sup> = 4.06 × 1017(Nm). It is assumed

that the fault rupture propagates concentrically from the source with the propagating speedV<sup>r</sup> = 2900(m/s). The limit of the uncertainty level is given by α = 1.1 and β = 0.95. The area of asperity, the seismic moment, the rise time are given by S<sup>a</sup> = 5(km<sup>2</sup> ), M0<sup>a</sup> = 1.79 × 1017(Nm), τ<sup>a</sup> = 0.43(sec),τ<sup>b</sup> = 1.03(sec). **Figure 12C** shows the comparison between the observation wave and the artificial wave by the recipe model. It can be observed that the frequency content is well simulated at OSK002. However, the amplitude is evaluated in a damped manner. On the other hand, at OSK003, the wave can be simulated well until t = 22 s. However, the amplitude does not correspond well after t = 22 s.

**Figure 13** illustrates the comparison of the box-and-whisker plot of the ratio of displacement response to artificial wave to that to the observation wave between the recipe model (**Figure 13A**) and the critical model (**Figure 13B**). It may be concluded that we can produce the critical ground motion under which the displacement response spectrum becomes larger than that under the observed ground motion in the rate of 50% in all natural period ranges and in the rate of 75% in the natural period range up to 13 (s).

# CONCLUSIONS

The finite difference method (FDM)-based critical excitation method has been proposed for building structures with rather long fundamental natural periods. The uncertain parameter is the fault rupture slip distribution described in terms of seismic moments at fault elements. Since the FDM is time-consuming, it is unfavorable to use it in a simple sensitivity algorithm where an independent response sensitivity is calculated for each design parameter (seismic moment at a fault element). To overcome this difficulty, the control points (representative points in the fault) have been introduced. Then the bi-cubic spline interpolation of uncertain parameter distributions (seismic moment distribution) and the response surface method for predicting the response surface from some sampling points have been used effectively. The robustness of building structures for varied uncertainty level of the fault rupture slip (seismic moments in the fault elements) has also been evaluated by using the robustness function. The main conclusions can be summarized as follows.

(1) The introduction of control points in the fault enabled efficient calculation of response sensitivity with respect to change of uncertain parameters (seismic moments in the fault elements). The response surface method also enabled efficient search of the critical distribution of seismic moments in the fault elements.

## REFERENCES


# AUTHOR CONTRIBUTIONS

KM formulated the problem, conducted the computation, and wrote the paper. KK conducted the computation and discussed the results. IT supervised the research and wrote the paper.

## FUNDING

We used the GMS software as the FDM solver and the observation wave of K-NET and KiK-net which were provided by the National Research Institute for Earth Science and Disaster Prevention (NIED). In addition, **Figure 12A** was generated using the Generic Mapping Tools (Wessel and Smith, 1998). Part of the present work is supported by KAKENHI of Japan Society for the Promotion of Science (No. 17K18922, 18H01584). This support is greatly appreciated.

Bouchon, M. (1981). A simple method to calculate Green's functions for elastic layered media. Bull. Seismol. Soc. Am. 71, 959–971.


Earthquake Research Promotion. Available online at: https://www.jishin.go.jp/ main/chousa/17\_yosokuchizu/recipe.pdf


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Makita, Kondo and Takewaki. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# APPENDIX

# Verification of the Present Computational Method

To investigate the accuracy and reliability of the GMS (software), the comparison with the benchmark test (Yoshimura et al., 2011) has been conducted. **Figure A1A** shows the fault plane and recording points in the benchmark test. The fault length is 8 km and the fault width is 4 km. The other parameters are strike = 90◦ , dip = 90◦ , rake = 180◦ . The base point is located at (0,0,2)(Km). The seismic size is M<sup>0</sup> = 1.0 × 10<sup>18</sup> (Nm). The fault rupture is assumed to propagate concentrically from the source H(Hx, Hy, Hz) = (0, 1, 4)(km) with the propagation speed V<sup>r</sup> = 3000(m/s). The number of division of the fault plane is 80 × 40 (3,200). The triangle function was used as the source time function.

The three-dimensional difference grids are set as 30km × 30km × 17km (−15 km ≤ X ≤ 15 km, −15 km ≤ Y ≤ 15 km, 0 km ≤ Z ≤ 17km) . **Figure A1A** shows 1/4 of the total model. The material properties of soil layer and source layer are shown in the reference. The effective frequency is 0–1 Hz. Since inhomogeneous grids are used in GMS, the grid interval in the soil layer, and the source layer are 50 m and 150 m, respectively. The absorbing zones are set at the side and bottom of the object region to damp the reflected wave as shown in **Figure A1A** (60grids). The time duration is 15 (s) and the time increment is 0.005 (s). The number of time steps is 3,000.

**Figure A1B** presents the comparison (the observation point +010) between the result due to the GMS [designated by "Makita (GMS)"] and the benchmark test. Although the amplitude of the present GMS is slightly smaller than the benchmark test result, the overall velocity wave exhibits a similar property. This indicates the validity of the used software GMS.

# Optimal Viscous Damper Placement for Elastic-Plastic MDOF Structures Under Critical Double Impulse

#### Hiroki Akehashi and Izuru Takewaki\*

*Department of Architecture and Architectural Engineering, Graduate School of Engineering, Kyoto University, Kyoto, Japan*

A new method for optimal viscous damper placement is proposed for elastic-plastic multi-degree-of-freedom (MDOF) structures subjected to the critical double impulse as a representative of near-fault ground motions. The double impulse consists of two impulses with different directions and the critical interval between those two impulses is characterized by the criterion on the maximum input energy that is expected to lead to the maximum deformation. The critical timing of the second impulse is the timing at which the sum of the restoring force and the damping force in the first story attains zero. The objective function and constraint in terms of the maximum interstory drift along the building height or the sum of the maximum interstory drifts in all stories are selected and the corresponding optimization algorithm based on time-history response analysis and sensitivity analysis is investigated. Since the objective function in terms of the sum of the maximum interstory drifts in all stories is superior to the objective function in terms of the maximum interstory drift along the building height, it is employed in this paper. A new concept of double impulse pushover (DIP) is proposed for determining the input velocity level of the critical double impulse. It is demonstrated that the combination of two algorithms, one for effective reduction of the overflowed maximum interstory drift via the concentrated allocation of dampers and the other for effective allocation of dampers via the use of stable objective function, is effective for finding a stable optimal damper placement.

Edited by:

*Nikos D. Lagaros, National Technical University of Athens, Greece*

#### Reviewed by:

*Nikos Pnevmatikos, University of West Attica, Greece Baki Ozturk, Hacettepe University, Turkey*

> \*Correspondence: *Izuru Takewaki takewaki@archi.kyoto-u.ac.jp*

#### Specialty section:

*This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment*

Received: *18 January 2019* Accepted: *12 February 2019* Published: *07 March 2019*

#### Citation:

*Akehashi H and Takewaki I (2019) Optimal Viscous Damper Placement for Elastic-Plastic MDOF Structures Under Critical Double Impulse. Front. Built Environ. 5:20. doi: 10.3389/fbuil.2019.00020* Keywords: earthquake response, critical excitation, critical response, elastic-plastic response, viscous damping, resonance, double impulse, multi-degree-of-freedom model

## INTRODUCTION

Viscous dampers such as oil dampers are passive dampers effective for broader amplitude ranges under the condition that the stiffness is not added. This property is advantageous in most structural design of framed building structures because the change of stiffness in building structures usually leads to the variation of design loads and the necessity of change of member size. Various theories of optimal damper have been proposed (see Takewaki, 2009; Domenico et al., 2019).

In an early stage, most research is limited to elastic problems. Zhang and Soong (1992) proposed a simple algorithm to insert dampers sequentially into the location exhibiting the maximum response. Tsuji and Nakamura (1996) presented an optimality condition-based sequential algorithm for optimal damper placement. Takewaki (1997) proposed an incremental inverse problem approach to the investigation of optimal damper placement by using the transfer function as an objective performance. Takewaki et al. (1999) introduced a sensitivity-based approach in the field of optimal damper placement. Garcia (2001) extended the approach by Zhang and Soong (1992) and compared the optimization performances of several algorithms. Singh and Moreschi (2001) investigated optimal design problems using the optimality conditions and the nonlinear programming. Uetani et al. (2003) presented a practical and general damper optimization method for framed structures based on the mathematical programing. Lavan and Levy (2006) investigated a methodology for the optimal design of added viscous damping for an ensemble of realistic ground motion records with a constraint on the maximum drift. They transformed the original problem into some equivalent problems. Silvestri and Trombetti (2007) proposed physical and numerical approaches for the optimal insertion of seismic viscous dampers in shear-type structures and compared the optimization performances of several algorithms. Aydin et al. (2007) investigated the optimal damper distribution for seismic rehabilitation of planar building structures. Whittle et al. (2012) compared some viscous damper placement methods for improving seismic building design.

As far as non-linear dampers or non-linear structures with dampers are concerned, a limited number of researches has been proposed. Lavan and Levy (2005) investigated a problem of optimal design of supplemental viscous dampers for irregular shear-frames in the presence of yielding. Attard (2007) investigated a problem of optimal viscous damping for controlling interstory displacements in highly non-linear steel buildings. Lavan et al. (2008) proposed a non-iterative optimization procedure for seismic weakening and damping of inelastic structures. Adachi et al. (2013) proposed a practical theory of optimal relief-force distribution for oil dampers by setting the maximum interstory drift and the maximum buildingtop acceleration as the objective performances. Murakami et al. (2013) treated a problem of simultaneous optimal damper placement using oil, hysteretic and inertial mass dampers and proposed a sensitivity-based algorithm. The cost ratio among oil, hysteretic and inertial mass dampers may be a key point in such problem. In addition, the irregular unstable sensitivity for hysteretic dampers and recorded earthquake ground motions seems to be a difficult but challenging issue. Pollini et al. (2017) dealt with a problem of optimal placement of nonlinear viscous dampers by using the adjoint sensitivity analysis method. Shiomi et al. (2018) investigated a problem of optimal hysteretic damper placement for elastic-plastic multi-degreeof-freedom (MDOF) shear building models under the double impulse as a representative of near-fault ground motions and proposed a sensitivity-based method. Their approach set the shear building model with uniform damper distribution as an initial design and reduces the unnecessary dampers based on the sensitivity information.

However, there is no method that enables an efficient analysis of optimal viscous damper placement for MDOF building structures experiencing rather large plastic deformation.

In this paper, a new method for optimal viscous damper placement is proposed for elastic-plastic MDOF shear building structures subjected to the critical double impulse as a representative of near-fault ground motions. The critical interval between two impulses of the double impulse is characterized by the criterion on the maximum input energy. The critical timing of the second impulse is proved to be the timing at which the sum of the restoring force and the damping force in the first story attains zero. The objective function and constraint in terms of the maximum interstory drift along the building height or the sum of the maximum interstory drifts in all stories are selected and the corresponding optimization algorithm based on time-history response analysis and sensitivity analysis is investigated. Since the objective function in terms of the sum of the maximum interstory drifts in all stories is superior to the objective function in terms of the maximum interstory drift along the building height, it is employed in this paper. A new concept called the double impulse pushover (DIP), an extended version of incremental dynamic analysis (IDA) by Vamvatsikos and Cornell (2001), is proposed for determining the input velocity level of the critical double impulse. It is demonstrated that the combination of two algorithms, one for effective reduction of the overflowed maximum interstory drift via the concentrated allocation of dampers and the other for effective allocation of dampers via the use of stable objective function, is effective for finding a stable optimal damper placement.

Although the proposed two algorithms are useful, the limitation is that both algorithms deal with the inter-story drift only and do not account for the inter-story velocity which is also critical for damping device placement. For further development, the related works (Hatzigeorgiou and Pnevmatikos, 2014; Papagiannopoulos et al., 2018) discussing this limitation should be investigated.

## DOUBLE IMPULSE AND ITS CRITICAL INPUT TIMING

It is well-known that earthquake ground motions are uncertain in occurrence and properties (see Abrahamson et al., 1998). On the other hand, it is also well-recognized that near-fault ground motions possess peculiar characteristics, e.g., pulse-type waves (Bertero et al., 1978; Mavroeidis and Papageorgiou, 2003; Ozturk, 2003; Kalkan and Kunnath, 2006). To model such peculiar characteristics of near-fault ground motions, Kojima and Takewaki (2015) introduced the double impulse as a representative of the pulse-type main portion of near-fault ground motions. The acceleration of the double impulse with the time interval t<sup>0</sup> of two impulses can be expressed by

$$
\ddot{u}\_{\xi}(t) = V\delta(t) - V\delta(t - t\_0), \tag{1}
$$

where V is the input velocity amplitude and δ(t) is the Dirac delta function. **Figure 1A** shows the acceleration, velocity and displacement of the double impulse together with those of the corresponding one-cycle sine wave of acceleration. While several investigations using the double impulse have been made for single-degree-of-freedom (SDOF) models (Kojima and Takewaki, 2015, 2016; Kojima et al., 2017; Akehashi et al., 2018a,b), researches on MDOF models are few (Taniguchi et al.,

2016; Saotome et al., 2018; Shiomi et al., 2018). This is due to the fact that the simple energy balance law for the simple derivation of the maximum response is difficult to apply for MDOF models because of the phase lag. Taniguchi et al. (2016) treated an undamped 2DOF elastic-plastic model under the double impulse and found by using the criterion of the maximum input energy that the timing of the second impulse at the zero restoring-force state becomes the critical timing. However, this critical timing cannot be used for damped models.

Consider an SDOF elastic perfectly-plastic model with viscous damping and a MDOF elastic perfectly-plastic model with viscous damping as shown in **Figure 1B**. Let us consider first the SDOF model subjected to the double impulse. The critical input timing of the second impulse can be defined as the timing which maximizes the input energy by the second impulse. To demonstrate this fact, consider the following equation of motion.

$$m\ddot{u} + c\dot{u} + f(k, u, d\_{\mathcal{Y}}, d\_r) = 0,\tag{2}$$

where m, c, f , k, u, dy, d<sup>r</sup> denote the mass, damping coefficient, restoring force, initial stiffness, displacement, yield displacement and residual displacement. The super dot indicates the differentiation with respect to time. Assume that this SDOF model is subjected to the double impulse and the free vibration occurs. Since the displacement does not change instantaneously at the second impulse, the strain energy does not change at the second impulse. Therefore, the input energy by the second impulse can be expressed by

$$E = \frac{1}{2}m(\dot{u} + V)^2 - \frac{1}{2}m\dot{u}^2 = (m\dot{u})V + \frac{1}{2}mV^2,\tag{3}$$

where u˙ is the velocity just before the input of the second impulse. Equation (3) indicates that, when the velocity attains the maximum, the input energy by the second impulse yields the maximum. The condition that u˙ attains the extremum is u¨ = 0. Substitution of u¨ = 0 into Equation (2) leads to cu˙ + f = 0. This means that, when the sum of the restoring force and the damping force becomes zero (i.e., the velocity becomes the maximum), the input energy by the second impulse becomes the maximum.

Consider next an N-story MDOF shear building model of mass m<sup>i</sup> in the i-th story. Let u<sup>i</sup> denote the horizontal displacement of mass m<sup>i</sup> . As in the SDOF model, since the displacements do not change instantaneously at the action of the second impulse, the strain energy does not change at the second impulse. Therefore, the input energy by the second impulse can be expressed by

$$E = \sum\_{i=1}^{N} \frac{1}{2} m\_i (\dot{u}\_i + V)^2 - \sum\_{i=1}^{N} \frac{1}{2} m\_i \dot{u}\_i^2 = V \sum\_{i=1}^{N} m\_i \dot{u}\_i$$

$$\begin{split} + \sum\_{i=1}^{N} \frac{1}{2} m\_i V^2 \text{ (4)} \end{split}$$

Equation (4) means that, when P<sup>N</sup> i=1 (miu˙i) attains the maximum, the input energy by the second impulse becomes the maximum. The condition that P<sup>N</sup> i=1 (miu˙i) attains the extremum is P<sup>N</sup> i=1 (miu¨i) = 0. Since P<sup>N</sup> i=1 (miu¨i) is equal to F<sup>1</sup> = c1u˙ <sup>1</sup> + f<sup>1</sup>

owing to the dynamic equilibrium, i.e., the sum of the damping force and the restoring force in the first story, the extremum condition becomes F<sup>1</sup> = c1u˙ <sup>1</sup> + f<sup>1</sup> = 0. This critical condition is very simple and can be used in the time-history response analysis in a simple manner.

**Figure 2** shows an example of the time history of F<sup>1</sup> = c1u˙ <sup>1</sup> + f<sup>1</sup> for the model subjected to the first single impulse and the variation of the input energy E by the second impulse with respect to t0.

It seems important to investigate the correspondence of the double impulse with recorded ground motions, the Rinaldi Station FN component (Northridge 1994) was used (see **Appendix**).

# PROBLEM OF OPTIMAL DAMPER PLACEMENT AND ALGORITHM OF SOLUTION

The direct problem of optimal damper placement may be stated as follows: To minimize the maximum interstory drift (or the sum of the maximum interstory drifts along height) under the condition on the specified total quantity of passive dampers. Another problem may be described as: To minimize the total quantity of passive dampers under the constraint on the maximum interstory drift. These two problems may be proved to be almost equivalent. It seems possible to deal with these problems by using the sensitivity-based approach that includes the time-history response analysis for the double impulse and the finite difference method. However, it was found that the direct application of this approach to the above mentioned problems leads to unstable results. To overcome this difficulty, a mixedtype approach including the following two problems (Problem 1 and 2) is proposed in this paper. The problems treated here will be explained next.

Consider first the following problem.

#### Problem 1

$$\begin{aligned} &\min \, \mathbf{c}^{\mathrm{T}}\_{\mathrm{add}} \cdot \mathbf{1} \\ &\text{subject to: } d\_{\mathrm{max},i} \le d\_{\mathrm{target},i} \text{, for all } i \end{aligned}$$

In Problems 1, **c**add is the damping coefficient vector for added dampers in all stories and **1** is the vector including one in all components. The superscript T indicates the transpose. dmax,<sup>i</sup> is the maximum interstory drift in the i-th story under the critical double impulse and dtarget,<sup>i</sup> is the target value of dmax,<sup>i</sup> . The solution algorithm for this problem may be described as follows.

## Algorithm 1


The flow of Algorithm 1 is shown in **Figure 3A**. Algorithm 1 is intended to implement the effective reduction of the overflowed maximum interstory drift (the stories i in Step 2) via the concentrated allocation of dampers. It was found after some attempts that Algorithm 1 sometimes encounters difficulties, i.e., inability to deal with the problem for a prescribed damper quantity.

Consider secondly the following problem.

#### Problem 2

$$\begin{aligned} \min f &= \sum\_{i=1}^{N} d\_{\text{max},i} \\ \text{subject to: } \mathbf{c}^{\text{T}}\_{\text{add}} \cdot \mathbf{l} &= \text{const.} \end{aligned} $$

The solution algorithm for this problem may be described as follows.

## Algorithm 2


Step 3 Find the story in which the largest reduction of f occurs. For that story i, update the damping coefficient as c<sup>i</sup> → c<sup>i</sup> + 1c. Put j→ j+1. If 1c· j = **c** T add · **1** is satisfied, then finalize the process. If 1c· j < **c** T add · **1** is satisfied, return to Step 2.

The flow of Algorithm 2 is shown in **Figure 3B**. It was found after some attempts that Algorithm 2 sometimes encounters difficulties in efficiency depending on models. The detailed explanation will be shown in numerical examples afterwards.

FIGURE 9 | Distribution of added damping coefficients, P*d*max,*i*/*d<sup>y</sup>* with respect to step number and the distribution of *d*max,*i*/*d<sup>y</sup>* under the critical double impulse for Model 1 (Problem 2): (A) Elastic limit, (B) *V* = 0.84 [m/s], (C) *V* = 1.23 [m/s].

Consider thirdly the following mixed problem.

# Problem 3: Mixed Problem of Problem 1 & Problem 2

Problem 1 is solved at first until some stage and then Problem 2 is solved.

The solution algorithm for this problem may be described as follows.

## Algorithm 3: Combination of Algorithm 1 & Algorithm 2

First of all, apply Algorithm 1. Set dtarget,<sup>i</sup> = dy,<sup>i</sup> for all i and obtain the model in which the largest interstory drift attains the elastic limit. Adopt this model as another initial model and repeat Algorithm 2 in **c** T add · **1**/1c steps.

Algorithm 3 is superior to Algorithm 1 and 2 because the combination of two algorithms, one for effective reduction of the overflowed maximum interstory drift via the concentrated allocation of dampers and the other for effective allocation of dampers via the use of stable objective function, is effective for finding a stable optimal damper placement.

# MODELS FOR NUMERICAL EXAMPLES

Consider three models of 12 stories with different story stiffness distributions. Model 1 has a uniform distribution of story stiffnesses. Model 2 has a straight-line lowest eigenmode. Model 3 has stepped distribution of story stiffnesses (upper four stories, middle four stories and lower four stories have uniform stiffness distributions with different values/ the ratios among them are 1: 1.5: 2). The undamped fundamental natural period of these three models is 1.2[s] and the structural damping ratio is 0.01 (stiffness proportional type). All the floor masses have the same value. The common story height is 4[m] and the common yield interstory drift ratio is 1/150. The story shear-interstory drift relation obeys the elastic perfectly-plastic rule.

**Figure 4** shows the eigenmodes multiplied by the participation factor (participation vectors) and the natural periods for three models.

# DYNAMIC PUSHOVER ANALYSIS UNDER AMPLIFIED CRITICAL DOUBLE IMPULSE

To determine the input velocity level of the critical double impulse, the incremental dynamic analysis (IDA) procedure

(Vamvatsikos and Cornell, 2001) is applied to the critical double impulse. It should be reminded that only the critical double impulse is treated here, i.e., the interval of two impulses of the double impulse is varied depending on the input velocity level (also depending on the maximum interstory drift). We call this procedure "Double impulse pushover (DIP)". DIP provides the relation between the maximum interstory drift and the input velocity level of the critical double impulse. While the conventional IDA includes multiple recorded ground motions for taking into account the uncertainty in predominant periods of ground motions, DIP adopts the critical input and enables an efficient analysis of the relation between the maximum response and the input level.

**Figure 5** shows the maximum interstory drift distributions by DIP. The velocity level is increased from V = 0.2 [m/s] to V = 1.6 [m/s] by 0.2 [m/s].

Since the input velocity level of the critical double impulse influences greatly the maximum interstory drift and the optimal damper placement, its determination appears very important. The determination process of the input velocity level of the critical double impulse is explained in the next.


Model 3 is treated as an example for determining the input level. The double of the yield interstory drift is taken as the specified target value of the maximum interstory drift. Then it is found that over 0.6 [m/s] is necessary. Although the input velocity level should be chosen for each structural model, V = 0.84 [m/s] and 1.23 [m/s] are employed in the following section.

#### NUMERICAL EXAMPLES

#### Examples for Problem 1 Using Algorithm 1

Consider first some examples for Problem 1. The amplitudes of the critical double impulses have been determined from the results for DIP explained in Section Dynamic Pushover Analysis Under Amplified Critical Double Impulse.

**Figure 6** shows the distribution of added damping coefficients, max(dmax,i/dy) with respect to step number and the distribution of dmax,i/d<sup>y</sup> under the critical double impulses with V = 0.84 [m/s] and V = 1.23 [m/s] for Model 1. The condition dtarget,<sup>i</sup> = d<sup>y</sup> (for all i) is employed here.

**Figure 7** presents the similar figures for Model 2 and **Figure 8** illustrates the similar figures for Model 3.

It can be observed from **Figures 6**–**8** that, as the input velocity level is increased, the ratios of damping coefficients of added dampers along height change and their allocation becomes smooth in the wide height range. This is because, as the input level becomes larger, the number of stories experiencing plastic deformation becomes larger. Secondly, the maximum interstory drift distribution in the final model becomes almost uniform by Algorithm 1. Furthermore, it is understood that the increase of max(dmax,i/dy) in the damper allocation process is allowed.

For Model 3, it can be observed that the dampers are not allocated in the 4th and 8th stories for the input level V = 0.84 [m/s], but those are allocated for the input level V = 1.23 [m/s]. However, the maximum interstory drifts in those stories are smaller than the elastic limit in the initial stage. This means that Algorithm 1 acts first so that the dampers are allocated to the 1, 5, 9 th stories experiencing large plastic deformation. This process helps the energy required for deformation distribute to the neighboring stories. As a result, among the stories neighboring to the 5, 9th stories, the deformations in the 6, 10th stories with relatively small stiffness becomes larger. Since the 4, 8th stories with relatively large stiffness go into the plastic range for the input level of V = 1.23 [m/s], the dampers are allocated so as to strengthen the model. A similar observation may be possible also for Model 1 and 2.

It may be concluded that Algorithm 1 is apt to allocate added dampers in a concentrated manner to the stories where the plastic deformation develops fast, then to allocate additional ones to rather weak stories after the strengthening is completed.

#### Examples for Problem 2 Using Algorithm 2

Consider next some examples for Problem 2. The parameter specification 1001c = **c** T add · **1** is given here.

**Figure 9** shows the distribution of added damping coefficients, Pdmax,i/d<sup>y</sup> with respect to step number and the distribution of dmax,i/d<sup>y</sup> under the critical double impulses with V = 0.84 [m/s] and V = 1.23 [m/s] for Model 1. The distributions for the elastic limit are also shown for reference (max(dmax,i/dy) = 1).

**Figure 10** presents the similar figures for Model 2 and **Figure 11** illustrates the similar figures for Model 3.

It can be observed that, when the elastic limit is employed for determining the input velocity level, the dampers are allocated so that the maximum interstory drifts become almost uniform for all models (Model 1–3). For Model 1, the dampers are allocated so that the maximum interstory drifts become almost uniform regardless of the input velocity level. On the other hand, for Models 2 and 3, the damper distributions are different depending on the input velocity level. Furthermore, for Models 2 and 3, the maximum interstory drifts in specific stories do not change between the initial model and the final model. This means that, since Algorithm 2 is aimed at finding the optimal damper allocation by using the steepest direction, it does not provide better damper allocation from the viewpoint of uniform reduction of the maximum interstory drifts in all stories.

## Examples for Mixed Problem (Problem 3) of Problem 1 and 2 Using Algorithm 3

Examples for Mixed Problem (Problem 3) of Problem 1 and 2 using Algorithm 3 are presented here. The parameter specification by 1001c = **c** T add · **1** is used for V = 0.84 [m/s] and the parameter specification by 2501c = **c** T add · **1** is used for V = 1.23 [m/s].

**Figure 12** shows the distribution of added damping coefficients, Pdmax,i/d<sup>y</sup> with respect to step number and the distribution of dmax,i/d<sup>y</sup> under the critical double impulses with V = 0.84 [m/s] and V = 1.23 [m/s] for Model 1.

**Figure 13** presents the similar figures for Model 2 and **Figure 14** illustrates the similar figures for Model 3.

It can be observed that, while the analysis in the elastic range is not easy by Algorithm 1 for Problem 1 owing to the inability to set the total damper quantity and the maximum interstory drift distribution of the final model obtained by Algorithm 2 for Problem 2 is unstable (not uniform) depending on the model, a stable damper allocation is possible by Algorithm 3 for Problem 3 regardless of the input velocity level of the critical double impulse and the final maximum interstory drift distributions are apt to become uniform. It may be said that Algorithm 3 enables the procedure to guarantee the minimum performance by using a small amount of added dampers and to reduce the structural response globally by using the additional amount of added dampers.

## REFERENCES


## CONCLUSIONS

A new method for optimal viscous damper placement has been proposed for elastic-plastic MDOF structures subjected to the critical double impulse as a representative of near-fault ground motions. The main conclusions can be summarized as follows.


# AUTHOR CONTRIBUTIONS

HA formulated the problem, conducted the computation, and wrote the manuscript. IT supervised the research and wrote the manuscript.

# FUNDING

Part of the present work is supported by the Grant-in-Aid for Scientific Research (KAKENHI) of Japan Society for the Promotion of Science (No.17K18922, 18H01584) and Sumitomo Rubber Industries, Co. This support is greatly appreciated.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbuil. 2019.00020/full#supplementary-material


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Akehashi and Takewaki. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Influence of the Brace Configurations on the Seismic Performance of Steel Concentrically Braced Frames

T. Y. Yang1,2 \*, Hediyeh Sheikh<sup>2</sup> and Lisa Tobber <sup>2</sup>

*1 International Joint Research Laboratory of Earthquake Engineering, Tongji University, Shanghai, China, <sup>2</sup> Department of Civil Engineering, University of British Columbia, Vancouver, BC, Canada*

#### Edited by:

*Izuru Takewaki, Kyoto University, Japan*

#### Reviewed by:

*Yutaka Nakamura, Shimane University, Japan Constantinos Repapis, University of West Attica, Greece Peng Pan, Tsinghua University, China Manolis S. Georgioudakis, National Technical University of Athens, Greece*

> \*Correspondence: *T. Y. Yang*

*yang@civil.ubc.ca*

#### Specialty section:

*This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment*

Received: *13 December 2018* Accepted: *18 February 2019* Published: *14 March 2019*

#### Citation:

*Yang TY, Sheikh H and Tobber L (2019) Influence of the Brace Configurations on the Seismic Performance of Steel Concentrically Braced Frames. Front. Built Environ. 5:27. doi: 10.3389/fbuil.2019.00027* Concentrically braced frame (CBF) is an effective and prevalent seismic force resisting system which is commonly used in low-rise buildings. This type of structural system utilizes steel braces to provide the stiffness and strength needed to dissipate earthquake energy. Several bracing configurations have been proposed in different building codes worldwide. These codes provide detailed design requirements for the structural members and connections, but no guidance is provided in selecting the best bracing configuration for the design. In this study, the impact of the bracing configuration on the seismic response of a five-story prototype office building located in Vancouver, Canada, is systematically examined. Five different bracing configurations were designed according to the National Building Code of Canada and CSA S16 standard. Detailed structural responses, initial costs, and life cycle costs of the prototype building with five different bracing configurations were systematically compared. The results show that the different bracing configurations play an important role in sizing the structural members, which impacts the initial material usage and the overall life cycle cost of the building.

Keywords: seismic performance evaluation, concentrically braced frame, time history analysis, life cycle cost, finite element model

## INTRODUCTION

Steel concentrically braced frame (CBF) is a seismic force resisting system (SFRS) commonly used in seismic zones around the world. This type of SFRS is effective in providing the stiffness and strength needed to resist earthquake forces. In the Canadian steel building code, CSA S16-14 (CSA, 2014), there are two types of steel CBFs: (1) moderate ductile concentrically braced frame (MD-CBF) and (2) limited ductile concentrically braced frame (LD-CBF). The MD-CBF is targeted to be used in high seismic zones, where the SFRS is designed to have enhanced ductility through yielding of the steel braces, while the beams and columns are capacity designed to resist the maximum load produced by the braces. On the other hand, the LD-CBF is targeted to be used for locations with less earthquake shaking, where the ductility requirement of the braces can be relaxed.

Multiple MD-CBF configurations have been pre-qualified by the CSA S16-14 (CSA, 2014). The configurations include: inverted V-braced (I-VBF), V-braced (VBF), X-braced (XBF), and Multistory-X-braced (M-XBF).

According to CSA-S16-14, the probable tensile and compressive resistance of bracing members shall be taken as Equations (1, 2), respectively.

$$\mathbf{T\_u = |A\_{\mathfrak{K}} R\_{\mathfrak{Y}} F\_{\mathbf{Y}}} \tag{1}$$

$$\mathbf{C\_u = \ min\left(A\_{\mathfrak{K}} \mathbf{R\_{\bar{Y}}} \mathbf{F\_{\bar{Y}}}, \ 1.2 \mathbf{C\_r}/\phi\right)}\tag{2}$$

Where A<sup>g</sup> = gross cross sectional area of the brace; R<sup>y</sup> = the expected strength factor; F<sup>y</sup> = specified yield stress; φ = resistance reduction factor and Cr is computed using RyFy.

The probable post-buckling compressive resistance of bracing members shall be taken as Equation (3).

$$\mathbf{C}'\_{\rm u} = \min \left( \mathbf{0}.2 \, \mathbf{A}\_{\rm \S} \mathbf{R}\_{\rm \Y} \mathbf{F}\_{\rm y}, \, \mathbf{C}\_{\rm r} / \phi \right) \tag{3}$$

In the V- or inverted V- bracing configurations, the beam must be designed to be continuous between columns and capable of resisting the maximum unbalanced vertical and horizontal loads when the braces below (inverted V) or above (V) the beam has buckled and yielded. This makes the beams in the V- or inverted V- bracing configurations of the MD-CBF very large. To mitigate the unbalanced force demand in the beam, the suspended zipper braced frame (SZBF) was proposed by Yang et al. (2009b). SZBF uses zipper columns placed between the braces to transfer the unbalanced vertical forces to higher stories. To prevent the structure from losing total stiffness, the top story braces of the SZBF are capacity designed to resist the combined unbalanced vertical loads when all the braces yielded and buckled. Similarly, the beams in the M-XBF system must be capacity designed for the yielding and buckling of the braces. Due to the symmetrical geometry of the M-XBF configuration, the unbalanced forces from the floor above and below cancel each other out (if the braces are of the same size). This results in smaller beams compared to the V- or inverted V- bracing configurations. Similarly, in the X-bracing configuration, there are no unbalanced vertical forces in the beams, but the columns are designed to resist large axial forces when the braces have yielded and buckled.

The building code provides detailed design requirements for each of the MD-CBF bracing configurations. However, there has been little discussion about the differences in seismic performance. Although there have been some previous studies comparing different bracing configurations (Mahmoudi and Zaree, 2010; Patil and Sangle, 2015; Ozcelik et al., 2016), the present study is the only one that shows a detailed cost comparison on a Canadian designed building.

In this study, a five-story office prototype building located in Vancouver, British Columbia, Canada was used to compare the seismic performance of the prototype building using each of the I-VBF, VBF, XBF, M-XBF, and SZBF bracing configurations. The prototype building was designed according to the National Building Code of Canada (NBCC, 2015) and the Canadian Institute of Steel Construction design standard S16-14 (CSA, 2014). Detailed finite element models of the prototype building with each of the bracing configurations were developed using OpenSees (McKenna et al., 2009). The finite element models were calibrated against available experimental data. The calibrated finite element models were then subjected to ground motions selected and scaled to three earthquake shaking intensities at the prototype site. Detailed initial construction costs, structural

**Abbreviations:** CBF, Concentrically Braced Frame; SFRS, Seismic Force Resisting System; MD-CBF, Moderate Ductile Concentrically Braced Frame; LD-CBF, Limited Ductile Concentrically Braced Frame; I-VBF, Inverted V-braced Frame; VBF, V-braced Frame; XBF, X-braced Frame; M-XBF, Multistory-X-braced Frame; SZBF, Suspended Zipper Braced Frame; ISD, Peak Inter Story Drift ratio; FA, peak Floor Acceleration; Std, Standard Deviations; PBEE, Performance-Based Earthquake Engineering; EDP, Engineering Demand Parameter; DS, Damage State; PG, Performance Group; Min\_Qty, Minimum Quantities; Max\_Qty, Maximum Quantities; MAL, Mean Annualized Loss.

response, and life cycle costs of the prototype building using each of the pre-qualified MD-CBF configurations were examined in this study.

# DESCRIPTION OF THE PROTOTYPE BUILDING

The prototype structure is a five-story office building located in Vancouver, British Columbia, Canada, which was adopted from Yang and Murphy (2015). The prototype building was designed according to the National Building Code of Canada (NBCC, 2015) and the CSA-S16-14 (CSA, 2014). The plan dimensions of the building are 45 m by 63 m. The story height is 4.25 m for the first floor and 3.65 m for other floors. **Figure 1** shows the global view and plan view of the SFRS. Five different bracing configurations: (1) Inverted V-braced (I-VBF); (2) Vbraced (VBF); (3) Suspended zipper column braced (SZBF); (4) Multistory-X-braced (M-XBF); (5) X-braced (XBF), were included in this study. **Figure 2** shows the member sizes of the SFRS. **Table 1** shows the summary of the gravity frame elements used for all bracing configurations.

## DESCRIPTION OF THE FINITE ELEMENT MODELS

The numerical models of the prototype building were generated using OpenSees (McKenna et al., 2009). Due to the symmetrical nature of the prototype building in the two principal axes, only half of the building was modeled. Because the SFRS was designed to be decoupled in the two principal axes, only the response in the East-West direction is presented in this paper. Masses were lumped at the nodes of the model at each floor based on the tributary area. An average value of 633,000 Kg was assigned at each floor. Rayleigh damping of 2.5%, based on the mass and tangent stiffness proportional damping, was assigned in the first and third modes. The first three vibrational periods of each MD-CBFs are shown in **Table 2**.

The steel braces were modeled following the approach developed by Yang et al. (2009a), where the braces were modeled using two flexibility-formulation non-linear beamcolumn elements with fiber sections. Uniaxial Menegotto-Pinto steel material (Steel02) in OpenSees was used to model the kinematic and isotropic hardening of the steel material. Rotational springs with rotational stiffness of EI/5Lb, where E is the modulus of elasticity, I is the moment of inertia about the plane of bending, and Lb is the total length of the brace, were added to the ends of the braces using zero-length elements to simulate the gusset plate response (Yang et al., 2009b). A geometric imperfection of 0.1% of the brace length was used according to the procedure presented in Uriz (2008). The corotational geometry transformation proposed by De Souza (2000) was used to simulate the second-order geometry effects (P-1 effects). **Figure 3A** shows the analytical model for the steel braces. **Figure 3B** and **Figure 3C** show the comparison of the brace hysteresis from the analytical simulation against the experimental data presented by Black et al. (1980). The results TABLE 1 | Size of the gravity frame elements.


TABLE 2 | Structural periods of systems.


show that the proposed numerical model can accurately simulate the non-linear behavior of the steel braces.

The beams and columns of the CBFs were modeled using non-linear fiber-based beam-column elements in OpenSees. The column base connections were modeled using zerolength element according to stiffness relationships proposed by Fahmy et al. (1998). Beam to column connections for shear tab connections were modeled with the rotational properties suggested by Astaneh-Asl (2005). **Figure 4** shows the finite element models developed for this study.

## SELECTION AND SCALING OF GROUND MOTIONS

To simulate the seismic response of the prototype building, the 5% damped spectra for the site was obtained from Natural Resources Canada (NRC) (2016) at three hazard levels: (1) 2% probability of exceedance in 50 years (2/50); (2) 10% probability of exceedance in 50 years (10/50); and (3) 40% probability of exceedance in 50 years (40/50). Ground motions were selected from the PEER NGA database [Pacific Earthquake Engineering Research Center (PEER), 2011]. The ground motions with magnitude (Mw) between 5.5 and 8, soil shear wave velocity (Vs30) between 360 and 760 m/s, and the limited distance to the fault (0–100 km) were selected. The ground motions were amplitude scaled such that the mean spectrum for the scaled ground motions does not fall below the target spectrum by 10% over the period range from 0.15 Tmin to maximum of 2.0 Tmax and 1.5 s (where Tmin and Tmax are the shortest and longest fundamental periods of the structures, respectively). To avoid over scaling, the scale factors were limited between 0.5 and 4 as

noted in the commentary of NBCC (2015). **Table 3** summarizes the selected ground motions used in this study. **Figure 5** shows the scaled response spectra and the target spectrum for the three hazard levels considered.

# RESPONSE QUANTIFICATION

A total of 57 non-linear dynamic analyses were carried out to examine the response of the prototype building using each of the bracing configurations. **Table 4** shows the summary of the structural responses. **Figure 6** shows the comparison of the peak inter story drift ratio (ISD) and peak floor acceleration (FA) and the corresponding standard deviations (Std).

At the 2/50 hazard level, all configurations have the highest ISD at the first floor. Among the bracing configurations, the SZBF, VBF and the I-VBF has the highest 1st floor ISD, while the M-XBF and XBF has the lowest 1st floor ISD. All bracing configurations have a significant reduction in ISD at the 2nd floor. The I-VBF, VBF, M-XBF, and XBF have some increase of ISD at the 3rd and 4th floor. This is likely caused by the higher mode effect. All systems show a reduction of ISD at the 5th floor. XBF has the highest ISD at the 5th floor, while the SZBF has the lowest 5th floor ISD. The variation in the ISD shows very similar trends as the median response, where the 1st floor has the largest variation. Peak floor accelerations tended to increase with height. XBF has slightly higher FA at all floors, while the other configurations all have very similar FA. The variation in FA is very similar among all configurations.

At the 10/50 hazard level, the ISD for the I-VBF, VBF, and XBF configurations are quite similar among the floors. On the other hand, the SZBF has lowest ISD at the 1st floor with increased ISD at the 2nd and 3rd floors and reduced ISD at the 4th and 5th floor. M-XBF has a similar trend as the SZBF, except 4th and 5th floor have less reduction as the SZBF. The variation in ISD is

similar among the floors. In terms of FA, all configurations have a similar trend, where acceleration increases linearly at the floor number increases. The variation of FA shows a similar trend for all systems.

At the 40/50 hazard level, I-VBF and XBF configurations have similar ISD at all floors. SZBF has a similar trend as the M-XBF at the 2nd and 3rd floors but has smaller ISDs at the 4th and 5th floor. Among the different configurations, the VBF has the highest ISD for almost all floors (except the 2nd floor). The variation in ISD is very consistent among all floors. The trend for the FA is very similar for all configurations, where the FA increases as the floor number increases. The variation of FA is very similar among all floors and bracing configurations.

# PERFORMANCE ASSESSMENT

To quantify the seismic performance of the prototype building with different bracing configurations, the performance-based earthquake engineering (PBEE) method developed by Yang et al. (2009a) was used. The main steps of the PBEE method can be summarized as the following: (1) seismic hazard analysis; (2) response analysis; (3) damage analysis, and (4) loss analysis. **Figure 7A** illustrates the four analysis phases. The goal of the seismic hazard analysis is to determine the ground motions that are most representative of the seismic hazard of the prototype building site. The response analysis is used to identify the engineering demand parameters (EDPs), such as drift and acceleration that could damage the structural and nonstructural components. Once the peak structural responses have been identified, the damage analysis utilizes fragility curves obtained from past research to determine the damage state (DS) of each structural and non-structural component in the prototype building. Because many structural and non-structural components could be affected by the same EDP, and the DS could be similar, these structural and non-structural components could be grouped into performance group (PG). The last step of the PBEE method is the loss analysis, where the analysis uses DS identified in the damage analysis to quantify the building repair cost.

In this study, 26 different PGs were defined. Each PG consists of one or more building components whose performance was similarly affected by the same EDP. The DS of each PG was determined in relation to the different types of repairs that may be required. For each DS, a fragility curve was defined to present the probability that the damage of the components of a PG will be equal or greater than the specified DSs for a given EDP. **Figure 7B** shows a sample fragility curve for the braces. In this fragility curve, there are four DSs. As shown in this figure, if the peak ISD equals to 1.5%, the component has a 5, 70, 25, and 0% probability of being in DS1, DS2, DS3, and DS4, respectively. The numerical values of the fragility relations, unit repair cost data, and repair quantities for all the structural and non-structural PGs were adapted from the FEMA P-58 project (FEMA, 2012). **Table 5** provides the summary of the PGs and fragility data used in this study. The first five rows of **Table 5** present the structural PG for the MD-CBFs system, and the rest of the table presents the non-structural components. The median and beta parameters of each DS used to determine fragility curve are also reported in **Table 5**. **Figure 7** shows the tri-linear curve used to represent the unit repair cost. Min\_Qty represents the minimum quantities, Max\_cost represents the maximum cost, Max\_Qty represent the maximum quantities and Min\_cost represent the minimum cost. As shown in **Figure 7**, the unit repair cost is expected to reduce when the repair quantity is increased. **Table 6** lists the data used to quantity the unit repair cost for each PG.

# RESULT OF PERFORMANCE EVALUATION

**Table 7** shows the initial structural costs for each of the MD-CBFs configurations. The construction cost of steel

#### TABLE 3 | Summary of ground motions.



components was calculated using typical steel prices in Vancouver Canada, where a unit cost of \$5.6 and \$7 USD per kg was used for members over and below 60 kg/m, respectively. The result shows the M-XBF uses the least column and brace material, while the XBF and SZBF have the lowest beam material. Overall, the M-XBF has the minimum structural cost, while the I-VBF has the maximum structural cost.

**Table 8** shows the summary of the initial construction costs, median repair costs, mean annualized cost and life cycle costs for the different bracing configurations. The initial construction cost of the prototype building was defined as the sum of the structural and non-structural components costs. The non-structural component costs were calculated using an average of \$300 USD per square foot which includes the building contents, partition walls, mechanical systems, finishes,

FIGURE 5 | Response spectra of scaled ground motions at (A) the 2/50 hazard level, (B) the 10/50 hazard level, and (C) the 40/50 hazard level (Target spectrum from Natural Resources Canada (NRC), 2016).

TABLE 4 | Median and standard deviation (Std) of peak EDPs.


*dui and ai represent the ISD at the ith story and the FA at the ith floor, respectively.*

office furniture, computers and other contents (Yang and Murphy, 2015). The result shows the difference in structural bracing configurations makes <1% difference in the overall construction cost.

deviation in right plots at (A) 2/50 hazard level, (B) 10/50 hazard level, and (C) 40/50 hazard level.

**Table 8** also shows the median repair costs. By comparing the median repair costs of different MD-CBFs configurations at 2/50 hazard level, the I-VBF, SZBF, and M-XBF have very close median repair cost in the range of \$1.8 million USD, while the XBF and VBF have the highest median repair cost of about \$2.3 million USD. At the 10/50 hazard level, the lowest median repair costs are for M-XBF and VBF which is about \$1.2 million and the highest one is for XBF which is about \$1.3 million USD. The median repair costs at 40/50 hazard level are very close for all configurations which are within the range from \$340,000 USD to \$370,000 USD.

The mean annual rate of the repair cost exceeding a threshold value can be determined by combining the repair cost and the seismic hazard relations. The annual exceedance

probability is obtained by multiplying the slope of the hazard curve at the corresponding ground motion intensity level by the complement of the cumulative repair distribution function then integrated the resulting curves across the seismic hazard levels. This process creates a loss curve that represents the mean annual rate of the repair cost exceeding a threshold value. The area under the loss curve represents the mean cumulative total repair cost for all earthquakes over the period of 1 year. Owners and stakeholders can use this information to quantify the mean annualized loss (MAL) when comparing the performance of different bracing configurations (Yang et al., 2009a). The results of this study show that the M-XBF has the lowest MAL of \$24,800 USD, while the XBF has the highest MAL of \$32,800 USD.

Using the initial costs and mean annualized loss, the total life cycle costs for 50 years with 3.5% annual interest rate is summarized in **Table 8**. Comparing the five MD-CBF bracing configurations, the M-XBF is the most economical system and has the lowest life cycle cost. Followed by the SZBF, VBF, I-VBF, and XBF.

# SUMMARY AND CONCLUSIONS

Steel CBF is a commonly used SFRS. This structural system is efficient in providing the stiffness and strength needed to resist the seismic load. In this type of SFRS, the braces are designed to be the main energy dissipation component, where


TABLE 5 | Summary of performance groups and their corresponding fragility curve data.

*dui and ai represent the ISD at the ith story and the FA at the ith floor, respectively.*

they are expected to buckle and yield during strong earthquake shaking. To ensure the system still has a stable force-deformation response after the braces yielded or buckled, the beams and columns of the CBF are capacity designed based on the maximum expected brace forces. Multiple bracing configurations have been presented in the building codes. However, the building codes do not indicate which brace configuration is best suited for different application. In this study, a 5-story office building located in Vancouver, Canada was designed using five different CBF bracing configurations according to NBCC (2015), Natural Resources Canada (NRC) (2016) and CSA S16-14 (CSA, 2014). Detailed finite element models of the prototype building with each of the bracing configurations were developed. The finite element models were calibrated against experimental test then subjected to ground motions selected and scaled to three hazard levels. Detailed structural response, initial construction costs and life cycle costs of the prototype building with each of the bracing configurations were systematically examined. The following findings were observed in this study:

1) The I-VBF uses the most steel material, while the M-XBF requires about 20% less steel material than the I-VBF, making the M-XBF the lightest system among the different bracing configurations considered. The XBF has the second lightest structural materials, followed by the VBF and SZBF systems.



*"w" is weight per unit length, and "PLF" is lb. per LF.*

TABLE 7 | Summary of the structural elements cost [USD].


TABLE 8 | Summary of the initial, median repair, and life cycle costs [USD].


except the 4th and 5th floors which have less reduction than the SZBF. The variation in ISD is similar among the floors. In terms of FA, all configurations show a similar trend, where acceleration increases linearly as the floor height increases. The variation of FA shows a similar trend for all systems.

4) At the 40/50 hazard level, I-VBF and XBF configurations have similar ISD at all floors. SZBF has the similar trend as the M-XBF at the 2nd and 3rd floors, with a reduction in ISD at the 4th and 5th floor. Among the different configurations, the VBF has the highest ISD for almost all floors, except the 2nd floor. The variation in ISD is very consistent among all floors. The trend for the FA is similar for all configurations, where the FA increases as the floor height increases. The variation of FA is very similar among all floors and bracing configurations.

5) XBF is the most expensive system as it has the highest life cycle cost. By contrast, the M-XBF is the most economical system having the lowest life cycle cost. After that, the SZBF, VBF, and I-VBF are ranked from the lowest life cycle cost to highest.

#### DATA AVAILABILITY

All datasets generated for this study are included in the manuscript and/or the supplementary files.

## AUTHOR CONTRIBUTIONS

HS contributed development of the proposed models, performed the analysis, and wrote the draft of the manuscript. LT wrote sections of the manuscript and revised the paper. TY supervised

#### REFERENCES


the research and revised the paper. All authors read and approved the submitted manuscript version.

#### ACKNOWLEDGMENTS

We would like to acknowledge the funding provided by the National Natural Science Foundation of China (grant number: 51778486), International Joint Research Laboratory of Earthquake Engineering (ILEE), National Science Foundation China, State Key Laboratory of Disaster Reduction in Civil Engineering. Any opinions, findings, and conclusions or recommendations expressed in this paper are those of the authors.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Yang, Sheikh and Tobber. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Experiments and Fragility Analyses of Piping Systems Connected by Grooved Fit Joints With Large Deformability

Tao Wang<sup>1</sup> \*, Qingxue Shang<sup>1</sup> , Xi Chen<sup>2</sup> and Jichao Li <sup>1</sup>

*<sup>1</sup> Key Laboratory of Earthquake Engineering and Engineering Vibration, Institute of Engineering Mechanics, CEA, Harbin, China, <sup>2</sup> Beijing Institute of Architectural Design, Beijing, China*

Pipes with a diameter of 150 mm, also called DN150, are often connected by grooved fit joints and employed as stem pipelines, which are used to transport water vertically to different building stories and distribute it horizontally to different rooms. A large deformability is often required for the grooved fit joints to accommodate the deformation concentrated between adjacent stories during an earthquake. To this end, the grooved fit joint is often improved with a wider groove to achieve such a large deformability. However, its seismic performance has not been thoroughly studied yet. This study conducted quasi-static tests on twelve DN150 grooved fit joints, including four elbow joints and eight DN150-DN80 Tee joints. The mechanical behavior, rotational capacity and failure mode were examined and discussed. The test results indicate that the fracture of the grooved fitting and the pull-out of pipes from the grooved fitting are the major damage patterns at deformations larger than 0.1 rad. At small deformations of <0.06 rad, although slight abrasions and wear were observed on the contact surface between the galvanized steel pipe and the grooved fitting, they would not result in significant leakage. Three damage states are defined accordingly, and the fragility models are developed for different grooved fit joints based on test results. Finally, seismic fragility analysis of DN150 stem pipeline system in a 10-story building was conducted. It is demonstrated that the improved joints survive under the maximum credible earthquake and the leakage is highly unlikely to occur.

#### Edited by:

*Antonio Formisano, University of Naples Federico II, Italy*

#### Reviewed by:

*Fabrizio Paolacci, Roma Tre University, Italy Huseyin Bilgin, Epoka University, Albania*

> \*Correspondence: *Tao Wang wangtao@iem.ac.cn*

#### Specialty section:

*This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment*

Received: *19 December 2018* Accepted: *28 March 2019* Published: *17 April 2019*

#### Citation:

*Wang T, Shang Q, Chen X and Li J (2019) Experiments and Fragility Analyses of Piping Systems Connected by Grooved Fit Joints With Large Deformability. Front. Built Environ. 5:49. doi: 10.3389/fbuil.2019.00049*

Keywords: piping system, grooved fit joints, large deformability, quasi-static test, leakage limit state

# HIGHLIGHTS


# INTRODUCTION

Recent earthquakes have demonstrated that the non-structural components of critical facilities such as power plants, hospitals, and industrial units suffered much more damage in comparison with structural components, which mitigated the functionality of these critical facilities (Filiatrault and Sullivan, 2014; Dhakal et al., 2016). The pressurized fire sprinkler piping system is one of the most critical non-structural systems in a building. Its operational function is expected to be maintained after a destructive earthquake to prevent a subsequent fire disaster. However, the piping system was reported to be one of the most vulnerable non-structural systems in past earthquake disasters.

During the 1994 Northridge earthquake, pipelines with a diameter of <25 mm in different piping systems such as HVAC systems, sprinkler piping systems, and water piping systems experienced widespread failures (Fleming, 1998). Leakage of fire sprinkler piping systems forced the temporary evacuation of several hospital buildings after the earthquake (OSHPD, 1995; Ayres and Ezer, 1996). After the 2010 Chile Earthquake, four hospitals in the central south region of the country were inoperable, and 12 hospitals lost almost 75% of their functionalities. Most losses were caused by damages to non-structural components such as suspended ceilings, light fixtures, and fire sprinkler piping systems. The two largest airports in Chile were closed as well because of non-structural damages and flooding from failed sprinkler piping systems (Ju, 2011; Miranda et al., 2012; Tian et al., 2014). One typical fracture failure of the piping tee joints at the Santiago Airport is shown in **Figure 1**, which resulted in significant water leakage and further mitigated the functionality of the airport. **Figure 2** shows the percentages of damaged piping systems over all investigated piping systems in previous earthquakes (Steinbrugge, 1982; McKevitt et al., 1995; Fleming, 1998; Global, 2001; Dillingham and Goel, 2002; LeGrone, 2004; Chock, 2006; Miranda et al., 2012). In the 2011 Tohoku Earthquake, about 8% of the investigated piping systems suffered damages, among which pipelines dominated with 49% of economic loss of the entire

fire sprinkler system, as shown in **Figure 3**. It was also reported that 42% of the piping systems with damaged components exhibited water leakage during this earthquake which caused great inconvenience to earthquake rescue efforts (Mizutani et al., 2012; Soroushian et al., 2015a).

Several studies have investigated the seismic performance of piping systems. Antaki and Guzy (1998) conducted four-point bending tests on 16 pressurized pipe specimens with different types of joints. He et al. (2018) conducted a total of 72 cyclic tests on common water supply pipes to study the failure modes and seismic fragilities of Tee joints. Similar tests were performed by Blasi et al. (2018) to examine the hysteretic response of castiron and copper tee joints. Tian et al. (2014) and Soroushian et al. (2015b) tested 48 Tee joints including grooved fit joints under reversal cyclic loading to determine their rotational capacities. The damage states such as water leakage and fracture were identified. The grooved fit joints featured had a diameter no larger than 100 mm. Although the mechanical performance was stable, the leakage deformation was about 0.02 rad, implying high possibility of leakage during a rarely occurred earthquake. The seismic performance of piping system is also often examined by shaking table tests apart from quasi-static tests. Zaghi et al. (2012) investigated the seismic characteristics of welded and threaded tee joints of hospital piping systems under various intensities of seismic loading using a biaxial shake table. Tian et al. (2015) tested three full-scale pressurized sprinkler piping specimens with different materials, joint arrangements, and types of seismic bracings under dynamic loading. Jochen et al. (2018) studied the swing of CPVC sprinkler branch lines with tee joints under strong earthquake motions. Soroushian et al. (2014b) used the results of two shaking table experiments and numerical simulation of a hospital piping system to evaluate the probability of seismic damage to the piping system. The damage to the piping system itself and the pounding effects with adjacent objects such as walls, floors, and suspended ceilings are involved in the analysis. Threedimensional dynamic response analyses on a suspended thin-wall water piping system under seismic floor motions were conducted by Tatarsky and Filiatrault (2019) to evaluate the fragility of thin-wall piping system.

Recently most public buildings in China have started to adopt pipes with large diameters, such as 150 mm (DN150), as the stem pipelines for fire sprinkler systems. The pipes are connected through grooved fit joints which are improved to accommodate larger rotations to avoid leakage at large deformations. Chinese style joints are quite different from American style ones and the seismic performance needs further study.

In this study, the configuration of the improved grooved fit joint is first introduced. Then, a series of quasi-static tests are conducted. The mechanical behavior, the rotational capacity and the failure mode are investigated, and fragility models are developed considering various damage states. Finally, a three-dimensional finite element model is constructed for a 10-story RC building where the DN150 pipeline system is installed to transfer water vertically along the height. Based on the probabilistic seismic demand model and fragility model obtained by experiments, probabilistic seismic fragility analysis is conducted for piping systems.

## QUASI-STATIC TEST OF GROOVED FIT PIPING JOINTS

# Grooved Fit Joints With Large Deformation Capacity

Grooved fit joints, as shown in **Figure 4**, are a relatively new piping construction product that has gained popularity in earthquake prone areas because of the improved flexibility they provide to fire sprinkler piping systems. The grooved fit joint is composed of pipes with grooves, lathedog-plumbings, rubber gasket and bolts. Dimensions of 90◦ Elbow joint and Tee joint are presented in **Figures 4B,C**. Nominal size of 90◦ Elbow joint is 150 mm while the actual size is 165 mm. Nominal size of Tee joint is 150 × 80 mm while the actual size is 160 × 90 mm. **Figures 4D,E** present the deformation mechanism of the grooved fit joint. The gaps in the grooves and between the pipes make it possible to accommodate earthquake induced rotations. The lathedog-plumbings and pipe joints are made of cast iron and the water pipes are galvanized steel pipes. Detailed information of the materials and components can be found in Meide (2016).

Although this type of joint is massively used in piping construction in China, there's no sufficient investigation on its mechanical behavior, deformation capability, and failure mode. This study conducted three sets of quasi-static tests to examine the seismic performance of two types of grooved fit joints. The elbow joint employs a 90◦ elbow joint connected with two DN150 pipes, as shown in **Figure 5A**. The length of each pipe is 840 mm. The Tee joint uses a triplet connecting to two DN150 pipes and one DN80 pipe. The length of DN150 pipes is 840 mm, while the DN80 is 350 mm. Eight Tee joints are divided into two groups. One tests the rotation of DN150 pipes and the other tests the behavior of DN80, as shown in **Figures 5B,C**, respectively. In each group, there are one monotonic loading test and three cyclic loading tests.

Installation procedure of grooved fit joints is presented in **Figure 6**. Before the installation, lubricants are applied to the surface of rubber gasket. It is necessary to ensure that there is no gap between pipe end and joint. When tightening the bolts, take turns on both sides until the torque meets the installation requirements (60 Nm) and then the installation is completed.

# Experimental Setup and Loading Protocol

The functionality examination is one of the objectives of the experiment. The joints shall be tested with the internal water pressure. The maximum working pressure of water pipeline is 1.2 MPa. According to requirements of GB 50084-2001 (2001) and GB 50261-2005 (2005), the test pressure shall be increased by 0.4 MPa when the designed working pressure is >1.0 MPa. Therefore, all specimens will be filled with water and pressurized to 1.6 MPa to capture the leakage during the test and simulate the working condition of pipes. To achieve a precise water pressure, each DN150 pipe is perforated and welded with a DN32 pipe. One DN32 pipe is used as the water injection orifice to the water pump, while the other for the air escape and equipped with a ball valve to shut off the orifice once all the air has escaped. Then the water pressure will increase. Once the target pressure, 1.6 MPa in this study, is reached, the valve close to the pump will be shut off to maintain the water pressure. The pipes full of pressurized water will be kept for 24 h before tests to examine the sealing performance.

The Tee joints are to be tested in two patterns, one examining the rotation capacity of the triplet connected to DN150 pipes, the other testing the capacity of the triplet to the DN80 pipe, as shown in **Figures 8B,C**, respectively. The test setup for Tee joints (DN150) is shown in **Figure 7**. One hydraulic-servo actuator is employed to realize the quasi-static tests at a low loading speed of 0.5 mm/s. One end of the actuator is fixed on the reaction beam which is securely fixed on the strong floor. The other end is attached to a load jig which is bolted on two linear bearing

Wang et al. Fragility Analysis of Piping System

sliders. The sliders restrain the movement of the load jig in the axial direction only. The end of each DN150 pipe is pinned to one reaction pier fixed on the strong floor, while the end of DN80 is connected to the load jig along the loading direction (see **Figures 7**, **8B**). Similar set-ups are used to test the Tee joints (DN80) and the elbow joints. When testing the DN80 connection of the Tee joints, the DN80 pipe is pinned to the load jig perpendicular to the loading direction (see **Figure 8C**). It is worth noting that the load jig is lengthened using a steel component as shown in **Figure 8C**. The elbow joints are pinned to the load jig at one end and to the reaction beam at the other, as shown in **Figure 8A**. The axes of the initial positions of the two DN150 pipes of the specimen are 45 degrees from the horizontal line. More details of the test setup are presented in Shang et al. (2018).

The specimens were subjected to a cyclic loading following loading protocol suggested by FEMA 461 (2007), as shown in **Figure 9**. This protocol is commonly used for the seismic performance assessment of non-structural components and equipment, as it allows all damage states to be quantified to develop the corresponding fragility models. The loading protocol consists of repeated cycles with an amplitude that steadily increases by 1.4 times at each step. Two cycles of loading are conducted at each amplitude level. At least six cycles must be performed before any damage can be observed.

#### Measurement Scheme

The bending moment of the Tee joints (DN150), M, was calculated using Equation (1), where F is the load applied by the actuator, and D<sup>i</sup> is the distance from the pin end of the pipe to the center of the triplet. The bending moment of the Tee joints (DN80) and the Elbow joints can be calculated in a similar way:

$$\mathcal{M} = F \cdot D\_i / 2 \tag{1}$$

Full field displacement and strain measurement with high speed 3D digital image correlation (DIC) system is used for the measurement of strain and rotation in this study. DIC method is a powerful digital image correlation technique that gives the full-field strain distribution across a specimen (Reedlunn et al., 2011). The pipe surface is sprayed with black and white speckle. The measurement system tracks the speckle changes and uses the optimized 3D DIC algorithm to calculate the displacement and strain. The measured strain using DIC method can cover the range from 0.005 (50 micro-strains) to 2,000%. During the loading process, water leakage will first happen at one end of the joints. Therefore, the rotations of all pipes connected to the joints are measured and the one with the first damage will be adopted to quantify the damage states. The rotation is defined as the relative rotation angle between the joint and the pipeline and can be obtained in the post-processing software of the high speed 3D digital image correlation (DIC) system.



Meanwhile, four strain gauges are pasted on each pipe of the specimen, as shown in **Figure 10**. The curvature of each section is measured by two strain gauges respectively attached to the top and bottom of the pipes, as formulated in Equation (2), where, κ represents the curvature of the measured section, ε<sup>1</sup> and ε<sup>2</sup> are the strains measured from the strain gauges, and s is the interval between two strain gauges (pipe diameter in this study). The moment M of the section is then calculated with Equation (3), where E is the young's modulus, and I is the moment of inertia of the concerned section. The distribution of bending moment along the pipe is assumed to be linear. The calculated bending moments of the two measuring points on the failure side are used to obtain that at the pipe joint section by linear assumption. The data measured from strain gauges will be used to calibrate the optical measurement and compared with the bending moment calculated using the method of Equation (1).

$$\kappa = \frac{\varepsilon\_1 - \varepsilon\_2}{s} \tag{2}$$

$$M = EI\kappa \tag{3}$$

#### EXPERIMENTAL RESULTS OF GROOVED FIT PIPING JOINTS

### Physical Damages of Grooved Fit Piping Joints

The observed damage at first leakage was consistent for all specimens tested. **Table 1** lists and shows photographs of the observed damages. The main failure modes were the fracture of coupling flanges and the pulling-out of pipes. Slight wear at interface between pipes and joint can be observed at DN150 pipes and heavy wear can be observed at DN80 pipes but the wear didn't cause any leakage and did no harm to the functionality of the pipes, which implies a good sealing performance of the improved joint.

#### Rotation Capacities of Grooved Fit Joints

Three damage states (DS) of grooved fit joints are defined based on the test results. DS<sup>1</sup> represents the case in which pressure in the pipe drops to 85%, DS<sup>2</sup> represents massive leakage when the water pressure drops to 0 MPa, and DS<sup>3</sup> represents complete damage of the joints, including the resistance force dropping to 85% of the peak force, visible cracks, and the pipe pull-out. The rotations corresponding to different damage states for each specimen are compared in **Figure 11**, both monotonic and cyclic test results are given. All types of joints exhibit large rotational capacities ranging from 0.054 to 0.274 rad before massive leakage. Rotational capacity of the DN80 Tee joint is much larger than that of the DN150 Tee joint and the Elbow joint. The monotonic rotational capacities of massive leakage (DS2) are larger than their corresponding cyclic rotational capacities. This result is consistent with the results of Tian et al. (2014).

## Moment-Rotation Cyclic Response of Grooved Fit Joints

The measurement scheme provides two methods to obtain the bending moment of the concerned sections. **Figure 12** shows the measured hysteretic curves from DIC and strain gauges, respectively. The moment of DIC curve is calculated using the

#### TABLE 2 | Summary of fragility parameters.


method of Equation (1), while that of strain gauge is calculated using the method of Equations (2, 3). They measured generally well. And in the following discussion, the data measured from DIC will be used. Typical moment vs. rotation hysteretic curve of each type of grooved fit joint is shown **Figure 13**. The occurrence of massive leakage (DS2) is indicated by a solid red dot on each plot. After massive leakage, the tests were continued up to complete damage (DS3) of the test specimens. It can be observed from the hysteretic loops that the moment values would keep increasing even at a large rotation around 0.1 rad. And a pinching effect was introduced by the gaps in the improved configuration. As soon as the load intensity was increased, during the load-inversion phase, the gap generated caused relative rotations with near-zero force variation. Furthermore, reduction of the reloading stiffness can be observed after each loading step because the response could be highly influenced by the cumulative damage.

# Experimental Fragility Analysis of Grooved Fit Joints

The experimental results of the cyclic tests were processed to develop the seismic fragility models for grooved fit joints. The cyclic behavior of the piping joints was governed primarily by joint rotation, and this is the only engineering demand parameter (EDP) considered. The three damage states were considered in the seismic fragility analysis, i.e., pressure dropping, massive leakage, and complete mechanical failure, corresponding to DS1, DS2, and DS3, respectively. Based on the method proposed by Porter et al. (2007), experimental fragility curves were defined for the three joint types based on the measured rotational capacities given in **Figure 11**. The median rotational capacity, θm, and associated logarithmic standard deviation, β<sup>t</sup> , were computed by Equations (4, 5).

$$\theta\_m = \exp\left(\frac{1}{M} \sum\_{i=1}^{M} \ln \theta\_i\right) \tag{4}$$

$$\beta\_t = \sqrt{\frac{1}{M-1} \sum\_{i=1}^{M} \left( \ln(\theta\_i / \theta\_m) \right)^2} \tag{5}$$

$$
\beta = \sqrt{\beta\_t^2 + \beta\_u^2} \tag{6}
$$

where θ<sup>i</sup> denotes the i-th measured rotational capacity and M is the number of cyclic tests conducted for this type of joints (M = 3 in this study). Considering the fact that the specimen number is <5 and all specimens experienced the same loading history, a correction factor β<sup>u</sup> = 0.25 is introduced in the calculation of logarithmic standard deviation and the modified logarithmic standard deviation β is calculated by Equation (6). **Table 2** summarizes the median rotational capacity and logarithmic standard deviation obtained for each type of joint. **Figure 14** compares all the fragility curves derived from the experimental data. The fragility curves of water leakage damage state (DS2) are further compared with the results of Tian et al. (2014) in **Figure 15**, where GFC50.8 (3.8) represents the pipes with diameter of 50.8 mm and the pipe wall thickness is 3.8 mm. It can be found that the rotation capacities of Chinese style joints are quite larger than those of American style joints. It is worth noting that more tests are needed for detail comparison to consider the influence of pipe diameter and pipe wall thickness.

#### BUILDING SPECIFICATIONS FOR THE CASE STUDY

A 10-story reinforced concrete (RC) frame was used for the case study of grooved fit piping joints. The story heights for the first two stories are 3.8 and 3.95 m, respectively, while the height of the penthouse is 3.9 m. The height of the remaining stories is 3.45 m. Uniformly distributed dead and live loads on each floor were 3.0 and 2.0 kN/m<sup>2</sup> , respectively and those for the roof were 4.0 and 0.5 kN/m<sup>2</sup> , respectively. It was designed following the Chinese Code for Seismic Design of Buildings (GB 50011-2010, 2010). The PGA of the design-basis earthquake ground motion was 0.2 g. The site was classified as type II and the design characteristic period was 0.35 s. The layout of the building was significantly irregular as shown in **Figure 16**. The concrete was C30 with a standard compressive strength of 20.1 MPa and the steel rebars were HRB400 with a yielding strength of 360 MPa.

A finite element model of the building was created in the Open System for Earthquake Engineering Simulation (OpenSees) software (2012). The beams and columns were simulated using non-linear beam-column elements, which were defined as forcebased elements with distributed plasticity. Five Gauss-Lobatto integration points were inserted along the length of each element with two at the ends to simulate the formation of plastic hinges. The composite RC cross-sections were conveniently simulated by the fiber formulation. Concrete material with a



linearly decaying tensile strength proposed by Hisham and Yassin (1994) was adopted. The input parameters of the confined concrete were determined using Mander's theoretical stressstrain model (Mander et al., 1988). A Giuffre-Menegotto-Pinto

uniaxial material model (Menegotto and Pinto, 1973) with isotropic strain hardening was used to model the longitudinal steel rebars. The story mass and the corresponding gravity were uniformly distributed to each beam-column joint and the rigid diaphragm constraint was used for each floor. P-Delta effect was considered for the geometric non-linearity. A Rayleigh damping of 5% was assigned to the first two modes of vibration. The first three fundamental natural periods of the building were 1.2362, 1.1599, and 1.0772 s, respectively. Detailed information can be found in Shang et al. (in press).

The far-field record set that includes 22 records was selected from the Pacific Earthquake Engineering Research Center Next Generation Attenuation (PEER NGA) database (Ancheta et al., 2013) and the criteria suggested by FEMA P695 (2009) were used. These ground motions were adopted to conduct non-linear time history analysis for evaluating the seismic performance of the building. The acceleration response spectra of the selected ground motions are presented in **Figure 17**.

## FRAGILITY ANALYSIS OF GROOVED FIT PIPING JOINTS

#### Probabilistic Seismic Demand Model

A probabilistic seismic demand model (PSDM) can be obtained through regression analysis using the computed responses. The relationship between the interstory drift ratio (IDR) of each story and ground motion intensity was estimated using a standard power function, as presented by Li and Ellingwood (2007):

$$S\_D = a(IM)^b \tag{7}$$

which can be rewritten in the logarithmically transformed space as:

$$
\ln\left(S\_D\right) = \ln\left(a\right) + b\ln\left(IM\right) \tag{8}
$$

where S<sup>D</sup> is the conditional median interstory drift ratio demand, a and b are unknown regression coefficients which can be found from linear regression in the log space, and IM is the ground

motion intensity measure [Sa(T1) in this study]. A conditional logarithmic standard deviation, βD|IM, was estimated based on the calculated data to consider the uncertainty associated with the demand model (Li and Ellingwood, 2008; Tavares et al., 2012):

$$\beta\_{D|IM} \cong \sqrt{\frac{\sum\_{i=1}^{N} \left(\ln\left(d\_i\right) - \ln\left(a(IM)^b\right)\right)^2}{N-2}}\tag{9}$$

where N is the number of simulations and d<sup>i</sup> is the calculated maximum interstory drift ratio. **Figure 18** shows the result of the probabilistic seismic demand analysis for first floor of the 10 story building. R 2 is the coefficient of determination which can be used to indicate the robustness of the regression. Parameters used to define the probabilistic seismic demand model of the case study building are given in **Table 3**.

#### Probabilistic Seismic Fragility Analysis

It is assumed that the pipe joints of vertical pipes are installed at the ends of the pipe, nearly at the same level of the floor slab, as presented in **Figure 19**. We assumed that the rotation on vertical pipes imposed by building deformation (drift) is only concentrated at end joint rotations, so that the joint rotation is equal to the interstory drift ratio (IDR) (Soroushian et al., 2014a). The following assumptions have been made in this study:

(1) flexural and shear deformation of pipe segments is negligible;

(2) the top and bottom of vertical pipe at each floor is fixed by solid sway braces.

With the developed PSDM of the interstory drift ratio of each story and the capacity as defined by component fragility curve (section Experimental Fragility Analysis of Grooved Fit Joints), the exceedance probability that the pipe joint demand (D) would be larger than the capacity (C) at a given earthquake intensity measure was computed by Equation (10) for the three damage states. It can be rewritten as Equation (11) by substituting the demand median S<sup>D</sup> in the form of Equation (8) (Tavares et al., 2012):

$$P\left[D \ge C\left|IM\right] = \Phi\left(\frac{\ln\left(\mathbb{S}\_{D}/\mathbb{S}\_{C}\right)}{\sqrt{\rho\_{D|IM}^{2} + \rho\_{C}^{2}}}\right) \tag{10}$$

$$P\left[D \ge C\left|IM\right] = \Phi\left(\frac{\ln\left(IM\right) - \frac{\ln(S\_C) - \ln(a)}{b}}{\frac{\sqrt{\beta\_{D|IM}^2 + \beta\_C^2}}{b}}\right) \tag{11}$$

**Figure 20** shows fragility curves of grooved fit joints (DN150) used for vertical pipes using Sa(T1) as an intensity measure, where F<sup>i</sup> (i = 1,2,. . . ,10) represents the fragility curve of i-th floor. The damage probability of grooved fit joints (DN150) used for the vertical pipe under 0.07, 0.20 and 0.40 g [corresponding Sa(T1) are 0.0498, 0.1424, and 0.2847 g] corresponding to the service level earthquake (SLE), design-basis earthquake (DBE), and maximum considered earthquake (MCE) in the Chinese codes (GB50011–2010) is nearly zero. The large deformation capacity of the improved grooved fit joints would ensure the functionality of fire protecting after a rarely occurred earthquake.

#### SUMMARY AND DISCUSSIONS

Monotonic and reverse cyclic testing were conducted on 12 grooved fit joints to evaluate their seismic performance. Three damage states were defined and fragility curves for different grooved fit piping joints were developed based on test results. A probabilistic seismic fragility analysis was conducted on a 10 story frame building equipped with DN150 vertical stem pipeline as the fire sprinkler system. Major findings are as follows:


# REFERENCES


considered earthquake, implying the fire protection would be functional.

This is a preliminary study on the functionality of fire sprinkler system, as part of the series studies of seismic resilience on important buildings. More realistic boundary conditions and configurations of this system shall be considered in the future. The damages, such as the pounding damage between the piping system and nearby objects, anchorage damage of fasteners, and anchors that connect the piping system to structural members, failure of pipe hangers, breaking of sprinkler heads, should be included to comprehensively investigate the seismic vulnerability of sprinkler piping systems in different types of buildings.

#### AUTHOR CONTRIBUTIONS

TW is in charge of the whole project including tests and analysis. QS completed the experiment and analysis. XC helped to prepare the test. JL helped to complete the fragility analysis.

#### FUNDING

This research was funded by the Scientific Research Fund of Institute of Engineering Mechanics, China Earthquake Administration (Grant No. 2017A02), the National Natural Science Foundation of China (Grant Nos. 51378478, 51678538), and the Heilongjiang Provincial Natural Science Fund for Distinguished Young Scholars (Grant No. JC2018018). Any opinions, findings, and conclusion or recommendation expressed in this paper are those of the authors and do not necessarily reflect the views of the sponsors.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Wang, Shang, Chen and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Application of Association Rules to Determine Building Typological Classes for Seismic Damage Predictions at Regional Scale: The Case Study of Basel

#### Lorenzo Diana1,2 \*, Julien Thiriot <sup>1</sup> , Yves Reuland<sup>1</sup> and Pierino Lestuzzi <sup>1</sup>

*<sup>1</sup> Applied Computing and Mechanics Laboratory (IMAC), School of Architecture, Civil and Environmental Engineering (ENAC), Swiss Federal Institute of Technology Lausanne (EPFL), Lausanne, Switzerland, <sup>2</sup> Department of Structural and Geotechnical Engineering (DISG), University of Rome 'La Sapienza', Rome, Italy*

#### Edited by:

*Katsuichiro Goda, University of Western Ontario, Canada*

#### Reviewed by:

*Viviana Iris Novelli, University of Bristol, United Kingdom Xinzheng Lu, Tsinghua University, China*

> \*Correspondence: *Lorenzo Diana lorenzo.diana@epfl.ch*

#### Specialty section:

*This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment*

Received: *31 January 2019* Accepted: *01 April 2019* Published: *17 April 2019*

#### Citation:

*Diana L, Thiriot J, Reuland Y and Lestuzzi P (2019) Application of Association Rules to Determine Building Typological Classes for Seismic Damage Predictions at Regional Scale: The Case Study of Basel. Front. Built Environ. 5:51. doi: 10.3389/fbuil.2019.00051* Assessing seismic vulnerability at large scales requires accurate attribution of individual buildings to more general typological classes that are representative of the seismic behavior of the buildings sharing same attributes. One-by-one evaluation of all buildings is a time-and-money demanding process. Detailed individual evaluations are only suitable for strategic buildings, such as hospitals and other buildings with a central role in the emergency post-earthquake phase. For other buildings simplified approaches are needed. The definition of a taxonomy that contains the most widespread typological classes as well as performing the attribution of the appropriate class to each building are central issues for reliable seismic assessment at large scales. A fast, yet accurate, survey process is needed to attribute a correct class to each building composing the urban system. Even surveying buildings with the goal to determine classes is not as time demanding as detailed evaluations of each building, this process still requires large amounts of time and qualified personnel. However, nowadays several databases are available and provide useful information. In this paper, attributes that are available in such public databases are used to perform class attribution at large scales based on previous data-mining on a small subset of an entire city. The association-rule learning (ARL) is used to find links between building attributes and typological classes. Accuracy of wide spreading these links learned on <250 buildings of a specific district is evaluated in terms of class attribution and seismic vulnerability prediction. By considering only three attributes available on public databases (i.e., period of construction, number of floors, and shape of the roof) the time needed to provide seismic vulnerability scenarios at city scale is significantly reduced, while accuracy is reduced by <5%.

Keywords: seismic vulnerability, urban scale, existing buildings, data mining, damage distribution

# INTRODUCTION

Seismic-risk assessment at urban scale is an essential step toward earthquake-resilient communities; not only in highly seismicprone regions but also in regions with low-to-moderate seismic hazard. The concept of resilience is not exclusively related to safety of inhabitants, but also to the capacity of systems to recover from a seismic event and to get back to previous levels of load-bearing capacity (Bruneau and Reinhorn, 2007; Cimellaro et al., 2010a,b; Alexander, 2013; Burton et al., 2015). Indeed, even slight damages can cause production interruptions and communication breakdowns with potential consequences on the prosperity of entire regions that last for months and even years. Architects, engineers and urban planners should be involved in the pre-earthquake phase, not only in towns with high seismicity, for the definition of scenarios regarding seismic vulnerability. Such scenarios are useful to estimate the levels of damage that urban systems are expected to suffer in case of seismic events. Thus, efficient pro-active actions can be implemented, which help increasing the general resilience.

Evaluating seismic risk at large scale is a complex and wide process of knowledge. Multiple domains are involved: seismic hazard, exposure and seismic vulnerability (Carreño et al., 2007). Several models and datasets exist in literature: WHE—World Housing Encyclopedia (EERI, 2004); FEMA—Federal Emergency Management Agency (ATC (Applied Technology Council), 2005; Pager—Prompt Assessment of Global Earthquakes for Response (Jaiswal and Wald, 2008); GEM—Global Earthquake Model (Brzev et al., 2013). The time demand and the expensiveness of evaluations depend on the type of detail required. In order to develop scenarios of urban vulnerability, several thousands of buildings would need to be evaluated. Lately, researches have been developed for large scale building-bybuilding evaluations (Yamashita et al., 2011; Xiong et al., 2018). In general, building-specific and detailed evaluation of the entire building stock is undermined by the economic and technical needs of such a process, even more so in regions with low-to-moderate seismic hazard. Therefore, simplifications are necessary. First, a reduction in the number of assessed buildings is performed: buildings are clustered into typological classes and the vulnerability of each class—rather than each building—is calculated in detail.

A fundamental starting point to urban assessment is thus the definition of an appropriate taxonomy. The taxonomy is particularly important because several typological classes are introduced to describe accurately the structural behavior of the existing building stock. In Europe, several studies have defined taxonomies: starting with the EMS-98 (Grünthal et al., 2001) with the definition of building classes and vulnerability classes and then with the Risk-UE project (Lagomarsino and Giovinazzi, 2006) who have implemented the EMS-98 classes by adding some specifications.

Once the taxonomy is defined, a specific typological class can be attributed to each building. This attribution process is essential to achieve meaningful seismic evaluations at urban scale, since the estimation of the expected damage is directly related to the building type. Attributing types to buildings is commonly based on surveys (ATC (Applied Technology Council), 2005). Several approaches exist to survey existing buildings: starting with the well-established building-by-building visual inspection, which at urban scale has the drawback of being highly time demanding (McEntire and Cope, 2004; Marquis et al., 2017). Over the last few years, surveys involving drones, satellites, and open maps have been introduced to reduce the time demand (Suzuki et al., 2010; Ehrlich et al., 2013). Deep-learning approaches, collecting data from pictures of buildings or other sources are then implemented for automatic attribution of types (Mallepudi et al., 2011) or directly for seismic-vulnerability estimates (Alizadeh et al., 2018).

At the current state of the art, an increase in the speed of the survey phase (i.e., attribution of typological classes to all buildings and consequently the estimation of the expected damage) can be obtained by using building data that is readily available in public databases (Riedel et al., 2014). By the selection of appropriate attributes (such as number of floors, material and year of construction) and the application of data-mining methods, correlations between building attributes and typological classes on small learning sets can be defined (Riedel et al., 2014). The association-rule learning (ARL) is a data-mining method (Agrawal et al., 1993) that is based on finding association rules between attributes of buildings that help discovering statistical links between building features. Once correlations are defined on the learning set, a class (or a probability to belong to a class) can be attributed to each building that is subsequently analyzed. The ARL method allows to use available building attributes to assign buildings to predefined classes and, by extension, to define their vulnerability. Vulnerability represents the intrinsic predisposition of a building to be affected and to suffer damage following the occurrence of a given event (Guéguen, 2013).

When dealing with vulnerability analysis at larger scales, the aim is to assess the impact of an earthquake on a set of buildings within an area of interest. Two main approaches exist for seismicvulnerability assessment of existing buildings at large scale: empirical (or macro-seismic) methods and mechanical methods (Lestuzzi et al., 2016). In empirical methods, the vulnerability of each class is measured in terms of a vulnerability index, V, that is calculated based on the observed damage of buildings of every class in past earthquakes. Mechanical methods are based on a model-based evaluation of the structural behavior of buildings: by the interaction between the structural behavior (identified for example by capacity curves) and the seismic demand (identified for example by response spectra), the expected damage reached by a typological class is determined.

In Europe, a fundamental research project for both empirical and mechanical method is the Risk-UE project. The Risk-UE project represents the first collaborative and comprehensive research program that studied territorial seismic risk focused on the European built environment (Lestuzzi et al., 2017).

Within the Risk-UE project (Mouroux et al., 2004; Mouroux and Le Brun, 2006), which proposes "an advanced approach to earthquake risk scenarios with applications to different European towns," the vulnerability of existing buildings is evaluated according to the two approaches: Level 1 or LM1, based on the empirical method; and Level 2 or LM2, based on the mechanical method. The empirical method, LM1, involves the macro-seismic Diana et al. ARL for Seismic Damage Predictions

intensity according to EMS-98 (Grünthal et al., 2001) and vulnerability and ductility indexes. Lagomarsino and Giovinazzi (2006) proposed a correlation between macro-seismic intensity of the European macro-seismic scale EMS-98 and building damage. This correlation is shown in terms of vulnerability curves for each building type. The vulnerability indexes are derived from the vulnerability curves of EMS-98 (Giovinazzi and Lagomarsino, 2001). In EMS-98, the expected damage of each type is defined using linguistic terms ("few," "many," "most") considering five damage grades. With the LM1 methods of the Risk-UE project, these terms have been transformed into numbers (parameter V) applying an implicit Damage Probability Matrix and using the fuzzy set theory (Dubois and Parade, 1980). For LM1 calculations, the European macro-seismic scale, EMS-98, defines the hazard and the damage-grade scale (from D1 to D5).

Mechanical models (LM2) are based on the structural response of the buildings, expressed by the force-displacement curve (capacity curve). Three parameters are needed to represent capacity curves in a simplified elastic-perfectly plastic model: dy—yield displacement; Ay—yield acceleration; du—ultimate displacement. In general, several mechanical methods exist in literature, many of them are based on the Capacity Spectrum Method of ATC-40 (ATC (Applied Technology Council), 1996). Within the framework of this paper, vulnerability is calculated using the macro-seismic approach (LM1).

In this paper, the application of the ARL method for attributing typological classes (for seismic evaluation) to buildings is proposed. The methodology is then applied to the city of Basel, in Switzerland, where several datasets are available containing information of all buildings. The analyzed datasets provide elementary attributes and characteristics of buildings. The main goal of the paper is to evaluate the performance in the attribution of typological classes on a learning set (containing more than 700 buildings), considering various combinations of attributes. In a second phase, the distribution matrixes obtained are applied to the entire city of Basel, which has been partially surveyed as part of a Master thesis at the IMAC Lab, EPFL (Thiriot, 2019).

# METHODOLOGY

When assessing the seismic vulnerability of existing buildings at large scale, the definition of typological classes and their attribution to the building stock of an urban system are elementary. The building-by-building evaluation is a process that can hardly be performed at urban scales, since it requires amounts of time and money that are not available, especially in regions with low-to-moderate seismicity that have reduced mobilization of resources (Riedel et al., 2015). Thus, methods are needed to attribute, in a simplified way, a specific typological class to each building in a region of interest. Such a list of typological classes (or types) that compose the urban system is called taxonomy (Porter et al., 2001). The classes proposed by Lagomarsino and Giovinazzi (2006) have been adopted by the Risk-UE project. This taxonomy is an improvement of the TABLE 1 | Typological classes that are considered for the city of Basel.


previous taxonomy proposed in EMS-98. In order to provide an example of a taxonomy, all building classes that are considered to form part of the city of Basel, which is considered as a case study in this paper (see section Case Study) are presented in **Table 1**. The typological class definitions can be retrieved in the Risk-EU project and in the EMS-98.

The aim of the presented methodology is to find correlations between building attributes, such as period of construction, number of floors (related to the building height), shape of the roof and footprint surface [as described in section Collection of Building Attributes (STEP 1)], which are available in large datasets, and typological classes (see **Table 1**). By doing so, a given class is attributed, with a certain probability, to a given building. The process of identifying patterns, defining correlations and attributing types is based on data-mining large datasets. The method used for the determination or correlations between attributes and typological types is the ARL (Association Rule Learning) method. In detail, understanding the influence of attributes and their combination on the general performance of the proposed class attribution is investigated. The flowchart presented in **Figure 1** describes the methodology applied in the present work explaining how typological classes are attributed and vulnerability is assessed for large-scale studies: in step 1 attributes of buildings are gathered from large-scale databases; in step 2 a data-mining process is applied on a training set of buildings to find correlations between attributes and typological classes; in step 3 typological classes are assigned to buildings; and finally in step 4 vulnerability assessment is performed. Two validation phases are included: one at the end of step 3 regarding the performance of class attribution; one at the end of step 4 regarding the performance of vulnerability evaluations. The validation process can only be performed on the parts of the city for which information concerning the real typological classes is available, for instance from visual inspection of buildings.

Therefore, the two validation phases that are involved in the cityscale steps (steps 3 and 4), are performed exclusively on a part of the learning set where knowledge is complete. As the distribution of building types is not homogeneous within a city, the case study analyzed in this paper involves a further validation set, which is lying outside the learning set (see section Evaluation of the Accuracy of Damage Distributions).

In the following sections, the four steps involved in learning correlations between building attributes and typological classes in order to determine city-scale seismic vulnerability are described.

#### Collection of Building Attributes (STEP 1)

The first step consists in the selection of the attributes needed to describe the analyzed building stock. While several other parameters could be possibly chosen from databases, such as number of inhabitants and monetary value of buildings, four building attributes are collected in this paper in order to apply the ARL method: the period of construction; the number of stories; the footprint surface; and the shape of the roof. These attributes are selected because they are considered to fulfill two criteria: having an influence on the structural behavior of buildings and being retrievable from city-wide databases. Thus, these four parameters are estimated to facilitate understanding of correlations between building attributes and typological classes. The period of construction, the number of stories and the shape of the roof have been chosen accordingly to the previous work of Riedel et al. (2015) while the footprint surface has been added in order to evaluate its influence on the final results. Other valuable attributes, such as shape and percentage of openings, are not considered as few current databases contain this type of attributes.

For applications in Switzerland, the date of construction can be taken from the Federal Register of Buildings and Dwelling (RegBL) of the Federal Statistical Office (BFS in German/OFS in French). This register includes all residential buildings in Switzerland and their dwellings. Data is updated on a trimester basis. Five intervals for the period of construction are defined: before 1900; 1900–1944; 1945–1969; 1970–1988; and after 1989. These intervals correspond to major construction periods: the year 1900 corresponds to a turning point in the development of brick constructions; 1945 and the end of the World War II, has led to a widespreading of new methods in the construction industry; the year 1970 corresponds to the first and very basic seismic considerations in Swiss building codes; and finally in 1989 real seismic considerations have been introduced in Switzerland (Lestuzzi and Badoux, 2013).

The number of floors and the footprint surface are derived from the same database (RegBL). The intervals chosen for the number of floors are: less than 4 stories; between 4 and 6 stories; 7 or more stories. As well other data, such as the official number of the building (EGID number) and geo-coordinates, that are useful for the combination of attributes and for figure editing, are available on the same database.

The shape of the roof (flat/sloped) is not available in the Federal Register of Buildings and Dwellings. This attribute provides information regarding the nature of the roof material and, as a consequence, regarding the construction materials of the entire building. Information concerning the shape of the roof is therefore taken from another source. For the case study of section Case Study, the 3D Model of the city of Basel is considered. This model contains all the buildings of the city in 3D with their roof. Each 3D model corresponds to a single building with its proper EGID (federal building identifier) reference number. Therefore, it is possible to link every 3D with the corresponding building in the RegBL database.

The four attributes selected have been collected for all the buildings composing the asset of the city of Basel.

### The Association-Rule-Learning Method (STEP 2)

An essential step of the methodology consists in learning correlations between building attributes and typological classes on a learning subset of the city, for which attributes and typological classes are both known. The method used within the framework of this paper for the attribution of typological classes to buildings is the Association Rule Learning (ARL) method (Riedel et al., 2014): a method for determining relationships between variables in large databases (such as databases of buildings in a city). It has been introduced by Agrawal et al. (1993) as a list of if/then statements that help discovering relationships between seemingly unrelated data. The combination of attributes of buildings (Yj) and the typological classes, X<sup>i</sup> , are linked together by a conditional matrix that is derived from the learning set. The probability that the building belongs to a specific class X<sup>i</sup> knowing that the combination Y<sup>j</sup> has a non-zero probability is written P (X<sup>i</sup> |Y<sup>j</sup> ) and is defined by:

$$P\left(\mathbf{X}\_{i}|\mathbf{Y}\_{j}\right) = \frac{P\left(\mathbf{X}\_{i}\cap\mathbf{Y}\_{j}\right)}{P\left(\mathbf{Y}\_{j}\right)}\tag{1}$$

This method has first been used for seismic applications in the assessment of the city of Grenoble (France) (Riedel et al., 2014). In the case of Grenoble, the vulnerability classes (A, B, C, D and E) of EMS-98 have been considered. Used attributes have been those available in the INSEE (French Institute of Statistics) database, such as the number of floors and the construction period. In this paper, the classes of the Risk-UE project are considered. As mentioned before, the object of the study is the city of Basel in Switzerland.

#### Definition of the Learning Set

For the learning phase, a small subset of the building stock of the city should be selected. On this learning subset buildings are evaluated individually in order to derive the typological class and the building attributes (as shown in step 2 of **Figure 1**). Ideally, classes are attributed based on original drawings. However, often such drawings cannot be retrieved and specific types are attributed based on expert-conducted visual inspection. A support for decision-making that takes the form of a decision tree is developed to help data collectors in the field to classify buildings and to guarantee uniform and correct attributions. When typological class attribution is based on visual inspection, misclassifications cannot be excluded. The number of floors—for seismic evaluations—is taken as the number of vibrating masses (number of floor masses in a multiple-degree-of-freedom model). At the end of the survey campaign, all the buildings composing the learning set correspond to a specific typological class (see **Table 1**) and are characterized by the four attributes collected in step 1.

As first step for the application of the ARL method and to enable a subsequent validation step, the learning set is divided into two parts. The first part is the training set and is composed of randomly selected 30% of the buildings forming the learning set. The training set is used to define correlations between building attributes and typological classes. The remaining 70% of the surveyed buildings are used to evaluate the accuracy of class attribution using the correlations that are defined from the training set. Riedel et al. (2014) have shown that using 30% of the learning set for training results in stable values for assessing the performance of the correlations for type attribution.

The data-mining process defining the ARL method consists in generating a distribution matrix that contains the probability for a building—defined by a combination of the selected attributes to belong to each typological class. Once the distribution matrix is derived from the training set, the attribute combination allows each building to be assigned with a typological class.

#### Typology Determination (STEP 3)

The third step consists in the attribution of typological classes to all the buildings composing the city. For the buildings that are not part of the training set, the real typological class is not available since performing a building-by-building survey campaign is hindered by the associated time and money demands. Therefore, for the whole city, only the attributes collected in step 1 are available.

Once correlations between building attributes and typological classes are derived, classes can be attributed following three approaches. The first two methods—called Pmax and Prandom– associate each building to one specific class. The Pmax method attributes each building to the class that has the maximum probability among the same attribute combination. The second method, Prandom, attributes a class to buildings according to random selection. A cumulative probability is calculated for each combination and then, a number from the interval [0, 1] is randomly generated. The random number defines the class of the building from the cumulative probability distribution. The third method—called Pdistribution–does not associate each building to a single class; it associates to each building the probability to belong to each class. For more specifications, see section Wide-Spreading the Typological Association to the Entire City of Basel (STEP 3).

An important step related to city-scale applications of association rules, which are learned on a small subset (training



set) of buildings, is validation on a set (validation set) of buildings that are not used to learn association rules and for which real typological classes are known. Depending on the attribution method, various approaches for validation are available, as discussed in the following section.

#### Validation of Typology Attribution

The accuracy of distribution matrixes can be assessed in two ways, which are related to the two approaches: attributing classes to buildings (association of one specific type to each building, using either Pmax or Prandom) or repartition of probability values corresponding to each type (Pdistribution). For both ways of attribution, the accuracy can be evaluated based on the total distribution of classes. Only for the case of associating one specific typological class to each building, the accuracy of the distribution matrix can also be evaluated using a confusion matrix. **Table 2** summarizes the approaches that are used to associate classes to buildings and to evaluate the accuracy on the validation set.

Within this paper, accuracy is defined with respect to both, false positives and false negatives. However, when a typology distribution is used to evaluate building attribution, only false negatives affect the accuracy (while false positives tend to increase the accuracy). On the other hand, a confusion matrix allows, in the case of direct association of a specific type to each building (Pmax and Prandom), to compare the attribution of buildings to typological classes by the ARL method with the attribution that is obtained by visual survey. An example of such a confusion matrix is shown in **Table 4** (section Evaluation of the Typological-Attribution Accuracy on the Learning Set). Columns of this matrix correspond to the "real" classes attributed, which are obtained by visual survey. Rows correspond to classes that are associated by the ARL method. The values of the diagonal correspond to correctly associated buildings: a building on the diagonal is assigned to the same class by the ARL method and by visual surveys. Buildings that are outside of the diagonal are assigned by the Pmax–ARL method to another class than based on visual surveys.

For confusion matrixes, the accuracy is defined as the number of buildings, for which classes are correctly attributed (elements of the diagonal, Aii), divided by the total number of buildings (Stehman, 1997). In other terms, the accuracy is the sum of the diagonal divided by the sum of all elements (correctly attributed elements are limited to true positives as in a confusion matrix true negatives for one class are equivalent to true positives for another class):

$$accuracy\_{conf.} = \frac{\sum\_{i} A\_{ii}}{\sum\_{i} \sum\_{j} A\_{ij}} \tag{2}$$

Thus, confusion matrixes provide the accuracy with respect to false positives and false negatives. An additional step can be done by separating the assessment of accuracy in columns or lines: on one hand, non-diagonal elements in one column provide false negatives and thus, allow calculating the recall score; on the other hand, non-diagonal elements in one row provide false positives and thus, allow quantifying the precision.

Evaluating classification accuracy using probability distributions of classes involves checking whether the general distribution of each class is conserved between the ARL method and visual surveys. The real distribution of classes is the outcome of visual surveying. This distribution is compared with the distribution provided by the ARL method. The distribution corresponding to the ARL method is obtained by multiplying the probability distribution for a given combination of attributes and the number of buildings having this given attribute combination. The accuracy is derived as the difference (in number of buildings) between the "real" and the "derived" distribution (see for an example **Table 5** at section Evaluation of the Typological-Attribution Accuracy on the Learning Set):

$$accuracy\_{distr.} = 1 - \frac{\sum\_{k} \left| class\_{k\_{REAL}} - class\_{k\_{ARL}} \right|}{tot.\,\text{build}.} \tag{3}$$

The main difference between the two validation methods is based on the fact that the confusion matrix allows a building-bybuilding attribution. As a consequence, all miscategorizations are considered as errors. The accuracy values obtained with Equation (2) are therefore lower (accuracy is checked with respect to false positives and false negatives) than the accuracy obtained with the probability distribution. Indeed assessing accuracy with respect to distributions, is a comparison of the sum of columns of the confusion matrix (real distribution) with the sum of rows of the confusion matrix (ARL-based type distribution) and thus, compensation between errors increases accuracy. Probability evaluation is commonly used for large-scale evaluations. It is a mean to understand the total distribution of types. The goal is to establish the total number of buildings in a certain typological class. Thus, errors (whether false positives and false negatives) can compensate each other without impacting the final results. Therefore, general accuracy is higher than the accuracy based on the confusion matrix. If the goal is to obtain an estimate of seismic vulnerability of a given region, rather than the vulnerability of a specific building, this approach is deemed acceptable. It is recalled that both accuracy evaluations can only be performed on the parts of the city where the real typology classification is available.

#### Seismic-Vulnerability Assessment (STEP 4)

The typological-class attribution (step 3) is the starting point for seismic-vulnerability assessment at urban scale. As introduced in section Introduction, the vulnerability of each type is defined TABLE 3 | Example of a St. Alban typology probability matrix (30% of the visual inspection | 213 buildings) with three attributes (number of stories, construction period, and roof shape).


TABLE 4 | Example of a confusion matrix obtained by Pmax-ARL method on the St. Alban district, considering 30% of the visual survey as training set.


*The association rule with three attributes (number of stories, construction period and roof shape) is subsequently applied to all the buildings that have been visually surveyed. Bold values are on the diagonal and provide the number of buildings that are correctly assigned by the method.*


TABLE 5 | "Real" distribution of buildings in St. Alban (obtained from the visually surveyed buildings) and the distribution given by ARL with the Pmax method.

*Bold values are the difference in absolute value between the real distribution and the Pmax-ARL distribution.*

according to appropriate vulnerability functions. Such functions are described by Equation (4):

> <sup>µ</sup>**d**=**2.5 <sup>1</sup>** <sup>+</sup> **tanh I** + **6.25V** − **13.1 Q** (4)

Where I is the seismic input provided in terms of macro-seismic intensity according to the European scale EMS-98; V is the vulnerability index according to the value of the vulnerability indexes for each type of construction; Q is the ductility index of the structure; µ<sup>D</sup> is the mean damage value that the building is expected to sustain.

Each typological class is linked to specific vulnerability and ductility indexes (V and Q), following previous studies (Lagomarsino and Giovinazzi, 2006). Vulnerability of buildings of the same typological class changes with the height category they belong to. The macro-seismic model (LM1) is used to perform vulnerability assessments. The mean damage that is expected for each class is determined using Equation (4). Thus, the accuracy achieved with the ARL method is assessed with respect to the predicted vulnerability. Classes, such as type M3 and M5, which are hard to distinguish in classification have similar vulnerability indexes, an improvement in accuracy can be expected with respect to typological-class attribution.

Starting from the mean damage (see Equation 4) obtained for each type and height category (related to the number of floors), damage distributions can be obtained. The probability, p<sup>k</sup> , related to each damage grade D<sup>k</sup> (k = 0 · · · 5), for a given mean damage µD, is evaluated according to the probability mass function of the binomial distribution (Lagomarsino and Giovinazzi, 2006):

$$p\_k = \frac{5!}{k!\left(5-k\right)!} \left(\frac{\mu\_d}{5}\right)^k \left(1 - \frac{\mu\_d}{5}\right)^{5-k} \tag{5}$$

Vulnerability can be assessed in two ways. The first one involves a deterministic evaluation: it is only based on the mean damage (µD) (Equation 4) for each class and height category. The second one is based for each class and height category on the distribution of damage grades obtained according to a probability mass function (using Equation 5).

#### Validation of Vulnerability Assessment

For both methods of vulnerability determination (deterministic and probabilistic), the accuracy is calculated as the sum of the differences in absolute value between the number of buildings attributed by the ARL to each damage grade and the "real" number of building for each damage grade (Diana et al., 2018), according to the following equation:

$$accuracy\_{vuln.} = 1 - \frac{\sum\_{k=0}^{5} |n^{\circ} 
u 
u 
u 
d.k\_{kreal} - n^{\circ} 
u 
u 
u 
d.k\_{katter.}|}{
to \text{t. } 
u 
u 
u 
d.} \quad (6)$$

In a similar way to typology attribution (see section Validation of Typology Attribution), the validation of vulnerability assessments is performed exclusively on the parts of the city for which real typological classes are available following a visual inspection campaign.

#### CASE STUDY

The objective of this paper is related to the seismic-vulnerability assessment of the city of Basel. Seismic hazard is not uniform throughout Switzerland. Some areas are characterized by higher hazard, such as the Basel region in northwestern Switzerland. In particular, the region of the city of Basel is classified in Zone Z3a, where the design value of horizontal acceleration of the ground is <sup>a</sup>gd <sup>=</sup> 1.3 - m/s 2 (SIA, 2014). The Basel region is characterized by the second-highest seismic hazard in Switzerland after the Valais region (Z3b), according to Appendix F of Swiss codes SIA 261 (SIA, 2014).

In 1356, the Basel region suffered the strongest earthquake ever recorded in the North of the Alps (Lestuzzi and Badoux, 2013). The Swiss Seismological Service (SED) estimates that this earthquake was of magnitude 6.6 on the Richter scale (Swiss Seismological Service, 2016). The return period of this earthquake is larger than the actual one of building codes (475 years) and is estimated to exceed 2000 years. In the present work, the Basel region is exposed to a macro-seismic intensity of IX according to the EMS-98, which corresponds to the shaking level of the 1,356 earthquake.

## The Collection of Building Attributes (STEP 1)

Three out of four attributes listed in methodology step 1 [section Collection of Building Attributes (STEP 1)] are collected in the Federal Register of Buildings and Dwelling (RegBL): the period of construction, the number of stories, and the footprint surface. The shape of the roof is collected in the 3D Model of the city of Basel. All these attributes are collected for the whole building stock of the city of Basel, composed of almost 21,000 buildings.

Concerning the number of floors, when performing the visual inspection phase that is essential to step 2, a slight difference between the attribute listed in the RegBL database and the data retrieved on site has been pointed out. For some buildings (32% of the buildings visually inspected) the number of stories listed in the RegBL database does not correspond to the number of stories derived from visual survey. These differences result from underground floors being included in the floor count of the Federal Statistical Office (BFS/OFS). In addition, "Attics and basements are only taken into account if they are designed, even partially, for housing. However, cellars are not taken into consideration" (Federal Statistical Office, 2017, p. 67). As the number of stories (seismic masses) is defined in groups (less than 4 stories, 4–6 stories, equal or more than 7 stories), errors are reduced. Only 8.2 percent of buildings composing the learning set [see section The Learning Phase (STEP 2)] are classified by the RegBL database in a different height category than by on-site visual inspections.

In order to derive the shape of the roof, data are collected from a vectorial 3D model of the entire city, which contains building-by-building features. In order to be able to process roof attribution in an algorithm, data are exported with a specific code. For each building, a mesh polyface of the building coordinates is extracted as a string of coordinates representing triangles of the surfaces.

An algorithm detects which triangles form the threedimensional surface of the roof. Then, the inclination of the roof is calculated as the angle between the normal vector of the triangle surface and the vertical. In case of multiple surfaces, if horizontal triangles represent more than 40% of the roof surface, then the roof is considered as flat. If not, the roof is considered as sloped. For sloped roofs, the highest slope of all surfaces is taken to define the inclination.

# The Learning Phase (STEP 2)

The learning set for the city of Basel is composed of 710 buildings, which are all located in the St. Alban district. The St. Alban district is chosen to be the learning set as it is deemed representative of the general distribution between the typological classes for the entire city of Basel. Indeed, the St. Alban district presents a miscellaneous building stock: buildings with different periods of constructions, heights (number of floors), materials and construction techniques (in aggregates or isolated), compose this district. The district is composed: for the 26% of simple stone masonry buildings (M3), for the 38% of unreinforced masonry brick buildings with flexible floors (M5) and for the 17% of masonry brick buildings with r.c. floors (M6). The remaining 17% is composed of r.c. buildings (RC1 + RC2 + RC3), mainly with resisting walls (RC2). As a consequence, applying ARL method on this district reduces the risk to bias the evaluation. Therefore, the typology correlations that are obtained for St. Alban can be expanded at city scale.

In addition to the four attributes collected in the public databases (step 1), the typological class for each building is defined for the 710 buildings composing the learning set of St. Alban. The attribution of typological classes to buildings is performed by analyzing archive drawings (on 30 buildings) and by performing either on-site or remote visual inspections. In total, among the 710 buildings, 53% have been surveyed online (using google maps or street view), 43% by on-site visual one-by-one inspection, and 4% based on archive drawings. The determination of the typological class for these buildings, following the classes defined in the Risk-UE project (**Table 1**), has a high certainty as the error in class attribution reduces with increasing amount of information.

To enable the successive validation phase, typological correlations between attributes have been performed on a subset of the learning set. This training set is composed of randomly selected 30% of the buildings composing the learning set (710<sup>∗</sup> 0.3 = 213 buildings) while the remaining 70% is considered as the validation set (710<sup>∗</sup> 0.7 = 497 buildings). As an example, **Table 3** shows a possible ARL matrix obtained from the combination of three attributes (number of stories, construction period and roof shape).

# Wide-Spreading the Typological Association to the Entire City of Basel (STEP 3)

The distribution matrix is the starting point to attribute typological classes using the ARL methodology. For the association of a specific typological class to each building, the Pmax and the Prandom methods are considered. For example, when considering the Pmax method attribution, all buildings built in 1929 with six floors and a sloped roof are classified as M5 buildings (see tenth row of **Table 3**). In the case of the Prandom attribution, a building built in 1929 with six floors and a sloped roof with a random instance equal to 0.93 will be classified as M6 (**Figure 2**). In the case when a specific typological class is not associated to a single building, the Pdistribution method is considered. In this case, for example, a building built in 1929 with six floors and a sloped roof will have the probability distribution reported in **Table 3** (tenth row) to belong to one of the classes M3 (11%), M4 (2%), M5 (76%), or M6 (11%).

The distribution matrix (**Table 3**) is applied to the entire city of Basel (almost 21,000 buildings). **Figure 3** shows that the three typological-class distributions (following the three attribution methods Pmax, Prandom and Pdistribution) yield similar predictions when the whole building stock of Basel is considered. Using the Pmax method, classes with low probabilities are ignored unlike for the Prandom and Pdistribution methods. Thus, classes M4, RC1, RC3, and S are not attributed to any building when the Pmax method is used. Prandom and Pdistribution distributions are highly similar, which is expected for large sets of buildings, such as an entire city.

# Seismic-Vulnerability Assessment for the Entire City of Basel (STEP 4)

With the typology distribution being defined in the previous section, the damage distribution for the entire city of Basel can be estimated for any seismic intensity. The damage distribution is achieved by the aggregation of the probability distribution multiplied by the number of related buildings according to the typology distribution (Diana et al., 2018).

Damage distributions for the city of Basel can be obtained for each ARL attribution method (Pmax, Prandom, Pdistribution). The three damage distributions derived from the three ARL attribution methods are similar. As the city of Basel has almost 21,000 buildings, the damage distribution for the Prandom and

FIGURE 2 | Example of the determination of the building class according to Prandom method for buildings built between 1900 and 1944, with 6 floors and a sloped roof. For random instances between 0 and 0.11, a building type M3 is attributed, for random instances between 0.11 and 0.13 type M4, between 0.13 and 0.89 M5 and over 0.89 type M6.

Pdistribution are identical. Moreover, only slight differences can be observed for damage distribution obtained using the Pmax-ARL.

Damage distribution, based on a probabilistic approach, for the city of Basel is shown in **Figure 4** for an earthquake of intensity I = IX. In the case of **Figure 4**, where buildings are assigned based on three attributes (construction period, number of stories and roof shape), no difference in the damage distribution is observed.

The distribution of damage can be calculated for each zip code (ZIP) of the city, in this way the vulnerability of districts

of the city can be evaluated. The Pmax-ARL method considering three attributes (construction period, number of stories and roof shape) is used for the whole Basel building stock.

**Figure 5** gives the distribution of damage for each region delimited by a shared ZIP code of the city of Basel after an earthquake of Intensity IX. This graph shows the vulnerability of specific districts in Basel. The most vulnerable district corresponds to ZIP 001, which is the historical city center. In the historic center, 68% of buildings are expected to sustain heavy damage or worse (DG3, DG4, and DG5) after an earthquake of intensity IX. This statement confirms the findings after the L'Aquila earthquake in 2009 about particular vulnerability of city-centers by Guéguen (2013).

The least vulnerable district corresponds to the zip code 002: the Universitätspital Basel (Basel university hospital). After an earthquake of intensity IX, there will be 44.6% of buildings with at least heavy damage (DG3, DG4, and DG5). If the Basel university hospital is ignored, the least vulnerable district has ZIP 011, which represents the suburban area. Fifty-six percentage of buildings in this area sustain or exceed heavy damage (DG3, DG4, and DG5) after an earthquake of Intensity IX.

# Evaluation of the Typological-Attribution Accuracy on the Learning Set

As discussed in section Validation of Typology Attribution, an important step is validation of the typological-class attribution based on the ARL. In **Table 4**, the confusion matrix obtained with the direct attribution of typological classes to each building (in this specific case with Pmax) is displayed. The confusion matrix is the mean for evaluate the accuracy of the attribution of types provided by direct attribution methods (Pmax and Prandom). In columns the "real" classes attributed by visual inspection are displayed while in rows the classes associated by ARL method. Therefore, on the diagonal the number of buildings that are correctly assigned are highlighted in bold while on the positions outside the diagonal the errors in attribution are displayed. The buildings on the last row (θ) of **Table 4** correspond to buildings that cannot attributed by the ARL method. For such buildings, there exists no building with the same attribute combinations within the buildings that form the training set. For **Table 4**, accuracy, in accordance with Equation (2), is equal to: 524/710 = 0.738.

The confusion matrix also allows to assess the performance of typological attribution ofspecific types. For instance, 74 buildings are wrongly attributed to be class M5 and thus, are false positives (precision of 261/(261+74) = 0.78). Also, 25 buildings of type M5 are false negatives as they are not recognized by the ARL method (recall of 261/(261+25) = 0.91). The general score of accuracy combines both types of errors: false positives and false negatives.

In **Table 5**, the distribution obtained applying the ARL method and the "real" distribution obtained considering visual inspections are displayed. The distribution evaluation is admissible for all the three attribution methods considered (Pmax, Prandom, and Pdistribution). It is an evaluation that is appropriate for large-scale evaluations since it considers exclusively the total number of items being part of a specific class, without considering building-by-building misclassifications. For the example of **Table 5**, accuracy, calculated using Equation (3), is equal to: 1–115/710 = 0.838. It is worth noting that in distribution evaluation, the accuracy is always higher since errors in class attributions compensate each other.

One of the main goals of the paper is to evaluate the accuracy of building-class attributions that is achieved using several attribute combinations. As stated in section Collection of Building Attributes (STEP 1), considered attributes are: the period of construction; the number of stories; the shape of the roof; and the surface footprint of the building. These four attributes are combined in multiple ways and possible improvements in the classattribution accuracy are evaluated. Combinations of two attributes (period of construction and number of stories), three attributes (adding either roof shape or surface footprint) and all four attributes are evaluated and compared in **Table 6**.

When using the confusion matrix as metric of accuracy, maximum accuracy is obtained by combining three attributes (construction period, number of stories and shape of the roof) and using the Pmax-ARL attribution method. The achieved accuracy is 73.9%. The Pmax method provides more accurate results than the Prandom method as each building is assigned the typological class of maximal probability. More often than not, this assignment is correct.

As can be seen in **Table 6** (first two rows), considering more attributes can result in decreasing attribution accuracy. For example, with three attributes (construction date, number of stories and footprint surface), the accuracy of the Pmax-ARL method drops to 68.9% (−2.1%) with respect to the accuracy that is obtained by considering two attributes. When considering four attributes, the accuracy is lower than the accuracy obtained with three attributes (72.1% as opposed to 73.9%). This results from overfitting the training set: the distribution matrix matches the particular data of the training set too closely (with four attributes) and thus, when wide spreading it to a larger set, the lower general validity reduces attribution accuracy. Indeed, random fluctuations in the training data is learned as a correlation between attributes and types. These fluctuations do not apply to new data (validation set) and thus, negatively impact the models ability to generalize.

When assessing classification accuracy using typological distribution (see **Table 3**), the best accuracy is obtained when three attributes (construction periods, number of stories and shape of the roof) are combined using the Pdistribution-ARL method. The maximum accuracy is 97.2% (see **Table 6**, last three rows). Unlike when accuracy is assessed using a confusion matrix, the classification obtained based on the Pdistribution (and the Prandom) methods are more accurate than Pmax. For


TABLE 6 | Accuracy and change in accuracy for confusion matrix and distribution of classes.

*Bold values are the most accurate results.*

Pdistribution and Prandom attribution methods, the classification follows probabilistic distributions, which explains why they outperform Pmax for which only the class with maximum likelihood is attributed. In general, the accuracy based on the typology distribution is higher than that based on the confusion matrix, since possible errors can be compensated by other errors at larger scale. From **Table 6** it can be seen that the shape of the roof seems to be particularly influent, improving substantially the distribution accuracy.

In **Figure 6** the number of buildings that are wrongly classified in the case of the two attributes approach (number of stories and period of construction) are represented. On the left of **Figure 6**, the "real" types for the misclassified buildings are shown, while on the right the types given by the Pmax-ARL method are represented. Thus, **Figure 6** can help in understanding the tendencies in mis-attribution of typological classes: as it is possible to notice, errors occur mainly for buildings with four to six floors. In this height interval, for the "before 1900" period, the Pmax-ARL method does not assign any building into the M5 class (bricks masonry with flexible floors), that is in the visually surveyed distribution the most common class. Between 1900 and 1969, the Pmax-ARL method assigns buildings exclusively to class M5 while in the "real" distribution M3- (stone masonry), M4- (massive stone), M6- (brick masonry with rigid floors), and RC2-classes (reinforced concrete wall building) are present.

As stated before, the accuracy of class attribution is improved by introducing the shape of the roof. A limit of 5.7◦ has been chosen as the inclination that separates flat roofs from sloped roofs. **Figure 7** shows the distribution of roof shapes for each class that has been visually surveyed in St. Alban. This graph shows that most M3, M4, and M5 buildings have a sloped roof while most reinforced-concrete (RC) and steel buildings (S) have flat roofs. For M6 buildings, 50% have a flat roof and the other 50% have a sloped roof. The shape of the roof seems to be closely related with the typological class, as shown by the improvement in terms of general accuracy in **Table 6**.

In order to improve the classification accuracy, addition of a third interval for roof slopes, has also been tested. The limit of this angle interval is set as an unknown variable first and the optimal value is derived as the value that maximizes the accuracy. The introduction of this limit angle on classification accuracy is irrelevant. No improvement in accuracy (neither based on the confusion matrix nor on the typological distribution) has been found. Therefore, taking into account two inclination intervals, either flat or sloped, is sufficient.

# Evaluation of the Accuracy of Damage Distributions

The accuracy of damage distributions are summarized in **Table 7**. The best accuracy is obtained when three attributes (construction periods, number of stories and roof shape) are combined and the Prandom-ARL method is used for typological-class attribution (98.4%). The Prandom-ARL method is not stable, as it involves random numbers, but studied random instances provide better results than the Pdistribution-ARL method. For larger amounts of buildings, the accuracy of the Prandom-ARL method equals the accuracy of the Pdistribution-ARL method. Without considering the Prandom-ARL method, the best accuracy is obtained with four attributes (construction periods, number of stories, shape of the roof and footprint surface) with the Pmax-ARL method. The accuracy is 96.3%, which also consists the highest improvement compared with the two-attribute results (+ 7.5%).

Regarding the accuracy of probabilistic vulnerability determination (obtained using Equation 5), the best accuracy is obtained using three attributes (construction periods, number of stories and roof shape) with the Prandom-ARL method (99.8%). Again, this method does not provide stable results and for a large amounts of buildings, accuracy of the Prandom-ARL and the Pdistribution-ARL methods, will yield similar results.

When the Prandom-ARL method is ignored, the best accuracy is obtained with three attributes (construction period, number of stories and shape of the roof) and with the Pmax-ARL method. The accuracy is 99.4%.

#### Evaluation of the Accuracy for Another District

The accuracy of the ARL method in attributing typological classes is also performed on 223 buildings that have been surveyed in Iselin, another part of the city of Basel. **Table 8** (first part) summarizes the accuracy based on the confusion matrix with the improvement of considering three and four attributes. The best accuracy (based on the confusion matrix) is obtained with four attributes (construction periods, number of stories, roof shape and footprint surface) with the Pmax-ARL method. The accuracy is 70.4%.

Considering additional attributes does not always provide better accuracy, as can be seen in **Table 8** first row. For instance, with three attributes (construction date, number of stories and footprint surface), the accuracy of the Pmax-ARL method is lower than with two attributes (−3.3%). For the Pmax-ARL method, considering four attributes (construction periods, number of stories, roof shape and footprint surface), instead of two, improves the accuracy by 5.4%.

**Table 8** (row three to five) summarizes the accuracy in the distribution of types with the improvement of considering other attributes as the roof shape or the footprint surface. The best accuracy (92%) of the typology distribution is achieved using all four attributes with the Pdistribution-ARL method. The Pdistribution-ARL method, considering four attributes (construction periods, number of stories, roof shape and footprint surface), instead of two attributes, improves accuracy by 4.4%.

**Table 8** (first three rows of the second part) summarizes the accuracy for the damage distribution based on the mean damage with the improvement of considering other attributes. The best accuracy for the distribution of mean damage is obtained with four attributes (construction periods, number of stories, roof shape and footprint surface) with the Pmax-ARL method. The accuracy is 96.1%. For the Pmax-ARL method, considering four attributes (construction periods, number of stories, roof shape and footprint surface), instead of only two improves the accuracy for the distribution of mean damage of 9.1%.

The accuracy for the damage distribution based on the probabilistic distribution is summarized in **Table 8** (last three rows). The accuracy of damage distribution (with probabilistic damage) is the best for the Pdistribution-ARL method considering all four attributes (construction periods, number of stories, roof shape and footprint surface). The accuracy of the blind check is 98.8%. For the Pmax-ARL method, the accuracy of the damage distribution is higher when only three attributes (construction periods, number of stories and the roof shape) are considered. This accuracy equals 98.6%. While considering three attributes (construction periods, number of stories and roof shape) instead of two improves the accuracy by almost 6%, considering all four attributes instead of three increases accuracy by 5.4%.

## SUMMARY AND CONCLUSIONS

Seismic vulnerability at large scale is time consuming and thus, simplifications are needed. Although a taxonomy should be initially defined to cover the typological classes that compose the building stock, there is potential to speed up the attribution of typological classes to each


TABLE 7 | Accuracy and change in accuracy for the seismic vulnerability assessment based on the deterministic and probabilistic approach.

*Bold values are the most accurate results.*

building using few attributes that are available in existing data-bases, rather than performing visual surveys for each building. In this paper, a typological class attribution based on the association rule learning (ARL) is proposed, based on various combinations of the following attributes: construction periods, number of stories, shape of the roof and footprint surface.

Through a case study, in which the seismic vulnerability is assessed for the entire city of Basel, located in the region with second-highest seismic hazard of Switzerland, the accuracy of a large-scale ARL-based typological-class attribution using attributes that are available in existing databases is assessed. Typological classes of several hundreds of buildings have been derived from visual surveys in a selected district and are used as training set in order to derive correlations between attributes and typological classes. Buildings of another district of the city, for which typological classes have been visually derived, are used as validation tests to evaluate the accuracy of the attribution method.

The following conclusions are drawn:

• in general, for large scale assessments, the loss in accuracy in ARL-based typology attribution is irrelevant when considering seismic damage predictions. Therefore, after defining a rather small learning set (around 5% of the entire building stock of a town), building classes can be assigned to all buildings by considering only three attributes that are readily available in


*Bold values are the most accurate results.*

public databases (i.e., period of construction, number of floors, shape of the roof). The reduction in terms of time demand for the preparation of reliable seismic vulnerability scenarios at city scale is particularly pronounced and only slightly reduces accuracy (<5% of error).


## AUTHOR CONTRIBUTIONS

LD: general organization of the research. JT: elaboration of the case study. YR: generation of distribution matrixes. PL: general overview. All authors proofread the paper and agree with the content.

## ACKNOWLEDGMENTS

This work was included in a more comprehensive project framework addressing seismic risk analysis for the city of Basel and was partially financed by the Canton Basel-Stadt.

## REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Diana, Thiriot, Reuland and Lestuzzi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

digital media

of impactful research

article's readership