A stochastic approach for the evaluation of the reliability of wood furniture in an industrial context: managing virtual standardization tests

– In this paper, a procedure is developed in order to manage virtual tests using numerical simulation in order to reduce time and to optimize design before the real validation tests on wood products. This procedure in 3 steps can be used in an industrial context: (i) characterization of the material properties, (ii) construction of probabilistic density function for each parameter, (iii) Monte Carlo simulation of the structural behavior of the wood products furniture. For this purpose, experimental tests associated with digital image correlation technique were managed using 3 point bending tests easy to develop in an industrial environment. A probabilistic density function was constructed for each mechanical parameter based on the test results using the principle of maximum entropy. The simulation is managed by a FEA developed in the Matlab environment and allows a Monte-Carlo simulation in a reasonable CPU time. We can then evaluate the eﬀect of the input data dispersions on the overall reaction of the structure. A real test is performed on a prototype product in order to conﬁrm the simulation results and to validate the approach.


Introduction
For the furniture industry, CEN the European committee for standardization has approved a number of tests on furniture to ensure its strength and durability before authorization can be given to sell the product.To pass through these tests, a prototype has to be manufactured for each new design and submitted to validation tests.On the basis of the results, the design is evaluated and modified.This trial and error process continues until a compromise is reached.Numerical simulation could be an interesting alternative to this expensive process.Mainly, three points have to be addressed to manage virtual tests on furniture: characterization of the furniture materiel property; characterization of the connections behavior between furniture components; focuses on strength characteristics according to the accepted levels of applied loads that had been determined for a desired category of use.More recently [3,4] focused on wooden furniture and the structural modeling of the front rails of sofas or the performance tests of school chairs constructed with round mortise and joints.
The problem of the link between elements of furniture was addressed and Kasal et al. [5][6][7] investigated the strength properties of glued-dowel jointed sofa frames constructed with solid wood and wood-based composite materials.As a result of his work, Kasal found that FEM gives reasonable estimates of the overall strength performance of the sofa frames; it was also concluded that the wood composite materials could be used instead of solid wood materials in the production of frame constructed furniture.He also worked on experimental tests to determine the ultimate shear and bending moment capacities of glued corner blocks under controlled laboratory conditions.Finally a model for estimating the moment resistance of L-type screw corner joints constructed of particleboard and medium density fiberboard was also published.All these contributions were essentially experimental and could not test all 3D loading possibilities on the furniture or the part of the furniture.
Meanwhile Nicholls and Crisan [8] analyzed the stressstrain state found in corner joints and box-type furniture by using FEA.FEA was beginning to be widely used in the furniture modeling context later than in automobile or aerospace applications.For example Mishra and Sain [9] used 3D finite element simulation of static loading on a chair base made up of wood thermoplastic composites by means of the commercially available software SolidWorks.Paoliello et al. [10] worked on defining loads that occur on a Eucalyptus wood chair while sitting different postures and he used ergonomic concepts.Lately, C ¸olakoglu and Apay [11] have presented the simulation and strength analysis of three chairs produced from different types of wood in free drop using ANSYS TM software.As a result of FEA, suggestions were developed for furniture designers and manufacturers, but accurate simulation needs a fine knowledge of the wood properties.
Wood is known as an inhomogeneous material with highly variable mechanical properties [12].Thereby, these properties can be considered as random variables that can be characterized by probability distribution.The general procedure for selecting the appropriate distribution of the parameter in question depends on selecting the statistical distribution that best fits to the empirical data set [13,14].Based on assignments for timber by [15], the normal, lognormal and Weibull distributions are considered as the most appropriate family of probability density function (pdf) to model the clear wood properties.However, it can be noted that strength and stiffness cannot attain negative value, so that normal distributions may be relatively unable to represent the physical sense of these parameters [16].
Consequently, managing virtual test by using the FEA that takes into account the dispersion of wood can provide accurate and numerous anisotropic material data via a Monte Carlo simulation and gives a confidence region for the furniture's rigidity and strength.This procedure has already been used successfully in different contexts: (see for example [17][18][19][20]).In order to manage fast simulations in the Monte Carlo process, we chose to model the structure of the bunk bed with anisotropic beam elements in order to reduce drastically the number of Degrees Of Freedom (DOF For short) in regard of a complete 3D FEA.Considering the important difference between the longitudinal and transverse properties, the natural beam theory is developed taking into account the shear effect.The use of beam element reduces the accuracy of the stress field especially in corners and connections between beams.
One way for approximating the effects of connection flexibility is to use extensional and rotational springs.The stiffness of each spring is the ratio of the transmitted force to the corresponding relative displacement within the connection.This system of springs can be represented with a zero length element 0D with two coincident nodes and thus assembled in a matrix form as a finite element which makes the integration of the connection behavior in FEA code relatively easy at programing level.Calculations showed that the stiffness of rotational components are significantly lower than the extensional ones so only the contribution of rotational springs is considered and the effects of extensional springs can be ignored.The identification of the connection stiffness matrix components is done by a multi-scale numerical approach, using 3D FE simulations.Details can be found in reference [21].
The studied bunk bed is made up of spruce.Measurements of elastic characteristics and strength of the material have required characterization tests; 3 point bending test associated with the technique of digital image correlation (DIC).This technique is well known since the early 80's [22][23][24][25] and applied in a lot of fields to identify material behavior, in particular for flexible material like polymers or rubbers where strain gages cannot be used.Usually pictures are taken during plane stress tests using 2D images [26,27] or even 3D structures [28].In this paper, characteristics are measured from 3 point bending tests.
The classical method for Young's modulus E and shear modulus G measurement by static bending test for wood beams needs to provide tests by varying the depth/span ratios of two or more beam specimens and is well known since early 1960's.Three point bending test results usually underestimate the value of G due to stress concentration imposed by the loading point which is not taken into consideration in beam theory.To overcome this problem and by means of finite element simulation, a correction factor on the Timoshenko's shear coefficient was proposed by [29] and plausible shear modulus was obtained.Murata and Tsubasa [30] managed to determine simultaneously both E and G by polynomial regression analysis of deflection curves.DIC was recommended as a powerful tool for full-field normal and shear strain measurement.Yoshihara and Tsunematsu [31] examined the feasibility of estimation methods for measuring Young's modulus of wood by 3 point bending test.They found that shear Fig. 1.Three principal axes of wood with respect to grain direction and growth rings [12].
deformation had no influence on the value of Young's modulus if it was obtained by the load-longitudinal strain relation.Based on these previous studies, we assume a transversely isotropic material behavior for spruce and we propose to measure the displacement field by DIC on a "short" beam in order to measure accurately the shear effect.We determine both E and G out of the loaddisplacement relations by means of least square method on the bottom part of the specimen.Thus, the localized effect in the region near the loading point has no influence and a single test leads to both E and G values.The same test also gives the limit of proportionality R el .
Putting together DIC, 3 point bending test analysis, Monte Carlo simulation and anisotropic beam FEA, we provide a complete example to prove the feasibility of managing virtual tests for wood furniture.To address these points, the paper presents anisotropic behavior characterization in a first section together with the probability density function construction for all variable parameters.In the second section, we present the stochastic simulation of the bunk bed and statistical treatment of the results that are finally validated by comparing with the measurements done on real tested bunk bed.
Using all results presented, one can first manage series of 3 point bending test on specimen cut out from wood sample that will be used to manufacture the studied furniture.From DIC, one can obtain the experimental results that are used to identify the probability density functions for each uncertain parameter of the mechanical model of the furniture.Thanks to the finite element beam code developed one may run the Monte Carlo simulation and obtain the confidence region for deflections and strength security factor for the furniture.This approach is quite simple to develop in an industrial context.

Anisotropic behavior of wood (spruce)
In account of the orientation of the wood fibers and the manner in which a tree increases in diameter as it grows, wood is generally considered as cylindrically orthotropic material (Fig. 1) that has its unique and independent mechanical properties along three mutually perpendicular axes: longitudinal (l), radial (r) and tangential (t).
However, it has been noticed that the differences in mechanical properties between the radial and tangential axes are relatively minor when compared to the difference between the radial or tangential axis and the longitudinal axis [12].Consequently, to simplify we assumed that spruce could also be considered as a transversely isotropic material where the properties are symmetrical around the preferred axis of fiber.Generally, 5 independent constants are required to identify the elastic behavior of transversely isotropic material; this is illustrated by the following stress-strain relations of generalized Hooke's law (Eq.( 1)) where l (fibre direction) is the axis of symmetry: (1) where: In this study, we focus on Young's modulus E l , shear modulus G lt , the elastic resistance of spruce and their dispersions.We also study the influence of the fluctuations of these parameters on the displacement values at critical points.For this purpose a set of 3-point bending tests and tensile tests has been carried out to define these characteristics.

Characterization tests of spruce: 3-point bending test
Clear specimens of European spruce (Picea abies, approximate density of 450 kg.m −3 ) were sawn from different elements of the bunk bed (3 safety rails and 3 lathes).The specimens were small in size with a length span of 130 mm (parallel to grain), 15 mm in width and 30 mm in depth and with a moisture content of about 12%.The main objective of this experimental study is to obtain appropriate statistical measurements of two simultaneous parameters for each test (Young's modulus E l and shear modulus G lt ).Therefore, we resorted to the displacement fields of each wood specimen.These fields can be measured by using the technique of digital image correlation (DIC) which is widely used recently and has shown efficacy in mechanical testing of wood and wood products [32].The 3-point bending tests, accompanied by DIC, were executed until failure occurred on 27 specimens using the universal testing machine Deltalab500 (Fig. 2).
The results containing the rough displacements fields went through a smoothening process to avoid the influence of experimental noise on the results and compared with the analytical solution of the displacement field coming from the beam theory.One can see in Figure 3 that an important effect is localized near the application of the load, so it is possible to manage accurate identification using only the lower part of the beam (i.e.y < 0).The analytical solution (see Appendix for details) is given in Equation (2) for 0 < x < L/2 where the origin is located at beam center: F is the applied load, I z the quadratic moment of the cross section of the beam; U 0 , V 0 and Ψ are the rigid body displacement components of the beam.The objective is to adjust the parameters (E l , G lt ) and the rigid body displacement from the displacement field measured by DIC using the least square method.For each of the 27 tests, an initial picture is taken and four pictures (F = 200 N, 400 N, 600 N and 800 N) are analyzed to get the displacement components U i and V i at every node (x i , y i ).The unknowns of the problem are obtained by minimizing the sum of squared residuals e 2 This is a direct application of the least square method.
Developing the 5 derivations leads to a 5× 5 linear system that writes A.  where N is the number of nodes and: Each test gives a mean longitudinal modulus and a transversal shear modulus.A typical comparison between the measured and the analytical displacements fields is shown in Figure 4.The final results yielded average on all tests of E l = 9170 MPa with standard deviation (sd El ) of 1570 MPa (17% dispersion), and average of G lt = 550 MPa with standard deviation (sd Glt ) of 95 MPa (17% of dispersion).The mean values of E l and G lt corresponded to the literature; for example in reference [33].The tests yielded a mean value of 24.4 MPa for the bending strength at the proportional limit R el with a fluctuation of 13%.The cor-relation between the parameters derived from the static bending tests is shown in Figure 5.The coefficient of correlation amounts to r ER = 0.62 between elastic modulus E l and bending strength R el .Lower correlation was reported between the shear modulus G lt and both E l and R el (r EG = 0.26, r GR = 0.21).Consequently, we continued the study assuming that R el , E l and G lt are mutually dependent variables.These correlations and dispersions are not surprising; both stiffness and strength properties of wood are strongly and positively correlated to wood density, which in turn varies from one tree to another and even within a tree where it varies with height and along the radial profile [34].

Stochastic modeling of spruce's behavior
In order to use the probabilistic approach, the dispersion of each parameter, considered as random variables, (b) (a) should be represented by their probability density function (pdf for short).For a bunk bed problem, the uncertain parameters were limited to the wood elastic characteristics i.e.Young's and shear modulus (E l and G lt ) and the bending strength R el .

Statistical dependence of the random variables
The parameters (elastic moduli and strength property) were identified from the previous experiments and considering the correlation reported, they must be considered as mutually dependent.Indeed the statistical correlation between these variables should require more data to be proved.The literature points out an effective correlation between them, see for example [35].The number of experiments in this kind of industrial approach will always remain too small to give rise to possible correlations (typically from 20 to 40 when more than 400 are needed).This small sample size is not sufficient to use the Maximum Entropy principle including the covariant matrix as an information, but it will help here to turn the 3 chosen parameters in a space where the variables will be independent.If − → X is the vector of the natural variables: and [C x ] the covariance matrix associated: The pdf for each independent parameter Y i can provide numerous realizations of the random variable.Turning back in the natural variables space with the inverse rotation one can obtain the realizations of the coupled natural variables.For each realization of (E l , G lt , R el ) a FEA, presented later, was executed.The construction of the pdf density function is briefly summarized below.

Probability density functions for the random parameters
The 27 experimental values of triplets (E li , G lti , R ei ) give 27 triplets of the Y i variables.This number of experimental values is not enough to consider that the standard deviations are helpful information to use with the principle of maximum entropy, a widely used technique, to construct the pdf [20].The available information necessary to build the pdf p Y i (y i ) for the variables Y i are limited to Equation (9).
The first assures, through the normalization conditions, that all characteristics are positive, the second that the pdf gives the experimental mean value and the last condition of Equation ( 9) assures that the inverses of moduli 1/E l and 1/G lt exist and are strictly positive.1/E l and 1/G lt must be 2nd order random variables so that the strain value remains finite.c i is a finite positive value but we don't have any experimental information on it.It can be adjusted a posteriori using the evaluation of the standard deviation of Y i .Using this evaluation of σ Y i and conditions of Equation ( 9) leads to a pdf Gamma that writes: where: is the Gamma function, m Y i and δ i standing for respectively the mean value and the coefficient of variation of the random variable Y i .Without any further information, δ i can be considered as model uncertainty.If the standard deviation σ Y i is evaluated from the experimental values, the δ i value is equal to the ratio σ Y i /m Y i .We managed 5000 realizations of the 3 independent variables Y i following Gamma law's pdf and coming back to the natural random parameters E l , G lt and R el one can see in Figure 6 the "cloud" (E, G) is slightly correlated; At the opposite, the "cloud" (E, R e ) highlights a strong correlation that has an important effect on the numerical simulations presented in the next part.

Numerical simulation of structural bunk bed behavior
A bunk bed is usually manufactured from long wood planks where one dimension is large compared to the two others.As we consider it as a structure designed to take loads, we propose the beam theory to model its structural elements.The choice of the beam theory depends on several factors, one of which is the anisotropy degree of the material.Generally wood has a much higher Young's modulus (up to 20 times) than its shear modulus (see the experimental section).Under these circumstances, the deformation of the beam due to shear, which is usually small and neglected in Euler-Bernoulli's theory, increases and should be taken into account.Therefore, two nodded Timoshenko beam element was chosen to model the bed structure.For this particular case, a Matlab code for finite element analysis was developed to simulate the normalized tests on a bunk bed with the principle aim of including the uncertainty of the mechanical parameters and its effect on the response of the structure as shown later.
Finally, a test was performed on a real bunk bed to confirm the methodology adopted in this work.

Mean problem and results
Based on the European standard [36], the mean problem we focused on is the normalized test in which we applied a vertical force of 1000 N on the middle of the long top safety barrier (node 29).
Considering a mass density of 450 kg.m −3 for spruce and the mean test results for E l and G lt , the FEA results showed a maximum displacement of 2.27 mm at the loading point as was expected (Fig. 7) and the associated maximum stress in this section was σ xx = 4.7 MPa and σ xy = 0.3 MPa.Whereas the maximum stress on the level of the structure was found at the node (13) which is the intersection point of the bed column and the long safety barrier and it was σ xx = 5.1 MPa, σ xy = 0.4 MPa.The contribution of shear stress in Hill criteria is only 15% of the tension stress contribution.Hence, we chose to ignore this contribution and evaluate the strength of the structure using strength criteria written for a beam problem as: where K s is the safety factor amounts to 4.8 and 5.2 for nodes 13 and 29 respectively.For the mean problem no node exceeds 1 and thus the structure is validated for this 503-page 7 test.Still, the strengths and the mechanical properties used have shown dispersions up to 17% according to the tests.Therefore, a Monte Carlo simulation, shown below, was performed to take into account the uncertainty of these parameters.

Monte Carlo simulation of bunk bed structural behavior
In this section we will examine the result of Monte-Carlo simulation [37] of the bunk bed structural behavior problem in the former normalized test.In relation to 2 uncertain variables (E l and G lt ) and their associated pdf, we managed 5000 finite element simulations.We investigated two results, the maximum displacement (here, the node where the load is applied) and the maximum stress.This procedure will lead to estimate a confidence region of the solution for both chosen results with respect to the measured dispersion of the uncertain parameters of the mechanical problem.We now present the results on maximum displacement with this approach.

Monte-Carlo numerical simulation: convergence
Convergence of the Monte-Carlo numerical simulation can be shown by calculating the mean value and the standard deviation of the maximum displacement.In Figure 8, both terms are plotted versus the number of realizations managed.We can see that convergence was managed for 200 realizations.If this was known, we could have reduced the number of simulations and consequently the CPU time used.
The mean value of the displacement appears to be equal to 2.36 mm which is slightly different from the maximum displacement value obtained by solving the mean problem.Let us recall that the maximum displacement was equal to 2.27 mm in the mean problem which is 4% less than the mean value of the probabilistic approach.This small difference is due to the asymmetry of the pdf of the unknown parameters around their mean value.

Mean displacement and confidence range
With a relatively small number of realizations, convergence can be obtained for the mean value and the standard deviation of U max .However, this number is not sufficient to manage the convergence for higher moments.Based on the 5000 realizations we have, we can perform an estimation of the probability density function of the random variable.As illustrated in Figure 9 a bar chart was constructed using a number of classes (50 classes of probabilities for example) along the range [U m Min , U m Max ], where U m Min and U m Max are minimum and maximum values of U max over the 5000 realizations.Counting the number of realizations that give a U max value for each class and reporting this number to the 5000 values permits to determine the probability of having a realization that belongs to this class.When the number of realizations increases, the chart bar takes a continuous form (superposed line on the bar chart of Fig. 9).
The constructed chart bar in Figure 9 makes it possible to define a confidence range from whose probability of finding a realization is higher than P c (P c = 95% for example).Equation (12) defines its upper and lower limits.
where ζ(p) is a fractile given for a random variable X with a distribution function F X (x), by: One can see that the distribution is not symmetrical around the mean value and certainly higher moments should be taken into account to build a more accurate pdf for Umax which can be shown in the pdf (full line) that takes into consideration additional information of the third and fourth order moments.
Apparently, this can be executed by counting the realizations ordered by displacement values until the density function reaches P c for the upper boundary and (1 − P c ) for the lower boundary.The boundaries are U + m = 3.38 mm and U m = 1.7 mm.The mean value 2.36 mm belongs to the confidence region but it is not located in the middle of it (the median value should be 2.54 mm).

Stress and bed strength probability
According to the Monte Carlo simulation, stress presented low dispersions whereas bending strength R el yielded a dispersion of 13% during the experimental tests.If we focus on one critical section, the safety factor K s can be statistically investigated using the Monte Carlo method.The confidence interval of K s gives [4-6.7] for node 29 and [3.5-6.2] for node 13.Apparently the safety factor for both critical nodes reported is higher than 1 for all realizations (Fig. 10) which indicates that no risk is reported at the critical nodes and the bed strength design is validated for the normalized test.

Bunk bed test results
To confirm this work methodology and the different assumptions made until now, we managed a validation test on a real bed given by the industry (only one bed was available) shown in Figure 11.The force was applied using a rigid wooden loading pad of 100 mm in diameter.
The loading machine was capable of applying an accepted load of 1003 N and no damage was noticed at the contact zone.The resulting amplitude was captured using a displacement comparator with a tolerance of ±0.05 mm.The force was applied ten times.It caused almost the same displacement of 2.6 mm each time.This value was included in the confidence region.

Conclusions
In this study, several tools were developed and different techniques were used to achieve a virtual normalized test on a wooden bunk bed.In the experimental section, the DIC process yielded reasonable mean values for modulus E l and G lt .The results highlighted important percentages of dispersion in the mechanical properties of spruce.This dispersion was taken into account by constructing probability density functions for each uncertain parameter using the maximum entropy principle.
The same bending test also enables to measure the strength (proportionality limit R el ) and correlations between these mechanical properties are taken into account in our pdf construction.Only 27 tests were performed and more accurate pdf could have been constructed if more samples had been tested.
Applying FEM and the Monte Carlo simulation using beam elements permitted to identify the confidence range for the maximum displacement and the strength safety factor for the suspected nodes.These results were confirmed by applying the test on a full size prototype.
This accurate and simple to use procedure is well adapted to the industrial context and is a great help in developing better furniture designs that could pass the standardized tests in a short cycle of product development.Both experimental evaluation of the mechanical properties and the numerical simulation of the furniture behavior can be provided in a design department and used to optimize the structure before managing a real size validation test.where c 1 is chosen to satisfy the free edge conditions on y = h/2 and -h/2 For elastic beam one can obtain the strain from: This leads to the two differential equations: From the first equation, one can obtain: Including this result in the second, one can obtain: In beam theory ε yy is naught.This last condition leads to: 12G lt Iz + yψ + U 0 g (y) = 0 ⇒ g(y) = V 0 U 0 , V 0 and ψ are constant values that can be interpreted as the rigid body motion components.Replacing all in U and V expression leads to the analytical expression of Equation (2):

Fig. 2 .
Fig. 2. Digital image correlation setup (speckle pattern applied on the surface of the sample).

Fig. 4 .
Fig. 4. Comparison between the displacement components' contour for a typical test (bottom contours) and beam model components (top contours).

Fig. 5 .
Fig. 5. Bending test results.Linear regression correlation between elastic modulus E l and bending strength R el (a), E l and shear modulus G lt (b).

r
EG = 0.26 r GR = 0.21 r ER = 0.62(7) Searching the eigenvalues C i and associated eigenvectors − → ϕ i of [C x ] leads to the independent variables via:

Fig. 6 .
Fig. 6. 5000 realizations of the correlated random variables E l and bending strength R el (a), E l and shear modulus G lt (b).

Fig. 8 .
Fig. 8. Evolution of mean value (dashed line) and standard deviation (solid line) of the maximum displacement versus the number of realizations.

Fig. 9 .
Fig. 9. Probability density function of the maximum displacement on the bed.(Bar chart) is the representation from the Monte-Carlo simulation with 50 classes.(Dot line) is the pdf identified from the mean value and standard deviation calculated from Monte-Carlo simulation.One can see that the distribution is not symmetrical around the mean value and certainly higher moments should be taken into account to build a more accurate pdf for Umax which can be shown in the pdf (full line) that takes into consideration additional information of the third and fourth order moments.