Multi-model collaborative finite element approach for service-state analysis of timber components

Abstract

Current numerical models for timber components are often limited to single constitutive theories, making it difficult to accurately simulate their complex multi-stage mechanical behavior under diverse service conditions. To overcome this limitation, this study proposes an innovative “multi-model collaborative finite element analysis method.” Guided by the principle of “service-condition matching,” this method dynamically selects and integrates appropriate mechanical models to achieve high-fidelity simulation throughout the entire service life of the component: an orthotropic elastic model is used to reveal the “strong longitudinal but weak transverse” stress distribution under normal loads; the Hill anisotropic criterion captures the evolution of plastic strain as loads approach the yield point; and a viscoelastic model describes the rate-dependent and stress-relaxation behaviors under long-term loading. Results show that the collaborative method effectively elucidates the respective mechanical mechanisms, and Digital Image Correlation (DIC) measurements validate the simulation accuracy. The proposed method provides an innovative and efficient cross-scale numerical tool for timber structures, enabling integrated simulation from short-term safety assessment to long-term performance evolution. This is of great significance for the conservation, performance optimization, and lifespan prediction of historic timber components.

Full Article

A Multi-Model Collaborative Finite Element Approach for Service-State Analysis of Timber Components

Le Zhou  , Xiaoyi Hu  ,* Hongchao Liu  , and Xin Hou 

Current numerical models for timber components are often limited to single constitutive theories, making it difficult to accurately simulate their complex multi-stage mechanical behavior under diverse service conditions. To overcome this limitation, this study proposes an innovative “multi-model collaborative finite element analysis method.” Guided by the principle of “service-condition matching,” this method dynamically selects and integrates appropriate mechanical models to achieve high-fidelity simulation throughout the entire service life of the component: an orthotropic elastic model is used to reveal the “strong longitudinal but weak transverse” stress distribution under normal loads; the Hill anisotropic criterion captures the evolution of plastic strain as loads approach the yield point; and a viscoelastic model describes the rate-dependent and stress-relaxation behaviors under long-term loading. Results show that the collaborative method effectively elucidates the respective mechanical mechanisms, and Digital Image Correlation (DIC) measurements validate the simulation accuracy. The proposed method provides an innovative and efficient cross-scale numerical tool for timber structures, enabling integrated simulation from short-term safety assessment to long-term performance evolution. This is of great significance for the conservation, performance optimization, and lifespan prediction of historic timber components.

DOI: 10.15376/biores.21.1.1497-1514

Keywords: Mechanical model; Finite elements; Simulation; Wood components; Applications; DIC

Contact information: School of Optical-Mechanical Engineering, Zhejiang A&F University, Hangzhou, Zhejiang 311300 China; *Corresponding author: jimhxy08@zafu.edu.cn

Graphical Abstract

INTRODUCTION

Timber has gained significant prominence in modern construction engineering due to its environmental friendliness and renewability. Its exceptional mechanical properties, particularly its high strength-to-weight ratio, make it highly advantageous for load-bearing and flexural structural applications. Wooden members serve as fundamental components within timber structural systems, encompassing processed elements such as beams, columns, rafters, and other primary load-carrying units (as shown in Fig. 1). Collectively, these elements form the structural skeleton of buildings, bearing loads and providing essential structural support. In-depth research on these critical components establishes an essential foundation for their effective conservation, accurate performance evaluation, and reliable restoration.

Fig. 1. Schematic diagram of partial wooden components in timber structures

China has the world’s most enduring and systematically sophisticated tradition of timber architecture (Pan 2003). Since pioneering research by scholars such as Liang Sicheng and Liu Dunzhen in the 1930s, this field has evolved into a comprehensive research chain spanning material properties, component behavior, and structural systems.

At the material level, research primarily focuses on  the damage mechanisms and mechanical responses of wood. Key research domains encompass: strain distribution and damage evolution in wood containing defects (Berg et al. 2015; Wen and Zhao 2016; Chen 2017); the influence of cracks on seasonal flexural behavior (Bisong et al. 2024); the effectiveness of strengthening damaged timber columns (Shi 2023); and the development of hygro-thermal coupled constitutive models (Fortino et al. 2019). At the component and connection level, studies concentrate on fundamental mechanical performance: Examples include the static characteristics of traditional Chinese timber frames and Dou-gong brackets (Li 2023), the performance degradation of cracked beams and columns (Wang 2023; Yu 2022), and the validation of long-term safety performance for critical components, such as the Caibu-jin beam in the Hall of Supreme Harmony (Zhou and Zhang 2007). In the area of strengthening techniques, research has confirmed the effectiveness of Carbon Fiber Reinforced Polymer (CFRP) (Hassan et al. 2020; Işleyen et al. 2021). Further investigations have been conducted on the shear performance of connection joints and the lateral stability of structural systems (Mansor et al. 2025; Wang et al. 2025). At the structural system level, emphasis is placed on the global mechanical behavior of timber structures. Zhang (2022) analyzed the static and dynamic responses of arch-beam systems; Gomon et al. (2023) developed computational methods for nonlinear deformation applicable to timber structures. Through numerical simulations, Yu et al. (2024) and Pan et al. (2024) revealed the lateral resistance mechanisms and progressive collapse behavior of traditional timber frames, respectively. Xue et al. (2021) proposed effective strategies for optimizing the lateral performance of glued laminated timber frame-infill wall systems. However, these ancient timber structural components, having endured hundreds or even thousands of years of weathering, face significant challenges: they are not only subjected to sustained static loads and material degradation due to environmental aging but may also experience extreme events, such as earthquakes and strong winds, leading to plastic deformation or cyclic loading. These complex service conditions and their interactions are key mechanisms contributing to performance degradation, damage accumulation, and even failure. Accurately simulating the mechanical behavior under these conditions is crucial for assessing their current safety status, predicting remaining service life, and formulating scientific reinforcement and repair strategies. Despite significant progress in characterizing the mechanical behavior of timbers, the inherent heterogeneity, geometric variability, and loading complexity of wood continue to present major challenges for precise mechanical simulation. Classical finite element analysis (FEA) approaches typically rely on a single, predetermined constitutive model for the entire analysis process. This “one-model-fits-all” strategy struggles to dynamically adapt to the different mechanical states (e.g., from elastic response, plastic yielding, to long-term creep) a component may experience throughout its service life. Consequently, a significant modeling gap exists in accurately simulating complex, multi-stage service behavior, characterized by insufficient model adaptability and reduced predictive accuracy across different stages.

To address these limitations, this study proposes a multi-model collaborative analysis approach guided by the principle of service-condition matching. This approach systematically integrates key mechanical properties of wood, including elasticity, elastoplasticity, anisotropy, and viscoelasticity, to establish a high-fidelity finite element modeling system. The model was validated under typical loading conditions using four-point bending experiments combined with Digital Image Correlation (DIC) technology, demonstrating high predictive accuracy. Ultimately, this approach provides a robust numerical toolkit for the safety design, service performance optimization, quantitative damage assessment, and repair strategy formulation of timber structures.

EXPERIMENTAL

Methodological Framework

This study proposes an innovative “multi-model collaborative finite element analysis approach”, whose core lies in the principle of “service-condition matching.” Based on the different mechanical states that timber components may experience during actual service, this approach dynamically selects and applies the most suitable mechanical model for simulation analysis. Specifically, the decision for model switching is primarily based on the evaluation of the component’s real-time mechanical response or anticipated load history.

1) When external loads are below the material’s proportional limit and long‑term creep is not a concern, an orthotropic linear elastic model is activated to rapidly evaluate the stress distribution and deformation of the component under normal working loads and to verify compliance with standard requirements.

2) When the load history or real‑time monitoring indicates that stresses are approaching or exceeding the yield threshold, or when the simulation of damage accumulation is required, an elastoplastic model based on the Hill anisotropic yield criterion is employed to accurately simulate plastic strain distribution, ultimate bearing capacity, and damage evolution. In both of the above stages, an isotropic model is used solely for comparative validation, highlighting the advantages of the anisotropic model in practical applications.

3) When analyzing long-term sustained load or stress relaxation effects, the Generalized Maxwell model should be employed; when analyzing cyclic energy dissipation, a hysteresis model incorporating hyperelasticity is required. The parameters of the former are calibrated by fitting a Prony series to stress relaxation data; the parameters of the latter need to be calibrated separately: its hyperelastic part is determined from the stress–strain curve, while its viscoelastic behavior can be determined from stress relaxation curves (Zhang et al. 2019), to predict the long-term deformation, stress relaxation, and energy dissipation characteristics (Fig. 2). This decision-making process based on quantified thresholds ensures that model application is highly synchronized with the actual physical state of the component, thereby achieving a balance between simulation accuracy and computational efficiency.

Fig. 2 . Multi-model collaborative analysis approach

The advantage of this method lies in its “condition-model” matching mechanism, which dynamically selects the appropriate model based on the current or anticipated key service conditions of the component, thereby overcoming the limitations of a single-model approach. Building on this, the study establishes a three-tier collaborative analysis system of “small-deformation elastic verification—elastoplastic damage analysis—viscoelastic time-dependent life prediction.” Subsequent sections will elaborate on the development and parameter determination methods for each sub-model, systematically validate the system through four-point bending tests combined with DIC technology, and ultimately demonstrate its exceptional effectiveness in simulating the complex mechanical behavior of timber components.

Elasticity Model

Integrated with advanced Finite Element Analysis (FEA) technology, the proposed timber component modeling approach is capable of simulating complex material behaviors of wood, including anisotropy, non-linear characteristics, and time-dependent effects. According to the General Code for Timber Structures GB 55005-2021(2021), hereafter referred to as the Code, the permissible deflection limit for load-bearing timber components (such as beams, columns, and trusses) is set as 1/250 of the calculated span or 3 mm, whichever is greater. This standard plays a crucial role in ensuring the structural integrity and serviceability of important historical timber constructions such as palaces, ancient bridges, and pagodas. The analysis and simulation of these load-bearing components form the basis for verifying structural safety, providing scientific support for conservation and restoration, and planning subsequent maintenance and reinforcement measures. Given that under normal service loads, the stress level in timber members is generally expected to remain below the proportional limit with recoverable deformation, the Code permits the use of linear elastic models for simulation. However, such models must accurately represent the inherent anisotropy of wood—that is, the significant variation in mechanical properties along different grain orientations (longitudinal, radial, and tangential). Within the FEA approach presented in this study, wood is modeled as a homogeneous orthotropic material, with its elastic properties defined by engineering constants. The validity of the modeling approach was verified by comparing simulation results under isotropic and orthotropic assumptions (see specific engineering constants in Table 1).

Table 1. Parameters of Anisotropic Wood Members

In structural building systems, wooden members serve as core load-bearing components, primarily supporting vertical loads from upper structures such as roofs and floors. Among these, horizontal members are predominantly subjected to transverse bending stresses under service conditions, while vertical members mainly endure axial compressive stresses. These components work synergistically through frame connections to form a stable combined compression-bending system.

To investigate the critical influence of wood anisotropy on its mechanical behavior, this study developed a three-dimensional numerical model of a standard rectangular clear wood specimen (15 × 20 × 200 mm³), which was validated through four-point bending tests. This testing method employs symmetric loading to effectively eliminate shear effects, producing an ideal pure bending stress state within the mid-span region, making it particularly suitable for high-fidelity characterization of material constitutive behavior and model validation. It is important to emphasize that the use of such clear small specimens is primarily intended for calibrating fundamental material constitutive relationships, validating models, and facilitating the isolated study of core mechanisms such as wood anisotropy, elastoplasticity, and viscoelasticity under controlled conditions. Excluding natural defects, such as knots and the direct influence of size effects, this approach allows the research to focus effectively on the material’s inherent nonlinearity, anisotropy, and time-dependent behavior. This provides an essential foundation for subsequently developing models that account for size effects and natural defects, and for extending the analysis to larger-scale, real-world components with imperfections. In the numerical model, simply supported boundary conditions were simulated using two cylindrical rigid supports (with a span equal to the specimen’s length), and rigid body constraints were applied to ensure displacement compatibility and system stability at the interfaces. The loading system was implemented via a rigid cylindrical indenter applying symmetric displacement-controlled quasi-static loading. This setup ensured uniform stress distribution within the pure bending segment at mid-span, providing an ideal condition for standardized analysis of the influence of wood anisotropy on flexural behavior. The slenderness ratio of the specimen satisfies the basic assumptions of classical beam theory, enabling a 2D model to achieve accuracy comparable to that of a 3D model in predicting mid-span deflection and reconstructing longitudinal stress fields, while significantly improving computational efficiency. To ensure simulation reliability, a mesh convergence study was conducted. By progressively refining the mesh size, it was found that when the element size was less than 1.5 mm, the stress variation was below 2%, indicating result convergence. All subsequent analyses therefore employed a mesh size of 1.0 mm to balance computational accuracy and efficiency. Based on the macroscopic elastic parameters of wood listed in Table 1 (referenced from He 2011 and Zhang 2017), and following the requirements of the code, a quasi-static displacement load (controlled according to a deflection limit of 1/250 of the span) was applied, effectively capturing the mechanical response characteristics of the wooden member under small deformation conditions. The orthotropic linear elastic model was adopted as the default starting model for analysis. It was employed when the evaluation indicated that the external load level was below the material’s proportional limit (e.g., taken as 80% of the yield stress), the loading duration was instantaneous, and the analysis did not involve long-term creep effects. This model was used to efficiently and accurately obtain the initial stress and deformation state of the component.

Elasto-plasticity Model

The limitations of elastic models become increasingly pronounced as deformations in wood components amplify, particularly when approaching or exceeding their elastic limit. This is especially critical for historic timber members in ancient structures subjected to heavy or sustained loads. These aged components— some having survived for centuries or even millennia— are often structurally compromised. Long-term environmental exposure degrades mechanical properties and induces embrittlement, rendering the material more susceptible to entering the plastic phase under stress. Additionally, lateral loads, such as wind or seismic forces, can induce plastic deformation. Upon entering the plastic phase, the stress-strain relationship becomes nonlinear (elastoplastic) and may involve permanent deformation, deviating from the linear recovery assumption of elastic models (based on Hooke’s law). Beyond the yield point, wood typically exhibits work hardening, where the yield stress increases with accumulating plastic strain. Therefore, when preliminary analysis based on the elastic model or monitoring data indicates that the equivalent stress in critical regions of the component approaches or exceeds the material’s yield threshold (e.g., reaching 95% of the yield stress), or when cumulative plastic strain has been detected, the analysis framework should switch from the elastic model to the elastoplastic model. The stress-strain response for isotropic elastoplastic constitutive behavior is described by incremental plasticity theory. To accurately simulate the post-yield behavior of wood, its elasto-plastic parameters must be determined through four-point bending tests. The tests were conducted at a loading rate of 5 mm/min to simulate quasi-static conditions, while the finite element analysis was performed directly under the quasi-static assumption.

In numerical simulations employing an updated Lagrangian formulation, where accumulated strain is equivalent to true strain, the engineering stress-strain curves obtained from the experiments must be converted into true stress-true strain curves. For simplified analyses, an ideal plasticity model may be adopted, which assumes that unconstrained plastic flow occurs at a constant stress level. This is shown in Eq. 1,

 (1)

where σ_e is the Equivalent stress (MPa) and σ_y is the Yield stress (MPa).

Wood exhibits pronounced anisotropic mechanical behavior due to its inherent growth characteristics, with noticeable directional variations in strength and elastic modulus, particularly between longitudinal and transverse grain directions. This anisotropy renders the classical von Mises yield criterion inadequate for accurately describing wood’s yielding behavior. In contrast, the generalized Hill yield criterion effectively simulates wood’s mechanical response along different material principal axes, specifically capturing its anisotropic yielding characteristics. The Hill plastic potential function is given by Eq. 2. Notably, when parameters are assigned the values F = G = H = 0.5 and L = M = N = 1.5, the generalized Hill criterion reduces to the isotropic von Mises yield criterion. Within this framework: Parameters F, G, and H define the anisotropic strength ratios for yield under normal stress along three mutually perpendicular directions. Parameters L, M, and N define the anisotropic strength ratios for shear yield on three mutually perpendicular planes. The Hill anisotropic yield model, widely implemented in finite element analysis using the associated flow rule (Eq. 3), provides an established approach for simulating the elastoplastic behavior of orthotropic materials like wood. Equation 2 is given as,

 (2)

where σ₁₁, σ₂₂, σ₃₃ are the positive stress components of the material in the three principal axis direction; σ₁₂, σ₁₃, σ₂₃ are the shear stress components of the material on the spindle plane; F, G, H represent the differences in yield strength of the material along three orthogonal directions, and L, M, N govern the yield behavior on the three corresponding shear planes. Together, these parameters describe the anisotropic yielding characteristics of wood under complex stress states. In Eq. 3, dε^p is the plastic strain increment, and λ is the plastic multiplier.

 (3)

Table 2. Yield Stress of Pinus sylvestris var. mongolica

 This study established a finite element model of a softwood element under four-point bending. The nonlinear mechanical behavior was analyzed by coupling the Mises isotropic criterion with the generalized Hill anisotropic criterion. Material parameters from Ma (2013, see Table 2) were used. Plastic deformation was simulated via displacement-controlled loading. Elastoplastic damage analysis was performed to assess the wood element’s plastic deformation and ultimate load-bearing capacity. The analysis of such deformation represents a key component of the second phase in developing a multi-model collaborative analysis approach.

Visco-elastic Model

Environmental factors (including wind, precipitation, and fluctuations in temperature or humidity) drive progressive deterioration of wood components under long-term exposure. Wood exhibits viscoelastic behavior under sustained load, characterized by creep and stress relaxation (Rong 2008). Whether to adopt a viscoelastic analysis model should be jointly determined by the load duration and the analysis objectives. This model should be employed when it is necessary to study the creep and stress relaxation behavior of a component under long-term sustained loading (e.g., load duration > 12 h), or to analyze its energy dissipation and path dependency under cyclic loading. To simulate this behavior, this study employs a Generalized Maxwell model (comprising multiple spring-dashpot units in parallel) to develop a finite element model for the viscoelastic response of wood elements (see the center panel of Fig. 3). This model effectively captures relaxation and creep behavior across different time scales. (Note: The model reduces to the classical three-element model when n=1). The constitutive equation can be expressed in integral form by Eq. 4. Given the experimental challenges and minor influence of volumetric relaxation, the research prioritizes shear, tensile, and compressive relaxation tests. Due to limitations in experimental data, Poisson’s ratio is assumed to be constant, implying that the normalized relaxation functions for the elastic modulus and shear modulus are identical. For the stress relaxation experimental data for Pinus sylvestris var. mongolica (adapted from Liu 2011, See Fig. 3 (a)), the experimental curves were first digitized using the image digitization tool in Origin software to extract characteristic data points, which were then normalized. The normalized data were then fitted using a Prony series to determine the model parameters. Finally, the obtained parameters were imported into the finite element system. This data-driven modeling approach effectively reduces the cost associated with complex experiments and provides an efficient numerical simulation pathway for assessing the long-term performance of timber structures.

Fig. 3. Generalized Maxwell model (a), Stress relaxation curve of Pinus sylvestris var. mongolica (b), and Nonlinear visco-elastic model (c)

In Fig. 3, E₀, E₁, E₂, etc. are springs, η₁, η₂, η₃, etc. are viscous dampers, E’₁, E’₂, etc. are hyperelastic, and η₀ is a viscous damper.

 (4)

In Eq. 4, e and ϕ are the mechanical bias strain and the bulk strain, respectively; G is the shear modulus, and K is the bulk modulus, which are reduced time τ.

Linear viscoelastic models effectively predict the relaxation behavior of wood under small deformations but fail to account for the influence of strain level. Wood exhibits nonlinear stress relaxation under large deformations, a limitation stemming from the linear dashpot in classical Generalized Maxwell models. Consequently, the theoretical framework based on the Boltzmann superposition principle is constrained in this regime. To address this limitation, a numerical approach incorporating hysteresis is employed to integrate nonlinear viscoelasticity into a computational approach, thereby enhancing the simulation accuracy of wood’s mechanical response under large deformations. The key model improvement involves replacing the linear spring within the Maxwell element with a hyperelastic body and substituting the linear dashpot with a nonlinear viscoelastic element. While normalized relaxation curves for linear viscoelasticity are consistent, their nonlinear counterparts demonstrate strain dependency. The modeling framework utilizes two parallel branches (A and B) sharing the deformation gradient. Branch A is purely hyperelastic (supporting various models such as Neo-Hookean). Branch A consists of Eq.5, while Branch B comprises Eqs. 6 through 8). Crucially, the hyperelastic models in both branches differ only by a scaling factor, K. Specifically, this work combines a Neo-Hookean model with hysteretic effects to simulate the nonlinear viscoelasticity of wood. The model captures path dependency and irreversible deformations. Its nonlinear nature is intrinsically reflected in the nonlinear relationship between equivalent stress and equivalent deformation (as shown in the Eq. 8),

 (5)

where J_A is the Jacobi determinant under branch A, ψ is the strain energy function,  is the right Cauchy isolation strain without volume deformation,  is the left Cauchy isolation strain without volume deformation. The material constants C₁₀ and D₁ are modulus-related parameters. In Eqs. 6 and 7,  is the velocity gradient due to visco-elasticity under branch B,  is the velocity gradient due to super-elasticity under branch B, and  is the gradient of deformation due to super-elasticity under branch B.

 (6)

 (7)

In Eq. 8,  is the strain rate, and σ_B is the equivalent stress.

 (8)

This study employed four-point bending simulations to assess the capability of a timber member model in characterizing the viscoelastic properties (creep and stress relaxation) of wood. Within the simulation approach, uniform loading was applied to wood elements at multiple rates until reaching specified displacements to investigate the relationship between modeled stress and loading rate. Analysis of the subsequent behavior during the displacement-hold phase was performed, with a focus on stress relaxation. Curve-fitting methodology was applied to quantitatively analyze the strain-rate effects (loading rate dependency) and time-dependent characteristics (specifically stress relaxation behavior) within the simulation results. Based on this quantification process, this study evaluated the performance of the viscoelastic numerical model for wood. The Neo-Hookean hyperelastic parameters in the nonlinear viscoelastic hysteresis model were calibrated using the longitudinal tensile stress-strain curve of wood, while the viscoelastic part was fitted through the stress relaxation curve (Li et al. 2019; Muliana et al. 2016). According to Eq. 5 and in line with the requirements of finite element simulation, the representative key parameters of the model were set as follows: C₁₀= 250, D₁ = 0.03; the creep parameter was 5E-08, the effective stress exponent was 2, and the creep strain exponent was taken as 0.3 (Wu 2021; Luo et al. 2021). Cyclic loading-unloading processes were simulated using a defined amplitude profile (0→1→0→2→0), with the resulting force-displacement responses at loading points analyzed and compared against a linear viscoelastic model. This integrated simulation approach systematically evaluates the model’s fidelity in replicating wood’s viscoelastic response, ultimately establishing a three-stage analysis Approach: small-deformation elastic verification → elastoplastic damage evolution analysis → viscoelastic time-dependent life prediction.

Four-point Bending Test

This study investigated the mechanical properties and failure mechanisms of Pinus sylvestris var. mongolica specimens using four-point bending tests. The specimens measured 15 × 20 × 240 mm³ and were made of wood sourced from the Greater Khingan Range in Heilongjiang Province. Under air-dried conditions, the average moisture content was approximately 12%, and the air-dry density (ρ₀) was 480 kg/m³. All specimens were prepared and conditioned in accordance with the ASTM D143 (2014) standard. The test setup employed a support span of 200 mm, with loading spans of 50 mm and 150 mm, respectively. A loading rate of 5 mm/min was applied to simulate quasi-static loading conditions until the specimens underwent progressive delamination and fragmentation failure. During testing, load and displacement data were synchronously recorded using the testing machine’s built-in load cell and an LVDT displacement sensor. Throughout the loading process, full-field principal strains (ε₁₁, ε₂₂), shear strain (γ₁₂), and displacement gradients were captured using three-dimensional Digital Image Correlation (3D-DIC) technology. The displacement field was calculated with a subset size of 29 pixels × 29 pixels and a step size of 10 pixels. Subsequently, the strain field was derived via a strain window based on least-squares fitting. The obtained stress-strain curves and DIC strain fields provide a high-fidelity data basis for the nonlinear mechanical simulation of timber structures.

RESULTS AND DISCUSSION

Simulation Results of the Elastic Model

Numerical simulations of four-point bending tests on 15 × 20 × 200 mm³ wood specimens revealed fundamentally different mechanical responses between orthotropic and isotropic material models (Table. 3). The orthotropic model demonstrated a substantially lower maximum von Mises stress, reaching only 44% of the isotropic value at maximum deflection. This stress reduction originates from wood’s inherent anisotropic stiffness distribution: the orthotropic model accurately reflects the fact that the stiffness along the longitudinal (fiber) direction is substantially higher than that in transverse directions. Under bending loads, the higher longitudinal stiffness promotes more efficient load transfer along the grain, thereby reducing transverse stress concentrations and optimizing the overall stress distribution. In contrast, the isotropic model ignores this stiffness difference, overestimating transverse stresses and leading to a systematic overestimation of von Mises stress.

Table 3. Simulation Results of the Linear Elastic Model

Both tensile and compressive peaks of longitudinal normal stress (S11) showed substantial reductions of 55.09% and 56.34%, respectively, under orthotropic conditions, highlighting the stress-relieving effect of higher stiffness along the fiber direction. Similarly, the compressive peak of transverse normal stress (S22) decreased from -21.35 to -10.12 MPa (52.6% reduction) under orthotropy, accurately reflecting wood’s inherently weaker transverse mechanical properties. These findings collectively validate that the orthotropic model successfully captured the direction-dependent stiffness (E₁>E₂) and characterized the “strong-longitudinal, weak-transverse” behavior of wood. This result provides a crucial and reliable analytical approach for rapidly identifying actual high-risk areas in timber components during early service stages, effectively preventing stress misjudgments arising from the use of isotropic models.

Simulation Results of the Elasto-plasticity Model

A comparative elastoplastic analysis of a timber component under four-point bending, utilizing both the generalized Hill anisotropic and Mises isotropic yield criteria (Table 4), revealed distinctly different strain evolution patterns despite similar strain localization near loading and support points. The Mises criterion produced a smoothly graded, symmetric strain distribution decaying from the pure bending region toward the neutral axis. In contrast, the Hill criterion resulted in a multi-polar, discontinuous strain pattern, with high-strain zones concentrated at support/loading points that subsequently propagate inward to form secondary yield regions.

Table 4. Simulation Results of the Elastoplastic Model

Notably, the Hill criterion generated a more uniform in-plane principal plastic strain (PE) distribution overall, with higher plastic strains at loading/support points but lower values at the maximum deflection location compared to the Mises model. This distribution discrepancy arises from wood’s anisotropic regulation mechanism: strain amplifies in weakened material directions under multiaxial stress, while constrained deformation occurs in strengthened orientations due to transverse effects in principal stress-dominated zones. These findings validate the decisive role of anisotropic constitutive modeling in simulating structural nonlinearity. Therefore, the Hill anisotropic criterion is recommended for accurately predicting failure modes and identifying vulnerable regions in timber components, as its predictions fundamentally determine the effectiveness of reinforcement strategies and the resulting structural safety margin.

Simulation Results of the Visco-elastic Model

Figure 4 (a) clearly demonstrates the clear sensitivity of stress to loading rate in the viscoelastic wood element. As the loading rate (descent rate) increases, the viscoelastic stress exhibits a nonlinear escalation(precisely fitted by the polynomial equation). This results in a stress difference of 7.25% between high and low rates (as given by the equation in Fig. 4). This phenomenon directly evidences the enhancement of the viscoelastic energy dissipation mechanism under high strain rates. Intermolecular frictional resistance leads to increased energy absorption capacity during dynamic loading. This rate dependence, absent in the control group of purely elastic models, highlights an inherent characteristic of wood. Complementing this mechanical behavior, Fig. 4(b) illustrates the time-dependent attenuation effect under constant load. The principal stress undergoes an exponential decay of approximately 5.9% over the course of 1000 time increments. This attenuation stems from the gradual release of internally stored elastic strain energy-driven microscopically by molecular chain reorganization (e.g., slippage of cellulose microfibrils) and phase adjustment and manifested macroscopically as characteristic stress relaxation.

Fig. 4. Comprehensive comparative analysis of linear elastic wood column models

Beyond the aforementioned intrinsic viscoelastic characteristics of the material, such as rate sensitivity and stress relaxation, the comprehensive mechanical behavior of timber components under actual service or testing conditions is considerably more complex. These behaviors are profoundly influenced not only by boundary conditions and large deformations but also serve as a key criterion for distinguishing between the responses of linear and nonlinear models. To investigate these aspects, this study employed four-point bending simulations, with a focused analysis on the cyclic loading response of the components. Four-point bending simulations revealed that contact behavior between the specimen and supports/indenters (e.g., friction or micro-separation/re-contact) introduced nonlinearities. This can cause the unloading path to deviate from the ideal state, leading to the formation of hysteresis loops. Even for linearly elastic materials, this contact-induced nonlinearity may generate “pseudo-hysteretic loops,” necessitating model calibration or the use of uniaxial compression experiments in order to obtain accurate hysteretic properties. For pure linear elastic materials, loading/unloading force-displacement curves coincide perfectly, exhibiting no rate-dependence. In contrast, linearly viscoelastic wood models exhibit elastic hysteresis due to time-dependency (Fig. 5), manifesting as separated loading/unloading paths that form characteristic hysteresis loops and the presence of residual displacement, reflecting irreversible energy dissipation. Considering the fact that large deformations in nonlinear viscoelasticity further complicates the response: the hysteresis loop shape becomes sensitive to loading rate/amplitude, and modulus exhibits nonlinear decay at maximum compression, revealing more complex energy dissipation mechanisms. Deeper understanding requires single/double-amplitude cyclic simulations: double-amplitude loops exhibit larger enclosed areas, signifying enhanced energy dissipation. Subsequent loading paths in single-amplitude cycles deviate from the initial curve due to historical viscous effects, demonstrating noticeable material memory of prior deformation. These characteristics confer unique engineering advantages upon nonlinear viscoelastic materials: their high energy dissipation capacity enables effective energy absorption and vibration damping, while stress relaxation properties are crucial for the design of long-term load-bearing structures. These viscoelastic characteristics provide unique engineering advantages, with the model offering critical full life-cycle analysis capability for timber components. In the short term, the rate-dependent behavior enables quantitative assessment of instantaneous energy dissipation and dynamic response under seismic or wind loads. For long-term performance, the stress relaxation model accurately predicts creep deformation and prestress loss in critical members such as beams and columns, thereby providing indispensable quantitative basis for evaluating long-term safety and service life.

Fig. 5. Force-displacement curves for linear and nonlinear viscoelastic loading-unloading cycles

Results of the Four-point Bending Test

This study investigated the strain field evolution and constitutive characterization strategies of 15 × 20 × 200 mm³ wooden elements under four-point bending using an integrated Digital Image Correlation (DIC) and Finite Element Method (FEM) approach (Fig. 6). The FEM simulations, based on an anisotropic material model, demonstrated high consistency with DIC measurements in transverse strain (ε_xx) distribution.

Fig. 6. DIC experimental and FEA simulation comparison

Fig. 7. A schematic of the multi-model collaborative finite element approach guided by service-condition matching

Both methods distinctly captured the mechanical characteristics of the neutral axis: a compressive strain gradient above transitioning to a tensile gradient below, with noticeable agreement on the spatial distribution of strain extrema. Longitudinal strain (ε_yy) analysis revealed synchronized full-field evolution patterns between DIC and FEM, including localized strain effects at loading points, linked to stress concentration and large deformations associated with wood’s viscoelastoplastic energy dissipation mechanism. For shear strain (γ_xy), both DIC and FEM successfully characterized typical symmetrically diffusing shear concentration zones near loading points, exhibiting high congruence in gradient evolution patterns. Notably, DIC enabled continuous recording of time-dependent strain fields under sustained load, providing crucial experimental data for creep characterization and supporting multi-physical field coupling parameter identification, though requiring consideration of long-term observation needs inherent to wood’s time-dependent behavior. A novel “service-state matching” principle for constitutive model selection was proposed: linear elasticity for rapid assessment, elastoplasticity for ultimate load capacity analysis, and viscoelasticity for long-term performance prediction. This underpins a three-tiered collaborative analysis approach: 1) initial screening of critical zones via small-deformation elastic analysis, 2) ultimate state analysis using an elastoplastic damage model, and 3) residual life prediction employing a viscoelastic model. This forms a complete, cross-validating methodology for timber structure assessment. Applied specifically to historical buildings, this integrated chain enables damage zone identification, quantification of ultimate load capacity, and long-term service performance evaluation.

CONCLUSIONS

1. The multi-model collaborative finite element approach, conceptually summarized in Fig. 7, establishes a dynamically adaptive simulation framework for timber structural analysis. Its core innovation lies in the service-condition-matched selection and integration of constitutive models, enabling high-fidelity simulation across the component’s entire life cycle—from instantaneous elastic response to long-term time-dependent behavior.
2. This approach systematically constructs a three-tiered collaborative analysis system: “small-deformation elastic verification → elastoplastic damage analysis → viscoelastic time-dependent life prediction.” It dynamically employs an orthotropic linear elastic model for stress and deformation assessment under service loads, it switches to an elastoplastic model based on the Hill anisotropic yield criterion to capture damage and ultimate capacity near yield limits, and it introduces a generalized Maxwell or nonlinear viscoelastic model to predict creep, stress relaxation, and energy dissipation under sustained loading.
3. The “condition-model” matching mechanism fundamentally overcomes the limitations of traditional single-constitutive models. It provides a systematic, quantitative, and cross-scale numerical tool for the safety diagnosis, damage localization, residual capacity assessment, and long-term performance prediction of timber structures.
4. This methodology offers crucial technical support for the scientific conservation and performance optimization of historic timber members, facilitating a paradigm shift in cultural heritage preservation from empirical practice to quantitative analysis. Furthermore, the research establishes a theoretical and numerical foundation for enhancing the safety and durability of modern timber structures, deepening the understanding of wood’s multi-scale mechanical behavior and providing a core simulation platform for whole-life-cycle performance design and safety management. Furthermore, this study explicitly defines quantitative switching criteria between models. The criteria are used to dynamically select constitutive models during the simulation process. This approach overcomes the arbitrariness of traditional single‑model methods and provides a transparent, reproducible numerical workflow for analyzing multi‑stage behavior.
5. The multi-model collaborative framework developed in this study provides an effective tool for the high-fidelity numerical simulation of timber components. Looking ahead, this framework can be further extended: on the one hand, by systematically incorporating natural defects such as knots, cracks, and voids into the model; on the other hand, by integrating it with 3D reconstruction technology for the analysis of larger-scale, actual damaged timber components. This will contribute to a more comprehensive and realistic assessment of the safety status and long-term durability of practical timber structures.

ACKNOWLEDGMENTS

The authors gratefully acknowledge the support from the National Natural Science Foundation of China, Grant No. 32171700.

Conflict of Interest

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

REFERENCES CITED

Berg, S., Sandberg, D., Ekevad, M., and Vaziri, M. (2015). “Crack influence on load-bearing capacity of glued laminated timber using extended finite element modelling,” Wood Material Science & Engineering 10(4), 335-343. https://doi.org/10.1080/17480272.2015.1027268

Bisong, S., Lepov, V. V., and Etinge, A. R. (2024). “Study of the mechanical behavior and multiscale simulation of the crack propagation in a bilinga wooden beam,” Zavod. Lab. 90(3), 52-61. https://doi.org/10.26896/1028-6861-2024-90-3-52-61

Chen, L. T. (2017). Study on the Degradation of Mechanical Properties of Damaged Wooden Beams, Ph.D. Dissertation, Xi’an University of Architecture and Technology, Xi’an, China.

Dassault Systèmes (2016). Abaqus Theory Guide, Dassault Systèmes, Providence, RI, USA. https://ceae-server.colorado.edu/v2016/books/stm/default.htm

Fortino, S., Hradil, P., Genoese, A., Genoese, A., and Pousette, A. (2019). “Numerical hygro-thermal analysis of coated wooden bridge members exposed to Northern European climates,” Construction and Building Materials 208, 492-505. https://doi.org/10.1016/j.conbuildmat.2019.03.012

GB 55005-2021 (2021). “General code for timber structures,” China Architecture & Building Press, Beijing, China.

GB/T 1927.9-2021 (2021). “Test methods for physical and mechanical properties of small clear wood specimens—Part 9: Determination of bending strength,” Standardization Administration of China, Beijing, China.

Gomon, P., Gomon, S., Pavluk, A., and Homon, S. (2023). “Innovative method for calculating deflections of wooden beams based on the moment-curvature graph,” Procedia Structural Integrity 48, 195-200. https://doi.org/10.1016/j.prostr.2023.07.148

Hassan, M. A., Usman, M., Hanif, A., Farooq, S. H., and Ahmed, J. (2020). “Improving structural performance of timber wall panels by inexpensive FRP retrofitting techniques,” Journal of Building Engineering 27, article 101004. https://doi.org/10.1016/j.jobe.2019.101004

He, S. (2011). Finite Element Model Analysis of Mongolian Scotch Pine Finger-Jointed Timber, Dissertation, Nanjing Forestry University, Nanjing, China.

Karagöz İşleyen, Ü., and Kesik, H. İ. (2021). “Experimental and numerical analysis of compression and bending strength of old wood reinforced with CFRP strips,” Structures 33, 259-271. https://doi.org/10.1016/j.istruc.2021.04.070

Li, Z., Wang, Y. L., Li, X., Li, Z. R., and Wang, Y. (2019). “Experimental investigation and constitutive modeling of uncured carbon black filled rubber at different strain rates,” Polymer Testing 75, 117-126. https://doi.org/10.1016/j.polymertesting.2019.02.005

Li, Y. H. (2023). Study on Static Characteristics of Typical Chinese Dou-Gong Bracket, Dissertation, Inner Mongolia Agricultural University, Hohhot, China. https://doi.org/10.27229/d.cnki.gnmnu.2023.000035

Liu, S. D. (2011). Study on the Drying Characteristics of Russian Larch Timber and Its Square Timber Drying Technology, Ph.D. Dissertation, Inner Mongolia Agricultural University, Hohhot, China.

Liu, W. Q., and Yang, H. F. (2019). “Research progress on modern timber structures,” Journal of Building Structures 40(2), 16-43. https://doi.org/10.14006/j.jzjgxb.2019.02.002

Luo, Y. H., Zhang, R. B., Liang, J. X., Lu, C. X., and Diao, H. L. (2021). “Study on the creep characteristics of Eucalyptus dunnii wood,” Hubei Agricultural Sciences 60(7), 80-83. https://doi.org/10.14088/j.cnki.issn0439-8114.2021.07.013

Ma, W. L. (2014). The Experimental Research of Dynamic Loading the Mongolian Scotch Rift Grain and Cross Grain, Dissertation, Northeast Forestry University, Harbin, China.

Mansor, M., Mohareb, M., and Doudak, G. (2025). “Elastic lateral torsional buckling of two-ply built-up wooden beams connected with discrete fasteners,” Engineering Structures 331, article 119944. https://doi.org/10.1016/j.engstruct.2025.119944

Muliana, A., Rajagopal, K. R., Tscharnuter, D., and Pinter, G. (2016). “A nonlinear viscoelastic constitutive model for polymeric solids based on multiple natural configuration theory,” International Journal of Solids and Structures 100-101, 95-110. https://doi.org/10.1016/j.ijsolstr.2016.07.017

Pan, G. X. (2003). A History of Chinese Architecture, 5^th Ed., China Architecture and Building Press, Beijing, China.

Pan, Y., Li, T. Y., Yang, Q. S., Shi, X. W., Meng, X. J., and Chen, J. Y. (2024). “Numerical modelling strategy and load resisting mechanism of palace-type traditional wooden frame,” Engineering Structures 319, article 118851. https://doi.org/10.1016/j.engstruct.2024.118851

Rong, X. W. (2008). Wood Elasticity and Viscoelasticity, Dissertation, Anhui Agricultural University, Hefei, China. https://doi.org/10.7666/d.y1394034

Shi, B. X. (2023). Research on Degradation Mechanism of Mechanical Properties of Damaged Timber Columns and Composite Retrofitting Techniques, Dissertation, Southwest University of Science and Technology, Mianyang, China. https://doi.org/10.27415/d.cnki.gxngc.2023.000069

Wang, A. S. (2023). Research on Mechanical Behavior and Simulation of Damaged Wood Element with Oblique Crack, Ph.D. Dissertation, Liaocheng University, Liaocheng, China. https://doi.org/10.27214/d.cnki.glcsu.2023.000991

Wang, S., Lin, J. K., Dong, S. W., Chen, Z. Y., Kong, F. X., Ma, P. P., Wang, F. B., and Que, Z. L. (2025). “Lateral resistance performance of wood-frame shear walls with wooden nail connections: Experimental and finite element analysis,” Journal of Bioresources and Bioproducts 10(4), 410-424. https://doi.org/10.1016/j.jobab.2025.04.002

Wen, Y. X., and Zhao, D. (2016). “Determination of neutral axis position in wooden beams with holes using digital speckle correlation method,” in: Proceedings of the 16th Academic Conference of the Northern Seven Provinces and Autonomous Regions Mechanics Society, Chinese Society of Mechanics, Beijing, China, pp. 218-220.

Wu, Q. Y. (2021). Study on the Duration of Load Effect in Wood, Dissertation, Harbin Institute of Technology, Harbin, China. https://doi.org/10.27061/d.cnki.ghgdu.2021.005288

Xue, J., Ren, G., Qi, L., and Zhang, X. (2021). “Lateral behavior of glued-laminated timber frame infilled with light-wooden-frame wall hybrid system: Experimental and numerical analysis,” Structures 30, 352-367. https://doi.org/10.1016/j.istruc.2021.01.022

Yu, H., Chen, B., Huang, Y., and Zhou, L. (2024). “Experimental and numerical investigation on collapse performance of wooden drum-towers in southwest China,” Engineering Failure Analysis 159, article 108098. https://doi.org/10.1016/j.engfailanal.2024.108098

Yu, Y. Z. (2022). Nondestructive Detection of Defects in Wooden Columns of Ancient Building Walls and Safety Numerical Simulation Study, Dissertation, Beijing Forestry University, Beijing, China. https://doi.org/10.26949/d.cnki.gblyu.2022.000523

Zhang, L. (2017). Experimental Study and Numerical Simulation Methods for Mechanical Properties of Wood, Dissertation, Beijing Jiaotong University, Beijing, China.

Zhang, S. (2022). Analysis of Bearing Capacity and Dynamic Response of Wood Structure Components, Ph.D. Dissertation, Inner Mongolia Agricultural University, Hohhot, China. https://doi.org/10.27229/d.cnki.gnmnu.2022.001346

Zhou, Q., and Zhang, X. Q. (2007). “Analysis of the current status of the eaves purlin structure of the Xishan pavilion in the hall of supreme harmony, Forbidden City,” in: Proceedings of the Chinese Forbidden City Society, Palace Museum, Beijing, China, pp. 692-702.

Article submitted: November 12, 2025; Peer review completed: December 13, 2025; Revised version received and accepted: December 21, 2025; Published: January 2, 2026.

DOI: 10.15376/biores.21.1.1497-1514