The theoretical modeling of the dynamic compaction process of forest soil

1.   Introduction

Tree felling is a very common method in forestry. Mechanical systems make this process easy, more efficient, and safe. At the same time, they have adverse effects on the environment, such as causing air pollution from gas emissions and compacting the forest soil ^[1,2,3].

Forwarders and skidders with the wheeled or tracked chassis cause damage to the soil's physical properties ^[3]. For instance, vehicle passage leads to compaction in its upper and lower layers. Repeated passages aggravate this effect and result in soil compaction that persists for a long time, threatening soil productivity ^[4,5,6]. Soil compaction refers to the structural rearrangement of soil, which increases soil density and reduces pore space size ^[7]. It causes, among other things, the shift of the soil air and alteration in the mineral and microbial composition of the soil. By reducing soil porosity, decreasing air and water permeability, and restricting root expansion, soil compaction increases the risk of soil erosion ^[7,8]. Consequently, it slows or halts seed germination, plant growth, and the renewal of the forest ecosystem ^[6,9].

Factors affecting the degree of compaction include climatic conditions (such as temperature), soil parameters (such as texture, porosity, composition, water content), vehicle parameters (weight, cargo load, tire pressure, vibration, contact area), and more ^[8,9,10]. The compaction process is very complex, and it is almost impossible to assess how much damage the harvesting operations have done to the soil ^[11]. Therefore, it becomes popular to build mathematical models to evaluate the impact of skidding operations ^[12].

Modeling is a reliable technique that can be successfully assess human impact on nature ^[13]. Mathematical models of soil-wheel interaction help to create optimal designs of rubber wheels and tracks ^[14]. The process involves describing the load force applied to the soil surface, assessing the increase in the applied load, and quantifying deformation ^[15]. The resulting models can serve as a tool to investigate the propagation of stress through the soil profile caused by tire rolling ^[4].

The input data for models come from studies conducted under actual field conditions, because they will differ for different types of soils, loads and the type of transport used ^[13]. Researchers typically measure the degree of compaction achieved at various depths under various vehicles with varying soil and vehicle parameters (rolling speed, slope angle, etc.) ^[9]. Another focus is the effect of traffic frequency and angle of inclination on the soil's mechanical properties and nutrient content ^[7]. Their findings show that wheeled tractors damage the forest soil more than the tracked ones ^[3].

Mathematical models of the vibration response of the soil have multiple mathematical laws and principles behind them. One of the popular ones is the finite element method, because it allows one to note inelastic and nonlinear mechanical displacements in the soil, as well as to analyze the history of these reactions, which is a huge advantage ^[16]. Other techniques include, but are not limited to, the Lagrangian-Euler method ^[17], the discrete element method ^[18], the finite volume method ^[19], the Smoothed Particle Hydrodynamics model ^[17]. The discrete element method and the Smoothed Particle Hydrodynamics model are effective in modeling the process of interaction between the soil and structural elements of agricultural equipment, but they require large resources to work ^[18,20]. The Lagrangian-Euler method contains two phases that are responsible for different deformation domains: the Lagrange grid for small deformations and the Euler grid for large deformations, but it is difficult to apply it in the context of the three-dimensional problem ^[17]. Given that soil is a porous system, the finite volume method can be used to model compaction, where crack growth is considered separately from compaction ^[19].

The volume of timber harvest in Russia reaches 240 million m³ while showing a tendency to increase annually. By transporting the harvested wood, the timber industry contributes significantly to soil compaction ^[21]. In this regard, this paper aims to create a mathematical model to assess the mechanical response of soil to dynamic loading of a wheeled skidder, to optimize the process of determining soil compaction.

2.   Materials and methods

The study breaks down the tractor-bundle-soil system into two subsystems with the view to examine the dynamic soil compaction by a skidding system and the disturbing effect of tire surface irregularities. While the first subsystem signalizes a radial tire vibration, the second one indicates an elastic deformation of soil under dynamic loading. By analogy to models used in truck-focused research, this paper considers the first dynamic subsystem as an oscillating system resembling the oscillating behavior of a loaded tractor with a variable wheel configuration.

The masses and moments of inertia within the given system are fixed and constant. The equivalent scheme appears as a linear oscillating system with constant dynamic parameters. A schematic image of an oscillating system with a tractor having a wheel configuration of 6 x 6 is depicted in Figure 1.

Figure 2 displays a similar scheme for an 8-wheel drive tractor.

While the 6 x 6 tractor has an axle in the rear, the 8 x 8 unit has both front and rear tandems, balanced, with two wheels per axle. Unlike vehicles two axles connected by a balancing spring, the tandem bogie is rigidly connected to the frame. Due to the need to simplify the above schemes for oscillating systems, one can follow the principle of additivity to reduce the number of external forces (tire ground pressure) to one. The resulting system is similar to that with a 4 x 4 tractor.

The coupling scheme for the semi-suspended bundle of logs loaded onto the tractor with a pincer claw is shown in Figure 3, where: γ is a deflection angle of the suspension; βp is a deflection angle of the pack of wood.

With a chokerless skidding tractor, the top and butt parts of the logs are arranged in a bunk. In this case, the coupling scheme for the semi-suspended bundle of logs would be as the one depicted in Figure 4.

3.   Results

3.1.   Description of the subsystem 1

When examining the first subsystem, the focus was laid on the influence of vertical and angular vibrations of the half-frames connected with a rigid universal joint shaft on dynamic soil compaction. The vibrations occur in the longitudinal and transverse planes. Previous studies show that it is possible to consider vibrations of a tractor and a skidding system on wheels without giving the account to torsional vibrations of the transmission and engine. Consequently, vibrations of a loaded tractor can be considered as oscillations in a dynamic system with four degrees of freedom. The first degree of freedom is associated with the vertical vibrations of a tractor having the center of mass at point Z. The second degree of freedom is linked to longitudinal-angular vibrations; the turn of the loaded tractor's mass in reference to the center of masses is defined as φ. The third degree of freedom is associated with transverse-angular vibrations of the two half-frames connected by the hinge joint at point γ. The fourth degree of freedom relates to vibrations of the front half-frame relative to the rear half frame and vice versa.

The wheeled forestry tractors have suspension systems that are symmetrical with respect to the longitudinal axis. The longitudinal-angular vibrations of the suspended mass, symmetrical with respect to the longitudinal axis of the machine, can be considered independently from the lateral-angular vibrations and vise a versa.

It is customary to represent the supporting surface profile using the correlation function for the impact R (τi), which can be approximated by the formula:

where: D is the variance of surface irregularities; α and β denote the coefficients of approximation.

The lag time τ_i = t_i-t₁ is the time interval between the first wheel hitting the irregular surface and the i-th wheel hitting it, expressed as:

where: l_i denotes the distance between the axes and the center of mass of the vehicle moving with speed υ.

3.2.   Modeling results for an unloaded skidder

The motion of the system was considered using the inertial system of coordinates (Figures 5 and 6), where C_k is the universal drive stiffness.

The equation is derived by applying the D'Alembert's principle. The system of differential equations for the motion of the skidding system can be presented as follows:

where: relative angular deviations are expressed as ${a}_{1} = \frac{2n{\mu }_{t}}{М}$, ${a}_{2} = \frac{2n{C}_{t}}{\mu }$, ${b}_{1} = \frac{{\mu }_{t}}{{J}_{SS}}{\sum }_{i = 1}^{2n}{l}_{i}^{2}$, ${b}_{2} = \frac{{C}_{t}}{{J}_{SS}}{\sum }_{i = 1}^{2n}{l}_{i}^{2}$, ${С}_{1} = \frac{{\mu }_{t}}{{J}_{fh}}{\sum }_{i = 1}^{2n}{l}_{i}^{2}$, ${С}_{2} = \frac{{C}_{t}}{{J}_{fh}}{\sum }_{i = 1}^{2n}{l}_{i}^{2}$; ${d}_{1} = \frac{{\mu }_{t}}{{J}_{rh}}{\sum }_{i = 1}^{2n}{l}_{i}^{2}$; ${d}_{2} = \frac{{C}_{t}}{{J}_{rh}}{\sum }_{i = 1}^{2n}{l}_{i}^{2}$; and n is the number of axles in a multi-axle tractor.

When the system of differential equations (4) with a real variable t is transformed into Laplace space, the result is a system of algebraic equations with a complex variable s:

where: L is the Laplace transform; and s = σ+iω is the complex frequency.

By solving the system of equations (4), one can find the transfer functions of the support surface with respect to the suspension masses of the skidding system, as shown below:

where: ${C}_{2}^{*} = {C}_{2}+C(1-\Delta \gamma)$ and ${d}_{2}^{*} = {d}_{2}-C(1-\Delta \gamma)$ are coefficients accounting for the hinge joint parameters.

To determine the vibration frequency characteristics, s= iω was substituted (where ω is the vibration frequency and i is the imaginary unit) into the Eqs (5) to (8), as shown below:

The amplitude-frequency characteristics of vertical (A_z), longitudinal-angular (А_φ), and transverse vibrations for the front ${(А}_{{\gamma }_{1}})$ and rear ${(А}_{{\gamma }_{2}})$ half-frames are taken as modules (10) to (13). These are $\left|{W}_{z}\left(i\omega \right)\right|$, $\left|{W}_{\phi }\left(i\omega \right)\right|$, $\left|{W}_{{\gamma }_{1}}\left(i\omega \right)\right|$, and $\left|{W}_{{\gamma }_{2}}\left(i\omega \right)\right|$ :

By applying the Fourier transform of the correlation function for surface irregularities in the longitudinal and transverse planes, it is possible to obtain the spectral density, which describes the frequency composition of this function. Usually, the spectral density of the supporting surface unevenness is determined through the coefficients of approximation α and β, as shown below:

Equations (13) and (14) can help to determine the spectral frequencies of the vertical (${S}_{m}^{Z}\left(\omega \right)$), longitudinal-angular (${S}_{m}^{\phi }\left(\omega \right)$), and transverse (${S}_{m}^{{\gamma }_{1}}\left(\omega \right)$ and ${S}_{m}^{{\gamma }_{2}}\left(\omega \right)$) vibrations:

Based on Eq (15), one can determine the correlation functions as the primary sources of spectral planes:

If in Eq (16), τ = 0, then the correlation functions σ determine the variances and mean square deviations (MSD) of the corresponding vibrations:

The mean square deviations ${\sigma }_{Z}$, ${\sigma }_{\phi }$, ${\sigma }_{{\gamma }_{1}}$ and ${\sigma }_{{\gamma }_{2}}$ characterize the scattering of vibration amplitudes with respect to their grouping centers. The latter are mathematical expectations ${T}_{Z}$, ${T}_{\phi }$, ${т}_{{\gamma }_{1}}$ and ${т}_{{\gamma }_{2}}$, which correspond to the idle state of the system. The value of the dynamic coefficient k_d depends on the values of ${\sigma }_{Z}$, ${\sigma }_{\phi }$, ${\sigma }_{{\gamma }_{1}}$ and ${\sigma }_{{\gamma }_{2}}$. Because these are dimensional variables (${\sigma }_{Z}$ is a metric unit and other variables characterize angular dimensions), there is a need to convert them into dimensionless quantities characterizing the processes in an idle state. In other words, one has to align dimensional variables with mathematical expectations.

The dimensionless mean square deviations ${\stackrel{-}{\sigma }}_{Z}$, ${\stackrel{-}{\sigma }}_{\phi }$ , ${\stackrel{-}{\sigma }}_{{\gamma }_{1}}$ and ${\stackrel{-}{\sigma }}_{{\gamma }_{2}}$ characterize the degree of change within the oscillating system in the vertical, longitudinal-angular and transverse planes during the dynamic load as compared to the static one. The coefficient k_d is crucial for the complex assessment of the dynamic impact of vertical, longitudinal-angular and transverse-angular vibrations on the soil. The coefficient of dynamic amplification of the system k_d, is determined as the mean square deviation ${\stackrel{-}{\sigma }}_{П}$ :

3.3.   The first method of semi-suspended skidding

The scheme of the transverse-angular (γ) vibrations of the tractor with a semi-suspended load (β), which takes into account the inclinations caused by rolling over the uneven surface profile, is shown in Figure 3. The differential equations for oscillations would be:

where ${c}_{1} = \frac{{J}_{n}+G{\left(D+L\right)}^{2}}{9g}$; ${c}_{2} = \frac{D{G}_{T}\left(1+L\right)}{9g}$; ${c}_{3} = \frac{1}{2}{c}_{t}n{d}^{2}-{G}_{T}{H}_{k}-\frac{2}{3}GL$; ${d}_{1} = {c}_{2}$; ${d}_{2} = {J}_{nk}+\frac{Q{D}^{2}}{9g}$; and ${d}_{3} = \frac{2}{3}QD$. Here, J_n represents the moment of inertia of the system in the transverse plane; Q is the weight of the bundle; D is the distance between the suspension center point and the center point of the bundle; L is the distance between the roll center and the suspension center point; G_T is the weight of the tractor; Н_k is the distance between the center of the mass of the tractor and the base of the front half-frame; J_nk is the moment of inertia of the bundle relative to the axis passing through the center point of the suspension and through the center of the contact area between the second end of the bundle and the soil; f₁(t-τ_i) is the function of the impact taking into account the lag time τ_i.

Multiplying both sides of the Eqs (19) and (20) by ${e}^{-St}$ and integrating them from 0 to ∞ enables the use of Laplace transforms to transform the system of differential equations into a set of algebraic equations with respect to a new complex variable S:

The generalized coordinates γ (S) and β (S) after solving the system of Eqs (21) and (22) would be:

By dividing both sides of the Eqs (23) and (24) by F₁(S), one can obtain the transfer functions ${W}_{\gamma }\left(S\right) = \frac{\gamma \left(S\right)}{{F}_{1}\left(S\right)}$ and ${W}_{\beta }\left(S\right) = \frac{\beta \left(S\right)}{{F}_{1}\left(S\right)}$

The amplitude-frequency characteristics of the transverse vibrations of the tractor-bundle system were obtained with the Eq (25) taking $S = i\omega$, where ω is the frequency of forced vibration (expressed as 1/s)

The vertical and longitudinal-angular vibrations were considered by looking at a tractor of specific mass (M) to which a timber bundle of specific length (L_с) and mass (Q) was attached. While the end part of the semi-suspended bundle is flexible, the end part of the bundle that lies on the ground is considered fixed. Hence, tractor at point A with a mass of М_а = М+ $\frac{2}{3}Q$ experiences vertical (Z) and longitudinal-angular (φ) vibrations, described by two differential equations (in this case, the bending stiffness of the end part of bundle is neglected):

where: ρ_y is the radius of inertia of the mass М_а in the longitudinal plane; n is the number of axes.

The amplitude-frequency characteristics of the vertical (${W}_{Z}\left(i\omega \right)$), longitudinal-angular (${W}_{\phi }\left(i\omega \right)$) and transverse vibrations of the tractor (${W}_{\gamma }\left(i\omega \right)$) and the bundle (${W}_{\beta }\left(i\omega \right)$) were obtained by representing Eq (28) in the operator form and making the necessary Laplace and Fourier transformations with respect to Eqs (26) and (27):

According to the theory of the stationary random functions, spectral densities of the amplitudes of corresponding oscillations are determined by multiplying the squared values of Eqs (30) to (33) by the spectral density S(ω) of surface unevenness. The energy spectra of vibrations are (Table 1):

3.4.   The second method of semi-suspended skidding

The oscillatory scheme in Figure 4 differs from that for the first semi-suspended skidding method. Firstly, the bundle loaded on the tractor forms a single system with the tractor's half-frames. The roll angle of the bundle is β= 0; hence, the roll of the system at the angle γ is considered.

The differential equations for the vertical (Z) and longitudinal-angular (φ) vibrations are similar to Eq (28). The differences apply to transverse vibrations of the front (γ₁) and rear (γ₁) half frames:

where: ${\bar{c}}_{1} = \frac{{\mu }_{t}}{{J}_{t}\left[{\sum }_{i = 1}^{2n}{l}_{i}^{2}+\frac{1}{2}{d}^{2}\right]}$; ${\bar{c}}_{2} = \frac{{c}_{t}}{{J}_{fh}}\left[{\sum }_{i = 1}^{2n}{l}_{i}^{2}+\frac{1}{2}{\rho }_{xn}^{2}n-\frac{2G{H}_{n}}{3{c}_{t}}\right]$; ${\bar{d}}_{1} = \frac{{\mu }_{t}}{{J}_{rh}}\left[{\sum }_{i = 1}^{2n}{l}_{i}^{2}+\frac{1}{2}{d}^{2}\right]$; ${\bar{d}}_{2} = \frac{{c}_{t}}{{J}_{rh}}\left[{\sum }_{i = 1}^{2n}{l}_{i}^{2}+\frac{1}{2}{\rho }_{x3}^{2}n-\frac{{H}_{n}G{H}_{3}}{{c}_{t}}\right]$; ρ_xfh and ρ_xrh are inertia radii of the front and rear half frames, J_fh and J_rh are moments of inertia in the transverse plane; H_fh and H_rh represent the distance between the center of the mass of the bundle and the base of the half-frames.

After necessary transformations analogous to the semi-suspended skidding scheme we obtain transfer functions and frequency response characteristics of the vertical ${W}_{Z}\left(i\omega \right)$, ${W}_{Z}\left(i\omega \right)$ transverse oscillations of the front ${W}_{{\gamma }_{1}}\left(i\omega \right)$ and rear ${W}_{{\gamma }_{2}}\left(i\omega \right)$ half frames:

The energy spectrum of oscillations is determined by the squared moduli Eqs (36) to (39) and the spectral density S(ω) of unevenness by analogy with Eq (33).

3.5.   Description of the subsystem 2

The second subsystem refers to elastic soil deformations due to vibrations from a loaded skidder. However, the use of this methodology is limited by the following factors: the modulus of elasticity depends on the magnitude of the repeating load ^[22], soil properties during compaction ^[23], soil moisture ^[24], the hardness and texture of the soil, which in turn depends on its type ^[25]. These limitations can be partially accounted for by the acoustic hardness of the soil. The modeling is based on the Kelvin-Voigt (Gogenesmer-Prager model) model with the addition of one elastic element. The analysis of the existing elastic-viscous models (mechanical, mathematical and physical) of soil behavior revealed that the dynamic load q_d on the soil, caused by dynamic vibrations of the tractor, can be assessed through the transfer of acoustic stiffness from one medium to another.

The coefficient of load transfer К_λ is expressed as:

where: λ₁ is the acoustic stiffness of the protector's material; λ₂ is the acoustic hardness of the soil.

The acoustic stiffness of tires is determined by the formula:

where: ρ_t is the density of the tire material, kg/m³; υ is the velocity of the longitudinal compression wave, m/s.

The acoustic hardness of the soil is determined by the formula:

where: ρ is the soil density, kg/m³.

The study of the viscoelastic soil was carried out taking into account that deformations ε did not occur instantly, but were a function of time. The reconstruction of the soil structure is an internal and inter-crystalline process of particle rearrangement. To create a dynamic model of soil compaction, static and dynamic load diagrams were used.

Soil deformation ε_d corresponds to the dynamic compaction, and soil deformation ε_s corresponds to the static one. The compaction of the medium is defined as the sum of these two deformations: $\varepsilon = {\varepsilon }_{d}+{\varepsilon }_{s}$. The relative soil compaction $\bar{\rho } = {}^{\rho }\!\!\diagup\!\!{}_{{{\rho }_{0}}}\;$, where ρ₀ is the initial soil density, is related to the deformation ε by the following formula:

Given the linear relationship between the value of deformation and soil density, the compaction value can be taken as the sum of dynamic and static compaction values. In contrast, soil deformation during shock compression is defined just by the dynamic compaction curve, as no static deformation occurs. We accept, as a first approximation, the linear nature of the loading dependences of stresses during dynamic and static processes:

where: Е_d and Е_s- dynamic and static deformation moduli.

Dynamic soil compaction can be expressed as follows:

The Eq (45), however, does not take into account the impact of the load.

Knowing the initial contact parameters (tractor weight, material weight, wheel size, soil density and others) it is possible to determine the degree of influence of the system on the soil. This will make it possible to select such an amount of material so that the final impact on the forest floor and the resulting negative effects are minimal.

With each pass of the skidding system, the value of dynamic soil compaction will increase. The maximum compaction value depends on many factors. Among them are the tire surface irregularities, parameters of the skidding system, and the rheological properties of the soil.

Studies on soil compaction involve a variety of methods, Wang ^[26] analyzed soil compaction under multi-location tamping with the help of the finite element method. They calculated the effective stress rate using the elastic coefficient and the deformation rate, decomposed into elastic and plastic parts. The present work relies on the Kelvin-Voigt model, the transformation of which enabled to find the dynamic soil compaction as a function of the ratio of the load to the deformation modulus without taking into account the static component of deformation. Khitrov and Andronov ^[27] offered a mathematical model of the interaction between soil and timber trucks. They set the value of the compressive deformation at different depths. Compaction was taken as the sum of the productions of the coefficient ς (it depends on the deformation modulus E and the value of B) by the integral value of stress. Both studies did not take into account the effect of tractor vibration on the soil.

Some authors modeled the impact of a vehicle in vibrated state on the rheological parameters of the soil ^[28,29]. For this, they created a set of differential equations. The present study also applied a system of differential equations to model the soil behavior, but this system takes into account the number of axles in a multi-axle tractor. The said system was then transformed using the Laplace transform.

Edwin et al. ^[14] went for the systems of differential equations proposed by Solis and Longoria. They sought to model the impact of tracked vehicles on soil properties. The authors took into account the tangential forces acting on the track due to shear resistance of the soil and thus determined shear deformation. The present study used a similar approach to model the second subsystem, represented as a set of disturbing factors associated with the soil and wheel.

4.   Conclusions

Theoretical modeling of the soil compaction process under the influence of forest machines during their operation was carried out in this work.

To study breaks down the tractor-bundle-soil system into two subsystems to increase the accuracy of the mathematical model of dynamic soil compaction. While the first subsystem is associated with vibrations of the tractor in different planes and radial tire vibrations. The second one indicates an elastic deformation of soil under dynamic loading. The effect of tractor vibrations was depicted with the help of differential equations obtained using the D'Alembert principle, the method of Laplace transforms, and a system of equations for calculating vibrations. A modified Kelvin-Voigt method was used to model changes in the rheological parameters of the soil. The study developed models to determine vibration effects on the soil under an unloaded tractor and two skidding systems using different methods of semi-suspended skidding.

The results of the work can be applied to establish the degree of impact of forest machines with different levels of loading on the forest soil, which will make it possible to select the optimal weight of wood for transportation. As a result, this will reduce the load on the soil and reduce the negative impact on the root system of trees.

Future work in this direction is required to develop new mathematical models or improve the existing ones to create an optimal model that is versatile for different soil types and tractor systems. Future research may also try to find other ways to address the limitations, such as considering changes in soil moisture before and after compaction.

Acknowledgments

The work was carried out within the confines of the scientific school "Advances in lumber industry and forestry". The underlying content of this paper was supported by the grant of the Russian Science Foundation No 22-26-00009, https://rscf.ru/project/22-26-00009/.

Conflict of interest

The authors declare that they have no conflict of interests.