New Perspectives for LVL Manufacturing from Wood of Heterogeneous Quality—Part 2: Modeling and Manufacturing of Variable Stiffness Beams

: This paper presents a new strategy in the use of wood of heterogeneous quality for composing LVL products. The idea is to consider veneers representative of the resource variability and retain local stiffness information to control panel manufacturing fully. The placement of veneers is also no longer random as in the ﬁrst part of this group of papers but optimized for the quality of veneers according to the requirement of bending stresses along the beam. In a four-point bending test arrangement, this means the high-quality veneer is concentrated in the center of the beam in the area between the loading points where the bending moments are the most important, and the low quality is located at the extremities. This initiates the creation of variable stiffness beams. This is driven by an algorithm developed and tested on representative veneer samples from the resource. Four LVL panels were manufactured by positioning the veneers in the same positions as in an analytical calculation model, which allowed the calculation of beam mechanical properties in four-point bending. The proposed optimization of LVL manufacturing from variable quality veneers should help for more efﬁcient usage of forest resources. This optimization strategy showed notable gains for modeled and experimental mechanical properties, whether in terms of stiffness or strength. The analytical calculation of the local modulus of elasticity from modelized beams was satisfactory compared to the tests of the manufactured beams test results, allowing the reliability of the model for this property to be conﬁrmed.


Introduction
The laminated veneer lumber (LVL) is a lamellar composite material made of veneers from the rotary-peeling process. This material has several advantages. First, unlike sawn timber, it allows for a better distribution of defects within the material to avoid defects superimposition and thus homogenizes them in the product. Second, it can allow optimizing the placement of veneers according to their quality depending on the product's performance needs. This is a very important advantage, particularly on very heterogeneous woods in terms of density and number and size of defects, such as large Douglas-fir woods.
In the peeling industry, the sorting of veneers to differentiate their quality is usually performed according to appearance criteria (knot size, slots, resin pockets, etc.) based on EN 635-3 standard [1]. This classification is generally practiced by the intervention of cameras integrated on the panel-making line. In the first part of the series of articles [2], it was shown that local modeling of elastic properties from local fiber orientation and average density could establish quality sorting criteria, giving sorting results significantly different from the standard sorting method. Now, distributing this high-quality heterogeneity of veneers using their local mechanical properties to obtain the best possible performance of the material is a major issue. The scientific literature presents extensive work on the consideration of veneers of different quality in LVL. Until now, a number of works have been aimed at upgrading inexpensive or low-grade veneers with better quality veneers while maximizing the bending properties of the material. Different LVL stacking sequences mixing beech and poplar plies were also tested in edgewise bending by Burdulu et al. [3]. It is not surprising that the configuration with the greatest number of beech plies, the denser species, provided the best performance in terms of bending stiffness and strength. Chen et al. [4] also presented similar findings in a hybrid bamboo-poplar LVL. They also performed flatwise bending tests that highlighted the importance of veneer quality depending on their distance to neutral fiber in this configuration. Hybrid LVL was also tested in [5] with acacia mangium and rubberwood. This was exclusively done by varying the species according to the plies, showing improvements in flatwise bending stiffness and strength when more dense species were used. All these literature results dealt with sequencing different species or qualities of veneers in the thickness of LVL, but to the authors' knowledge, the optimization of the distribution of veneers of heterogeneous qualities along the product length has never been done before. It can be useful when LVL is used as slender beams loaded in edgewise bending.
This article belongs to a group of papers, whose first one is [2], which is dedicated to the use of heterogeneous species, such as large French Douglas-fir, in LVL manufacturing. In this first paper, the determination of local veneer properties by fiber orientation measurement integrated into the peeling line and veneer density enabled the establishment of a criterion to classify them according to their expected mechanical quality. Batches of veneers from the same classes were used to compose virtual panels with a random positioning of veneers. Panels were then virtually cut into beams, which were analyzed according to their equivalent modulus of elasticity to conclude the influence of sorting on the mechanical properties of the lamellar product. The aim was to propose an innovative grading method based on physical-mechanical properties and not just defects observed visually as it is performed for plywood industry (EN 635-3 standard), in order to manufacture reliable structural products.
However, the disadvantage of this random positioning of veneers is if the mechanical properties are rather low on average, consequently the mechanical properties of the LVL material will be low too. The sorting of the veneers can solve this issue but could lead to poor yields. However, in a beam subjected to bending, the stress is not homogeneous; thus, homogenized beams are partly oversized. Therefore, beams could be actually produced with heterogeneous mechanical properties to optimize utilization of raw material: using rather low-quality veneers where stress is low and high-quality veneers where stress is high. Mixing veneer qualities in a LVL beam is a challenge in the development of structural materials, as it means in this present context to consider the variability of the Douglas-fir resource to create beams with the highest possible mechanical properties.
The objectives of the present paper are multiple. First, a new strategy in the use of wood of heterogeneous quality for composing LVL products will be presented, enabling the quality of veneers to be optimized according to the requirement of bending stresses along the beam. This optimization is driven by a declination of a random veneer positioning panel manufacturing algorithm presented in [2]. Second, optimized panels were manufactured from a sample of veneers representative of the considered resource and were tested experimentally according to applicable standards [6,7]. The aim was both to draw conclusions on the model accuracy and to be able to estimate the performance of optimized beams realistically. In addition, to estimate the gains and losses in the mechanical properties of the optimized beams. Random beams from the same veneers were made numerically, simulating manufacturing without the view to optimizing panels. Finally, comparisons with beams from the different sorting methods from [2] were made so that the effects of grading could be estimated against the effects of optimized veneer placement.

Materials and Methods
First, this part explains how the veneers' longitudinal modulus of elasticity was calculated in relation to Part 1 [2], taking into account possible variations in the relationship between modulus of elasticity and density. Then, the algorithm developed to optimize the longitudinal position of the veneers within the laminate according to their quality is presented. Finally, the experimental manufacturing of optimized panels is described, as well as the equations involved in the calculation of mechanical properties, determined experimentally and analytically. Table 1 contains the nomenclature listing the main symbols used.

Pollet Equations
The Pollet Equation [8] was the basis of the LVL beam analytical modeling process. Determined from a Douglas-fir mature, clear wood specimen, it allows the calculation of a bending elastic module from density information. In [2], the longitudinal modulus of elasticity in clear wood of the veneers was calculated using the average veneer density at 12% of moisture content (MC) according to Equation (1). where: • E 0 is the theoretical elastic of modulus for straight grain (MPa) at 12% MC; • ρ veneer is the veneer average density (kg m −3 ) at 12% MC. Then, the local fiber orientation regular grid was resized by applying a 12% MC tangential shrinkage coefficient [9]. Thus, the calculated mechanical properties of randomly composed beams were given at 12% MC to comply with convention.
In this paper, comparing the performance of manufactured LVL beams with modeled beams requires that veneer dimensions be consistent with their local fiber orientation regular grid. Indeed, material moisture content impacts these physico-mechanical properties (size, density, modulus of elasticity . . . ) since the veneers were scanned in a green state, weighed after drying at 9% MC. Moreover, the beam moisture content was determined by double weighing small pieces of them during the destructive test [10], giving 7.5% MC for heartwood beams and 8% for sapwood ones. The Pollet equation (Equation (1)) requires 12% MC corrected density. Local longitudinal modulus of elasticity at 12% MC E x (x, y) was expressed from local fiber orientation regular grid resized in veneer width direction with a 7.5% (respectively 8%) MC tangential shrinkage coefficient for heartwood (respectively sapwood). Thus, the mechanical properties calculated from analytical tests were given at 12% MC. Concerning experimental moduli, they were corrected at 12% MC according to [11].
The experimental scatter used to calculate Equation (1) from this linear regression was very sparse, which raises the question of the accuracy of the formula. Thus, to consider the entire domain of potential modulus/density relationships, some declinations of it were calculated. The slope coefficients were all taken equal to the 36.605 m 2 s −2 value of Equation (1), and several Y-intercept values were recalculated based on the percentage of points below each linear equation, as showed in Figure 1. direction with a 7.5% (respectively 8%) MC tangential shrinkage coefficient for heartwood (respectively sapwood). Thus, the mechanical properties calculated from analytical tests were given at 12% MC. Concerning experimental moduli, they were corrected at 12% MC according to [11]. The experimental scatter used to calculate Equation (1) from this linear regression was very sparse, which raises the question of the accuracy of the formula. Thus, to consider the entire domain of potential modulus/density relationships, some declinations of it were calculated. The slope coefficients were all taken equal to the 36.605 m 2 s −2 value of Equation (1), and several Y-intercept values were recalculated based on the percentage of points below each linear equation, as showed in Figure 1. It means that, in the example of the 75% Pollet equation, the Y-intercept value was calculated in order that 75% of points were under the line. Table 2 presents the coefficients of these equations. Thus, for example, the 75% Pollet equation implies +1052.6 MPa to add to the value of the longitudinal elastic modulus of clear wood. The results from destructive tests presented afterward will be put into perspective with these equations to consider the sensitivity of Equation (1) in the mechanical properties prediction.

Favorable, Unfavorable, and Random Optimization Positioning
The presence of knots between loading points of the neutral axis creates zones of weakness that initiate failure. In addition, low density leads to low static bending strength It means that, in the example of the 75% Pollet equation, the Y-intercept value was calculated in order that 75% of points were under the line. Table 2 presents the coefficients of these equations. Thus, for example, the 75% Pollet equation implies +1052.6 MPa to add to the value of the longitudinal elastic modulus of clear wood. The results from destructive tests presented afterward will be put into perspective with these equations to consider the sensitivity of Equation (1) in the mechanical properties prediction.

Favorable, Unfavorable, and Random Optimization Positioning
The presence of knots between loading points of the neutral axis creates zones of weakness that initiate failure. In addition, low density leads to low static bending strength [8]. The main objective was, therefore, to minimize the presence of zones of mechanical losses due to high knottiness and/or low density between loading points in a four-point bending configuration. The presence of knots in the beam area under tensile solicitation is critical, but this optimization was done in an undifferentiated manner between areas of compression and tensile in order to not limit by the direction of the beam in its use. A model for optimizing the distribution of veneers was developed to improve the mechanical properties of beams according to the bending stress imposed. The principle of optimizing the distribution proceeded by placing the veneers from the center of the panel towards the extremities and selecting the best performing veneer available on a defined criterion. The expected mechanical consequences were, therefore, better strength and higher local elastic moduli due to lower local grain angles in the beam area between the loading points. Figure 2 is a flowchart describing the veneer-positioning algorithm, whose similarities can be observed with the functionality of the random placement algorithm presented in [2]. Boxes representing new functions are framed in red.
As in Part 1 [2], ( → X, → Y) was the coordinate system of the panel. Each of these positions was initially considered as available. All vacant positions were first classified in ascending order according to their eccentricity distance (noted X i,eccentricity,theo in [2]). The first vacant position to be assigned was the one with the lowest absolute value of the eccentricity distance. Then, the highest-criterion veneer of the highest class was selected and allocated to this position. At this stage, the term criterion needed to be defined.
Here, the hierarchy between veneers of the same class was made according E veneer value, defined in [2]. Each of the veneers was then picked to fill each of the vacant positions progressively. A feature of nonsuperposition of areas of low elastic modulus between plies was developed and added to avoid creating local mechanical weaknesses in the width of the material. The principle is explained in Figure 3. Each time a veneer was placed in a vacant position, a test was performed to choose between 0 • or 180 • to minimize defect superimposing. A sum of all stiffness profiles across the 15 plies between the X abscissa of the edges of the considered position (X min and X max ) was performed. A flag system allowed this localized profile to be divided by the number of veneers allocated at this stage on the 15 plies on each X-axis coordinate between X min and the X max : its designation was E x,panel (X). Then, the highest value found between the 2 cases dictated the positioning angle.
As in Part 1 [2] random algorithm, each time a veneer was allocated to a position, the theoretical abscissa of the vacant positions was recalculated from its real length. To simplify the algorithm, a veneer attributed to an end of ply was cut in half as required, and the second part completed the opposite ends. Again, for each one, the condition of no superimposing of defects was applied. Once all veneers had been assigned, and therefore all plies of the material panel were filled, numerical sawing of the beams in the panel was done according to the input setting dimension.
In order to put the mechanical results of the favorably optimized beams into perspective, the working principle of the algorithm previously described was reversed in order to also be able to create the opposite of favorably optimized beams. The principle of assigning positions in ascending order of eccentricity distance from the center of the panel was retained, but the choice by veneer performance was reversed. Veneers were selected and placed from worst to best quality according to prioritization criteria, the E veneer value.  1 2 3 4 5 6 7 8 9 10 11 12 13 14   The Part 1 [2] numerical panel manufacturing algorithm by random veneer positioning was used too in order to evaluate the performance of optimized beams by generating panels and random beams. However, it was not possible to generate all the combinations of possible panels. Indeed, in this study, a panel was composed of 60 full veneers, divided into 72 positions considering veneers at extremities cut in half. Each veneer had two possible positions: 0 • and 180 • . In other words, there were 60! × 2 72 = 3.93 × 10 100 possible combinations. Obviously, not all these configurations could be computer-generated in a reasonable time as it takes one minute per panel. For each batch of veneers, which were used to compose a panel, a sample of 150 combinations was generated per panel to allow some conclusions. A statistical test presented later will show that this hypothesis is relevant.

Composition of Representative Resource Batches
To allow for comparison of the mechanical results of the beams from heartwood and sapwood panels, a batch sampling of the 60 necessary veneers was conducted based on their individual E veneer criterion. The objective was to obtain 2 batches of heartwood and 2 batches of sapwood veneers as similar as possible between them according to this criterion.
For this study, 4 optimized panels were composed: -Heartwood batch n • 1 veneers were used to manufacture a favorable veneer placement optimization without the "defect nonsuperposition" option: dense clear wood was concentrated in the middle of the panel; -Heartwood batch n • 2 veneers were used to manufacture a favorable veneer placement optimization, with the "defect nonsuperposition" option: veneer positioning was driven by their quality and the homogenization of defects by the algorithm portion shown in Figure 3; -Sapwood batch n • 1 veneers were used to manufacture a favorable veneer placement optimization without the "defect nonsuperposition" option: all available good quality wood veneers were concentrated in the center of the panel; -Sapwood batch n • 2 veneers were used to manufacture an unfavorable veneer placement optimization: all available low-quality wood veneers were concentrated in the center of the panel, with the "defect superposition" option: in the final step of the algorithm in Figure 3, the lowest value of E x,panel (X) was chosen to maximize the accumulation of areas of low mechanical property.
Once optimized (favorably or unfavorably), veneer positioning algorithms were used to build the panels, and a virtual cut-out of each beam in their actual dimensions and respecting placement in the real panel was performed.
In order to compare the sapwood and heartwood beams, respectively, the sampling of the veneers had to meet two requirements: the batches, representative of the initial resource, must have an overall equivalent population distribution by wood type. In practice, it is important to remember that this is impossible since each veneer is unique. To compose the veneer samples, the averaged local modulus of elasticity on all veneer surfaces E veneer was used as a quality criterion.
However, this criterion alone was not sufficient to attest to the proper equivalency of batches. The two parameters involved in the determination of the local modulus of elasticity of veneer, which were fiber orientation angle and density (as specified in [2]), were also observed. Table 3 shows these parameters for the complete veneer resource and the 2 sampled batches for each type of wood. According to the Tukey HSD test, when comparing sapwood and heartwood veneers independently, there were no significant differences between batches and no statistically significant difference between batches compared to samples, either on the overall individual veneer density, the average standardized Hankinson value, or on the average veneer modulus of elasticity. This, therefore, ensured that the batches used to create optimized panels were globally equivalent according to these criteria and thus allowed for the comparison of properties between them.

Panel Manufacturing
Four panels were produced according to veneer positioning as identical as possible to the models. A measuring tape was used to check the positioning of the veneers in the → X length of the panel, according to the ( → X, → Y) coordinate system. A stop along the length of the panels enabled positioning the veneers at Y=0. The accuracy of this positioning after gluing was estimated to +/−2.5 mm. Because of this and because the veneer surfaces were not always rectangular, areas of mounted joints (superposition of two veneers that locally cause a double veneer thickness) and cages (local gap between veneers' borders that cause a hole in the panel) were observed, not exceeding 4 mm in length. The gluing process was made with a vacuum press (woodtec Fankhauser GmBh) [12], shown in Figure 4a The glue used was a polyurethane glue (PU) with a spread rate of 200 g·m −2 (Figure 4b). Panels were manufactured two by two (Figure 4c). The glued veneers had a nominal thickness of 3 mm. Their nominal length was between 700 or 750 mm, according to the length of each log, and their cutting width was 850 mm or 1000 mm, adjusted according to the external aspect of the log (visible knottiness, radial cracks). They were pressed together by a diaphragm vacuum press combined with a vacuum pump producing a pressure of 0.08 N mm −2 for 8 h (Figure 4d). After gluing and pressing, the LVL were stacked and stabilized for one week before cutting.
The bending beams were cut from the panels with a panel saw. The nominal dimensions of each beam were 2710 mm × 145 mm × 45 mm. Five beams were cut from each panel. The positioning of beams in the panels was measured to match with the coordinate entered in the algorithm. Because the edges of the panels were not perfectly conspicuous due to the imprecision of veneer placement, an accuracy between 5 mm and 10 mm was estimated in the (X,Y) coordinate positioning of beams in the panel. The bending beams were cut from the panels with a panel saw. The nominal dimensions of each beam were 2710 mm × 145 mm × 45 mm. Five beams were cut from each panel. The positioning of beams in the panels was measured to match with the coordinate entered in the algorithm. Because the edges of the panels were not perfectly conspicuous due to the imprecision of veneer placement, an accuracy between 5 mm and 10 mm was estimated in the (X,Y) coordinate positioning of beams in the panel.

Bending Test Mechanical Properties Assessment
All the specimens were tested in the edgewise direction only. This choice was made to guarantee a sufficient number of bending tests in this direction which is the one used for LVL beams as a slender flexural product. In addition, this test configuration presents two advantages for a fair comparison of the material properties: all the plies were subjected substantially to the same mechanical loading, and the glue joints between plies were globally less subjected to stresses than in flatwise bending. The four-point bending tests were made with a dedicated testing machine composed of an electric actuator equipped with a 100 kN load cell. The upper and lower supports were made of 4.8-cm wide metal plates, fixed on a pivot, allowing the rotation of the beam supports coated with PTFE material.
To get mechanical properties, a four-point bending test was performed on every specimen, modelized and manufactured, based on the EN 14374 [6] and EN 408+A1 [7] standards. A distance equal to 2610 mm, as 18 times the specimen height between the lower supports of 145 mm, was set, as shown in Figure 5. The distance between the loading head and the nearest support (a) was set to 870 mm, which was 6 times the height, in order to prevent possible shear failures. The distance between upper support was also 6 times the height to conform to standards. The destructive test local deflection was measured with a LVDT sensor in 5 times height length portion, between 6.5 h and 11.5 h length from the first lower support.

Bending Test Mechanical Properties Assessment
All the specimens were tested in the edgewise direction only. This choice was made to guarantee a sufficient number of bending tests in this direction which is the one used for LVL beams as a slender flexural product. In addition, this test configuration presents two advantages for a fair comparison of the material properties: all the plies were subjected substantially to the same mechanical loading, and the glue joints between plies were globally less subjected to stresses than in flatwise bending. The four-point bending tests were made with a dedicated testing machine composed of an electric actuator equipped with a 100 kN load cell. The upper and lower supports were made of 4.8-cm wide metal plates, fixed on a pivot, allowing the rotation of the beam supports coated with PTFE material.
To get mechanical properties, a four-point bending test was performed on every specimen, modelized and manufactured, based on the EN 14374 [6] and EN 408+A1 [7] standards. A distance equal to 2610 mm, as 18 times the specimen height between the lower supports of 145 mm, was set, as shown in Figure 5. The distance between the loading head and the nearest support (a) was set to 870 mm, which was 6 times the height, in order to prevent possible shear failures. The distance between upper support was also 6 times the height to conform to standards. The destructive test local deflection was measured with a LVDT sensor in 5 times height length portion, between 6.5 h and 11.5 h length from the first lower support.
The 4-point bending analytical calculation model was adapted from the one developed in [13] in order to switch from flatwise bending to edgewise bending. The 4-point bending analytical calculation model was adapted from the one developed in [13] in order to switch from flatwise bending to edgewise bending.
From the beam equivalent stiffness curve (defined in [2]), the deflection at mid-span of the modeled beams for a four-point bending test of the modeled panels was calculated using the Müller-Breslau's principle (see Equation (2)).
where:  Mf is the bending moment during a four-point bending test induced by the application of unitary load at x position of loading points (MPa);  Mv is the bending moment induced by a unitary load at a given x position (mid-span in this case) (MPa);  EIeff is the effective bending stiffness calculated previously, which is dependent on the local modulus of elasticity (N mm²);  nx is the number of discrete elements along ⃗ direction;  Δx (1 mm) corresponds to the resolution of the interpolated images along the ⃗ direction.
To be consistent with the experimental measurements, the deflection at the 2 support points used to measure local deflection was also calculated according to Equation (2). The mean value of the 2 deflections . and .
were subtracted from the mid-span deflection to get the local relative deflection, as shown in Equation (3).
The experimental and analytical local modulus of elasticity was calculated according to the beam theory in 4-point bending using Equation (4) of Table 4 given in the NF-EN-408+A1 standard [7]. From the EI e f f (x) beam equivalent stiffness curve (defined in [2]), the deflection at mid-span v 1 2 an of the modeled beams for a four-point bending test of the modeled panels was calculated using the Müller-Breslau's principle (see Equation (2)). where: • M f is the bending moment during a four-point bending test induced by the application of unitary load at x position of loading points (MPa); • M v is the bending moment induced by a unitary load at a given x position (mid-span in this case) (MPa); • EI eff is the effective bending stiffness calculated previously, which is dependent on the local modulus of elasticity (N mm 2 ); • n x is the number of discrete elements along → x direction; • ∆ x (1 mm) corresponds to the resolution of the interpolated images along the → x direction.
To be consistent with the experimental measurements, the deflection at the 2 support points used to measure local deflection was also calculated according to Equation (2). The mean value of the 2 deflections v 6.5 18 an and v 11.5 18 an were subtracted from the mid-span deflection to get the local relative deflection, as shown in Equation (3).

(3)
The experimental and analytical local modulus of elasticity was calculated according to the beam theory in 4-point bending using Equation (4) of Table 4 given in the NF-EN-408+A1 standard [7]. where: • l 1 is the gauge length for the determination of modulus of elasticity (mm) • ∆w = w 2 − w 1 is the increment of deformation in millimeters corresponding to ∆F, measured by the local deflection LVDT sensor (mm).
• ∆w = v l,an the mid-span local deflection term is the one calculated by Equation (2).
Contrary to the requirements of the calculation bounds of the modulus for the destructive test, the assumption of the pure elastic behavior of the material considered in the analytical model allowed the local modulus to be calculated on a unitary loading increment.
During the experimental testing of the beams, the bending strength of each beam was measured. The bending strength was calculated according to EN 408+A1 [7] and Equation (5).
It is usual to assume that low stiffness in timber is associated with low strength. Olsson et al. did this in [14] for Norway Spruce timber, also Viguier et al. in [15]. Thus, the lowest bending stiffness computed from the analytical model as defined in Part 1 [2] can be used as a criterion to correlate to the experimental bending strength. The criterion, symbolized E min , which consisted of calculating the minimum value of EI e f f (x) in the between the upper-support portion of beam length (6 times the height of the beam), and then redividing by the nominal moment of inertia to enable an easier reading of results (see Equation (6) of Table 5). • F max is the maximum bending effort given by the load sensor (Newton); • a is the distance between a loading point and the nearest support (mm); • h is the nominal beam height (mm); • a is the distance between a loading point and the nearest support (mm); • I Gz is the moment of inertia for a rectangular cross-section beam (mm 4 ) defined as: • EI e f f (x) is the effective bending stiffness profile along the length of the beam (N mm 2 ); • x ll p is the left loading point x-position in the beam system of axis; • x rl p is the right loading point x-position in the beam system of axis.
It has already been used in literature as indicating property to bending strength [14]. For timbers, a moving average is usually calculated. In the case of the LVL, this value was considered already averaged because the E beam (x, y) regular grid used for the EI e f f (x) calculation was already an average over the 15 plies of the material.

Effect of Beam Optimization on Modeled Mechanical Properties
First, observations can be made on the impact of the controlled distribution of veneers based on their quality within the different panels on their effective bending stiffness.
The colored EI e f f (x) curves presented in Figure 6a are the bending stiffness profiles for each modeled and manufactured beam. They allow for a very explicit observation of the effects of the veneer positioning processes. Favorably optimization led to observing maximum bending stiffness values in the central part of the beams and a regularly decreasing bending stiffness at the ends of the beams. The opposite was also observable for unfavorably optimized beams. The dotted vertical lines represent the loading points of the four-point bending test arrangement. The black curves in Figure 6b,c,f,g represent the EI e f f (x) curves of beams obtained from random positioning of veneers in the panel. The 150 cases per batch simulated represent 750 beams (5 beams cut for each panel). All displayed black curves in Figure 6b,c,f,g, and showed the possible "area of performance". This allowed showing how the optimized bending stiffness profiles provide the best possible stiffness to fit the applied bending moment, while the random arrangements had up and downs randomly all along the beams (shown by individual grey curve examples). Figure 6d,e,h,i are examples of colormaps of the local modulus of elasticity averaged over the 15 plies of LVL beam generated from different batches (E beam (x, y), as defined in Part 1). Figure 6d,h,i shows a concentration of green between loading points, corresponding to higher local modulus than in the beam end areas. There was also a variation in the color range between the sapwood beams (Figure 6d,e) and heartwood ones (Figure 6h,i), explained by a density effect between the two types of wood.
For cases of favorably optimized beams, an EI e f f (x) profile was observed, with higher local values at the center of the beams than at their extremities, which corresponded well to a positioning of the veneers with lower quality at the ends and those of higher quality at the center (Figure 6b,f,g). On the contrary, for unfavorably optimized beams, a concentration of local high EI e f f (x) values at extremities and lower-quality veneers at the center was observed in Figure 6c.
Comparing heartwood and sapwood, a general gap between profiles was observable on the optimized profiles. This difference was mostly imputable to a difference in the veneers' density. To estimate this difference, an overall average of all the points in the profiles of all the sapwood and heartwood beams was performed. The average value was 1.54.10 11 N mm 2 for heartwood and 1.74.10 11 N mm 2 for sapwood (difference of 13.0%).
There was a difference between the positioning of veneers leading to the favorably optimized heartwood curves in Figure 6f,g. For the first one, there was a high local maximum between loading points (slightly above the "area of performance"), but also a low local minimum between loading points compared to the curves of the second batch. This was due to the placement algorithm, modified between the creation of the two batches to improve strength. In the first case, the requirement to place the higher-quality veneers in the center was preponderant, while in the second case, the effect of maximizing the plateau value between the loading points was sought.

Modeling Result Analysis by Comparison with Destructive Results
To analyze the accuracy of the mechanical properties modeling of beams, a comparison was made between the results of destructive testing on the optimized beams and analytical tests on their numerical equivalent. Beam EI e f f profiles: (a) optimized beams from the four batches, (b) batch n • 1 sapwood beams among random beams, (c) batch n • 2 sapwood beams among random beams, (f) batch n • 1 heartwood beams among random beams, (g) batch n • 2 heartwood beams among random beams. Beam E beam (x, y) grid: (d) batch n • 1 sapwood beam, (e) batch n • 2 sapwood beam, (h) batch n • 1 heartwood beam, (i) batch n • 2 heartwood beam.

Modeling Result Analysis by Comparison with Destructive Results
To analyze the accuracy of the mechanical properties modeling of beams, a comparison was made between the results of destructive testing on the optimized beams and analytical tests on their numerical equivalent. E m,l,ex experimental results were recalculated to 12% MC according to EN 384 standard, to be compared with E m,l,an . These results are visible in Figure 7. Whisker plots enabled us to observe the domain of mechanical properties, whose edges represent the value resulting from the application of the 25% and 75% variants of the Pollet equation from Table 2. experimental results were recalculated to 12% MC according to EN 384 standard, to be compared with , , . These results are visible in Figure 7. Whisker plots enabled us to observe the domain of mechanical properties, whose edges represent the value resulting from the application of the 25% and 75% variants of the Pollet equation from Table 2. The relative differences between modelized results from Pollet [8] mature Douglas wood equation (Equation (1)) according to experimental results were on average of −0.9% for heartwood, with a minimum deviation of −7.1% and a maximum deviation of +2.5%. For sapwood, the difference was, on average, −3.8%, with a minimum deviation of −7.5% and a maximum of +1.2%. Thus, despite the fact that this equation is applicable for mature wood, the analytical bending modeling was found to have a tendency to undervalue the local modulus of elasticity, especially for sapwood, and more clearly for unfavorable optimization. Moreover, the model exhibited difficulty perceiving the singularities of the different beams from the same panel. Hence, it tended to be relatively consistent in its predictions (R² = 0.73), but the conclusions cannot be drawn any further given the reduced number of specimens.
The experimental tests showed the optimization effect of the panels on the , , and , results of the four batches (shown in Table 6).  The relative differences between modelized results from Pollet [8] mature Douglas wood equation (Equation (1)) according to experimental results were on average of −0.9% for heartwood, with a minimum deviation of −7.1% and a maximum deviation of +2.5%. For sapwood, the difference was, on average, −3.8%, with a minimum deviation of −7.5% and a maximum of +1.2%. Thus, despite the fact that this equation is applicable for mature wood, the analytical bending modeling was found to have a tendency to undervalue the local modulus of elasticity, especially for sapwood, and more clearly for unfavorable optimization. Moreover, the model exhibited difficulty perceiving the singularities of the different beams from the same panel. Hence, it tended to be relatively consistent in its predictions (R 2 = 0.73), but the conclusions cannot be drawn any further given the reduced number of specimens.
The experimental tests showed the optimization effect of the panels on the E m,l,exp and f m,exp results of the four batches (shown in Table 6). Indeed, the average relative deviation for E m,l,exp between sapwood favorably and unfavorably optimized beams was +11.8% (16,599 MPa vs. 14,847 MPa), while the relative deviation between the E m,l,an mean value was +15.5% (16,205 MPa vs. 14,036 MPa). The relative deviation in between number 2 and number 1 batches of heartwood-favorably optimized beams E m,l,exp mean value was +1.0% (14,665 MPa vs. 14,527 MPa ) while the relative gap between the E m,l,an mean value was −0.4% (14,413 MPa vs. 14,472 MPa). The absolute values were also very close, highlighting the great calibration of analytical model parameters, as noticeable in Figure 7.
Despite unavoidable variations between beams, these average results provided some confidence in the representativeness of the modeling method, allowing for deeper exploitation of modeling results in the following section.
Regarding bending strength in sapwood batches, there was an average relative deviation of +28.1% between the favorably optimized and unfavorably optimized batches. This value is remarkable given that only the positioning of veneers according to their quality differs between the batches. Considering the small sample of beams tested, it was difficult to conclude the predictive aspect of the model by the E min criterion on the strength of beams tested in four-point bending, as shown in Figure 8. between sapwood favorably a unfavorably optimized beams was +11.8% (16,599 MPa vs. 14,847 MPa), while the relati deviation between the , , mean value was +15.5% (16,205 MPa vs. 14,036 MPa). T relative deviation in between number 2 and number 1 batches of heartwood-favorab optimized beams , , mean value was +1.0% (14,665 MPa vs. 14,527 MPa ) while t relative gap between the , , mean value was −0.4% (14,413 MPa vs. 14,472 MPa). T absolute values were also very close, highlighting the great calibration of analytical mod parameters, as noticeable in Figure 7.
Despite unavoidable variations between beams, these average results provided som confidence in the representativeness of the modeling method, allowing for deeper expl tation of modeling results in the following section.
Regarding bending strength in sapwood batches, there was an average relative de ation of +28.1% between the favorably optimized and unfavorably optimized batches. Th value is remarkable given that only the positioning of veneers according to their qual differs between the batches. Considering the small sample of beams tested, it was diffic to conclude the predictive aspect of the model by the criterion on the strength beams tested in four-point bending, as shown in Figure 8. However, the effect of favorable and unfavorable optimizations on the beams of t batches was noticeable, as for local modulus of elasticity. For batch n°2 sapwood unfavo ably optimized beams, the , bending strength was between 40 MPa and 52 MP while for batch n°1 favorably optimized beams, it was between 55 MPa and 67 MPa. Co cerning the criterion, batch n°2 optimized beams averaged 13,091 MPa while bat number 1 optimized beams averaged 15,257 MPa (+16.5%). For heartwood, the two favo ably optimized batches presented similar strength between 39 MPa and 59 MPa. It seem that, even if an increasing global trend was observed, the proposed criterion was not ab to differentiate precisely the observed variation of bending strength for all favorable unfavorable LVL composition.

Comparison with Sorted and Unsorted Random Beams
As a point of comparison to evaluate the performance of optimized beams, rando panels created using the algorithm developed in [2] were generated from the same vene batches. Figures 9 and 10 show respectively , , and of beams with optimiz However, the effect of favorable and unfavorable optimizations on the beams of the batches was noticeable, as for local modulus of elasticity. For batch n • 2 sapwood unfavorably optimized beams, the f m,exp bending strength was between 40 MPa and 52 MPa, while for batch n • 1 favorably optimized beams, it was between 55 MPa and 67 MPa. Concerning the E min criterion, batch n • 2 optimized beams averaged 13,091 MPa while batch number 1 optimized beams averaged 15,257 MPa (+16.5%). For heartwood, the two favorably optimized batches presented similar strength between 39 MPa and 59 MPa. It seems that, even if an increasing global trend was observed, the proposed criterion was not able to differentiate precisely the observed variation of bending strength for all favorable or unfavorable LVL composition.

Comparison with Sorted and Unsorted Random Beams
As a point of comparison to evaluate the performance of optimized beams, random panels created using the algorithm developed in [2] were generated from the same veneer batches. Figures 9 and 10 show respectively E m,l,an and E min of beams with optimized veneer placements compared to random beams made from the same veneers. Violin plots allows observing the distribution law of random values. Dots representing optimized beams enabled us to observe the optimization of the beams in relation to the random sample.
veneer placements compared to random beams made from the same veneers. Violin plots allows observing the distribution law of random values. Dots representing optimized beams enabled us to observe the optimization of the beams in relation to the random sample.
(a) (b) For each sample of 150 random beam results, a Shapiro-Wilk test was performed. For , , and random beams results, a Gaussian distribution was observed (pvalue < 5%). This enables the use of 150 random beam sample trends to be relevant to conclusions in general terms. For results, the distribution was crushed over the higher values. This was due to the nature of the parameters. The local modulus elasticity was calculated for a portion of beam (5 ℎ length). This gave a kind of averaged value. The value was by nature a minimum value for substantially the same portion of beam (6 ℎ length). Table 7 summarizes the main results of the study, including Part 1 results. It compares the results of beams from optimized panels with random panels, from the same batches, and random beams of [2], made of veneers of different sorting classes. The values are presented with mean values, but also by their 5th percentile value, to continue the analogy with the strength property usually used in the construction design. In the case of the manufactured beams, since optimized batches were composed of only five beams maximum, the minimum of values was used.  veneer placements compared to random beams made from the same veneers. Violin plots allows observing the distribution law of random values. Dots representing optimized beams enabled us to observe the optimization of the beams in relation to the random sample.
(a) (b) For each sample of 150 random beam results, a Shapiro-Wilk test was performed. For , , and random beams results, a Gaussian distribution was observed (pvalue < 5%). This enables the use of 150 random beam sample trends to be relevant to conclusions in general terms. For results, the distribution was crushed over the higher values. This was due to the nature of the parameters. The local modulus elasticity was calculated for a portion of beam (5 ℎ length). This gave a kind of averaged value. The value was by nature a minimum value for substantially the same portion of beam (6 ℎ length). Table 7 summarizes the main results of the study, including Part 1 results. It compares the results of beams from optimized panels with random panels, from the same batches, and random beams of [2], made of veneers of different sorting classes. The values are presented with mean values, but also by their 5th percentile value, to continue the analogy with the strength property usually used in the construction design. In the case of the manufactured beams, since optimized batches were composed of only five beams maximum, the minimum of values was used. For each sample of 150 random beam results, a Shapiro-Wilk test was performed. For E m,l,an and E min random beams results, a Gaussian distribution was observed (p-value < 5%). This enables the use of 150 random beam sample trends to be relevant to conclusions in general terms. For E min results, the distribution was crushed over the higher values. This was due to the nature of the parameters. The local modulus elasticity was calculated for a portion of beam (5 h length). This gave a kind of averaged value. The E min value was by nature a minimum value for substantially the same portion of beam (6 h length). Table 7 summarizes the main results of the study, including Part 1 results. It compares the results of beams from optimized panels with random panels, from the same batches, and random beams of [2], made of veneers of different sorting classes. The E min values are presented with mean values, but also by their 5th percentile value, to continue the analogy with the strength property usually used in the construction design. In the case of the manufactured beams, since optimized batches were composed of only five beams maximum, the minimum of E min values was used.
First, a comment on the results of the random heartwood and sapwood batches beam test results can be made. For heartwood, there was a relative deviation of −0.7% between E m,l,an of batch number 2 random beams (13,362 MPa) and batch number 1 random beams (13,462 MPa). For E min , the relative deviation was on average −1.2% (12, (14,359 MPa vs. 14,360 MPa). This allowed drawing conclusions on the correct equivalence of the veneer batches in terms of mechanical property. Concerning E min , it was not proven that this indicator allows an accurate prediction of beam bending strength. A simple proportionality trend was shown. The results involving the E min were therefore indicative.  Table 7 gives information to compare results of favorable and unfavorable optimized beams quantitatively. For sapwood, the favorably optimized beams from batch number 1 had an E m,l,an of 16,205 MPa and an E min value of 15,257 MPa. It was respectively +15.5% and +16.6% than results of batch number 2 unfavorably optimized beams. The results of batch number 1 optimized beams could also be compared with random beams from the same batch. Thus, relative gaps of +7.6% between optimized beam E min minimum value and random beam E min 5th percentile value (13,779 MPa) were noticed. In terms of modulus of elasticity, the gap between means values was +6.5% (16,205 MPa vs. 15,209 MPa). For heartwood, the favorably optimized beams of batch number 2 had a E m,l,an of 14,413 MPa and a E min of 13,658 MPa. The batch number 1 heartwood favorably optimized beams had a E m,l,an of 14,472 MPa and a E min value of 12,829 MPa. This 6.5% difference in E min is explained by the fact that, as seen in Figure 6f,g, these beams concentrated excellent local properties at the center of the beams. However, they also had a local minimum in the upper support area that was lower than batch number 2 heartwood favorably optimized beams, mainly due to the no-overlay defect function. For local modulus of elasticity, it is normal to see that this effect had very little impact (−0.4% between 14,413 MPa and 14,472 MPa) compared to the random beams generated from the same veneers.
On the basis of this better homogenization, batch number 2 favorably optimized beams, E min was +12.5% higher than batch number 2 random beams E min 5th percentile value (13, It was also possible to compare the results of optimized beams composed of the veneers representative of the resource with the random beams composed of the strength classes created in [2], i.e., compare the proposed optimized use of the veneer material to the more classical strength. It can be noted that sapwood favorably optimized beams were stiffer in terms of E m,l,an than all the random beam classes, and more especially than Class A (+2.5% between 16,205 MPa vs. 15,803 MPa). The same trend was true for E min (a minimum value of 14,826 MPa vs. a 5th percentile value of 14,360 MPa). Similarly, heartwood batch number 2 was superior to all classes (E m,l,an , 14,413 MPa vs. 14,127 MPa for Class A, E min minimum value of 13,399 MPa vs. E min 5th percentile value of 12,615 MPa for Class A). For heartwood batch number 1, it is true for E m,l,an (14,472 MPa vs. 14,127 MPa), but E min minimum value was slightly lower than Class A E min 5th percentile value (12,313 MPa vs. 12,615 MPa). This contributes to show the beneficial effect of avoiding superposition defects, which permitted a better distribution of the defects in the zone between supports and thus improved the E min value in heartwood batch number 2 over heartwood batch number 1. Generally, these results prove that the proposed optimization can allow obtaining better bending stiffness from the full heterogeneous resource of peeled veneers than with a high-quality classical sorting and random LVL composing.

Discussion
To the best of authors' knowledge, this study is the first to deal with the optimization of the LVL material by the consideration and exploitation of a highly heterogeneous species, not by the arrangement of the different qualities through the plies, but along the length of the material, according to an edgewise bending testing configuration. This was possible thanks to a new approach based on a numerical model fed with veneer density and local fiber angle in the veneer plan. This last parameter was measured using a new apparatus specifically designed to measure online local veneer fiber orientation.
The effects of the customized were highlighted by favorable or unfavorable optimization and were observed in the analytical and experimental results of the beams of both batches. For tests on the sapwood LVL beams, variations of +11.8% for experimental local modulus of elasticity and +28.1% for bending strength were found on average for favorably optimized batch compared to unfavorable batch.
The results presented in terms of local modulus of elasticity calculation from modelized beams were satisfactory, which means that the prediction of local properties and the analytical calculation are relatively reliable, despite sources of inaccuracies. In terms of pure modeling, the use of veneer global density, neglecting the localized effects of latewood and earlywood, is an important perspective of improvement.
The gain of the favorably optimized beams over random beams was quantifiable thanks to the analytical results, and it was more than 7.5% in terms of stiffness. However, the manufactured beams of this study had a relatively low span. The effectiveness of this method, which consists of coinciding the veneer quality with the curve of the bending moment, remains to be tested in cases of larger beams. Moreover, this method takes advantage of the variability of the veneer resource used by positioning low-quality veneer to low-stressed beam area, and conversely, thus the main advantage is to allow more efficient usage of low-quality resources.
Concerning the use of the nonsuperposition option for heartwood batch number 2, the homogenization effect was not noticeable in experimental results compared to batch number 1. Indeed, in terms of elastic mechanical property, the local modulus averaged the behavior between the loading points of the beam and did not allow for a precise conclusion to be drawn on a phenomenon evolving along this portion. Second, the number of beams was not sufficient to identify trends.
Finally, the modeling of LVL beam strength is a very complex issue. The use of E min as a bending strength indicator did not show very high reliability. It did, however, make it possible to distinguish, on average, the differences in performance between the batches. It was calculated with a very large number of simplifying hypotheses, intrinsically linked to analytical modeling. For example, the conditions and quality of gluing were not taken into account at all. The need to determine an indicator based on other criteria, such as local defect density or average knot spacing, remains. Funding: This study is funded by the region of the Burgundy Franche-Comté, the TreeTrace project subsidized by ANR-17-CE10-0016-03. This study was performed thanks to the partnership built by BOPLI: a shared public-private laboratory built between the Bourgogne Franche-Comté region, LaBoMaP and the Brugère company.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to industrial application confidentiality.
Acknowledgments: Thanks to Roxane Dodeler for her investment in the experimental task during her Master project period. Thanks to Jean-Claude Butaud for the experimental task preparation and execution.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.