Direct Displacement-based Seismic Design of Glulam Frames with Buckling Restrained Braces

ABSTRACT

This paper presents a direct displacement-based design (DDBD) approach for the buckling restrained braces (BRBs) braced glue-laminated timber (glulam) frame (BRBGF) structures. First, the critical design parameters of the DDBD approach were derived for BRBGFs. Then, using experimentally verified numerical models, pushover analyses and nonlinear time-history analyses (NLTHAs) were conducted on a series of one-storey BRBGFs to calibrate the stiffness adjustment factor λ for BRB-timber connections and the spectral displacement reduction factor η. Finally, the DDBD approach was verified as a prospective approach for the seismic design of multi-storey BRBGF buildings by NLTHAs of the case study buildings.

KEYWORDS:

• Glulam frames
• buckling restrained braces (BRBs)
• direct displacement-based design (DDBD)
• BRB-timber connections
• nonlinear time-history analyses (NLTHAs)

1. Introduction

Multi-storey mass timber buildings have gained immense popularity worldwide due to the development of engineered wood products with enhanced structural performance compared to sawn timber. These timber buildings also hold potential for generalizing climate-change mitigation and preservation (Himes and Busby 2020). However, timber has a relatively low elastic modulus compared to other construction materials such as steel and reinforced concrete. In addition, timber members are prone to brittle failure under tension with a limited energy dissipation capacity. Larger member sizes are often needed to satisfy stiffness requirements, and limited ductility is usually assumed for multi-storey mass timber buildings (Kirstein, Siracusa, and Smith 2018). In earthquake-prone countries such as New Zealand, seismic considerations usually govern the design of lateral force resisting systems (LFRS) of multi-storey buildings. The limited ductility assumption may result in uneconomical member sizes and increase the amount of LFRS (e.g. shear walls and braces). Excess shear walls or braces may restrict the flexibility of architectural plans and reduce space efficiency.

Several new solutions have been proposed to overcome these limitations. One is to introduce hybrid systems that can take advantage of the strengths of different construction materials and provide suitable solutions for LFRS (Brown et al. 2021; Dickof et al. 2014; Dong et al. 2020a; He et al. 2014, 2018). Quintana Gallo, Carradine, and Bazaez (2021) conducted a comprehensive review of timber-based hybrid LFRS, such as steel frames with infilled cross-laminated timber (CLT) shear walls (Dickof et al. 2014), steel frames with light-timber frame shear walls (Luo et al. 2021), post-tensioned CLT shear walls with energy dissipators (Brown et al. 2021, 2022) and timber frames with energy dissipators (Gilbert and Erochko 2019; Song et al. 2022). The strength, stiffness, ductility and energy dissipation capacity were improved through the hybridisation. However, timber shear walls might be more suitable for residential buildings with fewer demands on open spaces. In addition, the post-tension loss of mass timber buildings during long-term service still requires more investigation (Granello et al. 2020). Timber frames with energy dissipators, providing larger open spaces and relatively accessible replacements for dissipators, are still more favoured for many projects. Recently, a new timber-steel hybrid system that integrates buckling restrained braces (BRBs) into glued laminated timber (glulam) frames was proposed and experimentally investigated (Dong et al. 2020b). As shown in Fig. 1, dowelled connections with inserted steel plates were used to connect the glulam frame with the gusset plates, and BRBs were connected to the gusset plates by pinned connections. BRBs were designed as the energy dissipating members in the system. All timber members and connections were protected by the capacity design to remain elastic before BRBs reached the maximum designed deformation (Dong et al. 2020b). Two full-scale specimens with two different BRB–timber interface connections were tested under cyclic loading. The test results showed that the BRB-braced glulam frames (BRBGFs) significantly increased energy dissipation capacity compared with conventional timber-braced glulam frames. In addition, the damage in the connections and glulam members was minimised as the capacity design approach was followed. As shown in Fig. 2, initial slips were observed on the load-drift curves of the BRBGFs because the specimens used pin-end BRBs (Metelli, Bregoli, and Genna 2016) and dowelled connections that required manufacturing tolerances to be installed. The experimental test data were then used to calibrate component-based numerical models in OpenSees and the modelling details are explained in Dong et al. (2021). The influence of connection stiffness and initial slips on the cyclic performance of the system was investigated through parametric studies. It was shown that the connections had a negligible moment-resisting capacity, and the connections designed with a connection relative overstrength factor γ = 1.5 were suitable for engaging the BRBs and protecting the glulam members and connections. In addition, the initial slips caused by the manufacturing tolerances had a negligible impact on the ultimate strength and overall energy dissipation of the BRBGFs under cyclic loading. Further information regarding previous experimental design and numerical studies can be found in Dong et al. (2020b) and Dong et al. (2021), respectively.

Figure 1. BRBGF specimen with the dowelled connections (Dong et al. 2020b).

Figure 2. Comparison of BRBGF hysteresis curves (Dong et al. 2021).

Fig. 2 shows that the numerical models of the BRBGFs can accurately predict the cyclic performance. However, the complicated modelling process and lack of necessary design parameters of this hybrid system may still restrict its practical application. In addition, the influence of initial slips under seismic ground motions requires further investigation. This study fills the research gaps for applying a direct displacement-based design (DDBD) approach (Priestley, Calvi, and Kowalsky 2007) to the seismic design of BRBGF structures. It provides suitable displacement profiles, stiffness adjustment for timber–steel interface connections and hysteretic modelling for DDBD. In addition, it verifies the suitability of the approach via nonlinear time-history analyses (NLTHAs) of case study buildings, and thus facilitates the application of this hybrid system.

2. Fundamentals of the DDBD Approach

The DDBD approach was first proposed by Kowalsky, Priestley, and MacRae (1994) and further developed by Priestley (1995, 2003) to provide an alternative design solution for the traditional force-based design (FBD) approach widely used by modern design standards (e.g. Eurocode 8 (2005a) and NZS 1170.5 (2004)). Previous research (Maley, Sullivan, and Corte 2010; Priestley, Calvi, and Kowalsky 2007; Sullivan 2013) has shown that the FBD approach may have limitations for seismic design. One limitation is that the displacement response, usually directly correlated with building damage levels (Erduran and Yakut 2004; Ghobarah, Abou-Elfath, and Biddah 1999; Medina and Krawinkler 2005), is given only secondary importance. Alternative methods have been developed to address some of the limitations of traditional FBD approaches. For example, in the New Zealand standard NZS 1170.5 (2004), the designer is prompted to satisfy structural and non-structural deformation limits and check that the design ductility is compatible with the expected ductility demands. However, these standards do not require engineers to address all the issues raised by Priestley, Calvi, and Kowalsky (2007). For example, elastic analysis is used to calculate the inelastic force distributions in the FBD approach and the unique force reduction factor for a given structure type may be invalid (Priestley, Calvi, and Kowalsky 2007). Previous research has shown that the DDBD approach might benefit some conditions more than the FBD approach for nonlinear seismic design (Salawdeh and Goggins 2016; Terán-Gilmore and Ruiz-García 2011). Consequently, the DDBD approach has been used in a variety of structural systems, including reinforced concrete structures (Belleri 2017; Pennucci, Calvi, and Sullivan 2009; Priestley and Kowalsky 2000; Sullivan, Priestley, and Calvi 2006; Yang and Lu 2018), steel structures (Malla and Wijeyewickrema 2021; Nievas and Sullivan 2015; Roldán, Sullivan, and Della Corte 2016; Sahoo and Prakash 2019; Sullivan 2013; Wijesundara 2012) and timber structures (Chen et al. 2021; Filiatrault and Folz 2002; Hashemi, Zarnani, and Quenneville 2020; Li et al. 2019; Pang and Rosowsky 2009; Thiers-moggia and Málaga-Chuquitaype 2021; Zonta et al. 2011). A model code for the DDBD approach was also developed for practical design guidelines (Sullivan, Priestley, and Calvi 2012). For BRBGF structures, the initial slips, as shown in Fig. 2, may cause difficulty in estimating the initial stiffness and fundamental period essential for the FBD approach. Therefore, the FBD approach may not be well suited to particular cases of BRBGFs. The DDBD approach implicitly assumes that a low initial stiffness will not affect the displacement demands in high-intensity earthquakes by utilising the secant stiffness concept (a hypothesis that will be verified in this work). In this regard, the DDBD approach is developed for the BRBGF structures.

The primary process of the DDBD approach (Priestley, Calvi, and Kowalsky 2007) is illustrated in Fig. 3. The first step is to substitute a multi-degree-of-freedom (MDOF) structure with an equivalent single-degree-of-freedom (SDOF) system (Fig. 3a). A displacement profile of the MDOF structure under seismic loads is assumed based on the knowledge of the design displacement profile. Therefore, the equivalent SDOF system characterised by its design displacement Δ_d, effective mass M_e, and effective height H_e can be calculated by Eq. 1(1) ${\mathrm{\Delta }}_d=\frac{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i^2}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}$(1) -Eq. 3(3) $H_e=\frac{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i{H}_i}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}$(3) .

(1) ${\mathrm{\Delta }}_d=\frac{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i^2}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}$(1)

(2) $M_e=\frac{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}{{\mathrm{\Delta }}_d}$(2)

(3) $H_e=\frac{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i{H}_i}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}$(3)

Figure 3. DDBD design approach (Priestley, Calvi, and Kowalsky 2007).

where Δ_i is the storey drift at the i-th storey (mm); m_i is the seismic mass at the i-th storey (tons); H_i is the height of the i-th storey from the ground (m); n is the total number of storeys.

The second step is to calculate the equivalent SDOF system ductility demand μ at the design displacement Δ_d according to Eq. 4(4) $\mu =\frac{{\mathrm{\Delta }}_d}{{\mathrm{\Delta }}_y}$(4) , with the design and yield displacement points shown in Fig. 3b. To quantify the effects of nonlinear behaviour on the displacement demands, Priestley, Calvi, and Kowalsky (2007) recommended characterising the equivalent SDOF system with an equivalent viscous damping ratio ξ_eq, which is a function of μ, as shown in Fig. 3c. ξ_eq is computed using an empirical ξ_eq-μ relationship calibrated to the results of NLTHA, as will be discussed further in Section 3.3. The design displacement spectrum of an SDOF system is scaled to ξ_eq, as shown in Fig. 3d, using the spectral displacement reduction factor η, a function of ξ_eq. Thus, the required T_e for the equivalent SDOF system can be obtained according to Δ_d. Then, K_e and the design base shear force V_d can be determined using Eq. 5(5) $K_e={\left(\frac{2\pi }{T_e}\right)}^2{M}_e$(5) and Eq. 6(6) $V_d=K_e{\mathrm{\Delta }}_e$(6) , respectively, when the P–Δ effects are not significant.

(4) $\mu =\frac{{\mathrm{\Delta }}_d}{{\mathrm{\Delta }}_y}$(4)

(5) $K_e={\left(\frac{2\pi }{T_e}\right)}^2{M}_e$(5)

(6) $V_d=K_e{\mathrm{\Delta }}_e$(6)

The design base shear force V_d for the equivalent SDOF system is also the design base shear force for the MDOF structure, so the third step is to distribute V_d along the height of the structure as storey forces F_i using Eq. 7(7) $F_i=\frac{m_i{\mathrm{\Delta }}_i}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}{V}_d$(7) . These forces are used to design the structural members in MDOF structures, with possible adjustments to the force profile for taller buildings to mitigate higher mode effects (as clarified later in Section 3.2).

(7) $F_i=\frac{m_i{\mathrm{\Delta }}_i}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}{V}_d$(7)

3. Extending the DDBD Approach to the BRBGF Structures

To apply the DDBD approach to BRBGF structures, the following knowledge is required: (1) the displacement profile and the design displacement at the performance limit state for BRBGF structures when substituting BRBGFs into an equivalent SDOF system (Fig. 3a); (2) μ for BRBGF structures (Fig. 3c), requiring estimation of the yield displacement; and (3) the relationship between ξ_eq and μ to obtain T_e (Fig. 3c and d). In this section, the determination of these design parameters is discussed.

3.1. Displacement Profile and Limit State Displacement

The design displacement Δ_d of the equivalent SDOF system depends on the assumed displacement profile and the limit state displacement of the MDOF structure, as shown in Eq. 1(1) ${\mathrm{\Delta }}_d=\frac{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i^2}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}$(1) . Past research has shown that the empirical displacement profile for moment-resisting frame (MRF) structures proposed by Priestley, Calvi, and Kowalsky (2007) is also suitable for concentrically braced frame (CBF) structures (Pettinga and Priestley 2005; Wijesundara and Rajeev 2012). Rajeev et al. (2017) verified the feasibility of the empirical displacement profile based on numerous NLTHAs of steel CBF structures; Maley, Sullivan, and Corte (2010) also found that it worked well for dual BRB-MRF systems. Therefore, the empirical displacement profile shown in Fig. 4 was used for BRBGF structures. The first storey has the largest inter-storey drift ratio, so the performance limit state displacement is governed by the first storey design inter-storey drift ratio θ_d. The empirical design displacement profile can be expressed as a function of θ_d as shown in Eq. 8a(8a) ${\mathrm{\Delta }}_i=\left\{\begin{array}{cc}{\theta }_d{H}_i& n\le 4\\ {\omega }_{\theta }{\theta }_d{H}_i\frac{\left(4H_n-H_i\right)}{\left(4H_n-H_1\right)}& n\ge 4\end{array}\right.$(8a) according to Sullivan, Priestley, and Calvi (2012). Suppose more accurate displacement profiles of multi-storey BRBGF structures from shake table testing or building monitoring are available in the future. In that case, the displacement profile assumption can be improved, but the entire process of the DDBD approach presented in this study could remain the same.

Figure 4. Assumed displacement profile.

Previous research illustrated that higher mode effects might increase the peak storey drift demands (Gardiner, Clifton, and MacRae 2013). Sullivan (2013) suggested reducing the first mode design displacement profile to account for the additional displacements caused by the higher mode effects. Therefore, the design drift reduction factor ω_θ (Eq. 8b)(8b) ${\omega }_{\theta }=\{\begin{array}{cc}1.0& n\le 6\\ 1-0.015\left(n-6\right)& 6<n\le 16\end{array}$(8b) adapted by Sullivan, Priestley, and Calvi (2012) was included in the displacement profile to account for the higher mode effects of BRBGF structures.

(8a) ${\mathrm{\Delta }}_i=\left\{\begin{array}{cc}{\theta }_d{H}_i& n\le 4\\ {\omega }_{\theta }{\theta }_d{H}_i\frac{\left(4H_n-H_i\right)}{\left(4H_n-H_1\right)}& n\ge 4\end{array}\right.$(8a)

(8b) ${\omega }_{\theta }=\{\begin{array}{cc}1.0& n\le 6\\ 1-0.015\left(n-6\right)& 6<n\le 16\end{array}$(8b)

The design inter-storey drift ratio θ_d depends on the performance requirements of both structural and non-structural elements. Recent research on low-damage non-structural elements showed that the improved design of these non-structural elements could sustain a 2.0–2.5% drift ratio with minor damage (Dhakal et al. 2016). As the performance limits of non-structural elements are out of the scope of this study, the performance limits were determined by the structural elements. The serviceability performance limit of θ_d was set to 0.33% based on New Zealand engineering practice (Hashemi, Zarnani, and Quenneville 2020). The BRBGF tests (Dong et al. 2020b) showed that the residual drift ratio could be over 0.5% after loading to a 1.0% drift ratio, and the residual drift ratio of 0.5% was suggested as the permissible residual drift ratio for safety by McCormick et al. (2008), so the repairable damage limit state of θ_d was set to 1.0%. The ultimate limit state (ULS) drift ratio was set to 2.0% conservatively because the BRBGF tests proved that the BRBGF structures could sustain this drift limit without collapse or significant damage to timber members and connections (Dong et al. 2020b). A higher θ_d may be adopted in the future if it is verified through further research and testing.

3.2. Ductility Factor

The ductility factor μ of the equivalent SDOF system is used to estimate the equivalent viscous damping ξ_eq, as shown in Fig. 3c, before obtaining the effective period T_e. The yield displacement Δ_y is required to determine μ, as shown in Eq. 4(4) $\mu =\frac{{\mathrm{\Delta }}_d}{{\mathrm{\Delta }}_y}$(4) . The following sub-sections discuss the determination of Δ_y and links the ductility factor of each storey of BRBGF structures with its equivalent SDOF system.

3.2.1. Yield Drift of One-storey BRBGFs

Numerical models of the one-storey BRBGFs with dowelled connections were built in the OpenSees, as shown in Fig. 5. The modelling methods in Dong et al. (2021) was followed and briefly explained below. The experimental tests showed that the BRB–timber interface connections had negligible moment-resisting capacity (Dong et al. 2020b), so their rotational stiffness was set to zero. The translational stiffness of the interface connections was modelled using Pinching4 models in OpenSees with ZeroLength elements between nodes 5 and 6, nodes 1 and 7, and nodes 3 and 8. As shown in Fig. 2, the initial slips were simulated by ElasticMultiLinear model in OpenSees as a spring in series between nodes 5 and 6. The stiffness of the ElasticMultiLinear model was low when the drift was smaller than the specified initial slips, but the stiffness would be significantly higher than that of BRBs when the drift exceeded the initial slips. All glulam beams and columns were simulated by elasticBeamColumn model because they were protected by the capacity design. The detailed parameters for the interface connections and the experimental verification process can be found in Dong et al. (2021). Several methods have been used to define the yield drift (Priestley, Calvi, and Kowalsky 2007). The definition from Park (Park 1989) was used here because it could include the influence of the initial slips that might exist in BRBGF structures. In addition, the two well-defined linear parts of the backbone curves of the BRBGF specimens fitted the definition well (Dong et al. 2020b). The yield drift, δ_y,s, is defined as the lateral displacement corresponding to the yield strength, F_y, calculated using Eq. 9(9) $F_y=2{\varphi }_m{f}_y{A}_{c}cos\alpha$(9) and is equal to the lateral force when both BRBs yield. Pushover analyses were conducted for the one-storey BRBGFs with and without the initial slips; Fig. 6 shows the load-drift curves. The curves illustrate that the yield drift with initial slips δ_y,s was the yield drift without initial slips δ_y,0 plus initial slips δ_s, as shown in Eq. 10(10) ${\delta }_{y,s}={\delta }_{y,0}+{\delta }_s$(10) .

(9) $F_y=2{\varphi }_m{f}_y{A}_{c}cos\alpha$(9)

(10) ${\delta }_{y,s}={\delta }_{y,0}+{\delta }_s$(10)

Figure 5. Numerical model of the one-storey BRBGF.

Figure 6. Yield drift definition.

where ϕ_m is the material overstrength factor for the steel core of BRB; f_y is the expected steel yield strength (MPa); A_c is the area of the yield zone of steel core of BRB (mm²); α is the inclination angle of the BRBs as shown in Fig. 5.

3.2.2. Yield Drift of Multi-storey BRBGFs

Past research (Della Corte 2006; Maley, Sullivan, and Corte 2010; Sullivan, Priestley, and Calvi 2012; Wijesundara and Rajeev 2012) has illustrated that the inter-storey yield drift at the i-th storey δ_y,i in a multi-storey building may need to consider the inter-storey drift caused by the BRB deformation δ_y,i,BRB, and the inter-storey drift caused by the column axial deformation δ_y,i,col, as shown in Fig. 7. Similar to the one-storey BRBGF, the inter-storey drift contributions of the initial slips δ_s,i need to be included, so δ_y,i was calculated using Eq. 11(11) ${\delta }_{y,i}={\delta }_{y,i,BRB}+{\delta }_{y,i,col}+{\delta }_{s,i}$(11) .

(11) ${\delta }_{y,i}={\delta }_{y,i,BRB}+{\delta }_{y,i,col}+{\delta }_{s,i}$(11)

Figure 7. Yield drift components (a) the inter-storey drift contribution of the BRB deformation δ_y,i,BRB; (b) the inter-storey drift contribution of the column axial deformation δ_y,i,col.

The component δ_y,i,BRB was estimated using Eq. 12a(12a) ${\delta }_{y,i,BRB}=h_i{\theta }_{y,i,BRB}=h_i\frac{2{\epsilon }_y}{\lambda sin2\alpha }$(12a) based on research from Sullivan, Priestley, and Calvi (2012). In BRB steel frames, the connections between BRBs and steel frames are usually translationally rigid. The frame lateral stiffness is equal to the lateral stiffness of BRBs, whereas, in BRBGF structures, the connection stiffness should be considered. Stiffness is considered by introducing a stiffness adjustment factor λ, as shown in Eq. 12b(12b) $\lambda =\frac{k_{BRBGF}}{k_{BRBGF,rigid}}$(12b) . λ defines the lateral stiffness ratio between the BRBGF with translationally semi-rigid dowelled connections and the BRBGF with translationally rigid connections; therefore, λ should be between 0 and 1. The BRBs were modelled by Steel4 model in OpenSees as a truss, as shown in Fig. 5; thus the stiffness of BRBs was represented by an effective stiffness, K_eff,BRB, as shown in Eq. 13a(13a) $K_{eff,BRB}=E_s\frac{A_c{A}_{tr}{A}_e}{A_c{A}_{tr}{l}_e+A_c{A}_e{l}_{tr}+A_{tr}{A}_e{l}_c}=E_s{f}_{sm}\frac{A_c}{l_{wp}}$(13a) dependent on the geometry of the yield zone, transition zone and elastic zone of BRBs according to Vigh and Zsarnóczay and Balogh (2017). To simplify the calculation, the additional stiffness caused by the transition zone and the elastic zone was considered by a stiffness modification factor f_sm that amplified the elastic modulus of steel core E_s as shown in Eq. 13b(13b) $E_{eff,BRB}=f_{sm}{E}_s$(13b) according to Vigh, Zsarnóczay, and Balogh (2017). The detailed parameters for Steel4 can be found in the Appendix of Dong et al. (2021). In terms of that, the yield strain ε_y in Eq. 12a(12a) ${\delta }_{y,i,BRB}=h_i{\theta }_{y,i,BRB}=h_i\frac{2{\epsilon }_y}{\lambda sin2\alpha }$(12a) was calculated using Eq. 12c(12c) ${\epsilon }_y=\frac{{\varphi }_m{f}_y}{f_{sm}{E}_s}$(12c) to include the influence of steel core material (f_y,ϕ_m, and E_s) and BRB geometry (f_sm).

(12a) ${\delta }_{y,i,BRB}=h_i{\theta }_{y,i,BRB}=h_i\frac{2{\epsilon }_y}{\lambda sin2\alpha }$(12a)

with

(12b) $\lambda =\frac{k_{BRBGF}}{k_{BRBGF,rigid}}$(12b)

(12c) ${\epsilon }_y=\frac{{\varphi }_m{f}_y}{f_{sm}{E}_s}$(12c)

where h_i is the storey height of the i-th storey (m); and θ_y,i,BRB is the inter-storey drift ratio of the i-th storey caused by BRB deformations; ε_y is the BRB yield strain; f_sm is the stiffness modification factor as defined in Eq. 13; E_s is the elastic modulus of steel core (GPa).

(13a) $K_{eff,BRB}=E_s\frac{A_c{A}_{tr}{A}_e}{A_c{A}_{tr}{l}_e+A_c{A}_e{l}_{tr}+A_{tr}{A}_e{l}_c}=E_s{f}_{sm}\frac{A_c}{l_{wp}}$(13a)

(13b) $E_{eff,BRB}=f_{sm}{E}_s$(13b)

(13c) $l_{wp}=l_e+l_{tr}+l_c$(13c)

where K_eff,BRB is the BRB effective stiffness (kN/mm); A_c,A_tr, and A_e are the area (mm²) of the yield zone, transition zone and elastic zone of the BRB, respectively; l_c,l_tr, and l_e are their corresponding length (mm) (Zsarnóczay 2013); f_y is the expected steel yield strength.

Figure 8 shows a spring analogy model for the BRBGF, similar to the model proposed by Mahjoubi and Maleki (2016). As the top connection carried the load approximately twice as much as the bottom connections in the BRBGFs, the fastener number in the bottom connections was assumed to be half of that of the top connections. In this regard, the translational stiffness of the bottom connection k_con,b was approximately half of that of the top connections k_con.t in the model (i.e. k_con.t = 2k_con.b). To simplify the spring analogy model further, a stiffness ratio ν between k_con,b and the BRB lateral stiffness k_BRBGF,rigid is defined as Eq. 14a(14a) $\nu =\frac{k_{con,b}}{k_{BRBGF,rigid}}$(14a) . k_BRBGF,rigid represents the BRBGF stiffness with the assumption of translationally rigid connections; k_BRBGF can be calculated using Eq. 14b(14b) $\frac{1}{k_{BRBGF}}=\frac{2}{k_{con,b}}+\frac{1}{2k_{con,b}}+\frac{1}{k_{BRBGF,rigid}}=\left(\frac{2\nu +5}{2\nu }\right)\frac{1}{k_{BRBGF,rigid}}$(14b) considering the contributions from the translationally semi-rigid top and bottom connections. By comparing Eq. 12b(12b) $\lambda =\frac{k_{BRBGF}}{k_{BRBGF,rigid}}$(12b) and Eq. 14b(14b) $\frac{1}{k_{BRBGF}}=\frac{2}{k_{con,b}}+\frac{1}{2k_{con,b}}+\frac{1}{k_{BRBGF,rigid}}=\left(\frac{2\nu +5}{2\nu }\right)\frac{1}{k_{BRBGF,rigid}}$(14b) , λ is a function of ν, as shown in Eq. 14c(14c) $\lambda =\frac{2\nu }{2\nu +5}$(14c) .

(14a) $\nu =\frac{k_{con,b}}{k_{BRBGF,rigid}}$(14a)

(14b) $\frac{1}{k_{BRBGF}}=\frac{2}{k_{con,b}}+\frac{1}{2k_{con,b}}+\frac{1}{k_{BRBGF,rigid}}=\left(\frac{2\nu +5}{2\nu }\right)\frac{1}{k_{BRBGF,rigid}}$(14b)

(14c) $\lambda =\frac{2\nu }{2\nu +5}$(14c)

Figure 8. BRBGF spring analogy.

The component δ_y,i,col was calculated using Eq. 15a(15a) ${\delta }_{y,i,col}=h_i{\theta }_{y,i,col}=h_i\frac{2{\sum }_{j=1}^{i-1}{\epsilon }_{col,j}{h}_j}{L}=2h_i{\epsilon }_y\frac{{\sum }_{j=1}^{i-1}{\rho }_j{h}_j}{L}=2h_i{\epsilon }_y{\rho }_{avg}\frac{{\sum }_{j=1}^{i-1}{h}_j}{L}$(15a) assuming that the columns exhibit similar deformations in tension and compression (Wijesundara and Rajeev 2012). As it is difficult to determine the timber column strain before choosing the member sizes for BRBs, beams and columns, a strain adjustment factor ρ was used to convert the axial strain of timber column ε_col to the yield strain of yield zone of BRB ε_y, as shown in Eq. 15a(15a) ${\delta }_{y,i,col}=h_i{\theta }_{y,i,col}=h_i\frac{2{\sum }_{j=1}^{i-1}{\epsilon }_{col,j}{h}_j}{L}=2h_i{\epsilon }_y\frac{{\sum }_{j=1}^{i-1}{\rho }_j{h}_j}{L}=2h_i{\epsilon }_y{\rho }_{avg}\frac{{\sum }_{j=1}^{i-1}{h}_j}{L}$(15a) . An average strain adjustment factor ρ_avg along the entire height was empirically assumed as 0.4 to simplify the preliminary design. This value is verified by NLTHAs in Section 4.

(15a) ${\delta }_{y,i,col}=h_i{\theta }_{y,i,col}=h_i\frac{2{\sum }_{j=1}^{i-1}{\epsilon }_{col,j}{h}_j}{L}=2h_i{\epsilon }_y\frac{{\sum }_{j=1}^{i-1}{\rho }_j{h}_j}{L}=2h_i{\epsilon }_y{\rho }_{avg}\frac{{\sum }_{j=1}^{i-1}{h}_j}{L}$(15a)

(15b) ${\rho }_j=f_{sm}{E}_s\frac{{\sum }_{k=j+1}^{k=n}{A}_{c,k}sin\alpha }{A_{col,j}{E}_{GL}}$(15b)

where ε_col,j is the axial strain of glulam columns at the j-th storey below the i-th storey (i.e. j < i); ρ_j is the strain adjustment factor between the glulam columns and BRB at the j-th storey; L (= 8 m) is the span of BRBGF; A_c,k (mm²) is the BRB yield zone area at the k-th storey above the j-th storey (i.e. j < k < n); E_GL (= 10 GPa) is the glulam column elastic modulus; A_col,j is the glulam column cross-section area (mm²) at the j-th storey.

3.2.3. Ductility Factor of the Equivalent SDOF System μ_sys

The ductility factor of the equivalent SDOF system μ_sys is defined by Eq. 16a(16a) ${\mu }_{sys}=\frac{{\sum }_{i=1}^n{V}_i{\mathrm{\Delta }}_i{\mu }_i}{{\sum }_{i=1}^n{V}_i{\mathrm{\Delta }}_i}$(16a) according to Maley, Sullivan, and Corte (2010). Similar to Eq. 4(4) $\mu =\frac{{\mathrm{\Delta }}_d}{{\mathrm{\Delta }}_y}$(4) , the ductility factor for the i-th storey μ_i is defined by Eq. 16b(16b) ${\mu }_i=\frac{({\mathrm{\Delta }}_i-{\mathrm{\Delta }}_{i-1})}{{\delta }_{y,i}}$(16b) . Although the base shear force V_base is initially unknown, Eq. 16a(16a) ${\mu }_{sys}=\frac{{\sum }_{i=1}^n{V}_i{\mathrm{\Delta }}_i{\mu }_i}{{\sum }_{i=1}^n{V}_i{\mathrm{\Delta }}_i}$(16a) contains the shear force in the numerator and denominator. For the initial design, μ_sys can be obtained by assuming a total base shear V_base = 1.0 and recognising that the strength proportions are a design choice. The shear force at the i-th storey V_i was calculated using Eq. 16c(16c) $V_i=\sum _{j=i}^n{F}_j$(16c) . Different from Eq. 7(7) $F_i=\frac{m_i{\mathrm{\Delta }}_i}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}{V}_d$(7) , Sullivan, Priestley, and Calvi (2012) recommended Eq. 16d(16d) $F_j=\left\{\begin{array}{cc}0.9\frac{m_j{\mathrm{\Delta }}_j}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}{V}_{base}& j<n\\ \left(0.1+0.9\frac{m_n{\mathrm{\Delta }}_n}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}\right)V_{base}& j=n\end{array}\right.$(16d) to distribute the base shear force over the height of frame systems, where an additional 10% of the base shear force was allocated at the roof level to consider the higher mode effects.

(16a) ${\mu }_{sys}=\frac{{\sum }_{i=1}^n{V}_i{\mathrm{\Delta }}_i{\mu }_i}{{\sum }_{i=1}^n{V}_i{\mathrm{\Delta }}_i}$(16a)

(16b) ${\mu }_i=\frac{({\mathrm{\Delta }}_i-{\mathrm{\Delta }}_{i-1})}{{\delta }_{y,i}}$(16b)

(16c) $V_i=\sum _{j=i}^n{F}_j$(16c)

where F_j is the distributed load at the j-th storey:

(16d) $F_j=\left\{\begin{array}{cc}0.9\frac{m_j{\mathrm{\Delta }}_j}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}{V}_{base}& j<n\\ \left(0.1+0.9\frac{m_n{\mathrm{\Delta }}_n}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}\right)V_{base}& j=n\end{array}\right.$(16d)

3.3. Equivalent Viscous Damping ξ_eq and Spectral Displacement Reduction Factor η

The relationships among the ductility factor μ, equivalent viscous damping ξ_eq, and spectral displacement reduction factor η of the equivalent SDOF system are required for the DDBD approach, as shown in Fig. 3c and d. ξ_eq has traditionally been calculated by Eq. 17(17) ${\xi }_{eq}={\xi }_{el}+{\xi }_{hyst}$(17) as the summation of elastic viscous damping ξ_el and hysteretic damping ξ_hyst. Traditionally, ξ_hyst has been calculated using Eq. 18a(18a) ${\xi }_{hyst}=\frac{1}{2\pi }\frac{A_{hyst}}{F_m{\mathrm{\Delta }}_m}$(18a) based on an area-based approach proposed by Jacobsen (1960). However, several studies have indicated that the area-based approach can be inaccurate, especially for systems with a high energy dissipation capacity (Dwairi, Kowalsky, and Nau 2007; Grant, Blandon, and Priestley 2005). As such, ξ_hyst was calibrated using the results of NLTHAs and Eq. 18b(18b) ${\xi }_{hyst}=\frac{C\left(\mu -1\right)}{\pi \mu }$(18b) for different types of hysteretic models (Blandon and Priestley 2005; Dwairi, Kowalsky, and Nau 2007; Liu et al. 2015) and structures (Ghaffarzadeh, Jafari, and Talebian 2014; Khan, Kowalsky, and Nau 2016; Landi, Diotallevi, and Tardini 2012; Loss, Zonta, and Piazza 2012; Mazza and Vulcano 2014; Sullivan 2013; Sullivan and O’Reilly 2014; Wijesundara, Nascimbene, and Sullivan 2011; Yahyai and Rezayibana 2015; Yan, Gong, and Zhang 2018) including timber-steel hybrid structures (Bezabeh, Tesfamariam, and Stiemer 2016) and BRB steel frames (Sullivan, Priestley, and Calvi 2012). Even though expressions for timber and BRB systems exist, the hysteretic behaviour of BRBGFs is different from that of other timber-steel hybrid structures. In addition, ξ_hyst for BRB steel frames was derived by using the bi-linear hysteretic shape for BRBs (Sullivan, Priestley, and Calvi 2012), which might not accurately represent the performance of BRBs because the isotropic and kinematic hardening in tension and compression could be significant (Vigh, Zsarnóczay, and Balogh 2017) and the initial slips and semi-rigid connections in BRBGFs might reduce ξ_hyst.

(17) ${\xi }_{eq}={\xi }_{el}+{\xi }_{hyst}$(17)

(18a) ${\xi }_{hyst}=\frac{1}{2\pi }\frac{A_{hyst}}{F_m{\mathrm{\Delta }}_m}$(18a)

(18b) ${\xi }_{hyst}=\frac{C\left(\mu -1\right)}{\pi \mu }$(18b)

where A_hyst is the dissipated energy in an entire hysteretic loop (kJ); F_m and Δ_m are the maximum force (kN) and displacement (mm) for the hysteretic loop, respectively; C is a constant.

In addition to the ductility-dependent equivalent viscous damping expression, an expression is required to scale displacement spectra to different levels of equivalent viscous damping. Eq. 19a(19a) ${\eta }_1=\sqrt{\frac{7}{2+{\xi }_{eq}}}$(19a) and Eq. 19b(19b) ${\eta }_2=\sqrt{\frac{10}{5+{\xi }_{eq}}}$(19b) provide equations for the spectral displacement reduction factor η in the previous Eurocode 8 (1998) and current Eurocode 8 (2005a). This factor is used to scale the displacement spectrum based on ξ_eq and then obtain the required effective period T_e, as shown in Fig. 3d. However, past research (Grant, Blandon, and Priestley 2005; Pennucci, Sullivan, and Calvi 2011) has demonstrated that the ξ_eq-μ relationship is dependent on site seismicity, and this dependency is not explicitly expressed in Eq. 18b(18b) ${\xi }_{hyst}=\frac{C\left(\mu -1\right)}{\pi \mu }$(18b) . Using Eq. 18b(18b) ${\xi }_{hyst}=\frac{C\left(\mu -1\right)}{\pi \mu }$(18b) with Eq. 19a(19a) ${\eta }_1=\sqrt{\frac{7}{2+{\xi }_{eq}}}$(19a) and Eq. 19b(19b) ${\eta }_2=\sqrt{\frac{10}{5+{\xi }_{eq}}}$(19b) may result in inconsistent designs at different sites (Pennucci, Sullivan, and Calvi 2011). Therefore, Pennucci, Sullivan, and Calvi (2011) suggested combining the two steps in Fig. 3c and d into one step and deriving the relationship between η and μ directly because the η-μ relationship did not show an apparent dependency on the site seismicity. A more detailed explanation on using the form of Eq. 20(20) ${\eta }_3=\sqrt{\frac{\pi \mu }{11.04\mu -7.9}}$(20) over Eq. 19 can be found in Dong (2021). Improved accuracy was observed from numerous NLTHAs by Pennucci, Sullivan, and Calvi (2011). In this study, an attempt was made to establish a direct relationship between η and μ for the BRBGF structures.

(19a) ${\eta }_1=\sqrt{\frac{7}{2+{\xi }_{eq}}}$(19a)

(19b) ${\eta }_2=\sqrt{\frac{10}{5+{\xi }_{eq}}}$(19b)

Figure 9 shows the comparison between the numerical modelling results of the BRBGF with dowelled connections in Fig. 2 and the fitted Takeda fat hysteresis model (r_t = 0.05, β_t = 0.6 and D = 0.0) in Fig. 10. Although the Takeda fat model is often used to represent the hysteresis loops of concrete frames (Priestley, Calvi, and Kowalsky 2007), the comparison shows that it can also provide an approximate fit to the hysteresis loops of the BRBGF. The maximum difference of the total energy dissipation was about 10% (291 kJ versus 261 kJ). The η-μ relationship for the Takeda fat model, given by Eq. 20(20) ${\eta }_3=\sqrt{\frac{\pi \mu }{11.04\mu -7.9}}$(20) , was calibrated by Pennucci, Sullivan, and Calvi (2011) with numerous NLTHAs. In this study, we determined whether Eq. 20(20) ${\eta }_3=\sqrt{\frac{\pi \mu }{11.04\mu -7.9}}$(20) could also be applied to BRBGF structures.

(20) ${\eta }_3=\sqrt{\frac{\pi \mu }{11.04\mu -7.9}}$(20)

Figure 9. Comparison between the component-based model and Takeda fat model.

Figure 10. Takeda fat model (Priestley, Calvi, and Kowalsky 2007).

3.4. Calibration of Design Parameters

The stiffness adjustment factor, λ (Eq. 12b)(12b) $\lambda =\frac{k_{BRBGF}}{k_{BRBGF,rigid}}$(12b) must be determined and the applicability of the η-μ relationship (Eq. 20)(20) ${\eta }_3=\sqrt{\frac{\pi \mu }{11.04\mu -7.9}}$(20) for BRBGF structures must also be assessed. Pushover analyses and NLTHAs were conducted to determine λ and assess Eq. 20(20) ${\eta }_3=\sqrt{\frac{\pi \mu }{11.04\mu -7.9}}$(20) , respectively, using numerical models of the one-storey BRBGF as shown in Fig. 5.

3.4.1. Numerical Models

A series of one-storey BRBGF models were built in OpenSees. The models had an 8 m span and 3.6 m height with BRBs installed at an inclination angle α = 42°. All BRBs were made of S235 (British Standard Institution (BSI) 2005b) steel with a material overstrength ϕ_m = 1.2 (Sullivan, Priestley, and Calvi 2012), elastic modulus E_s = 210 GPa and a BRB overstrength factor γ_BRB = 1.5, i.e. ωβ = 1.5, where ω is the BRB strain hardening adjustment factor, and β is the BRB compression strength adjustment factor (American Institute of Steel Construction (AISC) 2016). Three design variables were considered: (1) cross-section of the yield zone in BRBs, (2) stiffness modification factor f_sm, and (3) initial slips due to manufacturing tolerances. Table 1 lists the configurations of the BRBGFs considered. Three cross-sections of BRBs were considered for implementation in the lower, middle and upper storeys of a BRBGF structure (Dong et al. 2020b), which corresponded to a lateral design load of 226 kN, 436 kN and 629 kN, respectively. The BRBGF yield load was calculated by Eq. 9(9) $F_y=2{\varphi }_m{f}_y{A}_{c}cos\alpha$(9) . All connections were designed with a connection relative overstrength factor γ = γ_BRB = 1.5, i.e. the design strength of the connections was 1.5 times the load transferred to the connections when the BRBs yielded (Dong et al. 2021). Although f_sm depends on the BRB geometry (Aguaguiña, Zhou, and Zhou 2019), it typically varies within ± 10% if the same steel grade is used. In this regard, three different values of f_sm were included to consider the influence of BRB stiffness. In addition, three different levels of initial slips caused by manufacturing tolerances were considered for ideally tight, medium and maximum allowable manufacturing tolerances in practice (Dong et al. 2021). Thus, 27 BRBGF configurations were considered in the simulations; each of them was denoted according to the combination of design variables. For example, “75 × 20–1.22-2.5” represents a BRBGF with 75 mm × 20 mm yield zone of BRBs, f_sm = 1.22, and ± 2.5 mm initial slips.

Table 1. Parameters of one-storey BRBGF models.

Cross-section of the yield zone (mm×mm)	fsm	Initial slips (± mm)
45 × 12 65 × 16 75 × 20	1.10 1.22 1.34	0.0 2.5 4.0

3.4.2. Analyses of One-storey BRBGF Models

Following the procedure shown in Fig. 11, numerical simulations were conducted on the 27 BRBGF configurations. Pushover analyses were performed first and the system and component stiffness values were evaluated. Table 2 lists the results for λ which are between 0.70 and 0.74. Therefore, the average λ = 0.72 is used in Eq. 12a(12a) ${\delta }_{y,i,BRB}=h_i{\theta }_{y,i,BRB}=h_i\frac{2{\epsilon }_y}{\lambda sin2\alpha }$(12a) to calculate the BRBGF drift δ_y,i,BRB when the BRBs begin to yield. λ = 0.72 also means that the translationally semi-rigid dowelled connections between the BRBs and glulam frames causes a 28% reduction in the system stiffness compared with the BRBGF with the translationally rigid connection assumption.

Table 2. Results of stiffness adjustment factor λ.

BRB cross-section	45 × 12			65 × 16			75 × 20
fsm	1.10	1.22	1.34	1.10	1.22	1.34	1.10	1.22	1.34
λ	0.74	0.72	0.70	0.74	0.71	0.70	0.74	0.72	0.70

Figure 11. Procedure of parameter verification.

NLTHAs were then conducted using ten ground motion records selected by Maley et al. (2013) from the Pacific Earthquake Engineering Research Centre database (PEER 2005) and GeoNet (Houtte et al. 2017), as listed in Table 3. The near-fault effect was not considered in the records. The average acceleration spectra were scaled to match the design acceleration spectra of soil type D in New Zealand seismic load standard NZS 1170.5 (2004). The average acceleration spectra were scaled to an intensity level of 0.23 g with a return period of 25 years, as shown in Fig. 12, ensuring that T_e of the equivalent SDOF system mostly fell into the T_e > 1.0 s region. Research by Dwairi, Kowalsky, and Nau (2007) and Grant, Blandon, and Priestley (2005) has shown that η is lower for T_e < 1.0 s than that for T_e > 1.0 s. Additionally, T_e for most structures falls into the T_e > 1.0 s region (Priestley, Calvi, and Kowalsky 2007), so scaling into T_e > 1.0 s region will result in more realistic predictions for most structures. The scale factors are denoted as SF1 in Table 3. The design acceleration spectrum S_a(T) in Fig. 12 are transferred to the design displacement spectrum S_d(T) using Eq. 21(21) $S_d\left(T\right)={\left(\frac{T}{2\pi }\right)}^2{S}_a\left(T\right)$(21) to calculate η.

(21) $S_d\left(T\right)={\left(\frac{T}{2\pi }\right)}^2{S}_a\left(T\right)$(21)

Table 3. Ground motion records and scale factors for NLTHA.

No.	Event	Station	RSN	Component	Magnitude	Vs30 (m/s)	SF1	SF2
1	Chi-Chi	TAP042	1430	E	7.62	273	0.71	2.92
2	Landers	Desert Hot Springs	850	090	7.28	345	0.86	3.00
3	Hector	USGS 5295 North Palm Springs Fire Sta #36	1816	270	7.13	345	1.35	5.51
4	Darfield	Westerfield (WSFC)	-a	N00E	7.10	-	1.90	7.23
5	Loma Prieta	CDMG 47179 Salinas-John & Work	800	160	6.93	271	1.08	4.72
6	Kobe	OSAJ	1113	090	6.90	256	0.77	3.43
7	Superstition Hills-02	USGS 5210 Wildlife	729	090	6.54	208	0.41	1.89
8	Imperial Valley-06	Delta	169	352	6.53	275	0.40	1.47
9	Chi-Chi Taiwan-03	CHY055	2477	W	6.20	226	1.54	6.65
10	Chalfant Valley-02	CDMG 54171 Bishop-LADWP South St	549	180	6.19	271	0.83	3.02

aNote: From GeoNet database (https://www.geonet.org.nz/)

Figure 12. Acceleration spectra scaling process.

The NLTHAs were conducted for each one-storey BRBGF model from μ = 1.5 to 6.0 with an increment of 0.5, resulting in 270 equivalent SDOF systems. 27 of them fell outside the matching zone in Fig. 12 and were discarded. The 243 remaining equivalent SDOF systems were used to calculate η by Eq. 22(22) $\eta =\frac{{\mathrm{\Delta }}_d}{S_d\left(T_e\right)}$(22) . The three variables in Table 1 had a minor influence on the η–μ relationship. For example, the results of η from equivalent SDOF systems with different initial slips are plotted with different symbols in Fig. 13 and no significant difference was observed among them. The reason might be that these variables primarily affected the magnitude of the yield drift δ_y,s, while μ included the influence of δ_y,s. Figure 13 shows that the best-fitted curve from the 243 equivalent SDOF systems (denoted as η₄) was similar to η₃ for the Takeda fat model and slightly smaller than η₃, that is, on the conservative side. This agrees with the earlier observation in Fig. 9 that the Takeda fat model tends to slightly underestimate the energy dissipated by BRBGF structures at large displacements. Therefore, η₃ (Eq. 20)(20) ${\eta }_3=\sqrt{\frac{\pi \mu }{11.04\mu -7.9}}$(20) can be used to represent the spectral displacement reduction factor for BRBGF structures conservatively.

(22) $\eta =\frac{{\mathrm{\Delta }}_d}{S_d\left(T_e\right)}$(22)

Figure 13. Verification of η–μ relationship.

4. Case Studies on the DDBD Approach

The DDBD approach with calibrated parameters was used for the ULS seismic design of case study BRBGF buildings with three, six, and nine storeys. Numerical models of these buildings were established in OpenSees and analysed using real earthquake ground motion records. The analysis results were compared with design values to verify the DDBD approach.

4.1. Case Study Buildings

The three BRBGF buildings with three, six and nine storeys were denoted as BRBGF-3, BRBGF-6, and BRBGF-9, respectively. They shared the same floor plans, and BRBGF-6 is shown in Fig. 14 as an example. These buildings are located in Christchurch, New Zealand. Continuous glulam columns for every three storeys were used to reduce splice joints of columns and for transportation convenience, which is also a common practice for timber construction. CLT floors and roofs with concrete topping were used with a rigid diaphragm assumption. The weights of the floors and roof were obtained from the CLT manufacturer (XLam NZ Limited 2016) and all the loading information is listed in Table 4 according to NZS 1170.5 (2004). Appendix I provides a design flowchart of the procedure for the extended DDBD approach mentioned in Section 3.

Table 4. Loading information of case study buildings.

Item	Value	Item	Value
Importance level	2	Return period factor R	1.0
Design working life	50 years	Near-fault factor N	1.0
Annual probability of exceedance	1/500	Dead load on floor	1.8 kPa
Site subsoil class	D	Dead load on roof	1.6 kPa
Hazard factor Z	0.3	Live load on floor	3.0 kPa

Figure 14. BRBGF-6 design example.

4.2. Design Results of the DDBD Approach

The design information for BRBGF-6 is listed in Table 5 as an example. Appendix II shows a detailed step-by-step calculation procedure for parameters in Table 5. The elastic viscous damping ξ_el was assumed to be 2% (Hashemi, Zarnani, and Quenneville 2020). Pennucci, Sullivan, and Calvi (2011) proposed Eq. 23a(23a) ${\eta }_{\xi =2\mathrm{\%}}=\frac{{\mathrm{\Delta }}_d}{S_{d,2\mathrm{\%}}\left(T\right)}={\gamma }_{2\mathrm{\%}}{\eta }_3=0.83\times 0.59=0.49$(23a) that uses γ_2% to adjust η₃ from ξ_el = 5% to ξ_el = 2%, which is denoted as η_ξ=2%. The corresponding design displacement spectrum at ξ_el = 2% (denoted as S_d,2% (T)) was first obtained using η₂ according to Eurocode 8 (2005a), as shown in Eq. 23c(23c) $S_{d,2\mathrm{\%}}\left(T\right)={\eta }_2{S}_d\left(T\right)=\sqrt{\frac{10}{5+{\xi }_{el}}}{S}_d\left(T\right)$(23c) , and then reduced by η_ξ=2% to obtain T_e, as shown in Fig. 15. Subsequently, K_e and V_d were calculated using Eq. 5(5) $K_e={\left(\frac{2\pi }{T_e}\right)}^2{M}_e$(5) and Eq. 6(6) $V_d=K_e{\mathrm{\Delta }}_e$(6) , respectively. Furthermore, the design base shear V_base should include the P–Δ effects using Eq. 24 according to Sullivan, Priestley, and Calvi (2012).

(23a) ${\eta }_{\xi =2\mathrm{\%}}=\frac{{\mathrm{\Delta }}_d}{S_{d,2\mathrm{\%}}\left(T\right)}={\gamma }_{2\mathrm{\%}}{\eta }_3=0.83\times 0.59=0.49$(23a)

(23b) ${\gamma }_{2\mathrm{\%}}={\left(1-0.25\frac{5\mathrm{\%}-{\xi }_{el}}{5\mathrm{\%}}\right)}^{1.5\frac{\mu -1}{\mu }}=0.83$(23b)

(23c) $S_{d,2\mathrm{\%}}\left(T\right)={\eta }_2{S}_d\left(T\right)=\sqrt{\frac{10}{5+{\xi }_{el}}}{S}_d\left(T\right)$(23c)

(24a) $V_{base}=V_d+V_{P-\mathrm{\Delta }}$(24a)

(24c) $V_{P-\mathrm{\Delta }}=C_{P-\mathrm{\Delta }}\frac{{\sum }_{i=1}^n{m}_{i}g{\mathrm{\Delta }}_i}{H_{eff}}$(24c)

(24d) $C_{P-\mathrm{\Delta }}=\left\{\begin{array}{cc}0.0& \frac{M_{eff}g}{K_{eff}{H}_{eff}}<0.05\\ 1.0& \frac{M_{eff}g}{K_{eff}{H}_{eff}}\ge 0.05\end{array}\right.$(24d)

Table 5. Design parameters of BRBGF-6.

Storey	hi	Hi	mi (ton)	Δi (mm)	δy,i,BRB (mm)	δy,i,col (mm)	δs (mm)	δy,i (mm)	μi	Vi	Fi (kN)
1	3.6	3.6	65.6	72.0	11.1	0.0	2.5	13.6	5.3	1.00	31
2	3.6	7.2	65.6	137.7	11.1	1.7	2.5	15.3	4.3	0.94	59
3	3.6	10.8	65.6	197.2	11.1	3.5	2.5	17.0	3.5	0.84	84
4	3.6	14.4	65.6	250.4	11.1	5.2	2.5	18.8	2.8	0.68	107
5	3.6	18	65.6	297.4	11.1	7.0	2.5	20.5	2.3	0.49	127
6	3.6	21.6	39.0	338.1	11.1	8.7	2.5	22.3	1.8	0.25	138

Figure 15. Adjusted displacement spectrum.

where V_P-Δ is the additional base shear owing to the P–Δ effects; C_P-Δ is the P–Δ force adjustment factor; g (= 9.8 m/s²) is the gravitational acceleration.

A similar process was conducted for BRBGF-3 and BRBGF-9, and the design information of the three BRBGF structures is listed in Table 6.

Table 6. DDBD information of case study buildings.

Parameters	BRBGF-3	BRBGF-6	BRBGF-9
Δd (mm)	132.4	243.7	341.1
Me (ton)	151	310	466
He (m)	7.6	14.4	21.4
μsys	4.5	3.9	3.4
ηξ=2%	0.48	0.49	0.50
Te (s)	1.38	2.46	3.4
Vbase (kN)	441	544	622

4.3. Design BRB and Glulam Members

The base shear V_base was distributed over the height according to Eq. 16d(16d) $F_j=\left\{\begin{array}{cc}0.9\frac{m_j{\mathrm{\Delta }}_j}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}{V}_{base}& j<n\\ \left(0.1+0.9\frac{m_n{\mathrm{\Delta }}_n}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}\right)V_{base}& j=n\end{array}\right.$(16d) . The shear force at the i-th storey V_i was resisted by two BRBs together, so the BRBs were designed using Eq. 25a(25a) $A_{c,i}=\frac{N_{i,BRB}}{{\varphi }_m{f}_y}$(25a) and Eq. 25b(25b) $N_{i,BRB}=\frac{V_i}{2cos\alpha }$(25b) according to the force equilibrium. Subsequently, glulam beams and columns were designed considering the BRB overstrength factor γ_BRB (= 1.5). The glulam columns were designed using the maximum axial load. The glulam beams were designed by combining the maximum axial load and corresponding moment caused by the uneven force of the two BRBs at the mid-span (López and Sabelli 2004). The member sizes of the three BRBGF structures are listed in Table 7-9.

(25a) $A_{c,i}=\frac{N_{i,BRB}}{{\varphi }_m{f}_y}$(25a)

(25b) $N_{i,BRB}=\frac{V_i}{2cos\alpha }$(25b)

Table 7. Member size information of BRBGF-3.

Storey	BRB cross section (mm)	Middle span connection (nr × nc)	Beam (mm)	Column (mm)
1	79 × 12	3 × 8	360 × 315	315 × 315
2	64 × 12	4 × 5	360 × 270
3	51 × 8	2 × 5	360 × 270

Table 8. Member size information of BRBGF-6.

Storey	BRB cross section (mm)	Middle span connection (nr × nc)	Beam (mm)	Column (mm)
1	82 × 16	4 × 8	405 × 315	360 × 360
2	77 × 16	5 × 6	405 × 315
3	68 × 16	4 × 7	405 × 270
4	74 × 12	3 × 8	405 × 270	270 × 270
5	53 × 12	3 × 5	315 × 270
6	42 × 8	3 × 3	315 × 270

Table 9. Member size information of BRBGF-9.

Storey	BRB cross section (mm)	Middle span connection (nr × nc)	Beam (mm)	Column (mm)
1	74 × 20	6 × 6	450 × 315	540 × 540
2	72 × 20	5 × 7	450 × 315
3	69 × 20	5 × 7	405 × 315
4	63 × 20	4 × 8	405 × 315	405 × 405
5	70 × 16	4 × 7	405 × 270
6	80 × 12	3 × 8	405 × 270
7	64 × 12	3 × 6	405 × 270	270 × 270
8	69 × 8	3 × 5	315 × 225
9	39 × 8	3 × 3	315 × 225

where A_c,i is the yield zone area of BRB at the i-th storey (mm²); N_i,BRB is the BRB axial load caused by V_i at the i-th storey (kN).

4.4. Verification by NLTHA

4.4.1. Multi-storey BRBGF Models

The DDBD approach was verified using NLTHA in OpenSees. The rigid diaphragm assumption was applied, so the seismic load was assumed to be equally distributed among the four BRBGFs in each direction of the buildings. For simplification, only one bay of the BRBGF structures was modelled and the gravity frames were simulated as a leaning column, as shown in Fig. 16, with the BRBGF-6 model as an example. One-fourth of the seismic weight of the entire building was concentrated and added to the top of each storey as shown in Fig. 16. The total seismic weight m_total was allocated to the BRBGF and the gravity frames according to their tributary areas. The seismic weight of floors and roofs in BRBGFs were calculated using Eq. 26(26) ${\mathrm{S}}_{d,earthquake}=\left(G+{\phi }_{E}Q\right)A_{tributary}$(26) according to NZS 1170.5 (2004) and the load information in Table 4, as shown in Appendix II. The tributary area A_tributary in the X-direction of Fig. 14 was 32 m² which came up with the seismic weight on floor node m_f,f (= 8.7 tons) and roof node m_r,f (= 5.1 tons) for the BRBGFs. The remaining m_total was added to gravity frames, modelled as a leaning column. The seismic weights on floor and roof node were m_f,l (= 48.1 tons) and m_r,l (= 27.8 tons), respectively.

(26) ${\mathrm{S}}_{d,earthquake}=\left(G+{\phi }_{E}Q\right)A_{tributary}$(26)

Figure 16. Six-storey BRBGF OpenSees model.

where S_d,earthquake is the seismic weight (kN); G and Q are the dead load, and live load (kPa), respectively; A_tributary is the corresponding tributary area (m²).

The damping model from Lee (2020a, 2020b, 2021) was used for the NLTHAs. In the current implementation of the Lee damping model in OpenSees, the damping coefficient matrix C_d is assumed to be proportional to the structural tangent stiffness matrix (Mazzoni et al. 2006). The performance and issues of this approach have been discussed by Lee (2022) and some details are briefly provided in Appendix III. The Lee damping model targets a constant elastic damping ratio ξ_el over an assigned frequency range, as shown in Fig. 17, which complies with the basic assumption in the earthquake dynamics of structures better than Rayleigh damping (Chopra and McKenna 2016). The general drawbacks of Rayleigh damping have been discussed in detail (Charney 2008; Chopra and McKenna 2016). Figure 17 shows ξ_el at the first natural frequency (ω_f,s,1) and the second natural frequency (ω_f,s,2) of a six-storey BRBGF structure with 2.5 mm initial slips as well as ξ_el at the first natural frequency (ω_f,ws,1) of a six-storey BRBGF structure without initial slips. The BRBGF structure with initial slips had lower initial stiffness than that without initial slips and thus had a lower first natural frequency ω_f,s,1 compared with ω_f,ws,1. As shown in Fig. 17, the Rayleigh damping model may slightly underestimate ξ_el when the slips of BRBGF structures are overcome and significantly overestimate ξ_el at higher modes compared with the Lee damping model. Therefore, the Lee damping model can avoid unrealistically high damping ratio predictions at higher modes. It also avoids spurious damping forces during the inelastic response (Lee 2020a).

Figure 17. Comparison of damping ratio.

4.4.2. Ground Motion Selection

The same ten ground motions listed in Table 3 were scaled to match the design acceleration and displacement spectra at the ULS with a return period of 500 years. The matching period is from 0.3–3.5 s, as shown in Fig. 18, which is enough to cover the T_e of the buildings studied under this case. The scale factors are listed in Table 3 as SF2.

Figure 18. Acceleration spectrum and displacement spectrum matching.

4.4.3. Displacement and Inter-storey Drift Ratio Response

Figure 19 shows the maximum displacement and inter-storey drift ratio (IDR) under each ground motion EQ-i from the NLTHA of the three BRBGFs compared to the mean and design target values. The mean values of the maximum IDR for BRBGF-6 and BRBGF-9 were 2.0% and 1.7% respectively, which are close to the design drift ratio θ_d (= 2.0%). The rest of the IDR results were smaller than θ_d. The DDBD approach provided reasonably accurate and conservative predictions for the maximum displacement and IDR.

Figure 19. Maximum displacement profile and inter-storey drift ratio response. (a) BRBGF-3. (b) BRBGF-3. (c) BRBGF-6. (d) BRBGF-6. (e) BRBGF-9. (f) BRBGF-9.

Figure 19 also shows that the IDR of the first storeys agreed better with the DDBD target drift than those of the upper storeys. One reason for this might be the actual shear force distribution along the height could be different from Eq. 16d(16d) $F_j=\left\{\begin{array}{cc}0.9\frac{m_j{\mathrm{\Delta }}_j}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}{V}_{base}& j<n\\ \left(0.1+0.9\frac{m_n{\mathrm{\Delta }}_n}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}\right)V_{base}& j=n\end{array}\right.$(16d) because BRB frames are more sensitive to the formation of a soft storey (Erochko et al. 2011). The first storey was designed to have a higher ductility μ₁, as listed in Table 5, so it was expected to enter the yielding stage early. Once the first storey yielded, the stiffness of the storey decreased significantly. Thus, more deformations tended to concentrate on the first storey. Although the soft storey is a potential issue for BRB frames, the design examples illustrate that the DDBD approach can avoid excessive deformation concentration in one specific storey, and the displacement of BRBGFs is well controlled in general. It is also worth mentioning that the similarity between the actual displacement profile and the assumed profile is essential to ensure a good prediction of the DDBD approach (Priestley, Calvi, and Kowalsky 2007). Continuous columns at lower storeys are suggested for BRBGF systems to help avoid excess drift concentration on the first storey, so that the actual displacement profile of BRBGFs is close to the assumed displacement shape. Another reason for the slightly conservative predictions is that the cross-section design of BRBs (Eq. 25(16b) ${\mu }_i=\frac{({\mathrm{\Delta }}_i-{\mathrm{\Delta }}_{i-1})}{{\delta }_{y,i}}$(16b) ) neglects the strain hardening effect after yielding. As shown in Fig. 19, the DDBD approach is more conservative for BRBGF-3 than for BRBGF-6 and BRBGF-9. The reason is that the average strain adjustment factor ρ_avg are 0.17, 0.35 and 0.32 for BRBGF-3, BRBGF-6 and BRBGF-9, respectively. The assumption of ρ_avg = 0.4 for BRBGF-3 overestimated δ_y,i,col and underestimated μ_sys. Because the BRBGF structures are more likely to work as LFRS for buildings with more than three storeys, it is suggested to use ρ_avg = 0.4 for the preliminary design. After the initial designs, further adjustments to ρ_avg along height can be made to optimise the configurations.

4.4.4. Glulam Member Strength Check

The glulam members were all designed by considering the BRB overstrength factor γ_BRB = 1.5 to ensure they remained elastic under ULS loads. The maximum moments and axial forces in the columns and beams from the NLTHA were extracted. Equation 27(16d) $F_j=\left\{\begin{array}{cc}0.9\frac{m_j{\mathrm{\Delta }}_j}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}{V}_{base}& j<n\\ \left(0.1+0.9\frac{m_n{\mathrm{\Delta }}_n}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}\right)V_{base}& j=n\end{array}\right.$(16d) from the New Zealand Timber Structures Standard NZS 3603 (1993) were used to check the beam and column strengths. Table 10 lists the glulam member sizes in the BRBGFs and the maximum combined strength factors (CSFs). The maximum CSF was 0.87 for columns and 0.69 for beams, so all glulam members satisfied the strength requirement.

(27a) $CS{F}_1=\frac{M_{t,mean}}{{\varphi }_{GL}{M}_n}+\frac{N_{t,mean}}{{\varphi }_{GL}{N}_{nt}}\le 1.0$(27a)

(27b) $CS{F}_2=\frac{M_{c,mean}}{{\varphi }_{GL}{M}_n}+\frac{N_{c,mean}}{{\varphi }_{GL}{N}_{nc}}\le 1.0$(27b)

Table 10. Combined strength factors of glulam members.

	BRBGF-3				BRBGF-6				BRBGF-9
Storey	Beam	CSF*	Col**	CSF	Beam	CSF	Col	CSF	Beam	CSF	Col	CSF
1	360 × 315	0.55	315***	0.52	405 × 315	0.69	360	0.87	450 × 315	0.64	540	0.61
2	360 × 270	0.43	315	0.33	405 × 315	0.52	360	0.61	450 × 315	0.58	540	0.56
3	360 × 270	0.30	315	0.03	405 × 270	0.44	360	0.34	405 × 315	0.53	540	0.40
4					405 × 270	0.42	270	0.56	405 × 315	0.45	405	0.58
5					315 × 270	0.49	270	0.44	405 × 270	0.46	405	0.48
6					315 × 270	0.27	270	0.03	405 × 270	0.45	405	0.29
7									405 × 270	0.42	270	0.57
8									315 × 225	0.52	270	0.49
9									315 × 225	0.30	270	0.03

Note: *CSF is the maximum of CSF₁ and CSF₂ in Eq. 27; **Col = column; ***all columns are square.

where N_t,mean and N_c,mean are the mean values of the maximum tension and compression loads under ten ground motions (kN), respectively; M_t,mean and M_c,mean are the corresponding moments (kN∙m) at N_t,mean and N_c,mean, respectively; N_nt and N_nc are the nominal strength in tension and buckling (kN), respectively; M_n is the nominal bending strength (kN∙m); ϕ_GL (= 0.8) is the strength reduction factor for glulam.

It was also observed that the CSFs were primarily smaller than 0.7, and CSFs of the beams were generally smaller than those of the columns at the same storey. This is because the actual overstrength that the BRB reached during the NLTHA was smaller than γ_BRB (= 1.5). The actual axial deformation of BRBs was smaller than the ultimate deformation; therefore, the loads transferred to the glulam members were smaller than the design loads with γ_BRB. The beams had a large span (L = 8 m) and the moment load in the beam was proportional to L², which made the beams more sensitive to the design load level than columns. The member sizes can be optimised by using different γ_BRB along the height based on μ_i as suggested by Lopez and Sabelli (2004), which is beyond the scope of this study.

5. Conclusions

This study extended the DDBD approach to a new timber-steel hybrid multi-storey system that consists of glulam frames and BRBs. The design steps following the DDBD approach are discussed and critical design parameters are derived. The DDBD approach was then used to design three case studies of BRBGF structures with three, six, and nine storeys. NLTHAs were conducted for the three buildings to verify the DDBD approach. The main conclusions drawn are as follows:

1. The inter-storey yield drift of BRBGFs must include the contribution of the BRB deformation, column axial deformation and possible initial slips due to manufacturing tolerances for an accurate prediction of inter-storey yield drifts.

2. In contrast to BRB steel frames, the stiffness of connections between the BRBs and glulam frames in BRBGFs must be considered. A stiffness adjustment factor λ (=0.72) was verified for the yield drift prediction of the BRBGF structures by pushover analyses.

3. The relationship between the spectral displacement reduction factor η and ductility factor μ for the Takeda fat model was also suitable for the BRBGF structures. It provided a slightly conservative prediction based on the results of NLTHA. The strength and stiffness of BRBs and initial slips caused by manufacturing tolerances had a negligible influence on the η-μ relationship.

4. The displacement profile for moment-resisting frames was also suitable for the BRBGF structures and provided relatively conservative predictions. The similarity between the actual displacement profile and the assumed profile is essential to ensure a good prediction of the DDBD approach.

5. The maximum displacement and inter-storey drift response of the three BRBGF structures agreed reasonably well with the design target drift on the slightly conservative side. A soft-storey issue was not observed although the maximum inter-storey drift was on the first storey and were relarively larger than that of other storeys. Therefore, this study demonstrates that the DDBD approach can be a prospective methodology for the seismic design of BRBGF structures. With the DDBD approach, BRBGFs can be designed and used as an effective LFRS for mid-rise buildings up to nine storeys. Further investigations of the DDBD approach, such as the optimisation design and risk-based evaluation, are recommended for future study.

Acknowledgments

The authors acknowledge the Natural Hazards Research Platform, QuakeCore, University of Canterbury, and Shanghai Research Institute of Materials for sponsoring this project. The authors also appreciate the help of the lab technicians Russell McConchie, Alan Thirlwell, Michael Weavers, Peter Coursey and Dave Carney when preparing the manuscript.

Disclosure Statement

No potential conflict of interest was reported by the author(s).

Appendix I

Flowchart of the presented DDBD approach

short-legendFigure 20.

Appendix

The storey height information is listed in Table 5. The mass on each floor and roof (m_f and m_r) for BRBGF-6 can be calculated according to load information in Table 4 and NZS 1170.0 (Australia and New Zealand Standards 2002).

(A1) $m_f=\left(D+0.3Q\right)\times A_f/g=\left(1.8+0.3\times 3.0\right)\times \frac{\left(40\times 24\right)}{4\times 9.8}=65.6tons$(A1)

(A2) $m_r=\left(D+0.0Q\right)\times A_f/g=\left(1.6+0.0\times 3.0\right)\times \frac{\left(40\times 24\right)}{4\times 9.8}=39.0tons$(A2)

where D and Q are the dead load (kN) and live load (kN), respectively; A_f is the tributary area for each BRBGF (m²).

The storey drift Δ_i is calculated based on the design drift ratio θ_d = 2.0% and the displacement profile for MRF as shown in Eq. A3(3) $H_e=\frac{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i{H}_i}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}$(3) –Eq. A10(9) $F_y=2{\varphi }_m{f}_y{A}_{c}cos\alpha$(9) . Δ_i for each floor is listed in Table 5.

(A3) ${\mathrm{\Delta }}_i={\theta }_d{H}_i\frac{\left(4H_6-H_i\right)}{\left(4H_6-H_1\right)}$(A3)

The inter-storey yield drift contributions of the BRB deformation and the column axial deformation at the i-th storey are calculated by Eq. A4 and Eq. A5, respectively. The inter-storey yield drift for each storey along the height of BRBGF-6 is calculated by Eq. A6, and the results are also listed in Table 5. The ductility factor for the i-th storey μ_i is thus calculated by Eq. A7. Eventually, the shear force at the i-th storey V_i in Eq. A8 can be calculated by assuming V_base = 1.0.

(A4) ${\delta }_{y,i,BRB}=h_i\frac{2{\epsilon }_y}{\lambda \sin 2\alpha }=h_i\frac{2{\varphi }_m{f}_y}{\lambda \sin 2\alpha f_{sm}{E}_s}=3600\times 2\times 1.2\times 2350.72\times \sin 84°\times 1.22\times 210000=11.1mm$(A4)

(A5) ${\delta }_{y,i,col}=2h_i{\epsilon }_y{\rho }_{avg}\frac{{\sum }_{j=1}^{i-1}{h}_j}{L}=2\times 3600\times 1.1\times {10}^{-3}\times 0.4\times \frac{{\sum }_{j=1}^{i-1}{h}_j}{8000}$(A5)

(A6) ${\delta }_{y,i}={\delta }_{y,i,BRB}+{\delta }_{y,i,col}+{\delta }_{s,i}$(A6)

(A7) ${\mu }_i=\frac{({\mathrm{\Delta }}_i-{\mathrm{\Delta }}_{i-1})}{{\delta }_{y,i}}$(A7)

(A8a) $V_i=\sum _{j=i}^n{F}_j$(A8a)

(A8b) $F_j=\left\{\begin{array}{cc}0.9\frac{m_j{\mathrm{\Delta }}_j}{{\sum }_{i=1}^6{m}_i{\mathrm{\Delta }}_i}{V}_{base}=0.9\frac{m_j{\mathrm{\Delta }}_j}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}& j<6\\ \left(0.1+0.9\frac{m_j{\mathrm{\Delta }}_j}{{\sum }_{i=1}^6{m}_i{\mathrm{\Delta }}_i}\right)V_{base}=0.1+0.9\frac{m_6{\mathrm{\Delta }}_6}{{\sum }_{i=1}^n{m}_i{\mathrm{\Delta }}_i}& j=6\end{array}\right.$(A8b)

Appendix

The Lee damping model has a damping matrix C_Lee and damping ratio ζ in a form of Eq. A9 and Eq. A10, respectively. The damping matrix C_Lee is symmetric, positive definite or semi-positive definite and possesses classical normal modes. The derivation can be found in (Lee 2020b).

(A9a) $C_{Lee}=\sum _{j=1}^{n_b}\left[M_{Cj}-M_{Cj}{\left(M_{Cj}+K_{Cj}\right)}^{-1}{M}_{Cj}\right]$(A9a)

where

(A9b) $M_{Cj}=4{\zeta }_{pj}{\omega }_{pj}M$(A9b)

(A9c) $K_{Cj}=\frac{4{\zeta }_{pj}}{{\omega }_{pj}}K$(A9c)

and M and K are mass and stiffness matrices, respectively; ζ_pj,ω_pj are parameters that control the values and location of the peak.

(A10a) $\zeta =\sum _{j=1}^{n_b}N\left(\omega ,{\omega }_{pj}\right){\zeta }_{pj}$(A10a)

where N(ω, ω_pj) is a bell-shape based function as shown in Eq. A10 and Fig. 2111.

(A10c) $N\left(\omega ,{\omega }_{pj}\right)=\frac{2\omega {\omega }_{pj}}{{\omega }_{pj}^2+{\omega }^2}$(A10c)

and ω is the circular frequency.

Fig. 2111 shows the relationship between ζ and ω of the Lee damping model with five bell-shape based functions, which provides an almost constant damping ratio within desired frequency range ω_p1–ω_p5. The formulas for ζ_pj,ω_pj can be found in (Lee 2020a) and the configuration of the Lee damping model.

This model has been implemented into a local branch of OpenSees and suanPan (Chang 2018) by the authors for the design verification, and will be merged into the main OpenSees repository soon.

Figure 21. Damping ratio-frequency relationship of the Lee damping model