Probabilistic Incremental Dynamic Analysis for Seismic Isolation Systems through Integration with the NHERI-SimCenter Performance-Based Engineering Application

: In the wake of the 1994 Northridge and 1995 Kobe earthquakes, structural designers adopted performance-based engineering concepts instead of traditional deterministic design approaches. The primary change was to evaluate the design according to stakeholders’ interests instead of the engineering parameters. This evaluation process required a probabilistic estimation for the included variables at all design stages. The NHERI-SimCenter application framework provides generic solutions implemented in di ﬀ erent hazard simulation problems. Seismic isolation is an e ﬃ - cient, proven technique for improving seismic performance by limiting drift ratios and reducing damage. During the design stage, seismic isolation-system parameters must be carefully calibrated to control di ﬀ erent aspects of the response, and it is necessary to run multiple simulations and count for parameter uncertainties. This research introduced components for seismic isolation and incremental dynamic analysis procedures integrated with framework modules, such as uncertainty quanti ﬁ cations and damage and loss estimation. Furthermore, an illustrative case study was included to re ﬂ ect the impacts of this development.


Introduction
The traditional philosophy of seismic design of structures is bounded by two primary variables: earthquake intensity and corresponding damage for structural and non-structural elements. The acceptance criteria are: preventing any damage in low-intensity earthquakes; accepting only repairable damage in medium-intensity earthquakes; and preventing collapse in high-intensity earthquakes. The structural engineering community recognized many deficits in this methodology after the 1994 Northridge and 1995 Kobe earthquakes. Although structures were designed according to seismic codes inherited from this philosophy, the reported total damage and economic losses were extremely high [1]. The deterministic methodology describes the structural design only in binary terms, as either safe or unsafe. Additionally, it depends only on meaningful parameters for engineers, not other stakeholders. These facts initiated the concept of performance-based earthquake engineering (PBEE).
In 1995, The Structural Engineers Association of California (SEAOC) published its PBEE Vision 2000 report, which is considered the most remarkable document of first- The Pacific Earthquake Engineering Research Center (PEER) developed a more robust methodology fitting the complex, multi-disciplinary problem, thereby overcoming the deficits of earlier PBEE methods. Unlike first-generation PBEE frameworks, estimating the system performance includes variables that reflect the direct interests of various stakeholders, such as monetary loss, repair duration, and casualties. System performance assessment is based on a rigorous probabilistic framework consisting of four logical, sequential components: hazard analysis, structural analysis, damage analysis, and loss analysis. Each component defines and quantifies a set of variables according to two inputs probability of exceedance (POE) and inherent uncertainties. The relationships among element outcomes are expressed in a probabilistic form, i.e., ( | ). These groups of variables are respectively known as intensity measures (IM), engineering demand parameters (EDP), damage measures (DM), and decision variables (DV), as shown in Figure 2. Many publications have summarized and explained the PEER methodology, such as [4][5][6].
The selection of optimal intensity measures (IMs) in probabilistic seismic demand models depends on several criteria, including efficiency, sufficiency, practicality, effectiveness, and hazard computability [7,8]. The efficiency criterion requires that the selected IMs be capable of providing accurate and reliable predictions of ground motion intensity. The sufficiency criterion requires that the selected IMs capture the key aspects of ground motion that are relevant to the response of the structure or system under consideration. The practicality criterion requires that the selected IMs be measurable and readily available. The effectiveness criterion requires that the selected IMs be capable of capturing the effects of earthquake hazards on the structure or system. Finally, the hazard computability criterion requires that the selected IMs be computationally efficient and amenable to probabilistic hazard analysis. In summary, the selection of optimal IMs for probabilistic seismic demand models involves a careful balance of these criteria to ensure that the models are both accurate and practical for use in seismic risk assessment and engineering design.
The hazard analysis considers the sources of earthquakes and facility characteristics to evaluate the ground motions' IMs. It requires selecting nearby faults and collecting their descriptions, i.e., source-site distance, magnitude-recurrence frequencies, and fault mechanism. The facility data are defined in two terms: the location (O) and design characteristics (D), i.e., the fundamental period of vibration and the foundation type. Engineering seismology relationships, such as ground motion prediction equations, are employed to quantify the desired variables from these inputs, including issued uncertainties. The intensity measure variables might include: peak ground acceleration (PGA); peak ground velocity (PGV); spectral acceleration at the period of the first mode (Sa(T1)); spectral shape; and the duration of ground motions. The current attenuation relationships are based on interpolating empirical data from past earthquakes. For example, the PEER Strong Motion Database (http://peer.berkeley.edu, accessed on 5 April 2023) contains over 1500 records collected from 143 earthquakes [9]. Moreover, it contains associated site conditions and rupture properties. Although this estimating procedure is industry-accepted, it can cause a high bias error and subsequential uncertainty in extensive scattered data cases. PEER is currently investigating the expansion of this database by adding new records or providing more information about the included records. The next generation will include a higher range of data variety, less dispersion, and more confidence in expecting shaking intensity. The hazard curve shows hazard analysis results by plotting the selected IMs versus their mean annual frequency (M.A.F.) of exceedance [10]. These results can be transformed into IM versus probability of exceedance (POE) in "t" years by using Equation (1) for the Poisson model, as shown in Figure 2: where ( ) is the annual frequency of exceedance of IM, and t can be considered as the duration of the facility life cycle.
After quantifying IM variables, the next step is to perform structural analysis to quantify the structure's responses to these various ground motions. This process involves developing a computational model, performing nonlinear simulations, and calculating engineering demand parameters (EDP). This step's uncertainty stems from the numerical model parameters such as mass, stiffness, strength, etc. Generally, EDPs characterize the response in terms of deformations, accelerations, induced forces, or other appropriate quantities [11]. Additionally, EDPs can be determined for local structural components' behavior or the overall building behavior. For a structural component, element forces (i.e., axial and shear forces) and plastic deformations (i.e., rotations and deflections) are examples of practical EDPs. The global building behavior can be described in EDPs such as inter-story drift, floor velocity, floor acceleration, roof drift, etc. Global EDPs can indicate damage for different damageable components (both structural and non-structural elements). For example, inter-story drifts measure structural system damage [12], while floor acceleration is used for non-structural equipment, e.g., office or laboratory buildings [13]. PEER developed the incremental dynamic analysis (IDA) procedure to investigate the relationship between ground motion intensity and structural simulations. This procedure is based on repeating structural simulations with scaling input ground motions. The resulting series of values can establish the statistical relationships between IM and EDP or the POE of EDP at a particular IM. The structural analysis step requires several nonlinear time history analyses, as the total number of simulations equals the number of ground motions multiplied by different modeling variations. The PEER-PBEE formulation adopts a single value for each defined EDP, so the peak value from all simulations is selected. Additionally, performing multiple simulations demands excessive effort from structural engineers and consumes high computational time for traditional software. However, PEER has overcome this barrier by developing OpenSees (Open System for Earthquake Engineering Simulation; http://opensees.berkeley.edu, accessed on 5 April 2023) [14]. It represents an open-source, object-oriented, and parameter-based finite element simulation software. It contains several commands to model materials, sections, joints, structural elements, and liquefiable soils in a nonlinear behavior. Additionally, integrating structures and soil models allows the integration of soil-structure-foundation-interaction (SSFI) investigations. Unlike traditional methods, OpenSees can handle modeling parameter variations and multiple ground motions with less computational cost. Its efficiency can be attributed to the software architecture of well-designed interfaces by a substantial community of developers. The results of the structural analysis can be summarized as the determination of probability distribution functions (PDFs) of EDP (e.g., EDPj) for a given IM (e.g., IMm), as shown in Figure 2. The number of PDFs equals α multiplied by β, where α is the number of IM data points and β is the number of considered EDPs [11]. The damage measures (DM) are the quantitative variables describing the physical damage to structural elements, non-structural elements, and contents. These measures can be related to given EDPs through the damage analysis conditional probability relationships. The conditional relationships can be calibrated by the results of experiments or previous earthquake reconnaissance reports. PEER continuously develops conditional probability relationships for different components by collecting previous test data or performing other experiments. This procedure requires defining damageable groups, as each group includes facility parts that are similarly affected by certain EDPs. For example, Bohl [15] proposed a methodology for damage analysis of a steel moment frame building. It considered 16 damageable groups: the structural system, exterior enclosure, drift-sensitive and acceleration-sensitive non-structural elements, and office content on each floor. Multiple DM levels describing the types of required repair for a damageable group should be defined. For example, the three DM levels of structural elements defined by [16] are light, moderate, and severe (or collapse). They respectively represent repair with epoxy injections, repair with jacketing, and element replacement. A single EDP value cannot be mapped deterministically into a certain DM level because this value can occur in different structural response scenarios and with different behaviors. Instead, any DM level can occur according to this value but with different probabilities. Fragility curves are established to show the variations in POE of all defined DMs in different EDP values, as shown in Figure 2. The probability of kth DM level can be estimated from fragility curves according to the following algorithm for the jth damageable group of the facility, and the ith value of the EDP utilized for the fragility function: for k 1:# of DM levels The last step of the PEER PBEE methodology is the loss analysis to estimate decision variables (DVs). These variables reflect the direct interests of various stakeholders in the design. Some of the DVs are economic loss, repair duration, and the number of fatalities and injuries. A loss function is determined to calculate the POE of the losses for different damageable groups at different DMs. The total number of loss functions required for a facility is γ multiplied by λ, where γ and λ are the number of DMs and damageable groups, respectively. The accumulative loss curve can be calculated based on the total probability theorem. The probability of the nth value of the DV is according to Equation (3) [17] in the form of a triple summation.
The structural engineering community started to realize the impact of PBEE methodology. It was included in several benchmark studies, such as [12,13,18].

Research Significance
The research presented in this study holds significant importance in the field of structural engineering and seismic performance. The adoption of performance-based engineering concepts following major earthquakes has revolutionized the traditional deterministic design approach. By focusing on stakeholders' interests and incorporating probabilistic estimation in the design process, engineers can enhance the reliability and safety of structures. Seismic isolation, a proven technique for mitigating earthquake damage, plays a crucial role in improving the seismic performance of structures. However, accurate calibration of seismic isolation system parameters is essential to control various response aspects. This research addresses this critical requirement by introducing components for seismic isolation and incremental dynamic analysis procedures integrated with uncertainty quantification and damage and loss estimation modules within the NHERI-SimCenter application framework. The addition of new features, including solving different isolated models and updating the user interface, enhances the functionality and usability of the framework. By advancing the capabilities of the NHERI-SimCenter application framework, this research contributes to the development of more reliable and efficient seismic design practices, ultimately leading to improved structural resilience and reduced seismic risks.

NHERI-SimCenter Application Framework
The Computational Modeling and Simulation Center (SimCenter) [139] develops many software tools specialized for natural hazard engineering (NHE) problems. The SimCenter's vision contributes to The Natural Hazards Engineering Research Infrastructure (NHERI) program, funded by the National Science Foundation (NSF) in the United States. The main goal is to provide an application framework for researchers to use, build, and extend scientific workflows for natural hazard simulations. This framework implements integrated relationships among three units: model, scientific workflow, and application. The model is a conceptual component with a particular role or analysis procedure (refer to Figure 3). Practically, many software applications can simulate each model. However, a workflow is a collection of models integrated into a sequence to automate multistep solutions. This approach requires developing pre-and post-processors for existing applications for data parsing. The end-user can customize the workflow by selecting the in-use application for each model depending on the nature of the problem, available input format, and the required specifications. Finally, the NHERI-SimCenter application is defined by a workflow as a backend package connected with a wrapping frontend user interface (UI). The UI components define workflow applications and their parameters, generate input file(s), and present output results. The architecture of the NHERI-SimCenter application framework [140], shown in Figure 4, has an intelligent design because it considers the following principles: • The architecture is designed as a modular framework to integrate with external applications or datasets. This feature can save repeating the work of an existing application and guarantee futural scalability and interoperability. For example, the NHERI-SimCenter team developed pre-and post-processors to benefit from existing applications such as OpenSees [141], OpenFOAM [142], and PEER Strong Ground Motion Databases [143]. Moreover, they developed and linked additional applications with the framework, such as Building Recognition using AI at Large-Scale (BRAIL) [144], to automate the collection of building inventory data; • End-users can excuse the computational workflow in their local desktop environment or on the cloud through DesignSafe [145]. The second method utilizes highperformance computing resources and parallel workflow executions for large-scale datasets; • The framework allows for defining input uncertainties in models. It implements uncertainty quantification (UQ) by utilizing Dakota [146] in sampling models with random variables. It supports various UQ methods: forward propagation, sensitivity analysis, and reliability analysis; • The architecture provides a high order of flexibility and generality. It is designed for NHE problems with different hazards, e.g., earthquake, wind, tsunami, or assets classification, such as buildings, network pipelines, bridges, etc. The ability to customize the scientific workflow, in areas such as available inputs, required outputs, and defined specifications in some models may accommodate the different interests of researchers.  PELICUN [147], the Probabilistic Estimation of Losses, Injuries, and Community Resilience Under Natural Disasters, is a performance assessment framework developed by the NHERI-SimCenter. It includes three different workflows to estimate DVs from input IM, as shown in Figure 5. Each workflow determines a specific path and assumptions as follows: 1. The first path is the most rapid and straightforward calculation approach, skipping the structural and damage analysis steps. It consists of a single direct step from IM to DVs using vulnerability functions. These functions are calibrated using ground motion intensity maps, and insurance claims data from past earthquakes. This approach may be efficient for limited asset types, i.e., single-family timber houses; 2. The second approach breaks the process into two steps: estimating DMs through fragility functions and using consequence functions to identify DVs of repair. Fragility functions are calibrated based on previous earthquake databases of IMs and corresponding DMs for different building types. Similarly, consequence functions can be determined to describe the cost and time of repairs as a function of damages from standard construction practices; 3. The third path is the most detailed and computationally demanding approach, as it adds structure analysis to the second path. It requires a sophisticated estimation of structure response through a finite element model or measured with a structural health monitoring system. This response is used to identify the values of EDPs, e.g., peak inter-story drift (PID), residual inter-story drift (R.I.D.), and peak floor acceleration (PFA). Then, calibrated fragility curves measure the damage variables as a function in the EDPs. Uncertainty quantification is a necessary model in the loss assessment process. Its implementation requires the creation of a large number of samples with random variables for inputs. The accuracy of the resulting DVs distribution curves increases according to the number of samples. However, increasing this number demands more computational resources and time for execution. This increase could be acceptable for the first and second paths and extreme for the third path because of the complexity of the structure response simulation. This problem is handled by estimating a few dozen EDP samples, then resampling them by fitting a probability distribution. The re-sampling technique will generate more samples for DMs and DVs models because they require fewer calculations with calibrated functions.
As explained in Figure 6, the software architecture has a flexible design that allows for all three paths [147]. Moreover, it considers minimizing computational resources and future scalability. The source code is written in Python to provide cross-platform usability, meaning that it works on Windows, Macintosh, and Unix operating systems. Furthermore, this will facilitate its usability in high-performance computing (HPC) cluster environments without requiring administrator support, compilation, or the installation of dependencies. Workflow inputs are the Asset Information Model (AIM) and EVENT files; both should follow a standard JSON format. Structure-response simulations are out of PELICUN's scope, so this model integrates an exterior application to the framework, such as OpenSees [141]. The process of integration requires only developing interface classes for parsing data formats. The architecture includes four main models, and their objectives are shown in Figure 7 and determined as follows [20]:

•
The response model describes the outputs of structural simulations according to types of EDPs written in the raw EDP file. The peak values are selected in terms of (EDP:t-l-d) where t is the EDP type, l refers to the location (the floor), and d represents the direction index according to a predefined list of directions. For example, (PID-2-1) describes the maximum value of inter-story drift that occurred on the second floor and towards the first axis; • The performance model organizes the asset components (structural, non-structural, and contents) in a descriptive hierarchy for damage analysis based on the FEMA P58 method [148]. The hierarchy consists of three layers: 1. The highest level is responsible for defining fragility groups (FGs). Each FG contains all the components whose damage behavior is controlled by the same EDP type and leads to similar consequences. According to the FEMA P58 method [148], buildings have FGs sensitive to either PID or PFA; 2. The intermediate layer divides F.G.s into several performance groups (PGs).
Components in the same PG are sensitive to the same EDP; their location and direction must be identical; 3. The third layer breaks each P.G. into different component groups (CGs   The implementation of damage and loss models handles both the high-fidelity component-based FEMA P58 method [148] and the more approximate building-level approach used by HAZUS [149]. The framework utilizes parallel executions by dividing calculations into two threads: (a) structure response simulations, creating the response model and estimating EDPs, and (b) assembling the performance, damage, and loss models. PELICUN is an essential component in the backend calculations of the NHERI-SimCenter application workflow, and many released applications are based on its capabilities.
The launched NHERI-SimCenter applications simulate performance-based engineering but with different problem definitions. These differences can be regarding desired workflow architecture, the type of hazard studied, or the investigation scale (individual assets or regional simulations). The launched application can be classified into several layers of hierarchy according to their complexity, refer to Figure 4. This hierarchy starts with the quoFEM [150] application, which calls for quantified uncertainty with optimization for the finite element method. This application provides a probabilistic structural analysis by integrating finite element modeling (FEM) with uncertainty quantification (UQ) procedures and simple UI. The FEM component includes several applications, e.g., OpenSees [14], OpenFOAM [142], OpenSeesPy [151], and FEAPpv, while the UQ model implements Dakota [146]. The second level contains three applications: EE-UQ [152], WE-UQ [153], and HYDRO-UQ [154], respectively, for earthquake engineering, wind engineering, and water-borne natural hazards with uncertainty quantification. In addition to using user-defined simulation scripts, this layer provides customized models to describe a building asset numerically. Also, it provides many options for inputs or to describe hazards, such as selecting records from a database, using previous records, or setting the parameters of other specifications. The outputs are the sampled EDPs determined by locations and directions. Moreover, applying performance-based engineering (PBE) [147] introduces a more comprehensive workflow than EE-UQ by combining damage and loss estimations using PELICUN components. The Regional Resilience Determination Tool (R2D) [155] presents the highest level, which was developed for simulating risk assessment on a regional scale on buildings inventory. It uses graphical components to classify the risk levels on a map view. This broad vision still has details and areas that need further development. The architecture of the framework is designed to launch more applications flexibly.

Seismic Isolation
According to conventional seismic specifications, structural design is based on collapse prevention with expected damages from strong earthquakes. This approach could not be acceptable for some facilities, such as hospitals, fire stations, nuclear plants, and historical buildings. These cases require a control mechanism on the structural response and reducing inter-stories drift and floor accelerations to decrease component damages and injuries. The seismic isolation technique was developed for this purpose, and the technique involves physically uncoupling the superstructure from the foundation [156]. The decoupling concept is based on isolation devices that provide a high stiffness in the vertical direction and flexibility in the horizontal direction. These properties dissipate motion energy and concentrate deformed displacement at the isolation interface. The fundamental natural period of isolated structures is in the long period range, e.g., two to four seconds, so the superstructure remains elastic or nearly elastic. In some cases, isolation systems contain dampers in cooperation with isolation devices to increase the capacity of the overall dissipated energy. Many publications compared isolation systems and design procedures, such as [157,158]. Practically, the increase in seismic isolation systems can prove the efficiency and impact of this technology. Various iconic and critical building designs implement seismic isolation to preserve serviceability even after large earthquakes [159], such as: (1) Apple's new corporate headquarters, which opened in April 2017, and is a ring-shaped building with a circumference of 1512 sq ft and houses 12,000 employees. It contains four stories above the ground and three stories below. Seven hundred bearing isolators support the superstructure. Each isolator has a 7 ft diameter and weighs about 15,000 lbs; and (2)  Isolation devices are designed to possess vertical stiffness for required service loads and lateral flexibility calibrated for damping and isolation demands. Figure 8 shows a traditional finite element model to simulate a bearing element. It consists of two primary nodes (i and j) connected by six springs in the six primary directions (12 degrees of freedom). In Opensees [141], most isolator-bearing element commands are an extension of this model. Predefined uniaxial materials, i.e., P, T, My, Mz, determine the general behavior of this element in the axial direction, torsional direction, moment around the local yaxis, and moment around the local z-axis. However, the uniaxial materials of shear springs are customized according to the mechanical properties of the isolator device to calibrate its horizontal flexibility. Additionally, the location of shear springs from the ith node is expressed as a relative ratio to the isolator height, called the shear distance ratio (sDratio). The isolator devices can be classified according to their mechanism into two main categories: friction-based and elastomer-based [161]. The concept of a friction-based isolator is to support the weight of structures on bearing support with a sliding interface [162]. The horizontal strength or zero-displacement force intercept ( ) is calculated according to the following equation: where W is the weight carried by the bearing isolator, and µ represents the friction coefficient of the sliding interface. Isolators should have a low value of µ to permit in-plane displacements and achieve flexible behavior. According to the manufacturers' reports, the practical value could range from 0.03 to 0.2. The slider is usually manufactured from ductile iron bonded with polytetrafluorethylene (PTFE) material. Figure 9 shows different friction-based slider-bearing devices and their nonlinear shear behavior. The friction effect can be simulated numerically on OpenSees [141] through friction models: (Coulomb), (velocity dependent), (velocity and normal force dependent), (velocity and pressure dependent), and (multi-linear velocity dependent). The flat slider bearing element is the simplest device with a flat steel-bearing plate and a slider that is allowed to move horizontally. This model leads to elastic-perfectly plastic behavior in the shear force and displacement relationship. In Opensees [141] an encapsulated element command, flatSliderBearing, was developed for this type of isolator. It considers the primary bearing element, shown in Figure  8. Shear springs are assigned to an elastic-perfectly plastic uniaxial material, whose behavior is a function of the elastic slope (Kinit) and an assigned friction model. Specifically, the sDratio is determined by the position of the slider interface. The single friction pendulum (FP) and the RJ-Watson EQS bearing elements are evolving forms of the flat slider element, as their behavior has a linear increase in the force capacity through the plastic zone. The restoring force capability in the single FP element results from exposing the slider to a spherical concave dish. The second-slope stiffness (Kd) is estimated according to Equation (5), where R is the radius of curvature of the spherical concave dish: Also, Opensees contains an element command for modeling a single F.P. element, named singleFPBearing, which follows the same assumptions of the flatSliderBearing model. However, the nonlinear sliding hinge movements are assigned to a spherical surface instead of a horizontal plane. The effective radius of the concave sliding surface (Reff) is required as an additional parameter to generate a uniaxial material with strain-hardening behavior. For the RJ-Watson EQS bearing elements, MER springs along the y-and zaxes provide restoring forces on the slider interface. Equation (6) demonstrates how to calculate Kd for this system, where k2 and k3 are post-yield stiffness of linear and nonlinear hardening in MER springs, dy is the displacement as shown in Figure 9, and η is the exponent of nonlinear hardening.
Similarly, the RJWatsonEqsBearing element command was developed following previously mentioned assumptions. In this case, k2, k3, and η are the additional parameters for shear springs material. Instead of a single concave dish in the single F.P. isolator, the triple friction pendulum bearing (TPB) elements have three spherical surfaces associated with an inner slider, the top plate, and the bottom plate, as shown in Figure 10. It is considered the most effective friction-based isolation system for intense earthquakes [164]. Since the normalized backbone curve is complicated, the TripleFrictionPendulum element command requires parameters to simulate this complexity perfectly. The essential prerequisite parameters are (1) three predefined friction models for the three sliding surfaces; (2) the displacements limits (d1, d2, and d3), which are clarified in Figure  10 The elastomer-based bearing consists of alternating natural or synthetic rubber layers bonded to intermediate steel shim plates and an outer rubber cover to protect steel layers from corrosion and environmental degradation, as shown in Figure 11. After placing un-vulcanized rubber sheets and steel plates in a mold, this mold is subjected to high temperature and pressure to vulcanize and bond the rubber. Increasing the bearing's vertical stiffness corresponds to the number of steel shims and closing their spacing. However, reducing the horizontal stiffness requires a relative increase in the total thickness of rubber (Tr). The design of elastomeric bearing elements is a compromising process among these factors to produce determined properties [165]. Elastomeric bearings can be classified into three categories: (1) low-damping natural or synthetic rubber; (2) high-damping rubber (HDR); and (3) lead-rubber bearings. Natural rubber with a type A durometer hardness of 50 is the most used low-damping material for seismic bearings. It has an approximate shear modulus (G) that ranges from 0.65 MPa to 0.9 MPa, and an equivalent damping ratio ξ between 2% and 3% at 100% shear strain. These properties can be modified by adding carbon black and other fillers to the raw rubber to create highdamping rubber (HDR) materials. The equivalent damping ratio of HDR bearings can reach from 10% to 20% at 100% shear strain, with an almost equivalent shear modulus. The horizontal stiffness of the bearing (Kh) can be calculated as the following, where Ab is the total bonded rubber area: The spacing between steel plates (the thickness of the individual rubber layer) defines the compression modulus of the individual elastomeric layer (Ec) due to changes in the bulging around the perimeter. According to [166,167] and assuming the rubber is incompressible, the Ec of a single rubber layer can be calculated as follows: where S is a dimensionless geometric parameter called the shape factor and defined as:

=
Loaded Area Area free to bulge (10) The simplification assumption of rubber incompressibility overestimates the compression (vertical) stiffness (Kv). According to [166], this overestimation increases in the cases of large S values. Typically, this assumption is invalid for thin rubber layers (close steel shims). The Kv is defined as: Inspection of Equations (8) and (11) reveals that the relative ratio of horizontal to vertical stiffness corresponds to the shape factor (S).
According to [168,169], designing elastomeric bearings with a high shape factor provides some benefits, such as mitigating rocking motion and enhancing stability. Most guidelines recommend designing elastomeric bearings with shape factors from 15 to 25 [170] or 30 [171]. However, Equation (12) demonstrates that the vertical stiffness would be thousands of times greater than horizontal stiffness for these high shape factors. Additionally, the actual value is potentially more, as rubber compressibility is neglected. This design procedure could be effective in horizontal isolation, but it causes vertical isolation periods ranging from 0.03 to 0.15 s and significant amplification of the vertical acceleration. Lead-rubber bearings are an enhancement of low-damping natural rubber bearings by adding a lead plug in the central hole of the bearing, refer to Figure 11. This plug deforms plastically under shear deformation, which increases the energy dissipation capabilities. The resulting horizontal force-deformation relationship of a lead-rubber bearing is defined using bilinear behavior, as shown in Figure 12.
where and are the shear yield strength and the cross-sectional area of the lead plug, respectively. For a given horizontal displacement, d, and the second-slope stiffness, , the effective or secant stiffness, , can be defined by Equation (14). The effective shear modulus, , is estimated for a given , then Equation (11) is applied to define the vertical stiffness.
= + OpenSees [10] numerically includes various commands developed specifically for elastomeric bearings. As shown in Figure 13, the elastomericBearingPlasticity and elasto-mericBearingBoucWen elements adopt the traditional bearings elements with customizing shear components, where is determined according to following equation: where , are, respectively, the post-yield stiffness ratio of linear and nonlinear hardening, is the initial stiffness, dy is the displacement, and µ is the exponent of nonlinear hardening. The elastomericBearingPlasticity element requires defining and the associated parameters with to simulate horizontal isolation. However, the elastomer-icBearingBoucWen requires additional parameters to predict the hysteretic energy dissipation precisely, and they are typically the yielding exponent (sharpness of hysteresis loop corners, η), first hysteretic shape parameter (β), and second hysteretic shape parameter (γ). Furthermore, the elastomericX, LeadRubberX, and HDR elements are developed for modeling circular low-damping rubber, lead-rubber, and high-damping rubber isolators, respectively. These commands might be more accessible as no predefined materials are required. Instead, modeling requires geometric properties, e.g., internal diameter, outer diameter, the number and thickness of steel or rubber layers, and material properties, such as shear modulus, bulk modulus, and yield strength. These commands were developed by [172] to provide an automation procedure for generating mechanical properties with a broad verification investigation. For general modeling purposes of square or circular shapes, the KikuchiBearing element, which was developed by [173], provides an elastomeric bearing element. It consists of multiple shear spring (MSS) and multiple normal spring (MNS) models, as shown in Figure 14. The required parameters are the assigned material, the total number of shear and normal springs, and the geometric parameters, e.g., the cross-section, total rubber thickness, and the total height.  Some isolation systems, usually low-damping rubber bearings, adopt supplemental damping devices across the isolator interface to control displacements and increase the provided damping ratio. The most common response control dampers can be categorized into steel hysteretic dampers, viscoelastic dampers, and viscous fluid dampers, as shown in Figure 15. A steel hysteretic damper is installed in a brace configuration and composed of low-yield-strength steel inserted between a pair of steel channels that prevent the plate from buckling. The device mechanism depends on absorbing energy through the elastic deformations of the plate. A similar technique is used for viscoelastic dampers but uses a high-damping rubber to absorb the vibration energy and convert it into heat [175]. The viscous fluid damper is a cylinder filled with hydraulic fluid and a piston that can move in both directions. The hydraulic fluid flows from one chamber to the other through valves. The motion of a piston resists this flow and causes energy dissipation. An axial hinge assigned to hysteretic uniaxial material could be a straightforward modeling approach for damper elements. The viscous material can be used for this model, which adopts the stress-strain relationship expressed in Equation (16).
where is the instant stress, refers to the strain rate, is the power factor (when set to 1, it makes a linear damper), and is called for the damping coefficient. Also, the Vis-cousDamper was developed to simulate the hysteretic response of nonlinear viscous dampers. It implements the Maxwell Model (linear spring and nonlinear dashpot in series), demonstrated in Equation (17).
where is the damper's extensional velocity, and represent the applied force and velocity, respectively, is the elastic stiffness of a linear spring, is the viscous damper parameter, and α is the viscous damper exponent.

Software Description
This research proposes extending the NHERI-SimCenter application framework to implement seismic isolation (SI) modeling rather than incremental dynamic analysis (IDA). The NHERI-SimCenter EE-UQ tool provides the capability to perform simulations and estimate the EDPs of multi-degree-of-freedom (MDOF) system structures, considering input uncertainties. Moreover, the PBE tool completes these steps with damage and loss measures. Some effort was required in the frontend and backend components to integrate previously mentioned isolator models with these tools. Figure 16 shows the added panel in the user interface to model multiple isolation systems. Since some SI models require predefined objects of materials or friction models, it is divided into three tabs. The user can define a list of materials and friction models (if needed) and assign their elements as properties for an isolator object. Combining dampers with the bearing elements is a valid option, rather than assigning their axial materials. This interface contains wrappers for all previously discussed bearing elements, all friction models in Opensees [141], damping materials, i.e., viscous, and ViscousDamper, and general behavior uniaxial materials, i.e., elastic, elastic-perfectly plastic, and hysteretic. Finally, an isolator object should be assigned to the MDOF system, as shown in Figure 17, while allowing the default option of the fixed base case. However, the backend extensions include adding several commands in wrapping scripts to retrieve the newly defined objects, execute their commands, and define new EDP cases for the isolator's horizontal displacements and damper forces.

Illustrative Examples
A simple case study is discussed herein to present the impact of integrating the recent features, i.e., seismic isolation (SI) modeling and incremental dynamic analysis (IDA), with legacy workflow models, e.g., uncertainty quantification (UQ) and PELICUN. The benchmark structural system is a realistic MDOF system introduced by [176].
Additionally, it was adopted in a previous publication by the lead author study [118] for multiple regression analysis (MRA) of natural rubber bearings (NRB) combined with viscous fluid dampers (VFDs) with isolators. The SI system was evaluated using a realistic five-story base-isolated building and an array of feasible damper combinations. Two ensembles of seven NF earthquake records were utilized to represent different seismic hazard levels. The key response parameters, including total maximum displacement, peak damper force, and top story acceleration ratio, were examined, and the MRA models provided satisfactory results with less computation. Table 1 shows the material properties (the mass and stiffness for each story) and the period of vibrations for the first six modes of the adopted system. Three models were generated based on this system: fixed base (A-FixedBase), isolated with only NRB (B-NRB-SI), and isolated with NRB combined with VFDs (C-NRB-VFD-SI). The damping coefficient (C) and the damping exponent (α) of the mean values were assumed to be 350 kN-s/m and 0.7, respectively. Model uncertainties were defined with a normal distribution, with a coefficient of variation values ranging from 5 to 20. The forward propagation Monte Carlo method was selected for the distribution of random values. In addition, the interactive interface for PEER Strong Ground Motions Databases was utilized to select 10 earthquake records, and their properties are defined in Table 2. It is important to highlight that the focus of this study was on the nearfield ground motion records, given that SI systems are typically employed in highly active seismic zones. The exclusion of far-field records is a deliberate decision that allows for a more specific and targeted investigation of the behavior of the combined NRBVFD system under the conditions most relevant to its practical application. This approach enables the researchers to achieve a more accurate and realistic assessment of the SI system's performance in response to the seismic events that are most likely to occur in the regions where these systems are implemented. The filtration criterion was a magnitude ranging from 6.5 to 7.5 and up to 12 km distance to the fault, as mentioned in [118]. A performance-based assessment was conducted for all three models, including structure responses, damage measures, and scaling factors for ground motions in different intensities.
The effect of the isolation systems can be measured by comparing structure-response characteristics. Table 3 introduces this analogy based on several EDPs: the peak floor acceleration (PFA) on the top floor, the peak inter-story drift (PID) on the first floor, and the peak roof drift (PRD). Comparing the "A-FixedBase" and "B-NRB-SI" models illustrates the impact of this isolation system as follows: (1) PFA range values were decreased from (0.3816-2.5476) to (0.2570-0.9959) m/s 2 , typically a 47% reduction for the mean value, as shown in Figure 19a; and (2) the fixed base model reported mean values of 0.0649 and 0.0394 for PID and PRD, while the isolated models provided only 0.0003 and 0.0001 drift ratios.   In addition, adding VFDs to the isolation system includes controlling the horizontal displacements at the isolation level to protect the bearing device. As a result, the total maximum displacement (DTM) for the bearing element was decreased by 4% on average (refer to Figure 19b). The peak damper force (PDF) values were estimated from 129.14 to 244.24 KN. For the global behavior, the mean values of PID and PRD ratios enlarged to 0.0018 and 0.0007 due to increasing Kh. In addition, a slight enhancement occurred regarding minimizing PFAs. This variation reflects the importance of the relative properties between dampers and bearing elements to compromise acceptable response.
An incremental dynamic analysis was conducted by scaling the ground motions up to 5 with 0.25 increments for a better assessment. Damages and loss analyses were applied according to the HAZUS approximate approach [149] Arbitrary values were assigned to replacement cost, replacement time, and population. Damage measures (DMs) are classified into three damage states (DS): structural (S), non-structural drift-sensitive (NSD), and non-structural acceleration-sensitive (NSA), while losses are described according to the time and cost for repairs, in addition to four levels of injuries. Table 4 shows DVs for the fixed base model at different intensities of motion, while Table 5 compares the results among the three models at the maximum scale factor. Damage and loss values followed the same trends of structural responses. Moving from a fixed base to an NRB isolation, the variations of the results can be summarized as: (1) a substantial decrease occurred for S and NSD damage states as a result of minimizing drifts ratios. For instance, the reported value for "B-NRB-SI" at maximum estimated damage state is about 0.08 of its counterpart; and (2) the estimation of the NSA damage state was reduced by about 50% on average, which is close to the reduction ratio in acceleration, as shown in Figure 20. Naturally, minimizing damages is accompanied by reduced losses, so the time and cost for repairs become almost negligible. However, the tiny range of damages in the isolated models impacts the employed VDPs to be extended.    The tendency of reducing accelerations and maximizing drifts is reflected in corresponding damage states. Although calculating decision variables reflected response behavior, this transformation is necessary. It is easier to judge the design regarding monetary losses and injuries than the drift ratios or floor accelerations.

Conclusions
The structural engineering community started to recognize the limited advantages of the traditional, deterministic earthquake design philosophy. PBEE was established as a probabilistic framework containing hazard, structure, designing uncertainties, and shifting decision variables from engineering parameters to general variables of losses, injuries, monetary, and repair downtime variables. The PBEE Vision 2000 report and FEMA-356 presented the first generation of PBEE. Then, the PEER Center introduced a more comprehensive and rigorous methodology, which comparted the design process into four modules of hazard, structural, damage, and loss analyses. This conceptual evolution created a demand for more flexible and generic computer programs. The NHERI-SimCenter team has developed an intelligent framework of well-structured integrated mini-programs. This architecture is interactive with the different interests of researchers corresponding to the studied hazard, asset classification, and desired outputs. The integration of uncertainty quantification procedures with finite elements simulation could be the most significant feature of the program. The framework design also permits futural development for more customized and scalable cases. This research adopted the development of seismic isolation components and an incremental dynamic analysis user interface as addons to the design. Seismic isolation is a powerful technique to control response and minimize damages implemented in several critical structures. It is based on decoupling the superstructure and foundations and absorbing ground motion energy at the isolation level. According to their energy-dissipation mechanisms, seismic isolation bearings are classified into friction-based and elastomeric-based. Several finite element models are provided for each bearing element type. All these models were implemented in the developed add-on software components. A simple case study was included to reflect the impact of integrating the developed components with previous components. An incremental dynamic analysis was conducted for three cases of a benchmark five-story structure, as a fixed base, NRB isolation, and NRB-VFD combined isolation system. For this case study, the uncertainties of structural properties and multiple earthquake motions were applied. The corresponding results demonstrated the efficiency of the isolation technique by reducing the drift ratios and drift-sensitive components' damages into negligible values. Moreover, it produced about a 50% reduction in floor accelerations and corresponding component damages. In addition, VFDs protect the bearing element by controlling the horizontal displacements at the interface, affecting the global isolation behavior. This demonstrates how critical it is to calibrate dampers and the properties of bearing elements to produce optimal designs. This process requires several analysis simulations according to uncertainties in the model parameters and earthquake intensity.