© L. Chevalier et al., published by EDP Sciences 2019

1 Introduction: from wood product properties to furniture behavior

In furniture industry, the recommendation of the European committee for standardization suggests manufacturing a prototype for each new furniture design. Such a prototype then undergoes a series of mechanical validation tests to ensure its strength and durability [1]. On the basis of the collected experimental results, the robust design and optimization of the prototype is then performed by a trial and error strategy until a compromise is reached. In order to shorten the trial and error design loop, a growing interest has been recently devoted to the development of numerical simulation methods adapted to the furniture industry needs [2]. Basically, three main points have to be addressed to perform virtual tests on furniture:

• characterization of the furniture material properties;

• characterization of the connections behavior between furniture components;

• construction of a numerical model adapted to the geometrical features of furniture for accurate and efficient numerical simulations of timber (wood-based) structures.

Each of these points has already been studied by numerous authors despite the fact that furniture industry began to develop numerical modeling strategies (using computer-aided engineering (CAE) techniques including finite element analysis (FEA) for instance) quite recently compared to other industrial sectors such as the transport (automotive, aeronautic, aerospace, rail) and energy sectors. Such numerical strategies allow the cost incurred by experimental tests and also the time delay before marketing to be significantly reduced due to a lower number of product development cycles. However, many scientific challenges remain to be addressed in order to perform robust and reliable virtual validation tests by taking into account the heterogeneity, anisotropy and variability of wood-based materials and the mechanical behavior (strength and durability) of connections between the structural components of furniture. The interested reader can refer to [3] for a general review of the FEM applied to the analysis of wood-based structures. This work makes reference to 300 research publications on wood analysis including papers and conference proceedings that were published between 1995 and 2004. Since the middle of the 1990s, several scientific contributions on furniture modeling have been published: in [4], a series of performance tests were carried out on side chairs in order to evaluate their strength and durability characteristics according to some acceptance levels of applied loads that were previously determined for a desired category of use. More recently in [5], Eckelman et al. focused on wooden furniture and the structural modeling of the front rails of sofas through performance tests conducted on school chairs made up of round mortise and joints.

The problem related to the modeling of the connection joints between furniture elements has also been addressed in [6–8]. In [6], Kasal et al. investigated the strength properties of glued-dowel jointed sofa frames manufactured from solid wood and wood-based composite materials. As a result of this work, it has been shown that the FEM allows obtaining reasonable estimates of the overall strength performance of the sofa frames, and that the wood-based composite materials could be used instead of solid wood materials for the production of furniture frames. In [7], experimental tests were performed to determine the ultimate shear and bending moment capacities of glued corner blocks under controlled laboratory conditions. Finally, in [8], a model was also developed and compared to experimental results of static compression and tension tests for estimating the bending moment resistance of screw connected L-type corner joints made up of particleboard and medium density fiberboard. All these contributions were essentially based on experimental tests under given load conditions but without considering all 3D loading possibilities on a piece of furniture or even just a furniture part. Another interesting work on joints modeling can be found in [9] which dealt with the analysis of the stress-strain state developed in corner joints for box-type furniture by using FEA.

Later on, a growing interest has been devoted to furniture modeling: in [10] for example, Mishra and Sain carried out 3D finite element simulations of a chair base made up of wood thermoplastic composites under static loading conditions by means of the commercially available SolidWorks (SW) software. Another example can be found in [11] where Çolakoğlu and Apay presented the numerical simulation and strength analysis of three chairs produced from different types of wood in free drop using ANSYS^TM software. A recent collaboration between FCBA (the French Institute of Technology for Forest-based and Furniture Sectors) and the MSME laboratory of Université Paris-Est has produced both experimental and numerical results considering the three following aspects, namely the heterogeneity, anisotropy and dispersion of wood properties [12], with more recently a special focus on the mechanical behavior of connections under static loading [13]. The furniture product concerned in this collaboration is a bunk bed made up of spruce and whose junctions are designed with steel through bolt and dowel nut connectors. Two wooden pins allow maintaining the horizontal part orientation stable (see Fig. 1).

As a result of FEA, design recommendations and suggestions have been provided for furniture designers and manufacturers. Nevertheless, performing more accurate numerical simulations would require a better knowledge of the wood properties under cyclic loading. Furniture design under cyclic loading conditions has been addressed in very few works in the literature until now: several research papers deal with fatigue as in [14–16] for example, but without considering any remaining strain under cyclic loading.

Considering the bunk bed study, Figure 2 presents some numerical simulations carried out on SolidWorks software where a static external force F is applied on the bolt head of the connection. The normalized tests for this specific configuration consist in applying a cyclic loading and controlling the evolution of the gap that may increase with the number of cycles. Focusing on the connection, one can observe two regions of interest where local compression of wood in the transverse direction (that is orthogonal to the wood fibers) generates a permanent strain that must be taken into account in order to predict the gap evolution with respect to the number of cycles. It appears that an in-depth analysis of the cyclic compression behavior of the wood material in the transverse direction is required to model the mechanical behavior of furniture accurately.

Considering the fact that we want to provide a fast tool for evaluating the connections behavior under cyclic loading, we chose to develop a simple homogeneous model. Wood products highlight complex behavior due to heterogeneity, orientation and anisotropy that would require to develop a more complex heterogeneous model.

Nevertheless, the influence of material heterogeneities is taken into account by considering the dispersion via a probabilistic approach that also takes into account the anisotropy in the normal to fibers plane. In real wooden part of furniture it is rather difficult and even impossible to make the difference between the radial and the circumferential orientations. Consequently, the proposed approach is consistent with the degree of uncertainties on the knowledge of the wood properties.

The paper is organized as follows. Section 2 presents the experimental test campaign performed on spruce cubic specimens. Cyclic compression tests have been performed with a fixed value of the maximum pressure during ten thousand cycles. The gap evolution has been modeled by using a decreasing exponential law with respect to the number of cycles, and a proportional law with respect to the local contact pressure. In Section 3, as the dispersion due to the strong heterogeneity of wood material in the transverse orientation (in relation to the compression direction) is important, we consider a stochastic modeling of the cyclic behavior of wood, as already done in [17,18], through the construction of an ad hoc probabilistic model for each parameter modeled as a random variable. In Section 4, a simplified pressure distribution is proposed to model the mechanical behavior of the corner of a bunk bed that has already been considered and studied under static loading in [13]. In that way, one can evaluate the gap evolution with the number of cycles during the cyclic tests performed on the piece of furniture under consideration (bunk bed). The mean model for the gap evolution in the bunk bed corner has been numerically simulated and used to study the influence of geometrical parameters of the bunk bed corner model on the gap after ten thousand cycles. Taking into account the parametric probabilistic model of the cyclic behavior of wood constructed in Section 3, one can generate independent realizations of the underlying random variables (that are the uncertain model parameters) with the identified probability distributions and perform a stochastic analysis of the gap evolution in the bunk bed corner using a Monte Carlo stochastic solver. Finally, we compute the first- and second-order statistical moments of the random gap in the bunk bed model and we construct a confidence region corresponding to a given probability level for robust design purposes.

Fig. 1 Test bunk bed and focus on the junction and connecting elements.

Fig. 2 (Top) Numerical simulation of the bolt head of a representative structure of the side of the bunk bed submitted to a static load F = 300 N and (bottom) focus on the two regions of interest located under the bolt head and at the bottom of the contact zone.

2 Mechanical behavior of spruce under cyclic compression

Wooden materials, even in the proportionality domain, corresponding to a linear elastic behavior, highlight a singular response under cyclic compression loading. After each compression cycle, a small amount of residual penetration remains and consequently, the penetration increases as the number of cycles increases. Let us now describe the test specimens and detail the experimental configuration of the cyclic compression tests.

Small sized cubic specimens of spruce (Picea sp. with an approximate mean density of 440 kg·m⁻³ and a standard deviation of 27.8 kg·m⁻³) were cut from different wood parts especially prepared to be used in the manufacturing of bunk beds provided by a furniture company. The specimens were cut with dimensions of approximately 57 mm × 57 mm × 57 mm and were chosen to be as clear as possible (without any visible defect). The wood fibers direction was estimated to be almost parallel to four of the six sides of the cube for all specimens. The orientation of the compression is normal to the direction of the wood fibers and the compression load is applied directly on the wood surface using an ad hoc compression cylindrical tool. All specimens have been stored in the same environmental conditions at a temperature of 20 ± 2 °C and a relative humidity of 65 ± 5% for approximately 6 months. The moisture content of the specimens was approximately 12%.

The specimens were tested with a universal testing machine Deltalab500 provided with a load cell of 5 kN capacity. The force is applied with a displacement rate of ±0.2 mm·min⁻¹ during 10 000 cycles. The applied force F(t) varies from 0 to F_max with F_max= 300, 400, 600 and 800 N. The associated penetration u(t) corresponds to the displacement measured from the position of the cylinder compression head. Graph of Figure 3a shows the first 25 cycles of time function F(t) with F_max= 300 N, and Figure 3b represents the evolution of force F(t) as a function of penetration u(t) in order to characterize the wood behavior under cyclic local compression. First, one can see a preliminary accommodation stage that can also be observed in static compression tests, followed from the second cycle by a regular evolution with an increasing penetration at each cycle. No limit cycle is reached and the penetration increases regularly and gradually, but one can also observe that the behavior becomes stiffer when the load increases but becomes softer when the load decreases, thus highlighting a typical nonlinear elastic behavior.

After a first cycle with a penetration of nearly 1 mm, displacement u(t) does not come back to 0 and one can observe an accumulation of penetration as the cyclic test continues. In the following, we are interested in the evolution of residual penetration δ_r after the first cycle, that is starting from time t₀ corresponding to displacement u(t₀)= 0.34 mm (represented by a circle in Fig. 3b). For each experimental test we determine the envelope of function u(t)−u(t₀) by computing the minimal and maximal values of u(t)–u(t₀) reached during each cycle as shown in Figure 4. Thus, we construct the graph of residual penetration δ_r versus number of cycles N corresponding to the lower (red) curve in Figure 4 (the upper (cyan) curve will not be considered in the following). The residual penetration evolves almost linearly after several cycles as it can be observed in Figure 4 (where only the first 1700 cycles of a representative cyclic test are plotted for a better visualization of both phases). It can then be modeled as an affine function of the number of cycles after a transition stage.

The same cyclic compression test has been carried out with a maximum load of 600 N and shows a higher residual penetration than in the case of 300 N with a greater slope in the transition phase. Conversely, in the regular (almost affine) part, the evolution of penetration versus the number of cycles is similar in both test cases as shown in Figure 4a and b. A measure of the slope in the regular (almost affine) part would then provide more information on the cyclic behavior of wood.

Each test cube has been photographed, the orientation θ of the wood rings with respect to the load direction has been measured as shown in Figure 5a, and the density index d has been characterized by counting the number of rings divided by the face area. We performed around 15 cyclic compression tests for each value of maximum force F _max resulting in a complete campaign of 64 tests. Each compression test consists of performing 10 000 cycles and each cycle duration is about 3 sec, consequently each test is about 8 h long and the complete campaign lasts around 2 months.

It is worth noticing that, for the same series of tests, different typologies of behavior can be identified. Indeed, Figure 6 shows that for a fixed value of maximum load F _max = 300 N, most experimental tests give results similar to the curve displayed in Figure 6a showing an almost linear increasing behavior, but some of them present a flat or even decreasing behavior (see Fig. 6b and c), while others show an almost linear increasing evolution of the residual penetration with a change of slope sometimes beyond several thousand cycles (see Fig. 6d). These various evolutions are rather difficult to explain especially the one displayed in Figure 6b where the residual penetration starts to decrease after 3000 cycles, which means that the material recovers strength as the cycles continue. A viscoelastic behavior would probably explain such a material growth under compression loading.

Considering the most representative evolution of the residual penetration shown in Figure 6a, the resulting curve can be fitted with an affine function weighted by an exponentially decreasing function of the number of cycles. The residual penetration δ_r is then defined as a function of the number of cycles N under the following form:(1)where δ ₀ is a virtual residual penetration defined as the ordinate at the origin (i.e. N = 0) of the affine function modeling the regular (almost linear) second part of the curve, k _∞ = Λδ ₀ denotes the slope of the regular (almost linear) part of the curve, and N _ref is the number of cycles that characterize the transition (nonlinear) first part of the curve. Precisely, N _ref is the number of cycles where the residual penetration reaches 62% of δ ₀(1 + ΛN). The mean value of this parameter is N _ref = 50 cycles that is a very small value compared to the 10 000 cycles of the compression test. Consequently, the fluctuations of this parameter have a very low impact on the cyclic behavior and will not be studied in the following. N _ref is assumed to be equal to 50 for all simulations.

Figure 8 displays the experimental values identified for both parameters δ ₀ and k_∞ for the eleven tests performed with a maximum load F _max = 300 N. The two top (resp. bottom) figures correspond to the values of parameter δ ₀ (resp. of slope k _∞) with respect to ring orientation θ and density index d. Both parameters δ ₀ and k _∞ do not seem to be correlated with the ring orientation or the density index and present an important dispersion of 54% (resp. 68%).

Table 1 reports the second-order statistical moments (mean value, standard deviation and dispersion) of both model parameters δ ₀ and k _∞ for the four considered values of F _max. Mean values and standard deviations are given for each series of tests and one can see that the dispersion illustrated in Figure 7 for F _max = 300 N can reach even higher values for the series of tests performed at F _max = 400 or 600 N. Considering the mean values, Figure 8 shows that the evolutions of both parameters δ ₀ and k _∞ with respect to maximum load F _max can be modeled as linear functions of F _max: δ ₀ = αF _max and k _∞= βF _max, where α and β are some strictly positive constants independent of F _max. Since both parameters δ ₀ and k _∞ are assumed to vary linearly with the maximum applied load F _max during each cycle, the parameter Λ = k _∞/δ ₀ introduced in equation (1) can be considered as a constant value that is independent of maximum load F _max and equal to β/α. Note that to go further in the exploitation of the proposed model, it would be useful to study the evolution of the model parameters with respect to the local pressure and not to the global force.

Considering the compression cylinder dimensions (cylindrical tool of 15 mm diameter) and assuming that the contact pressure is uniform, the maximum contact pressure p during a cycle can be determined for all tests and parameter δ ₀ can be written as:(2) where S denotes the cylinder contact area. Considering the high dispersion of the values of δ ₀ and k_∞ between all tests, both parameters α and Λ will be considered as uncertain and modeled by real-valued random variables in Section 3 to take into account the variabilities in the cyclic behavior of wood. Using equation (2), one can easily derive the local relation of proportionality between pressure p and residual penetration δ_r at each cycle in the following form:(3)

This proportionality relation between the pressure and the residual penetration will be used in the bunk bed corner model presented in Section 4.

Fig. 3 Cyclic compression test on the bolt head. (Left) Evolution of force F(t) with respect to time t during the first 25 triangular loading cycles, where the force does not come back exactly at 0 N but at 10 N in order to avoid rupture and contact problems that could be not detected by the load cell. (Right) Evolution of force F(t) with respect to displacement u(t), showing a preliminary phase (similar to a static compression test) during the first cycle, followed by more regular cycles (the model focuses on the regular cycles).

Fig. 4 Global response of the cyclic compression test carried out on the bolt head: evolution of penetration u(t)–u(t 0) with respect to the number of cycles when F(t) varies from 0 to 360 N and back following a triangular signal for a maximum load F max = 300 N (left) and F max = 600 N (right). The upper (cyan) curve corresponds to the maximum penetration when the maximum load is applied, while the lower (red) curve corresponds to the residual permanent penetration δr when no more force is applied.

Fig. 5 (Left) Measurement of the density index d and rings orientation θ on a face for a given test cubic specimen. (Right) Cyclic local compression managed through the cylinder head contact surface. The load cell applies regular compression load cycles.

Fig. 6 Different evolutions of the residual penetration (or remaining depth) with the number of cycles for the series of tests with F max = 300 N after 4000 cycles: (a) an almost linear increasing response (corresponding to most experimental results); (b) a flat then decreasing response; (c) a constant flat response; (d) an almost linear increasing response with a change of slope after 3000 cycles.

Fig. 7 Model definition and parameterization of the residual penetration for a representative cyclic compression test of 10 000 cycles: the almost linear increasing part is modeled by an affine function δ 0 + k ∞ N, where k ∞ is the slope and δ 0 the ordinate at the origin N = 0, and the nonlinear transition part is characterized by a characteristic number of cycles N ref.

Fig. 8 Dispersion of the experimental results obtained from the eleven cyclic compression tests performed with F max = 300 N: the 2 top figures show the identified values of δ 0 versus the ring orientation (left) and the density index (right), while the 2 lower figures show the ones of slope k ∞.

3 Stochastic modeling of the cyclic behavior of wood in the direction normal to fibers

To take into account the uncertainties and variabilities in the experimental results, each model parameter is represented as a second-order real-valued random variable whose probability distribution must be defined and constructed. In our case, the model parameters are the ratio Λ = k _∞/δ ₀ between k _∞ and δ ₀ and the one α = δ ₀/F _max between δ ₀ and F _max. As these two parameters have very different origins and there is a priori no available information about some statistical dependence between them, we further consider that they are mutually independent random variables.

Figure 9 shows all the experimental data for parameters α and Λ obtained from the 64 cyclic compression tests as well as independent realizations of the corresponding random variables which have been generated from the probability distributions constructed in Section 3.2. Table 2 gives the mean value, standard deviation and dispersion of random variables α and Λ estimated by using the 2000 independent realizations represented on the right side of Figure 9. The probability distribution of each random variable is assumed to be represented by a probability density function (pdf) whose construction is briefly described in the next paragraphs.

Fig. 9 Dispersion of model parameters α and Λ on the 64 cyclic compression tests: (left) data obtained from the experimental tests; (right) 2000 independent realizations of random variables α and Λ generated from the identified probability distributions.

3.1 Maximum entropy principle

Considering a second-order real-valued random variable X defined by the pdf p_X (x), the entropy S(p_X ) of p_X is given by:(4) S(p_X ) represents a measure of the uncertainties associated to the pdf p_X of X. Initially introduced within the framework of information theory [19–21], the maximum entropy (MaxEnt) principle [22–25] is a general and efficient optimization procedure that allows a parametric representation of the pdf of random variable X to be constructed by maximizing the entropy S(p_X ) under a set of constraints defined by the available information on X [26–28] such as the support of the pdf, the existence of the mean value, the standard deviation or higher order moments related to X for example. A Lagrange multiplier λ_i is introduced and associated with each constraint defined by the available information. In addition to the normalization condition, these constraints are generally related to some statistical properties of X and typically given under the following form:(5)where E denotes the mathematical expectation, g_i are given mappings defined from ℝ into ℝ and f_i are given values in ℝ. For instance, if the mean value of X is assumed to be known, then g_i (x) = x and f_i corresponds to the mean value E{X} of X. It can be shown that the Lagrange multipliers λ_i are the solution of a convex optimization problem that consists in minimizing the strictly convex function H defined by(6)where λ ₀ is the Lagrange multiplier associated to the normalization condition, the interval [a,b] denotes the support of the pdf p_X of random variable X and 1_[a,b](x) is the indicator function of interval [a,b] that is equal to 1 if x belongs to [a,b] and 0 otherwise. The pdf p_X of X is then parameterized by the (m+1) Lagrange multipliers λ₀, λ₁, …, λ_m and explicitly given by(7)

For instance, assuming that X is a second-order real-valued random variable whose support is ℝ and whose mean value and standard deviation are given, the MaxEnt principle leads to a Gaussian pdf for p_X . In the next paragraph, we will address the construction of the pdf associated with each of the random variables Λ and α involved in the stochastic cyclic compression problem.

3.2 Probability density functions of the random variables in the stochastic modeling of cyclic behavior of wood

In this section, we construct parametric representations of the pdfs of both second-order random variables Λ and α by having recourse to the MaxEnt principle. The mean value of random variable Λ is denoted by m _Λ, and the dispersion of Λ around m _Λ is characterized by its standard deviation σ _Λ. From a mathematical viewpoint, random variable Λ cannot take negative values and there is a priori no upper bound for the values of Λ, so that the support of Λ is assumed to be ℝ ⁺. Finally, the constraints defined by the available information on Λ and integrated in the MaxEnt formulation for the construction of the pdf p _Λ(λ) of Λ are given by(8)

The Lagrange multipliers λ ₀, λ ₁ and λ ₂ are numerically computed by minimizing the strictly convex function H(λ ₀ ,λ ₁ ,λ ₂ ) defined by equation (6), leading to λ  ₀ = 4.49, λ ₁= 0.011 and λ ₂ = 1.86 × 10⁻⁶. The corresponding identified pdf p _Λ (λ) of random variable Λ is represented in Figure 10. Obviously, this pdf does not correspond to a Gaussian pdf as it is not symmetric and its support is not the real line ℝ since pdf p _Λ (λ) is strictly equal to 0 for negative values of λ.

Knowing explicitly the pdf of random variable Λ, one can easily construct the inverse cumulative density function (cdf) of Λ and generate independent realizations of Λ (see Fig. 9) using the inverse transformation sampling method and a uniformly distributed (pseudo-)random number generator: from independent realizations of a uniform real-valued random variable between 0 and 1 (drawn from the aforementioned uniform random generator), the inverse cdf of Λ allows for computing independent realizations of Λ that are consistent with the experimental results. We can also proceed in the same manner to construct the pdf of random variable α by making use of the MaxEnt principle with the same constraints as for random variable Λ. In this case, minimizing the strictly convex function H(λ ₀ ,λ ₁ ,λ ₂ ) yields the three following values for λ ₀, λ ₁ and λ ₂: λ ₀ = 1.09, λ ₁ = −12.6, λ ₂ = 19.7. Figure 10 shows the pdf of both random variables Λ and α.

Let us note that this stochastic modeling based on the MaxEnt formulation leads to pdfs for Λ and α that are very different. For random variable Λ, the dispersion is very high and the shape of the pdf is a strictly monotonic exponentially decreasing function that is similar to an exponential pdf and very far from a Gaussian pdf. For random variable α, the dispersion is smaller than in the case of random variable Λ and the shape of the pdf looks like a Gaussian function truncated on ℝ ⁺. In the following section, we will use the proposed stochastic modeling of uncertainties to perform Monte Carlo numerical simulations of a bunk bed corner problem in order to analyze the random gap evolution in the bunk bed corner and provide statistical information (such as statistical moments and confidence regions) related to the random gap evolution after completing ten thousand loading cycles.

Fig. 10 Probability density functions of random variables Λ (left) and α (right).

4 Modeling the gap evolution in the bunk bed corner

In this section we present a simplified model describing the pressure distribution in the two regions of interest of the bunk bed represented in Figure 2 and allowing the gap evolution to be predicted. First, the mean model for the gap evolution in the bunk bed corner under cyclic bending loading will be presented and then, the statistical properties of the random gap evolution are analyzed by using the stochastic modeling presented in Section 3 and the Monte Carlo numerical simulation method.

4.1 Mean model for the gap evolution in the bunk bed corner under cyclic bending loading

We consider the bunk bed corner problem illustrated in Figure 11, in which a vertical part with a square cross section (of side length 57 mm) is linked with only one or two fixing components corresponding to through-bolt and dowel-nut connectors. As previously mentioned, two wooden pins prevent the horizontal part from rotating around the y axis but they will not be considered in the model for the sake of simplicity.

When applying a downward vertical load at the right end of the horizontal part, one generates a bending momentum around the z axis that passes through the junction. The right side of Figure 11 illustrates the bending effect on the relative displacement of both wooden parts. This relative displacement can be described using only two parameters δ and θ, respectively equal to the penetration of the bolt head in the vertical part and the relative rotation between the two wooden parts. One can see that this relative displacement generates some interferences between the two components in a localized area, where strain occurs in both wooden parts and in particular in the vertical part where the contact pressure is applied perpendicularly to the wood fibers direction (just like in the cyclic compression tests presented in the previous sections). In these zones, the local pressure p(x) located at point M can be linked to the local residual penetration δ_r(x) using the simple law p(x)= Kδ_r(x) presented in Section 2.

If the bolt is assumed to be rigid enough compared to both wooden parts, one can postulate that the penetration δ of the bolt head is equal to the gap between the two parts at point O. Consequently, the local penetration δ_r(x) can be expressed in terms of the vertical position x as:

(9)

There is contact between the two wooden parts only in case of effective penetration (i.e. for δ_r(x) > 0) and consequently, the length Δ of the loaded contact zone goes from x = δ/θ to x = e. We then denote by a, the length defined by the ratio a = δ/θ. This leads to the following expression for the length Δ of the loaded zone:(10)Equation (10) indicates that a is the length of the free surface between the bolt and the contact zone. Considering the relation between residual penetration δ_r(x) and local contact pressure p(x), one can model the pressure by a linear distribution near C and by a uniform distribution near A. The force balance equation of the bolt and the momentum balance equation of the traverse (horizontal wooden part) lead to the following set of equations:(11) (12) (13)where S_b is the contact area under the bolt head. Bending momentum M can be determined by performing a global numerical simulation of the corner problem, while δ, θ, Δ and F are unknown parameters (Fig. 12).

Using equations (11) and (12) and introducing length a lead to:(14)

Then, considering equation (10) we obtain:(15)Solving this 2nd degree polynomial equation, one can derive the expression of length a as follows:(16)

It is worth noticing that length a depends neither on the rigidity K which varies during cycles, nor on the load amplitude F, hence a is fully deterministic. As all terms S_b , e and b are constant geometrical parameters, the ratio a between penetration δ and rotation θ remains also constant during cycles.

We then use equation (10) to obtain Δ, equation (13) to obtain F, equation (11) to obtain δ and finally the definition of a is used to calculate the rotation θ. This procedure can easily be implemented as a post-processing of the finite element computation of the bending momentum in each corner and allows the evolution of the gap δ in the connection between both wooden parts to be predicted.

For example, Figure 13 shows the results obtained by performing a global numerical simulation of a normalized cyclic test on the bunk bed corner with a bending momentum M = 20 N·m, leading to a tension F = 306 N in the bolt, and with the following geometrical features: 2e = 148 mm, b = 25 mm and S_b= 148 mm². One can notice that in this case, length a = 50 mm so that the length of the contact zone between the horizontal traverse and the vertical part of the corner is Δ = 24 mm. It is also worth noticing that the gap δ between the two wooden parts does not exceed 0.22 mm after 10 000 cycles, which allows for a safe design for the mean model of the bunk bed corner.

Considering the example presented in Figure 13 as the reference case, one could be interested in studying the influence of geometrical parameters, such as the height 2e or the thickness b of the traverse, on the gap δ. Figure 14 represents the evolution of gap δ with respect to thickness b (with a fixed reference value for e) and also with respect to half height e (with a fixed reference value for b). One can see that both geometrical dimensions have a reducing influence on the value of the gap after 10 000 cycles, which means that the gap δ decreases as the thickness b or the half height e increases. The variation of the gap with respect to b is small compared to the one with respect to e. But even reducing drastically this last geometrical parameter, the gap remains less than 1 mm which seems to be secure enough according to the technical specifications for the considered piece of furniture. Besides, the design of the connection could be changed by adding a clamping washer of diameter D under the bolt head, so that the contact area S_b under the bolt head would be bigger and the contact pressure lower. Figure 15 shows the evolution of the gap δ as the bolt diameter D increases from 15 to 50 mm. One can see that increasing the bolt diameter D (or equivalently, the contact area S_b under the bolt head) allows reducing the gap evolution, since the length Δ of the contact zone increases and therefore the length a decreases and the maximum contact pressure p _max becomes lower.

Comparing the graphs of Figures 14 and 15, one can conclude that the most part of this gap reduction is mainly due to the contact area S_b under the bolt head and the height 2e of the traverse. Increasing the thickness b has a lesser impact on the gap. Quantifying such influences on the gap is of great importance for the design of furniture. Furthermore, in order to enhance the prediction of the gap within the context of robust design under uncertainties, we will perform a stochastic analysis of the random gap in accordance with the stochastic model derived in Section 3 in order to compute estimates of the mean value and standard deviation and to provide a confidence region (corresponding to a given probability level) for the random gap.

Fig. 11 Connection between two parts of the bunk bed (a vertical part with square cross section and a horizontal part of thickness b in the z direction): (left) initial configuration and (right) deformed configuration due to the bending momentum in the junction resulting in a penetration near point C.

Fig. 12 Balance of the horizontal wood part gives: (i) a relation between the maximum contact pressure p max and the tension F in the bolt: F = p max bΔ/2 with p max = K(θ e −δ) = K θΔ and (ii) a relation between the bending moment M in the connection and the tension F in the bolt: F = M/(e−Δ/3).

Fig. 13 Evolution of the gap δ and the relative rotation θ with respect to the number of cycles during a cyclic loading at a constant bending moment M = 20 N·m corresponding to a tension F = 306 N in the bolt.

Fig. 14 Influence of geometrical parameters on the gap: half height e and thickness b of the traverse have a reducing influence on the gap δ.

Fig. 15 Influence of the external diameter under the bolt head on the gap: using a clamping washer of diameter D = 30 mm allows reducing the gap δ by a factor four.

4.2 Stochastic analysis of the random gap evolution

In this section, we consider the stochastic modeling presented in Section 3 to take into account the uncertainties in the cyclic behavior of wood and we analyze the propagation of these uncertainties through the furniture corner model under cyclic bending loading using the Monte Carlo numerical simulation method [29–33] as stochastic solver. For this purpose, we carried out 2000 numerical simulations of the stochastic corner model and focused on the random gap δ _10 000 corresponding to the random residual penetration under the bolt head after 10 000 cycles.

A convergence analysis of the mean value and standard deviation of the random gap δ _10 000 has been performed with respect to the number of independent realizations used by the Monte Carlo stochastic solver. Figure 16 shows that the convergence is reached after around 500 realizations for the first two statistical moments. The converged estimation of the mean value of the random gap δ _10 000 is found to be equal to 0.26 mm which is quite close but slightly higher than the deterministic value of 0.22 mm obtained for the deterministic mean model.

The set of 2000 independent realizations of the random gap δ _10 000 computed to perform the convergence analysis for the first- and second-order statistical moments can also be used to construct an estimation of the probability density function of this random variable.

Figure 17 shows the shape of the underlying pdf in the form of a histogram built from an appropriate number of bins (here, 30 bins) to cover the range of values in the interval [δ _min, δ _max], where δ _min and δ _max are the minimum and maximum values taken by the random variable δ _{10 000} over the set of 2000 independent realizations computed by using the Monte Carlo stochastic solver. Smoother representations of the pdf may be obtained by using kernel density estimation (smoothing) techniques [34–36]. Considering the histogram of Figure 17, it is possible to define a confidence region corresponding to a given probability level p_c (for example, p_c  = 90%). The upper and lower bounds of this confidence region are defined by:(17)where ζ(p) is the quantile function defined for a real-valued random variable X with cdf F_X(x) by:(18)

The computation of such bounds can simply be done by sorting the realizations in increasing order and counting the number of realizations until the cdf reaches p_c for the upper bound and (1–p_c ) for the lower bound. In the case of the pdf of the random residual penetration δ _10 000 after 10 000 cycles shown in Figure 17, the upper and lower bounds of the 90% confidence interval are such that δ _10 000 ⁺ = 0.65 mm and δ _{10 000} ⁻ = 0.05 mm. It appears that the mean value of 0.26 mm is included but not centered in the confidence interval.

Finally, if we consider the span of the 90% confidence interval in comparison with the value obtained by using the deterministic mean model, we can conclude that the gap after completing 10 000 cycles may reach 0.65 mm that is almost three times higher than the expected value of 0.22 mm predicted by the deterministic mean model. This conclusion highlights the relevance of the proposed stochastic approach for modeling the gap occurring in the connections of furniture undergoing mechanical cyclic loading.

Fig. 16 Convergence analysis of the mean value and standard deviation of the random residual penetration after 10 000 cycles with respect to the number of realizations.

Fig. 17 Probability density function of the random response δ10 000 represented by a histogram.

5 Conclusions

The static response of furniture under mechanical loading can be numerically simulated using FEA but the large dispersion due to the strong heterogeneity of wooden materials necessitates to perform a stochastic analysis of furniture behavior by taking into account the uncertainties in the experimental configuration and between the different test specimens. In case of cyclic loading, the study is even more complicated considering the nonlinear behavior of the material under a repeated compression loading.

In this study, we first performed a series of cyclic compression tests on spruce specimens to identify the evolution of the residual penetration with respect to the number of cycles under different maximum loads. The experimental results show that the global response of the cyclic compression tests is quite regular but exhibits a large dispersion from one specimen to another. Based on these experimental curves, we have proposed a simple analytical model to describe the evolution of residual penetration versus the number of cycles. To take into account the uncertainties in the cyclic behavior of wood, we then have constructed a probabilistic model for the uncertain parameters modeled by statistically independent real-valued random variables using the maximum entropy (MaxEnt) principle.

Finally, the proposed probabilistic model has been used to study the mechanical behavior of a typical connection between two wooden parts of a bunk bed. An analytical model of the bunk bed corner has been built and used to investigate the influence of geometrical features on the gap evolution under cyclic loading within the deterministic framework. Subsequently, a stochastic analysis of the random gap has been carried out using the Monte Carlo numerical simulation method as stochastic solver, thus allowing the mean value and the standard deviation as well as a 90% confidence region for the random gap to be estimated. It has been shown that the natural dispersion in the wood behavior under cyclic loading leads to realizations of the random gap that can reach three times the value predicted by the deterministic mean model.

References

All Tables

All Figures

	Fig. 1 Test bunk bed and focus on the junction and connecting elements.
In the text

	Fig. 2 (Top) Numerical simulation of the bolt head of a representative structure of the side of the bunk bed submitted to a static load F = 300 N and (bottom) focus on the two regions of interest located under the bolt head and at the bottom of the contact zone.
In the text

Fig. 3 Cyclic compression test on the bolt head. (Left) Evolution of force F(t) with respect to time t during the first 25 triangular loading cycles, where the force does not come back exactly at 0 N but at 10 N in order to avoid rupture and contact problems that could be not detected by the load cell. (Right) Evolution of force F(t) with respect to displacement u(t), showing a preliminary phase (similar to a static compression test) during the first cycle, followed by more regular cycles (the model focuses on the regular cycles).

In the text

Fig. 4 Global response of the cyclic compression test carried out on the bolt head: evolution of penetration u(t)–u(t 0) with respect to the number of cycles when F(t) varies from 0 to 360 N and back following a triangular signal for a maximum load F max = 300 N (left) and F max = 600 N (right). The upper (cyan) curve corresponds to the maximum penetration when the maximum load is applied, while the lower (red) curve corresponds to the residual permanent penetration δr when no more force is applied.

In the text

	Fig. 5 (Left) Measurement of the density index d and rings orientation θ on a face for a given test cubic specimen. (Right) Cyclic local compression managed through the cylinder head contact surface. The load cell applies regular compression load cycles.
In the text

Fig. 6 Different evolutions of the residual penetration (or remaining depth) with the number of cycles for the series of tests with F max = 300 N after 4000 cycles: (a) an almost linear increasing response (corresponding to most experimental results); (b) a flat then decreasing response; (c) a constant flat response; (d) an almost linear increasing response with a change of slope after 3000 cycles.

In the text

	Fig. 7 Model definition and parameterization of the residual penetration for a representative cyclic compression test of 10 000 cycles: the almost linear increasing part is modeled by an affine function δ 0 + k ∞ N, where k ∞ is the slope and δ 0 the ordinate at the origin N = 0, and the nonlinear transition part is characterized by a characteristic number of cycles N ref.
In the text

	Fig. 8 Dispersion of the experimental results obtained from the eleven cyclic compression tests performed with F max = 300 N: the 2 top figures show the identified values of δ 0 versus the ring orientation (left) and the density index (right), while the 2 lower figures show the ones of slope k ∞.
In the text

	Fig. 9 Dispersion of model parameters α and Λ on the 64 cyclic compression tests: (left) data obtained from the experimental tests; (right) 2000 independent realizations of random variables α and Λ generated from the identified probability distributions.
In the text

	Fig. 10 Probability density functions of random variables Λ (left) and α (right).
In the text

	Fig. 11 Connection between two parts of the bunk bed (a vertical part with square cross section and a horizontal part of thickness b in the z direction): (left) initial configuration and (right) deformed configuration due to the bending momentum in the junction resulting in a penetration near point C.
In the text

	Fig. 12 Balance of the horizontal wood part gives: (i) a relation between the maximum contact pressure p max and the tension F in the bolt: F = p max bΔ/2 with p max = K(θ e −δ) = K θΔ and (ii) a relation between the bending moment M in the connection and the tension F in the bolt: F = M/(e−Δ/3).
In the text

	Fig. 13 Evolution of the gap δ and the relative rotation θ with respect to the number of cycles during a cyclic loading at a constant bending moment M = 20 N·m corresponding to a tension F = 306 N in the bolt.
In the text

	Fig. 14 Influence of geometrical parameters on the gap: half height e and thickness b of the traverse have a reducing influence on the gap δ.
In the text

	Fig. 15 Influence of the external diameter under the bolt head on the gap: using a clamping washer of diameter D = 30 mm allows reducing the gap δ by a factor four.
In the text

	Fig. 16 Convergence analysis of the mean value and standard deviation of the random residual penetration after 10 000 cycles with respect to the number of realizations.
In the text

	Fig. 17 Probability density function of the random response δ10 000 represented by a histogram.
In the text