Application of Temperature and Process Duration as a Method for Predicting the Mechanical Properties of Thermally Modified Timber

This study aims to investigate the influence of thermal modification (TM) on the physical and mechanical properties of wood. For this purpose, the experimental part focused on selected influential parameters, namely temperature, residence time, and density, while the four-point bending strength is obtained as the output parameter. The obtained experimental data are stochastically modeled and compared with the model created by genetic programming (GP). The classical mathematical analysis obtained treatment parameters in relation to the maximum bending strength (T = 187 °C, t = 125 min = 0.780 g/cm3) and compared with the results obtained by genetic algorithm (GA) (T = 208 °C, t = 122 min, and = 0.728 g/cm3). It is possible to obtain models that describe experimental results well with stochastic modeling and evolutionary algorithms.


Introduction
Wood is one of the oldest construction materials, and new wood-based materials continue to develop and are being successfully introduced into the engineering and construction marketplace [1]. Wood modification can be defined as improving wood properties by creating new materials. In the last twenty years, there has been a notable increase in the use of modified wood for a variety of applications [2]. Among various modification options, thermal modification (TM) is reported the most important in terms of traded volumes in Europe. Up to 500,000 m 3 of TM wood is placed on the market annually [3].
Despite its properties, wood demonstrates some limitations in the exterior environment use. The main disadvantages of wood, such as poor dimensional stability and biological degradation or deterioration, are mainly due to the nature of the main polymers of the wood cell wall, primarily due to the abundance of free hydroxyl groups (OH) [12][13][14].
The study [15,16] states that thermal modification improved dimensional stability and fungal attack resistance, while the TM process reduced the mechanical properties of wood. Thermal modification can be presented as partial pyrolysis in a reaction chamber with almost anoxic conditions. This process results in a changed chemical composition of the modified wood. Among the three key wood components, the hemicelluloses are being the most affected [17]. The first signs of hemicellulose modification appear at rather low temperatures. Decomposition starts at 140 °C and α-cellulose at 150 °C [18]. Lignin is more resilient to elevated temperatures [19]. One of the key results of thermal modification is the reduction of readily available hydroxyl groups, which affect the sorption properties of wood [20]. The equilibrium moisture content of TM wood is therefore considerably lower than that of non-modified wood when determined under the same climatic conditions [21,22]. It should be considered that the water exclusion efficacy is one of the key factors that determine the wood performance in above-ground conditions [23]. The effect of the thermal modification depends on the modification parameters. Shuang-Yan et al. stated that the chemical components are the main parameters affecting the mechanical properties of wood fibers [24]. The thermal modification of wood at different temperatures and durations decreased mechanical properties and increased dimensional stability and biological durability [25][26][27][28][29][30][31][32].
Physical and mechanical properties, modulus of rupture (MOR), modulus of elasticity (MOE), impact bending, compression strength, hardness in tangential and radial direction were determined in research [33,34].
The thermal modification process is influenced by various parameters, is including temperature, wood density, and residence time. Wood processing and testing technologies, which have been applied for many years, can modernize the application of knowledge through modeling, simulation, optimization, process theory, and computer technology [35][36][37]. The experimental results were obtained according to relevant standards for determining some physical and mechanical properties of timber EN 408 [38] and EN 789 [39], and using generally accepted scientific methods based on stochastic modeling and optimization [40][41][42][43][44]. Regarding genetic programming versus stochastic modeling, GP with high selection pressure finds a solution earlier, i.e., it converges faster without compromising the quality of the resulting solution, since faster convergence causes a higher probability of getting trapped in the local optimum [45]. On the other hand, if the selection pressure of the model is too low, computer time is spent on useless iterations because the convergence is too slow. Thus, this state has to be avoided with proper programming and specification of limit conditions [45][46][47].
This study attempts to predict optimal parameters of wood thermal treatment (modification temperature and duration) of wood products based on bending strength using stochastic modeling and genetic programming.

Description of the Resources Needed for the Experiment
Experimental research related to this work is based on the thermal modification of samples in the chamber for thermal modification, performed at Biotechnical Faculty in Ljubljana. Unmodified wood samples of European beech (Fagus sylvatica), Linden (Tilia sp.), and Silver fir (Abies alba) for bending process, with the average density of 0.675 g/cm 3 , 0.472 g/cm 3 , and 0.366 g/cm 3 , respectively. For modified samples, the average density for Beech was 0.655 g/cm 3 , for Linden 0.455 g/cm 3 , and 0.353 g/cm 3 for Fir. The test samples had a minimum length of 19 times the depth of the section. Ten samples without visible defects, full cross-section, free from knots and resin pockets, with the dimensions of 380 mm × 50 mm × 20 mm, were prepared for bending strength test. Before the TM process, all samples were dried to a constant weight, over 24 h in a drying chamber from the "Kambič" manufacturer, at a temperature of 103 °C. Drying was carried out to determine the oven-dry mass for further mass-loss calculations. After drying, the samples were cooled in a desiccator and weighed in an unmodified state. They were then transferred to a vacuum reactor, where the thermal modification process was performed according to a commercial process (Silvapro ® , Silvaprodukt, Ljubljana, Slovenia) [48]. After the specimens were placed in a vacuum-pressure reactor, about 95% vacuum was achieved in the chamber, while the absolute pressure was about 5 kPa. The parameters of the heat treatment process for testing the mechanical and physical properties of European beech (Fagus sylvatica), Linden (Tilia sp.), and silver fir (Abies alba) are following the parameters used in commercial processes [48]. The heating of the samples goes through five phases in which the samples were heated to maximum temperatures (170 °C, 180 °C, 195 °C, 210 °C, 220 °C) and treated with different maximum durations of the process (78, 120, 180, 240, 276 mins). In the next 12 h, the samples were cooled and weighed.
Measurement of four-point breaking force of unmodified and modified specimens was performed on a testing machine SHIMADZU type SIL-50KNAG at the Faculty of Technical Engineering in Bihać. The specimens were made and tested as recommended in the standard BAS EN 408 + A1 ( Figure 1) [49]. Bending force (Fb) testing was performed at four points of ten specimens for each experimental point. The specimens were always positioned so that they were subjected to maximum load. The test specimens, having a minimum length of approximately 19 times the depth of the section, were simply supported and symmetrically loaded in bending at two points throughout approximately 18 times the depth. Figure 1 shows the maximum force measured within the loading points.

Experimental Procedure and Setup
The experiment was conducted by using a central composition design with five levels of the three main independent parameters: process temperature (X1 = T), process duration (X2 = t), and density (X3 = ), as shown in Table 1 [35]. The overall number of experiments conducted of bending strength: where: number of variables, number of symmetrically set points, number of repetition points in the plan center.

Polynomial Form of Bending Strength Function
Based on the experiment, can present the functional relationship between the response of bending strength and the investigated independent parameters, by the following polynomial form of mathematical model: (2) where states that: -dependent variable corresponding to the physical or experimental value of bending strength , , -independent variables corresponding to the physical values of the processing regime , , , , , , , , -regression coefficient of the mathematical model.
The formalized description of , , changes to mathematical , , , which describes the influence of their interactions , , on the output value in addition to the influence on all three sizes , , individually. Table 2 lists the experimental design matrix or coded form of input parameters.
Regression coefficients are calculated as: The mathematical analysis provided the model with ten significant coefficients in coded form , , , , , , , , , . The dispersion is homogenous; thus, the experiment can continue. The mathematical analysis provided the model with ten significant coefficients in coded form: After the transformation of Equation (7), the bending strength model as a function of the process temperature (T), process time (t), and density (ρ), has the following physical form:   A multiple regression coefficient R = 0.97 was obtained, indicating a good interdependence of the mathematical model parameters. The mathematical model describes the four-point bending of thermally modified wood accurately and reliably within the space covered by the experiment.

Modelling by Genetic Programming
A genetic model of the four-point bending force is developed using GPdotNET software [45]. The output variable is the bending force of thermally modified wood ( Table 3). The population size is expected to correlate positively with the quality and rate of convergence of the solution. Mutation probability is a key parameter that allows GP new genetic material to be introduced into the population and obtain better solutions. Therefore, a significant complex effect of mutation on the solution's quality is expected. The number of generations should be positively correlated with the quality of the solution and negatively correlated with the program execution rate. Figure 2 shows a model of the bending force of thermally modified wood in the form of an expression tree. Table 3 lists the experimental results and the results calculated with the obtained GP model of bending force.  GP model of bending force obtained for a specific GP configuration depends on the following parameters: population size, mutation probability, number of generations, and window size. These parameters can take on a wide range of values, which can significantly affect GP solutions' performance and quality. The influence of window size and other parameters on the operation of the genetic algorithm can be found in the literature [45]. The following GP parameters were used during modeling: Based on the obtained results, the experiment and the GP model are presented in Figure 3 and Table 3. The presented diagram shows the quality distribution of the solution of the set of parameters of thermally modified wood for bending. The diagrams suggest that the algorithm consistently finds solutions close to the optimum for these sets of parameters. In other words, by changing the input parameters, it is known up to what final breaking force the thermally modified wood can withstand the load. The correlation coefficient between the model and the validation set is 0.99. The data acquired in this way show that the model obtained by the genetic programming method is reliable and can be used in further research.

Optimization of Bending Strength of Thermally Modified Wood by Classical Mathematical Analysis Method
The obtained mathematical model (8) (11).
The following condition is used for the determination of the optimum function : By solving the system of equations, we can obtain the coordinates of a stationary point: These are nonlinear equations and solutions ( , , ). Iterative methods obtain stationary points. If there is a minimum, maximum, bending, or saddle point in the stationary point, it is necessary to examine the sign of the determinants ∆1, ∆2, and ∆3, as follows: Maximum bending strength function is determined with the expressions (12)-(15), as follows:  Figure 3 shows the maximum bending force that depends on the input parameters. In contrast, the whole diagram represents the behavior of the maximum bending force and the relationship between the two input variables. The extreme point Fs max (T, t, ) shows the intensity of the change of the bending force depending on the individual change of each of the analyzed bending process parameters. Mathematical models of the maximum force are obtained and determined with two variables through the optimal point, with an orthogonal cross-section:  Orthogonal cross-section of the objective function Fs = Y (8) through the optimal point, enable mathematical models of the maximum bending force, determined by one physical variable:  Figure 4 shows the curves of the thermal modification parameters for the maximum four-point bending concerning the different input variables and optimal points. The application of stochastic modeling, based on the statistical method of defining a mathematical model, confirms the determined dependence of the maximum four-point bending force and the optimal parameters of the thermal modification process. The obtained mathematical model is adequate for the optimization of the thermal modification parameters if the maximum bending force is taken as the criterion of optimality.

Optimization of Bending Strength of Thermally Modified Wood by Genetic Algorithm
Genetic algorithms use stochastic search methods, meaning that multiple population goodness in a sequence can remain approximately the same over generations. After a while, a generation with superior goodness will undoubtedly emerge. This way of working with genetic algorithms can cause a problem in defining the objective function. A common practice in applying genetic algorithms is to complete the genetic algorithm after a certain number of generations (iterations). The best chromosomes in the population are then studied. If no satisfactory solution is obtained, the genetic algorithm is restarted. The mathematical model (10) of thermally modified four-point bending wood obtained by genetic programming was derived according to a genetic algorithm. A gene representation over real numbers was used for the genetic algorithm. The GA parameters used for optimization are as follows:  population size 500,  number of iterations 272,  probability of mutation 5%,  crossing probability 90%,  probability of reproduction 20%,  rank selection method. The maximum value of the mathematical model was calculated for T = 208 °C, t = 121 min, and = 0.728 g/cm 3 , where the maximum bending force of the thermally modified wood was F = 5950 N. The comparison of this result with the maximum value obtained by the regression analysis, where the maximum force is F = 6438 N, indicates that the forces are approximately equal.

Comparative Results of Bending Strength Experiment and Modelling
Genetic programming is a very efficient method of model making, and unlike classical regression analysis, there are no restrictions on the degree of polynomials. Based on the concept of computer programs, this method is an intuitive and simple way to build a model over experimental data that can then be used in future research in various wood tests as a starting point. The results showed that at maximum flexural fracture force, genetic programming can yield results that are closer to experimental data than classical regression analysis, that is, better results than stochastic modeling (classical mathematical analysis). Comparative results obtained by stochastic and genetic modeling of the maximum bending strength of thermally modified wood are presented, i.e., the reliability of using mathematical models without conducting expensive experimental research. Table 3 lists the comparative results of the bending experiment and the results obtained according to stochastic mathematical modeling and genetic modeling. Table 3 summarizes that the results of the maximum four-point bending force, obtained by the stochastic mathematical model and the genetic model, are almost indistinguishable, which further confirms the adequacy of the obtained models ( Figure 5). Figure 5 represents the results obtained by methods of mathematical modeling and GP are approximate values concerning the experimental results of the maximum fourpoint bending force. Not many models are developed that enable the optimization of wood properties based on density-independent of wood species. The majority of the current approaches are based on considering two parameters, modification temperature, and duration of the modification process [50]. These models are based on the mass loss of the timber after thermal modification and are predominately focused on durability. Up to our knowledge, this is the first approach to address the mechanical properties of wood with models based on the density of wood, modification temperature, and modification duration.

The Comparison of the Optimal Results
The obtained models (8) and (10) present a very good base for finding the optimal parameters to increase process efficiency, which is verified by the confirmation test. The presented optimization approaches provide the same optimal values for the TM parameters and give accurate results (as indicated by the confirmation test) with a slight deviation between each other (Table 4). The modeling provides a helpful tool for optimizing thermally modified wood properties. A considerable part of the chambers for thermal modification is somewhat of small capacity. This enables the producers to adopt the process based on the consumers' quality (density) and the expectation (needs). Thus, the respective solutions enable the higher added value of the thermally modified wood and better competitiveness on the market.

Conclusions
By mathematical modeling of the four-point bending process, have obtained a model in a coded form that has ten significant coefficients, which means that the relationship between (T), (t), and the wood ( ) plays a significant role in the adequacy of the mathematical model of bending strength.
The study shows that stochastic modeling and evolutionary algorithms can provide reliable results that are consistent with experimental data. By optimizing the obtained models of maximum bending force, using the method of classical mathematical analysis and evolutionary algorithm, optimal values for , , and were obtained. The optimal parameters of the thermal modification process by classical mathematical analysis and genetic algorithms have approximately the same value and represent an ideal treatment that can be applied to the wood.
There is the possibility of further research in finding a model for predicting the bending force of fracture which will be based on additional parameters such as the influence of humidity, or the behavior of the timbersʹ different conditions of exploitation, using GA in the input parameters optimization process.  Data Availability Statement: All data included in this study are available upon request by contact with the corresponding author.