Abstract The structural strength requirements for timber buildings have been significantly tightened in the second generation of Eurocodes (EN 1990:2023, EN 1991-1-7), which poses a particular challenge for solid timber frames with a beam-and-column structure, where the transfer of tensile forces via dowel connections is inherently limited. Existing multiscale frameworks for timber post-and-beam robustness lack operational detail at each scale, and no validated workflow currently bridges joint-level continuum damage mechanics and frame-level progressive failure analysis in compliance with the second-generation Eurocodes. This paper addresses this gap by proposing an effective two-scale finite element method (FEM) modelling framework for assessing the strength of such frames during column removal. Existing multiscale models describing the strength of timber structures with beam-and-column systems lack the operational details necessary to integrate failure mechanics at the joint level and progressive failure modelling at the frame level within a single, validated workflow. In this paper, this gap is addressed through three specific contributions: a physically modified quadratic Hashin-type failure criterion for timber, which eliminates the non-physical increase in shear strength under combined stress states perpendicular to the grain; a two-scale structure based on the finite element method (FEM), in which the results of continuous damage mechanics at the joint level directly parameterise non-linear joint elements with six degrees of freedom at the frame level, taking into account coupled directional wear and erosion of the elements; and quantitative validation of both scales against experimental data and the conversion factors for characteristic values of the second generation of Eurocode 5 (prEN 1995-1-1:2023). At the connection level, the simulated strength and stiffness values agree with the experiments to within an error of no more than 5%. At the frame level, the model correctly reproduces the non-linear ‘load–displacement’ relationship, the sequence of joint failure, and the axial forces in the chain line for vertical displacements up to 390 mm, which corresponds to experimental observations. 1. Introduction Although standard EN 1991-1-7 requires a dynamic non-linear analysis to be carried out for buildings in category CC3, the proposed concept provides a non-linear structural model—including constitutive properties, progressive failure of connections and chain deformation effects at large displacements—which serves as the basis for such an analysis. The next immediate development step is to extend this to a full dynamic implementation, taking into account physical inertia and damping, based on the quasi-static validation presented here. The main objective of this article is to address the gap caused by the lack of a fully functional, experimentally validated and Eurocode-compliant two-scale modelling strategy for assessing the stability of solid timber frames with columns and beams under column removal conditions—a gap that has been clearly identified but not addressed by previous comprehensive approaches. To this end, three specific contributions to the modelling have been made: (i) a physically modified quadratic Hashin-type failure criterion for timber, in which the non-physical term of the stress product σ 22 σ 33 is removed from the yield surface perpendicular to the fibres to eliminate the artificial increase in shear strength under combined stress states; (ii) a continuum fracture mechanics model at the joint level, which directly parameterises the non-linear joint elements with six degrees of freedom at the frame level, taking into account coupled directional damage and erosion of the elements, forming a complete and transferable workflow linking different scales; and (iii) quantitative validation of both scales against published experimental data and against the conversion factors for characteristic values of the second-generation Eurocode 5 (prEN 1995-1-1:2023), ensuring direct applicability in the regulatory design of structures in accordance with EN 1990:2023 and EN 1991-1-7. Compared with existing approaches, the proposed methodology makes the following unique contributions: (i) it provides a fully operational two-scale FEM scheme in which the fracture mechanics model at the joint level directly influences the global elements of the framework, thereby addressing the shortcoming of previous global methodologies associated with the lack of details ensuring the connection between scales; (ii) it introduces a modified Hashin-type failure criterion with a physically justified correction to the yield surface perpendicular to the fibre direction; (iii) it implements non-linear joint elements with six degrees of freedom, accounting for coupled degradation and erosion of the elements; and (iv) it is the first framework to have been directly verified against the characteristic values of stiffness and strength of the second generation of Eurocode 5, making it immediately applicable in the regulatory design of structures. Existing approaches to modelling the strength of timber structures fall into three categories, each of which has certain limitations: Comprehensive multiscale models (Voulpiotis et al. [ 7]): these justify the conceptual necessity of a multiscale approach, but clearly do not provide detailed calculation procedures for each scale, thereby depriving structural engineers of a working tool. Machine learning and stochastic approaches (Alshboul et al. [ 5], Cao et al. [ 8]): provide probabilistic results at the system level, but sacrifice the deterministic ‘stress–strain’ traceability required for regulatory design. None of the existing approaches simultaneously satisfies all three requirements of the modern regulatory framework: accuracy at the joint level, modelling of progressive failure at the frame level, and compliance with the definitions of characteristic values in the second-generation Eurocode. For structural joints modelling, the following aspects are covered: ▪ Joint’s non-linear FEM model with continuum damage mechanics for modelling up to failure load (no post-failure models). ▪ Discussion of joint load–displacement curves according to EC 5 [ 13] and EC 8 [ 14]. ▪ Hybrid model, where the initial part up to peak load is modelled with non-linear FEM and post-peak behaviour according to EC 8. For global FEM modelling, the following aspects are covered: Multifiber beam non-linear finite elements with plasticity and damage behaviour, including the element erosion technique. Nonlinear joint structural elements, where a joint-scale model defines their behaviour. Robust and straightforward moments, axial, and shear strength interactions for structural joints. The novelty of this work does not lie in the application of FEM to timber structures, which is well established [ 9, 10, 11, 12], but in four specific contributions that collectively address an unresolved gap in the design of timber post-and-beam buildings for structural robustness. First, a fully operational two-scale FEM workflow is developed and validated, in which joint-scale continuum damage mechanics results directly parameterise frame-scale non-linear joint elements—providing the detailed scale-bridging procedure that prior holistic frameworks identified as necessary but did not deliver. Second, a physically motivated modification to the MAT 143-type Hashin failure criterion is introduced, removing the non-physical stress-product term σ 22 σ 33 from the perpendicular-to-fibre yield surface to eliminate an artificial shear-strength enhancement under combined stress states. Third, a six-DOF non-linear joint element with coupled directional degradation and element erosion is proposed, enabling simulation of progressive and sequential joint failure at the frame scale. Fourth, the framework is the first to be explicitly validated against the characteristic value conversion factors of the second-generation Eurocode 5, establishing direct applicability in normative design practice under EN 1990:2023 [ 2] and EN 1991-1-7 [ 3]. Despite these contributions, none of the existing approaches simultaneously satisfies the three requirements imposed by the current regulatory environment—EN 1990:2023, EN 1991-1-7, and EU 305/2011 [ 1]—for robustness assessment of CC3 timber post-and-beam buildings. Their specific limitations are as follows. Holistic multiscale frameworks establish the conceptual necessity of treating tall timber buildings as complete systems and propose a multiscale philosophy, but explicitly acknowledge that no detailed calculation procedure is provided for each scale. Without operational detail at the joint scale and a defined transfer mechanism to the frame scale, such frameworks cannot be directly applied in structural design practice. Machine learning approaches produce probabilistic fragility curves efficiently but sacrifice the deterministic force-deformation traceability that EN 1991-1-7 requires for CC3 buildings. They are inherently data-driven and cannot extrapolate reliably to structural configurations or loading scenarios outside their training distribution, which limits their applicability to novel timber systems or non-standard joint typologies. Stochastic debris-tracking and whole-building collapse models operate at a scale too coarse to resolve individual joint behaviour and, therefore, cannot predict the sequence of joint failures or the catenary mechanism that governs load redistribution in post-and-beam frames during column removal. Existing detailed finite element guidelines for timber structures address monotonic loading of standard members and do not include post-peak joint behaviour, element erosion, or the large-displacement catenary effects that are essential for column-removal robustness assessment. The common limitation across all existing approaches is therefore the absence of a validated, operationally complete workflow that bridges joint-scale fracture mechanics and frame-scale progressive failure analysis within a single framework aligned with the second-generation Eurocode design parameters. This is the gap the present paper addresses. Table 1 summarises information about existing approaches and underlines their limitations. The novelty of current work lies not in the application of the finite element method to timber structures, which has long been standard practice, but in four specific and, taken together, unique contributions that address the existing gap in the strength analysis of timber buildings with column-and-beam structures: (1) a fully functional and validated two-scale FEA algorithm with an explicit procedure for switching between scales; (2) a physically based modification of the MAT 143-type Hashin failure criterion; (3) a non-linear connection element with six degrees of freedom, taking into account the coupled direction of degradation and erosion of the element; and (4) the first model, which has been explicitly verified for compliance with the second-generation characteristic value conversion factors specified in Eurocode 5. The practical engineering value of the proposed methodology lies in the fact that it constitutes a two-level design tool. Modelling at the joint level, which requires specialist knowledge of the finite element method and a minimum amount of calibration data (failure energy perpendicular to the grain and stiffness of the ‘wood-steel’ joint), allows for the derivation of verified joint characteristic curves. These directly parameterise the models at the structural frame level, enabling strength assessment in accordance with EN 1991-1-7, a regulatory requirement for CC3-class buildings, which none of the previously available tools could fulfil in a fully operational, experimentally verified and Eurocode-compliant form. 2. Materials and Methods 2.1. Experimental Reference Cases and Structural Configurations The validation of FEM modelling of dowel-type joints is performed using the experimental data shown in [ 15]. Two kinds of samples are analysed: one with 12 dowels ( Figure 1) and the other with eight bolts and 20 dowels. Perpendicular to the bolts, there are six glued-in rods (GIR) ( Figure 2). The timber cross-section was made of GL28h with size 315 × 405 mm 2. Two steel plates with a thickness of 10 mm and a yield strength of 700 MPa were installed within the timber cross-section and connected to the timber through predrilled holes and 20 mm diameter dowels/bolts. An axial tensile load was applied to steel plates, and the support reaction of timber was measured. Material properties of timber are given in Table 2. Additionally, contact elements were also used between steel and timber. The normal and shear stress–displacement curves are shown in Figure 3. It was assumed that the contact elements work only in compression, since there is no glue. The exception was the GIR elements, which had fully glued contact with the timber. Simulation of column removal in frame geometry is shown in Figure 4. Comparison of empirical and simulated strengths for designated samples is demonstrated in Table 3. The following describes the geometry of the two-span timber post-and-beam frame used for column removal validation. Column cross-section: 90 × 90 mm; beam cross-section: 150 × 31.5 × 2 mm (twin members); span: 2000 mm (each bay); bolted joints: 7 bolts d10 at each beam-column connection. Pin supports at outer column bases allow rotation but prevent translation in all directions. The central column (dashed outline) is instantaneously removed to simulate the column-loss scenario; the resulting load must be redistributed to the two outer columns through the beam-joint system. The joints immediately adjacent to the removed column (highlighted) are the primary load redistribution path and are expected to reach their ultimate rotation and displacement capacity first. Experimental data from Cheng et al. [ 4]. 2.2. Timber Constitutive Model (Joint Scale) The MAT 143-type quadratic Hashin criterion is selected over the nine alternatives reviewed in Table 4Table 5 on the basis of three criteria: physical appropriateness, parameter availability, and numerical stability. Physically, the quadratic formulation explicitly captures stress interaction under the multiaxial stress states that govern dowel-joint failure—a capability absent from maximum-stress models, which treat each stress component independently and predict premature failure at local stress concentrations without accounting for combined loading effects. Compared to WoodST, the MAT 143-type model avoids temperature-dependent parameters irrelevant to ambient-temperature mechanical loading, and requires no biaxial interaction parameters beyond those available in EC 5—while providing equivalent accuracy for the two-axial stress states encountered in dowel connections. Compared to nonlocal damage models, the MAT 143-type model achieves mesh-objective results through fracture energy regularisation using physically measured fracture energies ( Table 4) without introducing a nonlocal interaction radius that has no EC 5 equivalent and requires dedicated calibration. The modification proposed in Equation (2)—removing the stress product term σ 22 σ 33 from the perpendicular-to-fibre yield surface—further improves physical correctness by eliminating a non-physical artificial enhancement of shear strength under combined radial-tangential stress states not observed experimentally, while simultaneously improving Newton–Raphson convergence stability in the post-peak softening regime. Similar material models to “MAT 143” in LS-DYNA [ 16] are used. However, the following changes have been made to achieve a more robust and applicable model. More detailed explanations of the changes made are available in the publication [ 23]. Also, in addition to the changes mentioned in [ 23], an update is proposed regarding the perpendicular-to-fibre yield surface. f p = σ 22 + σ 33 2 f 90 2 + σ 23 2 − σ 22 σ 33 f r o l l 2 , (1) To the following version: f p = σ 22 + σ 33 2 f 90 2 + σ 23 2 f r o l l 2 , (2) where stress multiplication σ 22 σ 33 is omitted, since it can create unphysical effects; in case σ 22 σ 33 has a positive value, then it “artificially” improves shear strength (and also strength perpendicular to fibres), which is not experimentally observed and is, therefore, changed. The rest of the timber constitutive model is the same as shown in [ 23]. 2.3. Global-Scale Frame Numerical Models Multifiber beam finite element analysis and finite element model of post-and-beam using non-linear joint elements are shown in Figure 5Figure 6. The non-linear joint element local stiffness matrix is assembled in the following way: k 11 0 0 k 22 0 0 0 0 0 0 0 0 0 0 0 0 k 33 0 0 k 44 0 0 0 0 0 0 0 0 0 0 0 0 k 55 0 0 k 66 , (3) K l o c a l = k l o c a l 11 k l o c a l 11 k l o c a l 11 k l o c a l 11 , (4) where k 11 denotes axial stiffness [ N m ]; k 22 , k 33 denote shear stiffness [ N m ]; k 44 denotes torsion stiffness [ N m r a d ]; and k 55 , k 66 denote bending stiffness [ N m r a d ]. The stiffness parameters k i i are defined as an elastic–plastic model ( Figure 7). Direct interactions between different load directions are implemented. A factor multiplication rule is used, and each stiffness matrix coefficient is influenced by interaction factors: k 11 = k 11 * · f 21 · f 31 · f 41 · f 51 · f 61 … k 66 = k 66 * · f 16 · f 26 · f 36 · f 46 · f 56 , (5) where k i i is the initial stiffness when no failure has happened and f k i is the strength reduction factor from damage in the k-th load direction. The interaction matrix is shown below in Table 6. The simplified interaction matrix retains only the three physically dominant couplings for planar post-and-beam frames under column-removal loading. The neglected off-diagonal terms—particularly shear–bending ( f 25 , f 35 ) and torsion–bending ( f 45 , f 54 ) interactions—become relevant in short-span configurations with high shear-to-moment ratios, and in three-dimensional frames with eccentric or asymmetric loading. For such cases, the interaction matrix should be recalibrated using dedicated joint-level experiments or detailed 3D FEM simulations. The sensitivity of the results to the bending–axial coupling factor f 51 = f 61 = 0.5 was evaluated by bracketing with values of 0.25 and 0.75; only the adopted value of 0.5 reproduces the experimental catenary force measurements, confirming both the physical rationale and the empirical validity of this assumption. For beam and post timber frame with dowel-type joints, the following influence functions are used: If in axial tension load axial stiffness k 11 drops to 0 and it starts to yield, then stiffness in all other directions also drops to f 12 = f 13 = f 14 = f 15 = f 16 = 0 . If bending stiffness k 55 drops to 0, and it starts to behave as a plastic hinge, then axial stiffness is reduced by half f 61 = f 51 = 0.5 . Here, 0.5 accounts for the fact that tensile or compressive zone elements typically yield first, while the rest of the section remains elastic. If, in any direction, the critical displacement or rotation limit is reached and the material is fully damaged in that direction, then it affects all other directions, and all influence functions are set to 0. For different connection types, these influence parameters should be calibrated to capture all relevant influences for a specific problem. 2.4. Scale-Bridging Strategy The two-scale model proposed in this article operates on the basis of a systematic one-way relationship between the joint scale and the frame scale. This subsection describes the workflow by which the results of continuous damage mechanics simulations at the connection level are reduced to parameters that determine the behaviour of non-linear connection elements at the framework level, providing specific details linking the different scales—precisely those details that were considered necessary in previous holistic models, but were not implemented. Step 1—Joint-Scale Simulation For a given type of timber joint (timber species, quality class, steel plate configuration, and the diameter and arrangement of the dowels), a three-dimensional finite element simulation is carried out at the joint scale in accordance with the procedure described in Section 2.2. The simulation uses a modified Hashin failure criterion of the MAT 143 type with a corrected yield surface perpendicular to the grain (Equation (2)), contact elements for the cohesion zone at the wood-steel interface, and an explicit dynamic solver with kinetic damping. The main result of this simulation is the complete ‘force–deformation’ curve for each active load direction—axial tension, shear, and bending moment depending on the corresponding displacement or rotation—calculated from zero load through the peak load and up to the post-peak softening regime. Step 2—Extraction of Joint Element Parameters The following parameters are directly extracted from the ‘force–strain’ curves at the joint level to populate the formulation of the non-linear joint element at the framework level (Equations (3) and (4)): Initial elastic stiffness k i i in each active degree of freedom: axial stiffness k 11 [ [N/m], shear stiffnesses k 22 , k 33 [N/m], and bending stiffnesses k 55 , k 66 [Nm/rad]. These are read from the initial linear slope of each force–deformation curve, corresponding to the SLS stiffness range. Ultimate load-bearing capacity in each direction: the maximum force or moment that can be withstood before plastic deformation begins, corresponding to the ultimate limit state (ULS) of the connection. Ultimate displacement or rotation capacity: the deformation at which the force–deformation curve reaches a defined residual strength threshold (taken as zero for brittle splitting-governed failure modes, consistent with the element erosion criterion). Slope of the strength decay curve after peak load is reached: the rate at which strength decreases after peak load is exceeded, which determines the degree of plasticity available for load redistribution during progressive failure. In cases where the numerical behaviour after the peak load is reached is indeterminate or mesh-dependent, this slope is determined in accordance with the principles of fracture energy set out in Eurocode 8 (prEN 1998-1-2:2023 [ 14]), as described in Section 2.2. For the column-removal design validated in this study, these parameters correspond directly to the values, which were obtained from joint-level simulations of specimens C1 and C4 or derived from the experimental measurements reported earlier. Step 3—Assignment of Directional Interaction Rules Once the frame-level connection element has been assigned parameters for elastic stiffness, ultimate strength, and ductility, the rules for directional interaction set out in Table 6 are applied. These rules determine how plastic deformation or damage in one loading direction alters the stiffness in all other directions, and constitute the mechanism by which the physics of failure at the joint level—described in detail in the three-dimensional continuum damage model—is represented in a computationally efficient form at the framework level. For timber frames with beam-and-column floors and dowel-connected joints, the three rules of active interaction are defined as follows, based on the failure mechanisms identified during joint-level modelling: Axial tension governs all directions: when axial tensile yielding is reached—the critical failure mode in the column—removal catenary phase-all directional stiffnesses are simultaneously set to zero. This reflects the experimental observation that tensile splitting or net-section failure of the timber cross-section at the dowel row destroys the load-carrying capacity of the joint in all directions simultaneously. Bending plastic hinge reduces axial stiffness by half ( f 51 = f 61 = 0.5 ): when bending stiffness k 55 drops to zero, axial stiffness is reduced to 50% of its elastic value. This accounts for the partial yielding of the cross-section at plastic hinge formation: the tensile or compressive fibres yield first while the central fibres remain elastic, so approximately half the cross-sectional area continues to contribute to axial load transfer. Full damage triggers element erosion: when the critical displacement or rotation capacity is reached in any direction, all interaction functions are set to zero and the joint element is removed from the model. This represents the brittle, total loss of load-carrying capacity observed experimentally at 390 mm vertical displacement. For joint types other than the pin joints considered in this study, the interaction rules should be recalculated using specific joint-scale analyses that take into account the specific failure mechanisms of the joint in question. Step 4—Frame-Scale Progressive Failure Analysis Once the non-linear connection elements have been fully parameterised, a global model of the frame ( Section 2.3) is assembled and subjected to a loading sequence involving the removal of columns. Multistrand beam elements account for the axial and bending non-linearity of timber, whilst non-linear connection elements determine the redistribution of load between elements as connections gradually fail. The element erosion method allows the analysis to continue after the failure of individual joints, tracking the development of the chain line force and the sequential failure of the remaining joints until a stable equilibrium is reached or the structure collapses completely. 2.5. Distinction from Prior Holistic Approaches The procedure for scaling described above differs from previously existing multiscale models of timber frames in three aspects that are critical from a practical point of view. Firstly, the transition from the joint scale to the frame scale is quantitative and direct: the stiffness, strength, and deformation capacity values at the joint level, calculated using a three-dimensional finite element model, are used as numerical input data for the joint elements at the frame level without any intermediate empirical fitting or engineering estimates. Secondly, the procedure undergoes independent validation at both scales: results at the joint level are checked against experimental ‘load–displacement’ curves and EC 5 characteristic values before being transferred to the frame level, which ensures that no errors are inadvertently carried over between scales. Thirdly, the procedure is practically applicable: all input parameters at the connection level are either standard material properties according to EC 5 or minimally calibrated failure parameters (failure energy perpendicular to the grain, support stiffness of the contact element), and all parameters at the frame level are expressed in physically meaningful units that can be directly interpreted by structural engineers. This combination of quantitative scale transfer, independent validation and practical transparency of parameters constitutes an operational contribution that distinguishes the proposed framework from holistic but procedurally incomplete approaches previously proposed in the literature. 3. Results 3.1. Non-Linear Modelling of Dowel-Type Joints Simulation results and comparison with experimental measurements are given in Figure 10. Only minimal material property calibration was performed for sample C1 (the contact element stress–displacement curve and the fracture energy perpendicular to grains were calibrated). Sample C4 was simulated without any other calibration, and it also shows good agreement with experimental data. Failure mechanisms were also analysed and compared with experimental observations. For sample C1, there was initially tensile splitting failure perpendicular to grains, followed by partial axial failure of the timber section located at the last row of dowels (see Figure 11b). The section tended to open timber edge parts, which was also experimentally observed (see Figure 11b). Better results were obtained for the C4 sample since bolts and perpendicular GIR did not allow any significant tensile splitting damage; therefore, timber axial tensile damage at the cross section located on the last row of dowels was mainly observed (see Figure 12). The validity of the explicit dynamic solver for quasi-static problems was verified by monitoring the ratio of kinetic to internal strain energy throughout each simulation. This ratio remained below 5% at all load increments, confirming that inertial forces were negligible. The smoothness of the resulting load–displacement curves ( Figure 10) provides additional confirmation that no dynamic artefacts contaminated the static solution. Specimen C1 under peak load (1455 kN, displacement ~3.0 mm). In Figure 11a is shown deformed shape (displacement magnified 10-fold for clarity), demonstrating the opening of the wood edge zone in the last row of dowels—which corresponds to the experimentally observed splitting during tensile loading perpendicular to the grain. (In Figure 11b is shown distribution of damage severity (0 = undamaged, 1 = completely damaged), demonstrating Type I splitting, starting at the bearing surface of the last row of dowels and propagating perpendicular to the grain direction towards the edge of the timber. The damage zone is localised within approximately one pitch between the dowels of the last row of fasteners, confirming that the failure is due to stress concentration at the outermost dowel, rather than distributed plastic deformation along the entire length of the joint. The damage variable distribution in Figure 11b confirms that the governing failure mechanism for C1 is tensile splitting perpendicular to grain, initiating at the bearing zone of the last dowel row where the perpendicular-to-grain tensile stress first reaches the fracture energy-regularised damage threshold. The opening deformation visible in Figure 11a—where the timber edge zone separates from the main section—replicates the experimentally observed behaviour reported in [ 15], in which edge splitting was the primary failure mode for connections without transverse reinforcement. This failure mode is correctly captured by the modified MAT 143-type failure criterion (Equation (2)), specifically by the perpendicular-to-grain tensile damage component of the yield surface. Specimen C4 at peak load (2140 kN). In Figure 12a is shown deformed shape showing uniform axial elongation without visible edge opening, in contrast to the edge splitting observed in C1 ( Figure 11a). The perpendicular glued-in rods (GIR) provide transverse confinement that prevents splitting crack propagation. In Figure 12b is shown damage variable distribution showing axial tensile damage (parallel-to-grain) concentrated at the cross-section of the last dowel row, consistent with net-section tensile failure. The absence of perpendicular-to-grain splitting damage confirms that the GIR rods successfully redirect the failure mode from splitting-governed to strength-governed behaviour. Comparing Figure 11Figure 12 directly illustrates the fundamental effect of transverse splitting reinforcement on connection failure mode. In C1, damage is dominated by mode I perpendicular-to-grain fracture (splitting), producing localised edge separation at low energy dissipation. In C4, the GIR rods suppress this splitting mechanism entirely, forcing failure into the axial tensile net-section mode at a load 47% higher than C1 (2140 kN vs. 1455 kN). This failure mode transition—from splitting-governed to strength-governed—is the physical basis for the significant difference in rope effect contribution between the two specimens and is correctly reproduced by the proposed constitutive model without recalibration between specimens. 3.2. Analytical Load–Displacement Curves for Typical Joints According to EC 5 To better understand how numerical results align with the second generation of Eurocode 5 (EC5), comparisons are made among numerical, experimental, and EC5 results. The results of the C1 sample are given in Table 7, and for the C4 sample in Table 8. Experimental and simulation results for strength and stiffness differ by less than 5%. However, the characteristic strength of the C1 sample is approximately 10% lower than the experimental value. It should also be noted that the characteristic stiffness of C1 is around 2.45 times smaller than experimental measurements. This is in agreement with EC 5 Annex A.3.3.2 statement. Similar conclusions can be made for the C4 sample. It is also evident that for C4, the bolts have generated a significant rope effect, and perpendicular glued-in rods have prevented timber splitting, leading to a total strength similar to the net-tensile strength of timber. In this case, both analytical EC5 and numerical methods yield results that are reasonable compared to the experimental results. Stiffness at different load levels predicted by numerical simulations is in good agreement with the EC5 mean values. 3.3. Column Removal Modelling with FEM To see how the assembly of timber frame elements works in a column-removal scenario, the experiment from [ 4] is replicated in a numerical model. The frame displacement is shown in Figure 13. The model is created in the Versim4D software using special non-linear joint elements and multi-fibre beam elements for timber. The theory behind those elements is described in the previous section. A deformed frame with a central vertical displacement of 300 mm (scale 1:1, no magnification), corresponding to a rotation at the mid-span of approximately 8.5°—within the hinge’s maximum rotational capacity of 10° ( Table 9). The colour scale reflects the magnitude of the vertical displacement (from 0 mm at the outer supports to −300 mm in the middle of the span). At this level of deformation, the frame has transitioned from a mode dominated by bending to one in which the load is distributed between bending and torsional deformation: the beams act partly as tension members, generating significant axial forces in the connections. All four beam-to-column connections adjacent to the distant column have entered the plastic zone but have not yet reached their ultimate capacity for rotation or displacement. The 1:1 deformation scale illustrates the significance of geometric non-linearity: the displacement at the mid-span (300 mm) is 15% of the span length (2000 mm), confirming that small-displacement assumptions are entirely unsuitable for this loading regime. Figure 15 shows the ‘load–displacement’ curves for three cases: numerical simulation, a quasi-static experiment and a dynamic experiment. The ratio of peak dynamic displacement to peak quasi-static displacement at equivalent load levels serves as a direct empirical estimate of the dynamic amplification factor for this structural configuration. The ratio reduction from the tested frames can be explained by non-linear energy dissipation in the timber joints—damage, plasticity, and friction at the interface between the dowel and the timber—which provides effective damping that limits inertial amplification. The trajectory of the quasi-static simulation lies within the envelope bounded by the static and dynamic experimental curves across the entire displacement range up to 390 mm, confirming that the quasi-static model captures the essential behaviour of the structure, whilst the dynamic amplification remains within the limits demonstrated by the experimental comparison. Axial force compared to vertical displacement for the computational and experimental analysis of the frame column removal scenario is shown in ( Figure 16) 4. Discussion The proposed joint-scale model reproduces the experimental strength and stiffness of GL28h glued-laminated timber dowel joints with an error of less than 5% across two distinct configurations. One of the key quantitative results is a clear verification of compliance with the characteristic values of EC 5 (prEN 1995-1-1:2023). The characteristic stiffness predicted by EC 5 for specimen C1 is 226 kN/mm (SLS), whereas the average value obtained from the simulation is 565 kN/mm—giving a ratio of 2.50, which almost exactly matches the conversion factor Rm = 2.5Rk prescribed in Annex A4.2 to EC 5. Such direct normative traceability is absent in all the multiscale models of timber strength previously reviewed in the literature. In addition to the direct validation against specimens C1 and C4, the consistency of the proposed model with the characteristic value conversion factors of EC 5 (prEN 1995-1-1:2023) Annex A4.2 provides implicit validation against the extensive experimental database underlying the EC 5 design equations. Since EC 5 dowel strength predictions are calibrated against hundreds of connection tests across a wide range of timber grades and fastener configurations, the demonstrated agreement between numerical mean predictions and EC 5 characteristic values suggests that the framework captures the essential physics of dowel-type connections beyond the two directly simulated configurations. At the structural level, the multi-fibre beam model with non-linear connection elements accurately reproduces the full non-linear load–displacement relationship for vertical displacements of up to 390 mm, including: (1) the correct sequence of joint failure (complete failure of two joints at a displacement of 390 mm, which corresponds to experimental observations); (2) the dynamics of changes in the axial force in the chain, with the simulation results matching the experimental data across the entire displacement range from 0 to ~250 mm; and (3) non-linearity of the ‘load–displacement’ relationship under both quasi-static and dynamic experimental conditions. This represents a significant quantitative breakthrough compared to the stochastic debris tracking model proposed by Cao et al. [ 8], which focuses on the collapse of the entire building but does not account for the ‘load–displacement’ trajectory at the joint level, as well as the machine-learning-based fragility assessment approach proposed by Alshbul et al. [ 5] which provides probabilistic fragility curves but does not ensure deterministic ‘stress–strain’ traceability. A quantitative comparison of the proposed approach with bibliographic data is presented in Table 10. The differences between EC 5 and numerical results reflect three distinct sources: the characteristic-to-mean statistical conversion (factor 2.5, consistent with EC 5 Annex A4.2); failure mode representation (EC 5 governs on dowel yielding while the simulation captures splitting-governed failure for C1); and the absence of post-peak information in EC 5. For practical engineering design, these differences have the following implications. EC 5 characteristic values remain appropriate for ULS and SLS design of individual joints with standard failure modes and adequate safety margins. However, for robustness assessment under column-removal scenarios—where mean stiffness governs load redistribution, mean strength determines residual capacity, and post-peak ductility determines survival—the numerical model provides information that EC 5 cannot. Specifically, using EC 5 characteristic stiffness as input to a frame-level robustness model would overestimate flexibility by a factor of approximately 2.5 and could misidentify the failure sequence. The proposed numerical framework should therefore be seen as complementary to EC 5: EC 5 governs routine joint design, while the numerical model enables the advanced non-linear robustness assessment now required for CC3 buildings under EN 1991-1-7. From a practical engineering perspective, the proposed framework is best deployed as a two-tier tool: joint-scale simulations, requiring specialist FEM expertise and minimal calibration data (fracture energy perpendicular to grain and timber-steel interface bearing stiffness), produce validated joint property curves that parameterise frame-scale models; frame-scale progressive failure analysis, requiring expertise comparable to advanced non-linear structural analysis, then enables robustness assessment in compliance with EN 1991-1-7. The framework is validated for planar GL28h post-and-beam frames with dowel-type connections and should be extended and re-validated before application to three-dimensional frames, alternative timber products, or non-dowel connection systems. The use of mean—rather than EC 5 characteristic—stiffness and strength values as inputs to robustness models is strongly recommended, as characteristic values are intended for ULS design with partial safety factors and will systematically misrepresent load redistribution behaviour in progressive collapse analysis. 5. Conclusions An efficient numerical modelling framework is proposed and validated against experimental data at the structural joint and building frame scales. Structural joints with steel dowels were accurately modelled with continuum damage mechanics. Based on a literature review of published timber constitutive models, a modified Hashin-type quadratic criterion was selected, which proved to be numerically stable, efficient, and accurate for dowel-type joints. An effective two-scale finite element modelling scheme for assessing the stability of solid timber frames with columns and beams under conditions of column removal was developed and verified against experimental data. The choice of the MAT 143-type Hashin quadratic criterion over alternative criteria based on peak stress was motivated by the fact that it possesses a differentiable failure surface—which ensures convergence stability in post-peak stress regimes—as well as by the fact that its most important parameters correspond exactly to the design values specified in Eurocode 5. The proposed modification of the yield surface, perpendicular to the fibres, further enhances physical accuracy and numerical stability under combined radial-tangential stress states occurring in the support zones of dowel connections. The explicit comparison with EC 5 second-generation characteristic values demonstrates that numerical mean predictions are consistent with normative values via the statistical conversion R m = 2.5 R k , establishing the framework as a transparent and auditable extension of—rather than an alternative to—current design standards. Designers are advised to use EC 5 for routine ULS and SLS joint design, and the proposed numerical framework for robustness assessment, non-standard joint configurations, and large-displacement catenary analysis where EC 5 provides no applicable method. On the scale of the entire frame, multi-fibre beam elements combined with non-linear joint elements with six degrees of freedom accurately reproduced the full non-linear ‘load-displacement’ response of the column-and-beam frame upon removal of the columns. The model correctly predicted the sequence of joint failure, recording complete joint failure at a vertical displacement of 390 mm, which is consistent with experimental observations. The directional interaction functions implemented in the non-linear joint elements are calibrated for planar frames with dowel-type joints under dominant bending-to-catenary loading. Extension to three-dimensional spatial frames, short-span configurations with high shear forces, or joint types with significant torsional capacity will require recalibration of the off-diagonal interaction terms in Table 6, which is identified as a priority for future work. The catenary effect—a critical mechanism for load redistribution during large displacements—was quantified using curves showing the relationship between axial force and vertical displacement, which demonstrate close agreement with measured data across the entire deformation range up to 300–390 mm. The stiffness and strength of dowel-type joints were compared with those of second-generation Eurocode 5, and the results were in good agreement with experiments and simulations. At the whole frame scale, efficient numerical models are proposed that can capture the following effects: ▪ Timber axial and bending failure with multi-fibre finite element technique; ▪ Non-linear timber joint behaviour, including post-failure effects, such as bending–axial–shear strength interaction; ▪ Large displacement effects and catenary effect with accurate prediction of axial forces. The joint-scale model has been validated against two GL28h glued laminated timber dowel connection configurations representing the two physically distinct failure modes of splitting-governed and net-section tension-governed behaviour. Although the zero-recalibration transferability from C1 to C4, the consistency with EC 5 second-generation characteristic values, and the independent frame-scale validation against Cheng et al. [ 4] provide confidence in the framework’s physical correctness, broader generalisation to other timber grades, dowel diameters, connection geometries, and loading protocols requires additional validation. This is identified as a priority for future work, alongside the extension to three-dimensional spatial frames and cyclic loading scenarios. The frame-level model is verified quasi-statically, with dynamic amplification being assessed by comparison with data from dynamic experiments. For a complete assessment of a building in accordance with EN 1991-1-7, which requires explicit dynamic non-linear analysis taking into account structural dynamics and the interaction of loads with the structure, the proposed methodology must be extended to include physical mass and damping matrices in the formulation of multistrand beams and joint elements. The current quasi-static implementation is suitable for preliminary strength assessment and for configurations in which the dynamic amplification factor can be limited based on experimental data or energy-based analysis, but should not be used as a substitute for explicit dynamic analysis in cases where significant inertial amplification is expected—particularly for lightweight timber frames with low damping and high natural frequencies.” Author Contributions Conceptualization, J.S. and A.B.; methodology, J.S.; software, J.S.; validation, J.S., A.B., E.B. and D.S.; formal analysis, V.L.; investigation, J.S.; resources, V.L.; data curation, J.S.; writing—original draft preparation, J.S.; writing—review and editing, J.S., A.B. and V.L.; visualization, J.S.; supervision, J.S. and D.S.; project administration, J.S.; funding acquisition, J.S. All authors have read and agreed to the published version of the manuscript. Funding This research was funded by the long-term Latvian national research program project “Biomedical and Photonics research platform for innovative products”, grant name “Development of a Virtual Testing Laboratory for Wood Materials and Structures”, acronym VIWO-LAB. The part of research work was developed by A.Berzins within the framework of the EU ERDF-funded project “RTU Doctoral Grants for Supporting Scientific Excellence in Smart Specialization Areas” (No. 1.1.1.8/1/24/I/007) within the framework of a doctoral grant (ID 8004). Data Availability Statement The raw data on which the conclusions of this article are based will be provided by the authors upon request. Acknowledgments During the preparation of this manuscript, the authors used ChatGPT for the purpose of English grammar check. The authors have reviewed and edited the output and take full responsibility for the content of this publication. Conflicts of Interest The authors declare no conflicts of interest. Abbreviations The following abbreviations are used in this manuscript: FEM Finite Element Modelling EC 5 Eurocode 5 prEN 1995-1-1:2023 EC 8 Eurocode 8 prEN 1998-1-2:2023 CC3 Consequence class 3 according to Eurocode MAT 143 Material constitutive model number 143 in LS-DYNA software Annex A4.2 Annex to Eurocode 5 prEN 1995-1-1:2023 References Figure 1. Geometry of dowel-type joint C1. Figure 1. Geometry of dowel-type joint C1. Figure 2. Geometry of dowel-type joint C4. Figure 2. Geometry of dowel-type joint C4. Figure 3. Cohesive zone stress–displacement curves in ( a) normal and ( b) shear directions. Figure 3. Cohesive zone stress–displacement curves in ( a) normal and ( b) shear directions. Figure 4. Geometry of frame with column removal simulation. Experimental data given in [ 4]. Figure 4. Geometry of frame with column removal simulation. Experimental data given in [ 4]. Figure 5. Multifiber beam finite element. ( a) schematic representation of a beam finite element and cross-section with fibres. ( b) A bilinear stress–strain curve is used in this finite element. Figure 5. Multifiber beam finite element. ( a) schematic representation of a beam finite element and cross-section with fibres. ( b) A bilinear stress–strain curve is used in this finite element. Figure 6. Post (in blue) and beam (in green) finite element model with non-linear joint elements (in red). Figure 6. Post (in blue) and beam (in green) finite el
