Thermal power plant is an important life line structure^[1]. Maintaining its workability during and after earthquakes can be essential for people’s life and recovery^[1]. However, due to the functional requirements, thermal power plants are usually designed as complex structures with various irregularities^[2–6]. Besides, typically thermal power plants consists heavy coal scuttles at relatively high floors, which may generate significant inertial force and be detrimental to the structural seismic performance^[4]. An effective strategy to solve this problem is to convert the scuttles to sub-oscillators and tuning to the main structure, i.e.nonconventional multiple tuned mass damper (NC-MTMD)^[5–7].

During the real operational process, the coal storage in the scuttles may change^[6–7], which leads to random mass of the tuned mass dampers of NC-MTMD system. This can be classified as a mass uncertain NC-MTMD (MU-NC-MTMD) system. Jensen et al^[8]. first investigated TMD system with uncertain mass and proofed the necessities to consider uncertainty during design for moderate (coefficient of variation, $CoV $ =0.15) and high ( $CoV $ =0.30) cases.Mass uncertainty was considered in several subsequent studies^[9–10]. The mass uncertain tuned mass damper, as a case of isolated roof garden, was investigated in a previous study^[10]. The MU-NC-MTMD has multiple sub-oscillators, which is tuned from the original structure instead of tuning additional one, with mass varying from small mass ratio to large mass ratio^[6]. Despite NC-MTMD system contains multiple large mass ratio sub-oscillators and therefore has enhanced effectiveness and robustness^[6], those may not be a sure thing for MU-NC-MTMD systems. Therefore, performance and optimization of MU-NC-MTMD systems require further study. Besides, several random factors, such as stiffness and damping of sub-oscillators or parameters of main structures, may also lead to a different optimum design^[11] and therefore should be considered.

In this paper, a reliability-based optimization design framework for MU-NC-MTMD systems was proposed. Several parameters of hazard, main structures and scuttles were considered to be random. Unconditional failure probability considering multiple limit state bounds was adopted as objective function. A thermal power plant with multiple scuttles was used as a case to illustrate the application of the framework. Latin hypercube sampling (LHS) method was implemented to reduce the sample size needed and sampling size study was performed to determine the appropriate sampling size. Optimum design was obtained and discussed. A parametric study was further performed to study the influence of isolation mechanism and seismic gap. Pendulum MU-NC-MTMD and failure probability considering collision were investigated.

1 Structural model

For a NC-MTMD with a n_M DOF main structure (mass matrix ${{M}} $ , damping matrix ${{C}} $ , stiffness matrix ${{K}} $ ) and $n_ m $ TMDs (mass matrix $ {{m}}={\text{diag}}(m_j|j=1, 2, …, n_{\text{m}}$ ), damping matrix ${{c}}={\text{diag}}({c_j|j=1, 2, …, n_{\text{m}}}) $ , stiffness matrix ${{k}}={\text{diag}}({k_j|j=1, 2, …, n_{\text{m}}})) $ , the dynamic equation can be expressed as Eq. (1).

$\begin{aligned}[b]\left[ {\begin{array}{*{20}{c}}{{{{M}}_{{n_M} \times {n_M}}}} & {{O}}\\{{O}}& {{{{m}}_{{n_m} \times {n_m}}}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{{\ddot X}}}\\{{{\ddot x}}}\end{array}} \right] + \left[ {\begin{array}{*{20}{c}}{{{{C}}_{{n_M} \times {n_M}}} + {{{l}}^{\, \rm{T}}}{{{c}}_{{n_m} \times {n_m}}}{{l}}}& {{{{l}}^{\, \rm{T}}}{{{c}}_{{n_m} \times {n_m}}}}\\{{{{c}}_{{n_m} \times {n_m}}}{{l}}}& {{{{c}}_{{n_m} \times {n_m}}}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{{\dot X}}}\\{{{\dot x}}}\end{array}} \right] + \\\left[ {\begin{array}{*{20}{c}}{{{{K}}_{{n_M} \times {n_M}}} + {{{l}}^{\, \rm{T}}}{{{k}}_{{n_M} \times {n_M}}}{{l}}}& {{{{l}}^{\, \rm{T}}}{{{k}}_{{n_m} \times {n_m}}}}\\{{{{k}}_{{n_m} \times {n_m}}}{{l}}}& {{{{k}}_{{n_m} \times {n_m}}}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{X}}\\{{x}}\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{{{M}}_{{n_M} \times {n_M}}}}& {{O}}\\{{O}}& {{{{m}}_{{n_m} \times {n_m}}}}\end{array}} \right]{ l}{{a}}\end{aligned}$

(1)

where ${ X} $ and $ x$ are displacement vector of the main structure and TMDs; ${ l} $ is the location matrix with all 0 elements except the unit elements at $j $ th row and $p $ th column, where $ j$ takes 1, 2, $\cdots$ , $ n_{m}$ and $p $ is the label representing the order of the main structure DOFs that interacts with the sub-oscillators. $ {{\iota}}=[{{ I}}_{3\times 3},\;{{ I}}_{3\times 3},\; \cdots ,\;{{ I}}_{3\times 3}]^{\text{T}}$ and $ {{ I}}_{3\times 3}$ is a $ 3\times 3$ unit matrix; $ {{ a}}=[a_{\text{X}}, a_{\text{Y}}, a_{\text{Z}}]^{\text{T}}$ is accelerations; T is transpose.

By transforming only the main structure to its modal space, i.e., using the transformation in Eq. (2), and performing mode truncation^[12], one has the dynamic equation for the simplified NC-MTMD system (Fig. 1) as Eq. (3).

$\left[ {\begin{aligned}{{X}}\\{{x}}\end{aligned}} \right] = {{\hat \varPhi }}\left[ {\begin{aligned}{{q}}\\{{x}}\end{aligned}} \right] = \left[ {\begin{aligned}& {{{{\varPhi }}_{{n_M} \times {n_{\overline M}}}}} & {{O}}\\& {{O}} & {{{{I}}_{{n_m} \times {n_m}}}}\end{aligned}} \right]\left[ {\begin{aligned}{{q}}\\{{x}}\end{aligned}} \right] $

(2)

where ${{{\varPhi }}_{{n_M} \times {n_{\overline M}}}} $ is the mode matrix that composed of $ {n_{\overline M}}$ modes of the main structure; ${{{I}}_{{n_m} \times {n_m}}} $ is unit matrix; ${{q}} $ is generalized coordinate of the main structure.

$\begin{aligned}[b]\left[ {\begin{array}{*{20}{c}}{{{{{\overline M}}}_{{n_{\bar M}} \times {n_{\bar M}}}}}& {{O}}\\{{O}}& {{m}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{{\ddot q}}}\\{{{\ddot x}}}\end{array}} \right] + \left[ {\begin{array}{*{20}{c}}{{{{{\overline { C}}}}_{{n_{\bar M}} \times {n_{\bar M}}}} + {{{\varPhi }}^{\text{T}}}{{{l}}^{\text{T}}}{{cl\varPhi }}}& { - {{{\varPhi }}^{\text{T}}}{{{l}}^{\text{T}}}{{c}}}\\{ - {{cl\varPhi }}}& {{c}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{{\dot q}}}\\{{{\dot x}}}\end{array}} \right] + \\ \left[ {\begin{array}{*{20}{c}}{{{{{\overline K}}}_{{n_{\bar M}} \times {n_{\bar M}}}} + {{{\varPhi }}^{\text{T}}}{{{l}}^{\text{T}}}{{kl\varPhi }}}& { - {{{\varPhi }}^{\text{T}}}{{{l}}^{\text{T}}}{{k}}}\\{ - {{kl\varPhi }}}& {{k}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{q}}\\{{x}}\end{array}} \right] = {{{{\hat \Phi }}}^{\text{T}}}\left[ {\begin{array}{*{20}{c}}{{M}}& {{O}}\\{{O}}& {{m}}\end{array}} \right]{{\iota a}}\end{aligned}$

(3)

where $ {{\overline M}} = {{{\varPhi }}^{\text{T}}}{{M\varPhi }} = {\text{{\text{diag}}}}\left\{ {{{\overline M}_1}, {{\overline M}_2}, \cdots , {{\overline M}_{{n_{\overline M}}}}} \right\}$ , $ {{\overline { C}}} = {{{\varPhi }}^{\text{T}}}{{C\varPhi }}$ and ${{\overline K}} = {{{\varPhi }}^{\text{T}}}{{K\varPhi }} = {\text{{\text{diag}}}}\left\{ {{{\overline K}_1}, {{\overline K}_2}, \cdots , {{\overline K}_{{n_{\bar M}}}}} \right\} $ . ${{\overline C}} = {\text{{\text{diag}}}}\left\{ {{\overline C}_1}, \right. $ $\left.{{\overline C}_2}, \cdots , {{\overline C}_{{n_{\overline M}}}} \right\}$ if the main structure has proportional damping.

The uncertainty ubiquitously exists in the structural system and therefore the structural random variable set can be expressed as Eq. (4).

${{{\psi }}_S} = \left\{ {{{M}}, {{C}}, {{K}}, {{m}}, {{c}}, {{k}}} \right\}$

(4)

2 Earthquake hazard model

Stationary Kanai-Tajimi (KT) model^[13] was adopted to model the ground motion excitation. A three-dimensional filter, with one DOF in each direction, was assembled to the NC-MTMD system. For simplification, the filter was assumed to be isotropic and expressed as Eq. (5),

$\left\{\begin{aligned}& {{{I}}_{3 \times 3}}{{{\ddot x}}_{\text{f}}} + 2{\xi _{\text{f}}}{\omega _{\text{f}}}{{{I}}_{3 \times 3}}{{{\dot x}}_{\text{f}}} + \omega _{\text{f}}^2{{{I}}_{3 \times 3}}{{{x}}_{\text{f}}} = - {{w}}\\& {{a}} = {{{\ddot x}}_{\text{f}}}{\text{ + }}{{w}} = - \left( {2{\xi _{\text{f}}}{\omega _{\text{f}}}{{{{\dot x}}}_{\text{f}}} + \omega _{\text{f}}^2{{{x}}_{\text{f}}}} \right)\end{aligned}\right.$

(5)

where $x_{\text{f}} $ is the displacement vector of the KT filter; ${\xi _{\text{f}}} $ and ω_f are the damping ratio and circular frequency of the KT filter, respectively; $ {{w}} = {\left[ {\begin{array}{*{20}{c}}{{w_X}}& {{w_Y}}& {{w_Z}}\end{array}} \right]^{\rm{T}}}$ is the bed rock white noise excitation.

Substitute Eq. (5) into Eq. (3), one has Eq. (6),

${{{M}}_{{\text{ft}}}}\left[ {\begin{array}{*{20}{c}}{{{\ddot q}}}\\{{{\ddot x}}}\\{{{{{\ddot x}}}_{\text{f}}}}\end{array}} \right] + {{{C}}_{{\text{ft}}}}\left[ {\begin{array}{*{20}{c}}{{{\dot q}}}\\{{{\dot x}}}\\{{{{{\dot x}}}_{\text{f}}}}\end{array}} \right] + {{{K}}_{{\text{ft}}}}\left[ {\begin{array}{*{20}{c}}{{q}}\\{{x}}\\{{{{x}}_{\text{f}}}}\end{array}} \right] = {{{f}}_{{\text{bw}}}}$

(6)

where $ {{M}}_{\text{ft}}$ , $ {{C}}_{\text{ft}}$ , ${{K}}_{\text{ft}} $ and ${ f}_{\rm bw} $ are the mass, damping, stiffness matrix and load vector for the NC-MTMD system with the incorporation of the KT filter, respectively, and they are given by:

${{{M}}_{{\text{ft}}}} = \left[ {\begin{array}{*{20}{c}}{{{\overline M}}}& {{O}}& {{O}}\\{{O}}& {{m}}& {{O}}\\{{O}}& {{O}}& {{{{I}}_{3 \times 3}}}\end{array}} \right]$

(7)

${{{C}}_{{\text{ft}}}} = \left[ {\begin{array}{*{20}{c}}{{{\overline{ C}}} + {{{\varPhi }}^{\text{T}}}{{{l}}^{\text{T}}}{{cl{ \varPhi} }}}& { - {{{\Phi }}^{\text{T}}}{{{l}}^{\text{T}}}{{c}}}& { - 2{{{\varPhi }}^{\text{T}}}{{M}}{{{\iota }}_{n \times 3}}{\xi _{\text{f}}}{\omega _{\text{f}}}}\\{ - {{cl{ \varPhi}}}}& {{c}}& { - 2{{m}}{{{\iota }}_{{n_m} \times 3}}{\xi _{\text{f}}}{\omega _{\text{f}}}}\\{{O}}& {{O}}& {2{\xi _{\text{f}}}{\omega _{\text{f}}}{{{I}}_{3 \times 3}}}\end{array}} \right]$

(8)

${{{K}}_{{\text{ft}}}} = \left[ {\begin{array}{*{20}{c}}{{{\overline K}} + {{{\varPhi }}^{\text{T}}}{{{l}}^{\text{T}}}{{kl\varPhi }}}& { - {{{\varPhi }}^{\text{T}}}{{{l}}^{\text{T}}}{{k}}}& { - {{{\varPhi }}^{\text{T}}}{{M}}{{{\iota }}_{n \times 3}}\omega _{\text{f}}^2}\\{ - {{kl\varPhi}}}& {{k}}& { - {{m}}{{{\iota }}_{{n_m} \times 3}}\omega _{\text{f}}^2}\\{{O}}& {{O}}& {\omega _{\text{f}}^2{{{I}}_{3 \times 3}}}\end{array}} \right]$

(9)

${{{f}}_{bw}}{\text{ = }}{\left[{\begin{array}{*{20}{c}}{{{{O}}_{1 \times \left( {n + {n_i}} \right)}}}& { - {{{w}}^{\text{T}}}}\end{array}} \right]^{\text{T}}}$

(10)

Three components of earthquake excitation are not necessarily correlated in its principle directions. With $ {S\!\!_0}$ denoting the main direction intensity of earthquake excitation, Penzien and Watabe^[14] assigned the ratio of variance in three directions as 1: 0.75: 0.5, and therefore the excitation power spectrum density matrix $ { S}$ can be defined as Eq. (11).

${{S}}{\text{ = }}E \left( {{{w}}{{{w}}^{\text{T}}}} \right){\text{ = }}\left[ {\begin{array}{*{20}{c}}{{S\!\!_0}}& 0& 0\\0& {0.75{S\!\!_0}}& 0\\0& 0& {0.5{S\!\!_0}}\end{array}} \right]$

(11)

where E(•) is expectation operator.

Amplitude and duration are important characters of ground motion.In this paper, the amplitude was expressed as a function of peak ground acceleration ( $PGA $ ).Assuming that $PGA $ takes 3 times of acceleration variance $ {\sigma _w}$ ^[15], one has the expression of $ S_0$ as Eq. (12).

${S_0} = \frac{{2{\xi _{\text{f}}}{{\left( {PGA} \right)}^2}}}{{{3^2}{{\pi }}\left( {1 + 4\xi _{\text{f}}^2} \right){\omega _{\text{f}}}}}$

(12)

$ PGA$ is a function of earthquake magnitude ${M_{\rm s}} $ and focal distance $R $ ^[16–18] and was assumed to be expressed as Eq. (13).

$PGA{\text{ = }}{b_1}{{\rm{e}}^{{b_2}M_{\rm s}}}{\left( {R{\text{ + }}{R_0}} \right)^{ - {b_3}}}$

(13)

where $ b_1$ , $b_2 $ , $ b_3$ and $ R_0$ are parameters obtained from estimation over rather broad geographical regions.

The probability model of moment magnitude ${M_{\rm s}} $ can be described using truncated Gutenberg-Richter relationship^[19] in the interval of $ [{M_{\rm s}}_{{\min}}, {M_{\rm s}}_{{\max}}]$ , as Eq. (14).

$p\left( {M_{\rm s}} \right){\text{ = }}\frac{{{\beta _{M_{\rm s}}}{{\rm{e}}^{ - {\beta _{M_{\rm s}}}M_{\rm s}}}}}{{{{\rm{e}}^{ - {\beta _{M_{\rm s}}}{{M_{\rm s}}_{{\min} }}}} - {{\rm{e}}^{ - {\beta _{M_{\rm s}}}{{M_{\rm s}}_{{\max} }}}}}}$

(14)

where $\, \beta_{M_{\rm s}} $ is the regional seismicity factor and $ p(M_{\rm s})$ is the annual frequency.

The ground-motion duration $ t_{\text{d}}$ was assumed to be composed of the source duration and path-dependent duration $T_{\text{p}} $ ^[20], as Eq. (15).

${t_{\text{d}}}{\text{ = }}{T_{\text{p}}}{\text{ + }}{{0.5} /{{f_{\text{a}}}}}$

(15)

where $f_{\text{a}} $ is corner frequency.

Considering the uncertainty that associated in $\omega_{\text{f}} $ , $\xi_{\text{f}} $ , $ R$ , $ Ms$ , the hazard random variable set can be expressed as Eq. (16).

${{{\psi }}_H} = \left\{ {{\omega _{\text{f}}}, {\xi _{\text{f}}}, M_{\rm s}, R} \right\}$

(16)

Strictly speaking, the random variables should also include the bed rock white noise process. However, the structural response under random process can be solved by stochastic dynamics and eventually be expressed as a deterministic response variance. Therefore, it was not included in the random variable set.

3 Reliability Analyses

The state vector for the NC-MTMD system with the KT filter is ${{y}} = {[ {\begin{array}{*{20}{c}}{{q}}& {{x}}& {{{{x}}_{\text{f}}}}& {{{\dot q}}}& {{{\dot x}}}& {{{{{\dot x}}}_{\text{f}}}}\end{array}} ]^{\rm{T}}}$ and one can easily obtain the state function of the structure, as Eq. (17).

${{\dot {{y}}}} = {{Ay}} + {{B}}{{{f}}_{{\text{bw}}}}$

(17)

where A and B are state matrices given by

${{A}} = \left[ {\begin{array}{*{20}{c}}{{O}}& {{I}}\\{ - {{M}}_{{\text{ft}}}^{ - 1}{{{K}}_{{\text{ft}}}}}& { - {{M}}_{{\text{ft}}}^{ - 1}{{{C}}_{{\text{ft}}}}}\end{array}} \right], {{B}} = \left[ {\begin{array}{*{20}{c}}{{O}}\\{{{M}}_{{\text{ft}}}^{ - 1}}\end{array}} \right]$

(18)

The stationary response can be obtained by solving the Lyapunov equation, as Eq. (19).

${{AR}} + {{R}}{{{A}}^{\rm{T}}} + 2{{\pi }}\left[ {\begin{array}{*{20}{c}}{{O}}& {{O}}\\{{O}}& {{S}}\end{array}} \right] = {{O}}$

(19)

where $ {{R}}$ is the covariance matrix of $ { y}$ .

The response corresponding to certain state variable ${\textit{z}}_i (i=1, 2, …, n_{\textit{z}}) $ can be calculated as Eq. (20).

$\left\{\begin{array}{l}\sigma _{{{\textit{z}}_i}}^2 = {{n}}_i^{\text{T}}{{GR}}{{{G}}^{\text{T}}}{{{n}}_i}, \$6pt]\sigma _{{{\dot {\textit{z}}}_i}}^2 = {{n}}_i^{\text{T}}{{GAR}}{{{A}}^{\text{T}}}{{{G}}^{\text{T}}}{{{n}}_i}\end{array}\right. $

(20)

where $ \sigma _{{{\textit{z}}_i}}^2$ and $\sigma _{{{\dot {\textit{z}}}_i}}^2 $ are variance of variable ${\textit{z}}_{\text{i}} $ and its first derivative of time, respectively; $ { G}$ is observation matrix; ${{n}}_i $ is unit vector at the dimension of variable $ {\textit{z}}_{\text{i}}$ , therefore $ {{\textit{z}}_i} = {{n}}_i^{\rm T}{{z}}$ .

Take structure performance space $ \varPi $ as a $n_{\textit{z}} $ dimensions space with every dimension corresponding to a state variable and assume that every variable has a certain limit state bound $ \beta_i$ , one therefore has a hypercube space $ {\varPi _{\text{s}}}$ named safe polygon, in which all structure variables satisfies the limit state condition, as Eq. (21).

${\Pi _{\text{s}}} = \left\{ {{{\textit{z}}} \in R: \left| {{{\textit{z}}_i}} \right| < {\beta _i}} \right\}, i = 1, 2, \cdots, {n_{\textit{z}}}$

(21)

The first passage failure probability was expressed as Eq. (22).

${P_{\rm{f}}}\left( {{t_{\rm{d}}}} \right) = \int_0^T {P\left[ {{\textit{z}}\left( \tau \right) \notin {\Pi _{\text{s}}}} \right]} {\rm{d}}\tau {\text{ = }}1 - \exp \left( { - v_{\textit{z}}^ + {t_{\rm{d}}}} \right)$

(22)

where $v_{\textit{z}}^ + $ is the out-crossing rate and can be calculated from Eq. (23).

$\left\{\begin{aligned}& v_z^ + \approx \sum\limits_{i = 1}^{{n_z}} {{w_{{{\textit{z}}_i}}}\left( {{\lambda _{{{\textit{z}}_i}}}r_{{{\textit{z}}_i}}^ + } \right)}, \\ & {w_{{{\textit{z}}_i}}} = \int_{{{\text{B}}_i} \cap {\text{F}}} {p\left( {{{{z}}_ \bot }|{{\textit{z}}_i} = {\beta _i}} \right)} {\rm{d}}{{{z}}_ \bot }, \\ & r_{{{\textit{z}}_i}}^ + = \frac{{{\sigma _{{{\dot {\textit{z}}}_i}}}}}{{\pi {\sigma _{{{\textit{z}}_i}}}}}\exp \left\{ { - \frac{{\beta _i^2}}{{2\sigma _{{{\textit{z}}_i}}^2}}} \right\}, \\ & {\lambda _{{{\textit{z}}_i}}} \approx \frac{{1 - \exp \left\{ { - {q^{0.6}}{{\left( {\displaystyle\frac{2}{{\sqrt \pi }}} \right)}^{0.1}}\displaystyle\frac{{2\sqrt 2 }}{{{n_b}}}\displaystyle\frac{{{\beta _i}}}{{{\sigma _{{{\textit{z}}_i}}}}}} \right\}}}{{1 - \exp \left( { - {{\beta _i^2} / {2\sigma _{{{\textit{z}}_i}}^2}}} \right)}}, \\ & q = \frac{{\sigma _{{{\textit{z}}_i}}^5}}{{4\pi \int_{ - \infty }^{ + \infty } {\left| \omega \right|{S_{{{\textit{z}}_i}{{\textit{z}}_i}}}\left( \omega \right){\rm{d}}\omega \int_{ - \infty }^{ + \infty } {S_{{{\textit{z}}_i}{{\textit{z}}_i}}^2\left( \omega \right){\rm{d}}\omega } } }}\end{aligned}\right.$

(23)

where ${w_{{{\textit{z}}_i}}} $ is the correlation weighting factor and can be obtained by Eq. (23)^[21], ${{{z}}_ \bot } = {{z}} - {{{z}}_i}{{{n}}_i} $ , $ B_i$ is the limit state hyperplane bound of ${\textit{z}}_i $ and F is the surface of hypercube space ${\Pi _{\text{s}}} $ ; Eq. (23) is the Rice function^[22–23] and yields unconditional out-crossing rate; $ {\lambda _{{{\textit{z}}_i}}}$ is the out-crossing correction factor that calculated by Eq. (23)^[24], $ n_{\text{b}}$ is the number of limit state bound side; factor $q $ can be calculated by Eq.(23)^[21], $\omega $ is circular frequency and ${S_{{{\textit{z}}_i}{{\textit{z}}_i}}}\left( \omega \right) $ is the power spectrum density of structural response ${\textit{z}}_i$ .

Noting that the duration of earthquake is a function of hazard variables^[20], the unconditional failure probability can be obtained by the integration as Eq. (24).

${P_{\rm{f}}} = \int {\int {{p_{\rm{f}}}\left( {{t_{\rm{d}}}\left( {{\psi _{\rm{H}}}} \right){\rm{ /}}{t_d}\left( {{\psi _{\rm{H}}}} \right)\left\{ {{\psi _{\rm{S}}},{\psi _{\rm{H}}}} \right\}{\rm{ }}} \right)p\left( {{\psi _{\rm{S}}}} \right)p\left( {{\psi _{\rm{S}}}} \right){\rm{d}}{\psi _{\rm{H}}}{\rm{d}}{\psi _{\rm{H}}}} } $

(24)

Assuming failure event follows engineering Poisson distribution of independent occurrences^[25], the failure probability after a time period $ t_{\text{l}}$ can be obtained by Eq. (25).

$\left\{\begin{array}{l}{P_{\rm F}}\left( {{t_{\text{l}}}} \right) = 1 - {{\rm e}^{ - {t_{\text{l}}}{v_{\rm f}}{P_{\rm f}}}}, \\{v_{\rm f}} = {{\rm e}^{{\alpha _{M_{\rm s}}} - {\beta _{M_{\rm s}}}M_{\rm s}{_{{\min} }}}} - {{\rm e}^{{\alpha _{M_{\rm s}}} - {\beta _{M_{\rm s}}}M_{\rm s}{_{{\max} }}}}\end{array}\right.$

(25)

where $\alpha _{M_{\rm s}} $ and $\beta_{M_{\rm s}} $ are regional seismicity factors.

4 Optimization

With given limit state bounds, the objective function that adopted in this paper was defined as $ P_{\text{F}}$ .With design variables as ${{\theta}} $ , the optimization problem can be described as a constrained optimization problem as Eq. (26).

$\left\{\begin{array}{l}{\rm min}\;{P_F({\theta})}, \\{{\theta }}{\text{ = }}\left\{ {\begin{array}{*{20}{c}}{{\omega _1}, {\omega _2}, \cdots , {\omega _{{n_m}}}}& {{\xi _1}, {\xi _2}, \cdots , {\xi _{{n_m}}}}\end{array}} \right\}, \\{{\theta }} \in \mathop {I}\limits_{j = 1}^{{n_m}} {\left\{ {\left\{ {{\omega _j}: {\omega _j} \in \left[ {{\omega _{j{\min} }}, {\omega _{j{\max} }}} \right]} \right\} \cap \left\{ {{\xi _j}: {\xi _j} \in \left[ {{\xi _{j{\min} }}, {\xi _{j{\max} }}} \right]} \right\}} \right\}} \end{array}\right.$

(26)

where ${\omega _j}{\text{ = }}\sqrt {{{{k_j}} / {{m_j}}}} $ and ${\xi _j}{\text{ = }}{{{c_j}} /{\left( {2{m_j}{\omega _j}} \right)}} $ .

5 Case study

The case considered here is a concentrically braced steel thermal power plant building (Fig. 2). The structure consists of boiler frames, air heater houses and a scuttle bay. In its scuttle bay, 7 scuttles with each weight 1 098.3 tons locate at 32.2 m height of the structure.

5.1 Deterministic model

The structure was located at a Chinese site with relatively high hazard level. The expected bounds of earthquake magnitude are ${M_{\rm s}}_{\rm min} $ =5.5 and ${M_{\rm s}}_{\rm max} $ =8, respectively. The regional seismicity factors $\alpha _{M_{\rm s}} $ and $ \beta_{M_{\rm s}}$ are 4ln(10) and 2.16^[16], respectively.In Eq.(13), $b_1 $ =6.63; $ b_2$ =1.17; $b_3 $ =-1.43; $R_0 $ =14 for $PGA $ in cm/s^2[18]. In Eq. (15) $ f_{\text{a}}=10^{2.181-0.496M_{\rm s}}$ ^[26] and $T_{\text{p}}=0.05R $ ^[27]. The site was 40 km away from the focal center. The soil condition fits class III (stiff soil) with predominant period of 0.55 s. The KT filter parameters were assigned as $\xi_{\text{f}} $ = 0.6^[28] and $\omega _{\text{f}} = 3.63{\text{π}}$ rad/s (the corresponding circular frequency of the predominant period).

The main structure was modeled by its first 12 modes, at which the cumulated participation factor at two horizontal directions is larger than 90% (Tab. 1). Rayleigh damping was adopted to model the inherent damping of the main structure, with 2% at 0.1 s and 1 s^[29], respectively. Scuttles were modeled as a 2DOF lumped mass each, with a DOF at each horizontal direction.

5.2 Probabilistic model

Parameters in random variable sets $\varPsi_{\text{S}} $ and $\varPsi_{\text{H}} $ were modeled as random variables with probabilistic distributions in Tab. 2. Random elements in each matrices ( $ {{\overline M}}$ and ${{\overline K}} $ ) were assumed to be perfectly correlated to reduce the prohibitive computational burden. Scuttle mass is modeled by Beta distribution for its bounded character.Mass, stiffness and damping of scuttles were assumed to be independent with each other. Besides specific probabilistic distributions, all the rests were using the Lognormal distribution for its wide application and non-negative domain (Fig. 3).

Latin hypercube sampling (LHS) strategy was adopted in this paper to reduce the needed sample size. Appropriate sample size investigation was further carried out to determine the necessary sample size.With design variables $\omega_j=\{0.3\omega_{j{\max}}+0.7\omega_{j{\min}}, 0.7\omega_{j{\max}}+0.3\omega_{j{\min}}\} $ and $ \xi_j=\{0.3\xi_{j{\max}}+ 0.7\xi_{j{\min}}, 0.7\xi_{j{\max}}+0.3\xi_{j{\min}}\}$ , 128 combinations/cases were generated, with the consideration of computational burden and representativeness. Normalized probability, which is defined as the ratio between the probabilities $P_{\text{F}}^{\text{e}}$ estimated by using certain number and $ P_{\text{F}}^\infty $ by all of the samples, was used to assess the estimation accuracy of using the certain number of samples. The appropriate sample size is required to have $P_{\text{F}}^{\text{e}}/P_{\text{F}}^\infty $ of all cases in [90%, 110%]. From Fig. 4, it can be found that 4 699 samples are enough. Therefore, 4 699 samples were generated in the following studies.

5.3 Statement of optimization problem

Considering that a too high dimensional structural output may cause a prohibitive computational burden, only eight structural outputs that expected to be the most critical, i.e., drifts of corner columns at the 1st and 3th floor (shown in Fig. 5), were considered, which are expected to have higher value^[30].

Denoting $j={1, 2, \cdots, n_{{m}}} $ in Eq.(26) as {A_x, A_y, …, G_y} (Fig.2(b)), the optimization problem can be expressed as Eq.(27).

$\left\{\begin{array}{l}\!\!\!\! {\min}\;{P_{\text{F}}({{\theta}})}\\\!\!\!\! {{\theta }}{\text{ = }}\left\{ {\begin{array}{*{20}{c}}{{\omega _{{A_ x}}}, {\omega _{{ A_y}}}, \cdots , {\omega _{G_y}}}& {{\xi _{{A_x}}}, {\xi _{{A_y}}}, \cdots , {\xi _{{G_y}}}}\end{array}} \right\},\\\!\!\!\! {{\theta }} \in \!\!\!\!\!\!\!\!\! \mathop {I}\limits_{j = \{ {{\text{A}}_x}, {\text{ }}{{\text{A}}_y}, \cdots , {\text{ }}{{\text{G}}_y}\} }\!\!\!\!\! {\left\{ {\left\{ {{\omega _j}\!: \!{\omega _j} \!\in \!\left[ {{\omega _{j{\min} }}, {\omega _{j{\max} }}} \right]} \right\}\! \cap \!\left\{ {{\xi _j}\!: \!{\xi _j} \!\in\! \left[ {{\xi _{j{\min} }}, {\xi _{j{\max} }}} \right]} \right\}} \right\}} \end{array} \right.$

(27)

where ${\omega _{j{\min} }} $ , $ {\omega _{j{\max} }}$ , $ {\xi _{j{\min} }}$ , $ {\xi _{j{\max} }}$ takes 0, 17 rad/s, 0, 0.3, respectively^[6].

5.4 Optimum design

The problem contains uncertainty and deterministic algorithms (e.g., the gradient-based algorithms and direct search algorithms) are not fit for this situation. The genetic algorithm, which is one of the most popular stochastic optimization algorithms, was adopted in this paper to solve the optimization problem.

The optimum design obtained by the optimization process was summarized in Tab. 3.

It can be seen that $ \omega_{\text{opt}}$ varies over the range of natural periods of the main structure, which highlights the importance of higher modes. There is no obvious trend with $\xi_{\text{opt}} $ . But comparing with the results in^[6], which assumed no uncertainty, most $ \xi_{\text{opt}}$ here in Tab. 3 are higher. This is because that higher damping ratio enhances robustness of MTMD^[31] and therefore benefits its performance when uncertainty presents.

6 Parametric studies

Uncertain mass of the scuttles caused by changing coal storage cause large variation on oscillator frequencies.The variation can further cause problem on tuning. Pendulum system has an independent period with system mass, because of the perfect correlation between stiffness and mass. Consequently, the pendulum strategy is potentially a solution for the mass uncertainty. Another important practical problem is collision between coal scuttles and its surrounding structural members. These two important practical aspects were therefore investigated in this section.

6.1 Influence of oscillator type

For translational MU-NC-MTMD system (Fig. 6), stiffness of sub-oscillators was induced by elastic potential energy of isolators. As an alternative, the pendulum MU-NC-MTMD, in which stiffness of sub-oscillators is the function of gravity potential energy of oscillators, can possibly be effective to enhance the system robustness when mass uncertainty presents.

For pendulum MU-NC-MTMD, all modeling methods are same as translational MU-NC-MTMD, but stiffness. For pendulum MU-NC-MTMD, stiffness was expressed as Eq. (28).

$k = {\rm{diag}}\left\{ {{m_1}g/{\rho _1},{m_2}g/{\rho _2}, \cdots ,{m_{{n_m}}}g/{\rho _{{n_m}}}} \right\}$

(28)

where $g $ and $\rho $ denote gravity acceleration and length of pendulums, respectively.

1/ $\rho $ of each scuttle was assumed to follow lognormal distribution, with the $CoV $ of 0%, 5% and 10% for three different cases, respectively.

It can be observed from Tab. 4 that the translational system has slightly smaller $P_{\text{F}} $ than the pendulum system. Furthermore, $P_{\text{F}} $ of pendulum system is not sensitive to the uncertainty of 1/ $\rho $ .

6.2 Influence of seismic gap

The gap between the coal scuttle and the surrounding structural members is 130 mm. Considering possible collision in {A_x, A_y, $ \cdots $ , G_y} directions, 14 additional relative displacement limit state bounds were included.

For the optimum design in Tab. 3, probability of failure with the consideration of collision $ P_{\text{Fc}}$ was calculated, which equals 8.11%. The consideration of collision slightly increased the failure probability, which suggests that collision problem is not that critical. It is also possible to optimize $P_{\text{Fc}}$ to obtain a design that mitigates both $P_{\text{F}} $ and collision problem.

The author adopted multi-objective method^[5] in a previous engineering problem. The multiple output unconditional PF integrates multiple objectives into one failure probability and therefore avoids cumbersome multi-objective optimization.

7 Conclusion

In this paper, a multiple output unconditional reliability-based design framework for nonconventional multiple tuned mass damper for a complex structure was proposed. With the description of simplified structure model and earthquake hazard model, a multiple output reliability method was presented. Optimization problem was then formulated. A case of coal scuttles isolation for thermal power plant was solved to illustrate the method. From the study, conclusions as follows can be drawn:

1) The framework adopts failure probability with multiple limit state bounds as the objective function and therefore avoids the cumbersome multi-objective optimization. Multiple responses are integrated into one failure probability, which is clearer.

2) The pendulum system does not perform better than translational system in the case with equal uncertainty. Performance of pendulum system is not sensitive to the degree of its curvature uncertainty.

3) Collision problem is not critical for this case but it does affect the failure probability. Failure probability with the consideration of collision can be further optimized to search a design solution that mitigates both structural failure probability and collision problem.