MODELING THE EFFECT OF WHEELED TRACTORS AND SKIDDED TIMBER BUNCHES ON FOREST SOIL COMPACTION

Igor Grigorev1*, Olga Kunickaya1, Albert Burgonutdinov2, Evgeniy Tikhonov3, Valentin Makuev4, Sergey Egipko5, Edward Hertz6, Maksim Zorin7 1Yakut State Agricultural Academy, Department of Technology and Equipment of Forest Complex, Yakutsk, Russia 2Perm National Research Polytechnic University, Department of Roads and Bridges Building, Perm, Russia 3Petrozavodsk State University, Department of Transport and Technological Machinery and Equipment, Petrozavodsk, Russia 4Moscow State Technical University N.E. Bauman (Mytishchi branch), Department of Technological and Equipment LT 7, Mytishchi, Russia 5Don State Agrarian University, Department of Machines Environmental Engineering, Novocherkassk, Russia 6Ural State Forest Engineering University, Department of Technologies and Equipment of Timber Industry, Yekaterinburg, Russia 7Voronezh State University of Forestry and Technologies named after G.F. Morozov, Department of Industrial Transport, Construction and Geodesy, Voronezh, Russia


INTRODUCTION
Forests are considered not only a source of a timber and alternative energy source, but also ensure the environmental security of our planet. Forest is one of the main components of the human environment promoting biodiversity, which acquires undisputable importance for the full enjoyment of human rights [1]. Therefore, forest management, conservation, and regeneration methods must meet contemporary social, environmental, and economic human needs. Therefore, in improving the efficiency of timber harvesting and reforestation processes, such multiple factors as physical and mechanical properties of the soil, parameters of tires, the vehicle body, loading capacity, etc., should be considered [2]. Soil compaction is considered the main harm caused by timber skidding [3]. Besides, value of the correlation between economic value added and environmental harm is of high importance as well [4]. The integral indicator of soil compaction at a felling area is one of the most important factors for estimating the environmental efficiency of harvesting processes [4,5].
As the tire travels, it applies uneven pressure to the soil within the contact area [6]. This leads to the deterioration of the basic physical and physical-mechanical properties of the subsoil and topsoil layers, reducing soil fertility and increasing energy consumption. Numerous experimental studies present wheel and track designs for minimization of pressure on the soil [7][8][9]. A number of mathematical models of soil pressure distribution caused by tire have been proposed, which are based on the finite element model of soil compaction [10], smoothed-particle hydrodynamics [11], etc. Another important aspect in studying the environmental impacts of forest machinery on the soil is the interaction of skidding system with permafrost soils of cryolithic zone forests, which represent particularly susceptible forest ecosystems [12,13], and, according to [14], significantly subjected to commercial exploitation. Skidding process is the most environmentally harmful operation in wood harvesting [15,16], which may cause significant degradation of the soil at the interaction of the tractor`s skidding and the dragging timber bunch [17]. T problem of skidding systems and soil interaction. Also, the impact of forest machinery maneuvers on the change of soil compaction in side strips of the trail when transporting a bunch of tree-length logs is studied insufficiently. Оne of the main tasks for assessing the impact of skidding methods on the soil is to define the features of soil deformation and compaction under the action of individual elements of the skidding system like forwarder, limbs, butts, and tops [18].
Since new methods of more gentle impact on the soil have not yet been found, establishing such methods are relevant and in high demand to date. The main aim of this work was to develop a mathematical model for estimating processes of soil deformation and compaction under tires and certain parts of the skidding system. The method of mathematical modeling includes the description of physical processes and the conclusion of formulas for determining the parameters of compaction of different types of soil and the size of the shelterbelt under the influence of different components of the system. The obtained theoretical modeling data will be compared with the results of experimental measurements. It is assumed that the methodology developed in this study will allow to perform qualitative and quantitative assessment of technological impact on the soil during the projecting of maps for logging operations.

Methods of modeling the impact of tractor wheels on the soil compaction
Rheological models show a relationship between soil stress σ and relative soil deformation ε. A general linear equation for elastic-viscous-plastic soils takes the following form [19]: dε dp p p Eε μ t dt dt A nonlinear model of rut formation by tractor tires on soils with elastic-viscous rheology, which takes into account the change of properties with the track depth, can be expressed as follows: where: p -soil pressure; c, с х , μ o , μ, р n -variables characterizing track formation; t -time, x -compaction depth. When studying the performance of wheeled forest machines, one should look at the elastic properties of tires. Let us assume that the deformation of the tire at each point is proportional to the pressure it exerts on the soil: = ⋅ p β y where: β -tire elasticity, у -tire deformation. Figure 1 illustrates the process of rut formation by a wheel with an elastic tire after the (n-1)-th pass. This process can be described with the following expression: where: δx -rut formation after the (n-1)-th pass, h n -1rut depth after the (n-1)-th pass. The pressure acting on the soil in equal to: Figure 1: Rut formation by a wheel with an elastic tire after the (n-1)-th pass: ξ, а1, а -distance between the contact point and the wheel axis (measurements taken from different sides of the wheel); ω -angular speed of the wheel, R -wheel radius Here, the integration, according to Figure 1, is performed from -а 1 to а. The approximation is thus а 1 ≈а. After integration, we obtain the following formula: where the axis of the wheel shifts down by a certain value, which can be found using the following equation: The equation for wheel-induced rutting in elastic-viscous soils has the form Decrease in the rut depth is associated with the magnitude of the maximum elastic tire deformation. The depth of the rut after elastic tires is determined by applying the principle of superposition for the sequential solution of equations for a rigid rim under the initial condition h o =0 with respect to tire elasticity, as shown below: The relationship between the time of rutting in elastic-viscous soils, the rut depth, and elastic tire deformation can be expressed as follows: At the beginning of the soil compaction process, the viscous component therefore, dt = p -1 µ 0 dx and the rut depth is equal to: If μ о =0 and the initial conditions of rutting correspond to The expression (12) suggests that the equation for rut formation by elastic tires can be written as follows: From equation (14), the initial phase of rutting depends on both the elastic properties of the tire and the viscous properties of the soil. The asymptotic value of the rut depth can be found from the condition t → ∞, which corresponds to the maximum rut depth: It can be seen that the maximum depth of the rut under elastic tires linearly depends on the elasticity of the tire. Hence, the maximum sinkage depth of the wheel with the elastic tire in elastic-viscous soils can be determined by the following formula: whilst the first rut depth can be found using the formula and the depth of the following ruts is to be found with the help of the equation The presence of elastic deformability of the tire can be considered a degree of freedom manifestation, which, to a certain extent, makes it possible to control the performance of wheeled forest machines.

Methods of modeling the impact of skidding on the soil compaction
Most harmful effect on the soil during skidding is caused by turning and maneuvering operations while loading and transportation of timber bunches, which leads to additional shearing stress in a direction perpendicular to that of a normal forwarder's load [19]. Same great interest was aroused to studying the impact of a bunch weight on the soil which can cause an additional compaction during skidding [20]. Total forces acting while butt length immersion can be presented as follows: where Q is the gravity force of tree-length, P is the weight of tree-length acting on the tractor and amounting to 0.3Q,G is the weight of tree-length acting on the soil when it is skidded at the top and amounting to 0.7Q. At this, some part of a tree-length log can go into the soil to a certain value of primary immersion area h о . Тhe value hо can be determined as: where d is the trunk diameter, L is the length of a treelength log, the coefficient β≥1 is the relation of L to the maximum possible hanging part of the tree-length log, h к is the distance from top point of the tree-length log support on the bolster to the supporting surface. A and n are the parameters of soil and press tool of weight G with contact area d in exponential functions q=Ah n , i.e., in dependencies of press tool pressure q on the depth of its immersion h within the range from 0 to the depth H of the deformation area ( Figure 2).  The initial uniform pressure averaged over the area operates The process of soil deformation outside the contact zone with depth hо and radius а will be considered in a Cartesian coordinate system. As seen on Figure 2, the stress σ is determined as follows: Tangential stress is determined as follows: where u is the positive root of quadric equation; v is Poisson's coefficient According to the scheme presented in Figure 2, the ratio can be simplified to σ у = σ r . Then, from Eq. (23)- (25) follows that the stress state of the soil depends on the coordinates of the calculation contact point and the radius of the contact area а. The combined action of the vertical force G and tangential force F T of the tractor's traction initiate the force of resistance to tree-length log sliding F n , which results in in emergence of horizontal stress τc is formed perpendicular to the action of vertical stress σ z : = + c z τ σ tgφ C At the moment when the tree-length log deviates from the given traffic direction to the angle θ, the components of the stress tensor in the cylindrical system of coordinates zrθ can be determined as: The resultant value of tangential stress τ=τ с +τ rθ characterizes the total shear stress and should meet the following condition: If the given condition is fulfilled, the limit value of the contact approach or the maximum possible depth of the tree-length immersion area h s corresponds to the value Z on Figure 2. The dependence of τ on h̅ with a high determination coefficient (R 2 >0.93) and considering the angle parameter θ for a single tree-length can be expressed as follows:  29) and (30) show that the depth limit of the immersion area h s can be defined as: The value of a relative soil compaction within the trail is defined as: where ρ is the achieved soil density and ε is the deformation of compression. The radial stress q r along the radius r≥a calculated from the boundary of a contact area is described by the ellipse equation: The maximum radial distance r s , at which value q r will exceed the soil bearing capacity q s , can be interpreted as the size of the maximum compaction area in the radial direction from the track or as the width of a shelterbelt of undergrowth root system or cut-to-length wood:

Experimental procedure
For estimation of initial parameters of contact and other parameters of soil compaction, the basic values and constants in the course of modeling were used, which are presented in Table 1.
The soil densification experiments on the side strands of a trail were performed under laboratory conditions using a metrologically verified DOS-3-I electronic compression dynamometer, which includes 101WN strain sensor and indicator station R320. The analysis results were processed by software and Excel application [21,22]. For the tests, cylindrical wooden blanks were used, fabri-   Table 2. The sample of a soil was placed into the tray so that the distance from the soil surface to the bottom of the tray was 0.16, 0.08, and 0.06 m for the first, second, and third categories of soil, respectively. A press tool was attached to a circular plate of 0.4 m in diameter and manually im-  Note: ρ o is the initial density of soil samples and p is the density value of compacted soil samples mersed into the soil placed in the tray. At that, the plate with the attached press tool tightly fitted to the soil surface. Afterward, the press tool and the plate were rotated to the specified angle θ, wherein the measuring equipment of the sensor recorded the value of force applied to the indicator part of the strain sensor from the side of the soil under deformation. The main independent parameters to be monitored during the research were soil density used in the tests and consistency index. As variables were applied the die rotation angle θ and the distance from the die to the adjoining soil layer, for which the compressing force was measured at rotating the die. The output parameter was the maximum compressive force P occurring at the side soil layer adjacent to the die at its rotation. Based on the analysis, a three-level full-factor experiment was performed.
In order to test the model of rutting, three different categories of soil were sampled and then compacted using a standard compaction instrument. The plasticity limits of the soil were determined using the standard balancing cone of A.M. Vasilyeva. The soil elasticity and deformation moduli were determined using a tabletop lever press. The following data were found during the tests: a) diagrams of vertical soil stress vs relative soil deformation; b) the ultimate stress at failure of the soil sample and the corresponding relative deformation limits. Next, the influence of the cyclical vertical loads on the soil deformation process was evaluated. A load equal to 70-80% of the ultimate load was taken as a constant (fixed) one. Because the compacted soil samples had 1.8 to 2 times higher density than before the application of load (Table 2), the following assumptions were made: (1) this increase in soil density is associated with the application of the cyclic load; and (2) the true value of deformation is a combination of elastic deformation and residual (viscous-plastic) deformation.

RESULTS AND DISCUSSION
The results from the analysis of experimental data were converted to a graphical form to show the relationship between soil compaction at different moisture levels and variables such as the load acting on the soil specimens and load time (Figures 3 and 4). With these data, it is possible to determine the dependence of rut depth and runing speed (v) at a given load (P). Let us take the moving wheeled tractor as an example. For this, it is necessary to determine what kind of deformations occur due to the tractor wheel passes. To begin with, let us determine the duration of wheel-soil interaction. Even with a rut depth equal to the radius of the tractor wheel, the length of the contact ares will be less than 1.5 m. The tractor travels at the speed of 1-20 km/h. Therefore, the wheel-soil interaction will last 7 ÷ 0.45 sec. The relationship between rut depth and load time is depicted in Figure 3. As it can be seen in the Figure 3 below, the viscous component is associated with the load time of 10 sec. Therefore, the mathematical model may be simplified by considetring the viscous deformation alone. The model of rutting in elastic soils derived from the general equation for soil deformation will take the form . The pressure acting on the soil will be In this case, the rut depth after the first wheel pass (h1) will be = 1 0 P hV bμ , and, conse- . The graph will show a hyperbolic relationship between rut depth and wheel speed for the n-th wheel pass . The plastic soil deformation will be The calculation of total shear stress τ and stress tensor components showed that the soil deformation through shearing does not occur at skidding of a single treelength log over the given soil without any maneuvers of a skidding system outside the primary immersion area  Figure 5. As can be seen, the value ρ̅ varies from 1.28 to 1.67, i.e., by 30.5%, the respective values of r s differ by more than 2.5 times. This indicates the need to consider detailed characteristics of the soil in designing the skidding routes in order to determine the size of undergrowth shelterbelts. Mathematical modelling allowed identifying the dependencies of ρ̅ and r s final values on the skidding parameters for all soil categories. Figure 6 presents the calculation results of values ρ̅ and r s at N x changes from 1 to 10 and other parameters being constant. As can be seen, the dependencies r s (Nx) have logarithmical nature of curves. The values ρ̅ and r s at the variation of other skidding parameters h к , θ, and N were estimated based on these data. The analysis showed that when h к is changed by 100%, the change of ρ̅ does Noteworthy is that the 20% and lower variation is acceptable condition of skidding. Besides, the calculation analysis revealed that in milder soils, the variability of θ and q s parameters determines wider limits of allowable variation ranges of ρ̅ and r s values. Thus, given mathematical model allows calculating parameters and indicators of the skidding system and allowable variation range for stabilization of indicators [23,24]. The results of this research allows recommending the distance between the skidding trails, i.e., the areas between the skidding trails on which the machines can effectively use the boom, to be at least 20 m. Similar conclusions are drawn by other researchers [18,25]. The 4-meter-wide skidding roads at a distance of 20 meters imply that at least 20% of the forest stand covered by forest machinery. This percentage may be reduced by laying access roads at a greater distance from each other. Trees outside the reach of a skidder can be cut down manually in the direction of the machine so that the trunk would be within the reach of harvesting arm. Moreover, machines should move without turning within one straight line as much as possible to prevent the soil wounding on the roadsides. It is important to consider the balance between access roads and wood production areas. Narrower distance T between the tracks may lead to a reduction in the volume and value of production due to the soil compaction in the wood stands. However, wider spacing would increase the cleaning costs. Frutig et al. [25] found that a wider interval (optimum: 30-50 m) is advisable in the long-term perspective, as losses due to slower growth and less valuable wood are of great importance in calculating the entire turnover period. This conclusion coincides with the results obtained in this study. With the shrub wood covers, the weight of the machine is distributed over a larger contact area than the landing area of machinery, thereby reducing the contact pressure on the soil. However, this method requires the treelength logs to be carried, not skidded, out of the forest so that the brushwood cover remains intact [26]. However, the question of whether these wood covers provide sufficient protection in all situations is disputable. Another argument to consider is the displacement of nutrients in forest plantations, since, all organic matters usually spread throughout the plantation are now concentrated on the skidding trails.

ACKNOWLEDGMENTS
The work was carried out within the confines of the scientific school "Advances in lumber industry and forestry".

CONCLUSIONS
The results of the study allow concluding that the developed mathematical model for estimating the processes of soil deformation and compaction under the influence of the wheeled tractor and the skidded part of the timber bunch is quite accurate and effective. The behavior of soils under forest machines is influenced by a large number of factors, ranging from the rheological properties of soil to the forest machine design. The elasticity of tires and the running speed of forest machines can help to control the performance of forest machines. This can be done by reducing the pressure exerted on the soil and increasing the number of skidder passes 1.5-2fold. Predicting the size of the shelterbelt at skidding a bunch of tree-length logs considering the results of the experiments, the calculation data differ from experimental by not more than 10%. Besides, the results of modelling showed that the deformation of the soil by the bunch of tree-length logs occurs due to shearing processes, wherein the initial vertical (compressing) stress exceeds the radial (stretching) stress by 30-40%. Performed analysis allowed determining correlation between the intensity of total shear stress reduction and the size of the soil compaction area. Estimating the dependency of the shelterbelt width on the number of tree-length logs in a bunch for three different soil categories (mild, medium, and solid) showed, that a limit values for logs amount to 4 till 6 units depending on the soil category, accordingly. The results of the study may be useful in predicting the behavior of soil under forest machines. Developed models can be applied to draw up a technological map for projecting a logging area based on soil and skidding system specifics.