Determination of Elastic Properties of Beech Plywood by Analytical, Experimental and Numerical Methods

: This research article examines the application of various methods to determine the e ﬀ ective elastic properties of beech veneer-wood composites. Using laminate theory, the theoretically calculated e ﬀ ective values of the in-plane and out-of-plane modulus of elasticity as well as shear modulus are compared with the values determined from the natural frequencies of ﬂexural, torsional and longitudinal vibrations of samples having di ﬀ erent orientations and numbers of composite layers. The samples are also modelled using the ﬁnite element method, and their natural frequencies are calculated by the modal analysis. Research has shown that the laminate theory, which is well established and applied in the world of synthetic composites, can also be applied to beech plywood composites, where the theoretically calculated e ﬀ ective values can be up to 15% higher. Similarly, due to the higher calculated e ﬀ ective elastic properties, higher natural frequencies of ﬂexural, torsional and longitudinal vibrations are also calculated by the ﬁnite element method. element model of the 7A sample with the ﬁrst layer orientation of 22.5 ◦ , together with the relative deformation of the ﬁrst ﬂexural, torsional and longitudinal modes obtained by modal analysis using the Ansys software. The image clearly shows non-uniform ﬂexural and longitudinal vibrations as a result of the tissue orientation. The ﬂexural vibration image shows that the nodal line is not perpendicular to the long edge, which can lead to certain problems if perpendicular nodal lines are assumed, and specimen is perpendicularly supported at 0.22 and 0.78 L.


Introduction
Plywood is a widespread and long used building material. It is usually made from peeled veneer, where poor quality veneers can be used on the inside and higher quality veneers of various species are used on the visible outside. The good aspect of plywood is undoubtedly the reduction in anisotropy, as different combinations of veneer layers can affect the mechanical properties of the panel. One of the first who experimentally analysed the relationship between the number of layers, fibre angle, modulus of elasticity and strength, was undoubtedly Norris [1,2]. He was later followed by many other researchers. Kallakas et al. [3] investigated the influence of different layer structures and the use of various wood species on the mechanical properties of plywood, while Biadala et al. [4] carried out similar research trying to develop water-resistant plywood with combinations of different wood species. The influence of fibre reinforcement and wood species on the physical and mechanical properties of veneer plywood as well as on laminated veneer lumber was also investigated [5][6][7], as well as the influence of veneer densification on the mechanical properties of plywood [8]. Surface layer polymerisation by active monomers proved to have a positive effect on the modulus of elasticity and strength of wood-based flooring panels [9]. Kral et al. [10] investigated the use of cork in the middle layer of composite plywood panels, while Hao et al. [11] experimented with wood-based sandwich panels with honeycomb structure, and Zhou et al. [12] with the use of bamboo in laminated veneer lumber.
As different combinations of the factors listed can lead to different properties, many authors have dealt with the determination of the properties of plywood already produced. Some have determined the modulus of elasticity, where the comparison between the calculated and experimentally determined modulus of elasticity, using static bending tests has been carried out [13]. A similar study was conducted by Yoshihara [14], in which he investigated the influence of the dimensions of industrially produced plywood on the measured modulus of elasticity using flexural vibration. In another study [15], he compared the shear modulus determined by square-plate twist method with the values calculated from flexural vibrations. He found that the calculated shear modulus from the vibration test was higher than the values obtained with the static square plate twist method. The application of the same technique to full-size composite wood panels was also investigated [16,17] and new equations were developed to calculate the equivalent shear modulus of plywood from the shear moduli of the individual lamellae [18]. Since plywood can also be given spatial forms, the springback effect must be taken into account [19].
Between the methods, static and dynamic tests are the most popular to determine the mechanical properties of wood and wood composites. Among the static tests, the three-and four-point bending tests are used to determine the modulus of elasticity [20] and the plate-twist method is used to determine the shear modulus [21]. In the dynamic bending tests, the specimen can be clamped as a cantilever [22][23][24], whereby the boundary conditions, i.e., the stiffness of the clamping, must be taken into account, or the specimen can be free at both ends. As with both cantilever and free-free specimens, inertia and shear effects must be taken into account for smaller length-to-height ratios [25][26][27][28]. The influence of shear also depends on the flexural vibrational mode number [29] and the ratio between the modulus of elasticity and the shear modulus [29,30]. For example, the influence is greater in spruce than in beech, because in the first case the ratio can be up to 25, while in the second case it is 13 [31]. In most of the tests performed, the mechanical properties of the dynamic tests were higher than those of the static tests [32][33][34]. However, even more obvious were the differences between the mechanical properties obtained at very high strain rates and those obtained in the quasi-static tests; the former differences are much higher [35].
Today, many artificial materials [36] are known, such as carbon fibres, Kevlar, and glass fibres, which are characterised by very good mechanical properties in the main direction of orientation. For their modelling, the laminate theory [37,38] is used, with the help of which the mechanical properties of a composite material made of the mentioned materials can be calculated. In addition, the application of the finite element method is very widespread or indispensable, especially for complicated shapes and loads. While these approaches are already well established and used for the abovementioned synthetic materials, their application in modelling wood composites is not so widespread [39,40].
Plywood is nowadays a very common material in the furniture industry, especially in contemporary plywood furniture design [41], where it is used as a planotropic material due to its reduced anisotropy. It can therefore be used in the manufacture of a chair [41], where the front and back legs of one side of the chair, together with the crossbar and the support for the backrest, are cut from one piece of plywood. Such an element requires knowledge of the material properties of the plywood panel at various angles, both in the direction of orientation of the face material and at certain angles up to right angles. If a designer wants to calculate the dimensions of such an element, they must know the material properties in all mentioned directions and not only in the direction of orientation of the outer veneer layer. There are also well known examples of the use of plywood in the manufacture of table legs [41], where a single block has a certain dimension both in the direction of the outer veneer layer and in the rectangular direction. In case of insufficient stiffness in the cross direction or at a certain angle to the longitudinal direction of the outer veneer layer, the final product does not provide satisfactory stiffness, so that the end user of the product has no feeling of furniture stability.
The literature presented describes many individual investigations that were carried out to explain the influence of various factors on the mechanical properties of plywood. What they all have in common is that the results are based on experimental findings and are usually not supported by analytical calculations, such as the application of laminate theory. To the authors' knowledge, there is no research work that deals comprehensively with the determination of the elastic properties of plywood of different compositions, as this study attempts to show. The aim of this work is therefore to investigate the application of laminate theory and to combine individual approaches in order to model the elastic properties of beech plywood as a composite material. Constructive factors such as the number of layers and the orientation of the tissue in individual layers as well as the method for determining and calculating the elastic-mechanical properties are included. A 7-and 11-layer plywood panel made of beech veneer with uniform known mechanical properties were produced (Figure 1). Samples with different combinations of tissue orientation were cut from the board, the mechanical properties of which were then determined on the basis of free vibration of the samples. The mechanical properties were also determined theoretically by the laminate theory, and compared to experimentally determined values. Each sample was modelled by the finite element method using the composite principle, and then its natural frequency was calculated by modal analysis. The calculated frequency was then also compared to the measured frequency. In this way, different approaches such as the method of laminate theory and the finite element method, in the case of using wood in composite materials, were verified. The aim of this work is therefore to investigate the application of laminate theory and to combine individual approaches in order to model the elastic properties of beech plywood as a composite material. Constructive factors such as the number of layers and the orientation of the tissue in individual layers as well as the method for determining and calculating the elastic-mechanical properties are included. A 7-and 11-layer plywood panel made of beech veneer with uniform known mechanical properties were produced ( Figure 1). Samples with different combinations of tissue orientation were cut from the board, the mechanical properties of which were then determined on the basis of free vibration of the samples. The mechanical properties were also determined theoretically by the laminate theory, and compared to experimentally determined values. Each sample was modelled by the finite element method using the composite principle, and then its natural frequency was calculated by modal analysis. The calculated frequency was then also compared to the measured frequency. In this way, different approaches such as the method of laminate theory and the finite element method, in the case of using wood in composite materials, were verified. The weakness of the approach shown lies in the possible lack of knowledge of the input data, such as elasticity and shear modulus in the main directions of the wood tissue. In this case, the user could use the mean values known from the publicly available literature [42]. An additional weakness is also the variability of material properties, where the properties within a tree species can change by up to several 10%, which means that the calculated material properties for a given plywood combination vary accordingly. For a quantitative uncertainty determination, it would be necessary to perform an uncertainty analysis, whereby for a given case, if the mean values and standard deviations of the input material properties were known, the output material variations would be calculated under a given confidence interval. However, since the abovementioned analysis would go beyond the scope of this article, it is not performed in the present work. The weakness of the approach shown lies in the possible lack of knowledge of the input data, such as elasticity and shear modulus in the main directions of the wood tissue. In this case, the user could use the mean values known from the publicly available literature [42]. An additional weakness is also the variability of material properties, where the properties within a tree species can change by up to several 10%, which means that the calculated material properties for a given plywood combination vary accordingly. For a quantitative uncertainty determination, it would be necessary to perform an uncertainty analysis, whereby for a given case, if the mean values and standard deviations of the input material properties were known, the output material variations would be calculated under a given confidence interval. However, since the abovementioned analysis would go beyond the scope of this article, it is not performed in the present work.

Plyboards Fabricaton
For the purpose of the experiment, a custom-made peeled beech (Fagus sylvatica) veneer with tangential texture and nominal thickness of 1.5 mm, made from a single log with uniform annual growth, was purchased. Veneer sheets measuring 1400 × 600 mm were first conditioned in the laboratory at a constant temperature of 22°C and a relative humidity of 45% to an average equilibrium moisture content of 6.7%.
For the production of plywood, a veneer without visual defects and homogeneous structure was chosen. Seven-and 11-layer plywood boards of 600 × 600 mm in size, with different as well as identical veneer orientations, were produced (Table 1, Figure 2). The latter was used to determine the material properties in the principal material directions. Because the veneer was peeled from a cylindrical log and had a strictly tangential orientation in width, it can be described as an orthotropic material with a rectangular coordinate system of principal material properties [43,44]. Table 1. Tissue orientations of individual layers for 7-and 11-layer beech veneer plywood.

Plyboards Fabricaton
For the purpose of the experiment, a custom-made peeled beech (Fagus sylvatica) veneer with tangential texture and nominal thickness of 1.5 mm, made from a single log with uniform annual growth, was purchased. Veneer sheets measuring 1400 × 600 mm were first conditioned in the laboratory at a constant temperature of 22 ℃ and a relative humidity of 45% to an average equilibrium moisture content of 6.7%.
For the production of plywood, a veneer without visual defects and homogeneous structure was chosen. Seven-and 11-layer plywood boards of 600 × 600 mm in size, with different as well as identical veneer orientations, were produced (Table 1, Figure 2). The latter was used to determine the material properties in the principal material directions. Because the veneer was peeled from a cylindrical log and had a strictly tangential orientation in width, it can be described as an orthotropic material with a rectangular coordinate system of principal material properties [43,44].
The melamine-urea-formaldehyde adhesive Meldur H97 was used for production of the panels, to which 1% NH4Cl catalyst and 5% filler (rye flour) were added to increase the viscosity. The mixture was then stirred with an electric mixer for 15 min until a homogeneous mixture was obtained. The adhesive application was 180 g/m 2 , and the compression pressure and temperature were 1.6 MPa and 130 ℃, respectively, while the compression time for 7-and 11-layer plywood panels was 10 and 13 min, respectively. Afterwards the boards were stacked, weighed and conditioned for 1 week. In this way nine 7-layer plywood boards and twelve 11-layer plywood boards were produced.    The melamine-urea-formaldehyde adhesive Meldur H97 was used for production of the panels, to which 1% NH 4 Cl catalyst and 5% filler (rye flour) were added to increase the viscosity. The mixture was then stirred with an electric mixer for 15 min until a homogeneous mixture was obtained. The adhesive application was 180 g/m 2 , and the compression pressure and temperature were 1.6 MPa and 130°C, respectively, while the compression time for 7-and 11-layer plywood panels was 10 and 13 min, respectively. Afterwards the boards were stacked, weighed and conditioned for 1 week. In this way nine 7-layer plywood boards and twelve 11-layer plywood boards were produced. After conditioning, the boards were cut to the size of 550 × 550 mm, and their thickness was measured. The 7-layer board was 9.9 mm thick, while the 11-layer board was 15.6 mm thick, corresponding to 1.42 mm per layer.

Determination of Elastic Properties
This section is divided into four parts. The first part presents the experiments that were carried out to determine the material properties used in the finite element and analytical models using modal analysis and laminate theory, respectively. The second part describes the experimental determination of the effective elastic properties of plywood with different wood tissue orientations of the individual layers based on the natural frequencies of flexural, torsional and longitudinal vibrations. In the third part, the effective elastic properties are calculated using the laminate theory, while the fourth part describes the modal analysis using the finite element method, which is used to calculate the natural frequencies in the vibration modes mentioned above.

Determination of Elastic Properties of Beech Veneer in Principal (L-Longitudinal, T-Tangential, R-Radial) Directions by Flexural and Torsional Vibrations
From 11-layer boards (11E) with the same orientation of all veneer layers, 3 samples with dimensions 410 × 40 × 15.6 mm were cut in L × T × R direction and 3 samples were cut in T × L × R direction to determine the elastic and shear moduli in the principal tissue directions.
Since care was taken during the production of the boards to ensure that all veneers were oriented in the same direction and also had uniform mechanical properties, as they were made from a log with uniform annual growth, we assumed that 3 samples of the same veneer-orientation combination would be sufficient to determine the principal properties. In case of significant differences in the test results of the three samples, additional samples would be made from the boards.
The modulus of elasticity and shear modulus were determined from flexural and torsional frequencies of specimens with both ends free. Since the vibration amplitudes are the highest in the first flexural vibration mode, the specimen should be supported at the corresponding nodal lines, which are for the first flexural vibration mode at 0.22 and 0.78 L. To reduce the influence of the supports, the specimens were supported on elastic ropes in the lengths of 0.22 and 0.78 L and initially excited in the transverse direction by means of a hammer so that the specimens were subjected to flexural and torsional vibrations at their natural frequencies ( Figure 3). Then, the samples were further excited in the longitudinal direction so that they vibrated in the longitudinal direction at their longitudinal natural frequencies.
After conditioning, the boards were cut to the size of 550 × 550 mm, and their thickness was measured. The 7-layer board was 9.9 mm thick, while the 11-layer board was 15.6 mm thick, corresponding to 1.42 mm per layer.

Determination of Elastic Properties
This section is divided into four parts. The first part presents the experiments that were carried out to determine the material properties used in the finite element and analytical models using modal analysis and laminate theory, respectively. The second part describes the experimental determination of the effective elastic properties of plywood with different wood tissue orientations of the individual layers based on the natural frequencies of flexural, torsional and longitudinal vibrations. In the third part, the effective elastic properties are calculated using the laminate theory, while the fourth part describes the modal analysis using the finite element method, which is used to calculate the natural frequencies in the vibration modes mentioned above. Since care was taken during the production of the boards to ensure that all veneers were oriented in the same direction and also had uniform mechanical properties, as they were made from a log with uniform annual growth, we assumed that 3 samples of the same veneer-orientation combination would be sufficient to determine the principal properties. In case of significant differences in the test results of the three samples, additional samples would be made from the boards.
The modulus of elasticity and shear modulus were determined from flexural and torsional frequencies of specimens with both ends free. Since the vibration amplitudes are the highest in the first flexural vibration mode, the specimen should be supported at the corresponding nodal lines, which are for the first flexural vibration mode at 0.22 and 0.78 L. To reduce the influence of the supports, the specimens were supported on elastic ropes in the lengths of 0.22 and 0.78 L and initially excited in the transverse direction by means of a hammer so that the specimens were subjected to flexural and torsional vibrations at their natural frequencies ( Figure 3). Then, the samples were further excited in the longitudinal direction so that they vibrated in the longitudinal direction at their longitudinal natural frequencies.
The vibration of the samples was measured with a Bruel & Kjaer Type 4939 microphone, NI USB 6361 measurement card and LabVIEW software at a sampling rate of 100 kHz, with each measurement repeated 10 times. From each time signal, an FFT (Fast Fourier transformation) was performed and linearly averaged to obtain the averaged frequency spectrum from which the natural frequencies fn for different vibration modes were extracted.   linearly averaged to obtain the averaged frequency spectrum from which the natural frequencies f n for different vibration modes were extracted.
From the measured flexural natural frequencies of 1 to 4 vibration modes, the out-of-plane (OP) modulus of elasticity in the longitudinal direction (E L ), corresponding for the flexural stiffness, and the in-plane (IP) shear modulus (G LR ) were determined from L × T × R (length × width × thickness) samples, while modulus of elasticity in the tangential direction (E T ) and the shear modulus (G TR ) were determined from T × L × R samples, using the Timoshenko equation [34] where E x and G xz are the modulus of elasticity and shear modulus, in x-z direction of the sample coordinate system according to Figure 3, respectively; z is transverse displacement; I, S and ρ are moment of inertia, cross-sectional area and material density of the beam, respectively. s is the shear factor, which depends on the geometry of cross-section and the material properties, where the value of 1.2 was used in this case. The modulus of elasticity and the shear modulus were calculated by the linear regression of the equation [30,45] fi is the ith flexural natural frequency, and parameters m, F 1 (m) and F 2 (m) are calculated based on index i for the ith natural frequency.
The out-of-plane (OP) shear modulus G xy (G LT ), corresponding for out-of-plane torsional vibration, was calculated from an L × T × R sample with a rectangular cross-section with the length, width and thickness of L, b and h, respectively, according to the following equation [46][47][48] f t,j is the jth torsional frequency, ρ specimen density and G xz in-plane shear modulus. Since the shear modulus G xy is already involved in the calculation of K t (Equation (5)), the approximate shear modulus from the literature can be used to calculate K t the first time. The result is then used in Equation (3) to calculate the provisional value of G xy , which is then used again in Equation (5). The procedure is repeated until G xy converges to a constant value. Usually 5 to 10 iterations are required to obtain the converged value of G xy . Equation (5) also needs the shear modulus G xz (G LR ), which is determined from flexural vibrations (Equation (2)). The in-plane (IP) modulus of elasticity in the x direction (E L ), corresponding to in-plane longitudinal stiffness, was also determined from the first natural frequency of longitudinal vibration using the equation [33] where f x,j is jth natural frequency, and ρ and L are specimen density and length, respectively.
where , is jth natural frequency, and ρ and L are specimen density and length, respectively.  As in the case of the specimens with the same tissue orientation of all layers, the specimens were freely supported on elastic ropes in lengths of 0.22 and 0.78 L and excited in the transverse direction with a hammer, so that the specimens were subjected to flexural and torsional vibrations at their natural frequencies (Figure 3). The samples were then further excited in the longitudinal direction so that they vibrated in the longitudinal direction at their natural longitudinal frequencies. Elastic properties of each sample were then determined using Equations (2)- (6).
Since the properties in the x-direction vary from layer to layer, the constant value of 1.2 for the shear factor s in Equation (2) is not valid, since it can deviate significantly from the value of 1.2 [29,49] and a new shear factor was then calculated for each combination of layers using the equation [50] = , ( , ) where S is the cross-section of the sample, b is the width, Gxz,eff is the effective shear modulus of the whole sample, Ex,i is the modulus of elasticity of the ith veneer layer in the x direction, and Ii is the moment of inertia of the ith veneer layer, while Js can be written As in the case of the specimens with the same tissue orientation of all layers, the specimens were freely supported on elastic ropes in lengths of 0.22 and 0.78 L and excited in the transverse direction with a hammer, so that the specimens were subjected to flexural and torsional vibrations at their natural frequencies (Figure 3). The samples were then further excited in the longitudinal direction so that they vibrated in the longitudinal direction at their natural longitudinal frequencies. Elastic properties of each sample were then determined using Equations (2)- (6).
Since the properties in the x-direction vary from layer to layer, the constant value of 1.2 for the shear factor s in Equation (2) is not valid, since it can deviate significantly from the value of 1.2 [29,49] and a new shear factor was then calculated for each combination of layers using the equation [50] Forests 2020, 11, 1221 8 of 21 where S is the cross-section of the sample, b is the width, G xz,eff is the effective shear modulus of the whole sample, E x,i is the modulus of elasticity of the ith veneer layer in the x direction, and I i is the moment of inertia of the ith veneer layer, while J s can be written G xz,i is the shear modulus of the ith layer, n is the total number of layers of the sample, and y j dS is the static moment of area for the jth layer.
G xz,eff can be calculated using equation [51] 1 where for G xz,i and E x,i calculations, Equation (A8) (Appendix A) can be used h and h i are the total thickness of plywood and thickness of ith layer, respectively.
Since it was assumed that the modulus of elasticity, the shear modulus and Poisson's ratios in the principal directions were not known prior to determination of the elastic properties of plywood, literature data were used to calculate the shear factor as shown in Table 2 [31]. Table 2. Elastic properties for beech wood [31], taken for shear factor calculation. Samples with different tissue orientations of individual layers, for which the effective elastic properties were determined experimentally on the basis of flexural, torsional and longitudinal vibration, were also determined analytically using laminate theory. The effective out-of-plane (OP) laminate modulus of elasticity and out-of-plane shear modulus in x-y direction can be calculated using the following equation [37,38]  where constants a * 11 , a * 22 and a * 66 are elements of the matrix [A], also explained in detail in Appendix A (Equation (A21)).

Modal Analysis Using the Finite Element Method
In the fourth part of the experiment, all samples were modelled with the finite element method as composite material with 7 and 11 layers using the Composite module of Ansys software where the materials properties were taken from first part of the experiment. The length of the samples was 410 and 270 mm for 11-and 7-layer samples respectively, the width of the samples was 40 mm and the thickness of each layer was 1.42 mm. The size of the element in xand y-direction of the sample was 5 mm, while in z-direction the size was equal to the layer thickness, i.e., 1.42 mm. A test with a smaller element dimension was also carried out, but the results were the same as the given dimension. The flexural natural frequencies for 1-4 mode, the torsional frequency for the first mode and the first natural frequency for the longitudinal mode were calculated by modal analysis and then compared to the measured natural frequencies.

Results and Discussion
The average density of all samples was 755.3 kg/m 3 with a standard deviation of 12.9 kg/m 3 . Since the standard deviation is small, the same density was considered for all samples. Figure 5a shows an averaged frequency spectrum of the flexural vibration of the 11-layer sample with a tissue orientation of the top layer of −45 • , where the peaks for the first four flexural vibration frequencies and the first torsional frequency at transverse excitation of the sample are visible. The flexural frequencies have very high peaks, while the torsional vibration frequency is less pronounced and increases when the sample is excited at the corner. For the same specimen, the longitudinal vibration spectrum is shown in Figure 5b, where two distinct peaks belonging to the first and second mode of the longitudinal vibration are clearly visible. Only the first longitudinal natural frequency of the vibration was used to calculate the effective in-plane modulus of elasticity in the x direction.
Forests 2020, 11, x FOR PEER REVIEW 9 of 20 thickness of each layer was 1.42 mm. The size of the element in x-and y-direction of the sample was 5 mm, while in z-direction the size was equal to the layer thickness, i.e., 1.42 mm. A test with a smaller element dimension was also carried out, but the results were the same as the given dimension. The flexural natural frequencies for 1-4 mode, the torsional frequency for the first mode and the first natural frequency for the longitudinal mode were calculated by modal analysis and then compared to the measured natural frequencies.

Results and Discussion
The average density of all samples was 755.3 kg/m 3 with a standard deviation of 12.9 kg/m 3 . Since the standard deviation is small, the same density was considered for all samples. Figure 5a shows an averaged frequency spectrum of the flexural vibration of the 11-layer sample with a tissue orientation of the top layer of −45°, where the peaks for the first four flexural vibration frequencies and the first torsional frequency at transverse excitation of the sample are visible. The flexural frequencies have very high peaks, while the torsional vibration frequency is less pronounced and increases when the sample is excited at the corner. For the same specimen, the longitudinal vibration spectrum is shown in Figure 5b, where two distinct peaks belonging to the first and second mode of the longitudinal vibration are clearly visible. Only the first longitudinal natural frequency of the vibration was used to calculate the effective in-plane modulus of elasticity in the x direction.  Table 3 shows the mechanical properties determined from the flexural, torsional and longitudinal natural frequencies together with the coefficient of variation. The modulus of elasticity differs very little between individual specimens, as do the shear moduli. Likewise, the mean principal out-of-plane modulus of elasticity from flexural vibration differs from the in-plane modulus of elasticity from longitudinal vibrations in the L and T directions by only 3.5 and 2.5%, respectively, and it is comparable to values from the literature [52][53][54][55]. Due to the small differences, it was decided that three samples are sufficient for determination of the basic mechanical properties in the main tissue directions (L, T, R), therefore the number of samples was not increased.   Table 3 shows the mechanical properties determined from the flexural, torsional and longitudinal natural frequencies together with the coefficient of variation. The modulus of elasticity differs very little between individual specimens, as do the shear moduli. Likewise, the mean principal out-of-plane modulus of elasticity from flexural vibration differs from the in-plane modulus of elasticity from longitudinal vibrations in the L and T directions by only 3.5 and 2.5%, respectively, and it is comparable to values from the literature [52][53][54][55]. Due to the small differences, it was decided that three samples are sufficient for determination of the basic mechanical properties in the main tissue directions (L, T, R), therefore the number of samples was not increased.  The effective out-of-plane modulus of elasticity (E x,OP,eff ) is greatest in the orientation of the tissue of the outer layer of 0 • , as it contributes most to the stiffness due to the longitudinal modulus of elasticity (E L ) of the wood tissue, while inwardly the influence of the layers decreases. However, the situation is different with the effective in-plane modulus of elasticity in x direction, where all layers have the same effect.

Effective Elastic Properties in Specimen with Various Layer Orientations
The effective theoretical values calculated using the laminate theory are, in all cases, greater than the measured moduli, i.e., the ratio between the calculated and the average measured value is between 1.03 and 1.25, or the average values are 1.15 and 1.12 for the out-of-plane and in-plane modulus of elasticity of 11-ply samples, respectively, and 1.12 and 1.08 for 7-ply samples, respectively. The latter means that when dimensioning the plate, it must be considered that the calculated effective values will be, on average, 8 to 15% higher than the actual values.
The situation is different for the effective shear moduli, where the measured shear moduli between samples with the same combinations scatter very little in some cases, and in other cases very much. The reason for this is that for certain combinations of tissue orientation, it was difficult to determine the correct flexural and torsional frequency because they differed very little and it was not entirely clear which belonged to the flexural mode and which was for the torsional mode. In some cases, it was also very difficult to induce torsional vibration. This leads to an error in the calculation of both the G xz and G xy shear modulus. In some combinations of layer orientations, where the ratio of E x,OP,eff /G xz,IP,eff is small, the accuracy of the determination of G xz,IP,eff is also lower [29,30,46] which consequently produces a lower accuracy of G xy,OP,eff . This is the case of the first ply orientation of −45 • for 11-layer samples of P-type plates, where the average out-of-plane modulus of elasticity amounts 3059 MPa and the average in-plane shear modulus 937 MPa. In this case the ratio is only 3.26, and the accuracy of determined G xz,IP,eff is very low, which is seen as high scatter in Figure 6. However, in cases where the frequencies were clearly identifiable and the ratio between elasticity and shear modulus was correspondingly high, such as in the case of the first ply orientation of 0 • , where the average out-of-plane modulus of elasticity amounts to 10,500 MPa and the average in-plane shear modulus is 600 MPa, the ratio increases to 17.5, and correspondingly also increases the accuracy of the determined shear modulus.
The calculated effective shear moduli are comparable to the measured moduli and are slightly higher in some places and slightly lower in others. The ratios between the calculated and average measured values for the 11-layer plate vary between 0.59 and 1.25 for G xz,IP,eff and between 0.62 and 1.54 for G xy,OP,eff and are on average 1.07 and 0.94 for G xy,OP,eff and G xz,IP,eff , respectively. Similarly, the ratios for 7-layer plates vary between 0.63 and 1.47 for G xz,OP,eff and between 0.17 and 1.17 for G xy,OP,eff , or, on average, 0.97 and 0.76 for G xy,OP,eff and G xz,IP,eff , respectively. The high variation of ratios in some tissue orientations can be explained by the fact that at lower values of the E/G ratio, the accuracy of the shear modulus determination is lower, as already explained.     Figure 8 shows the calculated shear factor s for specimens with various combinations of tissue orientation, considering the modulus of elasticity, shear modulus and Poisson values from the literature. The figure shows that the shear factor is in the range of 0.9 to 1.65, which means that if the variability of the cross-sectional composition is not considered, a large error can be made. In our case, the shear factor was also calculated with the measured data, but the difference between the value calculated from literature data and the value calculated from measurement data for a given sample was only a few percent. This confirms the assumption that in the absence of knowledge of the material  Figure 8 shows the calculated shear factor s for specimens with various combinations of tissue orientation, considering the modulus of elasticity, shear modulus and Poisson values from the literature. The figure shows that the shear factor is in the range of 0.9 to 1.65, which means that if the variability of the cross-sectional composition is not considered, a large error can be made. In our case, the shear factor was also calculated with the measured data, but the difference between the value calculated from literature data and the value calculated from measurement data for a given sample was only a few percent. This confirms the assumption that in the absence of knowledge of the material properties required to determine the shear factor, values from the literature can be used. This results in an error of only a few percent, which is significantly smaller than in the case of using a constant value of 1.2. Figure 9 shows a composite finite element model of the 7A sample with the first layer orientation of 22.5 • , together with the relative deformation of the first flexural, torsional and longitudinal modes obtained by modal analysis using the Ansys software. The image clearly shows non-uniform flexural and longitudinal vibrations as a result of the tissue orientation. The flexural vibration image shows that the nodal line is not perpendicular to the long edge, which can lead to certain problems if perpendicular nodal lines are assumed, and specimen is perpendicularly supported at 0.22 and 0.78 L.
Forests 2020, 11, x FOR PEER REVIEW 13 of 20 properties required to determine the shear factor, values from the literature can be used. This results in an error of only a few percent, which is significantly smaller than in the case of using a constant value of 1.2.    in an error of only a few percent, which is significantly smaller than in the case of using a constant value of 1.2.   . Composite finite element model using Ansys software for a specimen 7A with first ply orientation of 22.5°: (a) first layer with specimen coordinate system and designation of the main direction (yellow arrow) and first ply orientation (green arrow); (b) relative deformation of the first flexural vibration mode in specimen coordinate system; (c) relative deformation of the first torsional vibration mode in specimen coordinate system; (d) relative deformation of the first longitudinal vibration mode in specimen coordinate system. Table 4 shows the average measured and calculated natural frequencies using the finite element method with the Ansys software for vibration mode 1-4 of the flexural vibration and the first mode Figure 9. Composite finite element model using Ansys software for a specimen 7A with first ply orientation of 22.5 • : (a) first layer with specimen coordinate system and designation of the main direction (yellow arrow) and first ply orientation (green arrow); (b) relative deformation of the first flexural vibration mode in specimen coordinate system; (c) relative deformation of the first torsional vibration mode in specimen coordinate system; (d) relative deformation of the first longitudinal vibration mode in specimen coordinate system. Table 4 shows the average measured and calculated natural frequencies using the finite element method with the Ansys software for vibration mode 1-4 of the flexural vibration and the first mode of the torsional and longitudinal vibration. The table shows that the frequencies of the 11-layer sample with the same tissue orientations (11E-0 • and 11E-90 • ) are in very good agreement between the measured and calculated values for the flexural as well as for the longitudinal and torsional vibration.
However, the difference is somewhat greater for different combinations of tissue orientations of individual layers, where the values calculated using the finite element method are somewhat higher.
The reason for this is that the Ansys program calculates a slightly higher effective modulus of elasticity and shear modulus, similar to the calculation of effective values using laminate theory.

Conclusions
The research showed the application of the experimental determination of the modulus of elasticity and shear modulus from the natural frequencies of flexural, longitudinal and torsional vibrations of plywood samples with various combinations of the orientation of the individual layers. It has been shown that the method is particularly useful for slender and wide elements where it is possible to excite and accurately determine both flexural and torsional frequencies from which the modulus of elasticity and shear moduli can be calculated.
It was also found that the accuracy of the determination of the shear modulus is strongly dependent on the ratio of modulus of elasticity to shear modulus, since the ratios vary greatly for different combinations of tissue orientations of the individual layers. For example, the ratio varied in the experiment from a value of 3.26, where the accuracy of the shear modulus determined was very low, as evidenced by a high scatter of the individual values, to 17.5, where the accuracy of the modules was much higher and the scatter was negligible.
The research also confirmed that the laminate theory, which is already well established and widely used for composites of synthetic materials, can be used for modelling plywood in all tissue combinations. In this way, it is possible to determine the effective out-of-plane flexural and in-plane longitudinal modulus of elasticity, as well as in-plane and out-of-plane shear moduli. Research has shown that the calculated effective modulus of elasticity can be on average up to 15% higher than the actual one, which must be considered when using such composites in practice.
When applying the laminate theory to plywood, the weakness of the theory could be the lack of input data and also the variability of the data for a particular tree species, which must also be taken into account. As mentioned in the introduction, in the absence of precise material properties and taking into account the data from the literature and the fact that the variability of wood properties in the main directions can be up to several 10%, the calculated values can also vary accordingly.
The situation is similar for modal analysis using the finite element method, where the software calculates higher effective values and thus higher natural frequencies. So, it can be said that the applicability of the finite element method for planar elements is equal to the applicability of the laminate theory, but the applicability for spatial elements, where the classical analytical methods can hardly provide an accurate calculation, is greatly increased. The linear elastic stress-strain Hook's law constitutive relationship can be written as [37,38] where C ijkl is a fourth-order tensor with 81 elastic constants, and the inverted form of Hooke's law in the case of orthotropic material Coefficients S ij are called compliance coefficients which is the inverse of the stiffness matrix The constitutive Equation (A2) for an orthotropic material can then be written in matrix form as where individual compliance coefficients in the terms of the engineering constants are E 1 , E 2 , E 3 are the modulus of elasticity in the principal material (1-2-3) coordinate system ( Figure A1), G 12 , G 13 , G 23 , shear modulus in and ν ij Poisson's ratios. When the principal material coordinate system does not coincide with the global x-y-z coordinate system, but it is rotated through an angle θ about the common z-(3) axis, as indicated in Figure A1, the transformed compliance matrix can be defined as and u= cos(θ) and v= sin(θ). From transformed compliance coefficients the desired modulus of elasticity and shear modulus in x-y-z directions for each separate layer can be calculated as follows It the case of thin composites for which a condition of plane stress exists or is a very good approximation, an orthotropic material with principal material coordinates, plane stress constitutive Equation (A4), has the following simplified form for reduced stiffness Figure A1. Multidirectional laminate with coordinate notation of individual plies, subjected to forces and moments per unit length.
When the principal material coordinate system does not coincide with the global x-y-z coordinate system, but it is rotated through an angle θ about the common z-(3) axis, as indicated in Figure A1, the transformed compliance matrix can be defined as where the transformation matrices [T 1 ] and [T 2 ] are and u= cos(θ) and v= sin(θ).
From transformed compliance coefficients the desired modulus of elasticity and shear modulus in x-y-z directions for each separate layer can be calculated as follows It the case of thin composites for which a condition of plane stress exists or is a very good approximation, an orthotropic material with principal material coordinates, plane stress constitutive Equation (A4), has the following simplified form for reduced stiffness and for reduced compliance where the stiffness is the inverse of the compliance In the case of rotation of x axis about the 3 (z) axis using [T 1 ] transformation matrix for stress and [T 2 ] for transformation of strain, we get the plane stress transformed reduced stiffness matrix and The stresses at any z-location of the kth layer can now be determined by plane stress constitutive equation where Q k is the transformed reduced stiffness of the kth layer in the z-location, which varies with the tissue orientation of each layer. Acting