Multi-stage terramechanics simulation: Seamless analyses between formation of wind ripple pattern and wheel locomotion

Abstract In futuristic vehicle traveling simulation on a rough terrain, it is desirable to perform terramechanics analysis that reflected appropriate road surface information. This paper presents a multi-stage analysis method that seamlessly performs the terrain formation process and trafficability evaluation. Simulation of wind-blown ripple formation based on cellular automaton was adopted as an example, and fields were created for terramechanics analysis. Then, trafficability characterization was performed with a single rigid wheel on virtually created terrain fields based on resistive force theory. The analysis results showed that the traveling performance of a single wheel could be evaluated for various terrain fields and under various traveling conditions. The multi-stage analysis method presented in this paper can be employed for simulations in extreme environments such as planetary surfaces, deserts, and disaster sites, where sensing is difficult.


PUBLIC INTEREST STATEMENT
Terramechanics is an interdisciplinary field that deals with the interaction between the ground and vehicles, and is effective for trafficability characterization. In future more advanced vehicle traveling simulation on a rough terrain, it is desirable to perform analysis reflected appropriate ground surface information. However, the terrain fields in most conventional analysis assumed homogeneous and flat surface conditions. Off-roads have different surface geometrical conditions depending on the formation process, which should be properly reflected in simulations of vehicles. This paper proposes a multi-stage analysis method that seamlessly performs the terrain formation process and trafficability evaluation. Simulation of wind-blown ripple formation was adopted as an example, and fields were created for terramechanics analysis. Then, trafficability characterization was performed with a single rigid wheel on virtually created terrain fields. The proposed method can be employed for simulations in extreme environments such as planetary surfaces, deserts, and disaster sites.

Introduction
Terramechanics is an interdisciplinary field that deals with the interaction between the ground and vehicles (Bekker, 1956;Vantsevich, 2015;Wong, 2008;WONG, 1997). The target of terramechanics includes all off-road vehicles such as automobiles, construction machinery, agricultural machinery, and military vehicles. Recently, terramechanics has attracted renewed attention and is being applied to problems in extreme environments such as disaster response robots and lunar/planetary exploration rovers (Schäfer et al., 2010;Sharma et al., 2018;Sutoh et al., 2013;Tadokoro, 2019;Trease et al., 2011;Zhou et al., 2014). Terramechanics analysis for a wide variety of rough terrains is classified into multibody dynamics analysis using semi-empirical theories (Arvidson et al., 2017;Raymond & Jayakumar, 2015;Zhou et al., 2014) and numerical analysis using finite element method (Chiroux et al., 2005;Fervers, 2004;González Cueto et al., 2016;Nishiyama et al., 2016;Ozaki et al., 2015;Ozaki & Kondo, 2016) or discrete element method (Johnson et al., 2015;Knuth et al., 2012;Li et al., 2010;Nakashima et al., 2010Nakashima et al., , 2007Nishiyama et al., 2016;Smith & Peng, 2013). In particular, multibody dynamics analysis can be a powerful tool for a wide variety of rough terrains, not only to improve the performance of off-road vehicles and robots but also to improve the efficiency of various tasks and create safe working plans.
To carry out the multibody dynamics analysis for practical applications, it is essential to establish a field modeling method for rough terrain, in addition to improving the prediction accuracy of the state-of-the-art terramechanics theories. In the past, research efforts have focused on the advancement of semi-empirical terramechanics theories for the interaction between the ground and vehicles (Bekker, 1956;Du et al., 2018;Ishigami et al., 2007;C Li et al., 2013;Wong, 2008). On the other hand, it is necessary to prepare the target terrain field in advance before carrying out terramechanics analysis; however, few studies have been conducted on how to create such terrain fields. Although the terrain fields in conventional terramechanics analysis are prepared based on various assumptions and experiences, most of them have homogeneous and flat surface conditions. Furthermore, off-road conditions may have different geometries and soil properties depending on the formation process, which should be properly reflected in simulations of the terrain field.
Sensing data can be used to create realistic terrain fields. For example, high-precision topographic data can be obtained by aerial laser scanning such as LIDAR, and a point cloud of the topographic surface layer can be obtained by ortho mapping using satellite images or images taken by drones. A terrain field for terramechanics analysis can be created using point cloud. In addition, the mechanical properties of soil in the target field can also be determined by performing the Bevameter test or various penetration tests.
Meanwhile, numerous studies have been conducted on numerical simulations related to the formation process of the terrain such as the formation of wind ripples in a desert (Guignier et al., 2013;Katsuki et al., 2011;Nishimori & Ouchi, 1993), washboarding of dirt road due to repeated traffic (Bitbol et al., 2009;Taberlet et al., 2007), and soil compaction by rollers in construction work (Kenneally et al., 2015;Pietzsch & Poppy, 1992). Furthermore, in terms of disaster, numerical simulations for slope failure, landslide, and avalanche have also been reported (Bhandari et al., 2016;Di Gregorio et al., 1999;Miles & Ho, 1999;Norem et al., 1989;Sailer et al., 2008;Sampl & Zwinger, 2004;Wang & Sassa, 2010). Therefore, if terrain fields can be created for terramechanics analysis using the results of these terrain formation simulations, they can be used for trafficability evaluation of the terrain and path planning in various situations.
In this study, using the simulation of wind-blown ripple formation as an example, a modeling method of a terrain field, which uses the soil surface geometry, is proposed for terramechanics analysis. In addition, in order to propose a method for multi-stage analysis for achieving seamlessness with terramechanics analysis, we demonstrated the systematic analysis of single wheel locomotion using numerically created sand ripple fields. For single wheel locomotion analysis, we adopted the resistive force theory (RFT), which has attracted attention in recent years (Agarwal et al., 2019;Askari & Kamrin, 2016;C Li et al., 2013;Slonaker et al., 2017;Suzuki et al., 2019). RFT can be employed to evaluate the resistive stress generated in an arbitrarily shaped object moving in granular media. It has been demonstrated that RFT is suitable for the traveling analysis of a legged mobile robot (C Li et al., 2013), and this method has been employed to evaluate the force generated in wheels (Agarwal et al., 2019;Slonaker et al., 2017;Suzuki et al., 2019). We adopted the programming language MATLAB and its numerical computing environment to perform the simulations of wind ripple formation and wheel locomotion.
The multi-stage analysis method proposed in this study can be used for simulations in extreme environments such as planetary surfaces, deserts, and disaster sites, where sensing is difficult.

Simulation of ripple formation by wind-blown sand
This section presents the simulation of the formation of ripple patterns by wind-blown sand based on the cellular automaton model proposed by Nishimori and Ouchi (Nishimori & Ouchi, 1993) as the first stage of seamless analysis. The model is a discrete model in space and time with a continuous field variable representing the averaged surface height at each cell. The simulation corresponds to a small-scale model, and a ripple pattern is spontaneously formed by the jumping process of sand grains, known as saltation, when the wind force exceeds a critical value.

Model
It is assumed that the wind ripple pattern is created by two processes, namely, saltation and creep of sand grains (Nishimori & Ouchi, 1993). When sand grains on the surface are released into the air under strong wind, the grains are accelerated by the wind. Further, the grains impact on the surrounding grains, releasing them into the air. As a result, sand in a certain area jumps and moves to another area. The dynamics of the saltation process is expressed as follows: where h is the height of the sand surface at each area (cell of the ith row and the jth column), and q is the transferred height. Eq. (1) indicates that when the wind blows in the ith direction, sand with a height of q in cell position (x i , j , y i , j ) is released, and is added to another cell position (x i + Δi , j , y i + Δi, j ). The subscript Δi is the flight length at one saltation, and is expressed using parameters l 0 and b. Using the floor function, Δi is given as follows: where l 0 is a control parameter proportional to the wind force, b is a constant related to the average wind velocity a grain experiences in flight, and dx is the lattice length in the ith direction (x direction) of a single cell.
Meanwhile, moving (rolling) to a nearby region when grains cannot maintain their position due to a steep surface gradient is called creep. The dynamics of creep is described by the following equation.
Eq. (3) has the following meanings: The sandhill in a certain cell is relaxed by gravity at a speed proportional to the convexity of the sand surface and receives inflow from adjacent cells at the same time. Here, r is the rate of relaxation.
The distribution of inflow is rh/6 from perpendicular crossing cells and rh/12 from diagonally adjacent cells.
The inflow term from the perpendicular directions, H cross , and the inflow term from the diagonal directions, H diag , are given as follows: It was assumed that creeping action occurred when the height difference between adjacent cells exceeds the angle of repose.

Analysis results
In the numerical simulation, the dynamics of saltation and creep are evaluated at each time step, and a ripple pattern is formed by repeated calculation. In this study, analysis was performed under five conditions labeled A-E, as listed in Table 1. Here, dy is the lattice length of a single cell in the jth direction and was set to be the same as dx. The number of repeated calculation steps is 30000 for each condition. The length of the analyzed soil surface is 1500 mm (row i direction), and the width is 400 mm (column j direction). Figure 1 shows the variation of the soil surface geometry with the number of repeated steps k, where condition C in Table 1 was used. The soil surface shape is irregular at k = 1 as the initial configuration of the height of the cells was set using random numbers. At k = 1000, the soil surface is close to being flat due to creeping action, and the formation of wind ripples is initiated due to saltation. Then, the wind ripple pattern starts becoming clear (k = 5000), and eventually steady state is reached (k = 30000). In this study, the wind ripple pattern at k = 30000 was used as the steady state and was reflected in the terrain field modeling for terramechanics analysis, which will be presented later.
Next, the effect of each parameter on wind ripple formation is verified. As presented in Table 1, we examined the effects of l 0 and q, which affect the dynamics of saltation. Figure 2 shows the wind ripple patterns under conditions A-E in steady state. As shown in conditions A, B, and C of Figure 2, l 0 affects the wind ripple interval because of the distance of saltation changes. Meanwhile, it was verified that variations in q affect the height of the ripples. In the case of condition D with a small q, the ripple pattern is unclear, whereas it becomes quite clear in condition E, which has a large q. It should be noted that the distributions of the initial height using random numbers were the same for conditions A-E.

Resistive force theory (RFT)
In this study, RFT (C Li et al., 2013) was adopted to analyze rigid wheel locomotion on a rough terrain. Although it is difficult to consider the movement and shearing of soil (Suzuki et al., 2019), RFT can be employed to evaluate the resistive force generated in an object with an arbitrary shape and motion at a low calculation cost.
According to RFT, a horizontal stress, σ x and a vertical stress, σ z on a small plate inside granular media are proportional to the sinkage z, and can be defined as follows: where subscripts x and z represent the coordinate axes. The stiffnesses α x and α z of the granular media depend on the orientation angle β of the plate and the velocity vector angle (intrusion/ Figure 1. Variation of soil surface geometry with a number of repetition steps k, where the condition C in Table 1 was adopted. The directions of row i and column j correspond to the x and y coordinates, respectively. The unit of each axis is in m.
extrusion angle) γ. The relationship between β, γ, and α can be obtained from the plate intrusion/ extrusion test of the target granular media (C .
To express the data of α x and α z for discrete β and γ as a continuous distribution, an approximation was adopted based on the following Fourier transformation.
where the coefficients A, B, C, and D are obtained by Fourier transformation. The values of the Fourier coefficients used in this study were obtained from Li et al. (C Li et al., 2013), and are listed  Figure 3, which corresponds to the stiffness distribution, can be adjusted based on the scale factor.
When applying RFT to a rigid wheel, the wheel surface is discretized into small sections (plates), and β and γ of each section are evaluated during locomotion. As shown in Figure 4, each plate has a rotational angular velocity ω and a translational velocity v. γ can be determined from these combined velocity vectors, whereas β can be determined from the posture of each plate.
The stress distribution of the contact surface can be obtained by substituting β and γ of each plate into Eq. (6). Further, the horizontal force (drawbar-pull), F x and the vertical force, F z are calculated by performing area integration of the stresses of the piecewise plates. Meanwhile, the sinkage z can be evaluated from the balance between the prescribed wheel load, W and the vertical force, F z calculated from RFT.
In the wheel locomotion analysis, the wheel radius, width, and load were set to R = 100 mm, B = 100 mm, and W = 120 N, respectively, unless stated otherwise. Further, the rotational angular velocity, ω, which is a source of driving for the wheel, was set to 0.5 rad/s. The translational velocity, v is calculated by solving the following equation of motion.
where m and a are the mass and acceleration of the rigid wheel, respectively. P is the towing load in the wheel backward direction and was fixed at 15 N. In addition, an important indicator of offroad trafficability is the slip ratio, s, which is defined as follows: Figure 4. Application of RFT to a rigid wheel, where da is the area of a small plate. Quadrangular prisms in the background of the (x, z) plane are the distribution of the cell height, which constitute the ripple pattern.
Here, s = 1.0 corresponds to the wheel being stuck.

Terrain field modeling
In this study, terrain field models were created using the analysis results of wind ripples under A-E conditions shown in Figure 2. The terrain fields for terramechanics analysis were discretely expressed since the simulation of wind ripples formation uses a discrete model in the spatial domain. In other words, the field was discretely modeled using a quadrangular prism corresponding to the cell size (dx, dy) of the ripple pattern analysis, as shown in Figure  4. This makes it possible to model a surface topography using the simulation results of the wind ripples.
It should be noted that to evaluate the sinkage in RFT, it is necessary to calculate the height difference between the ground surface and each plate, which constitute the wheel surface. Therefore, in this study, the quadrangular column where the center coordinate (x, y) of the small plate belongs were determined; then, the sinkage was evaluated using the height of the target quadrangular column. Here, the rigid wheel was discretized with small plates using intervals of 0.01 rad.
On the other hand, in terms of the terrain field modeling, it is necessary to set the soil parameters in addition to the surface topography. In RFT, the scale factor is the only soil parameter, and thus we set ζ = 1.0, assuming soft sandy soil (Suzuki et al., 2019). Figure 5 shows the relationship between the drawbar-pull and the slip ratio when the wheel specifications already described and ζ = 1.0 were used. Here, to determine the characteristics of the analysis conditions, we performed the wheel traveling analysis on a flat road surface in advance. It can be confirmed that the drawbar-pull of the wheel increases as the slip ratio increases, which is consistent with results in existing reports (Ding et al., 2011;Ishigami et al., 2007).

Trafficability characteristics of flat road surface
First, the analysis results of wheel locomotion on a flat road surface based on the equation of motion given in Eq. (7) are presented. Here, the wheel travels straight in the x-direction, which is achieved by using a constant rotational angular velocity, ω. Figure 6(a) shows the trajectory of the wheel center, whereas Figure 6(b) shows the variation of the slip ratio and the travel distance with the elapsed time. It can be observed from the relationship shown in Figure 5 that the slip ratio corresponding to a towing load of P = 15 N is approximately 0.38, indicating that the wheel shows steady traveling while balancing. In addition, after the drawbar-pull, F x exerted by the wheel balances with the towing load, P, a constant velocity linear motion is achieved.

Trafficability characteristics of wind ripples
Next, the wheel behavior when traveling in the x-direction on terrain fields obtained by wind ripple simulation is presented. Here, the traveling position of the wheel in the y-direction is the center of the field.
Figure 7(a) shows the trajectory of the wheel center obtained by terramechanics analysis using field C (Figure 2(c)). Figure 7(b) shows the variation of the slip ratio and the traveling distance with the elapsed time. It can be observed from the figures that the wheel repeatedly moves up and down according to the undulation of the terrain field. In addition, the slip ratio and the traveling speed of the wheels also vary according to the surface undulations. Figure 8 shows the stress distributions at t 1 and t 2 depicted in Figure 7. The lines represent the stress vector on the wheel surface. The vertical blue lines represent σ z ; the red and orange lines represent σ x acting in the traveling direction and its opposite direction, respectively. At time t 1 when the slip ratio reaches its maximum value, the traveling resistance due to the front slope is large, and the slip ratio increases to produce drawbar-pull to overcome this (Figure 8(a)). On the other hand, time t 2 is when the slip ratio reaches its minimum value. At this time, there is minimal traveling resistance, as the descent along the slope of the ripple is initiated, and the horizontal stress required for driving is small (Figure 8(b)). Therefore, it is possible to move forward even at a low slip ratio.   In addition, we will examine the trafficability characteristics of terrain fields A, B, D, and E using different parameters in the wind ripple simulation. Figure 9(a,b) shows variations of the slip ratio and the traveling distance with the elapsed time for each terrain field condition. It can be observed from the figure that the characteristics of the slip ratio and the traveling distance vary according to the terrain fields. In particular, as shown by the blue and black lines in Figure 9(a,b), the traveling characteristics (x-displacement and slip ratio) are not periodic on the irregular fields with undeveloped wind ripples such as the conditions A and D. On the other hand, the wheel becomes stuck in terrain field B (Figure 9(c)). The reason is that the traveling resistance due to ripples in front of the wheels is quite large, and it was difficult to produce drawbar-pull to drive forward under the wheel specification used in the analysis.

Effects of wheel specification and traveling direction
As described in Section 4.2, trafficability characterization of the wheel can be performed by reflecting the terrain geometry through terramechanics analysis of wind ripples formed under various conditions. The representative applications of terramechanics analysis are in wheel design and path planning. Therefore, this section presents a case study of trafficability characterization when the wheel specifications and motion conditions are changed for wind ripples formed under Figure 10. Trajectory of the wheel center for field C while varying the wheel radius.
the specified conditions (condition C). Here, the wheel load, W and the towing load, P were fixed at 120 N and 15 N, respectively.
First, the effect of the wheel radius, R was examined. In the analysis, the wheel radius was set to R = 50 mm, 100 mm, and 200 mm, and the trafficability data were evaluated. Here, ω = 1.0 rad/s, 0.5 rad/s. and 0.25 rad/s were set for each radius to obtain the same theoretical vehicle speed (corresponding to s = 0.0). Figure 11. Terramechanics analysis results for field C while varying the wheel radius: (a) variation of the slip ratio with the elapsed time; (b) variation of the x-coordinate of the wheel center with the elapsed time. Figure 10 shows the trajectory of the wheel center. Figure 11(a,b) shows variations of the slip ratio and the traveling distance with the elapsed time. It can be observed from the figures that the larger the wheel radius, the higher the traveling performance. In general, the larger the wheel radius, the smaller the contact pressure and the lower the traveling resistance. Moreover, the relative wheel radius with respect to the ripple height increases, producing low traveling resistance during climbing. As a result, the slip ratio required for traveling decreased and the wheel speed increased. On the other hand, when the wheel radius is small, the slip is required to exert the traction necessary for traveling, and the wheel speed is reduced. It should be noted that increasing the wheel radius increases the weight of the vehicle; thus, there is a trade-off relationship between energy consumption and payload during transportation.
Next, the effect of the wheel load, W is examined. In the analysis, the wheel load was set to W = 60 N, 120 N, and 240 N, and the trafficability data were evaluated. Figure 12(a,b) show variations in the slip ratio and the traveling distance with the elapsed time. It can be observed from the figures that the lower the wheel load, W, the lower the traveling performance; at W = 60 N, the wheel becomes stuck. This is because the maximum drawbar-pull that can be exerted at low wheel loads is reduced. Thus, the drawbar-pull required to overcome the slope cannot be exerted. On the other hand, a comparison of results obtained at W = 120 N and those of 240 N shows that there is minimal difference in the traveling speed, although there is a difference in the slip ratio under field condition C.
Finally, analysis results when the direction of travel is changed are presented. In the above analysis, although the wheel traveled perpendicular to the wind ripples (x-direction), it traveled straight at an angle, θ tr around the z axis. Figure 13 shows the trajectory of the wheel center. Figure 14(a,b) shows variations of the slip ratio and the traveling distance with the elapsed time, Figure 13. Trajectory of the wheel centre for field C while varying the transverse angle.
respectively. Here, the transverse angle, θ tr was set to 0 deg, 5 deg, and 10 deg. As shown in the figures, the periods of fluctuation of the slip ratio are different because the time to reach the slope and the contact angle are different. In addition, the effect of the transverse angle gradually increases with the traveling distance. Further, the traveling speed is also affected due to variation in the slip ratio. From the above, it can be concluded that terramechanics analysis using the wind ripple formation process under various conditions is effective for the robust design of vehicles (Sharma et al., 2018;Sutoh et al., 2013) and mobile robots (C Li et al., 2013) and optimal path planning (Arvidson et al., 2017;Schäfer et al., 2010;Trease et al., 2011;Zhou et al., 2014). The multi-stage analysis method, from terrain formation to vehicle locomotion, presented in this paper can be employed for simulations in extreme environments such as planetary surfaces, deserts, and disaster sites, where sensing is difficult.

Conclusion
The conclusions of this study are summarized below.
• In this study, we proposed a multi-stage analysis method that seamlessly performs the process of ground surface formation and trafficability evaluation. In particular, we used wind ripple simulation for rough terrain modeling as an example and performed the simulation of single wheel locomotion on virtually created ground surfaces using RFT. Then, we performed a terramechanics analysis of wheels on terrain fields created under various conditions. The proposed analysis method can be employed for an advanced off-road vehicle traveling analysis in the future because various terrain fields can be considered.
• The proposed method of multi-stage terramechanics simulation can be used with other numerical analyses and measurement technologies. For example, information on ground surface deformation can be analyzed in detail for the finite element method or smoothed particle hydrodynamics, and the strain (density) information obtained from the analysis can be used for terrain field modeling and subsequent terramechanics analysis.
• It should be noted that the extension of terramechanics theory for rough terrain needs to be continued, because it is difficult to consider the movement and shearing of soil in conventional theories.