Comparative analysis of FEM tire-soft snow interaction and theoretical model based on Bekker’s coefficients

Tire-deformable terrain interaction is a complex phenomenon that has proven difficult to characterize accurately and comprehensively using numerical models. While researchers have attempted to develop various numerical models to estimate tire behavior on deformable terrain, the most accurate models to date rely on empirical data and are only applicable to specific rolling conditions. In this research paper, a comparison between theoretical mechanics and a finite element method (FEM) model of a rigid tire rolling on a layer of soft snow is presented. The goal is to compare the data obtained from the FEM tire-soft snow interaction analysis with the data generated by a theoretical numerical model. The FEM snow model is developed using the Drucker-Prager cap material model while the tire is assumed rigid. The assumption of the rigid tire can provide reasonably accurate results in the case of the soft snow interaction model. The theoretical pressure-sinkage coefficients of the Bekker equation (kc, kϕ and n) for soft snow are calculated by simulating vertical penetration tests in the FEM snow model. The results of the pressure distribution and rolling resistance due to compaction obtained by the FEM and the theoretical model results are then compared.


Introduction
The interaction between tires and soil is a highly intricate phenomenon, posing significant challenges in accurately and comprehensively characterizing it through current numerical models.Numerous researchers have made attempts to develop various numerical models to assess the behavior of tires rolling on soft soil.However, to date, the most precise models rely on empirical data and are limited in their usability to specific rolling conditions.
The study on the motion resistance of two truck tires of different sizes on three distinct soil types was carried out near Buenos Aires, Argentina, at an elevation of 22 meters above sea level, and featured a slope gradient of 0.5%.The soil at the location belonged to drainage class 4 [1].The "Shear Interface Imaging Analysis Tool," a novel technique for experimentation and analysis, was introduced by Moreland [2].This technique is specifically designed to explore the fundamental principles of Terramechanics in great detail, enabling the visualization and analysis of soil motion both at and below the interface between the wheel and soil [2].Pressure-sinkage experiments were performed for utilizing penetration plates of varying dimensions to characterize a model based on Bekker's formula [3].In order 1303 (2024) 012046 IOP Publishing doi:10.1088/1757-899X/1303/1/012046 2 to simulate the soil conditions, a mixture of bentonite and water was selected as a substitute for seabed sediments, based on experimental data collected from China's ocean poly-metallic mining contract area [3].
Our research motivation is represented by the idea of integration between the theoretical mechanics of tire-soil interaction with the experimental data obtained in the field for a set of operation conditions.The primary objective of this research is to compare data acquired from a finite element method (FEM) analysis of tire-soft snow interaction with the data generated by a theoretical numerical model.For the FEM snow model we employ the Drucker-Prager cap material model, while assuming the tire to be rigid.This assumption is expected to yield reasonably accurate results for the soft snow interaction model.The coefficients utilized in the FEM snow model will be sourced from our comprehensive literature review, primarily relying on experimental measurements [4].Under typical conditions, snow displays intricate properties characterized by non-linear viscoelastic behavior and significant deformations, encompassing both volumetric and deviatoric strains [5], [6].The literature presents various approaches for modeling the interaction between snow and tires, including analytical and semianalytical methods [7], finite element methods [8], and, more recently, particle-based methods [9].Finite element methods (FEM) employ diverse material models to represent snow, specifically chosen to accurately capture its highly compressible nature [10].To determine the theoretical pressure-sinkage coefficients (k c , k ϕ , and n) for soft snow soil, vertical penetration tests will be conducted, using our FEM model.By calculating these coefficients, our objective is to compare the results of theoretical models, which account for pressure distribution and rolling resistance due to compaction, with the outcomes obtained from the FEM model.Furthermore, we aim to evaluate the thrust capabilities as a function of wheel slip, using the shear stress distribution in the case of a rigid wheel from the FEM analysis.Through these comparisons, our aim is to incorporate relevant experimental data to further validate the accuracy and reliability of the FEM analysis and numerical modeling.

Research Methodology
This section presents the methodology for developing the FEM snow model and the Bekker-Wong coefficients.Furthermore, this section presents the analytical approach to estimate the sinkage and the motion resistance due to terrain (snow) compaction of a rigid wheel.

FEM snow model development
The snow model was implemented using the explicit version of the commercial software ABAQUS, incorporating a methodology influenced by the research of Shoop [10].The Arbitrary Lagrangian Eulerian (ALE) meshing technique was utilized to model the snow.To capture the behavior of soft snow, the Drucker-Prager cap model was employed to define its material properties.Table 1 provides the specific values utilized to comprehensively describe the properties of the snow model.The Drucker-Prager failure criterion is a model that considers the pressure-dependent behavior of materials and determines if they have failed or undergone plastic yielding.This material model is applicable to various pressure-dependent materials such as rock, concrete, polymers, foams, and compacted snow.It provides a three-dimensional, pressure-dependent model for calculating the stress state at which a material reaches its ultimate strength [11].An important physical parameter of snow is represented by the snow density.According to Shoop [12], the snow density can vary from 100 kg/m 3 (fresh snow) to 600 kg/m 3 (snowpack subjected to melt-freeze cycles), however the compacted snow that automobiles can experience is typically classified as having a density ranging from 370 to 560 kg/m 3 .In order to verify the accuracy of the FEM snow model, a load-sinkage convergence test was conducted.Figure 1 presents the mesh convergence study of the soft snow model.
In the mesh convergence study, the elements size varied from 15 mm (714 nodes) to 2 mm (38152 nodes).The load-sinkage tests are performed for a snow model that has a height of 500 mm and a diameter of 600 mm.A rigid disc with a diameter of 230 mm is deforming the snow in the vertical direction with a prescribed sinkage of 200 mm.The output of the tests is the resulting vertical force.It can be observed from the plot that the last 3 iterations offer similar results for the resulting vertical force, therefore the element size of 2 mm is selected for further study.

Bekker-Wong coefficients calculations based on the FEM pressure-sinkage results.
This section introduces a methodology for calculating the Bekker-Wong coefficients.Equation (1) provides a suitable description of the pressure sinkage relationship for a homogeneous terrain [13]: Where p is the normal pressure, b is the radius of a circular plate, z is the sinkage, and n, k c and k ϕ are the coefficients which need to be calculated for every type of soil.Some researchers consider that the coefficient k c is related to the cohesive property of the soil while k ϕ is related to the frictional properties of the soil.However, later research proves that these parameters don't necessarily have a physical meaning and their values are obtained through fitting equation (1) to the measured data.The values of k c , k ϕ , and n can be calculated from the results of a minimum of two tests with two sizes of plates having different radii (or widths of rectangular plates).The tests produce two pressure equations (p 1 and p 2 ), for each one of the plate dimensions (b 1 and b 2 ): If natural logarithm function is applied, then equation (2) becomes: For sinkage z = 1, the natural logarithm will have the value of zero, therefore the expected shape of the equation ( 3) is presented in Figure 2 [13].
To find the values for the Bekker-Wong's coefficients, two different axial compression tests have been performed on the soft snow modelled in ABAQUS.Figure 3 presents the test conditions.
In the first test, a circular rigid plate with a diameter of 230 mm was utilized.For the second test, a circular rigid plate with a diameter of 340 mm was employed.Both tests involved a quasit-static loading condition, which was applied by imposing a vertical displacement on the rigid disc.The sidewall of the model was subjected to roller boundary conditions, while the bottom face of the model was constrained in all directions.The snow mesh was created using C3D8R solid elements with reduced integration.The rigid plates were represented as analytical rigid surfaces and modeled using a planar analytical surface combined with a rigid body constraint.The density of the snow was assumed to be constant throughout the simulation, with a value of 200 kg/m³.The analysis has been carried out with Abaqus Explicit with double precision.

Estimation of Motion Resistance due to snow compaction using Bekker-Wong's coefficients
The total motion resistance encountered by a tire rolling on soft terrain comprises three components: motion resistance due to terrain compaction, motion resistance due to tire deflection, and motion resistance caused by the bulldozing effect.However, in our project, we focus solely on estimating the  motion resistance attributed to terrain compaction.This is because we assume that the snow possesses low density, and we consider the tire to be rigid.
To estimate the motion resistance due to terrain compaction, we intend to utilize the Bekker-Wong coefficients, as suggested by Wong [13].Based on these coefficients, the estimation involves assuming that the terrain reaction at all points within the contact patch is purely radial.It is further assumed to be equivalent to the normal pressure experienced beneath a horizontal plate positioned at the same depth during a pressure-sinkage test.A simplified wheel-soil interaction model is illustrated in Figure 4, which depicts this estimation approach.
The equilibrium equations for this model are: If it is assumed that the normal pressure σ acting on the wheel rim is equal to the normal pressure p beneath a plate at the same depth z, the equation to calculate R c becomes: Where b is the width of the tire, σ is the normal pressure, W is the vertical load, r is the radius of the wheel, and k c , k ϕ and n are the Bekker -Wong coefficients.

Results and discussion
This section presents the results obtained for the wheel sinkage, motion resistance and thrust capabilities of a rigid wheel rolling on soft snow.To determine the Bekker-Wong's coefficients for various snow depths, load-sinkage tests were carried out on the FEM soft snow model.Based on these coefficients, analytical equations were employed to predict the wheel sinkage and motion resistance.Subsequently, the results obtained from the analytical study were compared with those obtained from an FEM analysis of a rigid wheel rolling on the soft snow model.

Estimating Wheel Sinkage and Motion Resistance utilizing Bekker-Wong's coefficients
As presented in section 2.2, vertical load-sinkage tests were performed on the FEM snow model using two rigid discs with diameters of 230 mm and 340 mm for different snow depths.The diameters of the disks were chosen to have the same value as those used by Shoop [10] in her laboratory compression tests.Figure 5 presents the FEM compression tests results for different snow depths.
The output of the FEM analysis is the resulting vertical force for a prescribed sinkage of the circular disc.The vertical load data is converted into pressure by dividing the force value with the area of the circular plates.The curves presented in Figure 5 are resulted after post-processing the raw data from the FEM analysis using the smoothing algorithm offered in MATLAB.As presented in Figure 2, in order to find the required curves to determine Bekker-Wong's coefficients, the data needs to be further processed by applying a natural logarithm function.From Figure 5 it can be observed that the loadsinkage curves are different for each snow depth, therefore the Bekker-Wong coefficients will have different values.Table 2 presents the values for the Bekker-Wong coefficients for the different snow depth values.The coefficients found by our algorithm are different than the ones presented in the literature for snow [13].However, the values of the coefficients might not have a physical meaning and are used only to curve-fit the experimental data with the Bekker-Wong pressure-sinkage equation.The coefficients from the literature are calculated from experimental data obtained in the field where snow has different compaction levels and different densities.Our FEM snow model assumes that the snow is soft (low compaction level), homogeneous, and has a constant density.As presented in section 2.3, the Bekker's coefficients obtained from the FEM snow model will further be used to estimate the sinkage and motion resistance, considering the rigid wheel case.
To analyze the accuracy of our analytical model, we conducted a comparison between our model's output and the empirical data collected in the field by Shoop [10].In order to do that, we employed identical wheel dimensions, D = 29.1"(0.74 m) and b = 10.7"(0.27 m).The vertical load is considered to be W = 1402 lbs (6236 N). Figure 6 presents a comparison between the theoretical and the empirical data for wheel sinkage and motion resistance coefficient.The motion resistance coefficient was calculated as a ratio between the longitudinal motion resistance force and the vertical load.
As it can be observed from the empirical data [10], at a snow depth of 0.4 m, the measured sinkage was 0.3 m (11.8") while the motion resistance coefficient was 0.2.At the same snow depth, our theoretical model output is given a sinkage of 0.29 m and a motion resistance coefficient of 0.47.Therefore, we can conclude that the accuracy of the theoretical model regarding the sinkage results is fairly good with respect to the measured values; however it can be observed that the motion resistance coefficients obtained using the theoretical model are significantly higher than the ones measured in the field.

Estimating Wheel Sinkage and Motion Resistance utilizing FEM snow model
Using the snow model described in section 2.1 (ALE method with C3D8R elements), we conducted an ABAQUS Explicit Analysis to examine the behavior of a rigid wheel rolling on soft snow.The purpose of this analysis was to investigate the changes in sinkage and motion resistance at different snow depths.
For this FEM analysis, the same wheel dimensions and properties have been used as for theoretical analysis in section 3.1.The ABAQUS Explicit software does not offer an input for the desired slip ratio, therefore the input longitudinal velocity and angular speed must be carefully chosen in order to simulate a desired slip ratio.For the analysis of wheel sinkage and motion resistance, a longitudinal speed of 2 m/s and an angular speed of 5.41 rad/s was applied to the rigid tire, therefore the simulated slip ratio was 0%.The boundary conditions for the snow model was selected as follow: the bottom face of the model was constrained in all directions and for the sidewalls roller boundary conditions were used to allow displacement only in the vertical direction.The simulated snow depth values were 0.2, 0.3, 0.4, and 0.5 m. Figure 7 presents the raw FEM output data for the longitudinal resistance force.

Figure 7. Effect of Snow depth on Longitudinal resistance force
Due to computational resources, for this type of explicit analysis, the mesh is coarse, therefore the output reaction force in the longitudinal direction varies significantly with respect to the simulation time.The resistance force data is post-processed by averaging the raw output data from the FEM model.After that, the motion resistance is calculated for every snow depth by dividing the resistance force to the vertical load.The raw output data regarding the wheel sinkage is relatively constant and it is processed in a similar manner.Figure 8 presents a comparison between the output of our own FEM model, the FEM model and the experimental results presented by Shoop [10], and the results from the analytical model of sinkage and motion resistance based on Bekker-Wong's coefficients.
It can be observed from Figure 8 that the sinkage output of our FEM model and the analytical model is fairly similar to the experimental data recorded in the field.This indicates that the rigid wheel approximation does not impact in a significant manner the accuracy of the predicated sinkage.The motion resistance coefficient of our FEM model is different from the FEM model and the experimental data presented by Shoop; however, it is similar to the analytical model.As presented before, the theoretical data is obtained by using the Bekker-Wong's coefficients, which were calculated based on the FEM snow model.Therefore, it is expected that the output of our FEM model and the output of the analytical model to be similar.In order to improve the output of our FEM model, further refinement of the mesh should be considered for the 3D explicit analysis.

Estimating Thrust capabilities utilizing FEM snow model
According to Wong [13], by integrating the horizontal component of tangential stress over the entire contact area, the total tractive effort F can be determined using equation (10).Also, the total motion resistance can be determined by integrating the pressure values over the entire contact area.If the total motion resistance is subtracted from the total thrust, the total drawbar pull force can be estimated using equation (11).The data for the pressure and shear stress distribution over the contact area is subtracted from the FEM Explicit analysis and then it is further post-processed using MATLAB.For this analysis the snow depth was 0.5 m. Figure 9 presents the evolution of pressure and shear stress distribution over the contact length of the rigid wheel for different slip ratios.It can be observed from Figure 9, that the variation of pressure and shear stress is similar to the one presented in the literature [13].By applying equations (10) and (11) to the data presented in Figure 9, the total thrust capabilities and the total resistive force can be predicted.Figure 10 presents the evolution of total thrust and resistive force with respect to slip ratio.
From Figure 10 it can be observed that the maximum thrust has a value in the close range of 4 kN for a slip ratio in the range of 10-20%, while the resistive force is estimated to have the same value for the same slip ratio range.This means that the drawbar pull has a small value, close to the limit where the wheel cannot develop enough thrust to overcome the motion resistance force.As presented in section 3.2, our FEM output produces considerable higher values for the motion resistance coefficient than those presented in literature, therefore the drawbar pull might have a higher value than the one presented in this section.Also, for this analysis, the wheel was estimated to be rigid.It is known that for snow rolling conditions, the pneumatic tire generates most of the thrust from the hysteresis component of the total longitudinal force, not from the adhesion component.Also, the rigid wheel has no tread attached.Therefore, the true values of the total thrust and the drawbar pull will be higher for a real-world scenario.

Conclusions
The implementation of the snow model using the explicit version of the commercial software ABAQUS followed a methodology similar to that of Shoop [10].To ensure accurate results, a mesh convergence study was conducted by varying the element sizes from 15 mm (714 nodes) to 2 mm (38,152 nodes).
The mesh convergence study ensures that subsequent analyses are performed with an appropriate element size, leading to reliable and precise simulations of load-sinkage behaviour in soft snow conditions.A comprehensive approach was employed to determine the values for Bekker-Wong's coefficients using the load-sinkage results from a soft snow FEM model.Two axial penetration tests were conducted, each utilizing a circular rigid plate with different diameters: 230 mm for the first test and 340 mm for the second test.The purpose of these tests was to evaluate the behaviour of the snow under quasistatic loading conditions.These coefficients play a crucial role in modelling the mechanical behaviour of the soft snow and are essential for accurate simulations and predictions.
The analysis aimed to investigate the behaviour of a rigid wheel rolling on soft snow with different depths, using both an FEM model and an analytical model.The primary focus was on examining the variations in sinkage and motion resistance at different snow depths.The results obtained serve as a basis for further refinement and validation of the analytical model, enhancing our understanding of wheel-snow interactions in practical scenarios.By comparing the empirical data [10] with theoretical model's result, it is observed that at a snow depth of 0.4 m, the measured sinkage was 0.3 m (11.8"), while the motion resistance coefficient was 0.2.In contrast, our theoretical model predicted a sinkage of 0.29 m (11.4") and a motion resistance coefficient of 0.47 at the same snow depth.Consequently, based on the findings of the study, it can be concluded that the theoretical model demonstrates a good level of accuracy in predicting the sinkage results when compared to the corresponding measured values.
However, it is noteworthy that the motion resistance coefficients obtained using the theoretical model are slightly higher than those measured in the field.On the one hand, the analysis of the sinkage output revealed that all models, including both FEM models and the analytical model, produced fairly similar results.This suggests that the approximation of a rigid wheel does not significantly impact the accuracy of the predicted sinkage.On the other hand, the output for the motion resistance coefficient exhibited variations between the two FEM models, although it was similar to the analytical model.As previously mentioned, the analytical model's data was obtained using Bekker-Wong's coefficients, which were calculated based on the FEM snow model.Consequently, it is expected that the output of our FEM model and the analytical model would be relatively similar.
To improve the performance of our finite element method (FEM) model, we plan to enhance the mesh refinement for the 3D explicit analysis.By incorporating this refinement, a higher level of detail regarding the behavior of the snow can be captured, resulting in enhanced accuracy in predicting the motion resistance coefficient.The optimization of the mesh size allows for a more accurate representation of the snow behavior, thus increasing the fidelity of the results obtained from the FEM model.

Figure 1 .
Figure 1.Mesh convergence study of the soft snow model.

Figure 3 .
Figure 3. Axial compression tests conditions performed on the soft snow model.

Figure 5 .
Figure 5. FEM compression tests results for different snow depths

Figure 6 .
Figure 6.Wheel sinkage and motion resistance coefficient comparison

Figure 9 .
Figure 9. Pressure and Shear Stress distribution over the contact length of the rigid wheel

Figure 10 .
Figure 10.Variation of Total Thrust and Resistive Force with Slip Ratio

Table 1 .
Coefficients used to characterize the soft snow material properties.

Table 2 .
Bekker-Wong coefficients for the different snow depth values.