Geometrically nonlinear behaviour of actively twisted and bent plywood

The current work aims to investigate the geometrically nonlinear behaviour of flexible plywood strips that undergo large displacements during the formation process of a structural element. Such deformations allow bending-active structures to reach curved geometries from initially planar elements. As our literature review points out, firstly, twist for plywood materials is not well studied, and secondly, the detailed determination of the stress state after the formation process of such structures, is usually oversimplified or even ignored. To this end, a finite element model using laminated shell theory is employed and calibrated against experimental data retrieved from physical prototypes, by paying particular attention in the accurate tuning of its boundary conditions. Furthermore, destructive tensile tests of twisted strips were implemented and their ultimate strength was compared with the numerical model. Our aim is to explore the material’s stress state during the strip’s 90 ◦ twist coupled with bending in two different directions, which is crucial for understanding the strip’s further behaviour as a structural element. In addition, we conduct parametric studies to showcase the non-negligible impact of geometric aspects on stress distribution along the strip. The results show that the utilisation ratios, using the Tsai–Hill failure criterion in terms of strength, for the case of the constructed beam element, are quite uniform throughout the strip and below 50%, apart from local stress concentrations. Finally, the parametric studies highlight the need to comprehend the material response to inform geometry within a structurally sound environment in order to unlock new possibilities for using thin plywood strips in lightweight structural and architectural applications, as demonstrated by the built beam elements.


Introduction
Bending-active structures, which reach their curved geometries through the deformation of initially planar elements, have recently been attracting growing interest from both architects and engineers [1].Such structures, that undergo large displacements, within the geometrical nonlinear spectrum during their formation processes, can be placed within the larger historical trajectory of studying the nonlinear aspects of elasticity.
For almost three and a half centuries, beginning with the experimental derivation of Hooke's Law [2], and after going through milestones such as Euler's elastica [3], the nonlinear aspects of the theory of elasticity have continually excited the scientific community, see for example the works: ''Elasticity and Geometry'' by Audoly and Pomeau [4] and ''W.T. Koiter's Elastic Stability of Solids and Structures'' edited by van der Heijden [5].The nonlinearities within elasticity can be On the other hand, application-wise, utilising the flexibility of the materials in the geometrically nonlinear range, to create curved shapes from flat elements, dates back even further and has a long history in vernacular architecture.Examples like Mongolian Yurts and Mudhif houses demonstrate the effectiveness of this approach by shaping the locally available and non-stiff materials, like straw or reed, to reach structurally sound, functional and inhabitable structures [16].Those structures combine spatial ideas and making, informed by material abilities, often resulting in lightweight, strong, functional, and aesthetically pleasing structures.Prominent examples that utilise elastic bending from the 20th century include Vladimir Schukov's lattice Towers [17], Buckminster Fuller's PlyDome from thin plates [18], and the elastic grid-shell of Mannheim multihalle by Frei Otto [19].Such architectural typologies, in the field of lightweight structures, have been distinguished by their intricate interplay of form, material behaviour, and forces.Their inherent need for an integrated design framework has called into question the conventional sequence in architectural design practice [20], where form, structure, and materialization follow a hierarchical order [21].
Despite the fact that the term ''bending-active'' typically refers to systems where bending is the main driver for elastically shaping a structure, torsion is occasionally regarded as falling within the scope of bending-active structures [22].It is not uncommon to observe the emergence of torsion as an outcome in elastic gridshells [23], prompting investigations into various form-finding and erection methodologies [24][25][26].Yet, the deliberate use of twisting as an input to initiate deformational behaviour is a less prevalent practice.For instance, the Focaccia Pavilion [27] employs flexible plywood strips, relying on twisting as the primary deformation mechanism to establish its final form.Schling et al. [28] utilise twisting of the slender steel lamellas to follow asymptotic curve networks on surfaces.In another case, the Timberfabric project achieves a braiding pattern by twisting plywood strips [29].On the other hand, the sun shading system of Softhouse [30] shows the potential of using elastic torsion for adaptive facade applications by allowing individual strips to responsively trace the sun through twisting.Slabbinck et al. [31] explored the possibilities of using elastic torsion in bending-active tensile structures, focusing primarily on its geometric stiffening effect.The Radical wood pavilion [32], designed by an Aalto University team, incorporates double-layered twisted plywood components, which demonstrates another application of twisting at an architectural component scale.Finally, the cantilevering beams of the Zero Gravity pavilion, designed and built by the authors, employ twisting as a main design driver where it serves to lock four thin adjacent plywood strips into one component [33].
In the early 1900s, the versatility of plywood for achieving curvilinear forms was acknowledged across a range of applications, including car bodies, canoes, and aircraft fuselages [36,37].Within the domain of building construction, the development of curved stressed skin panels leveraged the inherent flexibility of thin plywood sheets.By bending and adhering, thin sheets to a curved sub-structure formed curved surfaces, which may also serve as structural skin [38].Perkins [39] highlights the practicality of curving thin plywood plates to achieve tubular forms, referencing examples like grain silos during World War II.Expanding on these insights, Fuller's PlyDome, constructed in 1957, introduced a doubly-curved surface by bending the corners of individual plywood plates, today recognized as bending-active plate structure [22].Suitable materials for bending-active structures should possess the capacity to undergo large deformations without experiencing failure.Therefore, they require a high enough flexural strength-tostiffness ratio [40].Plywood is one of the most commonly used materials for such structures, not only due to its mechanical properties, but also due to its availability, ease of processing, renewability, and economic reasons [41].Over the past decade, diverse configurations of thin plywood plates and strips, featuring varying numbers of plies, have been used, contingent upon the specific application.While in some projects, like the Ongreening Pavilion [42], the Bend 9 project [22], the ReciPlyDome [43], and Geodesic winding [44], birch plywood was used in standard market sizes, the AA/ETH Pavilion [45] opted for custom-produced larger size spruce plywood panels.Furthermore, the ICD-ITKE Research Pavilion 2015/16 [46] utilised different grain orientations by creating custom-laminated beech plywood plates to achieve different bending stiffness across its components.
On the one hand, limited studies utilise large torsional deformations as a design driver as it has been reported in the literature, in less than 25% of the reviewed cases.On the other hand, plywood stands out as a prevalent choice for bending-active structures.Identifying the residual stress state in materials resulting from the formation process is particularly important for bending-active elements subjected to large deformations since it will remain as a residual stress state, in their further applications.Yet, in computational simulations, the investigation into residual stresses, due to formation process, in plywood is often limited, with exceptions in a few studies (e.g.[47]).To the best of authors' knowledge, no study applied laminated shell theory for analysing the behaviour of plywood strips within bending-active structures.The laminated shell theory can consistently incorporate both in-plane and bending stiffnesses, given the layer structure and their elastic properties.The classical laminated shell theory is often deemed sufficient for thin shells.Additionally, it facilitates in-depth investigations of stress states within individual layers for further analyses, such as predicting first-ply failure [48].
In our study, which lies at the crossroad between structural engineering and architecture, we use geometrically nonlinear finite element analysis to investigate the strength of a plywood strip, where large rotational and translational displacements were deliberately imposed as design drivers of its final twisted geometry.To investigate the residual stress distributions and overall behaviour, in detail, we simulate the plywood strip as a laminated shell, considering each veneer layer (ply), and we calibrate the model against the experimental data.In this regard, a beam element, previously employed in two different pavilion applications, see Fig. 1, is used as a case study.The beam element is composed of four identical twisted birch plywood strips.For our investigations, one of its four identical strip's deformation is simulated.In Section 2, the geometry of the beam element is described along with its full-scale applications.Section 3 describes the computational approach to simulate the strips.Section 4, discusses the results in terms of the element's stress state, along with parametric studies of varying: (i) surface veneer angle, (ii) the number of strip's plies, and (iii) and strip's width-to-length ratio.Additionally, in the same section, the uncertainties of the whole process are discussed, while Section 5 presents the conclusions of our study and the possible future steps.

Beam element from 4 twisted plywood strips
As the main objective of the current work is to study a hollowsection beam element that is actively formed by four flexible birch plywood strips, subsequent sections discuss its formation, its geometric digitization and finally its full-scale applications.

Physical prototypes
The assembly process of the beam element is initiated by forming a square tube from four individual thin-walled, 5-ply, 6.5 mm thick birch plywood strips.At mid-span, each of the four strips is screwed to a 20 mm thick square wooden block, which is used as distance keeper for practical reasons, see Fig. 2.1.In this way, the position of the strips can be adjusted relative to each other in the longitudinal direction, and inaccuracies in the manufacturing process due to the alignment of the strips can be prevented, while the connection between the strips is more easily facilitated.Two additional wooden blocks, as distance keepers, are inserted at quarters of the length of the beam element.Nevertheless, even without the aforementioned blocks, the final shape of the beam element would stay the same, see Fig. 2.2 to 5. As a Fig. 1.Full-scale applications using the beam element: (i) Zero Gravity Pavilion exhibition [34] and (ii) frame configuration from 20 × 4 Twist exhibition [35] (right), both built with 6.5 mm thick birch plywood strips.

Table 1
Geometrical properties of single birch plywood strip used for the fabrication of the beam element.

Length 𝐿 [mm]
Width  [mm] Thickness ℎ [mm] No. plies 3050 290 6.5 5 next step, the strips are interconnected loosely through lacing along a series of pre-drilled holes, to generate four common seamlines, see Fig. 2.2.At the same time, the two free ends of the tube are forced to close, so that the quadrilateral mechanisms at the ends of the hollow tube, transform from squares through rhombuses to straight seamlines, which are perpendicular to each other, see Fig. 2.2 to 5. Furthermore, the central axis of the initial tube (see Fig. 2.1) maintains its original position as the central axis of the final self-equilibrium geometry.At both ends of the beam, 2 mm steel plates are inserted in between the short edges of the plywood strips, see Fig. 2.3 to 5. Bolting the strips with the steel plate, locks the final shape of the beam element, see Fig. 2.4.The distances from the edge of the strip to the first and second rows of bolts are respectively 30 mm and 60 mm.During the closing process of the short edges, the already laced seamlines are adjusted by tensioning.Starting lacing before clamping the short edges makes the assembly process simpler, by limiting the individual spatial deformation tendency of individual strips and maintaining the reciprocity among them.The geometrical characteristics of the birch plywood we used for the strips are given in Table 1.For our application, the grain direction of the face veneer was chosen perpendicular to the longitudinal edge of the strips ( = 90 • ).In this way, the grain direction of only two out of the five plies is aligned with the longitudinal axis of the strip, see Fig. 3.This allows for decreased resistance to bending and torsion during the (de-)formation process.
During the described formation process, each strip deforms under (i) a twist of 90 • around x-axis (torsion), (ii) displacement equal to the strip's half-width along y-axis (bending around the strong zaxis) and (iii) displacement equal to the strip's half-width along z-axis (bending around the weak y-axis), see Fig. 4. Furthermore, during the (de-)formation process, the linear edges of rectangular strips transform into planar curves on either the xz-plane (vertical) or xy-plane (horizontal), see Fig. 3.While in the early stage scaled models made from thin cardboard, crease lines perfectly define planar curves on such planes [33], in beam elements made of plywood, the seamlines connecting the strips tend to separate near the end of elements, see Fig. 5 (Top).This phenomenon, observed due to the collision of material thickness, is incorporated into the calibration of the simulated model, later described in Section 3.

Digitizing the geometry: photogrammetry
The final geometry of the beam element is digitized through photogrammetric reconstruction, which will be used for the calibration of the numerical simulation.To this end, raw images are captured with a DSLR camera at regular spatial intervals around the beam element.The collected images are then processed in RealityCapture software to generate 3D point cloud consisting of 13.7 million points and mesh with variable densities.The conducted photogrammetric digitization process is based on the methods described by [49] and is explained in more detail by the authors in a previous work [50].

Full-scale applications
The beam elements were employed in two distinct architectural applications: initially as cantilever beams within the roof of a kinematic, lightweight pavilion, namely ''Zero Gravity'' [51], and subsequently as the main components of a lightweight kit-of-parts structure, namely ''20 × 4 Twist'' [35], see Fig. 1.Both structures are developed by a geometry-driven structural design approach, which necessitates the disciplinary knowledge of architectural design and engineering.The resulting pavilions are lightweight, easy to assemble, disassemble and reassemble into new configurations to meet the changing demands of their users in the future [52].
Even though the width-to-length ratio of the used strips for making the beam elements was consistent in both applications, namely 1/10, the lengths of the strips differed.Specifically, the ''Zero Gravity'' application utilized a strip length of 2000 mm, whereas the ''20×4 Twist'' had a strip length of 3050 mm.Nevertheless, in both cases, the beams were assembled from planar, 5-ply birch plywood strips with face veneer grain direction perpendicular to the longitudinal edge of the strips.

Numerical model
The numerical model aims to accurately reproduce the deformation of the strips in the beam element during its formation process, so that the residual stresses (i.e. the remaining stresses before further application) can be realistically estimated.In the model, only a single strip is simulated due to symmetry.The rest of the assembly is taken into account by appropriate boundary conditions, imposed by the adjacent plywood strips.To this end, the strip's two long edges are restrained to horizontal and vertical planes, as shown in Fig. 3, which correspond to mirrored deformations in the adjacent strips during the formation process of the beam element.However, as noted in Section 2, closing the gap between the strips in the middle of the element generates separation between the strips near the ends of the beam, see Fig. 5.To allow this observed behaviour in the real beam element, constraints in the corresponding plane, either vertical or horizontal, are partially released along the strip's long edges   , see Fig. 3.A 90 • rotation is introduced to the short edge of the strip, while the other edge is constrained against rotation.It is, however, not possible to capture the local shape of the twisted strip near the end realistically, by applying simple fixed or pinned conditions along the short edges.Therefore, in the simulation, the end conditions are introduced by using additional connection plates that are coupled with the surfaces of the strip by contact.This approach was essential to capture locally the shape at the ends of the strips, and by extension, to provide a realistic estimation of the stress state, which is one of the main objectives of the current work.The stiffnesses were calibrated against the photogrammetric data, so that the best matching between numerical simulation and the actual beam is achieved.Furthermore, in the beam element, the longitudinal elongations of the strips are not restricted.In the modelling approach, this is achieved by using only a small shear stiffness between the connection plates and the strip.
All numerical simulations presented herein have been performed as finite element (FE) analyses with Abaqus software [53].For verification of the FEM results derived from Abaqus, please refer to Appendix, where these results are benchmarked against a well-established nonlinear plate problem.An overview of the simulation model, including the mesh, the boundary conditions, and the strip-connection plate contact detail, is presented in Fig. 6.All simulations were performed as quasi-static, geometrically nonlinear analysis using the software's large deformation analysis option.Both the plywood strip and the connection plates have been modelled using 4-node conventional shell elements with full integration (S4).The restraints on the long edges have been applied by fixing the nodes in y-and z-directions on the bottom and top edges, respectively, except for the released parts of the edges.In the fixed end of the strip, the connection plate has been fixed along the x-direction and against all rotations, while in the loaded end, the displacement has been fixed along the x-direction and the rotations about the y-and z-axes.Note that ends are not fixed in the y-or zdirection, since this would introduce over-constraints and the current boundary conditions are sufficient to reach the required deformation.The load itself is applied as prescribed rotation about the x-axis up to  x = 90 • .
The section behaviour for the strip is defined via composite shell layup, which is used to define the section in terms of individual layers with corresponding thicknesses, material properties and layer orientations.The sanding process applied to the plywood panels during their fabrication is incorporated into the analysis.In line with this, the face veneers are attributed a thickness of 1.15 mm, whereas the inner veneers are set at 1.4 mm. Simpson's rule with three integration points per layer is used for calculating the section stiffnesses.For each layer, the material is defined as linear elastic, orthotropic, while the connection plates are modelled using homogeneous shell section with linear elastic isotropic material.The contacts between the connection plates and the strip are defined by surface-to-surface contact formulation with combination of normal contact and cohesive behaviour.This allows different normal stiffness under separation (cohesive behaviour) and compression (normal behaviour), while shear stiffness is always cohesive.In summary, in the resulting contact formulation, the interface normal stress and interface shear stresses,  N ,  T,1 and  T,2 , respectively, can be expressed in terms of relative interface displacements as where  N and  T, are the relative normal and tangential displacements in the interface, and  N (  N ) and  T are the normal and tangential interface stiffnesses respectively.To delineate the local interface coordinate system, refer to Fig. 6.The normal stiffness depends on the sign of the relative normal displacement such that where  N,t is the tensile stiffness and  N,c is the compressive stiffness.
Utilization of the plywood, in terms of strength, is evaluated by applying failure criterion for each veneer layer individually.Throughout the years, different multi-axial phenomenological failure criteria have been developed to predict the behaviour of wood and wood products [54].Among them, the Tsai-Hill failure criterion [55], which is interactive, between different stress components, quadratic criterion and is an extension of the von Mises-Hencky distortion energy hypothesis [56], is widely used for anisotropic composite materials.More specifically, Oh [57] studied the tensile strength of laminated veneer lumber (LVL) of various grain directions, under different failure criteria, and showed that Tsai-Hill criterion predicted accurately the strength for 0 • and 90 • , while underestimated the strength for angles 15 • , 30 • and 45 • .Pramreiter et al. [58] studied the tensile strength of Finnish birch veneers and showed that Tsai-Hill performed accurately.Wang et al. [59] showed that Tsai-Hill criterion performed well for the where  11 is the stress parallel to the grain,  22 is the stress perpendicular to the grain and  12 is the in-plane shear stress, and  1 ,  2 and  v,12 are the corresponding strength parameters.To account for the differences between tensile and compressive strengths of the material, the criterion is defined separately for each axial stress quadrant by defining the axial strength components as where subscripts t and c refer to tension and compression, respectively.

Model parameters
The elastic stiffness and strength properties of birch veneers are based on the information provided by the manufacturer (UPM) [60] and they are summarized in Table 2.In wood materials, there is natural variation between the material properties of individual samples.To this end, the elastic properties provided by the manufacturer, represent mean values.However, the strength properties were characteristic values (5%-quantiles), which were used to estimate corresponding mean values by multiplying with a factor of 1.17.This multiplication factor corresponds to the ratio between the mean and characteristic values as stated by the manufacturer.For the connection plates, typical properties of steel, elastic modulus  = 210 GPa and Poisson's ratio  = 0.3, were used.
The interface properties for the contacts were chosen such that (i) experimentally observed curvatures at the end of the strip are closely captured and (ii) axial deformations are not restrained, which is also the case in the physical prototype.The first point is achieved by setting compressive stiffness  N,c to a high value ( N,c = 1000 N,t , virtually allowing no penetration in compression) and tensile stiffness  N,t was calibrated based on observed strip shape (see Section 3.3).To allow almost free axial deformations of the strip, interface shear stiffness was set to a low value of  T = 0.25 MPa/m.

Model calibration
The calibration of the numerical model was conducted by aligning the scanned data retrieved from physical prototypes through photogrammetry and the simulated geometry in Rhinoceros 3D [61] environment.The two main parameters in question were the edge release length  rel and the normal contact stiffness  N,t , which were identified  through the calibration process, see Fig. 3.Moreover, the edge release length was estimated by measuring the opening lengths (separation between strips, see Fig. 5) from the beam elements and the length  rel = 0.15 L, corresponding to 457 mm, was chosen for the model, which is in good accordance with the measurements.The contact stiffness was calibrated by multiple rounds of comparison of the simulated geometry and the beam's photogrammetric reconstruction.The range of contact stiffnesses that was explored was  N,t = [10, 50, 250, 1250, 6250] MPa/m and finally the best fit was obtained with  N,t = 50 MPa/m, which was used for the subsequent simulations.Considering the possible uncertainties in the photogrammetric model, and its alignment with the simulation, investigating a larger set of stiffnesses was not meaningful.
As a result, the calibrated FEM model shows a good agreement with the scanned geometry of the beam prototype, as shown in the overlay between the numerical model and photogrammetric reconstruction, see Fig. 5 and in the resulting generated deviation map, see Fig. 7.The discrepancy in deviation values observed between the region close to the strip's vertical edge (approx. 1 mm) and the region close to the strip's horizontal edge (approx.4 mm) is attributed to the imperfections (tipping down) of the physical prototype, as outlined in Section 4.3.More specifically, the horizontal edge of the strip in the physical model is tipping down, and therefore the calibration is conducted based on its vertical edge.For the same reason, in Fig. 5, the cutting planes for sections of the beam element are located starting from the beam's vertical edge.

Preliminary strength verification through destructive tensile test
Tensile destructive tests were performed to verify the simulated strength of the twisted plywood strips.Fourteen samples were cut from 5-ply birch plywood panels, which belonged to the same batch used for the twisted beam elements, all maintaining a face grain angle of 90 • .Before testing, the strips were conditioned in a climate-controlled room set to a temperature of 20 • C and a relative humidity of 65%.
Before applying tension, the universal testing machine's upper head was pre-adjusted to a 90 • rotation, relative to its lower head, to facilitate testing of the twisted strips, with the following dimensions:  = 50 mm, length  = 940 mm, and thickness ℎ = 6.5 mm.The strips were clamped at both ends (Fig. 8 left), resulting in a clear test length of  = 844 mm.A constant displacement of 2 mm per minute was applied until the sample failure.All samples demonstrated failure at   due to cracks that aligned with the direction of the face veneers, see Fig. 8 right, while notably, no strips exhibited failure within the grip zone.
The numerical simulation followed the procedure outlined in Section 4.1.In order for the ultimate load   to be determined, the applied tensile load to the twisted strip is incremented in steps, and the stresses are computed for each ply, while the ultimate strength is identified when the Tsai-Hill failure criterion is fulfilled.The ultimate load   was extracted from the region exhibiting constant stress distribution.Fig. 9 presents the ultimate load result derived from the numerical simulation alongside the full dataset of ultimate loads obtained from the destructive tests.The mean value  = 11.1 kN of the physical tests is also highlighted, accompanied by a standard deviation of 2.2 kN.Accordingly, while   = 9.73 kN closely aligns with the experimental test results, it slightly underestimates the material's ultimate strength.The simulation indicates that shear failure precedes tensile failure, which might account for the observed underestimation, as elaborated by Wang et al. [59] regarding the shear strength of birch plywood.However, more extensive tests and related discussions are outside the scope of this study.

Stress analysis
The deformed state of the strip and the residual stresses in the plywood are of particular importance in our study, as they provide a starting point for the structural understanding of the aforementioned beam element.As a first step towards this direction, the three individual displacements, namely rotation 90 • , displacement along the weak and strong axis equal to  ∕2, are applied to the strip as explained in Section 2, and the preliminary results from the envelope of the utilisation ratios using the Tsai-Hill failure criterion are presented in Fig. 10.Inhere, and in all the following results, the term ''envelope'' denotes the value across the thickness of the strip.The results show a wide range of variation in the utilisation ratios between the case of bending in the strong axis and the other two cases.Nevertheless, the actual formation process of the element's strips, conceptually resembles the case of lateral torsional buckling of the strip bent around its strong axis.As it can be seen in Fig. 11 the deformation process of the beam's actual strip, does not lead to high utilisation ratios in comparison to the strong axis bending.The overall result from the strip's simulation presented in Fig. 11, shows an almost uniform utilisation ratio along the strip apart from local stress concentrations near the ends of the strips.While the utilisation ratios of the strip remain at around 30% to 40% throughout its length, it reaches a maximum of around 70% (point B in Fig. 11), which is most probably caused by the specific selection of the boundary conditions, and around 50% at its tip (point A in Fig. 11) of the material's capacity in local areas.
Besides assessing the overall strength of the strip through the Tsai-Hill failure criterion, stress values  11 parallel to the grain direction,  22 perpendicular to the grain direction, and the  12 representing the in-plane shear stresses were extracted, from the mid-plane of each individual ply, see Fig. 12.While ply 5 represents the outermost layer of the plywood, ply 1 represents the innermost layer, as illustrated in Fig. 3.It is worth noting, that the stress directions vary for each ply due to the changing grain direction.For example, in ply 5, the grain direction is parallel to the short edge of the strip, and so is the extracted stress  11 .On the other hand, in ply 4, the grain direction is parallel to the long edge of the strip, and so is the extracted stress  11 .The results show that for the case of parallel to grain stresses ( 11 ), plies 2 and 4 carry the maximum compressive and tensile stresses, see Fig. 12 left.More specifically, the maximum compressive stresses  11 occur in the central part of the plies 2 and 4 with a maximum compressive stress to strength ratio of the order of 33%, see Section 3.2.On the other hand, the maximum tensile stresses  11 occur in the longitudinal edges of the plies 2 and 4 with maximum tensile stress to strength ratio of the order of 68%.For the case of plies 1 and 5 the maximum compressive stresses occur at the ends of the strip with a maximum stress to strength ratio of the order of 18%, while the maximum tensile stress to strength ratio at the same region is of the order of 15%, see Section 3.2.In the case of stresses  22 perpendicular to the grain direction, similar stress patterns for plies 1, 3, and 5 are observed with the stresses  11 for plies 2 and 4, see Fig. 12 middle and left.The maximum stresses  22 occur in the inner and outer plies, namely plies 1 and 5. On the one hand, the maximum compressive stresses  22 occur in ply 5 and the two ends of the strip with maximum stress to strength ratio of the order of 15%, while on the other hand the maximum tensile stresses  22 occur in the same region in ply 1 with maximum stress to strength ratio of the order of 23%, see Section 3.2.Finally, for the case of in-plane shear stresses  12 the maximum positive shear stresses occur in inner ply 1 at the two ends of the strip with maximum stress to strength ratio of the order of 40%, while the maximum negative shear stresses occur at the ends of the strip in the outer ply 5 with maximum stress to strength ratio of the order of 38%, see Fig. 12 right and Section 3.2.Additionally, Fig. 13 displays the stress distributions across the thickness of the plywood, examined from the location where local stress concentrations are observed, at point A in Fig. 11.Fig. 10.Envelopes of Tsai-Hill utilisation ratios using from geometrically nonlinear FE analysis of clamped plywood strip, from left to right: (i) twisted 90 • around x-axis, (ii) displaced half of its width along y-axis (bending around the strong z-axis) and (iii) displaced half of its width along z-axis (bending around the weak y-axis).

Parametric studies
In the following sections, three independent parametric studies are presented by varying: (1) face grain angle, (2) number of plies and (3) strip's length-to-width ratio.It is worth mentioning that the model was parameterised by using Python scripting in Abaqus.

Parametric study 1: varying the face grain angle
In the first parametric study, which is inspired by observations made in nature [62][63][64], we evaluate the dependency between the face grain direction angle and the strength of the beam element's strip, by varying the face grain angles ranging from 0 • to 90 • .We assume firstly, that 0 • orientation is parallel and 90 • orientation is perpendicular to the strip's longitudinal edge, secondly, that the geometrical characteristics of the strip, namely width, length, and thickness, along with the number of plies are kept constant and thirdly, that the grain direction between two consecutive plies always remains equal to 90 • .The envelopes of the utilisation ratios are presented for face grain angles:  = [0 • , 30 • , 45 • , 60 • , 90 • ] in Fig. 14.The results show that if the grain orientation of three ( = 0 • ) instead of two ( = 90 • ) plies, would have been along the longitudinal direction of the strip, larger regions with higher utilisation ratios would have occurred at the end of the strips, but still with the same order of magnitude, namely 64%, see Fig. 14.On the other hand, for the cases of face grain angle  = [30 • , 45 • , 60 • ], the maximum utilisation ratios appear to be lower compared to  = [0 • , 90 • ], namely of the order of 60%, but with larger regions with higher ratios than 40%.Furthermore, a more detailed observation can be made from Fig. 15 for the case of face grain angle  = 45 • , where the stresses  11 ,  22 and  12 for the individual plies, are presented.The maximum tensile stresses parallel to the grain direction  11 are observed in plies 2 and 5 with a maximum stress to strength ratio equal to 38%, while the maximum compressive stresses are observed in ply 4 at the ends of the strip with maximum stress to strength ratio equal to 37%, see Fig. 15 left.For the case of the stress  22 perpendicular to the grain direction, the maximum tensile stresses are observed in ply 1 at both ends of the veneer with maximum stress to strength ratio equal to 16%, while the maximum compressive stresses are observed at both ends of ply 5 with maximum stress to strength ratio equal to 15%, see Fig. 15 middle.The maximum negative shear stresses  12 occur at the longitudinal edges of plies 1, 3, and 5, while the maximum positive shear stresses occur at the longitudinal edges of plies 2 and 4, both with max stress to strength ratio equal to 45%, see Fig. 15 right.Finally, the graph in Fig. 16 depicts how the maximum utilisation ratio varies with changes in the face grain angle .The upper line represents the trend on the tip of the strip, including the critical stress concentrations at point A in Fig. 11, while the lower line shows the trend in the central third section of the strip.

Parametric study 2: varying number of plies
The alteration in the behaviour of the strips with different thicknesses was investigated by modelling them as 3-ply, 5-ply, and 7-ply laminates, having respective nominal thicknesses of 4 mm, 6.5 mm, and 9 mm.The comparison range for the number of plies is determined according to the market availability of Finnish birch plywood.The envelopes of the utilisation ratios of the strips are presented in Fig. 17.The results show a slight increase in the utilisation ratios using the Tsai-Hill failure criterion, by increasing the number of plies from 5 to 7. For the case of 3 plies a triangulated pattern of utilisation is observed with uniform values and lower compared to the other two cases, see Fig. 17.Similar triangulated patterns resembling cone-like shapes have also been observed in the case of twisted ribbons, see [65].

Parametric study 3: varying strip's length-to-width ratio
Before the fabrication process of the beam elements, we conducted a series of physical experiments to determine the most suitable size of plywood strips for deformation.Plywood strips with small lengthto-width (∕ ) ratios, as expected, exhibit higher resistance to twist, whereas strips with high ∕ ratios, were not resistant enough after twisting.Therefore, an optimum ratio for building the prototypes was chosen, which was close to ∕ = 10, as it was the most practical.The insights gained from these experiments determined the range for our investigations in the digital environment.In order to be realistic, we prefer to present three different ∕ ratios, namely 8, 10, and 12, for stress investigations, by keeping the width of the strips constant while adjusting the length accordingly.Furthermore, the release length is kept proportional to the length of the strips such that  rel = 0.15.The results of the utilisation ratios are presented in Fig. 18.The results, as expected, show higher utilisation ratios for ∕ = 8 and lower ratios for ∕ = 12.More specifically, the maximum ratios for the case of ∕ = 8 reaches 94% at the longitudinal edges of the strip.Furthermore, there is a reduction of 20% between the case ∕ = 10 and ∕ = 12, which would result in an element with lower residual stresses.

Uncertainties and sensitivity
There are a number of uncertainties that propagate through the overall process and can lead to discrepancies between the model output and the response of the physical system of interest.In this section, a brief summary of the identified uncertainties is given.The sources of uncertainties are subdivided into ones related to: (i) manufacturing  physical prototype, (ii) digitizing and calibration process, and (iii) simulation model.

Physical prototypes
The final geometry of the physical prototypes is influenced heavily by imperfections that arise during the assembly process, as they are handmade.Differences in the by-hand-applied tension force during the lacing process of the neighbouring strips can result in variations (i) in the seamlines and (ii) in the opening (separation between strips) near the closing part of each element, see Fig. 5 top right.Even though the calibration of the FEM model was based on the scanned element, the separation between the strips, was varying among the physical prototypes.Additionally, due to a varying angle between two neighbouring strips in their common seamline, sliding between the strips can locally occur [50], which is prevented by the inserted wooden blocks within the beam.Finally, the overall geometry of the beam element is affected by the variations in the assembly setup, which are used to close its short edges (see Fig. 2).These variations can cause mismatches between the strips, as it is observable in Fig. 7, where the short edge of the strip on the right side of the image is pressed more than it should, causing the tip of the beam to curve downwards.

Photogrammetry and digitizing physical data
During the digitizing process of the physical prototype, a specific set-up was prepared to hold the beam element from its both ends only, which allows to capture the images all around it.Additionally, spotlights were located around the beam to prevent inconsistencies that can occur due to shadows.To test the accuracy, scale bars with control points are distributed on the floor and the distances between them are used as references during the reconstruction process, as described [49], and the average error was 0.063 mm.However, the rope that is used to lace the strips partially, covers the connection between them, which causes local discontinuities on the long edge of the strips in the digitized model.On the other hand, the overlapping of the meshes retrieved from photogrammetry and the simulation was conducted manually.However, considering the inherent imperfections associated with the manual production of beam elements, it is assumed that the comparison conducted is sufficient for calibrating the digital model.

Simulation
There are numerous sources of uncertainties in the numerical model, which can be categorized into three groups: (i) idealizations related to geometry and boundary conditions, (ii) effects of finite element discretization, and (iii) material behaviour.
A notable geometry-related idealisation is that the strip is treated as a zero-thickness shell, due to which planar restraints (see Section 3.1), are applied to the mid-plane of the strip instead of the edges of the strip, as in the physical prototype, which can lead to discrepancies in the final geometry.In addition, the selection of the boundary conditions (release length   and contact stiffness   ) can lead to local stress concentrations, see point B in Fig. 11.
It is well accepted that finite element approximation leads to errors that are dependent on the mesh density, [66].Here, the element size has been chosen as a compromise between the computational cost and the accuracy (see mesh in Fig. 6).By performed tests, refining the element size to 1/3 of the current size increases the Tsai-Hill utilisation  ratio only by 3% in the critical section near the end of the strip, see point A in Fig. 11.Therefore, it is not considered a significant source of error.However, there is an additional stress concentration at the start of the edge releases (point B in Fig. 11), which is mesh-dependent (utilization at this point increases significantly with refining mesh), and the actual location of the critical section should be checked carefully when evaluating the simulation results.
The choice of the failure criterion, namely Tsai-Hill, has an effect on the utilisation ratios.In this regard, an important uncertainty relates to the fact that Abaqus built-in Tsai-Hill failure criterion does not account for out-of-plane shear stresses  13 (rolling shear) and  23 .Although the simulated maximum rolling shear stresses are of the order of  23 = 0.85 MPa for the strip, they do not appear in the Tsai-Hill criterion, see Eq. ( 3).
It should be also noted that wood, apart from other uncertainties regarding its material properties, has a tendency to creep under prolonged load [67][68][69][70], which in the case of the beam element studied, would lead to stress relaxation.Since the effects of creep are not taken into account in the analysis, the real residual stresses in the beam element are expected to be lower than the ones predicted by the model.

Discussions and conclusions
Our studies have shown that plywood strips' behaviour in deformed states, as a result of undergoing large and complex deformations, including coupled twisting and bending, is strongly influenced by changes in their initial geometry.Therefore, informed design processes are essential to utilise the material's capacity effectively.Overall, by following an integral design framework, we can understand and further manipulate/tailor the structural performance of a highly deformed strip through variations in its initial geometry and material lay-up.This reciprocity between architectural and structural limits and qualities, expands the design and application space of thin plywood strips in lightweight structural and architectural applications, as demonstrated by the built pavilions.
A detailed FE simulation approach for assessing the behaviour of the beam element consisting of four interconnected strips has been presented.The parameters of the numerical simulation were calibrated against an experimental prototype.The following points are worth highlighting: (i) the importance of investigating the stress state after the formation process of bending-active plywood elements with detailed FE simulations is overseen in the literature, (ii) the lack of investigations reported in the literature for twisted bending-active structures, and (iii) the importance of calibrating the numerical model with experimental data (photogrammetry).Based on the numerical results from the numerical simulations the following conclusions can be drawn: 1. Fairly uniform utilisation ratios, based on the Tsai-Hill failure criterion, of 40% for the beam element, apart from local stress concentrations that lead to utilisation ratios of 70% in point B and of 50% at point A in Fig. 11, has been observed.In terms of individual veneers, the maximum stress to strength ratio was 68%, for the case of tension parallel to the grain direction, which shows that thin birch plywood is sufficient material for large displacements.In addition, the residual stresses observed during

Table 3
Geometric dimensions of the annular ring alongside its associated material properties used for the numerical simulation.The results are presented without employing specific units, consistent with the reference case [71].

Appendix. Numerical verification: annular plate strip under edge shear force
The FEM results obtained from Abaqus software have been validated through a comparison with the results presented by Arciniega and Reddy [71], focusing on a well-established benchmark problem concerning an annular plate strip's nonlinear deformation.The choice of this particular benchmark for validation is due to its relevance and similarity to the deformation studied in this work.Notably, the benchmark problem encompasses both bending and torsion and is modelled as a multi-layered composite, offering a coherent foundation for comparative analysis.
In the problem setup, the geometrical dimensions of the plate along with the mesh size have been modelled in accordance with the reference case, specified in Fig. 19.The model employed S4, 4-node conventional shell elements with full integration, consistent with the modelling used in the previously discussed analysis of the plywood strip.It is noteworthy that the referenced benchmark employed higher-order elements.Four laminate configurations are considered: [90 • /0  ], with respect to polar coordinates.In this context, the radial direction corresponds to the 0 • orientation.Each laminate layer is presumed to maintain uniform thickness, specifically, 0.01 for the 3-ply (h/3) and 0.0075 for the 4-ply (h/4) configurations.The plate's dimensions and the material properties of the lamina are given in Table 3.The annular plate with a slit, depicted in Fig. 19, has one edge labelled 'CD' that is clamped, while the opposing edge, labelled 'AB', is subjected to a uniformly distributed load of  = 0.45.For each laminate configuration, the load-deflection curve at point B is extracted and subsequently compared with the benchmarked results, see Fig.

Fig. 3 .Fig. 4 .
Fig. 3. Simulation notations on the strip.: length of the strip,   : length of the released edge,  : width of the strip,   : width of the clamped region and : plywood face grain angle.Blue marked edges show the constrained parts within the vertical and horizontal planes.The supplementary figure provides a visual representation of the veneer orientations, starting from the outermost ply 5 and progressing inward to ply 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .Fig. 6 .
Fig. 5. Overlay of low-poly mesh from photogrammetric reconstruction and the simulated plywood strip.(Top) Top view with the detail of separation behaviour of the strips, (middle) right view from the overlay, (bottom) Cross-sectional deviations between the simulated plywood strip (in red) and physical prototype (in black).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Spatial deviations between the simulated plywood strip and the physical prototype, based on photogrammetric data.

Fig. 9 .
Fig. 9.The full dataset of ultimate loads   from the destructive tests, highlighting their mean value , in comparison to the results derived from the Abaqus simulations   .

Fig. 11 .
Fig. 11.Envelope of Tsai-Hill utilisation ratio for the simulated plywood strip with the face grain angle perpendicular (90 • ) to the longitudinal edge of the strip as employed in the prototypes, see Fig.5.

Fig. 12 .
Fig. 12. Simulated stresses from the mid-plane of each ply for the case of face grain angle perpendicular (90 • ) to the longitudinal edge of the strip as constructed in the prototypes, see Figs. 5, 11.Ply 5 denotes the outermost layer and ply 1 denotes the innermost layer.

Fig. 13 .
Fig.13.The stress distributions across the thickness of the plywood from the critical stress concentration in the strip (point A in Fig.11).

Fig. 15 .
Fig. 15.Simulated stresses from the mid-plane of each ply for the case of surface grain angle 45 • to the longitudinal edge of the strip.Ply 5 denotes the outermost layer and ply 1 denotes the innermost layer.

Fig. 16 .
Fig. 16.Maximum utilisation ratio of the plywood strip in relation to its surface grain angle.

Fig. 18 .
Fig. 18.Envelope of Tsai-Hill utilisation ratios of the simulated plywood strip with three different ∕ ratios (parametric study 3).

Fig. 20 .
Fig. 20.Vertical displacement at point B versus applied force  = 4 in different laminate configurations.