Predicting the Passability of Wheeled Tractors

The variety of forwarder parameters should be considered when developing a new formula in a new iteration in the case of either machinery type. Consequently, the new model should be universal, without significant improvements in the mathematical structure of formulas, but with only minor clarifications.

The principal soil parameters used in the mathematical modeling of tractor passage are rigorously standardized and tested indicators:

• Modulus (of deformation – E, MPa and shear – G, MPa);
• Specific characteristics (coupling – C, kPa);
• Poisson's ratio - ν;
• Internal friction angle – φ, °;
• The thickness of deformable layer – H, m.

Such a number of features may cause difficulties in calculations. Moreover, they become more complex given the variability of the physical and mechanical properties of the soil. Due to the forest soil surface heterogeneity that tractors may encounter while operating, it is advisable to add such complications.

2.1 Physical and mechanical properties of bearing surfaces

The complexity of working with dynamic systems on bearing surfaces with variable parameters has been well documented in review [24]. Therefore, some developments allow systematizing the approaches used to define physical and mechanical characteristics of forest soils acting as bearing surfaces and create conditions for the wheeled tractor operation.

For the estimation of soil parameters, the following dependencies were obtained:

$D=\left(X_{x D}+Y_{x D} I_L\right) e^{X_{y D}+Y_{y D} I_L}$       (1)

where, X_xD, Y_xD, X_bD, Y_bD are numerical factors (Table 1), in which D summarizes the calculations for E, С or φ; I_L is the soil consistency index, e is the soil porosity index.

Table 1. Numeric values of the factors used to calculate the soil's physical and mechanical parameters

Values of deformation modulus allow determining the shear module:

$G=\frac{5 E}{2+2 v}$      (2)

Due to the lack of a forest soils classification, the version prepared in this study will be considered (Table 2).

Table 2. Example of forest soil classification by strength level

The influence of soil moisture content, composition, and temperature is derived indirectly through changing parameters in Table 2.

2.2 Experimental modeling of the wheel impact on the soil

Figure 1 shows the schematic version of indenting the conical portion of the penetrometer into the semi-space part of the soil.

1.png

Figure 1. Cone penetrometer operation diagram

A numerical representation of this process can be expressed as follows:

$C I=-C \cot \phi+\Theta \cdot \frac{24 G^m(\tan \alpha+\tan \phi)(1+\sin \phi) \tan \alpha}{d^2 \gamma^2(m-2)(m-3)(3-\sin \phi) \tan ^3 \phi}$                 (3)

where, Θ and m are factors taking into account the angle of internal friction.

$\Theta=\{C+(Z+L) \cdot \gamma \operatorname{Tan} \varphi\}^{3-m}-\{C+Z.$$\gamma \operatorname{Tan} \varphi\}^{2-m} \cdot\{C+(Z+3 L-L m) \cdot \gamma \tan \varphi\}$      (4)

$m=\frac{4 \sin \phi}{3(1+\sin \phi)}$     (5)

where, d is the width of the penetrometer’s cylindrical part, L is the length of the penetrometer's cone part, Z is the depth of the penetrometer's cylindrical part, α is the angle of the cone's sharpness.

A series of random values such as soil porosity and soil consistency factors were generated to account for soil parameters and cone index. The E and CI values were then compared when the deformation modulus ranged from 0...3 MPa and 0...2 MPa for forest and wetland soils, respectively, in increments of 250 kPa.

The above equation is verified for partial cases - different forest and waterlogged soils, clays, and loams.

In this study, a penetrometer with the following readings was used:

• L=31 mm;
• d=35.7 mm;
• α=30°.

2.3 Study of soil parameters affecting recovery

This experiment used 100 soil samples from different forest zones in Arkhangelsk. Moisture content was measured by the thermostatically controlled weight method. To this end, a 20 cm long drill bit with a diameter of 4 cm was employed. Samples were placed in aluminum beakers and dried at 105℃. Soil temperature at this depth was measured using a TPV-50 mercury thermometer.

The time of forest soil recovery was analyzed based on some of its characteristics. The results showed that soil moisture content ranges within 30% to 60% and temperature within –5℃ to 25℃. Also, the following properties of soils were measured: Deformation modulus, shear modulus, internal friction angle, and specific coupling. It was established that the time of soil recovery depends on its temperature (T, ℃), moisture content (W, %) and deformation modulus (E, MPa) with a high correlation ratio (R²=0.8660):

$t_r=0.039 W^{1.42} T^{-0.27} E^{-0.21}$     (6)

Hence, the recovery time of the soil increases with the increasing percentage of moisture content and decreasing deformation modulus and temperature. The interaction between the forwarder and the forest soil was modeled using the above formula.

2.4 Modeling the forwarder indentation into the bearing surface

A schematic representation of this process, taking the soil as a half-space with a homogeneous texture, is provided in Figure 2.

2.png

Figure 2. Experimental diagram of the indentation process

Based on the above, the modulus of soil deformation can be determined:

$E=\frac{\sigma}{\varepsilon}$     (7)

where, σ is the regular soil condition and ε is the relative soil deformation.

Soil stresses are distributed according to the law of attenuation:

$\sigma=\frac{J p}{1+(z / a b)^2}$      (8)

where, J is the forwarder indentor factor relative to its size; p is the pressure exerted on the soil surface by indentor; z is the counted vertical coordinate; a is the criteria of deformed layer thickness; b is the indentor width.

Considering the initial height of the elementary soil layer as z₀, its height in compression is:

$d z=(1-\varepsilon) d z_0$     (9)

During compression:

$d h_L=\varepsilon d z_0$      (10)

Based on these equations, compression can be calculated as follows:

$d h_L=\frac{\sigma}{E-\sigma} d z$     (11)

Then overall deformation is calculated by the following formula:

$h_L=\int_0^{H-h_L} d h_L$      (12)

The deformation modulus varies for each soil level according to a linear law, which will be considered as the first iteration:

$E=a_0+a_1 z$      (13)

where, a₀, a₁ are linear function coefficients:

$\left\{\begin{array}{c}a_0=E_1 \\ a_1=\frac{E_2-E_1}{H}\end{array}\right.$       (14)

where, E₁and E₂ are the deformation modulus values at the soil surface and depth corresponding to the thickness of the deformed layer (H), respectively. All three variables are distributed according to the law of equal density in the range from 0.5 to 5 MPa for modulus and from 0.3 to 0.8 m for H.

Given the previous formula, the total deformation can be calculated as follows:

$h_L=\int_0^{H-h_L} \frac{(a b)^2 J p}{(a b)^2\left(a_0+a_1 z-J p\right)+a_1 z^3+a_0 z^2} d z$      (15)

This integral was solved using a sample of random values of deformation modulus for the soil surface at a certain depth corresponding to the thickness of the deformed layer and the thickness of this layer itself. The values ranged from 500 kPa to 5 MPa (for modulus) and 30 - 80 cm (for thickness).

Given that the modulus of deformation can vary according to a quadratic law, the following formula shall be considered:

$E=b_0+b_1 z+b_2 z^2$      (16)

where, b₀, b₁, and b₂ are coefficients of the quadratic function:

$\left\{\begin{array}{l}b_0=E_1 \\ b_1=-\frac{3 E_1-4 E_2+E_3}{H} \\ b_2=\frac{2 \cdot\left(E_1-2 E_2+E_3\right)}{H^2}\end{array}\right.$      (17)

Then the formula for calculating the total deformation for this case is as follows:

$h_L=\int_0^{H-h_L} \frac{(a b)^2 J p}{z^2\left(b_0+b_1 z+b_2 z^2\right)+(a b)^2\left(b_0+b_1 z+b_2 z^2-J p\right)} d z$       (18)

In addition, the possibility of changing the deformation modulus according to the cubic law was considered:

$E=c_0+c_1 z+c_2 z^2+c_3 z^3$      (19)

where, с₀, с₁, с₂, с₃ are coefficients of the cubic function:

$\left\{\begin{array}{l}C_0=E_1 \\ c_1=-\frac{1}{2} \cdot \frac{11 E_1-18 E_2+9 E_3-2 E_4}{2 H} \\ c_2=\frac{9}{2} \cdot \frac{2 E_1-5 E_2+4 E_3-E_4}{H^2} \\ c_3=-\frac{9}{2} \cdot \frac{E_1-3 E_2+3 E_3-E_4}{H^3}\end{array}\right.$       (20)

Then the total deformation for such variant is:

$h_L=\int_0^{H-h_L} \frac{a^2 b^2 J p}{z^5 c_3+z_4 c_2+z^3\left(c_3 a^2 b^2+c_1\right)+z^2\left(c_2 a^2 b^2+c_0\right)+z c_1 a^2 b^2+a^2 b^2\left(c_0-J p\right)} d z$      (21)

Pressure at compression is calculated as follows:

$p=\frac{G_W k_d}{b l k_f}$      (22)

where, l is the length of forwarder indentor; k_d and k_f are coefficients of load dynamics and the shape of the contact pattern, respectively:

$k_f=0.8949 E^{-0.12}$      (23)

$k_d=\frac{l}{l+v t_{\mathrm{P}}}$        (24)

where, v is tractor driving speed.

Given the relationship between the compressive strain and the total deformation, we can find the total indentation:

$h=h_L \frac{p_S}{p_S-p}$       (25)

where, p_S is an index of pressure that the soil can withstand without deformation:

$p_S=p_{S 0} \alpha_Z$       (26)

where, p_S0 is the index of pressure that a soil layer with unlimited thickness can withstand without deformation, which depends on several factors, including the width and height of the indentor; α_Z is the coefficient that takes into account the thickness of the deformable soil layer.

The value p_S0 can be calculated as:

$p_{S 0}=0.5 K_1 B_1 N_1 \gamma b+N_2 \gamma h+K_3 B_3 N_3 C$       (27)

where, K_iis the factor that considers the shape of the contact pattern; N_i and S are coefficients that depend on the internal friction angle of the soil, and B_i is the factor that considers the angle of load application.

Indentor’s geometric parameters allow calculating coefficients K_i.

$K_1=\frac{l}{l+0.4 b}$       (28)

$K_3=\frac{l+b}{l+0.5 b}$        (29)

The values of N_i and S coefficients are calculated relative to φ:

$S=\tan \left(\frac{\pi}{4}-\frac{\phi}{2}\right)$       (30)

$N_1=\frac{1-S^4}{S^5}$      (31)

$N_2=\frac{1}{S^2}$      (32)

$N_3=\frac{2\left(1+S^2\right)}{S^3}$       (33)

Indentor load angle can be calculated as:

$\beta=\arctan \left(\frac{\tau}{p}\right) \pm \alpha$       (34)

where, α is the angle of rut slope and τ is the shear stress:

$\tau=\frac{j G}{t} \cdot \frac{t p \tan \phi+C t-C j}{j G+t p \tan \phi+C t-C j}$       (35)

where, t is the grouser pitch and j is the shear strain, which can be calculated by the formula obtained earlier:

$j=1.33 \sqrt{s^3 l}$       (36)

The value of the deformable layer thickness factor can be obtained as follows:

$\alpha_Z=1+\frac{H^* h}{2 H \cdot\left(H-h-0.25 H^*\right)}$       (37)

where, H^* is the factor that depends on the internal friction angle and indentor width:

$H^*=\frac{\sqrt{2}}{2} \exp \left[\left(\frac{\pi}{4}+\frac{3 \phi}{4}\right) \tan \frac{3 \phi}{4}\right] b \cos \frac{3 \phi}{4} \tan \phi$       (38)

2.5 Passability conditions

Thus, considering all the parameters of the soil and forwarder, the conditions ensuring the tractor passability in all situations can be established:

(1) The rut depth must be below a specific critical value:

$h \leq h_{\text {crit }}$      (39)

The remaining conditions relate to forwarder characteristics.

(2) The factor value of resistance to tractor movement is equal to the ratio of rolling resistive force F_R to the shear modulus:

$\mu_R=\frac{F_R}{G_W}$      (40)

In turn, the force depends on the indentor width and the pressure on the soil at a certain depth:

$F_R=\int_0^h b p d h$       (41)

Thus, the second condition is as follows:

$\mu_R M \leq N v$      (42)

(3) Sufficient level of traction on the soil should be ensured, i.e., the value of the tractive force, in this case, is non-negative:

$\mu_P \geq 0$      (43)

The force itself can be calculated as:

$\mu_P=\mu_T-\mu_R$      (44)

where, μ_T is a coupling factor:

$\mu_T=\frac{R_X}{G_W}$      (45)

(4) The forwarder must be capable of traversing an elementary unit of roughness with a certain height, which is calculated based on the tractive force value:

$h_{\Pi}=\frac{D}{2} \cdot\left(1-\frac{1}{\sqrt{1+\mu_p^2}}\right)$       (46)

Given that the thickness of the deformable layer during passage is limited, the relative compaction of the soil can be determined:

$\xi=\frac{\gamma_P}{\gamma}=1+\frac{h_L}{H}$       (47)

where, γ_P is the specific weight of soil after tractor passage.

If uneven compaction of the soil shall be considered, the equation is reduced to the following form:

$d \xi=1+\frac{1}{H} d h_L$       (48)

2.6 Software implementation of models

All model codes were written and compiled in the Maple 2015 software language. The first step implies entering the known characteristics (γ, W, I_L, CI, α, e, T, H) of the soil and calculating the missing ones (E, C, φ, G) based on the CI value. Also, knowing the E value, CI, C, φ, G can be calculated.

The equivalent deformation modulus was used to calculate inhomogeneous surfaces for three cases of dependence on depth (linear, quadratic, and cubic) using two, three, and four values of given factors, respectively.

The parameters of different soil types (sandy, clayey, and loamy soils) were calculated according to the established formulas for calculating E, C, φ, G, CI.

Besides, the written code considers the characteristics of the forwarder: D, B, S, H_T, P_W, G_W, t, k_m, v, which are compared with critical values.

All calculations are repeated as many times as necessary to obtain optimal results. 

The value ranges of parameters used to run the code and generate indicators are shown in Table 3.

Table 3. Variances of working model parameters

The process of cyclic soil loading was also modeled by considering the data obtained when generating values in each subsequent iteration and generating new values.

To compare the model's performance in this study, the WES model was applied.

2.7 Result accuracy assessment

For defining the level of reliability of the data obtained, an approximation was made using the least-squares method. In fact, the accuracy is indicated by the value of R², which must be equal to 1 in the ideal condition of the system.

A mathematical model for predicting tractors' passability was developed during the study. Its efficacy is demonstrated when used in the field and relative to another model.

Modulus of deformation and cone index are closely related linear function indicators. They can be interchangeable when calculating other soil parameters, such as specific coupling, angle of internal friction, and shear modulus.

Observations have shown that in sandy, loamy and clayey soils, compaction characteristics change linearly in soil surface homogeneity. For heterogeneous soils, the value of the equivalent deformation module may be used, thus simplifying the calculations. Also, these soils are characterized by additional dependencies of the deformation module during the deepening of compaction: square and cubic.

Study findings have shown that drier and warmer soils are much faster to deform, as the modulus of their deformation is relatively high.

The developed model is consistent with the existing model and takes account of the following soil parameters: the type, modulus of deformation and shear, slope angle, specific weight, consistency, porosity, temperature, moisture, the thickness of the deformation layer; and the forwarder parameters: Width and height of the wheel, tractor driving speed, tire tread expression and internal pressure in them.

The model also predicts the extent of soil compaction and rutting depth after multiple tractors passages.

Therefore, the proposed model can be used to analyze the extent of the impact of forestry and agricultural tractors on soils.