Soft Soil Tire Model

If a terrain is assumed to be homogeneous, its pressure-sinkage relationship (the result of the plate-sinkage test) may take one of the forms which are shown in Figure 1, and it may be characterized by the following empirical equation proposed by Bekker. [1]

Where is pressure, is the width of the rectangular contact patch area, is sinkage, and , and , are pressure-sinkage related parameters.

There is a variety of criteria proposed for the failure of soils. One of the most widely used is the Mohr-Coulomb criterion which states that the maximum shear strength of the soil is:

Where is the apparent cohesion, is the normal stress, and is the angle of internal shearing resistance of the material. The aforementioned parameters can be derived with a shear test for different pressures, as shown in Figure 2:

In the case of a moving wheel, the contact patch area can be determined as a function of sinkage , the tire’s radius , and elastic sinkage using the following equations:

Using the pressure-sinkage relationship proposed by Bekker (Eq. (1)), the tire normal stress distribution can be calculated as a function of wheel angle [2]. Specifically,

where b is the tire’s width. Moreover, is the wheel angle at which the normal stress is maximized [2] and can be determined as:

In the above equation, and are parameters that depend on the wheel-soil interaction, is the longitudinal slip, and is tire entry angle.

Furthermore, Eq. (5) can be modified in order to account for the contribution of soil damping [3]. In that case, the normal stress distribution is given by:

Where is the soil damping, is soil’s compression rate/velocity, and is the contact patch area.

As already mentioned, the maximum shear strength of the soil can be calculated using Eq. (2) based on the Mohr-Coulomb failure criterion. In addition, the maximum shear strength between the tire and the soil can be approximated as a function of the pressure and the friction coefficient:

Consequently, the minimum shear strength (soil–tire and internal soil) is used for the shear stress calculation in order to account for the friction between the tire and the soil [4].

The friction coefficient normally has a value between 0.3 and 1.0 and largely depends on the material of the tire and the water content of the soil [4].

The shear stresses and are calculated using identical expressions [5], [6].

In the above equations, and represent the shear deformation modules, which are provided by the following equations:

Where is the slip angle of the tire. Moreover, the soil deformations and can be formulated as functions of the wheel angle [2], [7]:

Typical values for the parameters and can be found in the following publications:

1. ‘Terramechanics-based Model for Steering Maneuver of Planetary Exploration Rovers on Loose Soil’ by G. Ishigami and K. Yoshida, 2007 [5]
   • [m]
   • [m]
2. ‘Soil-tire interaction analysis for agricultural tractors: Modeling of traction performance and soil damage’ by A. Battiato, 2014 [8]
   • [m] for clay
   • [m] for clay loam
   • [m] for silty loam
   • [m] for loamy sand
3. ‘Terramechanics-based analysis for slope climbing capability of a lunar/planetary rover’ by K. Yoshida and G. Ishigami, 2004 [9]
   • [m] for dry sand
   • [m] for regolith simulant
4. ‘Analysis of off-road tire-soil interaction through analytical and finite element methods’ by H. Li, 2013 [3]
   • [m]
   • [m]

Moreover, according to Wong [10], based on experimental data collected, the value of ranges from 0.01 [m] for firm sandy terrain to 0.025 [m] for loose sand and is approximately 0.006 [m] for clay at maximum compaction. For undisturbed fresh snow, the value of varies in the range from 0.025 [m] to 0.05 [m].

Bulldozing resistance is developed when a soil mass is displaced by the wheel. Regarding lateral forces, due to tire’s sinkage, a bulldozing force which acts on the side of the wheel must be added to the shear force exerted on the contact patch due to the tangential stresses . In this tire model, the Hegedus resistance estimation method [11] is used in order to calculate the bulldozing force. As depicted in Figure 7, a bulldozing resistance is developed per unit width of a blade, as the blade moves toward the soil, given by the equation:

Where is the soil density. Moreover, the destructive angle , based on Bekker’s theory [12], can be approximated as:

In order to account for tire deformability, a larger substitute circle is used to describe the contact patch between the tire and the soil [1], [4], [13], as depicted in Fig. 8. In order to calculate the diameter of the substitute circle, an iterative procedure is followed until the soil vertical reaction force and the tire vertical force are balanced. The former is calculated from an integration of the normal and shear stresses in the contact patch area, while for the latter the vertical tire stiffness is used along with the tire’s deflection. More specifically,

The last three equations are solved iteratively until:

Where the tolerance is a function of the tire’s maximum vertical load. If the tire’s deformation is negligible, you can set the boolean flag RIGID_MODE= 'TRUE' in the tire’s parameter file in order to simplify and accelerate the calculation of the vertical force . Specifically, no iterations are required since only Eq. (18) is used in which the quantities and are replaced by and respectively. By default, the aforementioned boolean flag is set to FALSE, making use of the substitute circle concept and thus providing an increase accuracy.

For the multipass effect, the response of the soil to repetitive normal load needs to be established. Specifically, the mathematical description of the normal pressure distribution must be modified in cases of existing pre-compaction of the soil. As shown in Fig. 9, the normal pressure distribution will be comprised initially from an elastic part , which is equal to the elastic (unloading) sinkage that a previous tire has created. Then the pressure-sinkage relationship continues according to Eq. (1). Finally, an unloading elastic part is encountered.

As already mentioned, one part of the induced soil deformation is elastic (elastic sinkage), and the remaining part (plastic sinkage) is irreversible. The elastic part is provided by the equation:

Where is the soil elastic stiffness.

If the multipass effect is not of interest, you can set the boolean flag MULTIPASS='FALSE' in the road data file. By default, this flag is set to TRUE in order to enable updating the road data.

There are two available contact methods in the soft-soil tire model. The required outputs for both contact models are the soil elevation and the local inclinations of the road surface . The default option is to collect all the spring regions that belong to the contact patch area and use an averaging process to calculate a representative value for the soil elevation. This is accomplished through a Radial Basis Function (RBF) interpolation. Moreover, using this approach, the local inclinations of the road surface can be calculated in an efficient way [14].

The second option is to use the 3D Enveloping contact method [15], [16]. This contact model can provide improved accuracy in cases of short wavelengths (sharp steps) in road height. However, it is a more computationally expensive approach in comparison to the aforementioned RBF interpolation method. The core idea is to use a series of ellipses to scan the road profile and produce an effective road plane which is defined by three quantities. Namely, the modified effective height , the effective forward slope and the effective road camber angle . The simplest form of this contact model is depicted in Fig. 12, for which:

Where is the cam center height.

In general, two parallel tandems are not sufficient to obtain accurate results for sharp irregularities. Therefore, more parallel tandems as well as intermediate cams are added, as depicted in Fig. 13. Intermediate cams between the left and right sides are not necessary since their contribution cancel out when calculating the effective road camber angle [15]. The same equations as presented for the previous simple enveloping contact model apply again, but now extended for the case of multiple parallel tandems and intermediate cams. Specifically,

It should be emphasized that you can set which of the tires (if any) will use the 3D Enveloping contact method. This means that in a model both tires that use the Radial Basis Function (RBF) interpolation method as well as tires that use the 3D Enveloping contact model can co-exist leading in an optimal combination of accuracy and computational efficiency. This can be accomplished through different tire parameters files depending on the application’s needs.