Minimizing flooring strip weight: A shape optimization approach

Deteix, J., Djoumna, G., Blanchet, P., Fortin, A., and Cloutier, A. (2012). "Minimizing flooring strip weight: A shape optimization approach," BioRes. 7(2), 1931-1947.

Abstract

In North America, flooring strips are manufactured with grooves at the back. There are various reasons for these grooves but, historically, they were considered a strategy to reduce weight and transportation costs as well as improving dimensional stability. As no data are available to assess best practices in terms of performance, we have investigated methods to reduce flooring strip weight. One way to achieve this is to adjust the number and shape of grooves. Using warp as a comparison tool, we were able to analyze the merits of a finite number of designs. With this approach, however, we could not guarantee that the result was the most favourable. The search for a solution led to design optimization, i.e.: minimizing weight by acting upon a part of the strip’s shape, taking into account its warp resistance or stiffness. This paper describes an optimization strategy adapted to the calculation of the optimal design subjected to arbitrary mechanical and geometrical conditions (including the thickness of the wear layer). This approach is not limited to flooring strips, and it can be used in any situation where a linear hygromechanical model is relevant. This strategy involves two steps: global optimization with respect to admissible variations of the shape (or design) followed by a post-processing phase that takes into account various other mechanical and possibly geometrical conditions imposed on the strip.

Download PDF

Full Article

MINIMIZING FLOORING STRIP WEIGHT: A SHAPE OPTIMIZATION APPROACH

Jean Deteix,^a George Djoumna,^b Pierre Blanchet,^a,c,* André Fortin,^a and Alain Cloutier ^a

Keywords: Wood flooring; Design optimization; Finite element method; Weight, stiffness; Moisture content; Linear elasticity

Contact information: a: Université Laval, Québec, QC, Canada; b: University of Waterloo, Waterloo, Ontario, Canada; c: FPInnovations, 319 Franquet, Quebec City, Canada, G1P 4R4.

* Corresponding author: pierre.blanchet@fpinnovations.ca

INTRODUCTION

The hardwood flooring industry has a long manufacturing tradition, and back side grooves are part of this tradition. According to the National Wood Flooring Association (NWFA 2009), back side grooves have been used to: a) show which side is up in rough flooring strips; b) create air movement and exchange from the back side; c) help the adhesive grab the strip in glued-down installations; and d) insure that the strips are lying flatter. It is also said that, historically, back side grooves have been used to reduce strip weight, hence transportation costs. Today, many flooring manufacturers are using these fibers from the grooves for self-supply of energy or for pellets production. In many cases, the pellet production has become an important factor of the economic viability of the manufacturers as well as improving sustainable qualities of the flooring manufacturing. As there is no standard for these grooves, a wide variety of patterns can be found on the market. In fact, there may be as many grooving patterns as there are manufacturers.

Our purpose in this work was to reduce flooring strip weight and to optimize waste fiber production for energy while minimizing distortion. Lower weight leads to lower transportation costs and a better use of the resource. Analyzing the various designs of parts made of materials controlled by hygrometric (and possibly thermal) conditions can be a difficult and time-consuming task. The development of a simulation tool allowing for weight reduction at the research stage would therefore be useful.

Variations in hygrometric conditions can induce undesirable hygromechanical deformations. For appearance products such as parquet flooring, even slight deformation is unacceptable. Adsorption and desorption of water vapour may induce cupping, thereby decreasing product value. Consequently, reducing warp caused by variations in the thermo-hygrometric conditions is of critical importance. The use of warp as a comparison tool could serve to analyze the merits of a ﬁnite (discrete) number of designs relative to some hygrometric variations. Given that the number of designs is finite, a hygromechan-ical analysis of each design is possible. However, this approach gives biased results since we deal with a pre-established set of designs. There is no guarantee that the result is the most favourable. The search for discrete design optimization must avoid subjectivity. The development of a systematic and unbiased way of lowering weight while maintaining acceptable dimensional stability provides a strong incentive to optimize strip design.

This project was a first approach to the development of an analysis and optimiza-tion strategy adapted to the calculation of the optimal design of parts subjected to arbitrary static loading and unsteady thermal and hygrometric conditions. The hygromechanical analysis followed the works of Ganev (2005), Cloutier (2006), Blanchet (2006), Deteix (2008), and Belleville (2008). The numerical model relied on a three-dimensional finite element approximation of the solution of the hygromechanical model. For the shape optimization part, we first devised a strategy (loosely inspired by Mins 1999) allowing us to place this structural optimization problem with transient loading in the context of classical shape optimization. Shape optimization is now a mature theory (Duysinx 1996; Bendsøe and Sigmund 2002; Delfour and Zolésio 2002) that is used in various disciplines. For the optimization process, the numerical approach was based on the SIMP method (Solid Isotropic Microstructure with Penalization), as presented in Bendsøe (1989) with modification for the anisotropic nature of wood, and using the heuristic approach of Sigmund (1997) for the filtering of the checkerboard effect (numerical artefact produced by the choice of some degrees of interpolation of the displacements) and making the design independent of the finite element grid. The proposed strategy was implemented in an object-oriented language (C++), taking full advantage of a modern and powerful software library called MEF++, developed by the GIREF (Groupe Interdisciplinaire de Recherche en Éléments Finis) research centre over recent years.

Hygromechanical Model

The aim was to describe the mechanical behaviour of a wooden part, identified as , subjected to mass transfer (moisture movement). We assumed isothermal conditions and a variation of hygrometric conditions. Mass transfer occurs by free convection from the surfaces. Transient moisture movement is described by the following three-dimensional moisture conservation equation,

where K_M is the tensor of effective water conductivity (kg m^-1 s^-1 %^-1); D is the tensor of effective moisture diffusion (m2 s^-1); d_b is the basal density (kg m^-3); and M is the moisture content (%).

Wood is orthotropic and elastic. No mechano-sorptive effects were taken into account, as previous work demonstrated elastic behaviour under the conditions considered in this study (Blanchet et al. 2003). Since we assumed the material to be conditioned, the flooring strip was initially free of stress. The governing equation describing the mechanical aspect of the problem was a three-dimensional (incremental quasi-static state) equation of equilibrium,

where body forces were assumed to be negligible. The quantities σ_ij are the normal and shear stress components, expressed in a rectangular coordinate system. Because the wood material was considered elastic, Hooke’s law could be used to relate stresses to strains,

where C_ijkl is the stiffness tensor; ε_kl is the strain tensor; β_kl is the moisture shrinkage/swelling coefficients (%^-1); and ∆M is the moisture content change between two time steps (%). The stiffness tensor was constant (in time, no ageing), and we took hysteresis in adsorption/desorption into consideration for the shrinking/swelling coeffi-cients (Goulet-Fortin 1975). The strains were related to the displacements, u₁, u₂, and u₃, measured along the x₁, x₂, and x₃ directions, respectively.

Shape Optimization under Transient Conditions

The mechanical behaviour of a flooring strip is transient. The search for an optimal design (we refer indiscriminately to “shape”, “design” or “geometry” to describe the complete definition of the geometry of the part under study) should take this transitory nature into account. This implies the development of an optimization process to establish the most favourable design (shape) in relation to weight reduction and suitable stiffness over a fixed period of time [0, t_max].

As this approach was not feasible from a practical point of view, we adopted a suboptimal approach, inspired by (Min et al. 1999), which is mathematically sound and easy to implement. This approach consisted in determining the optimal shape at a fixed time (critical time), and then verifying a set of mechanical criteria on the design for the whole period. When the criteria were satisfied, we considered our design to be optimal. In case of failure, we repeated the process using the current shape and the time at which failure of the criteria had occurred as the new critical time (Fig. 1)

Fig. 1. Optimization process for transient behaviour or cost function

Obviously, such an approach has some drawbacks. The shape of the part generally affects mass transfer, which is needed to calculate stress. Any shape variation will there-fore lead to a new calculation of the moisture content up to the critical time. Considering that the stiffness tensor or shrinking/swelling coefficients vary with time, it is all the more necessary to consider displacement history (moisture content in particular). Moreover, stringent mechanical criteria can lead to numerous and costly resolution of the static shape optimization process.

By limiting our study to a stiffness tensor that is constant in time, we avoided some of the difficulties. However, due to the effect of hysteresis on shrinkage coeffi-cients, we had to take moisture content history into account. Even with hysteresis, other assumptions can make this strategy effective.

The simplest case occurred when the moisture content was indifferent to modifications to the domain resulting either from the nature of the boundary conditions added to or from the geometry being restricted to a region with little or no mass transfer. The problem was much simpler in this case, as the moisture content had to be calculated at most once for each value of the critical time. With flooring strips, the boundary condition on moisture content meant that these conditions applied. Even though this strategy could be relatively slow (depending on the mechanical criteria imposed) it yielded acceptable results.

Shape Optimization of the Unsteady Strain Energy

It will be remembered that we used to identify the domain of occupied by the part. This domain is included in the maximum domain and consists of two components and

where is a fixed component (possibly empty) and is the zone where variation of the geometry is permitted and contained in a larger and fixed research region (Fig. 2).

Fig. 2. Two-dimensional example of geometrical decomposition. On the left, the maximum domain (known as the hold-all). On the right, a possible design included in the hold-all with a variation of the geometry contained in the modifiable zone.

For given and values, we defined the set of admissible shapes as the set of domains satisfying .

(6)

Since we were looking for an admissible optimal shape, the optimization process was restricted to .

The weight of the part was considered at moisture content and could be expressed as a function of the basal density of the strip (FPL 1999):

If only cupping was considered, bad design could result, with potential increases in other types of warp such as bow. To avoid this problem, we used strain energy as a global indication of warp. When, as usual, we limited imposition of the displacements to a discrete number of points in the domain (Dirichlet boundary condition), since the moisture component of the strain related to an external force, we had:

where is the load applied on , a part of the boundary. If no load was applied ( ) and was constant in time on the interval , corresponding to a period of adsorption or desorption, then we had:

Two approaches were available to deal with the problem of reducing the weight of a part while taking into account its stiffness over a period of time . The first was to minimize weight within an acceptable range of displacements (or stiffnesses). The second consisted in minimizing displacements (through compliance with ) while imposing a restriction on weight. The latter model is classical in the sense that it is the approach usually described in the optimization literature (Delfour 2002; Bendsøe 2002; or Duysinx 1996). Of course, both models are acceptable. We chose the second, as it was easier to set a weight reduction target and try to achieve it while keeping the part as rigid as possible. We maximized stiffness for an imposed minimum weight reduction (at least % of the weight of ),

, such that

with M and u resulting from the solution of and on with appropriate initial and boundary conditions.

In order to solve this problem using the strategy illustrated in Fig. 1, we needed to refine the classical (i.e. steady objective function and state variable) shape optimization process.

Static Shape or Topological Optimization

Shape optimization in its most general setting should consist of a determination for every point in space whether there is material at that point or not. As the theory related to shape optimization evolved, a clear distinction has emerged between shape optimization (applying to perturbations at the boundaries of an established shape) and topological optimization (determination of the topology (holes, disjoint components) of the shape). These techniques are generally used at different stages of optimization and performed separately. Topological optimization is usually considered as a tool for finding efficient design concepts at an early design stage, whereas shape optimization is viewed as a tool for detailed design at a later stage. In either case, the theory is based on a finite element (FE) mesh used to discretize the design domain.

The boundary variations approach consists in studying perturbations of the boundary for a given topology (lay-out) of the structure. The literature on the subject is quite extensive (see Delfour and Zolésio 2001 for a review). This method can be implemented in a number of ways, for example by employing certain mesh moving schemes to define the shape of a given structure. Implementation of the boundary variations approach is not straightforward, and it normally requires repeated application of a re-meshing method during an iterative optimization scheme. This approach also has limitations, as it only allows for the prediction of the optimum shape of the boundaries of a given initial topology.

In topological optimization, the topology of the structure is not set a priori (it is possible to generate/eliminate holes and disjoint components). For finite element discreti-zation, every element is a potential void or structural member. Topological problems formulated this way are inherently discrete optimization problems, but there are various ways of solving them without resorting to discrete optimization algorithms. One way is to use continuous approximations based on heuristics. The most satisfactory approximations are obtained by introducing composites such as layered structures or porous, periodic media (Bendsøe and Kikuchi 1988). This homogenization approach leads to parts made of physically acceptable material; however, the determination and evaluation of optimal microstructures and their orientations are cumbersome. Furthermore the resulting structures are of no practical interest if microstructures are to be avoided.

The SIMP method (Solid Isotropic Material with Penalization) (Bendsøe 1989) provides an alternative approach to homogenization. Here, material properties are assumed to be constant within each element of the FE mesh, and the variables are the element volume fraction of solid. The material properties are modeled as,

where is a material property, the solid material property, the volume fraction of solid, and the power of the penalizing. This approach eliminates the use of complex rules of mixture for intermediate values of volume fractions. For reasonable values of p, it actually yielded solutions with very few, if any, intermediate values.

As a first approach in flooring strip shape optimization, we favoured topological optimization, which allowed us to produce more general designs and avoid the multiple re-meshing associated with boundary variations. We selected the SIMP method because this approach to topology optimization can be applied to problems with multiple constraints, multiple physics and multiple materials, and was therefore more general. It also allowed us to circumvent the difficulties related to intermediate values of the volume fraction.

The two major numerical difficulties with this approach were the checkerboard effect (regions with alternating solid and void elements) and mesh dependency (qualitatively different design for refined FE mesh). A review of these problems and their treatment can be found in Diaz and Sigmund (1995), Duysinx (1996), Sigmund and Petersson (1998), Bendsøe and Sigmund (2002), and Bruns (2005). In this study, we applied a mesh-independent algorithm that both eliminated the checkerboard problem and the mesh-dependency problem (Sigmund 1997). The method works by replacing the element-wise derivatives with respect to the volume fraction (sensitivities) with a weighted average of the sensitivities of their neighbours within a given radius.

Topological Optimization of a Flooring Strip

We considered a flooring strip made of sugar maple as shown in Fig. 3. Sugar maple was selected due to the fact that all the physical and mechanical properties were available for modeling and as it is an important commercial species in floor covering in North America. Physical and mechanical properties are the basis of this model; species with similar properties will have similar behavior according to this model. This basic flooring strip, geometrically composed of three rectangular prisms, was our hold-all (we were looking for shapes contained in this basic shape). The top surface was 107.9 mm wide by 609 mm long. The lay-up consisted of a 7.2 mm-thick surface layer, a 6.4 mm-thick core layer and a 5.1 mm-thick backing layer, for a total thickness of 18.7 mm and with 5.5 mm-deep tongue and groove. All three layers were oriented in the same direction. For modeling purposes, the directions X, Y, and Z correspond to the longitu-dinal, tangential, and radial direction of wood, respectively.

Fig. 3. The basic flooring strip used as hold-all for the model

The temperature was kept constant at 20ºC, and a variation of hygrometric conditions caused water vapour desorption from 8.6 to 6.3 percent moisture content at the top surface only. This corresponds roughly to a relative humidity (RH) decrease from 80 to 50 percent at 20ºC. The initial condition and boundary condition were as follows:

The moisture transfer boundary conditions specified no moisture flux through any surface other than the top surface, through which desorption was assumed to occur. For that surface, the moisture flux was given as follows,

where is the convective mass transfer coefficient (kg m⁻² s⁻¹ %⁻¹) and is the equilibrium moisture content (%).

The set of admissible shapes was defined on the premise that the top (apparent) surface had to be kept intact, as did the tongue and the groove portions. To satisfy these restrictions, we defined the research region as a region included in the last two layers (Fig. 4). And was the complement set of with respect to the original strip (Fig. 3).

Fig. 4. The flooring strip viewed upside down, view of (without any material in this case)

With this definition of the set , the moisture content solution to (1)-(12)-(13) was indifferent to variations of the domain in . It performed as a restriction to the solution on . Our shape optimization problem therefore fell into the category of cases where the above-described strategy was most favourable: a single resolution of — was adequate to determine the optimal solution.

In terms of the mechanical component, the strip was initially at rest with no stress

In order to have a freestanding strip, we imposed only sufficient conditions to avoid rigid motions and rotations:

For our purpose, we had to set a time interval, and chose to monitor strip behaviour over 42 days. Our goal was to determine the optimal shape of a flooring strip exposed to an ambient relative humidity of 80% for a period of 42 consecutive days. To implement our strategy (Fig. 1), we needed a value for the critical time. From a finite element analysis (Fig. 5) for a freestanding strip, we concluded that the first critical time should be 15 days, which approximately coincided with the time when extreme displacement values were observed.

Fig. 5. Typical displacements of a free standing flooring strip measured a three points: at one corner (0,0,0.0187), at one edge(0,0.05395,0.0187) and on the side (0.3045,0,0.0187)

In line with our strategy, we had to solve a static version of the optimization problem at the critical time,

where M is the solution to — over the 42 days on the hold-all . Obviously is based on the fact that the stiffness tensor is constant over time, and the boundary conditions imposed on the flooring strip. For the same reason the displacements u solution of and – is the solution of,

with through . In , given that basal density was constant and we chose to measure weight at the initial moisture content which was constant throughout the hold-all, the weight constraint was equivalent to a volume constraint,

where V identifies the volume. Hence, the following topological optimization problem:

, such that

where u is the solution of with through .

Numerical Approach

The numerical approximation was based on finite element modeling, and was performed using the finite element code MEF++. Finite element discretization of the Galerkin weak form of mechanical equilibrium (19) and moisture conservation (1) was performed using standard isoparametric and linear interpolation of the unknown displacements u and moisture content M.

Time discretization of the equation was performed by the standard Euler implicit time-marching scheme. The predicted values of M, and u depend on position and time. A single coupled system of discrete equations was solved for the displacement components and M at each time step. A user-specified initial time increment of 0.5 s was used. The following time increments were automatically adjusted between 0.1 and 100,000 s by the MEF++ software based on convergence rate.

On the basis of the finite element approximation of the displacement, and using the SIMP method, we produced a discrete optimization problem that could be solved with a standard optimizing scheme.

Discrete Topological Optimization

was considered the finite element mesh composed of Q elements and N nodes corresponding to the hold-all , and identified the element of the mesh. Following -(7) we imposed:

Since we were using the SIMP method to achieve topology optimization, we actually optimized material distribution by controlling material density in each finite element. Hence, the stiffness tensor was dependent on the relative density, , of the solid in element

with p the power applied to the volume fraction of the solid in element e acting as a penalty factor. The value of p was selected in such a way as to avoid designs containing intermediate values of the fraction volume. In this case p = 3 produced designs presenting almost no intermediate values.

With being used to identify the Q-vector containing the design variables, the finite element problem became,

where the global stiffness matrix K is a (3N 3N) matrix (three displacement degrees of freedom (d.o.f.) per node), U is the 3N displacement vector, F is the structural load vector (3N). The system matrix and vectors are dependent on the element relative density vector . The system matrix K and global load were assembled in the usual way, as sums over elements,

where is the element stiffness matrix and the element load vector, both functions of the material properties in the element. With B identifying the usual strain-displacement matrix, and using and , the element system matrix and load vector could be written as:

Given that the volume constraint only applied to the research region , and introducing V_e as the volume of the element e, we had:

We were then in a position to formulate the discrete optimization problem:

such that

with U: the solution of .

This discrete optimization problem can be solved by various methods. For simplicity’s sake, we used a standard optimality criteria (OC) method (see Bendsøe and Sigmund 2001). This method required the derivative of with respect to . Using we had:

We omitted the derivatives of , K, and F, which were also needed, as they were relatively easy to compute.

Checkerboard Effect and Mesh Dependency

The checkerboard effect is merely a discretization error of the FE method, which has nothing to do with SIMP. This error results in an overestimation of the stiffness. If computer constraints on memory or processing time are not a consideration, checkerboard can be controlled using higher order elements (Sigmund and Petersson 1998). However, as we were dealing with very fine meshes, we did experience time and memory problems, which dictated the use of lower order elements (linear in this case), with a potential risk of checkerboard effect. Moreover, mesh-dependency (qualitatively different design for refined FE mesh) could not be ignored.

In order to avoid the problems associated with checkerboard effect and mesh dependency in the solutions to through , some sort of restriction on the resulting design had to be introduced (see Sigmund and Petersson 1998, for an overview). We used a ﬁltering technique (Sigmund 1997; Bruns 2005). It must be emphasized that this ﬁlter has not yet been proven to ensure existence of solutions, but numerous applications in different fields of research have established that it can yield mesh-independent designs in practice.

The basic idea of the filter is to use blurring techniques borrowed from image processing. Most image processing filters modify the colour (density) of each pixel (element), based on the colours (densities) of its immediate neighbours. To ensure mesh-independency in the present problem, the filter had a fixed geometrical size, meaning that the modification of an element was dependent on the elements lying within a fixed radius of the considered element. The ﬁlter works by modifying element derivatives (or sensitivities) as follows:

The convolution operator (weight factor) is expressed as:

where the operator dist(e, f ) is deﬁned as the distance between the centre of element e and element f, and r_min is a fixed mesh-independent radius. The convolution operator is zero outside the ﬁlter area (determined by the radius r_min) and decays linearly with the distance from element f.

Instead of the original sensitivities , the modified derivatives are used in the optimal criteria (OC) method. Obviously filtered sensitivity converges towards original derivatives as r_min, approaches zero. Basically, the filter eliminates structural variation under a certain length scale. Typically an optimization problem solution is initiated with r_min, equal to 10% of the smaller dimension of the design domain.

Numerical Results

We constructed a finite element mesh sufficiently refined for the moisture content

and displacements approximations to be independent of the mesh. As per , this mesh was composed of two zones. The fixed zone, , consisted of 12,264 elements (11,604 8-node linear elements and 620 6-node linear elements), the prismatic elements being used to allow for a finer mesh for the research region. was restricted to the last layer of the flooring strip (0.594 m x 0.0719 m x 0.0051 m). This region was located 0.021 m away from the border on the tongue side and 0.015 m on the groove side. It consisted of 5,940 8-node linear elements with an average volume four times smaller than the volume of the elements of . The flooring strip thus consisted of 18,504 elements and contained 19,890 nodes.

Obviously, the size assigned to the research zone limited the maximum reduction of the flooring strip, which stood at 17%, since the research zone represented 17% of the total volume. A volume decrease of (100 ) % of the research zone corresponded to a (17 ) % decrease of flooring strip volume. For the physical and mechanical properties of the material composing the flooring strip, we used the following values (Table 1).

Table 1. Material Properties of the Model

¹Belleville (2008),²Goulet et Fortin (1975), ³Bodig et Jayne (1973),⁴Siau (1995). L : longitudinal ; R : radial ; T : tangential

We chose two values for weight reduction, =0.5 and = 0.75, corresponding to reductions of 8.5% and 12.75%, respectively, of total strip weight. In both cases, it is worth underlining the effect of such reductions. Based on the given flooring strip, we computed the approximate reduction in weight for a 1 square metre surface area of flooring strip (see Table 2, second column). The third column indicates weight reductions for a fixed volume of flooring strips (1 cubic metre). The last two columns relate to transportation: For a load of 12,500 kg, we have an increase of the volume of flooring strip; and, in the last column, we have the weight reductions achieved for a fixed volume of strips (16 m³ representing approximately 12,500 kg of maple strips).

Table 2. Impact of Weight Reduction for the Two Cases Treated

These approximations clearly illustrate the impact on transportation. Of equal importance is the impact on resource optimization, since we have recovered 66 and 100 kg respectively of material for every cubic metre of flooring strips, which can be used in different applications such as wood pellet for energy production.

CONCLUSIONS

Even though the application chosen to illustrate this approach was relatively simple and academic in its considerations, it contains in its formulation and implementa-tion the most serious problems that one faces when dealing with shape optimization for a material under hygromechanical conditions. The approach can be applied to a problem involving more complex conditions in terms of geometry, make-up (engineered wood flooring (EWF), composites, etc), physics (elaborate boundary conditions, multiple mechanical constraint), or cost function.

As a first approach in shape-topological optimization, the authors consider the results promising. Clearly, the method has to be refined for better integration into practical and industrial applications. However, the authors feel that it should not be seen as a mere theoretical exercise. As environmental protection and better resource management become paramount in our society, the use of tools such as the proposed approach appears inevitable.

REFERENCES CITED

Belleville, B., Blanchet, P., Cloutier, A., and Deteix, J. (2008). “Wood-adhesive interface characterization and modeling in engineered wood flooring,” Wood Fiber Sci. 40(4), 484-494.

Bendsøe, M. P. (1989). “Optimal shape design as a material distribution problem,” Struct. Optim. 1, 193-202.

Bendsøe, M. P., and Kikuchi, N. (1988). “Generating optimal topologies in optimal design using a homogenization method,” Comp.Meth. Appl. Mech. Engrg. 71, 197-224

Bendsøe, M. P., and Sigmund, O. (1999). “Material interpolations in topology optimization,” Arch. Appl. Mech. 69, 635-654.

Bendsøe, M. P., and Sigmund, O. (2002). Topology Optimization: Theory, Methods and Applications, 2nd ed., Springer Verlag, Berlin, Heidelberg

Blanchet, P., Cloutier, A., Gendron, G., and Beauregard, R. (2006). “Engineered wood flooring design by the finite element method,” Forest Prod. J. 56(5), 59-65.

Blanchet, P., Beauregard, R., Cloutier, A., Gendron, G., and Lefebvre, M. (2003). “Evaluation of various engineered wood flooring constructions,” Forest Prod. J. 53(5), 30-37.

Bodig, J., and B. A. Jayne (1993). Mechanics of Wood and Wood Composites, Krieger Publishing Company, Malabar, FL, USA. 712 pp.

Bruns, T. E. (2005). “A reevaluation of the SIMP method with ﬁltering and an alternative formulation for solid–void topology optimization,” Struct. Multidisc. Optim. 30, 428-436.

Cloutier, A., and Deteix, J. (2006). “Modélisation du procédé de lamination des panneaux de fibres,” Research report presented to Natural Resources Canada, Value for Wood, Projet no.UL4, 79 pp.

Delfour, M. C., and Zolésio, J. P. (2001). “Shapes and geometries: Analysis, differential calculus, and optimization,” Society for Industrial and Applied Mathematics, Philadelphia, PA, USA. 2001.

Deteix, J., Blanchet, P., Fortin, A., and Cloutier, A. (2008). “Hygromechanical modeling of the adhesive line in laminated appearance wood composites,” Wood Fiber Sci. 40(1), 132-143.

Diaz, A. R., and Sigmund, O. (1995). “Checkerboard patterns in layout optimization,” Struct. Optim. 10, 40-45.

Duysinx, P. (1996). “Topology optimisation: From continuum media to elastic structures,” Ph. D. Thesis, Université de Liège, Raport LTAS OF-37, (372 pages). In French.

FPL (1999). Wood Handbook – Wood as an Engineering Material, Gen. Tech. Rep. FPL-GTR-113. U.S. Department of Agriculture, Forest Service, Forest Products Laboratory, 463 pp.

Ganev, S., Cloutier, A., Gendron, G., and Beauregard, R. (2005). “Finite element modeling of the hygroscopic warping of medium density fiberboard,” Wood Fiber Sci. 37(2), 337-354.

Goulet, M., and Fortin, Y. (1975). “Mesures du gonflement de l’érable à sucre au cours d’un cycle de sorption d’humidité à 21ºC,” Note de recherches nº12, Département d’exploitation et utilisation des bois, Université Laval, Québec, QC, Canada.

Min, S., Kikuchi, N, Park, Y.C., Kim, S., and Chang, S. (1999). “Optimal topology design of structures under dynamic loads,” Struct. Optim. 17, 208-218.

NWFA. (2009) http://www.nwfa.org/member/. Accessed on November 23^rd, 2010.

Siau, J. F. (1995). “Wood: Influence of moisture on physical properties,” Department of Wood Science and Forest Products, Virginia Polytechnic Institute and State University, Blacksburg, VA, USA, 227 pp.

Sigmund, O., and Petersson, J. (1998). “Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima,” Struct. Optim. 16, 68-75.

Sigmund, O. (1997). “On the design of compliant mechanisms using topology optimization,” Mech. Struct. Mach. 25(4), 495-526.

Article submitted: July 22, 2011; Peer review completed: January 8, 2012; Revised version received and accepted: March 2, 2012; Published: March 12, 2012.