Macroscopic Model for Spatial Timber Plate Structures with Integral Mechanical Attachments

The integral mechanical attachment (IMA), which is used to join timber elements with only wood–wood interlocking mechanisms (i.e., without connectors such as screws, nails, adhesives, and metal fasteners), is considered one of the oldest joinery techniques (Messler 2006; Robeller and Weinand 2015). In the last decade, the introduction of digital geometry processing to the design of complex freeform structures and the application of computer numerically controlled (CNC) machines for fabrication and robotic arms for assembly have resurrected the IMA technique. Moreover, engineered timber products such as laminated veneer lumber (LVL) planar boards have addressed the deficiencies of solid woods through their enhanced mechanical properties (Rezaei Rad et al. 2019a).

This paper focusses on two specific types of IMAs: multiple tab-and-slot dovetail joints [Fig. 1(a)], and mortise-and-tenon or through-tenon (TT) joints [Fig. 1(b)]. The flexural behavior of both types of IMAs was studied by Roche et al. (2015), Roche (2017), and Nguyen et al. (2019a). Rezaei Rad et al. (2019a, b) investigated the performance of TT joints under tension and in-plane (IP) (edgewise) loads, respectively. Additionally, the influence of dovetail and TT joints on the global behavior of timber folded plate surfaces was investigated by Stitic et al. (2018, 2019). They concluded that the semirigid behavior of IMAs considerably affects the performance of spatial timber plates and therefore should be considered in the design process.

An interdisciplinary design framework that combines the knowledge of architects, engineers, and computer and robotic scientists is critical to the construction of novel timber plate structures. Therefore, the laboratory for timber construction (IBOIS) at École Polytechnique Fédérale de Lausanne (EPFL) recently designed and participated in the construction of the first timber folded plate theatre using only wood–wood IMAs. Gamerro et al. (2018) experimentally studied the flexural behavior of the IMA joints used in the structure and found that they can adequately satisfy the relevant design requirements.

The IMA technique has also been used to construct freeform surfaces. This is partly the result of a biologically inspired computational design strategy that was introduced by Magna et al. (2013) and Krieg et al. (2015a, b). The IMA finger joints, including pairs of parallel cross-screws, were used to achieve the desired mechanical properties. The framework addressed challenges in the fabrication, assembly, and service-life performance of the structure. With the support of IBOIS researchers, it was implemented in the design and construction of freeform integrally attached timber plate (IATP) structures using only wood–wood connections (Fig. 2). The lightweight construction system uses CNC routers to fabricate timber plates equipped with wood–wood connections along the plate edges. The computational framework used for prefabrication and assembly was described by Robeller et al. (2016) and the parts that are relevant to the current work were discussed further in the section in which the IATP system and assembly process are described.

The numerical models used to analyze plate structures typically lie somewhere along a spectrum of complexity. Finite-element (FE) models with three-dimensional (3D) solid elements, which can be viewed as existing on one end of that spectrum, are detailed, complex, and require a high level of sophistication. However, they do provide the benefit of enabling an evaluation of the state of stress across the plate thickness. In less complex models, the plate thickness is neglected, and two-dimensional (2D) shell elements can be used to simulate the behavior of the plates. Both of the aforementioned approaches provide enhanced simulation capabilities to capture stress–strain response and local failure modes. However, they are not without their challenges. First, a refined mesh discretization is required. Furthermore, a high level of expertise is needed to develop reliable models and to postprocess and interpret the results (Bao et al. 2008; Pantò et al. 2016). In other words, detailed FE models are not suitable for most industry applications (Bao and Kunnath 2010; Kareem and Pantò 2019) because of the associated computational expense, especially for large-scale structures (Caliò and Pantò 2014). Moreover, it has been shown that more-detailed models do not always provide improved results compared to simpler models in large-scale systems (Caliò et al. 2012; Pantò et al. 2016). On the other hand, more-simplified models, which utilize only beam and spring elements, can be used to simulate structural response. These models, which are generally referred to as macromodels, are computationally inexpensive and involve less complexity in their construction and analysis and in postprocessing the results (Pantò et al. 2019). The primary goal of the macromodeling approach is to simulate the global behavior of a structure (Pantò et al. 2018).

This research introduces a novel macromodeling approach for timber plate structures which can be implemented in structural engineering practice. The proposed macromodel is only valid for linear static analysis. Fig. 3 presents an overview of the contents of the paper. The mechanical characterization chapter describes the assembly of the IATP system, the kinematic degrees of freedom (DOFs) associated with timber plates and IMAs, the force flow, and equilibrium equations [Fig. 3(a)]. The macromodel is then introduced and the mechanical properties associated with the beam/spring elements are determined [Fig. 3(b)]. Next, the response of an individual plate and a full-scale subassembly is simulated using the macromodel and the results are verified using experimental data (only for the full-scale subassembly) and FE models using shell elements (for both the individual plate and the full-scale subassembly) [Fig. 3(c)]. Lastly, the conclusions of the study are discussed.

The IATP structural system consists of several planar elements with rectangular/parallelogram shapes, which are connected by TT and dovetail joints. In brief, the target surface of a freeform structure is transformed into a target geometry [Figs. 4(a and b)] and the configuration of the timber plates is determined [Fig. 4(c)]. Double-layered IATP systems include the assembly of multiple four-sided hexahedron-shaped boxes, denoted $C_i$ in Fig. 4(d), where $i=1,2,\dots ,8$. Each box consists of top and bottom plates [denoted $T_i$ and $B_i$ in Fig. 4(d), respectively] and two cross plates [denoted cross longitudinal, ${\mathrm{CL}}_i$, and cross transverse, ${\mathrm{CT}}_i$, in Fig. 4(d)]. Within each box, the tenons are located along the edges of the $T_i$ and $B_i$ plates, and their mates (slots) are located in the ${\mathrm{CL}}_i$ and ${\mathrm{CT}}_i$ plates. The ${\mathrm{CL}}_i$ and ${\mathrm{CT}}_i$ plates are connected via multiple dovetail-shaped joints. Each box is then connected to its neighbor with only the TT joints. Each timber box is labeled with a normal vector (${\mathbf{nor}}_i$, $i=1,2,\dots ,8$). The collective set of normal vectors within a structure can be algorithmically determined so that they comply with the target surface to achieve the desired geometry.

A 1-degree-of-freedom process is used to assemble the plates (Robeller et al. 2016), meaning that there is only one translation vector in a 3D space to (dis)assemble adjacent plates; and accordingly, the joint geometry reduces the assembly degrees-of-freedom to that single vector. Two primary steps are involved in assembling an IATP system. In the first step, the plates are assembled to form each box (intrabox assembly). The ${\mathrm{CL}}_i$ and ${\mathrm{CT}}_i$ plates first are connected along vector ${\mathbf{u}}_i$. Fig. 4(d) shows ${\mathrm{CL}}_8$ being connected to ${\mathrm{CT}}_8$ along vector ${\mathbf{u}}_8$. Next, the $T_i$ and $B_i$ plates are connected simultaneously to the corresponding ${\mathrm{CL}}_i$ and ${\mathrm{CT}}_i$ plates along the internal assembly vector, ${\mathbf{V}}_{i,\mathrm{int}}$ [Fig. 4(d), ${\mathbf{V}}_{8,\mathrm{int}}$]. The neighboring boxes are assembled as part of the second step (interbox assembly). Each box is simultaneously connected to its two nearest neighbors along the same external assembly vector (${\mathbf{V}}_{i,\mathrm{ext}}$). Fig. 4(d) shows $C_8$ being simultaneously connected to $C_5$ and $C_7$ along vector ${\mathbf{V}}_{8,\mathrm{ext}}$.

Given the primary focus on the global (structural) scale, the overall behavior is deemed to be governed by the in-plane and out-of-plane (OP) response of the IATP elements and the behavior of the TT joints. To formulate the equations governing the IP and OP behavior, only the relative displacements associated with the perimeter of the plate are considered. In other words, only the DOFs associated with nodes located at the IATP perimeter are used. This reduces the size of the system of equations to be solved in the stiffness method, which in turn lowers the overall computational expense. The governing equations are presented in the next two subsections and in Appendices I–VI.

Depending on the location of each timber box, the plate types (e.g., $T_i$ and $B_i$ or ${\mathrm{CL}}_i$ and ${\mathrm{CT}}_i$), and the geometry (parallelogram or rectangular), an IATP element can be subjected to three main IP load cases:

The kinematic formulation of the IP behavior of a single IATP element was adapted from (1) Christovasilis and Filiatrault (2013) and Casagrande et al. (2016), which focused on describing the in-plane behavior of sheathing panel elements in wood-frame shear walls; and (2) Pantò et al. (2018), which focused on characterizing the in-plane behavior of masonry vaults. Additional assumptions were made to ensure that the formulation is compatible with the mechanical constraints of the IATP system.

The equations governing the IP behavior of the $T_i$ and $B_i$ plates, the associated local and global coordinate systems, and the applied loads in Figs. 5(a and b), are formulated in this section and in Appendices I–IV. A similar formulation can be established for the ${\mathrm{CL}}_i$ and ${\mathrm{CT}}_i$ plates [Figs. 5(c and d)]. Five kinematic DOFs, including uniform shear (${\gamma }_{\mathrm{XZ}}$) [Fig. 6(a)] and uniform flexural deformations ($\varphi$) [Fig. 6(b)], two rigid-body translations ($U_X$ and $U_Y$) [Fig. 6(c)], and rigid-body rotation ($\mathrm{\Theta }$)[Fig. 6(c)], are defined for each IATP element. The vector ${\mathbf{D}}_{\mathrm{IP}}={\{\begin{array}{ccccc}{U}_X& U_Z& {\gamma }_{XZ}& \mathrm{\Theta }& \varphi \end{array}\}}^{\mathrm{T}}$ is used to describe these kinematic DOFs when the element responds within the linear elastic range. The corresponding force vector ${\mathbf{F}}_{\mathrm{IP}}={\{\begin{array}{ccccc}{P}_X& P_Z& P_{XZ}& M_{Y,\mathrm{rgd}}& M_{Y,\mathrm{flx}}\end{array}\}}^{\mathrm{T}}$ is also defined, where $P_X$ and $P_Y$ are the IP forces along the global $X$- and $Z$-axes, respectively, and $P_{XZ}$ is the IP shear force. The IP moment, $M_Y$, is decomposed into two components: $M_{Y,\mathrm{rgd}}$ which comes from the neighboring boxes and causes IP rigid-body rotations, and $M_{Y,flx}$, which is the IP moment that causes flexural deformations.

The IP displacement of node $i$, which represents the $i$th IMA located at the perimeter of the IATP element (Fig. 6), is described by the vector ${\mathbf{u}}_{\mathrm{IP},i}={\{\begin{array}{cc}{X}_{\mathrm{IP},i}& Z_{\mathrm{IP},i}\end{array}\}}^{\mathrm{T}}$ in the global coordinate system (GCS) $(X,Z)$. The vector includes contributions from shear deformation $u_{\mathrm{IP},i}^{{\gamma }_{XZ}}$, flexural deformation $u_{\mathrm{IP},i}^{\varphi }$, rigid-body translation $u_{\mathrm{IP},i}^{UX,UZ}$ and rigid-body rotation $u_{\mathrm{IP},i}^{\mathrm{\Theta }}$. The principle of superposition is used to compute $u_{\mathrm{IP},i}$ [Eq. (1)], which is deemed valid because only linear elastic responses are considered

To determine the location of node $i$ at each load–displacement step, a local coordinate system (LCS) $(x^{\prime },z^{\prime })$ and a coordinate vector, ${\mathbf{\xi }}_{\mathrm{IP},i}={\{\begin{array}{cc}{x}_i^{\prime }& z_i^{\prime }\end{array}\}}^{\mathrm{T}}$ are defined. The latter specifies the position of node $i$ in the LCS with respect to its origin at each load step. The origin of the LCS is located at the center of the panel (Fig. 6) and, in the undeformed and undisplaced condition, is assumed to be parallel to the GCS. As noted previously, the IP deformation is computed as the vector sum of the deformations from the various sources. Consequently, ${\mathbf{\xi }}_{\mathrm{IP},i}$ is updated at each load step and used for the calculations in the next step. The detailed mathematical formulation, including the relationship between (1) the ${\mathbf{u}}_{\mathrm{IP},i}$ and ${\mathbf{\xi }}_{\mathrm{IP},i}$ vectors, and (2) the ${\mathbf{u}}_{\mathrm{IP},i}$ and ${\mathbf{D}}_{\mathrm{IP}}$ vectors, is presented in Appendix I.

To compute the IP stiffness of an IATP element, a virtual displacement is associated with node $i$where ${\mathbf{B}}_{\mathrm{IP}}(D_{\mathrm{IP}})$ = matrix representing the partial derivative of IP displacements with respect to $D_{\mathrm{IP}}$. The calculation needed to obtain each entry of the ${\mathbf{B}}_{\mathrm{IP}}(D_{\mathrm{IP}})$ matrix is given in Appendix II.

Assuming that there are $m$ TT joints along the edge of the timber plate, the virtual work principle ($\delta W_{\mathrm{external}}=\delta W_{\mathrm{internal}}$) is applied. The associated mathematical formulation is provided in Appendix III. The global force vector is expressed in terms of the internal forceswhere ${\mathbf{F}}_{IP}$ = global force vector; ${\mathbf{B}}_{IP}$ = matrix of partial derivatives of the IP deformation (Appendix II); $f_{\mathrm{IMA},\mathrm{IP},i}={\{\begin{array}{cc}{p}_{\mathrm{IMA},\mathrm{IP},i}& q_{\mathrm{IMA},\mathrm{IP},i}\end{array}\}}^T$ = internal force vector of each IMA joint [Fig. 5(b)]; and ${\int }_{L}N(x)dx$, ${\int }_{L}V(x)dx$, and ${\int }_L{M}_z(x)dx$ = axial force in global $X$-direction, IP shear force, and bending moment, respectively (Appendix III).

Finally, the IP stiffness of the IATP, including the plate ($K_{\mathrm{IP},\mathrm{Plate}}$) and IMAs ($K_{\mathrm{IP},\mathrm{IMA}}$), is (with additional details provided in Appendix IV)

In summary, the IP behavior of a timber plate equipped with TT IMAs can be simulated while accounting for the axial ($K_{\mathrm{Plate}.ax,0}$), shear ($K_{\mathrm{Plate}.shr}$), and flexural ($K_{\mathrm{Plate}.flx}$) stiffnesses of the timber plate and the semirigid behavior of the TT joints under tensile ($K_{\mathrm{IMA},\mathrm{tns}}$) and edgewise shear ($K_{\mathrm{IMA},\mathrm{edg}}$) loads. The calculation of these stiffnesses are discussed in the macromodel section.

Fig. 7 shows the loads applied to the $C_8$ box. Considering (1) the loads imposed by the neighboring boxes [$C_5$ and $C_7$ in Fig. 4(d)] on the $C_8$ box, and (2) any external load applied directly to $C_8$, an IATP element can be subjected to two main OP load cases:

Because the plate thickness is small relative to the perimeter dimensions, the OP behavior of IATPs is characterized using multiple beam elements (i.e., to capture biaxial bending) (Choi and Park 2010). Fig. 8 shows the beam idealization of an IATP element. The beam elements, which represent a strip of the plate, are along the directions parallel [fiber-parallel (FPa)] and perpendicular [fiber-perpendicular (FPe)] to the fibers. The contribution to OP bidirectional bending from each fiber-parallel and fiber-perpendicular beam is ${\alpha }_i$ and ${\beta }_i$, respectively and is reflected in the stiffness matrix of the IATP element. Assuming that the fiber-parallel and fiber-perpendicular directions have $n$ and $m$ strips, respectively, it follows that ${\sum }_{i=1}^n{\alpha }_i+{\sum }_{j=1}^m{\beta }_j=1.0$, where $0<{\alpha }_i$ and ${\beta }_j<1.0$. In Appendix V, the contribution factors ${\alpha }_i$ and ${\beta }_i$ are analytically calculated for a single IATP element. However, for more complex assemblies with large numbers of DOFs, a numerical solution is required (Bozorgnia and Bertero 2004).

The OP behavior of each IATP is represented by that of the beam (strip) elements. To calculate the stiffness matrix of each beam (${\mathbf{K}}_{\mathrm{Strip}}$), the end nodes are indexed as 1 and 2 in Fig. 9(a) and the OP behavior is described using 6 DOFs in the GCS. The vector ${\mathbf{D}}_{\mathrm{OP}}={\{\begin{array}{cccccc}{u}_1& v_1& {\theta }_1& u_2& v_2& {\theta }_2\end{array}\}}^{\mathrm{T}}$ describes the kinematics associated with these DOFs when the element responds within the linear elastic range. The components of the vector ${\mathbf{D}}_{OP}$ are shown in Fig. 9(a). The terms $u_1$, $v_1$, $u_2$, and $v_2$ are the translations, and ${\theta }_1$ and ${\theta }_2$ are the rotations at the end nodes. Similar to the $D_{OP}$ vector, the element force components are described by the global vector ${\mathbf{F}}_{OP}={\{\begin{array}{cccccc}{p}_1& q_1& m_1& p_2& q_2& m_2\end{array}\}}^{\mathrm{T}}$ [Fig. 9(b)]. The details of the formulation of the beam element used to characterize the OP behavior and the corresponding force-deformation relationship are provided in Appendix VI.

The stiffness matrix for each beam element, ${\mathbf{K}}_{\mathrm{Strip}}$, iswhere ${\mathbf{B}}_{OP}$ = matrix representing the partial derivative of the OP displacements with respect to ${\mathbf{D}}_{\mathrm{OP}}$; $E$ = modulus of elasticity; $I$ = moment of inertia; and $L_0$ = initial length of the element.

Once ${\mathbf{K}}_{\mathrm{Strip}}$ is determined, the principle of superposition is used to obtain the stiffness matrix for the entire IATP. This is taken as the sum of the stiffnesses of the IMAs and the strip elements along the fiber-parallel and fiber-perpendicular directions while accounting for the contribution of each strip element to the OP behaviorwhere ${\mathbf{K}}_{\mathrm{Strip},i}$ and ${\mathbf{K}}_{\mathrm{Strip},j}$ = stiffness matrices for the beams parallel and perpendicular to the fiber orientation, respectively; ${\alpha }_i$ and ${\beta }_i$ = associated flexural stiffness contribution factors; and $K_{\mathrm{IMA},\mathrm{tns}}$, $K_{\mathrm{IMA},\mathrm{flt}}$, and $K_{\mathrm{IMA},\mathrm{flx}}$ capture the semirigid behavior of the TT joints under tensile and flatwise shear loads and out-of-plane moments, respectively.

According to Nilson et al. (2010) and Hillerborg (1996), computing the flexural stiffness contribution of the orthogonal strips analytically is approximate. In fact, the load path is determined by multiple parameters, including the boundary conditions, distribution of the surface and edge loads, and the material orientation. Therefore, numerical simulations, which are described in the verification section, are needed to calculate ${\alpha }_i$ and ${\beta }_i$, and, consequently, the OP stiffness matrix of the strips and their contribution to the overall OP behavior.

Using the fabrication contours of an IATP [Fig. 10(a)], the midsurface of each 3D plate is established [Fig. 10(a)]. Next, the macromodel is generated from the midsurface using four main steps: (1) using the corner nodes of the midsurface to create a planar polygon (hereafter called boundary elements) [Fig. 10(b)], (2) identifying the location of the IMAs in the boundary elements [Fig. 10(c)], (3) dividing the polygon into equal segments along the fiber-parallel and fiber-perpendicular directions and adding inner lines at the interface between segments (hereafter called inner beam elements) [Fig. 10(d)], and (4) splitting the boundary elements at each intersection and connecting the end nodes of the inner beam elements, TT joints, and the contour corners [Fig. 10(e)].

The proposed macromodel [Fig. 11(c)] is formulated to be compatible with different quadrilateral geometries and load cases. It consists of boundary elements [Fig. 11(d)], inner beam elements [Figs. 11(e and h)], two (uniaxial) shear springs [Figs. 11(j and k)], two-node link elements representing the joints [Fig. 11(g)], and uniaxial springs distributed along the boundary elements [Figs. 11(f and i)], which capture the tension–compression behavior of the plate. Overall, structural stability is maintained by connecting the inner and perimeter beam elements, the uniaxial springs, and the IMAs.

The boundary elements are pin-ended and free to rotate about their local $x^{\prime }$- and $y^{\prime }$-axes [Fig. 11(d) for the LCS]. However, the fiber-parallel [Fig. 11(e)] and fiber-perpendicular [Fig. 11(h)] inner beams have a flexural release about their local $x^{\prime }$ and $z^{\prime }$, respectively. In other words, the boundary elements [Fig. 11(d)] are free to rotate in the IP and OP directions, whereas the inner beams [Figs. 11(e and h)] provide OP flexural resistance and are only free to rotate in the IP direction. It is within this context that the macromodel incorporates the mechanical properties of the inner beam elements to simulate the OP kinematic behavior and is independent of the boundary elements. To simulate the IP kinematic response, Eq. (4) and two modeling strategies are used: (1) a pair of uniaxial springs in the fiber-parallel and fiber-perpendicular directions is used to simulate the shear behavior ($K_{\mathrm{Plate}.\mathrm{Shr}}$) [Figs. 11(j and k)]; and (2) a set of distributed uniaxial tension–compression springs along the edges of the quadrilateral macromodel [Figs. 11(f and i)] is used to simulate the combination of direct axial ($K_{\mathrm{Plate}.\mathrm{ax}}$) and flexural ($K_{\mathrm{Plate}.\mathrm{flx}}$) loads. These two modeling strategies were also adopted by Christovasilis and Filiatrault (2013), Casagrande et al. (2016), and Caliò et al. (2012).

In IATP structures, the joints and plates are one unit. In other words, because no additional connectors such as screws or fasteners are employed, the joints and the associated plate cannot be separated. To better understand their behavior, experiments were performed by Roche (2017) and Rezaei Rad et al. (2019a, b), which included both the joints and the neighboring region of the plates in the test specimen. With this type of experiment, the behavior of what is referred to as the joint region is examined and the initial stiffness from an idealized backbone curve is assigned to the two-node link element [Fig. 11(g)]. Consequently, the use of link elements to explicitly describe the behavior of the joint region can account for the local behavior of IATPs.

Because the IMA boundary element behavior is simulated using the aforementioned two-node link elements [Fig. 11(g)], the primary role of the boundary beams is to provide overall stability, and are therefore modeled as rigid. These rigid boundary elements are also pin-ended so that they can freely rotate about their local IP and OP axes [Fig. 11(d)]. Moreover, by using pin-ended rigid boundary elements, the actual geometry of the timber plate is represented in the macromodel. Additionally, because the end nodes have flexural releases [Figs. 11(e and h)], the behavior of the macromodel is independent of the cross-section size of the boundary elements. Therefore, the inner beam elements and uniaxial shear and tension–compression springs encapsulate the mechanical properties of the macromodel. This modeling strategy makes the macromodel compatible with different polygon geometries and load cases. A similar strategy was adopted by Caliò et al. (2012), Caliò and Pantò (2014), and Pantò et al. (2017).

To determine the cross-section properties of each inner beam element, the associated tributary area [Figs. 11(e, f, h, and i)] is computed, and consequently, the stiffnesses of the shear springs, tension–compression springs located around the edges, and the inner beams are determined. The IP flexural stiffness of the plate [$K_{\mathrm{Plate}.\mathrm{flx}}$ in Eq. (4b)] is implicitly represented using the aforementioned uniaxial tension–compression springs [Figs. 11(f and i)]. The flexural stiffness of the plate can be expressed by the uniaxial tension–compression springs (Fig. 12) (Caliò et al. 2012; Pantò et al. 2016). In other words, in addition to simulating the tension–compression behavior, these springs can capture the IP flexural deformation of the element (Fig. 12).

Beech BauBuche LVL (Blaß and Streib 2017) (Pollmeier, Thuringia, Germany) was used as the construction material. Each timber board consisted of 13 3-mm-thick peeled veneer layers, 15% of which (2 of 13 layers) were aligned perpendicular to the fiber direction (referred to as cross-plies). Table 1 summarizes the mechanical properties of the BauBuche LVL material. LVL boards are engineered timber products with homogeneous and orthotropic material properties (Roche 2017; Stitic et al. 2019; Hassanieh et al. 2017).

The detailed FE model was constructed using Abaqus/CAE version 6.14. The midsurface of each timber panel was represented by a polysurface element which was divided into meshes and represented using shell elements. Only linear elastic response was considered for the shell element. Orthotropic material properties (Table 1) were assigned to the shell elements. A quad configuration and finite strain S4R elements were used for meshing. The mesh size was smaller near the connections, with a seed size of 20 mm, and larger in the other areas of the IATP, with a seed size of 50 mm. The connectors used to simulate the TT joints were modeled using the Bushing element.

The macromodels were constructed in OpenSees version 2.5.0 (Mazzoni et al. 2013). A uniaxial elastic material with infinite stiffness was assigned to the boundary elements. An elastic beam–column element was used for the inner beams based on their associated tributary area and the timber properties for each orthogonal direction. The dovetail joints, which form a continuous connection between plates ${\mathrm{CL}}_i$ and ${\mathrm{CT}}_i$, were assumed to be rigid (Nguyen et al. 2019a). The connectors used to simulate the TT joints were modeled using a two-node link element. To simplify the modeling process, the same division used to define the inner beam elements was employed to determine the location of the uniaxial tension–compression springs, which were distributed around the perimeter. Consequently, these springs were placed at the end nodes of each inner beam element.

Fig. 13(a) shows the geometric configuration of the IATP plates inside a box, which is based on the average dimensions of the elements in the full-scale prototype (e.g., Fig. 2). Because of the assembly process for a four-sided box, the ${\mathrm{CL}}_i$ and ${\mathrm{CT}}_i$ plates are subjected primarily to IP loads, whereas the $T_i$ and $B_i$ plates are subjected to OP loads. Therefore, the macromodels representing the IP behavior of the ${\mathrm{CL}}_i$ and ${\mathrm{CT}}_i$ plates [Fig. 13(b)] and the OP behavior of $T_i$ and $B_i$ plates [Fig. 13(c)] were constructed and the load–deformation behavior was compared with that of the corresponding detailed FE models. For the macromodels, Fig. 14 shows the section size of the inner beam elements, as well as the stiffness values for the uniaxial tension–compression and shear springs. The inner beam elements were generated by dividing the ${\mathrm{CL}}_i$ and ${\mathrm{CT}}_i$ plates into three equal segments along the fiber-parallel and fiber-perpendicular directions. Furthermore, the $T_i$ and $B_i$ plates were divided into three and five equal segments along the fiber-parallel and fiber-perpendicular directions, respectively. The sensitivity of the load–deformation behavior to variations in the number and configuration of inner beam elements (Fig. 15) was evaluated, and the results are discussed in the next two subsections.

For the IP load case [Fig. 13(b)], the force corresponding to the ultimate limit state (ULS) is calculated using Eq. (7) (CEN 2008). The shear strength of the timber boards is the weakest resisting element and control the behavior (Table 1)where ${\tau }_d$ and $f_{v,d}$ = applied stress and design strength values, respectively $A=\mathrm{24,000}{\mathrm{mm}}^2$ = effective cross-section area; $V_d$ = shear force; $f_{v,0,k}$ = characteristic shear strength (Table 1); and $K_{h,v}$, $K_{cr}$, and $K_{\mathrm{mod}}$ = modification factors (Table 2) in accordance with CEN (2008).

The IP ULS load $V_d=176.8\mathrm{KN}$ is then computed and distributed along the element height at each IMA [Fig. 13(b)]. The boundary conditions and location of the applied loads are shown in Fig. 13(b).

To compare the response of the FE and macromodels, the deformed shape along the element height for the different configurations defined in the sensitivity analysis are presented in Fig. 16. The deformed shape obtained from the macromodel was similar to the one obtained from the FE model. However, the deflections in the macromodel were slightly higher in the lower portion of the element and slightly lower in the upper portion. These differences are attributed to the different element formulations in the macro and FE models. Specifically, the FE model takes advantage of plate theory, in which nonlinear partial differential equations with an order of 8 govern the plate kinematics (Reddy 2007). In the element formulation, the IP shear, IP bending, and axial deformation are coupled. This leads to differences in the curvature along the element for the FE and macromodels, in which the latter used separate (uncoupled) springs to model the shear, flexural, and axial behavior of the plates [off-diagonal components of matrix $K_{\mathrm{IP},\mathrm{Plate}}$ in Eq. (4)]. Furthermore, in the FE model, the rigidity of the element is distributed over the plate, whereas in the macromodel, the rigidity is concentrated in a specific number of discrete elements, i.e., the shear and tension–compression springs. Nonetheless, the average difference between the horizontal displacement in the FE and macromodels was 9.01%, which is an acceptable degree of accuracy. The horizontal displacement values for both the macro and FE models in the joint regions identified in Fig. 13(b) are provided in Table 3.

For the OP load case, both the TT joints and timber plates play a role in the load-transfer mechanism. Generally, joints are the weakest elements in IATP structures (Rezaei Rad et al. 2019a). Specifically, the ultimate strength of IMAs is lower than that of timber plates (Rezaei Rad et al. 2019b). Therefore, the IMAs control the ULS of the OP load case [Fig. 13(c)]. Both the macro- and FE models were subjected to the OP load corresponding to the ULS of the TT joints (Rezaei Rad et al. 2019a, b; Nguyen et al. 2019b). The OP loads applied at the TT joints are shown in Fig. 13(c). The mechanical properties of the two-node link elements used to simulate the TT joints are presented in Table 4.

The undeformed and deformed shapes of the macromodels, including the configurations defined in the sensitivity analysis, and the FE model are shown in Fig. 17. Similar to the IP load case, the behavior captured by the macromodel was similar to that of the FE model. The average difference between the displacement in the FE and macromodels was 8.3%. The displacement values for both the macro and FE models in six different joint regions [Figs. 13(c) and 17] are provided in Table 5.

The macromodel was verified against the results of a three-point bending test performed on a $1:1$ IATP subassembly, which was carried out at IBOIS. Three replicates of the prototype were tested by Nguyen et al. (2019a). The full-scale prototype, which represented a portion of an IATP structure, included 15 boxes. LVDTs were installed at the midspan of the assembly to record the displacements. The structural dimension, test setup, and the load–deformation curves are shown in Fig. 18.

Portions of the macro and FE models are shown in Fig. 19. The mesh size and element type noted in the description of the FE and macromodel section was adopted for the subassembly verification. For the macromodels, the section size of the inner beam elements and the stiffness values for the uniaxial tension–compression and shear springs are shown in Fig. 14.

The macro and FE models were subjected to the same loading protocol used in the physical test. Using the deflection-based criteria (span length/300), the responses were assessed at the serviceability limit state (SLS) to ensure a linear elastic response. Two LVDTs were installed on the north and south sides of the physical specimen (Fig. 18), and a force transducer was attached to the load actuator (hydraulic jack). Accordingly, the load–deformation curves were obtained from the macro and FE models and the physical experiments (including the average response from the three replicates). Figs. 20(a and b) show the load–displacement curves associated with the LVDTs installed on the north and south parts of the physical specimen. The results demonstrate that the responses obtained from both the macro and FE models were very similar to the experimental load–deflection behavior. The average linear stiffnesses computed from the load–deflection curve in Fig. 20 for the macro ($K_{\mathrm{Macro}}=2.95\mathrm{kN}/\mathrm{mm}$) and the FE ($K_{\mathrm{FE}}=3.64\mathrm{kN}/\mathrm{mm}$) models were similar in value to that of the physical specimen ($K_{\mathrm{Test}}=3.37\mathrm{kN}/\mathrm{mm}$). The stiffness values of the macro and FE models were $-12.4\%$ and $+8.0\%$ different, respectively, from the experimental stiffness value.

Considering the joints, plate types ($T_i$, $B_i$, ${\mathrm{CL}}_i$, and ${\mathrm{CT}}_i$), and the fiber orientation (fiber-parallel and fiber-perpendicular), 13 different element groups existed in the model domain (6 for the inner beams, 6 for the uniaxial springs, and 1 for the joints). To understand the contribution of the inner beams to the overall response, one element from each group and the associated local coordinate system is shown in Fig. 21(a). The local forces for each element group, including axial and shear and torsion and bending moments, were first calculated. Given the section properties (area and moment of inertia), the associated stresses were computed. Fig. 21(b) compares the absolute maximum stress for each element group. The overall response was governed primarily by the IMAs, and the role of the timber plates was less important (Fig. 21). Furthermore, among the DOFs associated with the response of the IMAs, the IP behavior (DOFs 1 and 2 in Fig. 21) governed the response of the full-scale prototype. The stress associated with these DOFs was 4.4 and 7.49 MPa, respectively. This state of stress is consistent with what was reported by Rezaei Rad et al. (2019a, b), and Nguyen et al. (2019b) for the SLS. In addition, the top and bottom timber plates ($T_i$ and $B_i$) had a greater influence on the response compared to the cross plates (${\mathrm{CL}}_i$ and ${\mathrm{CT}}_i$). This was directly related to the assembly process used for IATPs, in which each box is connected to its neighbor via the TT joints of the top and bottom plates. For DOF 6, which corresponds to the OP bending of the inner beams, the $T_i$ and $B_i$ plates had the greatest influence on the deformation, with stress values of 1.36 and 2.88 MPa for the fiber-parallel and fiber-perpendicular directions, respectively. In other words, the top and bottom plates were primarily responsible for resisting the OP loads. Regarding the cross plates, the ${\mathrm{CL}}_i$ plate was also subjected to the OP loads, with a stress of 1.65 MPa [DOF 6 of the fiber-perpendicular inner beam, CL-FPe, in Fig. 21(b)]. This is attributed to the fact that the four-sided box is geometrically asymmetric, and therefore, the cross plates also carry OP loads.

In addition to comparing the load-displacement response of the macro and FE models with that of the experiment, the time required to construct the model (data creation runtime) and run the analysis (analysis runtime) was evaluated. As stated previously, the TT joints govern the overall response of the system, and therefore, the FE model increases the computational complexity without improving the efficiency. In addition to the analysis runtime, the creation of the model database is very time consuming for the FE models. For the mesh size and configuration used in the verification studies, the model creation and analysis runtimes were 4 min 10 s and 8 min 36 s, respectively using an Intel Core i7-4800MQ CPU @ 2.7GHz with 16 GB of RAM. In contrast, the macromodel was constructed and analyzed in only 34 s.