Online estimation of terrain parameters and resistance force based on equivalent sinkage for planetary rovers in longitudinal skid

Mechanical Systems and Signal Processing 119 (2019) 39–54

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Online estimation of terrain parameters and resistance force based on equivalent sinkage for planetary rovers in longitudinal skid Zhen Liu a, Junlong Guo a,b,⇑, Liang Ding a, Haibo Gao a, Tianyou Guo c, Zongquan Deng a a

State Key Laboratory of Robotics and System, Harbin Institute of Technology, Harbin 150001, China Department of Mechanical Engineering, School of Naval Architecture and Ocean Engineering, Harbin Institute of Technology (Weihai), Weihai 264209, China c Mechanical Engineering, University of Michigan, Ann Arbor, MI 48109, USA b

a r t i c l e

i n f o

Article history: Received 9 June 2018 Received in revised form 10 August 2018 Accepted 12 September 2018

Keywords: Wheeled mobile robots Equivalent sinkage Analytical wheel-soil interaction modeling Online estimation of terrain parameters and resistance force

a b s t r a c t Wheel-soil interaction mechanics plays a crucial role for wheeled mobile robots (WMR) on rough and deformable terrains such as Martian and Lunar surfaces. Skid terramechanics is an essential component for WMRs and generates resistance force when a WMR brakes or on downhill slopes. The basis of classical terramechanics theories for WMRs – Bekker’s normal stress and Janosi’s shear stress equations – are so complex that the wheel-soil interaction force/torque equations are not amenable to closed form solutions, which seriously limits the application of terramechanics theories to WMRs. To establish analytical wheel-soil interaction expressions, the normal and shear stresses that can be characterized linearly by the proposed terrain stiffness and shear strength, respectively, are presented in this paper. Terrain stiffness and shear strength can be used to characterize terrain mechanical properties. Compared with the experimental data, the maximum relative error of the resistance forces estimated using these expressions at steady state is less than 7%. These validated expressions can be applied to estimate terrain parameters and resistance force online with high accuracy. Terrain’s stiffness and shear strength increase first, and then reach a constant. Before wheels entering steady state, the online estimated resistance force’s relative error is much higher, which can be explained using wheel’s vertical velocity. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction The analysis of wheel-terrain interaction mechanics has implications for the system’s design [1,2], sensing subsystem [3,4], and estimation and control algorithms [5–7]. In most cases, this interaction is assumed to follow a simple Coulomb friction law [8–11], and the effects of such phenomena as wheel skid and vertical sinkage are ignored [12]. Although such an approach may be sufficient for some applications, operation near a system’s performance limits - that is on challenging terrain - often requires more sophisticated analyses of wheel-soil interaction [13]. The images that were sent back by the YUTU Lunar Rover and NASA’s Mars Rovers show that the Lunar and Martian surfaces are covered with fine and soft soil which are difficult to traverse, and their access presents an ongoing challenge for WMRs [14]. In the coming 2020, CNSA (China National Space Administration) [15], NASA [16], and the European Space ⇑ Corresponding author at: Department of Mechanical Engineering, School of Naval Architecture and Ocean Engineering, Harbin Institute of Technology (Weihai), Weihai 264209, China. E-mail address: [email protected] (J. Guo). https://doi.org/10.1016/j.ymssp.2018.09.017 0888-3270/Ó 2018 Elsevier Ltd. All rights reserved.

40

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

Nomenclature A C FD FP FR FRr FRs FV FVr FVs K K1, K2 Kr Ks TB W b c1, c2 c d1, d2 f1(h) f2(h) fo(t) hL j k1, k2, k3 k11, k12 k21, k22 kc ku N n0, n1, n2 nL p r sd t v z z1 z2 zF zM zV zr

a x xO

h h0 h1 h2 hm kPRC

s

u /

r r1 r2

amplitude of oscillation intercept when oscillation frequency is plotted versus skid ratio Damp force caused by wheel vertical velocity (N) pushing force from the vehicle (N) resistance force (N) resistance force component integrated from the normal stress (N) resistance force component integrated from the shear stress (N) vertical force (N) vertical force component integrated from the normal stress (N) vertical force component integrated from the shear stress (N) shearing deformation modulus of soil (mm) coefficients used to calculate resistance force terrain stiffness (Pa/m) shear strength (Pa/m) braking torque (Nm) wheel vertical load (N) wheel width (mm) coefficients used to compute maximum normal stress angle soil cohesion (Pa) coefficients used to determine angular position of the transition point of shear stress function used to calculate equivalent wheel sinkage function used to calculate normal stress function used to describe oscillation lug height (mm) soil shearing deformation (mm) coefficients used during the analytical modeling process coefficients used to calculate equivalent wheel sinakge coefficients used to calculate normal stress cohesive modulus of soil (kPa/mn
1) frictional modulus of soil (kPa/mn) sinkage exponent of soil coefficients used to compute sinkage exponent lug number pressure between plate and soil (Pa) wheel radius (mm) skid ratio time (s) wheel forward velocity (mm/s) plate sinkage (mm) wheel vertical sinkage (mm) soil rebounding height (mm) wheel vertical sinkage caused by the vertical force (mm) measured wheel vertical sinkage (mm) wheel vertical sinkage caused by its vertical velocity (mm) equivalent wheel sinkage (mm) slope angle (rad) wheel angular velocity (rad/s) oscillation frequency (rad/s) wheel-soil contact angle (rad) angular position of the transition point of shear stress (rad) wheel entrance angle (rad) wheel leaving angle (rad) angular position of the maximum normal stress (rad) resistance force coefficient tangential stress (Pa) soil internal friction angle (°) phase shift (rad) normal stress (Pa) normal stress in the front region (Pa) normal stress in the rear region (Pa)

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

41

Agency (ESA) in cooperation with Roscosmos [17], are going to explore the Red Planet by landing WMRs, as illustrated in Fig. 1. The use of rovers in missions significantly increases the exploration area, thereby improves the scientific return of the mission. To prevent the robots from getting stuck on the deformable terrain like the Spirit Mars Rover, it is necessary to conduct thorough investigations of the contact mechanics between the wheel and the soil, and eventually, have a better understanding of the motion behavior of the wheel on the loose soil. Slip between wheel and soil generates tractive force, while skid produces resistance force. As the slip terramechanics directly affects the traveling performance, and obtained extensive attention in the past several decades, mainly including the wheel performance estimation [18,19], dynamic oscillations prediction [20], slip-sinkage evaluation and its effect on the tractive performance[21], arising drawbar pull using active lugs [14], evaluating the influence of configuration parameters on rover’s traveling performance [2]. Compared with slip terramechanics, skid terramechanics which can determine the ability of downhill a slope by computing the resistance force that balances the gravity component received little attention. Traditional wheel-terrain interaction models rely on terrain-specific parameters, which must be correctly estimated to accurately predict rover mobility. In classical terramechanics approaches, it is typically necessary to conduct a series of direct shear tests and a series of pressure-sinkage tests to obtain these parameters [22,23], which are not available for Mars/Lunar regolith. To estimate the Martian parameters, the Viking landers employed manipulator arms to perform soil trenching experiments [24], and the Sojourner rover used the rover wheel as a trenching device to identify soil cohesion and internal friction angle [25]. Both missions used visual cues and earth-based analysis techniques, which leads to lengthy communication time delays that reduce rover efficiency and limit its autonomy. To estimate the cohesion and internal friction angle online for WMRs, a simplified form of classical terramechanics equations was proposed [6]. However, even the classical terrain parameters (e.g., sinkage exponent, cohesion, and internal friction angle) can be obtained using specially designed apparatus, to reflect the slip-sinkage [26] and skid-sinkage [27] phenomena, the sinkage exponent cannot be utilized directly, but was assumed to be a function of the skid ratio. In this paper, the basis of wheel-soil interaction modeling - the normal and share stress - is characterized linearly with the proposed equivalent wheel sinkage using the proposed terrain mechanical parameters, which can be computed in real time by running a wheel on the terrain. Then, the online estimated terrain parameters can be used to compute the resistance force in real time. The remainder of this paper is organized as follows. Part 2 concentrates on the introduction of the limitations of traditional terramechanics models. The process of establishing the closed-form analytical equations that rely on the linear normal and share stress is presented in Part 3. The experimental validation results are discussed in Part 4. The online applications are performed in Part 5. Finally, the conclusions and discussions are presented in Part 6.

(a) Mars exploration rover designed by CNSA

(c) Mars Rover of the ExoMars program carried out by the ESA and Roscosmos (b) NASA’s Mars 2020 Rover Fig. 1. Future Mars exploration rovers.

42

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

2. Limitations of traditional theories of skid terramechanics In classical terramechanics, skid ratio is a key state variable to describe the wheel running condition, and its definition is shown in Eq. (1), where r denotes wheel radius, v and x denote wheel forward velocity and angular velocity, respectively.

sd ¼ 1


rx

v

0 < rx < v

ð1Þ

The force diagram and the normal and shear stress distribution beneath a driving wheel in skid is depicted in Fig. 2, where W is the wheel vertical load, FV is the vertical force, FP is the pushing force from the vehicle, FR is the resistance force from the soil, TB is the braking torque, z1 is the wheel vertical sinkage, z2 is the rebounding height of soil, h1 is the entrance angle that can be computed using the wheel vertical sinkage as shown in Eq. (2), and h2 is the leaving angle.


z1  h1 ¼ arccos 1
r

ð2Þ

The force diagram of a wheel driving downhill a slope is shown in Fig. 3, where a is the slope angle. The resistance force FR balances the gravity component Wsina as given in Eq. (3), while the vertical force FV balances the gravity component Wcosa as shown in Eq. (4). The ratio of the resistance force to the vertical force was defined as resistance force coefficient kPRC [28], as shown in Eq. (5). If the resistance force can be estimated accurately, the slope angle that the planetary rover can safely traverse can be determined explicitly.

F R ¼ Wsina

ð3Þ

F V ¼ Wcosa

ð4Þ

kPRC ¼

FR Wsina ¼ ¼ tana FV Wcosa

ð5Þ

The normal stress increases from h1 and reaches the peak value at h = hm, which is considered as a function of the skid ratio [27,29]. It decreases to zero at h = h2. The shear stress is positive fromh1 toh0, but negative between h0 and h2 [27]. And h0 is denoted as the angular position of the transition point of the shear stress. The normal stress can be computed using [27,29]:

8   < r1 ðhÞ ¼ kbc þ ku r n ðcosh
cosh1 Þn ðhm  h  h1 Þ 
 in   h
: r2 ðhÞ ¼ kc þ ku r n cos h1
h
h2 ðh1
hm Þ
cosh1 ðh2  h < hm Þ b hm
h2

TB

W ω FP

r

z2

θ2

θ1

FR

ð6Þ

W

TB

v

r

v

FP θ1 θ0 τ θm τ

θ2

z2

z1

ω

σ

τσ

z1

σ (b) Normal and shear stress distribution

FV (a) Force diagram

Fig. 2. Force diagram and stress distribution under a driving wheel in skid.

Wsinα TB Wcosα

W ω

FV

v FR α

Fig. 3. Force diagram of a driving wheel downhill a slope.

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

43

where kc is the cohesive modulus of soil, ku is the frictional modulus of soil, and n is the sinkage exponent. These three pressure-sinkage characteristics parameters need to be measured by a series of pressure-sinkage experiments, and are used to characterize terrain’s stiffness [13,18,27]. The expression of shear displacement used to compute the shear stress in the region from h0 to h1 is different from that between h0 and h2 as shown in Eq. (7) [27,29], where c is the cohesion of soil, u is the internal friction angle, and K is the shear deformation modulus of soil. These three shearing deformation characteristics parameters need to be measured using a series of shearing experiments, and are used to characterize terrain’s shear strength [27].

8 ( > > > < ½c þ rðhÞtanu 1
exp

sðhÞ ¼

sinh
sinh0 sinh1
sinh
1
s  d d 1
h0 Þ


r½ðh1
hÞð1
s 1Þðh

!) h0  h  h1

K

n
o > > > : ½c þ rðhÞtanu 1
exp
r½ðh0
hÞ
ð1
sd Þðsinh0
sinhÞ h2  h  h0 K

ð7Þ

Eqs. (8) and (9) are used to determine hm and h0 using the entrance angle h1, respectively [27,29].

hm ¼ ðc1 þ c2 sd Þh1

ð8Þ

h0 ¼ ðd1 þ d2 sd Þh1

ð9Þ

The vertical force FV, resistance force FR, and braking torque TB can be computed by integrating the distributed normal and shear stress as follows:

Z F V ¼ rb Z F R ¼ rb

h1 h2

h1 h2

Z T B ¼ r2 b

rðhÞcoshdh þ rðhÞsinhdh þ

h0

h2

sðhÞdh


Z

h1

h0

Z

h1

h0

Z

h0

h2

sðhÞsinhdh


sðhÞcoshdh


Z

h0 h2

Z

h1

h0

sðhÞsinhdh

sðhÞcoshdh

 ð10Þ

 ð11Þ



sðhÞdh

ð12Þ

To reflect the skid-sinkage phenomenon, the measured sinkage exponent could be assumed not to be a constant, but a function of the skid ratio [27].

n ¼ n0 þ n1 sd þ n2 s2d

ð13Þ

In conclusion, the limitations of these traditional models mainly contain the following three aspects: (1) 16 parameters and variables, i.e., terrain mechanical parameters (kc, ku, n0, n1, n2, K, c, u), wheel running condition variables (v, sd, z1 used to compute h1), wheel-soil interaction parameters used to compute hm and h0 (c1, c2, d1, d2), and computational expensive soil shearing deformation (j), are needed in the modeling process. Terrain classical mechanical parameters are typically necessary to be measured by conducting two series of tests in advance using two special designed apparatuses, which is not suitable for planetary exploration missions. (2) In practice, the required wheel forward velocity and skid ratio are very difficult to be properly measured. (3) Due to the complexity of the normal and shear stress’ expressions, these integrated equations are not amenable to closed-form integrations and highly coupled, which means that it is impossible to identify all the terrain parameters using physically measurable quantities in real time. 3. Closed-form analytical wheel-soil interaction modeling To overcome the above mentioned constraints, the foundation of wheel-soil interaction modeling, i.e., the normal and shear stress equations, need to be simplified. 3.1. Linear normal stress As shown in Fig. 4, if the undisturbed soil surface is utilized as the reference of wheel sinakge, the wheel sinkage z can be computed using Eq. (14). According to this equation, the maximum wheel sinkage zmax occurs at h = 0, i.e., the lowest point of the wheel. As indicated by Bekker’s normal stress equation, the normal stress increases with the wheel sinkage, and thereby the maximum normal stress should also occur at the same location. However, the measured maximum normal stress occurs at the maximum radial stress angle hm, and the corresponding point is not the lowest point but almost the midway between h1 and h2 [30].

z ¼ r½cosðhÞ
cosðh1 Þ

ð14Þ

44

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

W

TB

ω

v

FP

r θ2

Traditional reference

θ1 θ

z

zσ New reference

zmax

σ

Fig. 4. Proposed reference of wheel sinkage used to compute normal stress.

To overcome this drawback, the connection of the entrance point and leaving point was proposed to be another reference of wheel sinkage, as shown in Fig. 4 [31]. To be differentiated from the traditional one, the wheel sinkage relative to this new reference is named as equivalent wheel sinkage and denoted by zr as given in Eq. (15).

  h1 þ h2 h1
h2 zr ðhÞ ¼ r cos h

cos 2 2

ð15Þ

The angular position of the maximum equivalent wheel sinkage zrm and the corresponding maximum normal stress angle hm occur midway between h1 and h2. The advantages of this new reference are: (1) the maximum radial stress occurs midway between h1 and h2 instead of the wheel’s lowest point; (2) the normal stress in the rear region (from h2 to hm) can be computed directly using the equivalent wheel sinkage instead of referring to the front region (between hm and h1). For sandy terrains, the soil rebound quantity is so small that the leaving angle h2 can be neglected [6,7,27,32]. In this paper, the leaving angle h2 was not considered for simplification, and Eq. (15) can be simplified as Eq. (16).

  h1 h1 zr ðhÞ ¼ r cos
h
cos 2 2

ð16Þ

When the wheel (radius r 160 mm  width b 165 mm) was running on sand HIT-LSS1 with the vertical load of 65 N under the skid ratio of 0.4, the equivalent wheel sinkage zr(h) and the normal stress r(h) computed using the traditional terramechanics models are shown in Fig. 5. According to Fig. 5, it was found that both zr(h) and r(h) can be well fitted with a quadratic equation. It is assumed that zr(h) and r(h) can be expressed using Eqs. (17) and (18), respectively. The quadratic fitted results are shown in Table 1.

 2 2 k12 k zr ðhÞ ¼ f 1 ðhÞ ¼ k11 h2 þ k12 h ¼ k11 h þ
12 2k11 4k11 

k22 2k21

2

2




k22 4k21

ð18Þ

×103

6 5 4 3 2

Computed data Fitted curve

1 0 -5

0

5 10 15 20 25 30 35

Wheel−soil contact angle θ (º)

Normal stress σ (Pa)

Wheel sinkage zσ (mm)

rðhÞ ¼ f 2 ðhÞ ¼ k21 h2 þ k22 h ¼ k21 h þ

ð17Þ

8 7 6 5 4 3 2 1 0 -1 -5

Estimated data Fitted curve 0

5 10 15 20 25 30 35

Wheel−soil contact angle θ (º)

Fig. 5. Computed and fitted equivalent wheel sinkage and estimated and fitted normal stress.

45

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54 Table 1 Quadratic fitted results of the equivalent wheel sinkage and normal stress.

r(h)

zr(h) R-square

k11(10
2)

Std Err

k12(10
1)

Std Err

R-square

k21

Std Err

k22

Std Err

0.9999


2.42

9.01  10
6

7.62

2.20  10
4

0.9986


29.50

0.14

924.80

3.47

Besides the above mentioned quadratic approximation, the Taylor series expansion of the equivalent wheel sinkage at zero equilibrium was also used, as given by Eq. (19). The terms higher than second-order were ignored. When the entrance angle is equal to 60°, the comparison of the approximated and computed equivalent wheel sinkage and its relative error are shown in Fig. 6. According to this figure, the relative error of the Tylor approximation increases dramatically with the wheelsoil contact angle. The maximum relative error of the Tylor approximation is around 37%, while that of the quadratic approximation is about 0.50%. When the entrance angle is less than 30°, the maximum relative error of the Tylor approximation can be captured within 7%.

 

    h1 h1 h1 1 h1 2 1 h1 3 zr ðhÞ ¼ r cos
h
cos ¼ r sin h
cos h
sin h þ Oðh4 Þ 2 6 2 2 2 2 2

ð19Þ

In conclusion, if the entrance angle is less than 30°, both the quadratic approximation and Tylor series expansion can be used to approximate the equivalent wheel sinkage. However, if the entrance angle is larger than 30°, the quadratic approximation would perform better. Because both zrmax and rmax occur at h = hm, r(h) can be expressed linearly using zr(h) with a constant Kr named as terrain stiffness:

rðhÞ ¼

k21 zr ðhÞ ¼ K r zr ðhÞ k11

ð20Þ

The linearity between the normal stress and equivalent wheel sinkage zr has been demonstrated by the experimental data obtained using two wheels and on two types of planetary soil simulants. Moreover, the proposed terrain stiffness has also been used to compare the stiffness of various sandy terrains, i.e., it can be used to characterize sandy terrains’ pressure-sinkage mechanical property [31]. 3.2. Linear shear stress Even with the linear normal stress, closed-formed wheel-soil interaction equations still cannot be achieved due to the complexity of Janosi’s exponential shear stress equation. This motivates the simplification of the shear stress. The shear stress distribution has been measured and reported in [30], and it is found that the ratio of the negative area to the positive is so large that the positive portion in the front region is neglected in this paper. And it is assumed that s(h) is negative along the entire wheel-soil interface and the maximum shear stress smax occurs also at h = hm. In conclusion, the supposed stress distribution can be illustrated using Fig. 7. Similar to the linear normal stress, the simplified shear stress can also be expressed linearly using zr(h) with a constant Ks named as shear strength in this paper, as shown in Eq. (21).

sðhÞ ¼
K s zr ðhÞ

ð21Þ

To verify the accuracy of the simplified shear stress using physically measurable variables (FV, FR, TB, h1), the linear normal and shear stress is integrated along the wheel-soil interface, as detailed in the following subsection.

20 15 10 5 0 -50

Computed Tylor approximation Quadratic approximation

10 20 30 40 50 60 Wheel-soil contact angle (°) (a) Comparison of computed and approximated zσ(θ)

Relative error (%)

Equivalent sinkage (mm)

25

40 35 30 25 20 15 10 5 0 0

Tylor approximation Quadratic approximation

10 20 30 40 50 60 Wheel-soil contact angle (°) (b) Relative error

Fig. 6. Comparison of the approximated and computed equivalent wheel sinkage and its relative error.

46

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

TB

W ω

r

v

FP θ1

θm

τ

τ σ zσ

z σ

Fig. 7. Assumed stress distribution under a driving wheel in skid.

3.3. Analytical modeling Substituting Eqs. (20) and (21) into Eqs. (10)–(12) yields the analytical forms of the force balance equations. The contribution of r(h) to FV can be determined using Eq. (22):

Z F Vr ¼ K r br

h1

2



0

"    2 # h1 h1 h1 h1 h1 2 h1 cos
h
cos coshdh ¼ K r br cos
sin cos 2 2 2 2 2 2

ð22Þ

The shear stress contributes to the vertical force negatively as:

Z F Vs ¼
K s br2



h1

0

"    2 # h1 h1 h1 h1 h1 h1 cos
h
cos sinhdh ¼
K s br2 sin
cos sin 2 2 2 2 2 2

ð23Þ

The analytical expression of vertical force is given by:

" "  2 #  2 # h1 h1 h1 h1 h1 h1 h1 2 h1
K s br F V ¼ F Vr þ F Vs ¼ K r br cos
sin cos sin
cos sin 2 2 2 2 2 2 2 2 2

ð24Þ

The shear stress contributes to the resistance force as follows:

Z F Rs ¼ K s br2

h1

0



"    2 # h1 h1 h1 h1 h1 h1 cos
h
cos coshdh ¼ K s br 2 cos
sin cos 2 2 2 2 2 2

ð25Þ

The resistance force integrated from the normal stress is given as:

Z F Rr ¼ K r br

2

h1

0



"    2 # h1 h1 h1 h1 h1 2 h1 cos
h
cos sinhdh ¼ K r br sin
cos sin 2 2 2 2 2 2

ð26Þ

The analytical expression of resistance force can be obtained as:

" F R ¼ F Rs þ F Rr ¼ K s br2

"  2 #  2 # h1 h1 h1 h1 h1 h1 h1 h1 þ K r br2 cos
sin cos sin
cos sin 2 2 2 2 2 2 2 2

ð27Þ

The analytical braking torque equation is given as:

Z T B ¼ K s br3

0

h1

   h1 h1 h1 h1 cos
h
cos dh ¼ K s br 3 2sin
h1 cos 2 2 2 2

ð28Þ

In contrast with the classical terramechanics equations, these analytical expressions are functions of only the proposed terrain stiffness Kr, shear strength Ks, and the entrance angle h1. The running state variables (i.e., skid ratio sd and wheel forward velocity v), classical terrain mechanical parameters (kc, ku, n0, n1, n2, K, c, u), coefficients (c1, c2, d1, d2) used to compute hm and h0, and computationally expensive soil shear displacement (j) are no longer required.

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

47

4. Experimental validation The vertical load can be computed from a quasi-static force analysis of the rover, with the knowledge of the rover configuration and its mass distribution. The braking torque can be estimated from the current input to the motor. Therefore, to validate the accuracy of these analytical expressions, the system of Eqs. (24) and (28) was first used to compute Ks and Kr as follows:

Kr ¼ Ks ¼

FV þ

T B k2 r k3

ð29Þ

br2 k1 TB br 3 k3

ð30Þ

where

k1 ¼

 2 h1 h1 h1 h1 cos
sin cos 2 2 2 2

k2 ¼

 2 h1 h1 h1 h1 sin
cos sin 2 2 2 2

k3 ¼ 2sin

h1 h1
h1 cos 2 2

Then, the estimated Kr and Ks were used to calculate FR (see Eq. (31)) to compare with the experimental data. If the WMR was running on hard terrain, the resistance force expression can be simplified further into Eq. (32) [33]. 2

2

FR ¼

T B k1 þ k2 k2 þ FV r k1 k3 k1

ð31Þ

FR ¼

TB r

ð32Þ

4.1. Experimental set-up To measure FV, FR, TB, and z, a wheel performance test-bed (Fig. 8) was developed at the Harbin Institute of Technology. This test-bed contains a driving motor and a carriage motor, and the skid ratio was changed by actively controlling the driving and carriage motors simultaneously. FV and FR were measured using a force/torque sensor, while TB was measured by a torque sensor on the axle of the driving motor. The wheel vertical sinkage z was measured by a linear displacement sensor. Reference [34] describes the details of experimental procedure. The two experimental wheels used in this research were named as WH1 and WH2 (see Fig. 9) for convenience, and their parameters are listed in Table 2, where hL and nL denote lug height and lug number, respectively. Impurities were removed from soft sand, which was then sieved, ventilated, and dried to produce two types of experimental soil HIT-LSS1 and HIT-LSS2 to simulate planetary regolith. HIT-LSS1 was provided by the Shanghai Academy of

Fig. 8. Photograph of wheel performance test-bed.

48

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

(a) Wheel WH1

(b) Wheel WH2

Fig. 9. Photograph of experimental wheels.

Table 2 Parameters of experimental wheels. Wheel code

r (mm)

b (mm)

hL (mm)

nL

WH1 WH2

140.5 135

150 165

15 10

28 24

Table 3 Classical mechanical parameters of HIT-LSS1 and HIT-LSS2 [26]. Sand

kc (kPa/mn
1)

ku (kPa/mn)

n

c (kPa)

u (°)

K (mm)

HIT-LSS1 HIT-LSS2

15.6
20.7

2407.4 1594.8

1.10 0.79

0.25 0.46

31.9 38.1

9.7–13.1 4.0–16.6

Table 4 Experimental setup. Wheel code

v (mm/s)

W (N)

sd

WH1 WH2

10 10

110 85

0.1, 0.2, 0.3, 0.4, 0.5, 0.6

Spaceflight Technology and HIT-LSS2 was made from commercial standard sand. Their classical mechanical parameters are shown in Table 3 [26]. The experimental setup is illustrated in Table 4. Note that, when sd = 0, due to the influence of wheel grousers, it was found that the wheels were in slip. And the smallest experimental skid ratio was set as 0.1. 4.2. Quasi-static validation After the wheels entering steady state, the measured wheel-soil interactions were used to validate the accuracy of these analytical equations. The estimated mechanical parameters of HIT-LSS1 using WH1, and those of HIT-LSS2 using WH2 are shown in Fig. 10. According to this figure, the terrain stiffness firstly increases with the increase of skid ratio, then reaches a peak value, finally suffers a decrease. Like the terrain stiffness, the shear strength increases with the increase of skid ratio, but then suffers a slight decrease. To verify the accuracy of the estimated terrain parameters, the resistance force that was calculated by these soil parameters was used to validate against its experimental data. Fig. 11 shows the comparison between the measured and calculated resistance force, and its relative error. When wheel WH1 was running on sand HIT-LSS1, these analytical models can capture the maximum relative error within 4.26%. As for WH2 and HIT-LSS2, the maximum relative error is under 6.37%. Though the error is large in a relative sense when sd = 0.5, the corresponding absolute error of the estimated FR is 2.21 N (2.60% of the vertical force). 4.3. Comparison of analytical and traditional stresses According to the experimental validation results, the established closed-form analytical models has good performance. The base of these analytical equations is the simplification of the normal and shear stress equations. To further prove the

49

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

Mechanical parameters (Pa/m)

Mechanical parameters (Pa/m)

6

16 ×10 14 12 10 8 Stiffness 6 Shear strength 4 2 0 -2 0.1 0.2 0.3 0.4 0.5 Skid ratio sd (a) HIT-LSS1 using WH1

0.6

7 2.5 ×10

2 1.5 Stiffness Shear strength

1 0.5 0 -0.5 0.1

0.2

0.3 0.4 0.5 Skid ratio sd (b) HIT-LSS2 using WH2

0.6

Fig. 10. Computed terrain equivalent parameters after wheels entering steady state.

4.5

Measured Computed

50

Relative error (%)

Resistance force FR(N)

60

40 30 20 10 0

0.1

0.2

0.3 0.4 Skid ratio sd

0.5

4 3.5 3 2.5 2

0.6

0.1

0.2

0.3 0.4 Skid ratio sd

0.5

0.6

0.2

0.3 0.4 Skid ratio sd

0.5

0.6

(a) WH1 & HIT-LSS1 7 6

Measured Computed

30

Relative error (%)

Resistance force FR(N)

35

25 20 15 10

5 4 3 2 1

0.1

0.2

0.3 0.4 Skid ratio sd

0.5

0.6

0

0.1

(b) WH2 & HIT-LSS2 Fig. 11. Comparison of the measured and computed resistance force and its relative errors.

feasibility of this kind of simplification, the comparison of the normal and shear stresses computed using the analytical models and those by the traditional theories is shown in Fig. 12. The wheel-soil interactions used to estimate the normal and shear stresses were obtained by the wheel (radius r 160 mm  width b 165 mm) and the planetary soil simulant HIT-LSS1. Note that ‘‘Ana.” and ‘‘Tra.” in the legend stand for ‘‘Analytical models” and ‘‘Traditional models”, respectively, and the normal and shear stresses computed by the traditional models are plotted using solid lines. When the wheel-soil contact angle is 0, the normal stress is 0. However, due to the influence of the soil cohesion, the shear stress estimated using the traditional terramechanics theories is not 0, as given by Eq. (7). The quadratic approximation has little influence on the normal stress. When the skid ratio is equal to 0.6, the maximum relative error of the normal stress computed by the analytical models is less than 10% against the normal stress by the traditional models. While this approximation has some effect on the shear stress’s magnitude. Under various skid ratios, the shear stresses’ magnitudes obtained by the analytical models are smaller than those by the traditional models. However, for a specific skid ratio, the integrations of the shear stresses along the wheel-soil interface are the same.

50

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

Normal and shear stress (Pa)

10

×103 σ(θ),Ana.,sd=0.1 τ(θ),Ana.,sd=0.1 σ(θ),Ana.,sd=0.4 τ(θ),Ana.,sd=0.4 σ(θ),Ana.,sd=0.6 τ(θ),Ana.,sd=0.6

σ(θ),Tra.,sd=0.4

σ(θ),Tra.,sd=0.1

8 6

σ(θ),Tra.,sd=0.6

4 2

τ(θ),Tra.,sd=0.1

0 -2

τ(θ),Tra.,sd=0.6

τ(θ),Tra.,sd=0.4

-4 -5

0

5

10

15

20

25

30

35

40

Wheel-soil contact angle (°) Fig. 12. Stress comparison obtained using analytical and traditional models.

5. Online applications These validated analytical terramechanics models are then applied to estimate terrain equivalent parameters and resistance force in real time. To estimate Kr and Ks online, Eqs. (29) and (30) must be rewritten in matrix form:



Kr Ks

"

 ¼

1 br 2 k1

k2 1 br3 k3 k1

0

1 br3 k3

#

FV

 ð33Þ

TB

where

k1 ¼

 2 h1 ðtÞ h1 ðtÞ h1 ðtÞ h1 ðtÞ cos
sin cos 2 2 2 2

k2 ¼

 2 h1 ðtÞ h1 ðtÞ h1 ðtÞ h1 ðtÞ sin
cos sin 2 2 2 2

k3 ¼ 2sin

h1 ðtÞ h1 ðtÞ
h1 ðtÞcos 2 2

Further, the resistance force can be estimated using the matrix of terrain stiffness and shear strength as:

" F R ¼ ½ k2

k1 

br 2 K r br2 K s

#

" ¼

k2 k1

2

1 r

ðk1 Þ þ ðk2 Þ k1 k3

2

#

FV TB

 ð34Þ

When the wheel WH2 was running on sand HIT-LSS1 with the vertical load of 85 N under the skid ratio of 0.3, the experimental data was selected. The computed Kr and Ks is shown in Fig. 13 (a), and the comparison of the estimated and measured FR is illustrated in Fig. 13 (b). The terrain stiffness increases with time, then reaches a constant, and the shear strength also shears the same trend. Not only the measured and estimated resistance force, but also the estimated terrain stiffness fluctuates in the form of sinusoid, partly due to the difficulty in preparing a perfectly uniform soil for testing. More importantly, it is caused by the oscillation occurring in the wheel-soil interaction. Irani et al. proposed a dynamic model to quantitatively describe this phenomenon, as expressed in

f O ðtÞ ¼ AsinðxO ðsd Þt þ /Þ ¼ Asinððks sd þ CÞt þ /Þ

ð35Þ

where A is the amplitude of oscillation, xO is the frequency of oscillation, sd is the skid ratio, / is a phase shift that can be used to fit the simulation predictions to the experimental data, ks is a constant of proportionality, and C is the intercept when xO is plotted versus sd [20,35]. The amplitude of oscillation was found to be a function of the soil density, skid ratio, and vertical load. It increases with both the vertical load and the skid ratio [3,35]. However, dry sand is a typical dilative soil and has a limited range of volume fractions (ratio of particle to air volume). Hence, the wheel-soil interaction has a limited effect on the sand density, which means that dry sand’s density has a limited effect on the amplitude of oscillation [36].

51

×107

5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0

40

Resistance force FR (N)

Mechanical parameters (Pa/m)

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

Sand stiffness Shear strength

35 30 25 20 15 5 0

10

20

30

40

Time t (s)

50

60

Measured Computed

10

0

10

20

(a) Kσ and Kτ

30

40

Time t (s) (b) FR

50

60

Fig. 13. Online estimated terrain equivalent parameters and resistance force.

35 30

Relative error (%)

25 20 15 10 5 0

0

10

20

30

40

50

60

Time t (s) Fig. 14. Relative errors of the online estimated resistance force.

The relative errors of the estimated FR are shown in Fig. 14. In the beginning 20 s, the relative error decreases from about 18% to around 5%, and then it fluctuates around 5%. The amplitude of the estimated resistance force is much larger than that of the measured especially in the beginning 10 s, which can be explained by considering the wheel vertical velocity derived from the wheel vertical sinkage. The derived vertical velocity and acceleration are shown in Fig. 15. In the beginning 10 s, the wheel vertical velocity decreases from about 1.3 mm/s to
0.5 mm/s. The measured wheel vertical sinkage zM includes the sinkage zF caused by the vertical force and the sinkage zV by the vertical velocity (see Eq. (36)). And the measured wheel sinkage must be larger than that caused by the vertical force, which means that the computed entrance angle using the measured wheel vertical sinkage is larger than that caused by the vertical force.

zM ¼ zF þ zV

ð36Þ

The definitions of K1 and K2 are shown in Eqs. (37) and (38), respectively.

K1 ¼

k2 k1

K2 ¼

ðk1 Þ þ ðk2 Þ k1 k3

ð37Þ 2

2

ð38Þ

As shown in Eq. (34), the resistance force increases with both the K1 and the K2. Furthermore, these two coefficients increase with the entrance angle h1 as shown in Fig. 16. In conclusion, the larger the entrance angle, the larger the estimated resistance force. Therefore, the computed resistance force is larger than the measured in the beginning 10 s.

52

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

Vertical acceleration (mm/s2)

Vertical velocity (mm/s)

1.5 1 0.5 0 -0.5 -1

0

10

20

30

40

50

Time t (s)

60

2 1 0 -1 -2 -3 -4

0

10

20

30

40

50

60

Time t (s)

(a) Vertical velocity

(b) Vertical acceleration

Fig. 15. Derived wheel vertical velocity and acceleration using wheel vertical sinkage.

1.4 1.2

K1 and K2

1

K1 K2

0.8 0.6 0.4 0.2 0

0

0.2

0.4

0.6

0.8

Entrance angle θ1(rad) Fig. 16. Curves of K1 and K2 versus the entrance angle.

Measured vertical force FV (N)

110 100 90 80 70 60 50 40

0

10

20

30

40

50

60

Time t (s) Fig. 17. Measured vertical force.

In addition, the maximum vertical acceleration is about 4.0  10
3 m/s2, and the wheel vertical load is 85 N. The maximum unbalance force in the vertical direction is 0.034 N and can be neglected. The relative vertical velocity between the wheel and soil results in a vertical damp force. The force balance equation in the vertical direction can be written as Eq. (39), where FV is the vertical force, FD is the damp force caused by wheel’s vertical velocity.

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

W ¼ FV þ FD

53

ð39Þ

Before the wheels entering steady state, the damp force FD plus the measured vertical force FV equal the wheel vertical load W, which means that the measured vertical force FV during this period is much smaller than the experimental set vertical load W. After the wheels entering steady state, the measured vertical force FV can be approximated as wheel’s vertical load W. Therefore, the measured vertical force FV first increases, and then fluctuates around a constant, which can be demonstrated using the measured vertical force as shown in Fig. 17. 6. Conclusions and discussions The following summarizes the findings of this research: (1) Bekker’s normal stress and Janosi’s shear stress equations, which are the basis of conventional terramechanics theories for wheeled mobile robots, are so complex that the wheel-soil interaction force/torque equations are not amenable to closed form integrations, and it is impossible to estimate all the terrain classical parameters and resistance force online. (2) To simplify wheel-soil interaction models, the positive portion of shear stress is neglected and it is assumed to be negative along the entire wheel-soil interface and its maximum value occurs at h = hm. Both the normal stress and the simplified shear stress can be expressed linearly using the equivalent wheel sinkage relative to the connection of the entrance and leaving points with a constant named as equivalent terrain stiffness and shear strength, respectively. (3) The integrated analytical equations are validated using experimental data obtained with two wheels and on two sand types. Both the equivalent terrain stiffness and shear strength first increase with the skid ratio, and then reach a peak value, finally, suffer a decrease. The maximum relative error of the estimated resistance force is less than 7%, and the corresponding absolute error is 2.21 N (2.60% of the vertical force). (4) These analytical models can be applied to estimate terrain’s stiffness, shear strength, and resistance force in real time. Sandy terrain’s stiffness and shear strength increase first, and then reach a constant. The online estimated resistance force’s relative error is much higher in the beginning process, which is caused by wheel’s vertical velocity. The proposed online estimation approach was implemented on a single wheel performance test-bed, where the wheel vertical sinkage z was measured by a linear displacement sensor. On the test-bed, it is easy to find a fixed reference point to install the linear displacement sensor, which cannot be applied to a practical mobile robot. As for WMRs, the wheel sinkage z can be computed with vision-based techniques [37,38] or by kinematic analysis of the rover suspension [39]. 7. Funding This study was supported in part by the National Natural Science Foundation of China [grant number 61370033], Foundation for Innovative Research Groups of the National Natural Science Foundation of China [grant number 51521003], Foundation of Chinese State Key Laboratory of Robotics and Systems [grant number SKLRS201707A], Fundamental Research Funds for Central Universities [grant number HIT.NSRIF.2017027], and the ‘‘111 Project” [grant number B07018]. References [1] A.N. Ibrahim, S. Aoshima, N. Shiroma, Y. Fukuoka, The effect of assistive anchor-like grousers on wheeled rover performance over unconsolidated sandy dune inclines, Sensors 16 (2016) 1507. [2] M. Sutoh, J. Yusa, T. Ito, K. Nagatani, K. Yoshida, Traveling performance evaluation of planetary rovers on loose soil, J. Field Rob. 29 (2012) 648–662. [3] J. Guo, L. Ding, H. Gao, T. Guo, G. Liu, H. Peng, An apparatus to measure wheel-soil interactions on sandy terrains, IEEE/ASME Trans. Mechatron. 23 (2018) 352–363. [4] Y. Yang, Y. Sun, S.G. Ma, Paddle trajectory generation for accessing soft terrain by an ePaddle locomotion mechanism, The IEEE Int. Conf. Rob. Autom. (2013) 403–408. [5] Y. Li, L. Ding, Z. Zheng, Q. Yang, X. Zhao, G. Liu, A multi-mode real-time terrain parameter estimation method for wheeled motion control of mobile robots, Mech. Syst. Sig. Process. 104 (2018) 758–775. [6] K. Iagnemma, S. Kang, H. Shibly, S. Dubowsky, Online terrain parameter estimation for wheeled mobile robots with application to planetary rovers, Ieee T Robot 20 (2004) 921–927. [7] K. Iagnemma, S. Dubowsky, Traction control of wheeled robotic vehicles in rough terrain with application to planetary rovers, Int. J. Rob. Res. 23 (2004) 1029–1040. [8] Z. Zhao, D. Lei, J. Chen, H. Li, Optimal control of mode transition for four-wheel-drive hybrid electric vehicle with dry dual-clutch transmission, Mech. Syst. Sig. Process. 105 (2018) 68–89. [9] X. Zeng, G. Li, G. Yin, D. Song, S. Li, N. Yang, Model predictive control-based dynamic coordinate strategy for hydraulic hub-motor auxiliary system of a heavy commercial vehicle, Mech. Syst. Sig. Process. 101 (2018) 97–120. [10] R. Wang, C. Hu, Z. Wang, F. Yan, N. Chen, Integrated optimal dynamics control of 4WD4WS electric ground vehicle with tire-road frictional coefficient estimation, Mech. Syst. Sig. Process. 60–61 (2015) 727–741. [11] Y.H. Liu, T. Li, Y.Y. Yang, X.W. Ji, J. Wu, Estimation of tire-road friction coefficient based on combined APF-IEKF and iteration algorithm, Mech. Syst. Sig. Process. 88 (2017) 25–35. [12] W. Yu, O.Y. Chuy Jr, E.G. Collins Jr, P. Hollis, Analysis and experimental verification for dynamic modeling of a skid-steered wheeled vehicle, Ieee T Robot 26 (2010) 340–353. [13] K. Iagnemma, P. Martinet, D. Wang, Introduction: vehicle–terrain interaction for mobile robots, J. Field Rob. 27 (2010) 105–106. [14] Y. Yang, Y. Sun, S. Ma, Drawbar pull of a wheel with an actively actuated lug on sandy terrain, J. Terramech. 56 (2014) 17–24.

54

Z. Liu et al. / Mechanical Systems and Signal Processing 119 (2019) 39–54

[15] Xinhua News, CNSA announces the configuration of the first Mars explorer and rover of China. http://zhuanti.spacechina.com/n1411849/ c1412769/content.html, 2016 (accessed 02 October 2016). [16] D. Brown, NASA announces Mars 2020 rover payload to explore the red planet as never before. http://www.nasa.gov/press/2014/july/nasa-announcesmars-2020-rover-payload-to-explore-the-red-planet-as-never-before, 2014 (accessed 02 October 2016). [17] ESA, The ExoMars programme 2016-2020. http://exploration.esa.int/mars/46048-programme-overview/, 2016 (accessed 20 October 2016). [18] M.G. Bekker, Off-the-Road Locomotion: Research and Development in Terramechanics, University of Michigan Press, Ann Arbor, 1960. [19] J. Wong, A. Reece, Prediction of rigid wheel performance based on the analysis of soil-wheel stresses part I. Performance of driven rigid wheels, J. Terramech. 4 (1967) 81–98. [20] R. Irani, R. Bauer, A. Warkentin, A dynamic terramechanic model for small lightweight vehicles with rigid wheels and grousers operating in sandy soil, J. Terramech. 48 (2011) 307–318. [21] M. Lyasko, Slip sinkage effect in soil–vehicle mechanics, J. Terramech. 47 (2010) 21–31. [22] K. Iagnemma, C. Senatore, B. Trease, R. Arvidson, K. Bennett, A. Shaw, F. Zhou, L.V. Dyke, K. Iagnemma, C. Senatore, Terramechanics modeling of mars surface exploration rovers for simulation and parameter estimation, Am. Soc. Mech. Eng. (2011) 805–812. [23] M. Apfelbeck, S. Kuß, B. Rebele, B. Schäfer, A systematic approach to reliably characterize soils based on Bevameter testing, J. Terramech. 48 (2011) 360–371. [24] H. Moore, R. Hutton, R. Scott, C. Spitzer, R. Shorthill, Surface materials of the viking landing sites, J. Geophys. Res. 82 (1977) 27. [25] J. Crisp, J.R. Matijevic, Characterization of the Martian surface deposits by the mars pathfinder rover, sojourner, Science 278 (1997) 1765–1768. [26] L. Ding, H. Gao, Y. Li, G. Liu, Z. Deng, Improved explicit-form equations for estimating dynamic wheel sinkage and compaction resistance on deformable terrain, Mech. Mach. Theory 86 (2015) 235–264. [27] H. Gao, J. Guo, L. Ding, N. Li, Z. Liu, G. Liu, Z. Deng, Longitudinal skid model for wheels of planetary exploration rovers based on terramechanics, J. Terramech. 50 (2013) 327–343. [28] L. Ding, J. Guo, B. Yan, H. Gao, Z. Liu, N. Li, H. Yang, Z. Deng, Longitudinal skid experimental investigation for wheels of planetary exploration rovers, J. Mech. Eng. 51 (2015) 45–53. [29] J. Wong, A. Reece, Prediction of rigid wheel performance based on the analysis of soil-wheel stresses: Part II. Performance of towed rigid wheels, J. Terramech. 4 (1967) 7–25. [30] C. Senatore, K. Iagnemma, Analysis of stress distributions under lightweight wheeled vehicles, J. Terramech. 51 (2014) 1–17. [31] J. Guo, H. Gao, L. Ding, T. Guo, Z. Deng, Linear normal stress under a wheel in skid for wheeled mobile robots running on sandy terrain, J. Terramech. 70 (2017) 49–57. [32] K. Nagatani, A. Ikeda, K. Sato, K. Yoshida, Accurate estimation of drawbar pull of wheeled mobile robots traversing sandy terrain using built-in force sensor array wheel, IEEE/RSJ Int. Conf. Intell. Robots Syst. (2009) 2373–2378. [33] H. Zhang, W. Zhao, Decoupling control of steering and driving system for in-wheel-motor-drive electric vehicle, Mech. Syst. Sig. Process. 101 (2018) 389–404. [34] L. Ding, H. Gao, Z. Deng, K. Nagatani, K. Yoshida, Experimental study and analysis on driving wheels’ performance for planetary exploration rovers moving in deformable soil, J. Terramech. 48 (2011) 27–45. [35] R.A. Irani, R.J. Bauer, A. Warkentin, Modelling a single-wheel testbed for planetary rover applications, ASME 2010 Dyn. Syst. Control Conf. (2010) 181– 188. [36] G. Meirion-Griffith, C. Nie, M. Spenko, Development and experimental validation of an improved pressure-sinkage model for small-wheeled vehicles on dilative, deformable terrain, J. Terramech. 51 (2014) 19–29. [37] G. Reina, L. Ojeda, A. Milella, J. Borenstein, Wheel slippage and sinkage detection for planetary rovers, IEEE/ASME Trans. Mechatron. 11 (2006) 185– 195. [38] G.M. Hegde, C. Ye, C.A. Robinson, A. Stroupe, Computer-vision-based wheel sinkage estimation for robot navigation on lunar terrain, IEEE/ASME Trans. Mechatron. 18 (2013) 1346–1356. [39] B. Wilcox, Non-geometric hazard detection for a Mars microrover, AIAA Conf. Intell. Rob. Field, Factory, Service, Space (1994) 675–684.