Issue 
Mechanics & Industry
Volume 19, Number 3, 2018



Article Number  301  
Number of page(s)  12  
DOI  https://doi.org/10.1051/meca/2018004  
Published online  03 September 2018 
Regular Article
Modeling the influence of connecting elements in wood products behavior: a numerical multiscale approach
^{1}
Université ParisEst, Laboratoire Modelisation et Simulation Multi Echelle,
MSME UMR 8208 CNRS, Champs sur Marne,
77454
MarnelaVallée cedex 2, France
^{2}
Institut technologique Forêt Cellulose BoisConstruction Ameublement, FCBA,
10 rue Galilée,
77420
Champs sur Marne, France
^{3}
ESIPE, Université ParisEst MarnelaVallée,
5 bd Descartes, Champs sur Marne,
77454
MarnelaVallée cedex 2, France
^{*} email: luc.chevalier@univparisest.fr
Received:
8
August
2017
Accepted:
10
January
2018
Wood furniture is often composed of simple parts that may be modeled as beams or plates. These particularities allow using simplified approaches that reduces the number of degrees of freedom (dof for short) in a finite element simulation of the furniture's behavior. Generally, connections are not taken into account in such simulations but these connections are critical in the failure process of the furniture and it worth studying it precisely. Using a multiscale approach, this paper introduces a numerical procedure to identify the connection contribution in the furniture's stiffness. Comparing 3D finite element calculations with a Timoshenko's beam calculation on a corner of two wooden parts, we identify the specific behavior of the connection elements (pins, nut, screw… and local 3D effects) to introduce it as a punctual 0D element in the beam code. Two validations of the approach are presented here: (i) a numerical validation by comparing the result of the beam code with a complete 3D finite element simulation on a representative plane structure of wooden furniture; (ii) an experimental validation by managing a global bending test and measuring the displacement field using digital image correlation (DIC for short).
Key words: Wood product / connection / multiscale approach / digital image identification
© AFM, EDP Sciences 2018
1 Introduction: from dispersion of wood properties to the necessity of a fast modeling tool
Bunk bed furniture is widely produced for domestic use. Since furniture may be described as a mechanical structure submitted to various loadings [1], it can be classified into three groups depending on their construction design [2]:

frame type structure constructed of beamlike elements of wood or wood based materials, such as chipboard, Medium Density Fiberboard (MDF for short) or plywood;

can type structure constructed of platelike elements of wood based materials;

combined structure constructed of both beam and plate elements.
The different components of each structure are assembled with a variety of joining techniques. These techniques are either simple with only one or two fixing components or more complex.
In service, bunk bed structure is subjected to static loads such as punctual static forces and dynamic loads, such as cyclic or impact loads. It is mostly for safety reasons, that the European committee of standardization (CEN) has issued a particular norm [3] for this type of furniture that includes validation tests for both dynamic and static loads.
Generally, the strength and the stiffness of wood furniture are highly influenced by three main factors [4]:

the type of wood or woodbased material out of which furniture is manufactured;

the assembling technique used to join the elements;

the type of load applied.
Wood highlights a complex behavior, heterogeneous, anisotropic and sensitive to hygrometric conditions and a classical orthotropic elastic model is often used to carry out finite element simulations. To take into account the uncertainties of the model in our calculations, it is necessary to characterize all parameters with probabilistic laws so it is possible to carry out numerous simulations. In that way we can obtain a statistical response for the structure: mean displacement or stress, standard deviation, confidence region as it is presented in [5]. Making thousands of simulation is possible in this context thanks to the geometrical particularity of the bunk bed structure composed of long elements where the use of the beam theory is possible.
Nevertheless, 3D effects appear in the junctions between beams and one must take these effects into account to make sure of the accuracy of the simulations. In a beam finite element code, junctions are 0D punctual element that can be interpreted as local springs. In this paper we study the behavior of a specific assembly (bolt − dowel nut − pins) used to fasten the main corners in a wood bunk bed frames under static loads. The multiscale approach developed is quite classical but here, we chose to run a purely numerical procedure to identify the behavior of the 0D element: results of 3D finite element calculations done using the SolidWorks software are compare with a Timoshenko's beam calculation on a corner of two wooden parts, this enables the identification of the specific behavior of the connection elements (pins, nut, screw… and local 3D effects). The choice of SolidWorks instead of another finite element code is imposed by the furniture manufactory context.
Two validations of the approach are presented: (i) a numerical validation by comparing the result of the beam code with a complete 3D finite element simulation on a representative plane structure of wooden furniture; (ii) an experimental validation by managing a global bending test and measuring the displacement field using digital image correlation.
In the next section (Sect. 2), we present the multiscale approach. We develop the stiffness matrix of a perfect corner under beam theory as well as the 3D results and explain the identification process and justify the assumption made (neglecting the elongation effects in regards of rotational effects in the 0D model). We finally conclude that rotational stiffnesses are not constant values and depend on the load intensity.
In Section 3, we compare the deformed shape of a bunk bed structure using the SolidWorks 3D code (we focus on the mean line of each beam) with the displacement field obtained with the beam code where we included the nonlinear behavior of the 0D junctions. This comparison is very good despite the complexity of the bed structure in regards of the simplicity of the corner. Finally, an experimental section is devoted to validate this approach.
2 A multiscale analysis of the jointing element influence
A throughbolt with dowelnut connection is studied in this section. This system is usually used both as primary connectors and also to reinforce weaker joints [6]. The throughbolt with a dowelnut connection is of particular interest because of its strength and reliability [7] thus it's widely used to connect bed rails to bed posts. Figure 1 introduces the geometrical characteristics of the corner studied to manage the multiscale approach. First we manage calculation on a perfect corner using the Timoshenko's beam theory and we highlight a first stiffness matrix between the load applied at section B and the displacements and rotations of the section B. Second, we use a 3D finite element code (SolidWorks) to manage 3D calculations on the 3D model of the corner including bolt and other wood pins used to fix the two parts together. The difference between the displacements components of the 3D calculations and the beam ones leads to the influence of the connection on the corner behavior. This is the global procedure of the multiscale approach used in the following sections.
Fig. 1 Screen shot of the CAD model of the corner. Lengths are h = 447 m and d = 428.5 mm. Sectional areas are 57 × 57 mm^{2} edge for the vertical part and 20 × 140 mm^{2} for the horizontal one. All dof are naught in section A and loads are applied in section B. 
2.1 Perfect corner via beam theory
Considering that the beams' behavior is elastic, one can use superposition of the six independent loading cases: forces F_{X}, F_{Y}, F_{Z} and moments M_{X}, M_{Y}, M_{Z} applied on the section at B. These six cases are representative of any loading case. $$\left\{{F}_{ext}\right\}={\left\{\begin{array}{c}\hfill {\displaystyle \overrightarrow{F}}={F}_{x}{\displaystyle \overrightarrow{x}}+{F}_{y}{\displaystyle \overrightarrow{y}}+{F}_{z}{\displaystyle \overrightarrow{z}}\hfill \\ \hfill {\displaystyle \overrightarrow{{M}_{B}}}={M}_{x}{\displaystyle \overrightarrow{x}}+{M}_{y}{\displaystyle \overrightarrow{y}}+{M}_{z}{\displaystyle \overrightarrow{z}}\hfill \end{array}\right\}}_{B}\text{.}$$(1) In the following, the vertical part is denoted (1) and the horizontal one is denoted (2). Young's modulus and shear modulus are noted E_{i} and G_{i} (spruce is used for the two parts and subscript ‘i’ is not necessary here) and quadratic momentum and cross section are respectively denoted I_{ix} _{or z} and S_{i} When we apply forces F_{X}, F_{Y} or bending moment M_{z} the displacements field stays inplane and only two displacements u_{B} and v_{B} and a rotation θ_{z} are non naught. For example if we apply a force F_{X}, the displacements and rotation are: $$\begin{array}{ccc}\hfill {u}_{B}=\left(\frac{{h}_{\mathrm{}}^{3}}{3{E}_{1x}{I}_{1z}}+\frac{h}{{G}_{1xy}{S}_{1}}\right){F}_{x},\hfill & \hfill {v}_{B}=\frac{{h}^{2}d}{2{E}_{1x}{I}_{1z}}{F}_{x},\hfill & \hfill {\theta}_{z}=\frac{{h}^{2}}{2{E}_{1x}{I}_{1z}}{F}_{x}\hfill \end{array}\text{.}$$(2)Because of the great difference between E (11 580 MPa) and G (540 MPa) the shear effect is as important as the bending effect. Compression or tension effects in longitudinal direction are neglected because their contribution leads to small displacements in regards of the bending and shear effect. Consequently, if one applies a F_{y} force: $$\begin{array}{ccc}\hfill {u}_{B}=\frac{{h}^{2}d}{2{E}_{1x}{I}_{1z}}{F}_{y},\hfill & \hfill {v}_{B}=\left(\frac{h{d}^{2}}{{E}_{1x}{I}_{1z}}+\frac{{d}_{\mathrm{}}^{3}}{3{E}_{2x}{I}_{2z}}+\frac{d}{{G}_{2xy}{S}_{2}}\right){F}_{y},\hfill & \hfill {\theta}_{z}=\left(\frac{hd}{{E}_{1x}{I}_{1z}}+\frac{{d}^{2}}{2{E}_{2x}{I}_{2z}}\right){F}_{y}\hfill \end{array}\text{.}$$(3)For the bending moment M_{z} we have: $$\begin{array}{ccc}\hfill {u}_{B}=\frac{{h}^{2}}{2{E}_{1x}{I}_{1z}}{M}_{z},\hfill & \hfill {v}_{B}=\left(\frac{hd}{{E}_{1x}{I}_{1z}}+\frac{{d}^{2}}{2{E}_{2x}{I}_{2z}}\right){M}_{z},\hfill & \hfill {\theta}_{z}=\left(\frac{h}{{E}_{1x}{I}_{1z}}+\frac{d}{{E}_{2x}{I}_{2z}}\right){M}_{z}\hfill \end{array}\text{.}$$(4)These inplane results can be summarized in the compliance matrix where coupling between the 3° of freedom appears clearly. $$\left(\begin{array}{c}\hfill {u}_{B}\hfill \\ \hfill {v}_{B}\hfill \\ \hfill {\theta}_{z}\hfill \end{array}\right)=\left[\begin{array}{ccc}\hfill \frac{{h}_{\mathrm{}}^{3}}{3{E}_{1x}{I}_{1z}}+\frac{h}{{G}_{1xy}{S}_{1}}\hfill & \hfill \frac{{h}_{\mathrm{}}^{2}d}{2{E}_{1x}{I}_{1z}}\hfill & \hfill \frac{{h}^{2}}{2{E}_{1x}{I}_{1z}}\hfill \\ \hfill \frac{{h}_{\mathrm{}}^{2}d}{2{E}_{1x}{I}_{1z}}\hfill & \hfill \frac{h{d}^{2}}{{E}_{1x}{I}_{1z}}+\frac{{d}_{\mathrm{}}^{3}}{3{E}_{2x}{I}_{2z}}+\frac{d}{{G}_{2xy}{S}_{2}}\hfill & \hfill \frac{hd}{{E}_{1x}{I}_{1z}}+\frac{{d}^{2}}{2{E}_{2x}{I}_{2z}}\hfill \\ \hfill \frac{{h}^{2}}{2{E}_{1x}{I}_{1z}}\hfill & \hfill \frac{hd}{{E}_{1x}{I}_{1z}}+\frac{{d}^{2}}{2{E}_{2x}{I}_{2z}}\hfill & \hfill \frac{h}{{E}_{1x}{I}_{1z}}+\frac{d}{{E}_{2x}{I}_{2z}}\hfill \end{array}\right]\left(\begin{array}{c}\hfill {F}_{x}\hfill \\ \hfill {F}_{y}\hfill \\ \hfill {M}_{z}\hfill \end{array}\right)\text{.}$$(5)Outofplane loadings are force F_{z} and torques M_{x} and M_{y}. Some complementary torsion deformations superpose the shear and bending effects already observed in the in plane problem. The outofplane compliance matrix also shows coupling between rotations dof. $$\left(\begin{array}{c}\hfill {w}_{B}\hfill \\ \hfill {\theta}_{x}\hfill \\ \hfill {\theta}_{y}\hfill \end{array}\right)=\left[\begin{array}{ccc}\hfill \frac{h{d}_{\mathrm{}}^{2}}{{G}_{1xz}{J}_{1}}+\frac{{h}_{\mathrm{}}^{3}}{3{E}_{1x}{I}_{1y}}+\frac{{d}_{\mathrm{}}^{3}}{3{E}_{2x}{I}_{2y}}+\frac{h}{{G}_{1xz}{S}_{1}}+\frac{d}{{G}_{2xz}{S}_{2}}\hfill & \hfill \frac{{h}_{\mathrm{}}^{2}}{2{E}_{1x}{I}_{1y}}\hfill & \hfill \left(\frac{hd}{{G}_{1xz}{J}_{1}}+\frac{{d}^{2}}{2{E}_{2x}{I}_{2y}}\right)\hfill \\ \hfill \frac{{h}_{\mathrm{}}^{2}}{2{E}_{1x}{I}_{1y}}\hfill & \hfill \frac{h}{{E}_{1x}{I}_{1y}}+\frac{d}{{G}_{2xz}{J}_{2}}\hfill & \hfill 0\hfill \\ \hfill \left(\frac{hd}{{G}_{1xz}{J}_{1}}+\frac{{d}^{2}}{2{E}_{2x}{I}_{2y}}\right)\hfill & \hfill 0\hfill & \hfill \frac{h}{{G}_{1xz}{J}_{1}}+\frac{d}{{E}_{2x}{I}_{2y}}\hfill \end{array}\right]\left(\begin{array}{c}\hfill {F}_{z}\hfill \\ \hfill {M}_{x}\hfill \\ \hfill {M}_{y}\hfill \end{array}\right)\text{.}$$(6)3D finite element calculations are managed in the SolidWorks environment and the results are compared to these analytical expressions in order to identify the influence of the connection elements (screw, nut, pins) and the influence of local contact effects.
2.2 3D FEA of the real corner behavior and focus on rotational stiffness
To carry out the finite element simulation of the “real” corner including contacts, pins, screw and nuts, one must address several topics: (i) the material characteristics (ii) the meshing refinement (iii) the contact definition and (iv) the bolt preload. This last was the trickiest and we used the temperature difference between the bolt and the wooden parts to model the preload: (i) first assembly is created without any gap or preload, (ii) the temperature of the screw is imposed to a small value and (iii) the thermal problem is computed and the contraction of the screw generates an axial preload in the connection.
Finally, the determination of the displacements u_{B} and v_{B} of point B and the rotation θ_{z} of the section at B cannot be obtained directly from the finite element results. It is necessary to make an extraction of the nodes coordinates X_{i}, Y_{i} and Z_{i} and displacements u_{i} and v_{i} to obtain this angle using the least square method. The displacement field of a assumed rigid section can be represented by: $$\left\{{U}_{\text{section}B}\right\}={\left\{\begin{array}{c}\hfill {\displaystyle \overrightarrow{\mathrm{\Psi}}}={\theta}_{x}{\displaystyle \overrightarrow{x}}+{\theta}_{y}{\displaystyle \overrightarrow{y}}+{\theta}_{z}{\displaystyle \overrightarrow{z}}\hfill \\ \hfill {\displaystyle \overrightarrow{{u}_{{M}_{i}}}}=\left({u}_{B}+{\theta}_{y}{Z}_{i}{\theta}_{z}{Y}_{i}\right){\displaystyle \overrightarrow{x}}+\left({v}_{B}{\theta}_{x}{Z}_{i}\right){\displaystyle \overrightarrow{y}}+\left({w}_{B}+{\theta}_{x}{Y}_{i}\right){\displaystyle \overrightarrow{z}}\hfill \end{array}\right\}}_{{M}_{i}}\text{.}$$(7) The sum of squared difference between the rigid displacements and the finite element nodal displacements is given by: $$S={\displaystyle \sum _{i=1}^{N}}{e}_{i}^{2}={\displaystyle \sum _{i=1}^{N}}\left\{{\left({u}_{B}+{\theta}_{y}{Z}_{i}{\theta}_{z}{Y}_{i}{u}_{i}\right)}^{2}+{\left({v}_{B}{\theta}_{x}{Z}_{i}{v}_{i}\right)}^{2}+{\left({w}_{B}+{\theta}_{x}{Y}_{i}{w}_{i}\right)}^{2}\right\}\text{,}$$(8)where N is the total number of nodes at the face B. Minimizing S with respect to u_{B}, v_{B}, w_{B}, θ_{x}, θ_{y}, θ_{z} yields to a system of 6 linear equations that gives the needed rotations and displacements for the comparison. The linear system can be written in a matrix form as follows: $$\left[\begin{array}{cccccc}\hfill N\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\displaystyle \sum _{i=1}^{N}}{Z}_{i}\hfill & \hfill {\displaystyle \sum _{i=1}^{N}}{Y}_{i}\hfill \\ \hfill 0\hfill & \hfill N\hfill & \hfill 0\hfill & \hfill {\displaystyle \sum _{i=1}^{N}}{Z}_{i}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill N\hfill & \hfill {\displaystyle \sum _{i=1}^{N}}{Y}_{i}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\displaystyle \sum _{i=1}^{N}}{Z}_{i}\hfill & \hfill {\displaystyle \sum _{i=1}^{N}}{Y}_{i}\hfill & \hfill {\displaystyle \sum _{i=1}^{N}}({{Y}_{i}}^{2}+{{Z}_{i}}^{2})\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\displaystyle \sum _{i=1}^{N}}{Z}_{i}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\displaystyle \sum _{i=1}^{N}}{{Z}_{i}}^{2}\hfill & \hfill {\displaystyle \sum _{i=1}^{N}}{Y}_{i}{Z}_{i}\hfill \\ \hfill {\displaystyle \sum _{i=1}^{\mathrm{}}}{Y}_{i}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\displaystyle \sum _{i=1}^{N}}{Y}_{i}{Z}_{i}\hfill & \hfill {\displaystyle \sum _{i=1}^{N}}{{Y}_{i}}^{2}\hfill \end{array}\right]\left(\begin{array}{c}\hfill {u}_{B}\hfill \\ \hfill {v}_{B}\hfill \\ \hfill {w}_{B}\hfill \\ \hfill {\theta}_{x}\hfill \\ \hfill {\theta}_{y}\hfill \\ \hfill {\theta}_{z}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill {\displaystyle \sum _{i=1}^{N}}{u}_{i}\hfill \\ \hfill {\displaystyle \sum _{i=1}^{N}}{v}_{i}\hfill \\ \hfill {\displaystyle \sum _{i=1}^{N}}{w}_{i}\hfill \\ \hfill {\displaystyle \sum _{i=1}^{N}}({Y}_{i}{w}_{i}{Z}_{i}{v}_{i})\hfill \\ \hfill {\displaystyle \sum _{i=1}^{N}}{Z}_{i}{u}_{i}\hfill \\ \hfill {\displaystyle \sum _{i=1}^{N}}{Y}_{i}{u}_{i}\hfill \end{array}\right)\text{.}$$(9)
Thanks to the numerous bibliographic references [8–18] the mean values for the spruce properties are well known. The longitudinal Young modulus E_{ℓ} of the wood is greater than the moduli E_{r} and E_{t} in the cross section of the tree. E_{r} is the modulus in the radial direction of the expansion of the tree and is more or less 4% of the longitudinal modulus while E_{t} that is the modulus measured tangentially to the wood rings of growth is about 8%. Nevertheless, these two are significantly lower than E_{ℓ} and we will assume that the elastic behavior is transversally isotropic. Values of Table 1 are implemented in the finite element code for the simulations.
The finite element used is a quadratic tetrahedron with a characteristic length of 2 mm in regular parts and the mesh is refined near the connections where the size reduced to 0.5 mm. Junction needs also more than one hundred contact elements (see Fig. 2). This leads to more than 350 000 dof and, because of the nonlinearity of the contact element, one simulation needs 5–15 min to be completed on a PC using Intel i7 processor. This is already quite long for the corner only and one can easily understand that the entire bunk bed would lead to unreasonable CPU time. When furniture's elements are assembled, junctions are tightened and a preload appears between the assembled elements. To apply this preload on the finite element model there are actually no simple option in SolidWorks so we use the thermal properties and gave an initial temperature lower to the bolt than for all other parts that where at ambient temperature. In case of a classical assembly, the bolt force can be written in function of ΔT as: $${F}_{0}=\frac{{l}_{0}\alpha}{(\frac{1}{{k}_{b}}+\frac{{k}_{c1}+{k}_{c2}}{{k}_{c1}{k}_{c2}})}\mathrm{\Delta}T=\beta \mathrm{\Delta}T\text{.}$$(10)
k_{c}_{1} and k_{c}_{2} are stiffnesses of the assembled part 1 and 2 and k_{b} the bolt stiffness. In this relation α denotes the thermal expansion coefficient for the bolt material (α = 1.5 × 10^{−5}K^{−1}) and ΔT is the temperature change imposed on the bolt. Because k_{c}_{1} and k_{c}_{2} are global values not easily known, calculating the coefficient β that appears in equation (10) necessitates finite element calculation of the stiffness of two differently sized, compressed and oriented wood parts. Simulation gives β equal to a constant value of 8.74 N/K and equation (10) can be used to determine the temperature variation needed to generate the desired preload. Some experiences with different users were carried out and the maximum torque was measured for each tightening: the mean value was chosen to evaluate a reasonable preload value. Because of the low strength of wood in transversal direction the dispersion on all experiences is not so important.
For the inplane case, considering a nodal spring between the two wooden parts introduces 3 complementary terms in the matrix of equation (5). Our approach gives 5 potential relations from which one can identify the 3 values of k_{u}, k_{v} and k_{γ}: $$\left(\begin{array}{c}\hfill {u}_{B}\hfill \\ \hfill {v}_{B}\hfill \\ \hfill {\theta}_{z}\hfill \end{array}\right)=\left[\begin{array}{ccc}\hfill \frac{{h}_{\mathrm{}}^{3}}{3{E}_{1x}{I}_{1z}}+\frac{h}{{G}_{1xy}{S}_{1}}\begin{array}{c}\hline +\frac{1}{{k}_{u}}\\ \hline\end{array}\hfill & \hfill \frac{{h}_{\mathrm{}}^{2}d}{2{E}_{1x}{I}_{1z}}\hfill & \hfill \frac{{h}^{2}}{2{E}_{1x}{I}_{1z}}\hfill \\ \hfill \frac{{h}_{\mathrm{}}^{2}d}{2{E}_{1x}{I}_{1z}}\hfill & \hfill \frac{h{d}^{2}}{{E}_{1x}{I}_{1z}}+\frac{{d}_{\mathrm{}}^{3}}{3{E}_{2x}{I}_{2z}}+\frac{d}{{G}_{2xy}{S}_{2}}\begin{array}{c}\hline +\frac{1}{{k}_{v}}+\frac{{d}^{2}}{{k}_{{\theta}_{z}}}\\ \hline\end{array}\hfill & \hfill \frac{hd}{{E}_{1x}{I}_{1z}}+\frac{{d}^{2}}{2{E}_{2x}{I}_{2z}}\begin{array}{c}\hline +\frac{d}{{k}_{{\theta}_{z}}}\\ \hline\end{array}\hfill \\ \hfill \frac{{h}^{2}}{2{E}_{1x}{I}_{1z}}\hfill & \hfill \frac{hd}{{E}_{1x}{I}_{1z}}+\frac{{d}^{2}}{2{E}_{2x}{I}_{2z}}\begin{array}{c}\hline +\frac{d}{{k}_{{\theta}_{z}}}\\ \hline\end{array}\hfill & \hfill \frac{h}{{E}_{1x}{I}_{1z}}+\frac{d}{{E}_{2x}{I}_{2z}}\begin{array}{c}\hline +\frac{1}{{k}_{{\theta}_{z}}}\\ \hline\end{array}\hfill \end{array}\right]\left(\begin{array}{c}\hfill {F}_{x}\hfill \\ \hfill {F}_{y}\hfill \\ \hfill {M}_{z}\hfill \end{array}\right)\text{.}$$(11) If F_{X} ≠ 0 and F_{y} = M_{z} = 0 the first column gives a relation to identify k_{u} from the calculated displacement u_{B} obtained from the 3D finite element simulation. $${k}_{u}=\frac{1}{\frac{{u}_{B}}{{F}_{x}}\left(\frac{{h}_{\mathrm{}}^{3}}{3{E}_{1x}{I}_{1z}}+\frac{h}{{G}_{1xy}{S}_{1}}\right)}\text{.}$$(12)In the case of the k_{θz} stiffness, one can use three different relations: (i) F_{Y} ≠ 0 and F_{X} = M_{z} = 0 and using the angle θ_{z} coming from the finite element simulation; (ii) M_{z} ≠ 0 and F_{X} = F_{y} = 0 and using the displacement v_{B} coming from the finite element simulation or (iii) M_{z} ≠ 0 and F_{X} = F_{y} = 0 and using the angle θ_{z} coming from the finite element simulation. All three cases give quasiidentical values (difference less than 5%) and for sake of simplicity we will use the case (i) in the following, which leads to: $${k}_{{\theta}_{z}}=\frac{d}{\frac{{\theta}_{z}}{{F}_{y}}\left(\frac{hd}{{E}_{1x}{I}_{1z}}+\frac{{d}^{2}}{2{E}_{2x}{I}_{2z}}\right)}\text{.}$$(13)Finally, one can evaluate k_{v} from the case F_{y} ≠ 0 and F_{x} = M_{z} = 0 and using the displacement v_{B} coming from the finite element simulation. Stiffness k_{θz} appears in the relation so it must be evaluated previously: $${k}_{v}=\frac{1}{\frac{{v}_{B}}{{F}_{y}}\left(\frac{h{d}^{2}}{{E}_{1x}{I}_{1z}}+\frac{{d}_{\mathrm{}}^{3}}{3{E}_{2x}{I}_{2z}}+\frac{d}{{G}_{2xy}{S}_{2}}+\frac{{d}^{2}}{{k}_{\gamma}}\right)}\text{.}$$(14)Numerical applications clearly show that elongational stiffnesses k_{u} and k_{v} have high values and that their influence in the displacement components can be neglected in regard of the influence of k_{θZ}. Consequently, the 0D spring element is reduced to one inplane angular rigidity k_{θZ} and two outofplane angular rigidities k_{θX} and k_{θY}.
Transversal isotropic elastic characteristics for spruce.
Fig. 2 CAD and meshing of the “real” corner. 
2.3 A nonlinear model for rotational stiffness
Varying the force F_{Y}, the previous procedure gives different values of the spring stiffness k_{θz}: high value for low force and a decreasing evolution of the stiffness to an asymptotic value when the load increases. A constant value of k_{θz} would lead to a too smooth structure. Consequently, the junction behavior has to be nonlinear and one must compute the behavior law between the bending moment M_{z} across the junction and the relative rotation θ_{rz} between the two wooden parts.
Moment M_{z} is obtained from the load values F_{Y} by multiplying with the beam length d. Thanks to the relation M_{z} = k_{θz}θ_{rz} one can obtain the relative rotation θ_{rz} for each value of F_{Y}: $${\theta}_{rz}=\frac{{M}_{z}}{{k}_{{\theta}_{z}}}=\frac{d{F}_{Y}}{{k}_{{\theta}_{z}}}\text{.}$$(15)
This leads to the nearly linear plot of Figure 3 that we chose to model with a decreasing exponential function asymptotic to a linear one as in equation 16: $${M}_{z}=\left({M}_{0}+{k}_{\infty}{\theta}_{rz}\right)\left(1\mathrm{exp}\left(\frac{\left{\theta}_{rz}\right}{{\theta}_{\text{ref}}}\right)\right)\text{.}$$(16) k_{∞} is the slope of the linear part, M_{0} is the moment value when the linear part crosses the vertical line at the origin of the graph and θ_{ref} is an angular value that characterizes the length of nonlinear zone. In Figure 3, the line is the best fit of these parameters to represent the spots coming from the multiscale approach. The best fit leads to values summarized in Table 2.
The same identification procedure is carried out for the two other rotational stiffnesses k_{θx} and k_{θy} and, after transformation into a moment versus relative rotation graphs, it leads to the evolution shown on Figure 4. The evolution of momentum are quasi linear so the identification will be simple.
Apparently, the graphs of Figure 4 indicate a quasilinear behavior of the connection for both k_{θy} and k_{θx} but with stiffness 10 times higher in torsion (k_{θx}) then for bending around the Y direction (k_{θy}). However, the relative difference in flexibility between the rigid and semirigid system was found to be of 5.8% for torque around the X direction implying that the connection is stiff in torsion, whereas k_{θy} bestfit value is 1.615 × 10^{6} N.mm.rad^{−1}.
Fig. 3 Left, stiffness vs. moment and right, moment vs. relative rotation (results and model). 
M_{z} vs. θ_{rz} model parameters.
Fig. 4 Spring rotational behavior, k_{θx} and k_{θy} stiffnesses. 
3 Numerical and experimental validations
The nonlinear connection model has been implemented in the Timoshenko's beam code developed for furniture applications. Such code should give an identical deformed shape than a 3D simulation for a given structure. First, to validate the efficiency of the connection model to reproduce the junction behavior we carried out a complete 3D simulation on a plane structure representative of the bunk bed. Finally, it is also necessary to manage real tests on corner to validate experimentally the behavior of the jointing elements and its domain of validity.
3.1 Numerical validation of the approach on simplified 2D examples
At beam element scale, the connection is represented by fictitious zero length (0D) elastic element inserted at the intersection point between a beam and a column as used in [19,20]. For the inplane configuration, this element can be modeled by a multisprings system dividing the connection into equivalent shear, axial and rotational springs. The stiffness of each spring (denoted by k_{u}, k_{v} and k_{θ}, respectively) is the ratio of the transmitted force to the corresponding relative displacement within the connection. By satisfying the equilibrium requirements, the stiffness of the connection can be represented in the finite element beam code by the following stiffness matrix: $$\left[Kc\right]=\lfloor \begin{array}{cccccc}\hfill {k}_{u}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {k}_{u}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {k}_{v}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {k}_{v}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {k}_{\theta}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {k}_{\theta}\hfill \\ \hfill {k}_{u}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {k}_{u}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {k}_{v}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {k}_{v}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {k}_{\theta}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {k}_{\theta}\hfill \end{array}\rfloor \text{.}$$(17) The global structure stiffness matrix is obtained by assembling the frame elements matrices, including the spring's behavior, according to the conventional direct stiffness procedure programmed in a Matlab code which is also adapted to run a Newton–Raphson routine for nonlinear analysis.
A simplified frame of one side of a bunk bed was modeled at both scales. The semi rigid frame model includes a beam of sectional area (20 × 140 mm) and a varying length d and two columns of sectional area (57 × 57 mm) and a varying length h. The modulus of elasticity and the shear modulus are the same for all frame elements. The length of the beams is 1 and 2 m whereas the length of the different columns is 0.5, 2 or 3 m. These standard dimensions represent all subframes that may be found in a typical bunk bed. Two flexible connections are located at the intersection of beam and column and modeled as nonlinear rotational springs which properties are derived using equation 17.
The 3D model with the real configuration of the assembly technique was modeled in Solid works environment using the same solid element, materials properties, pretention value, mesh density and contact set as in the corner model. The frame was constrained at the base and loaded inplane at mid span with a vertical force of 1000 N. The nodal displacements of the 3D model (blue dots in Fig. 5) were extracted and implemented in the beam code to be compared with the deflected neutral axis of frame elements. The relative error in displacement ranged from 2% to 12% which falls in the expected range of error previously observed between the two scales.
In the corner case used to manage identification the maximum bending momentum is localized on the corner. In the case of the example, Figure 5, which is hyperstatic, the momentum value is maximum in the middle of the horizontal beam: the momentum distribution is really different from the reference case so the good coincidence between deformed shape at both scales is a first validation of the quality of the “beam + connection” modeling. It can be noticed that for the corner case, considering that stiffness varies with the applied momentum, the contribution of the spring in the total displacement goes from 13% for small loads to 56% for high loads. In the case shown in Figure 5, the contribution of the springs at max displacement is 23%.
Fig. 5 Comparison between the displacement field of the mean line of the 3D FEA and the “beam + spring” one. 
3.2 Experimental validation: nonlinear behavior and dispersion
Tests were performed on the corner case using a universal testing machine DeltaLab500 provided with 5 kN capacity loading cell. A special fixation system was designed and attached to the machine. Each specimen was first clamped to the fixation system via the post (column part). A vertical increasing static load was uniformly applied along the complete width of the rail (beam part) through a (20 × 10 × 2) mm metallic plate slicked to the rail free edge. Digital image correlation technique [21] was also employed where two digital cameras were installed 70 and 80 cm in front of the specimen face. The test arrangement is illustrated in Figure 6.
The tests were performed until a major failure occurred at a constant rate of 2.5 mm/min and was programmed to pause every 30 N of load increment in order to take photos. The primary requirements of the test setup were the measurements of:

the applied load F_{y} which is captured directly by the load cell;

the vertical displacement v_{B} at the free edge which is measured as the movement of the crosshead;

a direct measurement of the relative rotation between the post and the rail for each load increment using DIC.
The load displacement data was recorded every one second while the rotation was posttreated at each load increment which enables to reproduce the Momentum versus relative rotation curve as presented in the previous section (see Fig. 7).
One can see that most of experimental curves coincide with the 3D finite element simulation for low momentum values (i.e., M < 50 × 10^{3} N.mm). It is a first interesting result because simulation managed on complete bunk bed in normalized loading conditions showed that bending moment does not exceed 50 × 10^{3} N.mm in connections. Consequently, the finite element beam code with nonlinear connections is a perfect tool to manage validation before manufacturing real furniture.
Nevertheless, if one wants to identify the connection model from equation 16 it appears that (i) k_{∞} is much lower than from the finite element simulation; (ii) M_{o} highlights a large dispersion up to 22% and (iii) the initial slope k_{o} that can be identified to the ratio M_{o} over θ_{ref} highlights even higher dispersion: up to 30%! Where does all these dispersion come from is the question addressed in the next section.
Fig. 6 Presentation of the corner test experiment (a) details on the apparatus (b) angle gap θ_{rz} measurement with DIC. 
Fig. 7 M_{z} versus θ_{rz} curves compared with the FEA curve − dots are direct measurement obtained from the displacement of the testing machine and the line is obtained from 3D finite element simulation. 
3.3 Complementary test and identification of the dispersion's origin
When momentum increases the slope of the momentum versus rotation decreases and reaches a stable value lower than what the 3D simulation can estimate. More, the momentum value of the divergence between simulation and experimental tests is much dispersed: one can wonder why and answers are to be found in the low value of the limit of proportionality in radial or tangential direction of the wood. Since contact is localized between the head of the screw and the column the pressure can exceed the elastic limit and lead to nonlinear behavior. When the relative rotation becomes important there is no more contact on most of the area between the column and the beam: at this time the pressure is also high on the remaining surface between these two wooden parts and one must notice that the vertical part is compressed in the direction normal to fiber. Both effects conjugate to obtain this nonlinear global behavior.
25 compression tests have been managed to reproduce the contact between the head of the screw and the wooden part. The compression is managed on faces that are normal to the fiber direction. A hole of 7 mm was drilled across the center of each cube and an M6 bolt was cut to a length of 30 mm and was inserted in the hole. The bolt head was compressed through a metallic pad of the same diameter attached to the loading cell of the Deltalab machine. The load was applied with constant rate of 0.25 mm/min.
Figure 8 gives an illustration of the test and plots the different curves obtained. One can observe that initial slope K_{init} highlights dispersion as well as the F_{p} values (F_{p} is the load value obtained from the intersection of the linear part of the curve with the initial slope). The slope of the linear part is denoted K_{∞} and is more or less the same for all curves. The dispersion on F_{p} is the same that the dispersion on the M_{o} parameter of the behavior law of the connection, about 20% dispersion, so this confirms that the origin of dispersion has been identified. More, if we introduce the supplementary relative displacement in the behavior of the corner, one can reproduce the dispersion observed experimentally.
The results of the tests represented as forcedisplacement curves in Figure 8 exhibited an elasticplastic behavior of wood in compression perpendicular to grain. In the same manner than for the corner tests, we can fit a three parameter exponential model defined in equation (18) with numerical values defined in (Tab. 3) for each curve: $$F=\left({F}_{init}+{K}_{\infty}u\right)\left(1\mathrm{exp}\left(\frac{u}{{u}_{ref}}\right)\right)\text{,}$$(18) where F_{init} is the yintercept, K_{∞} is the slope in the second linear phase and u_{ref} is the scaling factor of displacement. Searching the bestfit parameters for each curve yielded the following average value of the curves characteristics and dispersions:
Assigning the mean, minimum and maximum values in equation (18) permits to plot the black dot curves of Figure 8 which represent the mean curve with two curves that envelope the test set quite well.
It is well known that mechanical properties of wood show dispersions related to the density. It is also well known that radial modulus E_{r} or tangential modulus E_{t} are significantly lower than the modulus E_{ℓ} in longitudinal direction but also very different from each other. In this case one can evaluate the modulus in a given direction in a cross section of the wooden post from the classical equation (19): $$E(\alpha )=\frac{{\sigma}_{\alpha \alpha}}{{\u03f5}_{\alpha \alpha}}=\frac{1}{\frac{(\mathrm{cos}{\alpha}^{4}{v}_{tr}\mathrm{sin}{\alpha}^{2}\mathrm{cos}{\alpha}^{2})}{{E}_{t}}+\frac{(\mathrm{sin}{\alpha}^{4}{v}_{rt}\mathrm{sin}{\alpha}^{2}\mathrm{cos}{\alpha}^{2})}{{E}_{r}}+\frac{\mathrm{sin}{\alpha}^{2}\mathrm{cos}{\alpha}^{2}}{{G}_{tr}}}\text{.}$$(19) Figure 9 shows the 25 cross sections of the post where the 25 different α_{i} angles have been measured. On pictures, the face shown is the normal to fiber face of the cubic specimen. The horizontal top face is the one being compressed and α_{i} is the angle between the vertical line (normal to compression) and the ring orientation.
Mean angle value on all α_{i} is quasi naught but because of the symmetry. For each angle measured, the positive value of the apparent E_{t} Young's modulus is then possible to calculate from equation (19) considering Poisson's ratios, G_{tr}, E_{t} and E_{r} values obtained from bibliographical references. Because of the low value of the shear modulus G_{tr} in the cross section of the tree, the mean value of E_{t} appears to be 220 MPa which is significantly lower than E_{t} or E_{r} of bibliographical references. The standard deviation of these measures is quite high and one obtains a dispersion of 23% on all specimens used to manage the compression tests. This geometrical uncertainty due to the orientation of the wood product coming from wood factory gives important dispersion that (i) can explain the dispersion observed on compression test but also on the corners tests and (ii) clearly indicates that numerical simulations of tests on wood products must take into account this dispersion to provide accurate information on the furniture's behavior.
This dispersion is truly related to the transversally isotropic behavior that is not a reasonable assumption. But, considering that it is also not reasonable to measure exactly the transversal orientation of all beams that are involved in wood furniture has a bunk bed for example, one cannot avoid this assumption in simulations. Taking into account dispersion due to this lack of knowledge on the orientation seems to be a good compromise.
Fig. 8 Screw compression in wood test and results. 
Fig. 9 25 pictures of cube section − influence of the radial and ring orientation on the elastic behavior. 
4 Conclusion: toward virtual test for wood product validation and nonlinear wood behavior identification
Using a multiscale approach, a numerical procedure to identify the connection contribution in the furniture's stiffness has been presented. Comparing 3D finite element calculations with a Timoshenko's beam calculation on a corner of two wooden parts, we can identify the specific behavior of the connection elements (pins, nut, screw… and local 3D effects) to introduce it as a punctual 0D element in the beam code. This procedure has been carried out for a specific connection technology, 1 bolt and 2 pins, and it has been showed that the low rotational stiffness must be taken into account in a global calculation. In comparison the elongational stiffness effects on the displaced shape can be neglected in regard of the rotational ones. At least one of the three stiffnesses depends on the bending moment value and must be modeled by a nonlinear law.
The approach has been validated first by a numerical simulation on a complete structure. The accuracy of the approach is proved by comparing the result of the beam code where the 0D nonlinear element has been implemented with a complete 3D finite element simulation on a plane structure of wooden furniture.
An experimental validation by a global bending test has also been carried out. The displacement field has been measured by using digital image correlation. Both experimental and numerical results agree for low value of the bending moment passing through the connection (i.e. M < 50 N.m). Considering the low values of momentum in corners during the normalized tests on bunk beds (about 30 N.m) the purely numerical multiscale gives an accurate behavior of the real assembly.
When bending moment is higher than 50 N.m, divergence appears between the simulation and the experiments with a severe slope decreasing and high dispersion from one test to another. This behavior is due to the head of the screw that generates high pressure on the wood in a transversal direction that is no more elastic. Dispersion appears to be due to the fluctuation of the orientation of the wooden part and the important difference between E_{t} and E_{r} in spruce properties.
To be able to completely reproduce the experimental behavior for bending moments higher than 50 N.m one needs to go further in the characterization of the nonlinear behavior of the wood in the direction normal to fiber. Three point bending tests are carried out actually using DIC to measure the displacement field and finite element model upgrade method is managed to identify these properties in the linear and nonlinear domain. This work is actually going on at MSME and first results have been presented in Wood GDR recently in Makhlouf et al. [22].
Acknowledgement
Authors would like to thank FCBA for financial support.
References
 A. Kasal, Determination of the strength of various sofa frames with finite element analysis, G.U.J. Sci. 19 (2006) 191–203 [Google Scholar]
 H.I. Demirci, The experimental and finite element analysis of diagonal tensile tests conducted on frametype constructed corner joints, Technol 14 (2011) 11–21 [Google Scholar]
 EN 747, Furniture − Bunk beds and high beds for domestic use − Part 2: test methods, 2012 [Google Scholar]
 T. Nicholls, R. Crisan, Study of the stress strain state in corner joints and boxtype furniture using Finite Element Analysis (FEA), Holz als Roh und Werkstoff. 60 (2002) 66–71 [Google Scholar]
 H. Makhlouf, L. Chevalier, E. Launay, E. Favier, A stochastic approach for the evaluation of the reliability of wood furniture in an industrial context: managing virtual standardisation tests, Mech. Ind. 2015a [Google Scholar]
 Y.Z. Erdil, J. Zhang, C.A. Eckelman, Withdrawal and bending strength of dowelnuts in plywood and oriented strandboard, FPJ 53 (2003) 54–57 [Google Scholar]
 C.A. Eckelman, Holding strength of through bolt with dowel nut connections, FPJ 39 (1989) 41–48 [Google Scholar]
 H. Carrington, The elastic constants of spruce, Philos. Mag. Lond. 45 (1923) 1055–1057 [CrossRef] [Google Scholar]
 J. Campredon, Contribution à l'étude des propriétés élastiques des bois, Ann. Éc. Natl. Eaux For. Stn. Rech. Exp. 3 (1935) 251–288 [Google Scholar]
 F.H. Neuhaus, Elastizitatszahlen von Fichtenholz, Diss. Ruhruniversitat, Bochum, 1981 [Google Scholar]
 D. Guitard, Mécanique du matériau bois et composites, Cepadues, Toulouse, France, 1987 [Google Scholar]
 J. Bodig, B. Jayne, Mechanics of wood and wood composites, Krieger, Malabar, Florida, USA, 1993 [Google Scholar]
 J. Charles, Wood properties, Encycl. Agric. Sci. 4 (1994) 549–561 [Google Scholar]
 E. Kahle, J. Woodhouse, The influence of cell geometry on elasticity of softwood, J. Mater. Sci. 29 (1994) 1250–1259 [Google Scholar]
 J.E. Winandy, Wood properties, USDAForest Service, Forest Products Laboratory, 1994, pp. 549–561 [Google Scholar]
 USDA, Wood handbook d wood as an engineering material, Forest Products Society, 1999 [Google Scholar]
 D. Keunecke et al., Determination of Young's and shear moduli of common yew and Norway spruce by means of ultrasonic waves, Wood Sci. Technol. 41 (2007) 309–327 [Google Scholar]
 P. Niemz, D. Caduff, Research into determination of the Poisson ratio of spruce wood, Holz. Roh. Werkst. 66 (2008) 1–4 [CrossRef] [Google Scholar]
 S.L. Chan, P.P.T. Chui, Non linear static and cyclic analysis of steel frames with semirigid connections, Elsevier, UK, 2000 [Google Scholar]
 L.M.C. Simões, Optimization of frames with semi rigid connections, Comput. Struct. 60 (1996) 531–539 [Google Scholar]
 L. Chevalier, S. Calloch, F. Hild, Y. Marco, Digital correlation used to analyse the multiaxial behavior of rubberlike materials, Eur. J. Mech. A/Solids 20 (2001) 169–187 [CrossRef] [Google Scholar]
 H. Makhlouf, L. Chevalier, T.T. Nguyen, E. Launay, Identification du comportement anisotrope du bois par comparaison calcul éléments finis et analyse d'images, journées scientifiques du GDR Bois, 46 novembre2015, Clermont Ferrand, 2015b [Google Scholar]
Cite this article as: L. Chevalier, H. Makhlouf, B. JacquetFaucillon, E. Launay, Modeling the influence of connecting elements in wood products behavior: a numerical multiscale approach, Mechanics & Industry 19, 301 (2018)
All Tables
All Figures
Fig. 1 Screen shot of the CAD model of the corner. Lengths are h = 447 m and d = 428.5 mm. Sectional areas are 57 × 57 mm^{2} edge for the vertical part and 20 × 140 mm^{2} for the horizontal one. All dof are naught in section A and loads are applied in section B. 

In the text 
Fig. 2 CAD and meshing of the “real” corner. 

In the text 
Fig. 3 Left, stiffness vs. moment and right, moment vs. relative rotation (results and model). 

In the text 
Fig. 4 Spring rotational behavior, k_{θx} and k_{θy} stiffnesses. 

In the text 
Fig. 5 Comparison between the displacement field of the mean line of the 3D FEA and the “beam + spring” one. 

In the text 
Fig. 6 Presentation of the corner test experiment (a) details on the apparatus (b) angle gap θ_{rz} measurement with DIC. 

In the text 
Fig. 7 M_{z} versus θ_{rz} curves compared with the FEA curve − dots are direct measurement obtained from the displacement of the testing machine and the line is obtained from 3D finite element simulation. 

In the text 
Fig. 8 Screw compression in wood test and results. 

In the text 
Fig. 9 25 pictures of cube section − influence of the radial and ring orientation on the elastic behavior. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.