Numerical Simulation of the Semi-Rigid Behaviour of Integrally Attached Timber Folded Surface Structures

Abstract

Timber folded surface structures assembled using semi-rigid multiple tab and slot joints (MTSJ) have been shown to form feasible structural systems with high load bearing potential. However, for their further development and use on large building scales, a pertinent model for prediction of their structural behaviour has yet to be developed. This paper focuses on simplified numerical methods for accurately modelling the semi-rigid structural behaviour of bidirectional timber folded surface structures with multiple tab and slot connections. Within this scope, the structure behaviour is considered to be in the elastic stage. Three practical methods of analysis for such structural systems are presented. The first two approaches use the Finite Element Method (FEM), where the theory of plates and shells are applied. In the first method, the MTSJs are modeled using strip element models, while, in the second strategy, spring models are used. The third modeling strategy elaborates on the new macroscopic mechanical models, referred to as macro models. Sets of one-dimensional (1D) elements are used to represent the mechanical behaviour of the entire system. Both linear and geometric nonlinear analysis are performed for all three modeling strategies. The numerical results are then validated against the large scale experiments. Comparison of the strip and spring element model results have shown that the strips represent more accurately the experimentally obtained values. Concerning the macro modelling approach, very good agreement with both detailed FE modelling approaches, as well as experimental results, were obtained. The results indicate that both linear and nonlinear analysis can be used for modelling the displacements within the elastic range. However, it is essential to include geometric nonlinearities in the analysis for accurate modelling of occurring strains as well as for displacements when considering higher load levels. Finally, it is demonstrated that including semi-rigidity in the numerical models is of high importance for analysing the behaviour of timber folded surface structures with MTSJ.

1. Introduction

Timber folded surface structures assembled using integral mechanical attachments, specifically multiple tab and slot joints (MTSJ), have been shown to form feasible structural systems with high load bearing potential [1]. However, for efficiently realising such structures on building scales, it is essential to provide a pertinent method for their structural analysis. The large scale tests performed by Stitic et al. [2] provided insight into timber folded structure behaviour, but only for a certain loading condition and a specific geometry, i.e., vertical loading and anti-prismatic geometry. For different forms as well as varying load combinations, the load paths and load resisting mechanisms will not necessarily be the same. Therefore, developing a pertinent model for prediction of their structural behaviour is crucial for the efficient use of timber folded surface structures on a building scale.

For structural analysis of complex folded surface structures, the commonly used theory for mathematical model definition is the theory of plates and shells. However, due to the complexities in the geometry, boundary conditions and heterogeneous material properties of timber, obtaining an analytical solution for such mathematical model is very difficult [3,4]. For these reasons, the solution is usually obtained using numerical modelling, mainly based on finite element methods (FEM). Previously developed numerical models for timber folded surface structures have shown that joint semi-rigidity representation has an important influence on the global structural behaviour [5,6,7]. The joints were either modelled as completely rigid or as hinges, leading to highly inaccurate estimations of structure displacements. The values of occurring plate strains were not addressed. Furthermore, experimental work performed on MTSJ connections [2,8,9] also demonstrated the importance of including semi-rigidity in the numerical models. However, defining the MTSJ connections semi-rigidity by modelling their complete complex geometry has proved to be too cumbersome and computationally expensive to be used on larger scale models with multiple panels [10]. Examples of simplified numerical models of connection semi-rigidity were found in [11] where joint strips of different stiffness parameters were used for connection representation in glass plate shell structures, and in [12,13] where springs were used for modelling finger joints in a segmental plate shell structure.

This paper studies two simplified FE models for MTSJ connection simulation within timber folded surface structures: (1) strip element model and (2) spring model. Additionally, a macro modelling approach is studied where one-dimensional (1D) elements are used to represent the mechanical behaviour of the entire system. Such a modelling strategy has been previously used to simulate the progressive collapse of reinforced concrete shear walls [14,15] as well as assessing the performance of timber framed structures [16]. Panto et al. [17] introduced a discrete modelling approach to analyse historical masonry structures where both in-plane and out-of-plane properties of masonry walls were taken into account. In this paper, we are inspired from the model proposed by [17] and adapt it to the mechanical, fabrication, and construction constraints of integrally-attached timber plate structures.

The test results used for validating the ABAQUS${}^{\mathrm{TM}}$-based FE models as well as the proposed macro model are taken from large scale timber folded surface structures tests presented in [2].

The paper is structured as follows: in Section 2, large scale tests performed on MTSJ closed slot folded surface structures are presented together with results relevant for numerical model validation; Section 3 presents the description of the numerical models, two proposed simplifications for the connection detail modelling and the timber folded surface structure macro model, as well as the obtained results; Section 4 includes the comparison of the results with the ones obtained from experimental testing as well as discussion on the found differences; finally, the main conclusions are summarized in Section 5.

2. Experimental Investigations

Large scale tests on folded surface structures with MTSJ closed slot connections are presented shortly hereafter. The material used was Kerto-Q structural grade Laminated Veneer Lumber (LVL) panels of 21 mm thickness. Figure 1 presents the global geometry of the structures, taken as a regular anti-prism based folded form with a single constant curvature.

Such curvature was chosen in order to have all plate elements of equal size and in this way simplify the load application setup. As uniform loading was simulated by applying load at discrete points at each relevant panel geometrical center, equal size panels enabled these loads to be of the same size as well. MTSJ closed slot detail parameters are shown in Figure 2.

Load application was done in a quasi-static rate using a specifically devised setup which enables simultaneous, continuous loading of the structures’ top surface. Boundary conditions that allow rotation about the y-axis were used along the two supporting sides. Force transducers, linear variable differential transformers (LVDTs) and inclinometers were placed as shown in Figure 3. Additionally, a three-dimensional digital image correlation (DIC) system was used for obtaining strain and displacement fields of the entire structure. The DIC system was calibrated for each test individually in order to ensure the accuracy of measured values. Due to the sensitivity of the DIC measurements to the light source and random speckle pattern quality, as well as difficulties related to obtaining reliable results near the boundaries and sharp edges, strain gauges were used for validating the DIC results at selected points. Data acquired from the tests were analysed using both VIC 3D^TM provided by Correlated Solutions (Irmo, SC, USA), and custom algorithms developed within Matlab^® (version R2017a, MathWorks, Inc., USA). For more exhaustive details on the performed tests, the reader is referred to [2].

Results

The total load vs. midspan deflection curves of tested large scale structures are shown in Figure 4a. A group of three experimental replicates was produced and each replicate was fitted with single curve. The results suggest that, when considering the proposed timber folded surface structures for building scale, the governing limit state in their design will be the serviceability one (SLS). The results further show that, for a span of 2.9 m, the SLS maximal allowed displacement, equal to 9.66 mm ($L/300$ according to [18]), stays well within the elastic stage for all tested structures. For this reason, only linear elastic material behaviour is considered here and in further presented numerical models. The representative structure stiffness in the elastic domain, k, was determined by fitting a linear regression model to the three replicates within the elastic region of their load–displacement curves (see Figure 4b). The upper bound of the elastic region of each replicate was determined by imposing $R^2>0.99$, where $R^2$ is the coefficient of determination of the linear regression. The displacements were extracted from the DIC measurements at the point marked with x (Figure 5), as the maximum deflection found around the loading ring.

Displacement field of the entire structure obtained from the DIC system is shown in Figure 5 (left), with respect to the structure global coordinate system (see Figure 1). Strain fields are shown for one selected central panel only, as the global system coordinates need to be transformed into local plate coordinates to obtain correct values of strain for each individual panel (see Figure 5 right). The replicate chosen for representing the results of an entire group is the one with most unfavorable results in terms of stiffness.

Validation of DIC obtained results was performed by using strain gauges for measuring strain at selected points. Low elastic strain gauges of LFLA-10-11 type and 10 mm length were used as recommended for wood based materials [19]. The position of strain gauges is shown in Figure 5 (right). The strain data obtained from DIC measurements was examined at the same positions and directions as those of respective strain gauges. Within Vic-3D^TM software, the image matching was performed on a subset size of 23 × 23 pixels with the correlation analysis performing at every 5th pixel within the area of interest (step size $=5$). With respect to the set step size the filter size was set to nine data points, in order to represent a virtual gauge with a physical smoothing length of ∼10 mm. Furthermore, the strain values were calculated with respect to the local plate coordinate system, with z-axis corresponding to the plate normal direction. For validating the DIC results by comparing them to the strain gauges data, a noise reduction algorithm using a moving average filter was applied for reducing the noise influence and enabling a more straightforward comparison [20]. The strain gauge and DIC obtained results that showed good agreement (see Figure 6), where DIC strain noise levels were within the typical prescribed resolution of $\pm 100$ microstrain [21].

3. Numerical Modelling

In order to create a pertinent numerical model for analysing the structural behaviour of timber folded surface structures, a mathematical model consisting of three fundamental parts is used: (1) the geometrical description of the structure; (2) specification of material properties; and (3) mathematical expressions of the basic physical laws which govern the structural behaviour.

The geometrical description of the structure form was identical to the one used in laboratory tests and was directly imported into the numerical analysis software from the CAD software used for its initial generation.

For modeling the chosen structural grade LVL panel material and its behaviour, an orthotropic model is used. Kerto-Q 21 mm thick panels consist of seven 3 mm thick spruce peeled-veneer laminates from which one fifth is glued crosswise in a lay-up $\mid -\mid \mid \mid -\mid$ (Figure 7). This kind of structure improves the lateral bending strength and stiffness of the panel. More importantly, in this way, very homogenous and mechanically strong panels are obtained which can be assumed as having orthotropic material properties [22]. The used characteristic values of elastic properties for Kerto-Q LVL material principal directions are shown in

Table 1. Elastic properties used for orthotropic material model of Kerto-Q 21 mm panels [22].

Property	Symbol	Value	Unit
Elastic modulus 11	$E_{11}$	10,000	N/mm${}^2$
Elastic modulus 22	$E_{22}$	2400	N/mm${}^2$
Elastic modulus 33	$E_{33}$	130	N/mm${}^2$
Shear modulus 12	$G_{12}$	600	N/mm${}^2$
Shear modulus 13	$G_{13}$	60	N/mm${}^2$
Shear modulus 23	$G_{23}$	22	N/mm${}^2$
Poisson’s ratio 12	${\nu }_{12}$	0.09
Poisson’s ratio 13	${\nu }_{13}$	0.85
Poisson’s ratio 23	${\nu }_{23}$	0.68
Density	$\rho$	510	kg/m${}^3$

. These moduli further define the elastic compliance matrix according to Equation (1).

$$\begin{array}{c}\hfill \begin{array}{c}\left(\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{33}\\ 2{\epsilon }_{12}\\ 2{\epsilon }_{13}\\ 2{\epsilon }_{23}\end{array}\right)=\left(\begin{array}{cccccc}{\displaystyle \frac{1}{E_{11}}}& -{\displaystyle \frac{{\nu }_{21}}{E_{22}}}& -{\displaystyle \frac{{\nu }_{31}}{E_{33}}}& 0& 0& 0\\ -{\displaystyle \frac{{\nu }_{12}}{E_{11}}}& {\displaystyle \frac{1}{E_{22}}}& -{\displaystyle \frac{{\nu }_{32}}{E_{33}}}& 0& 0& 0\\ -{\displaystyle \frac{{\nu }_{13}}{E_{11}}}& -{\displaystyle \frac{{\nu }_{23}}{E_{22}}}& {\displaystyle \frac{1}{E_{33}}}& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{1}{G_{12}}}& 0& 0\\ 0& 0& 0& 0& {\displaystyle \frac{1}{G_{13}}}& 0\\ 0& 0& 0& 0& 0& {\displaystyle \frac{1}{G_{23}}}\end{array}\right),\left(\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{12}\\ {\sigma }_{13}\\ {\sigma }_{23}\end{array}\right)\hfill \\ \mathrm{where}\frac{{\nu }_{ij}}{E_{ii}}=\frac{{\nu }_{ji}}{E_{jj}}\mathrm{for}1\le i,j\le 3,i\ne j.\hfill \end{array}\end{array}$$

(1)

The regarded timber folded surface structures are composed of discrete planar elements of continuous nature, where the plate thickness dimension is considerably smaller compared to its other dimensions. Therefore, the plates can be considered as “thin” (thickness/average side ratio: $t/L\le 0.05$ [3]). For this reason, appropriate “thin” plate theories are used for the mathematical model definition. The applied theory choice further defines the element type used in FE numerical solution procedures. For the presented problem, ABAQUS provides an element type called conventional shell element, where the element geometry is specified at the reference surface which is coincident with the shell’s mid surface (Figure 8). The desired thickness of the panels is defined by section property. In this case, a three-dimensional continuum is approximated with a two-dimensional theory. This reduction in dimensionality is achieved by using the fact that the element thickness is considered small when compared with its other characteristic dimensions [23]. For plates, FE formulation quadrilateral 4-node, six degrees of freedom, large-strain shell elements of type S4R with reduced integration were used. This is a robust, general-purpose quadrilateral element that is suitable for a wide range of applications [23]. Such elements normally use thick shell, Reissner–Mindlin plate theory, however, as the shell thickness decreases and the transverse shear deformation becomes very small, they converge to classical thin shell theory. The elements become discrete Kirchhoff thin shell elements where transverse shear strains are neglected. A detailed study on element type choice is presented in Appendix A. Both linear and geometric nonlinear analysis were performed, i.e., material nonlinearity is not considered. A geometrically nonlinear, updated Lagrangian shell formulation was employed for reviewing the influence of rigid body rotations and translations as well as to account for coupling between axial and bending action within the plates. A mesh density study was performed to ensure that reliable results are obtained (see Appendix B). The mesh was refined until the two consecutive solutions were found to be in good agreement and the results were converging. The total computation time was also taken into consideration. A mesh density corresponding to a seed size of 0.0075 m was found to be sufficiently refined and was used in all further presented models.

Since the plates were modelled as shell elements also the overall behaviour of the joint needed to be modelled on a two-dimensional level. In this regard, two different possibilities for detail modelling were explored:

• Introducing a connecting strip element between the panels with characteristics which represent the stiffness of the connection detail (Figure 8b); see Section 3.1.
• Assigning a set of springs along the adjacent edges with appropriate stiffness parameters (Figure 8a); see Section 3.2.

In addition to the detailed modeling of folded plate structures using conventional FE technique, a new modeling strategy is proposed where only 1D mechanical elements are used. Spring elements with rigid frame boundaries simulate the kinematic degree of freedom (DOF) of both thin shells and boundary conditions of folded surfaces; see Section 3.3.

Detailed experimental work performed on MTSJ connections [2,8,9], provided valuable input of bending and shear stiffness used for the simplified numerical modelling of the connections within the scope of this paper.

For realizing adequate boundary conditions within the experimental setup, heavy steel and LVL plates were used for fixing the structure (∼33 kg at each support side, for comparison the total structure weight was ∼65 kg) (see Figure 9). As the supports were positioned under an angle of ${27}^{\circ }$ in order to follow the curvature of the structure, their center of mass was not in line with the center of rotation of the steel rollers (see Figure 10a). This resulted in nonlinear secondary moment effects which needed to be taken into account in the FE model. In Figure 11, the order of magnitude of support weight influence on rotations is explained in a simple example of an equivalent size singly curved shell with the same support conditions as in the used models. Finally, the supports were represented with equivalent single plate of the same weight and the center of rotation determined with respect to the actual setup (Figure 10b). The modelled support plates were defined as rigid bodies, where only rigid body motion is allowed. Their rotation was constrained to the center of rotation, $RP$. The boundary conditions allowing rotation about the structure y-axis were prescribed at the center of rotations at the reference points $R{P}_1$ and $R{P}_2$, for left and right support side, respectively.

Furthermore, the degrees of freedom of the nodes along the plate edges fixed to the supports were constrained by using Multi-point constraint (MPC) definition within ABAQUS. As the local coordinate system of each plate differs from the globally defined one, the rotation of the above mentioned edges and their constituting nodes needed to be set in a global coordinate system to properly simulate the experiment conditions. Therefore, MPC type BEAM was used to provide a rigid beam between the edge nodes and the respective reference point. In this way, the nodal displacements and rotations were constrained to the ones of the reference point, corresponding to the presence of a rigid beam in between [23]. In the FE model, the load was applied in the same manner as in the actual test setup. In the first step, the gravity load was imposed on the entire model. In the second step, the 70 mm diameter circular surfaces in the plates geometrical center, obtained by partitioning the plates, were loaded by surface traction, the same as loading rings in the tested structures (see Figure 12).

It is noteworthy to mention that, due to complex geometry and a large number of plate elements, for defining material properties, sections, varying local plate material axes as well as individual spring and strip definitions, the entire model generation was done automatically using Python programming language and the ABAQUS Scripting Interface. This allowed for fast changes in element type, mesh size as well as spring and strip properties, for exploring different modelling techniques while avoiding tedious manual work.

3.1. Strip Element Model

In the first approach for modelling the connection details, a small gap of 17 mm was introduced between the connecting plate elements in order to represent the existing joints. The overall dimensions of the structure were kept the same, but the size of individual plate elements was reduced by $5\%$, with the origin point set at the geometrical center of each scaled plate respectively. The gap was then filled by adding a strip element for representing the connection between the neighbouring faces (Figure 12). The thickness and material characteristics of the connecting strip were assigned according to the stiffness of the actual connection detail used. This way of modelling joint details was used as in [11], where isotropic joint strips of different stiffness parameters were used for connection representation in glass plate shell structures and have shown satisfying results. Secondly, the 70 mm diameter circular surfaces in the plates geometrical center, obtained by partitioning the plates, were loaded by surface traction, the same as the loading rings in the tested structures (see Figure 12).

In the presented FEM, rotational, axial and shear stiffness were used for defining the strip material behaviour. The MTSJ connections representing strips were modelled as linear elastic and orthotropic. As opposed to [11] where the joint material was modelled as isotropic, in the presented model, the orthotropy needed to be introduced due to large over estimations of shear modulus parameters when an isotropic material stiffness matrix was used. The values for rotational and shear stiffness per strip unit length, $k_m$ and $k_{v,12}$ were taken based on small scale test results, presented in [2,9]. As MTSJ connection details are in fact an integral part of the plates themselves, the representing shear stiffness in the remaining two directions, $k_{v,13}$ and $k_{v,23}$, as well as axial stiffness, $k_n$, were calculated using the generalized Hooke’s law with the use of LVL material elastic properties (see Table 1). Concerning the axial stiffness, $k_n$, since the strips’ and the plates’ local coordinate systems do not coincide (see Figure 8a), the known plate elastic properties had to be rotated form the plate local coordinate system, where the z-axis corresponds to the panel normal, to the strip local coordinate system, where the z-axis corresponds to the respective strip normal orientation. The latter was performed by using Hankinson’s equation, Equation (2) [24], where $\theta$ was equal to $32.5^{\circ }$

$$E_{\theta }=\frac{E_{22}{E}_{33}}{E_{22}{sin}^2\theta +E_{33}{cos}^2\theta }$$

(2)

Furthermore, the rotational, axial and shear stiffness for determining the strip material properties are given by Equations (3)–(5), respectively, which were derived by assuming plane stress distribution and zero value for the Poisson ratio [11]:

$$k_m=\frac{1}{6}\frac{h^3}{w}E,$$

(3)

$$k_n=\frac{2h}{w}E,$$

(4)

$$k_{v,ij}=\left\{\begin{array}{c}{\displaystyle \frac{h}{w}{G}_{ij}\mathrm{for}i,j\in \{1,2\},i\ne j,}\hfill \\ {\displaystyle \frac{w}{h}{G}_{ij}\mathrm{for}i,j\in \{1,3\}\vee \{2,3\},i\ne j.}\hfill \end{array}\right.$$

(5)

Having a predefined strip width, w = 17 mm, and, according to Equations (3)–(5), the height, Young and shear modulus parameters for the orthotropic material strip are given by Equations (6)–(8), respectively:

$$h=\sqrt{12\frac{k_m}{k_n}},$$

(6)

$$E=k_m\frac{6w}{h^3},$$

(7)

$$G_{ij}=\left\{\begin{array}{c}{\displaystyle k_{v,ij}\frac{w}{h}\mathrm{for}i,j\in \{1,2\},i\ne j,}\hfill \\ {\displaystyle k_{v,ij}\frac{h}{w}\mathrm{for}i,j\in \{1,3\}\vee \{2,3\},i\ne j.}\hfill \end{array}\right.$$

(8)

The complete set of parameters used for modelling the strips are given in

Table 2. Strip element model, estimated stiffness parameters per strip unit length.

Property	Symbol	Value	Unit
Rotational stiffness	$k_m$	1684	N
Axial stiffness	$k_n$	126.5	N/mm${}^2$
Strip width	w	17	mm
Strip thickness	h	12.64	mm
Elastic modulus	E	85.13	N/mm${}^2$
Shear stiffness 12	$k_{v,12}$	31.35	N/mm${}^2$
Shear stiffness 13	$k_{v,13}$	108.76	N/mm${}^2$
Shear stiffness 23	$k_{v,23}$	39.88	N/mm${}^2$
Shear modulus 12	$G_{12}$	42.17	N/mm${}^2$
Shear modulus 13	$G_{13}$	80.85	N/mm${}^2$
Shear modulus 23	$G_{23}$	29.65	N/mm${}^2$

. The strip offset from the edge ends, marked as $l_e$ in (Figure 8a), was equal to the offset of the first tab, where the connection between the two plates is established. Within the MTSJ connections design process, the length of this unconnected edge part was limited to a maximum $10\%$ of the total edge length. Accordingly, the strip dimensions were modeled to correspond to the realized plate connection length.

3.2. Spring Model

In the second approach, the semi-rigidity of the joints is introduced by means of springs along adjacent edges. Similar to the strip model, each plate surface was reduced by 1%, introducing a small gap (∼3 mm) between the plates. This was needed in order to be able to define the direction in which they act. Each tab of the MTSJ was modelled by a group of sets containing three spring types, representing the rotational, axial and shear stiffness of the connection (Figure 8b). Appropriate springs stiffness values were determined according to the mechanical and material properties of the actual MTSJ, the same as in the strip model. However, their values were taken as such and divided only by the number of spring connectors of the respective spring type within each individual group (see

Table 3. Spring model, estimated stiffness parameters per tab of the MTSJ closed slot connection. The stated values are divided by the number of spring connectors of each spring type used within one group.

Property	Symbol	Value	Unit
Rotational stiffness	$k_m$	168,449.52	Nmm
Axial stiffness	$k_n$	12,656.63	N/mm
Shear stiffness 12	$k_{v,12}$	3135	N/mm

). It is important to note that, in the presented model, all springs were considered as being uncoupled.

Furthermore, assigning multiple sets within one group used for modelling the behaviour of each MTSJ tab was assessed for obtaining improved results. First, a group with one set of spring connectors was examined. The obtained results showed that such modelling generates an overly flexible system leading to overestimated vertical displacements. A higher number of sets consisting of three types of spring connectors, with their respective stiffness reduced according to the total number of sets used within a group, provided more accurate results. As shown in Figure 13, increasing the number of spring sets causes the obtained results to converge. A good agreement was found in using seven sets per group for representing each individual tab of MTSJ (as shown in Figure 12 right).

3.3. Macro Modelling Approach

In this approach, the plates are modelled by series of spring elements which try to simulate the behaviour of a finite portion of timber plates in tension/compression, bending, out-of-plane (OOP) shear and in-plane (IP) shear. Macro model representing a portion of a timber plate is illustrated in Figure 14. A comprehensive elaboration of the macro modeling strategy can be found in [17,25,26]. According to the orthotropic homogeneous properties of LVL timber plates, the mechanical characteristics of the macro model are calibrated using fibre discretization strategy using the properties provided by the manufacturer. In addition to the plates, the joints are modelled using the spring model, as discussed in Section 3.2.

The tension/compression axial stiffness of the macro model, $K_{axi}$, is shown in Figure 14a. Since two rows of axial springs are defined, the bending behavior of timber plates can be included in the kinematics of timber plates. Timber orthotropic properties are considered in the model by separately calibrating the interfaces according to the fibre orientation and corresponding mechanical properties mechanical properties. The axial stiffness, $K_{axi}$, is given by Equation (9):

$$K_{axi}=2\frac{EA}{L},$$

(9)

where E is the Young modulus (see Table 2), A is the tributary area and $L/2$ is the tributary length (see Figure 14a).

The transversal springs are introduced to control the out-of-plane response of the macro models. The elastic shear stiffness, $K_{OOP,shr}$ is computed using Equation (10):

$$K_{OOP,shr}=\frac{G_{OOP}A}{2L},$$

(10)

where $G_{OOP}$ is the out-of-plane shear modulus ($G_{13}$, see Table 2), $A/2$ is the tributary area and L is the tributary length (see Figure 14b).

Finally, the diagonal springs are calibrated using the in-plane shear properties of timber portion. Equations (11) and (12) show the IP shear stiffness:

$$K_{IP,shr1}=\frac{G_{IP,shr1}Lt}{h},$$

(11)

$$K_{IP,shr2}=\frac{G_{IP,shr2}ht}{L},$$

(12)

where $G_{IP,shr1}$ and $G_{IP,shr2}$ are the shear moduli of timber along the perpendicular and parallel direction of timber fibres corresponding to $G_{12}$ and $G_{23}$, respectively (see Table 2), L is the length of the timber portion, t is the thickness of the plate, and h is the height (see Figure 14c,d).

The size and number of the macro models are affected by the degree of robustness and confidence. The macro models can have different sizes in order to capture the precise response of the structure. Figure 15 shows two macro models having different mesh sizes. Finally, the macro models fare generated for each timber plate which are then assembled together. Figure 16 shows the macro model representation of the prototype. Macro models associated with the timber folded plates are modeled in an open source computational platform, OpenSees [27]. Two-node link elements are used to simulate the behavior of the already-defined springs.

4. Results and Discussion

In order to compare the numerical and experimental results, the numerically obtained data needed to be corrected for the values caused by the gravity load imposed in the first loading step. This was necessary because the measurement of the total load in the performed experiments was done by force sensors placed under the plates; therefore, the gravity load was not a part of the measured force. Additionally, all the measurements started to be recorded from the time when the structure was already put in position, meaning that all the initial effects of gravity load were also disregarded in the experimental results. On the other hand, in the experiments, both the displacement and the strain field were captured from above the structure with the DIC system, meaning that, in the subsequent loading steps, there existed a propagated influence of gravity on the resulting deflections and strains. For this reason, in the numerical model, the first loading step including gravity could not be simply excluded, as gravity propagated influence in the subsequent steps that needed to be simulated. Therefore, for correcting the results of the numerical model, the values obtained at each step were reduced for the initial gravity load influence of the first loading step, but the propagation of gravity load in subsequent steps was kept intact.

The elastic part of the total load vs. midspan deflection curve of the tested structures is compared to the numerically obtained results in Figure 17. Both strip and spring FE models as well as the marco model show very good agreement with the linear approximation of the three tested samples. On the other hand, the results of models where geometric nonlinearity was included seem to capture with better accuracy the slight tendency of experimental results to deviate from the linear approximation with increasing loads. In the case of the macro model, corotational geometric transformation technique is used to simulate the geometric nonlinearity. In both types of analysis, the strip shows slightly less stiff behaviour compared to the spring and macro model. This is attributed to its possibility to introduce additional properties for modelling the connections behaviour, such as shear stiffness in the remaining two directions, $k_{v,13}$ and $k_{v,23}$. When examining the displacement profile along two lines, as shown in Figure 18, the error of the spring model is more pronounced. The strip model and the macro model in both linear and nonlinear analysis represent the actual behaviour with more accuracy. The maximum absolute error with respect to the experimentally obtained results, amounted to $6.6\%$ for the spring, $2.2\%$ for the strip model, and $4.6\%$ for the macro model with linear analysis. Concerning the geometrically nonlinear analysis, the errors amounted to $4.8\%$ for the spring, $3\%$ for the strip model and $4\%$ for the macro model.

Furthermore, a comparison of the strain results show that not taking into consideration geometric nonlinearities gives highly erroneous values of strain, especially in some cases (see $S{G}_1$ in Figure 19). The stiffness matrix which is formulated for the undeformed geometry at the beginning of the linear analysis is not updated after each loading increment, as is in the case with geometrically nonlinear analysis, and cannot capture the change in strain values as the system deforms. Generally, the results show that the strains in the plates are fairly low in magnitude even for high total loads but have a strong nonlinear dependance on the load at certain positions, $S{G}_1$ and $S{G}_4$ (see Figure 19). However, even though the limit state design is governed by the values of displacements, for modelling accurate strain magnitude, including geometric nonlinearity seems to be essential. When comparing the strain results of strip and spring models, it can be seen that higher discrepancies are found along the skewed edges of the plates, strain gauges $S{G}_3$ and $S{G}_4$. This is due to the fact that the forces which appear at the straight edge connections are acting very much in the same directions as the defined bending, axial and shear springs within the spring model. However, on the skewed edges, where shear forces are more dominant, twisting of the mutually connected panels and shear which occurs in the other two in-plane directions, 13 and 23, are not taken into account with the assigned springs. The same modeling deficiency is also realized in the macro models since the models adopt the spring elements. This seems to be the main shortcoming of using the spring elements. As the plates in the models are represented by their mid surface, the springs can only be assigned to connect nodes on those mid surface planes, meaning that only one direction of shear in that specific plane can be represented. The numerically obtained results of the linear and geometrically nonlinear analysis show that the plates themselves undergo what can be considered small strain and deformations, i.e., the values of resulting engineering and logarithmic strains are equal, indicating that an important part of the structures displacements are in fact plate rigid motions as well as system rotations due to semi-rigid connections.

For examining the level of the MTSJ semi rigidity, a structure with completely rigid connections was modelled using high rigidity springs within the spring model (see Figure 20). At the load of 25 kN, when MTSJ structures reached the SLS maximal allowed displacement of 9.66 mm, the displacement of the rigid jointed structure shows to be 5.7 mm, only 60% of the maximal allowed value. The 9.66 mm displacement for the rigid model corresponds to the load of 42 kN, a 68% higher load then in the case of MTSJ structures, which clearly indicates that the semi rigidity indeed plays a decisive role in the structural behaviour of systems using closed slot MTSJ.

5. Conclusions

In this paper, three simplified methods for modelling the semi-rigid behaviour of timber folded surface structures with MTSJ connections were examined, namely, strip element, spring element and macro element modelling approaches. Based on the obtained results and observations, the conclusions are as follows:

• Compared to the spring, the strip element model presents more accurate results with respect to the experimentally obtained values. It is inferred that the main advantage of using strips is due to their possibility to include more connection properties into the model, such as shear stiffness $k_{v,13}$ and $k_{v,23}$. Respectively, the inability to do so was found to be the main shortcoming of the spring element model.
• The results showed that using the strip element seems to be a feasible way for modelling the MTSJ semi rigid behaviour within a folded plate system in its linear part. Both linear and nonlinear analysis could be used for this purpose as the displacements are a key parameter for modelling such structures. However, including geometric nonlinearities in the analysis proved to be essential for accurate modeling of occurring strains as well as for displacements when considering higher load levels.
• On the other hand, the potential advantages of using the spring element model could be found in their slightly more fast and straightforward generation. Total time needed for the automatic generation of the spring model takes only 2 min of total CPU time, while the strip element definition takes 4 min. This time includes geometry input, material and sections definition, assembly, step generation, defining the springs/strips and their local coordinate systems as well as specifying loads and boundary conditions (i.e., final analysis computation time was found to be approximately equal). When considering more complex geometries with possibly varying MTSJ geometrical parameters, for which the computation time would increase exponentially, using the spring elements could prove to be more efficient. Additionally, using spring elements would be more appropriate for geometry cases with nonparallel adjoining plate edges, as using strips would result in their inconsistent widths along a single edge (for example, timber plate-shell geometries as shown in [28]).
• The proposed macro modeling strategy takes the advantage of spring elements only and avoids the application of mesh-based FE elements. Taking the master DOFs that have the most influential contribution to the kinematics of such systems, the global stiffness matrix of the structure is considerably decreased and the computational overhead is incrementally reduced. While the accuracy remains in an acceptable level, total computational time takes less than 1 min for the prototype structure. Nevertheless, the main modeling drawback is the boundary conditions, where the interaction of DOFs cannot be realistically simulated. However, the macro models are proposed for fast-feedback analysis in practical works.

In conclusion, when comparing MTSJ structures to those with completely rigid connections, the presented results demonstrate that the MTSJ details have a significant influence on the structural behaviour of the global system. Therefore, it is suggested that rigorous models including the connection semi-rigidity are necessary for appropriately analysing the structural response of timber folded surface structures with closed slot MTSJ.