Seismic assessment of a heavy-timber frame structure with ring-doweled moment-resisting connections

The performance of heavy-timber structures in earthquakes depends strongly on the inelastic behavior of the mechanical connections. Nevertheless, the nonlinear behavior of timber structures is only considered in the design phase indirectly through the use of an R-factor or a q-factor, which reduces the seismic elastic response spectrum. To improve the estimation of this, the seismic performance of a three-story building designed with ring-doweled moment resisting connections is analyzed here. Connections and members were designed to fulfill the seismic detailing requirements present in Eurocode 5 and Eurocode 8 for high ductility class structures. The performance of the structure is evaluated through a probabilistic approach, which accounts for uncertainties in mechanical properties of members and connections. Nonlinear static analyses and multi-record incremental dynamic analyses were performed to characterize the q-factor and develop fragility curves for different damage levels. The results indicate that the detailing requirements of Eurocode 5 and Eurocode 8 are sufficient to achieve the required performance, even though they also indicate that these requirements may be optimized to achieve more cost-effective connections and members. From the obtained fragility curves, it was verified that neglecting modeling uncertainties may lead to overestimation of the collapse capacity.


Introduction
The performance of timber structures under intense earthquake ground shaking depends strongly on the type of failure modes, and in particular, their ductility (Blass et al 1994). The failure of timber elements is usually brittle, whereas the failure of connections between timber elements can be ductile. Design recommendations for timber structures focus mainly on the formation of inelastic deformations on connections by increasing the slenderness of connectors, guaranteeing that failure occurs after yielding of the connectors, and thus enhancing the capacity of joints to withstand large inelastic deformations without rupture. The common analyses methods for the design of multi-story timber buildings in seismic regions do not consider explicitly the inelastic behavior of joints. Since nonlinear dynamic analyses are time consuming and imply high computational efforts, it is common for practicing engineers to design timber framed structures through static analyses, simplified dynamic analyses (or equivalent static), and multi-modal response spectrum analyses (Foliente 1997).
During the design of timber structures, according to Eurocode 8 (CEN 2013), the capacity of a specific structure to dissipate energy is taken into account by considering a behavior factor q, which affects seismic loads by reducing the linear elastic seismic response spectrum. The q-factors considered in EC8 are prescribed considering three classes of ductility: low (DCL), medium (DCM), and high (DCH). Such levels are dependent on the structural typologies and their local and global ductility. The qfactor is a function of the ductility of the structure, which depends on its redundancy and the ductility of elements and connections. Local ductility is associated with the capacity of joints to undergo large deformations without failure. Large displacements can be developed in connections when slender dowels are used and adequate detailing is provided to prevent brittle failure modes such as splitting. Global ductility is related to the ability of the structural system to remain stable, exploiting its redundancy, after yielding of the joints (Jorissen and Fragiacomo 2001). These principles can lead to a reserve in strength and ductility while preventing significant structural damage.
According to Fajfar (1996) the behavior factor q is dependent on the ductility factor, R µ , which represents the dissipative capacity of a structure, and the so-called overstrength, R Ω , defined as the ratio of the actual strength to the design strength. Miranda and Bertero (1994) provided a basic definition of the ductility factor R µ while reporting the influence of specific elements on its value. In Vidic et al (1994) and Miranda and Bertero (1994) it is suggested that the ductility factor depends on the principal elastic period of the structure and the type of soil. On the other hand, several authors, including Mitchell and Paultre (1994), Humar and Rahgozar (1996), and Park (1996), provided discussions regarding the inclusion of overstrength values in the design of structures. These contributions are the basis of the overstrength factors used to compute the behavior factor q prescribed in the present seismic design codes.
Contrary to what happens in the case of steel and reinforced concrete structures, the q-factor values suggested in EC8 for timber structures do not account for the overstrength factor and structural specificities (e.g connection types). On the other hand, only the more common timber building systems are mentioned in EC8, which means that there is an absence of values regarding innovative systems such as the ones consisting of heavy-timber frame structures built with cross-laminated timber (CLT) or moment resisting frames. Nevertheless, the contributions of several authors, including Andreolli et al (2011), Cecotti et al (2013 and Gavric et al (2015), have provided studies on ductility and overstrength properties of innovative structural systems that make use of CLT. For example, in Cecotti et al (2013), the behavior factor q of a seven-story full-scale CLT building was computed. The approach used to calculate the q-factor was proposed in Cecotti and Sandhaas (2010), where the q-factor is given by the ratio of the peak-ground acceleration (PGA) observed in full-scale tests, for a previously defined near-collapse state, and corresponding PGA used in the design. A mixed analytical-experimental procedure was presented in Pozza (2013), by exploiting the results of cyclic tests, which was applied to several building typologies. The procedure proposed in Pozza (2013) combines quasi-static pushover tests of representative walls and analytical models to establish capacity curves. In this study, the q-factor is given by the ratio between the peak-ground acceleration of the ultimate spectra and the yielding spectra (Fajfar 1996).
Despite the remarkable results obtained so far, there is an absence of studies where material uncertainties are considered to compute the ductility properties of timber structures. Timber as a natural material presents a high variability in its mechanical properties (Faber et al 2011), which contributes directly to an increase in the uncertainties associated with the expected strength and peak drift capacities of timber structures. Therefore, it is crucial to develop studies that account for material uncertainties as well as the uncertainties associated with the seismic loading (Foliente 1997).

Numerical modeling of cyclic behavior of timber connections
The numerical models available in the literature for modeling steel connections in timber structures can be classified into three levels, ranging from micro-level to macro-level. In micro-scale models the hysteresis deformation of single connectors is modeled (Foschi 2000). The meso-models simulate the behavior of entire connections built with several fasteners (Rinaldin et al 2013). Macro-scale models can be used to represent the behavior of wood shear walls (Dean et al 1986), diaphragms (Falk and Itani 1989) or complete structures themselves (Foliente 1995). To analyze momentresisting timber frame structures under dynamic loads, meso-level models can be used to capture the behavior of connections and timber elements without disproportionate computational costs. These meso-level models include the use of nonlinear springs that capture the macro response of connections and linear beam-column elements for the structural members.
Laboratory test results are essential to calibrate numerical models and assess their adequacy to capture the degradation modes exhibited by moment-resisting joints. The experimental tests presented in the literature show that the cyclic response of such connections often present a pinched behavior with stiffness and strength degradation (e.g. Folz and Filiatrault (2004), Chui and Li (2005), Polastri et al (2013)). Consequently, the numerical model used to represent the joints must incorporate these features to reproduce the expected progressive weakening of timber structures due to strong ground motions. The degradation modes observed in timber moment-resisting joints, used in heavy timber structures, share many similarities with those observed in reinforced concrete and steel structures allowing the use of models originally developed for other materials, such as the ones in Foliente and Mohammad (1996) Bouc-Wen-Baber Noori (BWBN) model was used by Foliente (1995) to generate a general hysteresis model for wood joints and structural systems, incorporating stiffness and strength degradation with pinching effect. A more general phenomenological model that captures the pinched cyclic strength and stiffness degradation of plywood panels was developed in Ibarra et al (2005). Another model that can be used to represent the pinched hysteretic response of timber connections is a beam-column joint model proposed by Lowes et al (2003), which was implemented in OpenSees as Pinching4 force-deformation model. Although originally developed to simulate reinforced concrete beam-column joints, Pinching4 was also used to perform a seismic reliability analysis of a Timber-Steel-Hybrid-System (Zhang et al 2015) due to its ability to represent a pinched load-deformation response and to enable three modes of cyclic degradation: unloading stiffness degradation, reloading stiffness degradation, and strength degradation.
Despite the lack of a database of experimental results for moment-resisting connections between timber elements, existing experimental results indicate that these connections are able to dissipate energy during cyclic tests without significant strength reduction. One of possible moment-resiting connections in timber frame structures is known as the ring-type doweled connection (Bouchaïr et al 2007). Bouchaïr et al (2007) studied the distribution of loads among fasteners of a ring-type doweled connection. The authors compared theoretical results with experimental monotonic tests in terms of stiffness and strength properties. Ring-doweled joints were also the subject of an experimental campaign (Polastri et al 2013) that comprised monotonic and cyclic reversal tests to study the effectiveness of different joint patterns and their capacity to satisfy the criteria for DCH structures presented in EC8. One of the joint patterns studied in Polastri et al (2013) is used to calibrate the analytical model presented in this paper. Using a bilinear model, Polastri (2010) characterized the q-factor for structures built with ring-doweled moment connections. However, the work in Polastri (2010) did not account for strength and stiffness degradation nor material uncertainties, which could be important to account for uncertainty in the capacity of connections.

Fragility curves for timber structures
A fragility curve is a conditional distribution of the probability of exceeding a specific threshold (e.g. drift, damage, or collapse) as a function of one or more hazard intensity measures (e.g peak-ground acceleration, spectral acceleration at a fundamental period of the structure). In essence, a fragility curve F R is defined as the conditional probability of the structural demand parameter (D) exceeding the structural capacity (C) for a given level of intensity measure (IM ) of the hazard : Fragility curves can be developed using four different approaches: expert judgment, empirical methods, analytical methods, or hybrid methods that combine in some form the first three approaches (Porter 2015). Fragility curves are commonly assumed to follow a lognormal cumulative distribution function (e.g. Rosowsky and Ellingwood (2002), Porter et al (2007), Baker (2015)). The parameters of analytical fragility functions can be determined through different methods. Monte-Carlo simulation is the most widely used method due to its simplicity and robustness. For this, samples of the random variables are generated from their joint probabilistic distributions. The limit state function considered is evaluated for each realization that is treated as a deterministic quantity (Rubinstein and Kroese 2011). In order to achieve a robust fragility function, it is necessary to perform a considerable number of numerical simulations, which, when combined with nonlinear finite element models, may require very high computational efforts. Such drawbacks can be partially overcome by using efficient sampling methods such as Latin Hypercube Sampling (Melchers 1999) or other methods described in Baker (2015).
To assess the response of light-frame timber constructions, several seismic fragility analyses are available in the literature (e.g. Lee and Rosowsky (2006), van de Lindt and Dao (2009), Sutley and van de Lindt (2016)). Seismic fragility curves for shear wall systems, typically found in residential construction were developed in Li and Ellingwood (2007) for three performance levels, including Immediate Occupancy (IO), Life Safety (LS), and Collapse Prevention (CP). The authors considered both demand and capacity related uncertainties, accounting for record-to-record variability due to ground motion accelerogram amplitude and phasing and the effect of openings in the shear walls, known as SAWS. The variability of strength, stiffness, and various hysteretic parameters of a macro-model (Folz and Filiatrault 2001) were considered in Yin and Li (2010). In Yin and Li (2010), it was possible to understand the impact of aleatoric and epistemic uncertainties on the fragility functions. However, the study only considered typical one-story residential buildings with reduced collapse risk.
Even though extensive work has been performed to develop fragility curves for light-frame constructions, limited research has focused on the development of fragility curves for heavy-timber structures, and especially for those designed with ring-doweled moment resisting connections.

Objectives
The main objective of this paper is to evaluate the performance of heavy-timber structures designed with ring-doweled moment resisting beam-column connections. The performance is evaluated based on characterization of the behavior factor q using a large number of nonlinear structural analyses, which account for uncertainty in mechanical properties of members and connections. This paper involves the use of pushover analyses and multi-record incremental dynamic analyses (Vamvatsikos and Cornell 2002), where nonlinear analytical models account for the post-yielding behavior of calibrated moment-rotation connections models. Secondary objectives of this paper are the definition of interstory drift limit states, assessment of the distribution of q-factors, and the development of fragility curves for different damage state levels, including IO, LS, CP, and Global Collapse (GC) accounting for both material and seismic loading uncertainties.

Methodology
A three-story timber moment-resisting frame designed with ring-doweled connections is used as a case study in this paper. The structure had been designed in Callegari (2009) and was here re-checked using a modal response spectrum analysis, following EC8 and Eurocode 5 (CEN 2005a), for a site located in Lisbon, Portugal.
The structure was designed to verify the requirements for a DCH structure according to EC8.
The assessment of the q-factor and the fragility curves involves the use of non- To address the impact of modeling uncertainties, two structural models, in which median parameters and design parameters were assigned, were subjected to the same analysis, for comparison purposes. Henceforth, these structural models will be denominated as "median structure" and "design structure", respectively. It is worth noting that there is no spatial variability of members and connections properties on both median structure and design structure. Moreover, it is important to stress that the post-yielding rotations of the design structure were computed considering the characteristic value of the parameters used to define the envelope curve of the Pinching4 model.
Displacement controlled nonlinear static pushover analyses were performed in order to assess the impact of modeling uncertainties on the q-factor and on the seismic resistance of the timber moment-frame and to define damage state levels. The lateral load distribution, for each structural model is based on the the first vibration mode configuration of that specific structural model realization and a vector of lateral loads {p i } were applied at floor level i, as given by: where {M } is the lumped mass matrix, and {φ i } is the normalized first mode nodal displacements. The lumped masses are assigned to the points coincident with the intersections between beams and columns, at each floor level. The control node is positioned at the roof and a displacement increment of 0.001 m is considered at each analysis step.
The seismic loading uncertainties were combined with the modeling uncertainties performing an IDA study for the set of analytical models generated with LHS method. Henceforth, this set of analytical models will be designated as "structural set". An IDA study is accelerogram and structural model specific (Vamvatsikos and Cornell 2002). It consists in subjecting a structural model to an accelerogram of increasing intensity. Consequently, it is necessary to consider several accelerograms to account for ground motion record-to-record variability. The records used are selected in section 3. The intensity measure (IM) adopted was the "first mode" spectral acceleration (S a (T 1 , 2%)) and the engineering demand parameter (DMP) chosen was the peak interstory drift ratio. The selection of S a (T 1 ) as an IM follows suggestions by Shome (1999) and Baker et al (2008). Other authors have also provided discussions on the topic, including Barbosa (2011), Faggella et al (2013, Mollaioli et al (2013), and Donaire-Ávila et al (2015) and have shown that S a (T 1 ) corresponds to a good predictor of the structural response, both for linear and nonlinear response.
During an IDA, the global collapse of a structural model is related with dynamic instability (Vamvatsikos and Cornell 2002). In this work, it is assumed that an IDA curve reached dynamic instability when its slope is lower than 20% of the initial (elastic) curve.
The total number of IDA curves developed corresponds to the product of the number of models in the structural set (N Sim = 1000) and the number of ground motion records considered (N GM R = 24). Considering approximately 30 intensity levels, the total number of nonlinear dynamic analyses performed was approximately 720,000.
To reduce computational time, a sequential version of OpenSees and a batch-queue system called HTCondor (v7.8.0) was used (Ribeiro et al 2014). From the results obtained from multi-record IDA, the variability of structural responses due to both epistemic and aleatoric uncertainties was thus evaluated. The effect of considering modeling uncertainties was also assessed by comparing the results of the structural set with the results obtained for a structure in which the median parameters were assigned for mechanical properties of members and connections.
Damage state levels were defined using the peak interstory drift ratio as the limit state parameter, which is tied to the behavior of the moment-resisting connections or the building capacity curve, which is defined as the base-shear versus roof displacement. It is assumed that the IO damage state is reached when any connection exceeds the yielding point. In turn, the LS damage state is reached when any connection exceeds the deformation associated with the capping point. Finally, for the CP damage state, it is assumed that it is reached when the descending branch of the building capacity curve passes below 80% of the peak building strength capacity. The GC damage state is exceeded when an IDA curve reaches dynamic instability.
Fragility functions were defined for different damage state levels, based on the multi-record incremental dynamic analysis. These fragility functions were obtained by fitting a lognormal distribution to the values of intensity measure (spectral acceleration) that caused the exceedance of a predefined demand threshold values associated with different damage state levels (IO, LS and CP) considered. In addition, the impact of modeling uncertainties on the GC damage state fragility curves is evaluated. Moreover, incremental dynamic analysis and fragility curves were analyzed for different q-factor levels. A disaggregation process, which consists in aggregating the structural models according to their values of measured q-factors, was followed by the analysis of the common incremental dynamic analysis and fragility curves in each group.
3 Case study: A three-story heavy-timber frame with ring-doweled moment-resisting beam-column connections

Description of the structure
The structure under analysis is a residential three-story timber building designed with GL24h elements. The GL24h correspond to glued laminated timber with materials defined per EN14080, with characteristic bending strength equal to 24 MPa and a mean value of modulus of elasticity of 11.5 GPa. The plane frame structure has two 6-meter long bays in the X-direction and three 3-meter high stories (Figure 1(a)). The columns consist of two elements each with a rectangular cross-section of 160 mm by 600 mm. The beams consist of a single rectangular glulam element with a crosssection of 120 mm by 600 mm. The connections between columns and beams are executed with a ring-doweled joints with two layers of connectors, as shown in Figure   1(b). The first layer of connectors is composed by 10 dowels located 165 mm from the center, and the second layer by 16 dowels located at a radius of 240 mm. All dowels are M4.6 dowels with a diameter of 12 mm. At foundation level, columns are considered as hinged supports. Braced timber frames ensure lateral resistance in the perpendicular out-of-plane direction. To tie the structure together, these elements are fastened with screws to the main beams. Moreover, diagonal steel bars are connected to the columns at floor level, as presented in Figure 1(d).
The floors consist of a low weight solution frequently used in residential buildings in Europe. The top layer is parquet that is placed over a leveling 20 mm thick layer of mortar, followed by 18 mm thick Oriented Strand Boards (OSB) fixed to the GL24h joists with a cross-section of 120 mm by 200 mm. The joists are covered by a wood cover and plasterboard on the underside, and connected with screws to the main beams as shown in Figure 1(e). The model assumes that the floors respond as rigid diaphragms due to their large in-plane stiffness.
Their connection and in-plane stiffness of the floors allows to consider them as rigid diaphragms. For the sake of brevity, further information on the design can be obtained from Callegari (2009).
Glulam Beams (120 The design considered serviceability limit states and ultimate limit states from EC5 and EC8. Nevertheless, the governing loading for the beam column connections used was the seismic loading. The floor loads will be assumed as deterministic, and inertial effects of these loads were considered on the seismic assessment as prescribed in EC8. The combination coefficient (ψ Ei ) was considered for variable loads to account for a reduced participation of their mass, which depends on the quasi-permanent coefficient (ψ 2 ) for variable loads, and is given by: In equation (3), the parameter ϕ is related to the story occupancy. The value assumed is equal to 0.8, which is the value prescribed for stories with correlated occupancies. For residential buildings the quasi-permanent coefficient is equal to 0.3.
Consequently, the mass related with variable loads is computed considering an uniformly distributed area load of 0.48 kN/m 2 . There is no reduction for permanent loads and consequently the area load considered to compute the structural mass is 1.37 kN/m 2 .

Finite element modeling and analysis
The beams and columns were modeled using linear elastic frame elements connected with zero-length springs (Pinching4), as shown in Figure 2(a), to account for the rotational stiffness of moment-resisting joints. Geometric non-linearities are incorporated in the form of P − ∆ effects. Rayleigh damping is assumed and a damping ratio ξ = 0.02 is assigned to the first and second vibration modes. It is worth noting that the Rayleigh viscous damping ratio is assigned to the model to account for energy dissipation that cannot be directly captured in the analytical model (e.g. friction at steel connections and stressing of nonstructural elements). In addition, the energy dissipation due to the hysteretic response of the connections are explicitly accounted for in the Pinching4 model when performing the nonlinear analyses.
The Newton-Raphson method is used to solve the nonlinear algebraic system of equations. The tolerance adopted is 10 −8 on the inner product of the unbalanced load and displacement increments at each iteration (Chopra et al 1995). Newmark integration was used considering γ = 0.5 and β = 0.25, which results in the average acceleration method (Newmark 1959). The time step adopted for the transient analysis is 0.002 seconds.
A preliminary study comprised the evaluation of brittle failure occurrences. This assessment involved performing nonlinear dynamic analysis, where the structure was subjected to different time history records. After each analysis step, the occurrence of brittle failures was checked through a set of user implemented algorithms written in tcl/tk (Welch 1995 Figure 2(c), a symmetric behavior was assumed and consequently the same backbone-curve for both directions was considered. The response envelope obtained experimentally was fitted by defining the points I to IV , shown in Figure 2(b). The calibrated moment-rotation values used to define the backbone curve are presented in Table 1. The points where reloading begins, shown as points V I and IX in Figure 2(b), depend on the ratio between the rotation at which reloading begins and the maximum historic rotation demand (r Disp ), and the ratio between the strength at which reloading begins and the maximum historic strength demand (r F orce ). The moments corresponding to points V and V III are dependent on the ratio of strength developed upon unloading from negative (positive) load to the maximum (minimum) strength (u F orce ). The parameters r Disp , r F orce , and u F orce were assumed equal for both directions given the symmetric response of the connection observed in the testing results. The calibration of these parameters account for the energy dissipated per cycle and also the fitting of the position of points V to IX to match experimental results. Hysteretic damage is simulated through degradation of unloading stiffness, degradation in strength developed in the vicinity of the maximum and minimum rotation demands (strength degradation), and degradation in stiffness (reloading stiffness degradation). The form assumed for each damage rule is the same and represents a more general version of the damage index proposed in (Park and Ang 1985). Each damage index, δ i defined in the Pinching4 model is given by: where d max is given by: and where i refers to the current displacement increment, α i are parameters used to fit the damage rules to the experimental data, E is the hysteretic energy and E monotonic is the energy required to achieve failure under monotonic loading. The values def max and def min are, respectively, the positive and negative deformations that define failure, and d max,i and d min,i are, respectively, the overall maximum and minimum deformation demands achieved until increment i. In Mazzoni et al (2005) and Park and Ang (1985), each parameter of Pinching4 model is presented with more details.
To find the parameters that define pinching and degradation behavior, it was necessary to approximate the numerical model to the experimental results in terms of median strength degradation measured in three completed cycles at the same maximum rotation level. Based on the calibration performed, the parameters α i used to model strength degradation, reloading stiffness and unloading stiffness degradation, are presented in Table 2. Additional parameters obtained from the calibration are r Disp = 0.60, r F orce = 0.50, and u F orce = 0.05. Figure 2(c) shows a good approximation of the Pinching4 model adopted to the experimental results, both in terms of hysteretic response and strength degradation. Moreover, the approximation between experimental and analytical results can also be evaluated in terms of dissipated energy, as shown in Figure 2(d), where the maximum relative error determined was 18%, although it can be seen that the difference is mainly due to the energy dissipated in the initial cycles with is not as well captured by the model or could be associated with damping in the experimental setup.

Uncertainties in timber members
The inherent uncertainties of timber, as a material, are considered in this work by assigning different mechanical properties to each element. As presented in Table 3 Table 4 presents the intra-element correlation coefficients considered in this work.  It is considered that the properties of different elements are also correlated. A high inter-element correlation coefficient of 0.8, shown in Table 5, is assumed to simulate the assumption that laminations of different elements origin from the same consignment of sawn wood.   (2013), the elastic rotational stiffness K el can be estimated as: where n sp is the number of shear planes, r i is the distance of each fastener from the center, and k ser is the slip modulus per shear plane (CEN 2005a) given by: where ρ con is the wood density, in kg/m 3 , and d is the diameter of the dowels, in millimeters. It is worth noting that due to the variability of connected elements it is necessary to consider different densities for different members. Thus, ρ con used in the equation (7) is the geometric mean of the density of adjacent members given by ρ con = √ ρ den,1 · ρ den,2 , as proposed in EC5 for connections involving members with different densities. In this work, ρ den,1 and ρ den,2 are the densities of outer (columns) and inner (beams) elements respectively. The rotational stiffness of the beam column joint is then defined by the following equation: As reported in Polastri et al (2013), equation (6) overestimates the rotational stiffness. Considering the envelope experimental backbone curve, the value obtained through the experimental tests (K el,test ) is equal to 6936.6 kN.m/rad. Nevertheless, the rotational stiffness obtained through equation (8) is equal to 12572.0 kN.m/rad, with a measured wood density of 467 kg/m 3 . Thus, to consider this feature the rotational stiffness used to perform the nonlinear analysis is affected by a bias coefficient 55. According to Porteous and Kermani (2013), the yielding moment of a ring-doweled joint can be determined from equilibrium as where F R,i is the yielding force of fastener i, which is given by: where the parameter A is a constant that accounts for the nonlinear load-deformation response of isolated dowels and r max is the distance between the center and the farthest layer. For the moment-resisting connection studied in this work, A is assumed equal to 0.5 (Polastri et al 2013). The determination of dowel resistance force F V,R per shear plane is based on the Johansen's yielding theory (Johansen 1949). Considering moment resisting connections used in the case study, the resistance force F V,R of a dowel per shear plane is computed with the following equation: where t i is the thickness of the members, f h,α,i is the embedment strength of the members in the direction of the load applied (α) on the fastener; β is the ratio between the embedment strength of the beam and the embedment strength of the column (β = f h,α,2 /f h,α,1 ), and d is the diameter of dowels. Thus, f h,α,1 and f h,α,2 are determined considering the influence of the angle α relative to the grain (α = 0 corresponds to the direction parallel to grain). The embedment strength f h,α is computed through the following equation: where k 90 = (1.35 + 0.015d) for softwood (CEN 2005a) and f h,0 is the embedment strength in the parallel to grain direction given by: Failure modes given in equation (11 a) where f u is the ultimate yield capacity in tension. In this work, the ultimate capacity in tension of dowels f u is also assumed as a random variable varying from connection to connection. This variable follows a lognormal distribution with a median value of 427.2 MPa and a coefficient of variation equal to 0.04 (PMC 2011) corresponding to a characteristic value of 400 MPa.
A bias coefficient λ My has to be considered to affect the values of yielding moment obtained through equation (9). Such as in the case of rotational stiffness, the yielding moment is also over-estimated by the formulae proposed in EC5. Consequently, the coefficient λ My is the ratio between the yielding moment of the experimental backbone curve (M I = 83.93 kN.m) and the predicted yielding moment M test y,EC5 (= 101.04 kN.m) obtained through equation (9) using the mechanical properties of the experimental specimens (f u = 580 MPa). Thus, the deterministic value of λ My is considered in the Latin Hypercube samples is given by: Consequently, the yielding moment used for Latin Hypercube samples is defined by the following equation: From the set of structural model samples generated using the parameters presented in Tables 3 to 5, the yielding moment (equation (9)) and elastic rotational stiffness (equation (8)) were computed for each joint, and each structure in 9000 size samples for each property with parameters shown in Table 6. Due to the current requirements of EC8 for high ductility class structure, it is worth noting that the capacity of the ductile failure mode shown in equation (11 d) was the one that governed the strength capacity of the dowels, even when variability on timber properties and dowels ultimate strength was assumed. This is a consequence of the high slenderness of dowels (t ≥ 10d) and the mild steel quality assumed during design.
The parameters defining the post-yield behavior of the hysteretic models defined for the connections are also considered as random variables. The envelope curve of the cyclic test presented in Figure 2( ). The rotations necessary to define the backbone of moment-resisting joints are also assumed as random variables. Their values are defined considering the yielding rotation θ y and the parameters X II , X III , and X IV , as presented in the following equations: Lognormal distributions were assumed for X II , X III and X IV , as these parameters must be positive. The mean values for the parameters are illustrated in Figure   3, which was obtained from calibration of the backbone curve of the experimental results shown in Figure 2(c). Since there are no studies available in the literature regarding the variability of post-yielding rotations of ring-doweled connections, and, there are more factors (e.g. cracks on the vicinity of dowels) that may induce deviations on the idealized backbone, a high coefficient of variation was assumed for the variables X II , X III and X IV . The parameters used to represent post-yielding rotation properties are presented in Table 7 and illustrated in Figure 3.
Additional parameters used to define the hysteretic behavior of connections are assumed to be deterministic.

Distribution of the periods of vibration of the structural set
The random variables considered in this work affect mass and stiffness of the structure which in turn induce variability on periods of vibration of each generated structural model realization. It is worth noting that the fundamental period T 1 plays an important role on the development of multi-record IDA, since the time-history records considered are scaled to the same first mode spectral acceleration S a (T 1 ), and therefore its expected value needed to be determined in this section before performing the ground motion selection and before defining the reference intensity measure. The first natural period for the 1000 structural models generated through LHS has a median value equal to 0.97 s while the median of the second natural period is 0.14 s.
All the realizations of the structural set have in common the fact that 95% modal mass is reached with only two lateral modes of vibration. The probability distribution of computed natural periods can be approximated by lognormal distributions (p-value equal to 0.83 and 0.99 respectively, for a significance level of 5%): (T 1 ∼ LN (µ LN = −0.032, ξ LN = 0.057) and T 2 ∼ LN (µ LN = −1.991, ξ LN = 0.052)).

Ground motion selection
A set of 24 ground motion records was selected and scaled to the 2% damped linear elastic acceleration response spectra considered in EC8 for a site in Lisbon, Portugal.
According to the national document of application of EC8, it is necessary to consider two types of seismic action for any site in Portugal, including the Type 1 far-field ground motions which correspond to magnitudes higher than 7.0, and Type 2 nearfield ground motions characterized by magnitudes lower than 7.0. The set of ground motions were extracted from the Pacific Earthquake Engineering Next-generation Attenuation (PEER NGA) database (PEER 2002). The main characteristics of the timehistory records selected are presented in Table 8. The respective response spectra are shown in Figure 4 where the elastic response target spectra of EC8 is also presented.
In this figure, it can be seen that the response spectra of the selected records match the target spectra in the period range of interest, which was selected between 0.2T 1 and 3.0T 1 , where T 1 is the median fundamental period of vibration of the structural set. The values presented in Table 8 to characterize the time-history records used are the peak-ground acceleration (PGA), the "first-mode" spectral acceleration (S a (T 1 )), obtained for the median fundamental period of the structural set, and the scale factors (SF i ) used to scale the records to the target response spectra.

Performance assessment
In this section, the seismic capacity of the structural set of 1000 structural models obtained with LHS method, is evaluated through nonlinear static analysis and incremental dynamic analysis.

Estimates of equivalent reference system parameters
From the results of the capacity curves shown in Figure 5(a), the parameters shown in Figure 5(b) that define the equivalent bilinear inelastic were obtained, including the peak base shear F * max , the reference yield displacement d * y , and associated yield force F * y , and the ultimate displacement d * u . The initial stiffness is determined considering the secant stiffness to the point of the first yield of one of the moment resisting connections. The ultimate displacement d * u is defined as that corresponding to a decrease of 20% from F * max ("near -collapse" state). The yielding force F * y results from the application of the Energy Equivalent Elastic Plastic (EEEP) method (Foliente 1996), which consists in balancing the total energies below the obtained capacity curve and the equivalent bilinear model. The distribution of the parameters are listed in Table 9.

Estimation of the q-factors
According to Fajfar (1999) the q-factor is given by: where R µ is the ductility factor and R Ω is the overstrength factor. In this paper, R µ is computed with the formulae proposed by Vidic et al (1994): where the coefficient µ is the ratio between ultimate and yielding displacements of the equivalent inelastic model (i.e. µ = d * u /d * y ), T 1 is the fundamental period of a structure, and T C is the transition period that is equal to 0.8 s for the site considered in this study. All the fundamental periods computed, for the structural set, are higher than T C , which means that R µ = µ.
In turn, the overstrength factor R Ω considered in this work is equal to the ratio between the yielding force F * y , of each structural sample, and the base shear force F * 1,d . This value was taken from the design structure capacity curve, presented in Figure 5(a), when any connection reaches the yielding moment (F * 1,d = 59.8 kN).
The quantities used to evaluate the q-factor can be computed using the parameters listed in Table 9. The table shows that the median value of the q-factor obtained from the pushover analysis with modeling uncertainties is significantly higher than the q-factor recommended in EC8 for structural design which are q 0.50 = 7.1 and q EC8 = 4.0, respectively. Moreover, it is worth noting that the median value of the ductility factor R µ obtained is 4.0, which is similar to the q-factor used in the design of this structure. On the other hand, a linear correlation coefficient of 0.46 and 0.74 are obtained between the q-factor and R Ω and between the q-factor and R µ , respectively. A lognormal probability density function was fitted to the values of qfactor obtained, which is shown in Figure 6. It can be seen that the median value of the q-factor is approximately 7.1. Table 9 for the R µ and q-factor indicate that the q-factor considered in EC8 and the detailing requirements defined in EC8 and EC5 are adequate for the design of this moment resisting frame structure, even though they also indicate that the detailing requirements may be optimized to achieve more cost-effective connections and members in future studies.  16%[ −[16%,50%[ −[50%,84%[ −[84%,100%] Behaviour factor (q) Probability Density Function (PDF) Fig. 6 q-factor levels and fitted PDF

Estimation of interstory drift ratio limit states
The estimation of a relationship between the interstory drift ratios and the joint rotations is performed in this section. These values are obtained here from the capacity curves, where each damage state is associated with a certain structural damage. It is assumed that the Immediate Occupancy damage state (IO) is reached when an interstory drift ratio associated with the first yielding of any connection in the building.
It is considered that Life Safety damage state (LS) is reached when any momentresisting connection reaches its capping rotation θ III , shown in Figure 3. For the Collapse Prevention (CP) damage state, the threshold is given from the highest interstory drift ratio obtained when the structure reaches the "near -collapse" state (20% decrease from maximum capacity). Table 10 lists the lognormal parameters obtained for each limit state interstory drift ratio. The results in Table 10 show that the mean threshold values of interstory drift ratio are 1.2%, 4.9%, and 7.9% for the IO, LS, and CP damage states. The values of the coefficient of variation vary from 9% to 14.4%. The Kolmogorov-Smirnov (KS) goodness-of-fit test confirms that the lognormal distribution is an acceptable distribution for the limit state thresholds considered with a confidence level of 95%. The goodness-of-fit of the distributions is also shown in the lognormal probability plots presented in Figure 7.

Incremental Dynamic Analysis
In order to evaluate the effect of record-to-record variability, the single-record curves IDA for both design structure (T 1 = 1.13s) and median structure (T 1 = 0.97s) are shown in Figure 8(a) and 8(b), respectively. The design structure shows a lower capacity both in terms of spectral acceleration and peak interstory drift ratio at global collapse. Despite the differences regarding strength and stiffness between the two structures, the ability of joints to deform in the nonlinear range contributes considerably to these results. The parameters X II , X III and X IV , used to compute the backbone curve, differ considerably since the 5th percentile value was used for the design structure, instead of the median value used for the median structure. In Figure 9(a), the results for all multi-record IDA curves (N IDA = 24000) are summarized by the mean and median IDA curves for the structural set. In addition, the mean and median IDA curves of the median structure (N IDA = 24) are also shown in this figure. Figure 9(b) shows the coefficient-of-variation of the estimated capacity (S a (T 1 )) versus peak interstory drift ratio.
From the results shown it possible to see that, structural set and median structure present similar mean and median curves for peak interstory drift ratios lower than 5%.
The coefficient of variation of spectral acceleration is higher for the structural set, for peak interstory drift ratios lower than 5%, as shown in Figure 9(b). These results were expected as a consequence of including structural uncertainties. For peak interstory drift ratios greater than 5%, the median structure presents higher mean and median capacity, but similar coefficient of variation, when compared to the structural set.
Similar observations have been reported in the literature for RC structures (Liel et al 2009  The seismic capacity of heavy-timber frame structures depends strongly on the hysteretic energy dissipation of connections, which is related to the deformation capacity of the moment-resisting joints. In the following disaggregation of the results, the q-factor is subdivided into four intervals, as shown in Figure 6. The first interval comprises structures with q-factor lower than its 16 th percentile (q 0.16 = 6.5), while the second groups structures with q-factor between, the 16 th percentile q 0.16 and the 50 th percentile q 0.50 , respectively (q ∈ [6.5, 7.1[). The third interval is equal to q ∈ [q 0.50 , q 0.84 [ and the last one comprises structures with q-factor higher than its 84 th percentile (q 0.84 = 7.7).

Fragility Analysis
The fragility functions computed here result from fitting a lognormal distribution to the spectral accelerations extracted from the IDA curves for each limit state. In this section, the effect produced by modeling uncertainties and q-factor variability on the different fragility functions are presented. Table 11 summarizes the results for the fragility curves, determined for IO, LS, and CP damage states. For each damage state level, three different values were chosen to define exceedance of a certain capacity level, namely 16 th , 50 th and, 84 th percentile of each limit state. The influence of limit state threshold on the fragility functions can therefore be assessed for the median structure and the structural set. It can be seen that the expected values are similar for both structural set and median structure when IO damage state is considered. Nevertheless, the coefficient of variation is higher for the structural set. For example, in Table 11, it can be seen that, when the median value of interstory drift ratio is considered as threshold (θ max = 0.0118), the difference on expected values is negligible whereas the coefficient of variation of the structural set is 18% higher than the one obtained for the median structure. The expected values of the structural set are 2% to 4% higher when compared with the median structure for LS damage state. In terms of coefficient of variation, the values obtained for the structural set are higher, raging from 11% to 13%. For the CP damage state, the differences between the structural set and the median structure are negligible in terms of coefficient of variation. Never-theless, the expected value obtained for the median structure is 3% to 6% higher than the expected value of the structural set. From the IO and LS fragility curves, graphically represented in Figure 11(a) and Figure 11(b), respectively, it is observed that the limit state interstory drift ratios admitted as thresholds influence both the expected value and the coefficient of variation. Consequently, these values may be considered as random variables in further studies.    Table 12, the expected value of GC damage state peak interstory drift ratio, for the structural set, is lower than the expected value used for CP limit state value. Consequently, a considerable number of structural models reached instability for the interstory drift ratios considered for CP. Moreover, the coefficient of variation of global collapse peak interstory drift is 37% higher than the one obtained for CP damage state.
The impact of considering modeling uncertainties, on a GC fragility curve, is reflected in the expected values presented in Table 12, for both structural set and median structure. As shown in Figure 12(b), a similar dispersion is observed, but neglecting modeling uncertainties may lead to an overestimation of the collapse capacity.
A higher variability is observed for the structural set whereas the median structure presents a greater expected value. The cumulative distribution functions of peak interstory drift ratios are presented for GC in Figure 12(a).  The disaggregation of the multi-record IDA results by different levels of q-factor, allow to present fragility curves for different q-factor levels. Table 13 shows results for both LS and CP damage states, considering median values of interstory drift ratios as limit-state thresholds. The disaggregation of fragility curves in different q-factor levels shows that higher q-factors correspond to higher expected values and higher coefficients of variation of the structural capacity. These results are shown graphically in Figure 13. The capacity/demand ratio for timber members and the beam-column joints were checked for a sample of the ground motions with the 5% damped linear response spectral accelerations at the fundamental period of vibration of the structure ranging from 1.4g to 4.5g. It was observed that neither the members nor the beamcolumn joints reached their ultimate deformation capacity.  [0,16%[ [16%,50%[ [50%,84%[ [84%,100%]

Conclusions
The performance of a heavy-timber structure designed with ring-doweled moment resisting connections was evaluated in this work. It is important to note that the results obtained here are specific to one direction of a three-story building. The ringdoweled joints used to connect beams and columns had already been experimentally studied under cyclic testing in Polastri et al (2013). The results had shown that the connection could fulfill the requirements of EC8 for high ductility class structures.
To evaluate the effectiveness of this connection on a prototype structure, a three-story building designed following EC5 and EC8 codes was subjected to a comprehensive seismic performance assessment which included numerical nonlinear static and nonlinear dynamic analysis. OpenSees was used for the numerical analysis, in which the Pinching4 constitutive model was used to capture the moment-rotation behavior of ring-doweled moment resisting connections, which was calibrated based on testing data available in the literature. The inherent variability of the timber structural members was included in the analysis as modeling uncertainties, which influenced the structural capacity and the notable points used to characterize the hysteretic response of the connections. Using a set of 1000 structural models generated with the Latin Hypercube Sampling method, a probabilistic assessment was performed including spatial variability of strength and stiffness of timber elements and connections properties.
Nonlinear static analyses were first performed to evaluate the variability of limitstate interstory drift ratios and q-factors. The consideration of modeling uncertainties allowed the estimation of a probabilistic distribution of q-factors for the structural set.
The main observations from this part of the work were: -Parameters for the peak value of interstory drift ratio θ max associated with different damage states were obtained. Median values obtained were 1.2%, 4.9%, and 7.9% for the IO, LS, and CP damage states, respectively. In addition, the coefficient of variation of θ max ranged from 9% to 15%, which are relatively low values, specially since modeling uncertainties were explicitly considered.
-The median value of the q-factor obtained from analyses was considerably higher (q = 7.1) than the value prescribed in EC8 (q EC8 = 4.0). Moreover, the median value of the ductility factor R µ obtained was 4.0. These results indicate that the q-factor values considered in EC8 and the detailing requirements defined in EC8 and EC5 are adequate for design of this type of structure.
In a second part of the work, a multi-record incremental dynamic analysis (IDA) was performed assuming the first mode spectral acceleration as the intensity measure and peak interstory drift ratio as the damage measure. Main observations from the IDA curves were: -The results showed that modeling uncertainties have a slight influence on the expected values of the IDA curves for peak interstory drifts ratios lower than 5%.
Nonetheless, when modeling uncertainties are taken into account the coefficient of variation increases up to 43% more.
-From the disaggregation of the IDA curves according to four different q-factor levels, it was observed that structural models with higher q-factors are more likely to resist ground shaking with higher intensities. These results can be partially explained due to the fact that a positive linear correlation of 0.46 was observed between the q-factor and R Ω .
-No brittle failures were observed in the dynamic analyses conducted to compute the IDA curves, indicating that the sizing requirements in EC8 are adequate, but potentially too conservative, thus leaving room for improvements of the slenderness of the members and dowels. Several seismic fragility curves for different performance levels, including Immediate Occupancy, Life Safety, Collapse Prevention, and Global Collapse, were determined using multi-record incremental dynamic analysis. These were based on interstory drifts that were defined considering the limit-states of the moment-resisting joints and the building capacity curves, specifically for the global collapse estimates.
In this case, the main observations were: -Modeling uncertainties did not affect the median values of the fragility curves associated with IO and LS. However, the coefficient-of-variations increased by 18% and 13% when the modeling uncertainties were considered for the IO and LS damage states, respectively.
-When modeling uncertainties are neglected an overestimation of the capacity is obtained, both in terms of spectral acceleration and peak interstory drift for the CP and GC damage states, by approximately 3% and 5%, respectively. However, the coefficient-of-variation did not change for the CP and GC damage states, when modeling uncertainties were considered.
As described above, there is room to perform further experimental tests to evaluate how a reduction in connected elements thicknesses and slenderness of dowels would impact the ductility and strength degradation of ring-doweled joints. Such tests, along with the methodology proposed in this work, could contribute to propose new design values and detailing requirements to moment-resisting joints, in future works. Moreover, the variability on the responses observed during this study also indicates that more experimental campaigns need to be performed in the future to build a database of moment-resisting timber connections. This would also allow for characterization of the uncertainty of the expected model parameters used in design and their correlation with observed joint performance.

Acknowledgments
This work had financial support of the Portuguese Science Foundation (FCT), through PhD grant PD/BD/113679/2015 included in InfraRisk-PhD program. The results regarding the numerical model for the connections were developed during a Short-term Scientific Mission (STSM) within scope of COST Action FP1101: Assessment, Reinforcement and Monitoring of Timber Structures. Thus, the insightful discussions with Professor Roberto Tomasi and Doctor Andrea Polastri are gratefully acknowledged. The first author would also like to acknowledge the support of Oregon State University during the period in which he was a visiting Ph.D. student, at this institution.

A Notation
The following symbols are used in this paper: T 1 = First natural period; and T 2 = Second natural period;