Reliability-based optimal design of nonlinear viscous dampers for the seismic protection of structural systems

Abstract

Viscous dampers are widely employed for enhancing the seismic performance of structural systems, and their design is often carried out using simplified approaches to account for the uncertainty in the seismic input. This paper introduces a novel and rigorous approach that allows to explicitly consider the variability of the intensity and characteristics of the seismic input in designing the optimal viscous constant and velocity exponent of the dampers based on performance-based criteria. The optimal solution permits controlling the probability of structural failure, while minimizing the damper cost, related to the sum of the damper forces. The solution to the optimization problem is efficiently sought via the constrained optimization by linear approximation (COBYLA) method, while Subset simulation together with auxiliary response method are employed for the performance assessment at each iteration of the optimization process. A 3-storey steel moment-resisting building frame is considered to illustrate the application of the proposed design methodology and to evaluate and compare the performances that can be achieved with different damper nonlinearity levels. Comparisons are also made with the results obtained by applying simplifying approaches, often employed in design practice, as those aiming to minimize the sum of the viscous damping constant and/or considering a single hazard level for the performance assessment.

1 Introduction

Energy dissipation devices are widely employed to improve the performance of structural systems under earthquake input. In particular, the use of viscous dampers represents a valuable and cost-effective solution to achieve a structural performance enhancement (Christopoulos et al. 2006; Soong and Dargush 1997; Takewaki 2011). Viscous dampers provide a resisting force proportional to the relative velocity between the damper ends to the power of an exponent α. This exponent can assume values in a wide range, usually between 0.15 and 2 (Constantinou and Symans 1992). In structural engineering applications, values lower than 1 are preferred to limit the increase of forces for velocities beyond the design ones.

Many methodologies were proposed in the last decades for the optimal design of viscous damper properties in buildings subjected to seismic ground motion records (Constantinou and Symans 1992; Martinez-Rodrigo and Romero 2003; Uriz and Whittaker 2001). The main limitation of these methodologies is that they do not account explicitly for the uncertainties in the earthquake input or in the model. Thus, they cannot be employed to control directly the seismic risk of the protected structure and of the dampers.

More recent works have demonstrated the important effect of the exponent α on the propagation of these uncertainties in the seismic response assessment of buildings equipped with viscous dampers. Di Paola et al. (2007) and Tubaldi and Kougioumtzoglou (2015) proposed a stochastic dynamic approach based on the equivalent linearization of the damper behaviour to study the problem. Tubaldi et al. (2014) and Dall’Asta et al. (2016) analyzed the performances of respectively single degree of freedom (SDOF) and multi degree of freedom (MDOF) systems with nonlinear viscous dampers by developing response hazard curves for values between 0.15 and 1, whereas Guneyisi and Altay (2008) carried out fragility studies on systems equipped with linear viscous dampers. The work of Lavan and Avishur (2013) analyzed the influence of the uncertainties concerning model parameters on the seismic performance. Dall’Asta et al. (2017) studied the specific effects of the damper property variability by linking the tolerance allowed in quality control tests and production tests to the relevant variability of the risk of exceedance of the most interesting parameters controlling the seismic performance. Despite the importance of a probabilistic approach for the performance evaluation, only few methodologies have been developed so far (Marano et al. 2007; Taflanidis 2010; Tubaldi et al. 2015) that can be employed to control the risk of failure in buildings equipped with dampers, or the costs associated to it (Taflanidis and Beck 2009; Gidaris and Taflanidis 2015). None of these studies have however investigated in depth the influence of the velocity exponent α on the design results.

This study proposes a methodology for the optimal design of viscous dampers in building frames accounting for the effects of the uncertainty relevant to the intensity, duration and the frequency content of the seismic input. The methodology allows finding the optimal damper properties that minimize an objective function related to the dampers cost, while considering a constraint on the performance of the frame. In particular, the sum of the damper forces, a quantity dependent on the response, is considered as objective function, differently from other existing methodologies (Marano et al. 2007; Tubaldi et al. 2015) which consider the sum of the damper viscous constants, a deterministic quantity, to simplify the design problem. The proposed approach permits to compare the design solutions obtained for different levels of the damper velocity exponent.

The stochastic objective function and stochastic constraint are evaluated via Subset simulation (Au and Wang 2014; Au and Patelli 2016) and the auxiliary response method (Au and Wang 2014). OpenCossan (Patelli et al. 2016) is employed to solve the design problem through a direct approach (Patelli and De Angelis 2012). COBYLA algorithm (Powell 2007) is used in the optimization process (outer loop), while the structural analyses for the multiple earthquake records (inner loop) are carried out by using OpenSees (McKenna et al. 2000).

The rest of the paper is organized as follows: in Sect. 2 an overview of the Reliability-based design optimization (RBDO) problem is given, in Sect. 3 the methodology for the solution of the RBDO problem is described together with simplified approaches commonly employed in the literature; in Sect. 4 the RBDO methodology is applied to the design of the dampers for a case study consisting of a steel moment-resisting frame, and the results are compared with those obtained by using the simplified approaches illustrated in Sect. 3.

2 Reliability based design optimization problem formulation

In general, the reliability-based design of the viscous dampers properties for the mitigation of the seismic risk of a structural system should account for the uncertainties inherent to the seismic input, as well as to the model. In the present formulation, the uncertain parameters are represented by the vector Θ with joint probability density function $f_{{\varvec{\Theta}}} \left( {\varvec{\uptheta}} \right)$, collecting all the quantities necessary to completely describe the probabilistic properties of both the seismic action (Pinto et al. 2004) and the structural model (e.g., floor mass, geometry). The vector x collects the design variables, i.e., the damper properties that need to be designed.

In order to characterize the seismic input, it is necessary to specify a probabilistic model for the occurrence, intensity, frequency content, and duration of the earthquakes that may strike the building site. “Appendix” reports a complete description of the earthquake model considered in this study, which assumes that the earthquake occurrences follow a Poisson model with rate λ.

The seismic risk is expressed in terms of the probability of exceeding different limit states related to the performance of the components of the structural system during the design life time T _L . By assuming that the structural model is homogeneous and does not change at each earthquake occurrence, the risk for a generic component can be computed from the probability that the engineering demand parameter (EDP), used to monitor the response, exceeds the value edp at any earthquake occurrence. This probability, denoted as $P_{EDP} \left( {\left. {edp} \right|{\mathbf{x}}} \right)$, depends on the value assumed by the design variables, and its estimation corresponds to the solution of a direct reliability problem. The corresponding risk of exceedance in the design lifetime T _L can be expressed as (Soong and Dargush 1997):

$$P_{{EDP,T_{L} }} \left( {\left. {edp} \right|{\mathbf{x}}} \right) = 1 - \exp \left( { - \lambda \cdot P_{EDP} \left( {\left. {edp} \right|{\mathbf{x}}} \right) \cdot T_{L} } \right)$$

(1)

The design of the dampers is an inverse reliability problem aiming at finding the values of the design variables x ^*such that the system meets a pre-defined performance level. This problem can be cast in the form of an RBDO problem (Patelli and De Angelis 2012; Schueller and Jensen 2008; Jensen and Sepulveda 2012; Beck et al. 2014), aiming at identifying the optimal damper properties which minimize the dampers cost while satisfying the stochastic constraints on the probability of exceeding a prescribed global damage level. Finally, additional constraints are needed to ensure that the values assumed by the design variables are physically admissible.

In the problem at hand, the design variables and the objective function depend on the damper constitutive law:

$$f_{d} \left( {\dot{u}} \right) = C_{d} \left| {\dot{u}} \right|^{\alpha } \,\cdot\, sign\left( {\dot{u}} \right)$$

(2)

where $\dot{u}$ represents the relative velocity between the damper’s ends, C _d the damping coefficient and α the parameter that describes the nonlinearity of the damper response.

In particular, the design variables are the coefficients C _d of the N _d dampers to be added to the building frame at the various floors. These are collected in the vector x = [C _d,1 C _d,2 …C _d,Nd ]^T. In order to evaluate explicitly the influence of the damper exponent α on the design solution and performance enhancement that can be achieved with viscous dampers, this parameter is not included in x and the optimization process is repeated for different values of α. However, the optimization scheme is general and this parameter could also be included in the design process, if required.

In this study, the objective function is expressed as a function of the sum of the maximum forces observed in the dampers $\varXi = \sum\nolimits_{i = 1}^{{N_{d} }} {f_{d,i} }$. In this regard, it is important to observe that many works (Lavan and Avishur 2013; Marano et al. 2007; Tubaldi et al. 2015) employ the sum of the viscous damping constants $\sum\nolimits_{i = 1}^{{N_{d} }} {C_{d,i} }$ as objective function for the optimization problem. However, the cost of the dampers depends on the forces and the strokes the dampers have to withstand rather than on the values of the damper viscous constants (Bahnasy and Lavan 2013; Pollini et al. 2017). Considering the viscous damping constant leads to a deterministic objective function, whereas considering the forces leads to a stochastic objective function, since these parameters, according to the constitutive law, depend on the response through the velocity between the damper ends. Thus, the objective function needs to be formulated in terms of the values of the sum of the damper forces, which correspond to a prefixed probability of exceedance $\bar{P}$. If $P_{\varXi } \left( {\xi ,{\mathbf{x}}} \right) = P\left[ {\left. {\varXi \ge \xi } \right|{\mathbf{x}}} \right]$ denotes the probability of Ξ exceeding ξ, then the objective function is $\phi \left( {\mathbf{x}} \right) = \bar{\xi }\left( {\mathbf{x}} \right)$, where $P_{\varXi } \left( {\bar{\xi },{\mathbf{x}}} \right) = \bar{P}$. In (Altieri et al. 2017), the expected value of ξ is used as objective function.

In this study, the system performance, related to the stochastic constraint, is assumed to be controlled entirely by the inter-story drifts, i.e., the structure is assumed to fail if the maximum inter-storey drift among the various storeys ∆ exceeds a given limit related to the building damage $\bar{\delta }$. The same limit $\bar{P}$ assumed for the sum of the damper forces is considered for the probability of exceedance of the drifts $P_{\Delta } \left( {\bar{\delta },{\mathbf{x}}} \right) = P\left[ {\left. {\Delta \ge \bar{\delta }} \right|{\mathbf{x}}} \right]$.

By limiting the inter-storey drifts, the strokes of the dampers are indirectly controlled. It is noteworthy that other limit states of interest for monitoring the performance of the structural and non-structural building components could be included in the proposed formulation (Pinto et al. 2004; Freddi et al. 2013) and controlled by introducing additional constraints. The RBDO problem can be mathematically formulated as follows:

$$\begin{aligned} \mathop {\text{argmin}}\limits_{{\mathbf{x}}} \, \phi \left( {\mathbf{x}} \right) \, \hfill \\ {\text{subject to:}} \hfill \\ {\mathbf{c_{\it{i}}}}\left( {\mathbf{x}} \right) \le {\mathbf{0}} \, \quad \left( {i = 1,2, \ldots ,m} \right) \hfill \\ P_{\Delta } \left( {\bar{\delta },{\mathbf{x}}} \right) - \bar{P} \le 0 \hfill \\ \end{aligned}$$

(3)

where c _i (x) ≤ 0 = (linear and/or nonlinear) deterministic constraints specifying the feasible domain of the damper properties. This constraint can be used to introduce a lower bound for the viscous damping constant, thus avoiding excessively small dampers.

It is envisaged that in the case of large-scale structures the proposed design approach may lead to an excessive number of dampers with different properties. Thus, some post-processing of the design results may be required, e.g. grouping dampers with similar viscous constants at different locations to have the same damper property. More sophisticated design formulations such as the one described in Lavan and Amir (2014) could also be implemented in Eq. (3) to allow for a limited number of size groups of dampers in the design.

3 Proposed solution

3.1 SUBSET-RBDO method

The RBDO problem solution through a direct approach (Patelli and De Angelis 2012) entails performing, for each optimization loop (outer loop), a full reliability analysis (inner loop). This subsection describes a method, denoted to as SUBSET-RBDO method, which uses Subset simulation (Au and Wang 2014; Au and Patelli 2016; Zuev 2017) for performing the reliability analysis. Subset simulation is an advanced stochastic approach for simulating rare events and estimating the corresponding small tail probabilities. The basic idea is to express the rare-event probability as a product of larger conditional probabilities by introducing a sequence of less rare events.

In the problem at hand, the maximum inter-storey drift is assumed as the driving variable inside the Subset simulation and the probability to be computed is $P_{\Delta } \left( {\bar{\delta },{\mathbf{x}}} \right)$, i.e., the probability that the IDR exceeds the threshold $\bar{\delta }$. Subset simulation allows a drastic reduction in the computational costs associated with the reliability problem by considering this probability as a product of larger conditional probabilities for higher IDR thresholds $\delta_{1} > \delta_{2} > \cdots > \delta_{m} = \bar{\delta }$, which can be computed by considering a reduced set of samples of Θ.

The samples of Θ generated in the application of Subset simulation can also be used to obtain the statistics of the damper forces through the auxiliary response method described by Au and Wang (2014). This way, one can avoid executing another reliability analysis assuming the sum of the damper forces as driving variable. The auxiliary response method is applied by subdividing the sample space, Θ, generated with Subset simulation, based on the intermediate threshold values δ _i , i = 1,2,…,m. These threshold values define a sequence of mutually exclusive and collectively exhaustive events B _i , i = 0,1,…,m − 1, such that B _i  = [δ _i  ≤ ∆ ≤ δ _i+1]. These events form a partition of the sample space. Thus, by exploiting the Total Probability theorem, the exceedance probability for the sum of the damper forces Ξ can be computed as:

$$P\left[ {\varXi \ge \xi } \right] = \sum\limits_{i = 0}^{m - 1} {P\left[ {\left. {\varXi \ge \xi } \right|B_{i} } \right] \cdot P\left( {B_{i} } \right)}$$

(4)

where P(B _i ) is the probability of event B _i , which can be computed by dividing the number M _i of samples of ∆ contained in B _i by the total number of samples. For a given i, let Θ _ik , k = 1,2,…,M _i denote the samples of Θ conditional on B _i . Hence, the conditional probability P(Ξ > ξ|B _i ) is evaluated by:

$$P\left[ {\left. {\varXi \ge \xi } \right|B_{i} } \right] = \frac{1}{{M_{i} }}\sum\limits_{k = 1}^{{M_{i} }} {I\left( {\varXi_{ik} \ge \xi } \right)}$$

(5)

where I() is the indicator function, and Ξ _ik  = Ξ(Θ _ik ) is the value of the sum of the damper forces, corresponding to the sample Θ _ik .

Once the statistical distribution of Ξ is known, the objective function can be estimated by solving the problem: $P_{\varXi } \left( {\bar{\xi },{\mathbf{x}}} \right) = \bar{P}$. Figure 1 illustrates the application of the SUBSET-RBDO method.

Fig. 1

Description of the SUBSET-RBDO method

Using a simulation method for computing the failure probabilities introduces some noise in the functions ($P_{\Delta } \left( {\bar{\delta },{\mathbf{x}}} \right)$ and ϕ(x). This results in numerical difficulties and convergence problems if gradient-based algorithms are employed for the optimization process (Tubaldi et al. 2015). For this reason, the Constrained Optimization by Linear Approximation algorithm (COBYLA) (Powell 2007) is employed in this work. COBYLA is a gradient free algorithm for the solution of constrained optimization problems. The algorithm operates by evaluating the objective function and the constraints at the vertices of a trust region. If the optimization problem has a total of N design variables, then the trust region has a total of N + 1 vertices. With this information, linear approximations of the objective function and constraints are employed during the optimization process. In general, for smooth functions, the rate of convergence of COBYLA is slower than that of gradient-based algorithms, i.e. more function evaluations are required to find the optimum. However, one of the salient features of COBYLA is its robustness, which makes it suitable for noisy functions, and the low number of parameters that need to be tuned for performing optimization.

3.2 Simplified approaches for RBDO solution

Although the use of Subset simulation and the auxiliary response method allows estimating the exceedance probabilities of ∆ and Ξ with a relatively reduced number of samples, the application of the SUBSET-RBDO method remains still quite cumbersome for practical design applications. This sub-section describes some approximations commonly introduced to simplify the damper design problem, by reducing the number of simulations involved in the design process. These approximations concern the performance evaluation and the definition of the objective function.

In particular, with reference to the first aspect, many design codes and methodologies evaluate the seismic performance by considering only a single hazard level (SHL), corresponding to a target probability of exceedance during the design life. This approach is denoted as intensity-based demand assessment in (Bradley 2011), and employs the concept of intensity measure (IM), which is a scalar measure of the earthquake intensity. The target exceedance probability during the lifetime T _L for the IM can be assumed equal to $\bar{P}$, which is the target value of the risk in the RBDO problem, and the corresponding IM value, denoted as im*, can be obtained as the solution of the equation:

$$P_{{IM,T_{L} }} \left( {im} \right) = 1 - \exp \left( { - \lambda \cdot P_{IM} \left( {im} \right) \cdot T_{L} } \right) = \bar{P}$$

(6)

where P _IM (im) is the probability of IM exceeding the value im given an earthquake occurrence.

The latter can be evaluated again by employing the Total Probability theorem based on the seismic input model parameters Θ:

$$P_{IM} \left( {im} \right) = \int\limits_{IM} {P_{IM} \left( {\left. {im} \right|{\varvec{\uptheta}}} \right){\mathbf{f}}_{{\varvec{\Theta}}} \left( {\varvec{\uptheta}} \right)d{\varvec{\uptheta}}}$$

(7)

At each iteration, Subset simulation, or any other simulation technique, can be employed to generate a set of seismic input records characterized by the target im ^* value. The seismic inputs are then employed to evaluate, via time-history analysis, the constraint and the objective function, the former based on the mean value of the IDR, the latter on the mean value of the sum of the damper forces. The same set of records can be employed for evaluating the demand at the different iterations. This corresponds to implementing an exterior sampling approach (Taflanidis and Beck 2008), which reduces the noise in the objective function and in the stochastic constraint, helping to limit the number of iterations required to achieve convergence. Figure 2 illustrates the application of SHL-RBDO method for the performance assessment.

Fig. 2

Illustration of SHL-RBDO method

With regard the objective function evaluation, a simplification often introduced in design methodologies is to consider the sum of the viscous damping constants $\sum\nolimits_{i = 1}^{{N_{d} }} {C_{d,i} }$, rather than the sum of the damper forces $\varXi = \sum\nolimits_{i = 1}^{{N_{d} }} {f_{d,i} }$, as objective function. This approach simplifies the objective function, by making it independent of the stochastic response, and proportional to the model parameters. However, it should be noted that the damper forces depend on the inter-story velocities and the damping coefficients via Eq. (2). Thus, if the response is dominated by a single mode of vibration and if the inter-storey drifts are similar at all the floors, the peak inter-story velocities will also be similar, and using $\sum\nolimits_{i = 1}^{{N_{d} }} {C_{d,i} }$ rather than Ξ as objective function, should lead to similar results. However, in the case of buildings with nonlinear behavior, with the response dominated by higher modes, and with irregularities and different IDR demands at the various floors, the results obtained with the two approaches could be very different.

4 Case study

The case study considered for illustrating the application of the design methodologies consists of a 2d model of a 3-storey steel moment-resisting building frame. This structure has already been employed in several projects investigating the efficiency of seismic control (Barroso and Winterstein 2002; Ohtori et al. 2004). It has been designed in compliance with local code requirements and design practices for office buildings, by considering the gravity, wind, and seismic load. Figure 3 illustrates some geometrical properties of the system. The steel has a mean yield strength of 248 MPa for the beams, and 345 MPa for the columns. The total mass at the two lower levels is equal to 975 tons, whereas the total mass at the top floor is equal to 1040 tons. Further information about the structure can be found in Ohtori et al. (2004). A 2d model of the structural system has been developed in OpenSees (McKenna et al. 2000) by using distributed plasticity beam elements with nonlinear mechanical and geometrical behavior taken into account. A more detailed description of the finite element model is given in Dall’Asta et al. (2016). The first three modes of vibrations are characterized by a period of 0.99 s, 0.33 s, and 0.18 s.

Fig. 3

Geometrical properties of the case study

The seismic input for the building site is defined by the following values of the seismic hazard parameters: m _min = 5.5, m _max = 8, a = 4.35, and b = 0.9, corresponding to a mean annual frequency of seismic events $\lambda = \lambda_{M} \left( {m_{\hbox{min} } } \right) - \lambda_{M} \left( {m_{\hbox{max} } } \right) = 10^{{a - bm_{\hbox{min} } }} - 10^{{a - bm_{\hbox{max} } }} = \, 0.25\,{\text {years}}^{ - 1}$. The circular area from which earthquakes originate has a radius r _max  = 50 km. Generic soil conditions (Boore and Joyner 1997) are assumed for the soil at the site, corresponding to an average shear-wave velocity over the upper 30 m of soil layer of 310 m/s. The ground motion time histories are generated by using the Atkinson and Silva (2000) model, which is summarised in “Appendix”.

A diagonal damper brace is considered at each storey. The RBDO problem is solved by keeping the damper exponent α out of the optimization process, i.e., the design is repeated for different values of α. In particular, the optimal solution is sought for 8 discrete values of α, i.e., α = [0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]. Two different design scenarios are considered, one corresponding to uniform damper properties at the various storeys, the other corresponding to variable damper properties. The first case corresponds to having a single design variable, whose optimal value lies at the boundary of the reliability constraint and could also be found iteratively without resorting to sophisticated optimization tools. The second case, on the other hand, requires more advanced optimization tools such as COBYLA (Powell 2007), which allows to find the optimal damper placement.

The target probability $\bar{P}$ is equal to 0.10, for a design life time of 50 years. This corresponds to a mean annual frequency of exceedance of 0.0021 and a probability of exceedance for an earthquake occurrence of 0.0084, evaluated according to Eq. (1). In the following sub-sections the results of the application of the SUBSET-RBDO method are reported and compared with the results obtained by considering the previously discussed simplifications.

The numerical simulations in the inner loop have been carried out by exploiting the high-performance computing resources available at Liverpool University. In particular, a cluster of 31 CPUs has been employed to run the structural analyses in parallel with a final tolerance of 10⁻³ on the optimal solution.

4.1 SUBSET-RBDO results

This subsection illustrates the results of the application of the SUBSET-RBDO method for the case of uniform damper properties at all the storeys. In order to illustrate the application of Subset simulation method for the solution of the inner loop problem, Fig. 4a shows the samples of M and R, generated via Subset simulation, corresponding to different intermediate threshold levels reported in Fig. 4b, in the case of α = 0.3 and at the first iteration of the design process corresponding to a value of C _d  = 2000 kN^0.3 s^0.3 /m^0.3.

Fig. 4

Output of Subset simulation at the first iteration of the optimization process for α = 0.3 and C _d  = 2000 kN^0.3s^0.3/m^0.3 a samples of M and R corresponding to each Subset levels and b to the different IDR intermediate thresholds

It is observed that by increasing the IDR threshold level, the magnitude M of the simulated ground motion records increases, whereas the source-to-site distance R decreases. This trend is expected and is a consequence of the seismic attenuation law, described by the employed radiation spectrum, according to which the earthquake intensity increases for increasing M and decreasing R.

Figure 5a shows the samples of the target and auxiliary responses generated via Subset simulation at the optimal solution in the case of α = 0.3 and α = 1. The maximum interstorey drifts and the sum of the damper forces exhibit a significant direct correlation, i.e., by increasing the drifts also the forces increase. However, the trend of increase of the forces is highly nonlinear for α = 0.3, and is mainly a consequence of the damper constitutive law, reported in Eq. (2). Figure 5b shows the empirical complementary cumulative distribution function (CCDF) of the sum of the damper forces. This output is used to define the objective function $\phi \left( {\mathbf{x}} \right) = \bar{\xi }\left( {\mathbf{x}} \right)$, where $P_{\varXi } \left( {\bar{\xi },{\mathbf{x}}} \right) = \bar{P}$.

Fig. 5

Output from the auxiliary response method for α = 0.3 and α = 1 in terms of a relation between damper forces and IDR, b probability of exceedance of the damper forces

Figures 6a and 6b show the evolution of the objective function and the design variable respectively during the optimization process, carried out by using the SUBSET-RBDO method for α = 0.3 and α = 1, starting from an initial design variable of 2000 kN^0.3 s^0.3 /m^0.3 for α = 0.3 and 10,000 kNs/m for α = 1.

Fig. 6

Evolution of the objective function and design variable for a α = 0.3 and b α = 1

A maximum value of C _d  = 7000 kN^0.3s^0.3 /m^0.3 and C _d  = 15,000 kNs/m is considered for α = 0.3 and α = 1 respectively. After very few iterations both these quantities converge to a stable value equal to 1905 kN^0.3 s^0.3 /m^0.3(α = 0.3) and 10,030 kNs/m (α = 1) respectively. Similarly, there is convergence to 3746 kN (α = 0.3) and 5412 kN (α = 1) for ϕ(x ^*).

Figure 7 shows the values of the optimal solution x ^* = C ^* _d and the corresponding objective function ϕ(x ^*) obtained for the different α values. For each α value, Table 1 reports the starting C _d value and the number of iterations required before reaching convergence. On average, the optimal solution is found after 11 iterations. Each iteration corresponds to a full Subset simulation, involving about 800 model evaluations. Hence, 8800 structural analyses are required to obtain the optimum solution.

Fig. 7

Evolution of the optimal solution in terms of a design variable C _d and b objective function ϕ(x) for different α

Table 1 Starting design variable value and number of iterations required for each different α

The optimal solution increases with a nonlinear trend for increasing values of α, while ϕ(x ^*) increases with a linear trend for increasing α values. It can be observed that the use of non-linear dampers allows reducing the forces in the dampers, and thus corresponds to a more economical solution.

In fact, the value of ϕ(x ^*) for the case of α = 0.3 is about 70% of the value for α = 1. This can be seen also by looking at Fig. 5a, plotting the relation between the sum of the damper forces and the maximum IDR samples, in correspondence of the optimal solutions obtained for α = 0.3 and of α = 1. This relation reflects the damper constitutive law, and it is very close to being linear in the case of α = 1, and highly nonlinear with a maximum cap on the force level in the case of α = 0.3. Another effect of the damper constitutive law is that the CCDF of the forces in correspondence of the optimal solution is steeper for the nonlinear than for the linear case. Moreover, since the damper forces are further out of phase with the displacements and drifts in linear viscous dampers, the maximum absolute floor accelerations are generally higher for α = 0.3 than for α = 1. This can be observed in Fig. 8, plotting the samples of the maximum floor acceleration A _abs against the maximum IDR samples, for the two cases. The differences in the CCDF are notable for low values of the accelerations, and become negligible for larger values. The same trend has been observed in Dall’Asta et al. (2016) on a similar system, by considering a different earthquake model. In any case, the values of the maximum accelerations are found to be very low and the probability of damage of acceleration-sensitive non-structural building components is minimal (Taflanidis and Beck 2009).

Fig. 8

Output from the auxiliary response method for α = 0.3 and α = 1 in terms of a relation between maximum absolute acceleration and IDR, b probability of exceedance of the maximum absolute acceleration

Figure 9 shows how the optimal solution is influenced by the choice of the IDR limit $\bar{\delta }$. In particular, increasing the threshold, i.e., the allowable damage in the frame, results in a reduction of the optimal viscous constant and the objective function. The reduction is more significant for the nonlinear than linear viscous damper. For $\bar{\delta } = 1.48\%$, i.e., the maximum IDR for the bare frame, the forces are zero as no dampers are required.

Fig. 9

Evolution of the optimal solution in terms of a C _d and b ϕ(x) for increasing values of the maximum allowed interstorey drift δ, by fixing $\bar{P} = 0.0084$

In order to evaluate the robustness of the proposed methodology, the design procedure is applied again for the case corresponding to α = 0.3 by considering different starting points, i.e., values of C _d other than 2000 kN^0.3 s^0.3 /m^0.3. Figure 10 reports the evolution of C _d and ϕ(x) during the optimization process at different starting points, showing that the optimal solution is not affected by the initial value of C _d .

Fig. 10

Evolution of a design variable C _d and b objective function ϕ(x) for different starting points, case α = 0.3

Figure 11 shows the evolution of C _d and ϕ(x) by considering different numbers of samples for the intermediate steps of Subset simulation in the inner loop. The starting value is C _d  = 2000 kN^0.3 s^0.3 /m^0.3 and the maximum value allowed is C _d  = 3000 kN^0.3s^0.3/m^0.3. Sample sizes larger than 250 no longer have an effect on the results. Moreover, the design procedure has been repeated by considering different initial seeds for the random number generator used for sampling, no change in the optimal solution was observed, thus confirming again the robustness of the proposed procedure.

Fig. 11

Evolution of a design variable C _d and b objective function ϕ(x) for different numbers of samples at each intermediate level of the Subset simulation, case α = 0.3

Applying the SUBSET-RBDO method, by considering a design variable vector x = [C _d,1 C _d,2 C _d,3]^T instead of a scalar value C _d , allows evaluating the optimal solution for a variable distribution of the damping coefficients. This obviously leads to an increase in the computational costs of the optimization process. Figure 12 shows the variation, during the optimization process, of the viscous damping constant at the various storeys, and of the objective function, for the case corresponding to α = 0.3. It can be seen that the number of iterations increases drastically with respect to the uniform damper distribution case. In fact, the optimal solution is found after 29 iterations, for a total of about 24,900 structural analyses. This, together with the results of further analyses carried out by considering 6 design variables (corresponding to two different dampers at each storey), suggests that the number of iterations and thus of model evaluations increases close to linearly with the number of design variables. However, the choice of the initial point may also affect the number of iterations and further analyses should be carried out before drawing general conclusions.

Fig. 12

a Variation, during the optimization process, of the viscous damping constant at the various storeys and b of the objective function, in the case of variable damper distribution

Table 2 reports the results of the optimization process for the values of α = 0.3 and α = 1, which corresponds to the two extreme cases for the investigated range of variation of α. It can be observed that the value of C _d is highest at the second storey. This is also the storey undergoing higher drift demand when the undamped structures vibrate in its first mode of vibration. Moreover, it is interesting to observe that assuming a variable damper distribution yields lower values of the damper forces with respect to the case of uniformly distributed dampers. In particular, the relative reduction of the objective function is of 9% for α = 0.3, and of 17% for α = 1. However, it should also be pointed out that the reduction of the objective function may not represent an actual reduction in the total cost of the dampers, as employing dampers with the same characteristics at the different storeys may lead to some cost reductions (Lavan and Amir 2014).

Table 2 Optimal damper viscous constant and corresponding objective function for the case of dampers with variable distribution at the various storeys (α = 0.3 and α = 1)

4.2 Results for simplified RBDO methods

This subsection illustrates the results of the application of the simplified design approaches. First, the results of the application of the SHL-RBDO method are reported.

A total of 35,000 samples, generated via Latin Hypercube Sampling (LHS) by considering the seismic input model described previously, is employed to derive the hazard curve P _IM,TL (im), by assuming as IM the spectral acceleration at the fundamental building period and for 2% damping ratio, S _a (T,2%), and the peak ground acceleration, PGA. Figure 13 illustrates the corresponding hazard curves. In the same figure, the hazard curves obtained by considering Subset simulation with a reduced set of 3500 samples is considered. The two curves are almost identical.

Fig. 13

a Hazard curve in terms of probability of exceedance of S _a (T, 2%) in 50 years, b hazard curve in terms of probability of exceedance of the PGA in 50 years

The hazard curve has been then employed to identify the value of im ^*, which is characterized by the target risk of exceedance $\bar{P} = 10\%$. Successively, 20 records with the IM values closest to im* have been extrapolated from the sample space, and they have been employed to run the seismic analyses required in the inner loop.

Figure 14 reports the results of the optimization in terms of optimal values of the viscous constant $C_{d}^{*}$ and corresponding sum of the damper forces ϕ(x ^*). The results obtained by considering SUBSET-RBDO method are also reported for the sake of comparison. In general, it can be observed that the SHL-RBDO method leads to values of $C_{d}^{*}$ and of the corresponding objective function ϕ(x ^*) lower than the reference values obtained via the SHL-RBDO method. This result is expected, since a deterministic demand measure corresponding to a single hazard level is considered in the SHL-RBDO method, by disregarding the contribution from other hazard levels. This usually leads to underestimate the risk of exceedance (Dall’Asta et al. 2016; Bradley 2013). The maximum relative difference between the two methods is equal to 38% and 45% for C ^* _d and ϕ(x ^*) respectively (case α = 1).

Fig. 14

Results comparison between the RBDO and the simplified approach in terms of a design variable C _d and b objective function ϕ(x)

Figure 14 shows the results of the application of the SHL-RBDO method obtained by considering the PGA, rather than S _a (T,2%), as IM. It is observed that the optimal solution is affected quite significantly by the IM choice. This is a limit of an intensity-based demand assessment, as pointed out in Bradley (2013).

Finally, the effect of the objective function choice on the design results is investigated. In particular, the optimization is carried out once again with the SUBSET-RBDO method by considering the sum of the viscous damping constant as objective function. Figure 15 reports the variation with α of the optimal damping constant obtained by considering this simplified design approach and a uniform distribution of the damper properties. The values of the optimal damping viscous constants are very similar to those obtained by considering the sum of the damper forces as objective function. The two design solutions also return similar values of the total damper forces corresponding to a probability of exceedance equal to 0.0084 (Fig. 15b), with some scatter observed for high α values only.

Fig. 15

a Optimal C _d obtained via SUBSET-RBDO by considering the two alternative objective functions for the different α values, b damper forces sum for the different design approaches corresponding to a probability of exceedance equal to 0.0084

Table 3 reports the design results for the case of damping constants variable along the various storeys. Similarly to the case described in Table 2, the highest values of C _d are obtained at the second level. Moreover, the values of the objective function are lower than the corresponding ones obtained for the case of uniform damper distribution. The relative reduction is 20% in the case α = 1, and 8% in the case α = 0.3.

Table 3 Optimal damper viscous constant and corresponding objective function for the case of dampers with variable distribution at the various storeys (objective function equal to the sum of the damping constants)

5 Conclusions

A novel and rigorous methodology has been proposed for the optimal reliability based design of viscous dampers in building frames subjected to a stochastic earthquake input with uncertain intensity, duration, and frequency characteristics. The sum of the damper forces for the target exceedance probability is considered for the objective function definition. This parameter, accounting for the stochastic structural response, is more explicitly related to the dampers costs, differently from the other objective functions considered in the literature. An efficient reliability computational approach is used in the inner loop of the optimization process, and a robust gradient-free algorithm is employed for the optimization loop.

The proposed methodology has been applied to design the viscous dampers for protecting a steel moment-resisting frame, by considering both the cases of uniform and variable damper distribution along the building height. The following conclusions can be drawn from the design results:

• The combined use of Subset simulation and the auxiliary response method supports efficient evaluation of the quantities required by the inner loop. Thus, it can be conveniently employed for solving double-loop type RBDO problems.

• By reducing the damper velocity exponent α, the cost of the dampers, related to the sum of the damper forces for a target exceedance probability, decreases. In particular, a reduction of the sum of the damper forces of 30% is obtained for α = 0.3 compared to the case of α = 1, for a target failure probability of 10% in 50 years.

• Considering a variable distribution of the damper properties along the building height significantly reduces the cost of the dampers, while ensuring the same target performance of the system in terms of probability of exceeding the target value of the IDR.

• By changing the performance objective, the relative cost of the linear dampers with respect to the nonlinear dampers changes.

The methodology also provides reference results that can be used to evaluate the accuracy of some simplifications often introduced in design practice. In particular, the use of an intensity-based demand assessment has been investigated. Although this approach allows a significant reduction in the number of the model evaluations, it provides values of the optimal damper properties very different from those obtained by applying the proposed methodology, and this may lead to a violation of the pre-defined performance objectives.

A comparison with the results provided by the objective function commonly employed in other studies and consisting of the sum of the viscous damping constants has been carried out. The results of the comparison show that for the case study investigated, considering this simplified objective function yields optimal solution values very close to the reference ones obtained by considering the sum of the damper forces for the target exceedance probability. The proposed design methodology could be used in future studies to investigate the accuracy of other design approaches based on the use of metamodels for approximating the relation between the uncertain earthquake input parameters and the structural response.

Appendix

A point-source stochastic ground motion model is employed to describe the seismic input. This model is defined by the moment magnitude M and source-to-site (hypocentral) distance R of the seismic source, together with the ground motion radiation spectrum A(f) and the time envelope function e(t). The uncertainty of M is modelled by the truncated Gutenberg-Richter law, corresponding to the following probability density function (pdf) of M given an earthquake event:

$$p_{M} \left( m \right) = \frac{{\beta e^{ - \beta m} }}{{\beta e^{{ - \beta m_{\hbox{min} } }} - \beta e^{{ - \beta m_{MAX} }} }}m \in \left[ {m_{min} , m_{MAX} } \right]$$

(8)

where β = ln(10)b is a parameter related to the number of the expected earthquakes per annum with magnitude exceeding m. More precisely, it is assumed that the occurrence of an event with M > m is a Poisson process with exceedance frequency λ(m) = 10^a−bm and no event is expected for M > m _MAX . It is also assumed that no significant response is observed for M < m _min . The corresponding mean annual frequency of a seismic event is $\lambda = \lambda_{M} \left( {m_{min} } \right) - \lambda_{M} \left( {m_{MAX} } \right) = 10^{{a - bm_{min} }} - 10^{{a - bm_{max} }}$. For the uncertainty in event location, earthquakes of magnitude between M _min and M _max are assumed to occur equally likely in a circular area of radius r _max centered at the site where the structure is situated. This leads to a triangular pdf for the epicentral distance R confined to the interval [0,r _max ]:

$$f_{R} \left( r \right) = \left\{ \begin{aligned} 2r/r_{\hbox{max} }^{2} \, \quad r \in \left[ {0,r_{\hbox{max} } } \right] \hfill \\ 0 \quad {\text{ otherwise}} \hfill \\ \end{aligned} \right.$$

(9)

The ground motion for a given event of magnitude m and distance r is modelled as a time-modulated stochastic process, which is generated by starting from a white noise w(t), described by the N-dimensional vector w of values w _i assumed at the instant t _i  = i∆t, where ∆t is the finite time interval assumed for the numerical integration. Following the Atkinson–Silva model (Atkinson and Silva 2000), the ground motion is obtained by modulating in time the white noise by means of the function e(t), which yields the time-function z(t) = e(t)w(t). The amplitude and the frequency content are obtained by multiplying its Fourier transform z¯(f) (normalized to have a mean square amplitude of unity) by the radiation spectra A(f), where A(f) is a deterministic function of the frequency f. The final ground motion acceleration a(t) is obtained by the inverse Fourier transform of z¯(f)A(f). The time modulating function and the radiation spectra A depend on the moment magnitude, the distance and the local characteristic of soil. Further details on the ground motion model can be found in Taflanidis and Beck (2009).

The set of random variables Θ = [M, R, W] consists of the two scalar quantities (i.e., M and R), and the vector-valued quantity W, whose dimension N depends on the discretisation of the time interval.