# Predicting rolling shear failure of cross-laminated panels using a multiscale domain decomposition approach

Jorge A. Fernández González, Karin Saavedra, O. Allix, Pierre Gosselet

► To cite this version:

Jorge A. Fernández González, Karin Saavedra, O. Allix, Pierre Gosselet. Predicting rolling shear failure of cross-laminated panels using a multiscale domain decomposition approach. Advanced Modelling and Simulation in Engineering Sciences, 2026, 13 (1), pp.1. 10.1186/s40323-026-00318-9. hal-05519912

HAL Id: hal-05519912

https://hal.science/hal-05519912v1

Submitted on 20 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
RESEARCH

Open Access

<figure><img src="image_2.png" /></figure>

Check for
updates

# Predicting rolling shear failure of cross-laminated panels using a multiscale domain decomposition approach

J. Fernández<fnref n="1" />,<fnref n="2" />,<fnref n="5" />, K. Saavedra<fnref n="3" />*†, O. Allix<fnref n="2" />† and P. Gosselet<fnref n="2" />,<fnref n="4" />†

## Abstract

Cross-laminated timber (CLT) is gaining popularity as a sustainable building material due to its favorable structural performance and low environmental impact. This study introduces a multiscale domain decomposition method, known as the LaTIn strategy, to predict the out-of-plane behavior and failure of CLT panels. The proposed approach incorporates cohesive zone models to capture inter- and intralaminar damage using a cohesive parameter identification based on strength values reported in the literature. By aligning domain partitioning with principal stress directions and leveraging parallel computations, the method reduces degrees of freedom and the number of subdomains by up to 55% and 70%, respectively. Additional algorithmic enhancements allow for a 70% reduction of iterations. Validation with experimental data confirms the model's ability to predict rolling shear failure with high accuracy. These results demonstrate the potential of the strategy for efficient and reliable simulation of timber structures, in line with the needs of multiscale nonlinear analysis in structural engineering.

**Keywords:** Cross-laminated timber, Rolling shear strength, Domain decomposition, Cohesive zone model

## Introduction

The use of engineered wood products is gaining momentum in the construction sector due to the increasing demand for sustainable resources and the advantageous mechanical properties of these materials, such as the high strength-to-weight ratio, stiffness and suitability for prefabrication [1–3]. Among them, Cross-Laminated Timber (CLT), in particular, has emerged as a promising solution for mid- and high-rise buildings, with successful applications in countries such as Canada, Norway, France, and the United Kingdom.

CLT panels consist of multiple layers of lumber, each ply is typically oriented perpendicularly to adjacent layers and bonded with structural adhesives to form solid panels [2], typically in odd-numbered layers (three, five or seven). Despite extensive standardization efforts [4–7], several challenges remain, one of the most critical being the evaluation and prediction of rolling shear failure (RS) [8], largely due to the orthotropic nature of wood.

<footnote>†K. Saavedra, O. Allix and P. Gosselet contributed equally to this work.</footnote>

<footnote>*Correspondence: ksaavedra@utalca.cl</footnote>

<footnote n="1">¹Engineering Systems Doctoral Program, Faculty of Engineering, Universidad de Talca, Curicó, Chile</footnote>

<footnote n="2">²LMPS-Laboratoire de Mécanique Paris-Saclay, Université Paris-Saclay, Gif-sur-Yvette, France</footnote>

<footnote n="3">³Depto. de Tecnologías Industriales, Universidad de Talca, Campus Curicó, Curicó, Chile</footnote>

<footnote n="4">⁴Univ. Lille, CNRS, Centrale Lille, UMR 9013-LaMcube-Laboratoire de Mécanique, Multiphysique, Multiéchelle, 59000 Lille, France</footnote>

<footnote n="5">⁵Depto. de Ingeniería Mecánica, Universidad del Bío-Bío, Concepción, Chile</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Local and global coordinate systems and dimensions of a CLT panel</figcaption></figure>

As shown in Fig. 1, wood exhibits different elastic properties along its longitudinal (L), radial (R), and tangential (T) directions. The low shear strength in the transverse directions (R and T) leads to critical stress states in the cross-layers under out-of-plane loads, where RS stress can induce twisting and failure of the fibers (see Fig. 2) [9, 10]. RS failure has been widely studied using experimental methods such as the four-point short-span bending test [11, 12], although improvements have been proposed [13].

Experimental studies have provided valuable insights into RS failure, but the complex interaction of failure modes, especially under nonuniform stress states, highlights the need for advanced computational tools. Numerical models can complement experimental investigations and help reduce the cost and duration of physical testing, while also allowing the analysis of loading scenarios that are difficult to replicate experimentally [14, 15]. However, the predictive capacity of these models depends on rigorous validation, particularly when multiple fracture mechanisms are involved [16–18].

This study proposes a computational approach to predict RS strength through non-linear finite element analyses. A mixed multiscale domain decomposition method is introduced, integrating the LaTIn strategy [19] with an efficient partitioning approach specifically designed to capture RS and its interaction with other failure modes. This approach has been successfully applied in the simulation of laminated composites [20, 21] and other multiphysics problems [22]. To capture cracks due to RS failure, the cohesive zone model (CZM) is used at the parallel partitioning interfaces [23, 24]. In this work, simulations are focused on 3-layer panels and are compared with the experimental data available in the literature [11]. In summary, three key aspects are addressed here:

* a stress field analysis of CLT panels under out-of-plane loads, identifying the principal stresses and their orientations, in order to provide information to take into account to predict RS failure;

* the implementation of a mixed domain decomposition method incorporating cohesive zone models to capture RS failure. To reduce computational costs, it includes:
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Rolling shear failure observed in a 3-layer CLT panel</figcaption></figure>

(i) an adapted partitioning based on the preferred crack patterns and the size of the fracture process zone (FPZ), and (ii) algorithm modifications for initializing variables;

* the cohesive zone model parameters are identified from the wood's strength properties to capture the analysed failure mechanisms. The effect of the board's rolling shear resistance on the panel's overall response is also assessed.

The paper is organized as follows. "Mechanical modeling of CLT panels" section discusses the mechanical behavior of CLT panels and introduces the approach used to model failure. "Finite element analysis for predicting rolling shear failure" section details the mesomodel, the reference problem under study, and the domain decomposition strategy. "Assessment and optimization of the computational approach" section follows with three case studies, including the evaluation of RS failure in a CLT panel. Finally, "Conclusions" section summarizes the main conclusions of this study.

## Mechanical modeling of CLT panels

Wood can be characterized at the mesoscopic level as an orthotropic material in the local coordinate system (L, R, T). However, when considering a coniferous timber defined by EN 338 (see moduli of elasticity and strengths in Table 1), a transverse isotropy material can be assumed [25], while the global coordinates ($x_i, i \in \{1, 2, 3\}$) are adopted. In this case, the elastic moduli are expressed as fractions of the longitudinal elastic modulus $E_L$: $E_1 = E_L, E_2 = E_3 = E_L/30, G_{12} = E_L/16$, and the rolling shear modulus $G_{23} = E_L/160$ [4, 26].
<table><caption>**Table 1** Characteristic elastic properties and strengths for timber (MPa) [27]</caption><thead><tr><th>Class</th><th>E<sub>L</sub></th><th>E<sub>R</sub></th><th>E<sub>T</sub></th><th>G<sub>LT</sub></th><th>G<sub>LR</sub></th><th>f<sub>m</sub></th><th>f<sub>v</sub></th><th>f<sub>Lt</sub></th><th>f<sub>Lc</sub></th><th>f<sub>Rt</sub></th><th>f<sub>RC</sub></th><th>f<sub>Rt</sub></th><th>f<sub>Tc</sub></th></tr></thead><tbody><tr><td>C16</td><td>8000</td><td>270</td><td>270</td><td>500</td><td>500</td><td>16</td><td>1.8</td><td>10</td><td>17</td><td>0.5</td><td>2.2</td><td>0.5</td><td>2.2</td></tr><tr><td>C18</td><td>9000</td><td>300</td><td>300</td><td>560</td><td>560</td><td>18</td><td>2.0</td><td>11</td><td>18</td><td>0.5</td><td>2.2</td><td>0.5</td><td>2.2</td></tr><tr><td>C20</td><td>9500</td><td>320</td><td>320</td><td>600</td><td>600</td><td>20</td><td>2.2</td><td>12</td><td>19</td><td>0.5</td><td>2.3</td><td>0.5</td><td>2.3</td></tr><tr><td>C24</td><td>11,000</td><td>370</td><td>370</td><td>690</td><td>690</td><td>24</td><td>2.5</td><td>14</td><td>21</td><td>0.5</td><td>2.5</td><td>0.5</td><td>2.5</td></tr></tbody></table>

m, bending; v, shear; t, tension; c, compression

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Four-point bending test setup according to EN 16351:2021</figcaption></figure>

<table><caption>**Table 2** Geometrical parameters as a function of $h_T$ for a 4-point bending test according to EN 16351:2021</caption><thead><tr><th>Description</th><th>Variable</th><th>Value</th></tr></thead><tbody><tr><td>Individual thickness</td><td>h</td><td>–</td></tr><tr><td>Total thickness</td><td>h<sub>T</sub></td><td>Σ<sup>n</sup> h<sub>i</sub></td></tr><tr><td>Span</td><td>L</td><td>12 h<sub>T</sub></td></tr><tr><td>Overlap of support</td><td>a<sub>0</sub></td><td>0.5 h<sub>T</sub></td></tr><tr><td>Dist. support/load</td><td>a<sub>1</sub></td><td>3 h<sub>T</sub></td></tr><tr><td>Dist. load/load</td><td>a<sub>2</sub></td><td>6 h<sub>T</sub></td></tr><tr><td>Width</td><td>b</td><td>≥ 600 mm</td></tr></tbody></table>

This Section examines an analytical method to determine the stress distribution through the thickness of a panel under four-point bending. The idea is to identify the regions that are more prone to failure. Then, the approach to predict rolling shear failure is introduced.

## Stress analysis

Analytical methods are used to evaluate the equivalent bending and shear stiffness of CLT panels, enabling the observation of stress distributions across the thickness of the panel. Common approaches include the Shear Analogy Method (SAM) introduced by Kreuzinger [28], the γ-Method, and the modified Timoshenko method [29]. In this study, the SAM following Mestek's methodology [26,30] is applied to a four-point short-span beam, a standard configuration to evaluate bending strength. This methodology provides a more accurate prediction of the bending properties [25]. Panels with 3, 5, and 7 layers, made from the same wood grade, same layer width ($b_i = b$) and same layer thickness ($h_i = h$), are considered. Fig. 3 shows an example of such a test setup, and Table 2 provides their characteristic dimensions.

Applying SAM-Mestek's methodology to the four-point short-span beam, the stress distribution across the thickness can be determined, especially in the transverse layers, which have a lower shear strength. Fig. 4 illustrates the principal stress distribution through the
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Four-point bending test: principal stresses on the weakest layer for panels with $n \in \{3, 5, 7\}$ layers</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 5** Four-point bending test: experimental results for 3- and 5-layer panels [10]</figcaption></figure>

thickness in the layer $i = 2$ for panels with $n = 3$ and $n = 5$ layers and in the layer $i = 4$ for panels with $n = 7$ layers. It is observed that the principal stresses are predominantly influenced by rolling shear stress. The principal angle is close to 45° for panels with 3 and 7 layers, and slightly higher, approximately 52°, for 5 layers. This suggests that RS cracks are likely to propagate along these angles, which aligns with experimental observations [31] (see more details in Fig. 5a and b). This insight forms the basis for the numerical model proposed in this work, as it highlights the need to accurately capture the trajectories of potential cracks and to perform efficient computational simulations as well.

## Predicting rolling shear failure

Traditional approaches to evaluating the structural integrity of CLT panels often rely on strength failure criteria, comparing the stress state in areas of interest with characteristic values of wood strength [32,33]. Although these methods can be effective for ideal, defect-free materials, they can present drawbacks when dealing with areas containing stress concentrations or weak zones. For RS failure, several empirical equations have been developed to determine the RS strength $f_r$, and to improve current standards [11, 13, 34]. In particular, Mestek et al. [34] verified that the interaction between compressive and shear stresses enhances the RS strength. Although these criteria help predict the onset of failure, they do not provide insight into the post-failure behavior or the progression of damage. RS failure typically occurs in the central layer of CLT panels without a pre-existing crack, making traditional fracture mechanics approaches less suitable. Instead, a damage-based approach, in particular cohesive zone models, is more appropriate to model crack initial-
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Degradation of interlaminar and intralaminar interface</figcaption></figure>

tion and propagation. In this study, a damage mechanics approach is adopted using CZMs to model the deterioration of wood properties due to RS failure, considering that it occurs in a brittle but progressive manner. This approach allows for the evaluation of the load at which failure occurs and the structural response thereafter. CZMs have been successfully implemented in finite element frameworks to address delamination problems in composite, orthotropic, and brittle materials [35,36]. CZMs allow for a stress-displacement degradation relationship that simulates gradual deterioration of the material properties. This has already been used to model wood cracking under different boundary conditions [37,38], and specifically for RS failure [21]. In this work, the constitutive cohesive law proposed by Allix et al. [39] is employed (see details in “Mesomodel” section).

Finite element analysis for predicting rolling shear failure

A finite element approach for predicting RS failure is proposed in this Section, which is based on a mesomodel and a domain decomposition method.

Mesomodel

The mesomodel developed by Ladevèze and Allix [23,24] is used here, which describes the whole at the mesoscale using volumetric entities (layers) and 2D interfaces, each with their own behavior. This mesomodel has initially been formulated using only interlaminar interfaces to capture delamination, as shown in Fig. 6 (left). However, this work proposes the introduction of intra-laminar interfaces to cover other failure modes typically observed in CLT panels, such as rolling shear failure, transverse cracking, and normal failure due to bending stress, as illustrated in Fig. 6 (right).

The behavior of the interface must be defined on a local orthotropic basis ($e_1, e_2, e_3$), where $e_3$ is normal to the plane and ($e_1$ and $e_2$) lie in the plane of the interface. Therefore, the displacements and stresses of two adjacent laminates can be related through a CZM. The surface energy of the damaged interface can be defined on an orthotropic basis using Eq. (1). The damage evolution indicator $d$ is governed by thermodynamic forces $Y_i = -\partial e_d / \partial d_i$, representing the energy dissipated due to interface damage. The irreversibility of the process is represented by the relationship expressed in Eq. (2), while the evolution law proposed by Allix et al. [39] is provided in Eq.(3) (see Fig. 7). The surface energy of each mode $Y_i$ is calculated based on the relationship between the stiffness of the undamaged
<figure><img src="image_9.png" /><figcaption>Fig. 7 Potential cohesive zone model (normal direction)</figcaption></figure>

interface ($k_i^0$) and the displacement jump $[[u_i]]$: $Y_i = \frac{1}{2} k_i^0 [[u_i]]^2$.

$$e_d = \frac{1}{2} \left( \frac{h^+(-\sigma_{33})^2}{k_3^0} + \frac{h^+(\sigma_{33})^2}{k_3^0(1-d_3)} + \frac{\sigma_{31}^2}{k_1^0(1-d_1)} + \frac{\sigma_{32}^2}{k_2^0(1-d_2)} \right), \quad (1)$$

$$Y_t = \max_{\langle \tau \le t \rangle} (Y_3^\alpha |_\tau + (\gamma_1 Y_1 |_\tau)^\alpha + (\gamma_2 Y_2 |_\tau)^\alpha)^{\frac{1}{\alpha}}, \quad (2)$$

$$d(Y_t) = \min \left( 1, \left( \frac{n}{n+1} \frac{h_+(Y_t)}{Y_c} \right)^n \right). \quad (3)$$

The parameters $\alpha$ and $\gamma_i$ weigh modes I, II, and III for propagation in mixed mode, while $n$ controls the shape of the potential law. To ensure the energy equivalence between the cohesive model and the fracture mechanics, the Griffith criterion must be met for each mode, as follows:

$$\begin{cases} G_{3c} = Y_c, \\ G_{ic} = \frac{Y_c}{\gamma_i}, \quad \forall i \in \{1, 2\} \end{cases} \quad (4)$$

In order to propose a parameter identification based on wood strengths $f_i$, this work introduces a simplification of Eqs. (1)–(4), leading to the expressions of Eq. (5), as illustrated in Fig. 7. Thus, the cohesive law is derived from $f_i$, $k_i^0$, and $n$.

$$Y_c = \frac{\sigma_{3,c}^2 (2n + 1)^{1/n+2}}{8n(k_3^0)(n+1)} \quad (5)$$

## An adapted multiscale LaTIn strategy for rolling shear failure

Before presenting the strategy, it is important to introduce the reference problem to be solved. The structure under analysis, denoted E and illustrated in Fig. 8a, is a three-layer CLT panel spanning a domain $\Omega$ with boundary $\partial\Omega$. Each layer, labeled P, occupies a domain $\Omega_P$, satisfying $\bigcup_P \Omega_P = \Omega$, and has a boundary $\partial\Omega_P$ with outward normal $\bar{n}_P$. For
<figure><img src="image_10.png" /><figcaption>Fig. 8 An adapted multiscale LaTIn strategy: reference problem and partitioning</figcaption></figure>

each layer P, the displacement $\underline{u}_P \in \mathcal{U}_P$ and the stresses $\underline{\sigma} \in \mathcal{S}_P$ are considered, with $\underline{\epsilon}(\underline{u}_P)$ denoting the symmetric part of the displacement gradient. For each pair of adjacent layers P and P', the interlaminar interface $\Gamma_{PP'} = \partial\Omega_P \cap \partial\Omega_{P'}$ is defined, while the intralaminar interface $\Gamma_P$ is defined within each layer P. The set of layers adjacent to layer P is designated as $\mathbf{V}_P$. The surface force $\underline{F}_d$ is applied in $\partial\Omega_{F_d}$, and the displacement $\underline{u}_d$ is prescribed in $\partial\Omega_{u_d}$, both subsets of $\partial\Omega$. The volume force is denoted as $\underline{f}_{d'}$. Under the assumptions of small perturbations and quasistatic evolution, the problem is systematically solved incrementally at each time step using an implicit temporal integration scheme.

The multiscale LaTIn mixed domain decomposition strategy [20,23] is here used to solve the reference problem previously defined. As proposed in “Mesomodel” section, damageable interfaces can be defined within a layer and their partitioning should be chosen judiciously according to the potential crack patterns. Because in this paper a 3-layer CLT panel is being studied, intralaminar interfaces with a 45° orientation are taken into account (see “Stress analysis” section). To pose the partitioned problem, the entire domain $\Omega$ is divided into subdomains connected by interfaces, as schematized in Fig. 8b. The substructure E occupying the domain $\Omega_E$ is connected to the adjacent subdomain $E'$ by the interface $\Gamma_{EE'} = \partial\Omega_E \cap \partial\Omega_{E'}$. The surface entity $\Gamma_{EE'}$ applies force distributions $\underline{F}_E$, $\underline{F}_{E'}$ and displacement distributions $\underline{W}_E$, $\underline{W}_{E'}$ to E and $E'$, respectively. We note E the set of neighbors of the subdomain E (including the exterior when it applies) and $\Gamma_E = \cup_{E' \in E} \Gamma_{EE'}$.

The objective of the method is to determine the subdomain fields $\underline{u}_E$ (displacement) and $\underline{\sigma}_E$ (stress), as well as the interface fields $\underline{W}_E$ (displacement) and $\underline{F}_E$ (interface force). The scalability of the strategy is ensured by solving a global, linear and reduced problem, constructed from the interfaces and defined as the macroscopic problem (see Eq. (9)). Finally, after partitioning, the governing equations are divided into two groups:

* Linear equations in subdomains and macroscopic admissibility of interfaces whose solutions are elements of the manifold $\mathbf{A}_d$:

* Linear kinematic admissibility of each substructure E

$$ \text{At every point of } \Gamma_E, \underline{u}_{E|\partial\Omega_E} = \underline{W}_{E|\Gamma_E} \qquad (6) $$
* Static admissibility of E

$$ \begin{gathered} \forall (\underline{u}_E^*, \underline{W}_E^*) \in \mathcal{U}_E^0 \times \mathcal{W}_E^0 \quad \text{such that} \quad \underline{u}_{E|\partial\Omega_E}^* = \underline{W}_{E|\Gamma_E}^*, \\ \int_{\Omega_E} \underline{\sigma}_E := \underline{\epsilon}(\underline{u}_E^*) d\Omega = \int_{\Omega_E} f_d \cdot \underline{u}_E^* d\Omega + \int_{\Gamma_E} \underline{F}_E \cdot \underline{W}_E^* d\Gamma \quad (7) \end{gathered} $$

* Linear orthotropic behavior of the substructure E

$$ \text{At every point of } E, \underline{\sigma}_E = \underline{\underline{C}} := \underline{\epsilon}(\underline{u}_E) \quad (8) $$

* Macroscopic admissibility of the interfaces, which is a partial verification of the action-reaction principle

$$ \forall \underline{W}_E^{M*} \in \mathcal{W}_{E E'}^M, \quad \int_{\Gamma_{EE'}} (\underline{F}_E + \underline{F}_{E'}) \cdot \underline{W}_E^{M*} d\Gamma = 0, \quad (9) $$

where the space $\mathcal{W}_{E E'}^M$ (and its counterpart $\mathcal{F}_{E E'}^M$) of interface macroscopic admissible displacements is a parameter of the method described in [20]. These spaces are constructed using an orthogonal projector.

* Local (non-linear) equations in interfaces whose solutions belong to the manifold $L_Γ$:

* Behavior of interfaces between adjacent substructures

$$ \text{At every point of } \Gamma_{EE'} \text{ such that } E' \in V_E, \\ \mathcal{R}(\underline{W}_E, \underline{W}_{E'}, \underline{F}_E, \underline{F}_{E'}) = \underline{0} \quad (10) $$

* Behavior of interfaces transmitting boundary conditions:

$$ \text{If } \Gamma_{EE'} \cap \partial\Omega \neq \{0\}, \text{ at every point of } \Gamma_{E_d} = \partial\Omega_{F_d} \cup \partial\Omega_{u_d}, \\ \mathcal{R}_d(\underline{W}_E, \underline{F}_E) = \underline{0} \quad (11) $$

The behavioral relationship (10) represents the constitutive equation and the local equilibrium of the interfaces between adjacent substructures. This relation can be made explicit for:

* Perfect interface: $\begin{cases} \underline{F}_E + \underline{F}_{E'} = \underline{0} \\ \underline{W}_E - \underline{W}_{E'} = \underline{0} \end{cases}$

* Cohesive interface: $\begin{cases} \underline{F}_E + \underline{F}_{E'} = \underline{0} \\ \mathcal{A}_{PP'}(\underline{W}_E - \underline{W}_{E'}, \underline{F}_E) = \underline{0} \end{cases}$

* Unilateral contact interface: $\begin{cases} \underline{F}_E + \underline{F}_{E'} = \underline{0} \\ \underline{n} \cdot (\underline{W}_E - \underline{W}_{E'}) \geq 0 \quad \text{or} \quad \underline{n} \cdot \underline{F}_E \geq 0 \\ (\underline{n} \cdot (\underline{W}_E - \underline{W}_{E'}))(\underline{n} \cdot \underline{F}_E) = 0 \\ \underline{P} \underline{F}_E = \underline{P} \underline{F}_{E'} = \underline{0} \end{cases}$ where $\underline{n}$ is the unit normal and $\underline{P}$ is the tangential projection operator.

The scheme involves iteratively finding the solution $s_{ref} \in \mathbf{A}_d \cap \mathbf{L}_Γ$, alternating between the two manifolds. First within $\mathbf{A}_d$ where $s = (s_E)_{E \in E} = (\underline{W}_E, \underline{F}_E)_{E \in E}$ is determined,
<caption>**Algorithm 1 Multiscale LaTIn strategy**</caption>

**Construction of the operators**
Computation and assembly macroscopic operator $\mathbb{L}^M$
**for** $t = t_1, \dots, t_N$ **do**
**Initialization of interface quantities**
**for** $n = 0, \dots, N$ **do**
**Linear Stage**
Macro Stage
Micro Stage
**Relaxation**
**Local Stage**
**Error calculation**
**end for**
**end for**

and then within $\mathbb{L}_\Gamma$ by finding $\hat{s} = (\hat{s}_E)_{E \in \mathbb{E}} = (\hat{\mathbb{W}}_E, \hat{\mathbb{F}}_E)_{E \in \mathbb{E}}$. Both stages are connected through the search directions $k_E^-$ and $k_E^+$ as follows:

$$k_E^- (\mathbb{W}_E - \hat{\mathbb{W}}_E) + (\mathbb{F}_E - \hat{\mathbb{F}}_E) = 0, \\ k_E^+ (\mathbb{W}_E - \hat{\mathbb{W}}_E) - (\mathbb{F}_E - \hat{\mathbb{F}}_E) = 0.$$

The iterative schema is briefly presented in Algorithm 1. First, the meshes are read and the macroscopic operators are constructed. Then, for each time step, the interface variables are initialized. The iterative process begins with solving the local stage ($\mathbf{A}_d$), which is divided into a macroscopic stage (Eq. (9)), and then brought to the micro stage for each subdomain (Eqs. (8) and (7)). Once the linear stage is completed, the computed variables are relaxed. Next, the local stage ($\mathbb{L}_\Gamma$) is solved (Eqs. (10) and (11)). Finally, the residual associated with the convergence criterion is evaluated.

Extensibility is ensured by solving Eq. (9), whereas the convergence rate strongly depends on the local search direction $k_E^+$ and the global direction $k_E^-$. In order to ensure convergence, it is appropriate to choose $k_E^+ \to \infty$ for cohesive interfaces. The global search direction $k_E^-$ for perfect and cohesive interfaces is commonly chosen as $k_{E,i}^- = \frac{\mathbb{F}_i}{\mathbb{L}_c}$ and $k_{E,i}^- = 2k_i^0$, respectively [20]. For each new time increment, the method usually initializes the interface variables as $\hat{\mathbb{F}}^t = \hat{\mathbb{F}}^{t-1}$. Because a small perturbation framework is taken into account, in this work it is proposed to use an extrapolation to initialize them. The particular case of a regular monotone loading leads to the following:

$$\hat{\mathbb{F}}^t = \hat{\mathbb{F}}^{t-1} + (\hat{\mathbb{F}}^{t-1} - \hat{\mathbb{F}}^{t-2}). \qquad (12)$$

It is worth emphasizing that in this work, the stiffness of the cohesive zone model ($k_i^0$) associated with the traction-separation law (see “Mesomodel” section) remains constant and represents the physical stiffness of the interface. The parameters varied in the numerical study (see “Improving the convergence rate of the multiscale LaTIn strategy” section) are the search directions ($k_E^+$) and ($k_E^-$) of the LaTIn iterative scheme, which control the convergence rate behaviour of the algorithm but do not modify the physical cohesive law.

In this work, the LaTIn interfaces are placed at the same locations where the mesomodel introduces the physical interfaces (interlaminar adhesive planes and the 45° intralaminar planes for rolling shear). At each LaTIn iteration, each subdomain is solved with the current interface forces; on those same interfaces the cohesive law uses the displacement jumps to update the damage variables; and the resulting tractions are sent back to the macroscopic stage. Since damage is computed exactly on the interfaces that transmit the
<table><caption>**Table 3** Stiffness and strength values for timber (MPa) [40]</caption><thead><tr><th>grade</th><th>E<sub>L</sub></th><th>f<sub>m</sub></th><th>f<sub>L,t</sub></th><th>f<sub>L,c</sub></th><th>f<sub>T,c</sub></th></tr></thead><tbody><tr><td>MSG8</td><td>8000</td><td>14</td><td>4.0</td><td>18</td><td>8.9</td></tr><tr><td>MSG6</td><td>6000</td><td>10</td><td>4</td><td>15</td><td>8.9</td></tr></tbody></table>

kinematic and static information, any loss of stiffness at the mesoscale is seen by the neighbouring subdomains and by the global problem in the next iteration.

All calculations were performed using the parallel finite element code MULTI, developed in C++ at LMPS-ENS (Université Paris-Saclay). The code uses the MPI library for communication between processors and the METIS algorithm for automatic partitioning of the computational domain into balanced sub-structures with minimal connectivity between interfaces. Each processor solves the local problems associated with its substructures, including interfaces, while macroscopic coupling is performed through distributed reduction and assembly operations. The implementation is based on the LaTIn multiscale strategy and runs on high-performance computing environments. The code can be made available for research purposes upon request to the corresponding author.

## Assessment and optimization of the computational approach

This section is devoted to three case studies with increasing complexity, in order to optimize the numerical approach proposed in the previous section, with a focus on the initiation and propagation of cracks. The first goal is to accurately represent the stress state across the thickness of a short-span CLT beam under bending, using the analytical solution from [26] as a reference. Next, the influence of subdomain size, number of elements, and stiffness of the cohesive zone on crack initiation and propagation is analyzed. This step also serves to evaluate the independence of the results with respect to both mesh refinement and substructuring. Finally, an RS failure simulation of a short-span beam will be performed, comparing the numerical results to the experimental tests carried out in [11].

The material properties of the NZS 3603 structural wood grades [40] are taken into account (see Table 3), except as otherwise stated. All simulations are performed with 3D linear 6-node wedge elements.

## Case study 1: assessment of the stress field

To verify the precision of the multiscale LaTIn strategy, the stress field obtained from the MULTI code is compared to the analytical SAM method [26]. A three-layer CLT panel (length $L = 420$ mm, cross-sectional width $b = 50$ mm, layer height $h = 20$ mm and total height $h_T = 60$ mm) is subjected to a three-point bending condition with an applied load $P = 1154$ N. The material properties are as follows: layers 1 and 3 ($i = \{1, 3\}$) are MSG8, while layer 2 ($i = \{2\}$) is MSG6.

The numerical solution is reached using 480 subdomains with a size of $L_{sst} = 10$ mm and an element length of $l_{x1} = l_{x2} = l_{x3} = 2$ mm. This involves a total of 131,040 elements, resulting in 335,664 degrees of freedom (DOF). The normal and shear stress fields at $x = L/4$ are plotted in Fig. 9, showing a satisfactory representation.
<figure><img src="image_11.png" /><figcaption>**Fig. 9** Three-point bending test: analytical and numerical solution at L/4</figcaption></figure>

## Case study 2: representation of the 3D fracture process zone

Before addressing a large-scale structural problem, it is advisable to assess the proposed strategy in a controlled setting where crack initiation and propagation are well characterised. Such a preliminary analysis enables testing the algorithm under reproducible conditions and, if needed, incorporating adjustments that enhance convergence in wood fracture simulations.

For this purpose, we consider the double cantilever beam (DCB) test, a standard benchmark for delamination and to study the FPZ. In addition to the conventional finite element discretisation, we introduce domain substructuring as a variable intrinsic to the domain decomposition framework, thereby assessing its influence on the FPZ representation.

The DCB specimen (see Fig. 10a) has $L = 280$ mm, $b = 20$ mm, $h = 10$ mm ($h_T = 20$ mm) and $a_0 = 100$ mm. The elastic modulus is $E_L = 12000$ MPa. Cohesive interfaces are inserted in the middle plane, allowing crack propagation, with a potential model defined by $G_c = 0.264$ N/mm, $n = 0.5$ and $k_3^0 = E_L/h$. The load is applied through a prescribed displacement $u_d = (0, 0, \pm2.4)$ mm subdivided into 32 steps.

Preliminary results (see Fig. 10b) reveal pronounced nonlinearities and mesh sensitivity, along with a high iteration count to reach a fixed error, which motivates: (i) a mesh and partition convergence study; and (ii) targeted actions to improve the convergence rate.

## Influence of partition and element size on the FPZ

The numerical representation of crack initiation in the double cantilever beam (DCB) test is highly sensitive to the element size within the fracture process zone. Several analytical formulations have been proposed to estimate the characteristic FPZ length and, consequently, to define appropriate mesh sizes [41]. All these models, originally developed to
<figure><img src="image_12.png" /><figcaption>Fig. 10 DCB test: configuration and early results</figcaption></figure>

predict the cohesive zone length, share a general expression that relates the CZM length to the material properties and fracture energy [36,42–44], such as:

$$l_{cz} = M E \frac{G_c}{\sigma_c^2}, \qquad (13)$$

where M is the weighting constant (see [41]), E is the elastic modulus of the layer under study (E_L), G_c is the critical fracture energy, and σ_c is the strength in the i-direction (f_i). Following the domain decomposition strategy, the initial interface stiffness is set as k_i^0 = E_i/h [20]. Equation (13), together with the Eq. (5), provides a reference scale for the process zone, as illustrated in Fig. 11. For the chosen material parameters, the non-dimensional cohesive-zone length lies in the interval:

$$1 \lesssim \frac{l_{cz}}{h} \lesssim 2.6.$$

Since a better transfer of information is obtained when the subdomain size is of the same order as the FPZ, we select $L_{sst} \simeq l_{cz}$, i.e.:

$$1 \lesssim \frac{L_{sst}}{h} \lesssim 2.6.$$

Therefore, in this study, the lengths of the subdomains $L_{sst} = 5, 10, 20$ mm are considered, as shown in Fig. 12. Different substructure lengths along the $x_1$ direction are compared with the stress distribution within the FPZ and with a reference solution obtained from a beam-element model. The results indicate that, for the two configurations considered, the stress distribution is more accurately represented when $L_{sst} = 20$ mm.

Once the subdomain length has been fixed, the next step is to select an element size that allows the FPZ to be properly captured and the load-displacement response of the DCB test to be accurately reproduced. To this end, the engineering expression proposed in [41] is adopted, while following the common recommendation [45,46] of using approximately ten elements across the FPZ. According to these references, the resulting element size is slightly larger than 2 mm; therefore, $l_{x1} = 2$ mm is chosen for the present analysis (Fig. 13).

Finally, considering $L_{sst} = 20$ mm, $l_{x1} = 2$ mm, it is possible to achieve a satisfactory solution, as illustrated in Fig. 14a. An initial linear behavior (steps 1 to 13) is observed, followed by the onset of delamination starting at step 14. However, the huge number of iterations to reach the fixed error ($e = 1 \cdot 10^{-4}$), suggests that the domain decomposition method used could be further enhanced.
<figure><img src="image_13.png" /><figcaption>**Fig. 11** Effect of *n* and *M* on the theoretical cohesive zone length</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12** DCB test: surface tension distribution as a function of the substructure size</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13** DCB test: influence of the mesh and CZM initial stiffness ($L_{sst} = 20$ mm)</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 14** DCB test: final results ($L_{sst} = 20$ mm, $l_{x1} = 2$ mm, $k_3^0 = E_L/h$)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15** DCB test: influence of using $k_E^+ = 0$ (open contact update) for delaminated interfaces on the number of iterations</figcaption></figure>

Improving the convergence rate of the multiscale LaTin strategy

Different possibilities are explored to reduce the number of iterations. The efficient transfer of information between domains is essential to improve the rate of convergence. For this purpose, normally it is possible to enrich the macroscopic space or to enhance the choice of search directions. Here, the second alternative is addressed through: i) implementing $k_E^+ = 0$ over delaminated interfaces in open contact; and ii) updating the macro operators when an interface is fully delaminated. In addition, this paper proposes a new improvement by initializing the interface quantities (see Eq. (12)).

As shown in Fig. 15, the open contact update decreases the total number of iterations by 52%. Fig. 16 shows that the linear initialization of the quantities has a significant effect before the delamination begins, after which the update of macro operators plays a more critical role. Implementing the two last alternatives results in a 38% decrease in the overall number of iterations compared to the open contact update and a 70.5% reduction compared to the initial setting.
<figure><img src="image_18.png" /><figcaption>**Fig. 16** DCB test: influence of updating macro operators on delaminated interfaces and initialization of quantities on the number of iterations</figcaption></figure>

The DCB computations were run on a 64-nodes cluster (AMD OPTERON 6376, 2.3 GHz). Under these conditions, the reference case required 17.8 h (1.0). Enabling only the *open contact* update reduced the time to 8.12 h (0.46), and combining open contact with the update of macro operators further reduced it to 6.87 h (0.39). When open contact is used together with the new field initialization, the time increases to 9.32 h (0.52). Finally, combining open contact, new field initialization, and macro-operator updating brings the time down again to 6.75 h (0.38).

The analyses presented here establish the numerical reliability and scalability of the proposed multiscale strategy. By validating the cohesive behaviour and improving the convergence rate of the LaTIn solver in a well-known benchmark, the framework is now ready to be applied to a more demanding configuration: the simulation of rolling shear failure in CLT panels, characterised by strong nonlinearities and higher computational cost.

**Case study 3: prediction of rolling shear failure**

This case study focuses on the evaluation and prediction of RS failure in a short-span four-point bending test [27], as illustrated in Fig. 3. A three-layer panel is considered, allowing comparison with the experimental tests of [11]. The geometric parameters (see Table 2) are: $h_T$ = 105 mm, $h_i$ = 35 mm, $L$ = 1260 mm, $b$ = 600 mm, $a_1$ = 315 mm, and $a_2$ = 630 mm. The material properties are the same as in the experimental tests [11]: external layers are MSG8 and the central one is MSG6. The load is applied through a prescribed displacement $u_d$ = (0, 0, -22) mm subdivided into 50 steps.

To generate RS interfaces in a preferred direction, in this case at 45° as analysed in "Stress analysis" section, a partitioning with triangular subdomains is proposed, as shown in Fig. 17. Due to the symmetry of the problem, only one-quarter of the beam is modelled. The mesh results in 422,400 wedge elements (1,267,200 DOF), distributed across 1920 subdomains. The interfaces are represented by 166,464 quadrilateral elements (332,928 DOF).
<figure><img src="image_19.png" /><figcaption>**Fig. 17** Four-point bending test: partitioning</figcaption></figure>

This 45° partition was chosen because the analytical stress field in “Stress analysis” section (Figs. 4, 5) shows that the rolling-shear stresses in a three-layer CLT panel concentrate in the middle layer along a direction close to 45°. By aligning the subdomain interfaces with that direction, the LaTIn interfaces coincide with the expected crack path and with the cohesive interfaces of the mesomodel, so the RS crack can be captured with fewer subdomains than with an orthogonal partition [21].

Two types of interface are depicted in Fig. 17: (a) internal interfaces (blue), which may be perfect or damageable depending on their location, and (b) boundary condition interfaces (red), which transmit forces, displacements, and enforce symmetry conditions. The stratified structure, represented by volumetric elements, is shown in violet for the fiber direction and in yellow for the transverse orientation. Damageable internal interfaces allow interlaminar or intralaminar failures, as depicted in Fig. 6. In this example, all interfaces in the upper longitudinally oriented layer are treated as perfect. In contrast, in the transversely oriented layer, the 45° interfaces are referred to as rolling shear interfaces (RSI), while those with normal in the $x_3$ direction can capture adhesive damage (delamination interface, DI) or transverse wood fiber failure (transverse strength interface, TSI). In the bottom longitudinally oriented layer, vertical interfaces capture tension-induced damage due to bending stresses (bending strength interface, BSI). All simulations of this case study were run on the same computing environment described in “Improving the convergence rate of the multiscale LaTIn strategy” section (cluster with AMD OPTERON 6376, 2.3 GHz). The wall-clock time for the fully non-linear configuration, in which RSI, TSI and BSI are all allowed to activate, was about 7.8 h. Earlier test runs with only the rolling-shear interfaces enabled required a shorter time.

Unlike the DCB benchmark, where the standard choice $k_i^0 = E_i/h$ is adopted, the RS simulations use a single interface stiffness $k_i^0 = 300 \text{ N/mm}^3$ for the damageable interfaces of the transverse layer (RSI, DI, TSI). This value is of the same order as $E_L/h$ for the layer involved in rolling shear ($E_L = 6000 \text{ MPa}, h = 35 \text{ mm}, E_L/h \approx 171 \text{ N/mm}^3$) and prevents an overly compliant response. For vertical interfaces in the lower longitudinal layer (BSI), a slightly stiffer value, $k_i^0 = 800 \text{ N/mm}^3$, is used to preserve the initial bending stiffness of the panel, which was noticeably reduced when a lower value was adopted. This higher stiffness remains within a realistic range for the stiffer outer lamellae ($E_L = 8000 \text{ MPa}$) and, as shown in Fig. 18a, does not alter the linear part of the global response.

The cohesive zone models for these interfaces are configured following the parameter identification methodology described in “Mesomodel” section. DI properties are based on adhesive characteristics, whereas TSI and BSI properties are from characteristic wood strength data available in the literature. Using Eq. (5) and setting $n = 2$ to account for the quasi-brittle behaviour of wood, three different RSI are fitted considering a RS strength ($f_r$) equals to 1.6, 1.8 or 2.0 MPa, respectively. For all cases, $\gamma = 0$ is adopted because
<table><caption>**Table 4** Four-point bending test: cohesive interface properties for intra- and interlaminar interfaces</caption><thead><tr><th rowspan="2"></th><th colspan="3">Intralaminar</th><th rowspan="2">TSI</th><th rowspan="2">BSI</th><th>Interlaminar</th></tr><tr><th>RS-1.6</th><th>RS-1.8</th><th>RS-2.0</th><th>DI [47]</th></tr></thead><tbody><tr><td>σ<sub>c</sub> (MPa)</td><td>1.6</td><td>1.8</td><td>2.0</td><td>f<sub>T,c</sub> = 8.9; f<sub>r</sub> = 1.6</td><td>f<sub>m</sub> = 14</td><td>f<sub>DI</sub> = 2.7</td></tr><tr><td>k<sup>0</sup><sub>f</sub> (N/mm<sup>3</sup>)</td><td>300</td><td>300</td><td>300</td><td>300</td><td>800</td><td>300</td></tr><tr><td>γ<sub>c</sub> (N/mm) (Eqn. (5))</td><td>0.010<sup>a</sup></td><td>0.0126</td><td>0.016</td><td>0.31<sup>b</sup></td><td>0.285</td><td>0.40</td></tr><tr><td>n (-)</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>0.5</td></tr><tr><td>γ (-) (Eqn. (4))</td><td>0.0</td><td>0.0</td><td>0.0</td><td>31<sup>c</sup></td><td>0</td><td>1.0</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Fig. 18** Four-point bending test: effect of RS strength on overall panel performance</figcaption></figure>

the interface is perpendicular to the direction of the principal stress. TSI parameters are identified by evaluating the integrity of the central layer, combining rolling shear and transverse tensile strengths and using Eqn. (4) for the determination of $γ$, while the vertical interfaces of the lower layer are defined based on the wood's bending strength (BSI). All parameters derived from the aforementioned methodology are summarised in Table 4.

Firstly, based on the elastic properties of the wood presented in Tab. 3, the linear behavior of the panel can be represented and compared with the analytical solution described in “Finite element analysis for predicting rolling shear failure” section, as shown in Fig. 18a. Next, the effect of RS cracks can be incorporated in isolation (without interaction with other failure scenarios), with the three new curves showing the onset of RS failure and its effect on the overall load capacity. An RS crack is verified to reduce the global stiffness of the panel by 6.5%, a failure that could begin to interact with other modes. Fig. 18b provides a descriptive analysis of RS failure. Once the critical stress threshold is exceeded, the damage rapidly propagates towards the interface, as shown in the instant #2, for the three RS strength values ($f_r$) mentioned above.

After RS failure develops, it is crucial to know whether the crack will spread as delamination (DI) or due to transverse failure of the wood fiber (TSI). The evaluation depicted in Fig. 19a shows that the RSI does not induce adhesive failure. However, it induces a horizontal intralaminar crack. As this crack propagates, outward in the lower part and inward

<footnote><sup>a</sup> Obtained from Eqn. (5) for the prescribed rolling-shear strengths</footnote>

<footnote><sup>b</sup> Computed from Eqn. (5) using the transverse tensile strength $f_{T,c} = 8.9$ MPa</footnote>

<footnote><sup>c</sup> For TSI, the interface is horizontal and may fail either in opening (mode I) or in rolling-shear (mode II). $G_I$ and $G_{II}$ are therefore evaluated for both cases and $γ$ is set as the ratio $G_I/G_{II}$, which yields $γ ≈ 31$ for the adopted strengths</footnote>
<figure><img src="image_21.png" /><figcaption>**Fig. 19** Four-point bending test: simulating RS failure, delamination and bending failure</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 20** Four-point bending test: mid-span moment vs. displacement curves compared to the experimental data from [11]</figcaption></figure>

in the upper part, the global stiffness of the panel decreases, leading to bending failure in the lower layer near the load application point. The sequence of cracks is described in Fig. 19b.

The simulations are compared with the experimental tests of [11] in Fig. 20. With the onset of RS failure, a localized drop in slope is observed, which is visible in the experimental results and simulations. As the load increases, the damage propagates through the wood in the area adjacent to the interlaminar zone, resulting in a decrease in the overall stiffness, both in the experimental tests and in the virtual test. Finally, the tensile stress in the bottom layer causes a rupture in the fiber direction of that layer. From a descriptive point of view, the results can be compared using Fig. 21. The RS failure observed in the virtual test has an orientation and location similar to those in the experiments. Furthermore, crack propagation in wood occurs horizontally. Finally, the failure of the bottom layer near the load area is also consistent with the experimental observations.

In addition to accurately describing the RS failure and the subsequent failures it triggers, the number of subdomains required to represent these mechanisms is drastically reduced due to the partitioning method used. This is evident when comparing to the 6528 sub-domains and 2.9M DOF in a similar problem in [21], or, in other terms, when reducing
<figure><img src="image_23.png" /><figcaption>**Fig. 21** Four-point bending test: Descriptive representation of failure with respect to the experimental tests from [11] at the final instant before unloading</figcaption></figure>

from 8 to 4 subdomains through the layer's thickness, meaning a global reduction of 55% in DOF and 70% in subdomains.

## Discussions

The main findings derived from the previous sections are presented below:

* The results obtained using the domain decomposition method suggest that the linear behaviour of a CLT panel can be adequately captured, with a reasonable agreement observed when compared to analytical solutions available in the literature.

* In the analysis of non-linear behaviour, particularly regarding crack initiation and propagation, several influencing factors can be identified. Alongside the expected effects of mesh and element size, the type of substructuring adopted in the domain decomposition scheme appears to impact the quality of the solution. This indicates the relevance of a substructuring convergence study, similar to mesh convergence analyses commonly used in finite element analysis. In particular, substructures with a length comparable to that of the fracture process zone seem to promote improved information transfer between subdomains, enabling the representation of surface stresses without marked discontinuities.

* Once a reasonable representation of the behaviour was achieved, alternatives were explored to improve the computational efficiency of the algorithm. Some modifications that showed a positive impact on reducing the number of iterations include: consideration of open contact in already delaminated interfaces, updating the macro operators once an interface is fully delaminated, and reinitialising the force fields at the beginning of each load increment step. These adjustments appear to contribute to convergence with fewer iterations.

* Regarding the rolling shear phenomenon, the results suggest that both the mixed-domain decomposition strategy and the substructuring employed provide an adequate representation of the phenomenon, both qualitatively and quantitatively. The observed crack trajectories and stiffness loss show good agreement with experimental data reported in the literature.

* In addition, this methodology enables the investigation of how rolling shear may interact with other failure mechanisms, both intralaminar and interlaminar. Progressive loss of stiffness associated with rolling shear has been observed to trigger additional damage modes, compromising the panel's load-bearing capacity up to the point of fibre-parallel rupture—a phenomenon which, under real conditions, may lead to severe structural failure.
## Conclusions

This work provides a comprehensive numerical framework for predicting Rolling Shear failure in CLT panels, combining cohesive meso-modelling with a mixed multiscale domain decomposition approach. The adapted substructuring, informed by analytical stress field analysis in 3-layer panels, enable a reduction of at least 55% in the number of subdomains and 70% degrees of freedom, while accurately capturing dominant stress directions and failure zones.

The simulation results show strong correlation with experimental data, particularly regarding the initiation and propagation of Rolling Shear failure and its interaction with other failure mechanisms, such as transverse layer fracture and subsequent bending failure. In addition, improvements introduced in the iterative algorithm—such as macro operator updates and variable reinitialisation—led to a reduction of 70.5% in the number of iterations, significantly improving computational efficiency.

Overall, the proposed strategy offers a practical and reliable tool for evaluating the structural performance of CLT panels under complex loading scenarios. While it does not replace experimental testing, it contributes to a deeper understanding of failure mechanisms and supports the development of improved design criteria. Future work will extend this methodology to more complex configurations, including 5- and 7-layer panels and CLT-concrete composite systems.

## Acknowledgements

K. Saavedra acknowledges the financial support from the Chilean National Agency for Research and Development (ANID), FONDECYT regular research project No 1221793. J. Fernández acknowledges the financial support from ANID, PFCHA National Doctorat No 21171988.

## Author Contributions

JF: Simulation, writing—original draft, investigation. KS: conceptualization, supervision, writing—review & editing, investigation. OA: conceptualization, supervision, writing—review & editing. PG: conceptualization, investigation, supervision, writing—review & editing.

## Funding

This research is supported by Chilean National Agency for Research and Development (ANID).

## Data availability

No data was used for the research described in the article.

## Declarations

### Competing interests

The authors declare that they have no competing interests.

Received: 31 March 2025 Accepted: 26 January 2026

Published online: 18 February 2026

## References

1. Asdrubali F, Ferracuti B, Lombardi L, Guattari C, Evangelisti L, Grazieschi G. A review of structural, thermo-physical, acoustical, and environmental properties of wooden materials for building applications. Build Environ. 2017;114:307–32.

2. Brandner R, Flatscher G, Ringhofer A, Schickhofer G, Thiel A. Cross laminated timber (CLT): overview and development. Eur J Wood Wood Products. 2016;74(3):331–51.

3. Ahmed S, Arocho I. Analysis of cost comparison and effects of change orders during construction: study of a mass timber and a concrete building project. J Build Eng. 2021;33:101856.

4. Gagnon S, Bilek E, Podesto L, Crespell P. CLT: introduction to cross-laminated timber. Pointe-Claire: FPInnovations; 2013.

5. European Committee for Standardization. Eurocode 5: design of timber structures-Part 1–1: general-common rules and rules for buildings; 2004.

6. Association A-TEW. Standard for performance-rated cross laminated timber; 2011.

7. Teibinger M, Matzinger I. Construction with cross-laminated timber in multi-storey buildings: focus on building physics: guidelines. Vienna: Holzforschung Austria; 2013.
8. Ren H, Bahrami A, Cehlin M, Wallhagen M. A state-of-the-art review on connection systems, rolling shear performance, and sustainability assessment of cross-laminated timber. Eng Struct. 2024;317:118552.

9. Ehrhart T, Brandner R, Schickhofer G, Frangi A. Rolling shear properties of some european timber species with focus on cross laminated timber (CLT): test configuration and parameter study. In: International network on timber engineering research: proceedings of meeting 48, vol. 2015; 2015. p. 61–76. Timber Scientific Publishing, KIT Holzbau und Baukonstruktionen.

10. Ehrhart T, Brandner R. Rolling shear: test configurations and properties of some European soft-and hardwood species. Eng Struct. 2018;172:554–72.

11. Lim H, Tripathi S, Li M. Rolling shear modulus and strength of cross-laminated timber treated with micronized copper azole type c (mca-c). Constr Build Mater. 2020;259:120419.

12. Nero R, Christopher P, Ngo T. Investigation of rolling shear properties of cross-laminated timber (CLT) and comparison of experimental approaches. Constr Build Mater. 2022;316:125897.

13. Glasner D, Ringhofer A, Brandner R, Schickhofer G. Rolling shear strength of cross laminated timber (CLT)-testing, evaluation, and design. Buildings. 2023;13(11):2831.

14. Cheng Z-Q, Liu H, Tan W. Advanced computational modelling of composite materials. Eng Fract Mech. 2024;305:110120.

15. Okereke M, Akpoyomare A, Bingley M. Virtual testing of advanced composites, cellular materials and biomaterials: a review. Compos B Eng. 2014;60:637–62.

16. Flores ES, Saavedra K, Hinojosa J, Chandra Y, Das R. Multi-scale modelling of rolling shear failure in cross-laminated timber structures by homogenisation and cohesive zone models. Int J Solids Struct. 2016;81:219–32.

17. De Moura M, Dourado N, Morais J. Crack equivalent based method applied to wood fracture characterization using the single edge notched-three point bending test. Eng Fract Mech. 2010;77(3):510–20.

18. Gustafsson PJ, Henrik D. Perpendicular to grain stiffness of timber cross sections as affected by growth ring pattern, size and shape. Eur J Wood Wood Prod. 2013;71(1):111–9.

19. Ladevêze P, Nouy A, Loiseau O. A multiscale computational approach for contact problems. Comput Methods Appl Mech Eng. 2002;191(43):4869–91.

20. Kerfriden P, Allix O, Gosselet P. A three-scale domain decomposition method for the 3D analysis of debonding in laminates. Comput Mech. 2009;44(3):343–62.

21. Saavedra Flores EI, Saavedra K, Hinojosa J, Chandra Y, Das R. Multi-scale modelling of rolling shear failure in cross-laminated timber structures by homogenisation and cohesive zone models. Int J Solids Struct. 2016;81:219–32.

22. Ruda A, Louf F, Boucard P-A, Mininger X, Verbeke T. A parallel implementation of a mixed multiscale domain decomposition method applied to the magnetostatic simulation of 2d electrical machines. Finite Elem Anal Des. 2024;235:104136.

23. Ladevêze P, Allix O, Deü J-F, Lévêque D. A mesomodel for localisation and damage computation in laminates. Comput Methods Appl Mech Eng. 2000;183(1–2):105–22.

24. Allix O, Feissel P, Thévenet P. A delay damage mesomodel of laminates under dynamic loading: basic aspects and identification issues. Comput & Struct. 2003;81(12):1177–91.

25. Li H, Wang L, Wei Y, Wang BJ, Jin H. Bending and shear performance of cross-laminated timber and glued-laminated timber beams: a comparative investigation. J Build Eng. 2022;45:103477.

26. Winter S, Kreuzinger H, Mestek P. Flächen aus brettstapeln, brettsperrholz und verbundkonstruktionen. Technical report, Technische Universität München; 2008.

27. European Committee for Standardization. Structural timber-strength classes; 2016.

28. Platten KH. Scheiben und schalen-ein berechnungsmodell für gängige statikprogramme. Bauen mit Holz. 1999;1:34.

29. Timoshenko S. On the correction factor for shear of the differential equation for transverse vibrations of bars of uniform cross-section. Philosophical Magazine Series 1; 1921.

30. Mestek P. Punktgestützte flächentragwerke aus brettsperrholz (bsp)-schubbemessung unter berücksichtigung von schubverstärkungen. PhD thesis, Technische Universität München; 2011.

31. Hochreiner G, Füssl J, Eberhardsteiner J. Cross-laminated timber plates subjected to concentrated loading. Strain. 2014;50(1):68–81.

32. Norris C. Strength of orthotropic materials subjected to combined stresses. Forest products laboratory: technical report; 1962.

33. Tsai SW, Wu EM. A general theory of strength for anisotropic materials. J Compos Mater. 1971;5(1):58–80.

34. Mestek P, Dietsch P. Design concept for CLT-reinforced with self-tapping screws. In: Focus solid timber solutions-European conference on cross laminated timber (CLT); 2013. p. 103–18.

35. Barenblatt Gl. The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech. 1962;7:55–129.

36. Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids. 1960;8(2):100–4.

37. Coureau J-L, Morel S, Dourado N. Cohesive zone model and quasibrittle failure of wood: a new light on the adapted specimen geometries for fracture tests. Eng Fract Mech. 2013;109:328–40.

38. Danielsson H, Gustafsson PJ. Fracture analysis of glued laminated timber beams with a hole using a 3D cohesive zone model. Eng Fract Mech. 2014;124:182–95.

39. Allix O, Lévêque D, Perret L. Identification and forecast of delamination in composite laminates by an interlaminar interface model. Compos Sci Technol. 1998;58(5):671–8.

40. Standards New Zealand. Timber structures standard; 1993.

41. Turon A, Davila CG, Camanho PP, Costa J. An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Eng Fract Mech. 2007;74(10):1665–82.

42. Falk ML, Needleman A, Rice JR. A critical evaluation of cohesive zone models of dynamic fractur. Le J de Physique IV. 2001;11(PR5):5–43.

43. Rice JR. The mechanics of earthquake rupture, vol. 720. Providence: Division of Engineering, Brown University; 1979.

44. Hillerborg A, Modéer M, Petersson P-E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res. 1976;6(6):773–81.

45. Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Eng Fract Mech. 2002;69(7):813–33.
46. Carpinteri A, Cornetti P, Barpi F, Valente S. Cohesive crack model description of ductile to brittle size-scale transition: dimensional analysis vs. renormalization group theory. Eng Fract Mech. 2003;70(14):1809–39.

47. Oñate A, Fernández J, Saavedra K. Numerical simulation on the influence of delamination on the critical buckling load of CLT panels. Eng Struct. 2023;286:116029.

**Publisher’s Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.