© 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 beam-like elements of wood or wood based materials, such as chipboard, Medium Density Fiberboard (MDF for short) or plywood;

• can type structure constructed of plate-like 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 wood-based 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 multi-scale 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 multi-scale 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 non-linear 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 multi-scale analysis of the jointing element influence

A through-bolt with dowel-nut connection is studied in this section. This system is usually used both as primary connectors and also to reinforce weaker joints [6]. The through-bolt with a dowel-nut 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 multi-scale 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 multi-scale 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 mm2 edge for the vertical part and 20 × 140 mm2 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 in-plane 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}{3E_{1x}{I}_{1z}}+\frac{h}{G_{1xy}{S}_1}\right)F_x,\hfill & \hfill v_B=-\frac{h^{2}d}{2E_{1x}{I}_{1z}}{F}_x,\hfill & \hfill {\theta }_z=-\frac{h^2}{2E_{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}{2E_{1x}{I}_{1z}}{F}_y,\hfill & \hfill v_B=\left(\frac{h{d}^2}{E_{1x}{I}_{1z}}+\frac{d_{\mathrm{}}^3}{3E_{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}{2E_{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}{2E_{1x}{I}_{1z}}{M}_z,\hfill & \hfill v_B=\left(\frac{hd}{E_{1x}{I}_{1z}}+\frac{d^2}{2E_{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 in-plane 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}{3E_{1x}{I}_{1z}}+\frac{h}{G_{1xy}{S}_1}\hfill & \hfill -\frac{h_{\mathrm{}}^{2}d}{2E_{1x}{I}_{1z}}\hfill & \hfill -\frac{h^2}{2E_{1x}{I}_{1z}}\hfill \\ \hfill -\frac{h_{\mathrm{}}^{2}d}{2E_{1x}{I}_{1z}}\hfill & \hfill \frac{h{d}^2}{E_{1x}{I}_{1z}}+\frac{d_{\mathrm{}}^3}{3E_{2x}{I}_{2z}}+\frac{d}{G_{2xy}{S}_2}\hfill & \hfill \frac{hd}{E_{1x}{I}_{1z}}+\frac{d^2}{2E_{2x}{I}_{2z}}\hfill \\ \hfill -\frac{h^2}{2E_{1x}{I}_{1z}}\hfill & \hfill \frac{hd}{E_{1x}{I}_{1z}}+\frac{d^2}{2E_{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)Out-of-plane 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 out-of-plane 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}{3E_{1x}{I}_{1y}}+\frac{d_{\mathrm{}}^3}{3E_{2x}{I}_{2y}}+\frac{h}{G_{1xz}{S}_1}+\frac{d}{G_{2xz}{S}_2}\hfill & \hfill \frac{h_{\mathrm{}}^2}{2E_{1x}{I}_{1y}}\hfill & \hfill -\left(\frac{hd}{G_{1xz}{J}_1}+\frac{d^2}{2E_{2x}{I}_{2y}}\right)\hfill \\ \hfill \frac{h_{\mathrm{}}^2}{2E_{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}{2E_{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 non-linearity 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_c1 and k_c2 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⁻⁵K⁻¹) and ΔT is the temperature change imposed on the bolt. Because k_c1 and k_c2 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 in-plane 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}{3E_{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}{2E_{1x}{I}_{1z}}\hfill & \hfill -\frac{h^2}{2E_{1x}{I}_{1z}}\hfill \\ \hfill -\frac{h_{\mathrm{}}^{2}d}{2E_{1x}{I}_{1z}}\hfill & \hfill \frac{h{d}^2}{E_{1x}{I}_{1z}}+\frac{d_{\mathrm{}}^3}{3E_{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}{2E_{2x}{I}_{2z}}\begin{array}{|c|}\hline +\frac{d}{k_{{\theta }_z}}\\ \hline\end{array}\hfill \\ \hfill -\frac{h^2}{2E_{1x}{I}_{1z}}\hfill & \hfill \frac{hd}{E_{1x}{I}_{1z}}+\frac{d^2}{2E_{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}{3E_{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 quasi-identical 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}{2E_{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}{3E_{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 in-plane angular rigidity k_θZ and two out-of-plane angular rigidities k_θX and k_θY.

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

2.3 A non-linear 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 non-linear 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-\exp \left(-\frac{\left|{\theta }_{rz}\right|}{{\theta }_{\text{ref}}}\right)\right)\text{.}$$(16) k_∞ is the slope of the linear part, M₀ 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 non-linear zone. In Figure 3, the line is the best fit of these parameters to represent the spots coming from the multi-scale 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 quasi-linear 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 semi-rigid 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 best-fit value is 1.615 × 10⁶ N.mm.rad⁻¹.

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

Fig. 4 Spring rotational behavior, kθx and kθy stiffnesses.

3 Numerical and experimental validations

The non-linear 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 in-plane configuration, this element can be modeled by a multi-springs 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 non-linear 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 sub-frames that may be found in a typical bunk bed. Two flexible connections are located at the intersection of beam and column and modeled as non-linear 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 in-plane 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 hyper-static, 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: non-linear 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 set-up 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 post-treated 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³ 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³ N.mm in connections. Consequently, the finite element beam code with non-linear 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 Mz 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 non-linear 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 non-linear 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 force-displacement curves in Figure 8 exhibited an elastic-plastic 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-\exp \left(-\frac{u}{u_{ref}}\right)\right)\text{,}$$(18) where F_init is the y-intercept, K_∞ is the slope in the second linear phase and u_ref is the scaling factor of displacement. Searching the best-fit 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 }}{{ϵ}_{\alpha \alpha }}=\frac{1}{\frac{(\cos {\alpha }^4-v_{tr}\sin {\alpha }^2\cos {\alpha }^2)}{E_t}+\frac{(\sin {\alpha }^4-v_{rt}\sin {\alpha }^2\cos {\alpha }^2)}{E_r}+\frac{\sin {\alpha }^2\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 non-linear wood behavior identification

Using a multi-scale 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 non-linear 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 non-linear 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 multi-scale 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 non-linear 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 non-linear 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

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 mm2 edge for the vertical part and 20 × 140 mm2 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 Mz 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