Identifying mechanical property parameters of planetary soil using in-situ data obtained from exploration rovers

Planetary and Space Science 119 (2015) 121–136

Contents lists available at ScienceDirect

Planetary and Space Science journal homepage: www.elsevier.com/locate/pss

Identifying mechanical property parameters of planetary soil using in-situ data obtained from exploration rovers Liang Ding a,b, Haibo Gao a, Zhen Liu a, Zongquan Deng a, Guangjun Liu b a b

The State Key Laboratory of Robotics and System, Harbin Institute of Technology, Harbin 150001, Heilongjiang, China Department of Aerospace Engineering, Ryerson University, 350 Victoria Street, Toronto, Ontario, Canada M5B 2K3

art ic l e i nf o

a b s t r a c t

Article history: Received 11 February 2014 Received in revised form 8 May 2015 Accepted 2 September 2015 Available online 25 September 2015

Identifying the mechanical property parameters of planetary soil based on terramechanics models using in-situ data obtained from autonomous planetary exploration rovers is both an important scientific goal and essential for control strategy optimization and high-fidelity simulations of rovers. However, identifying all the terrain parameters is a challenging task because of the nonlinear and coupling nature of the involved functions. Three parameter identification methods are presented in this paper to serve different purposes based on an improved terramechanics model that takes into account the effects of slip, wheel lugs, etc. Parameter sensitivity and coupling of the equations are analyzed, and the parameters are grouped according to their sensitivity to the normal force, resistance moment and drawbar pull. An iterative identification method using the original integral model is developed first. In order to realize real-time identification, the model is then simplified by linearizing the normal and shearing stresses to derive decoupled closed-form analytical equations. Each equation contains one or two groups of soil parameters, making step-by-step identification of all the unknowns feasible. Experiments were performed using six different types of single-wheels as well as a four-wheeled rover moving on planetary soil simulant. All the unknown model parameters were identified using the measured data and compared with the values obtained by conventional experiments. It is verified that the proposed iterative identification method provides improved accuracy, making it suitable for scientific studies of soil properties, whereas the step-by-step identification methods based on simplified models require less calculation time, making them more suitable for real-time applications. The models have less than 10% margin of error comparing with the measured results when predicting the interaction forces and moments using the corresponding identified parameters. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Planetary exploration rovers Mechanical property Parameter estimation Wheel–terrain interaction Terramechanics

1. Introduction Scientists have long been interested in the mechanical properties of planetary soil, attempting to improve our scientific knowledge of the planets' geological properties and to provide the engineering knowledge required for planetary surface exploration or future settlement activities (Arvidson et al., 2004a; Perko et al., 2004; Carrier et al., 1991). Ever since NASA's Sojourner rover landed on Mars in 1997, there has been an upsurge in planetary exploration using autonomous wheeled rovers. The Spirit and Opportunity rovers endured several years of activity on Mars and made many significant discoveries (NASA JPL, 2013). The recently landed Curiosity rover of the Mars Science Laboratory (MSL) mission (NASA, 2013), and E-mail addresses: [email protected] (L. Ding), [email protected] (H. Gao), [email protected] (Z. Liu), [email protected] (Z. Deng), [email protected] (G. Liu). http://dx.doi.org/10.1016/j.pss.2015.09.003 0032-0633/& 2015 Elsevier Ltd. All rights reserved.

future rovers such as those in the projects of ExoMars (ESA, 2013) and SELENE-2 (Japan Aerospace Exploration Agency, 2013) are required to traverse more challenging terrains. To enhance a rover’s mobility and fulfill its exploration missions, the design of its control and planning systems needs to take into account both the physical characteristics of the rover and the mechanical properties of the soil that can be obtained through parameter identification or estimation (Iagnemma et al., 2004). The interaction mechanics model involves many soil property parameters. Identifying the model parameters of wheel–terrain interaction using in-situ data obtained by on-board sensors is important for both soil investigation and capacity enhancement of a rover. For example, the MSL mission of NASA was planned to use the Curiosity rover as a terramechanics instrument to learn about Martian soils. The operational purpose is to help Curiosity rover chart the best path across Martian terrain, avoiding sand traps and slippery slopes, and the scientific purpose is to study the soil properties (Lutz, 2013).

122

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

The mechanical properties of soils can be investigated using insitu experiments. However, the standard pressure–sinkage experiments and shearing experiments in canonical terramechanics (Bekker, 1969; Wong, 2010) are difficult to implement in unmanned missions. Interaction mechanics experiments between the exploration rovers and terrains are usually performed for in-situ parameter identification. Researchers from the Sojourner rover team conducted experiments by driving one wheel of the rover while keeping the other wheels stationary. The soil appeared to show little or no cohesion, and the friction angles were found to be between 32° and 41° (The Rover Team, 1997). Wheel–terrain interaction experiments were also performed during NASA's Mars Exploration Rover missions (Arvidson et al., 2003). The Spirit and Opportunity rovers were used to study the soil properties of Meridiani Planum (Arvidson et al., 2004a) and Gusev Crater (Arvidson et al., 2004b), respectively, by excavating the subsurface soil with their wheels for in-situ observations. Besides the off-line methods, on-line parameter identification methods are required to enhance the mobility of rovers. Shibly et al. (2005) linearized the Wong–Reece normal stress equation and Janosi shear stress equation to obtain a closed-form formula. Then, a linear least-squares estimator was applied to estimate the internal cohesion and friction angle online while setting the shearing deformation modulus to typical values (Iagnemma et al., 2004). Hutangkabodee et al. (2006) developed a method to identify the internal friction angle, shearing deformation modulus, and lumped pressure sinkage coefficient, while the internal cohesion was fixed at 3 kPa. The composite Simpson’s rule was employed to obtain an approximated model. Ray (2009) proposed a multiplemodel approach to estimate the terrain parameters for rigid-wheel rovers based on the Wong–Reece model. Setterfield and Ellery (2013) presented an approach of estimating terrain responses of a rover with rocker-bogie mobility system using on-board instruments. Cross et al. developed a method of estimating terrain parameters including internal friction angle and internal cohesion for a rigid wheeled rover using neural networks (Cross et al., 2013). The existing approaches can only identify some of the terrain parameters. Identifying all the terrain parameters in the interaction mechanics model is challenging because of the nonlinear and coupling nature of the functions. The main contributions of this paper lie in identifying all the terrain parameters for a high-fidelity interaction mechanics model that incorporates the effects of slip– sinkage, wheel lugs, dimensions and payloads, including an indepth model analysis of the coupling and parameter sensitivity, development of an iterative parameter identification method using the original integral model to enhance the accuracy of parameters that could be used for scientific soil research, development of a novel decoupled closed-form model as functions of the limited number of unknowns and the observed contact information, and realization of high speed parameter identification to support real-time control and simulation of wheeled rovers. The remainder of this paper is organized as follows. Section II introduces an integral form wheel–terrain interaction mechanics model with high accuracy that considers multiple effects. Section III presents an experimental investigation of the wheel–terrain interaction mechanics. Section IV analyzes the coupling degree and parameter sensitivity of the models. Section V proposes an iterative parameter identification method based on the original integral model. Section VI derives the closed-form decoupled analytical models, and parameter identification methods based on these models are presented in Section VII. Finally, the conclusions are drawn in Section XIII.

2. Integral-form wheel–terrain interaction mechanics model (M1) Planetary rovers are usually equipped with lugs of a certain height to improve their tractive ability in deformable soil. The free-body diagram of a lugged wheel on deformable terrain is shown in Fig. 1 (Ding et al., 2014a), where TD is the driving torque generated by the driving motor and gears; W and fDP are the vertical load and resistance force, respectively, that are exerted at the axle of a wheel by the rover's body; z is the wheel sinkage; θ1 is the entrance angle at which the wheel begins to contact the soil; θ2 is the leaving angle at which the wheel loses contact with the soil; θm is the angle at which the maximum normal stress and 0 shearing stress occur; θ1 is the angle where the soil begins to deform in the shearing direction; r and h are the wheel radius and lug height, respectively; v is the velocity of the wheel axle, and ω is the angular velocity of the wheel. The soil interacts with the wheel in the form of a continuous radial normal stress σ and tangential shearing stress τ, which are divided into a forward part (σ1 and τ1) from angle θ1 to θm and a rear part (σ2 and τ2) from angle θm to θ2. A lugged wheel is considered to shear with the soil at an equivalent shearing radius rs: r s ¼ r þ λs h;

ð1Þ

where λs is a coefficient ranging from 0 to 1 and can be calculated by (Ding, 2009a):

λs ¼ 1 

3ðr þ hÞ : 8½1 þ cotðγ L =2ÞcotX c h

ð2Þ

In (2), γL denotes the included angle between two adjacent wheel lugs and γL ¼2π/nL, where nL is the number of lugs. Xc ¼ π/ 4  φ/2, representing the slide angle of the soil. If λs ¼ 1, the shearing is considered to act at the maximum radius R ¼r þh; if λs ¼ 0, the shearing is considered to act at the minimum radius. An approximate value of λs is 0.5, indicating that the shearing occurs in the middle of the lugs. Let s denote the slip ratio of a wheel to describe the relationship between its actual linear velocity and theoretical circumferential velocity caused by rotation. It is defined using shearing radius rs (Ding et al., 2009b): ( ðr s ω vÞ=r s ω ðr s ω Z v; 0 r s r1Þ s¼ : ð3Þ ðr s ω o v;  1 r s o 0Þ ðr s ω vÞ=v

Fig. 1. Free-body diagram for lugged wheels (Ding et al., 2014a).

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

The normal stress is calculated with (4).   8 N kc N > < σ 1 ðθÞ ¼ b þkφ r ð cos θ  cos θ1 Þ   n oN h i > : σ 2 ðθÞ ¼ kbc þkφ r N cos θ1  θθ θθ2 ðθ1  θm Þ  cos θ1 m 2

ðθ m r θ r θ 1 Þ ðθ 2 r θ o θ m Þ

;

ð4Þ where kc is the cohesive modulus of the soil, kφ is the frictional modulus. According to Bekker (1969), in the pressure–sinkage experiments, b is the radius of a circular plate or the smaller dimension of a rectangular plate. For a wheel, b should be the smaller dimension between the wheel width and the length of contact patch, however, as the value of kc/b is much smaller than kφ for the sandy terrain in this study, and it is difficult to determine the variable contact patch length before the wheel sinkage is known. Therefore, the wheel width b is used in this study. N is the improved soil sinkage exponent. N is the linear function of the slip ratio (Ding et al., 2010a): N ¼ n0 þ n1 s;

ð5Þ

where n0 and n1 are coefficients for calculating N. The linearized exponent is effective for predicting the entire wheel sinkage, including the static sinkage and severe dynamic slip–sinkage, which cannot be reflected well by the conventional model. The sinkage exponent N is an increasing function of the slip ratio, whereas the normal stress is a decreasing function of N, because the wheel sinkage z is generally less than 1 m. Therefore, the wheel sinkage should increase with an increase in the slip ratio to compensate for the decrease in the normal stress. Eq. (5) has been used by Zhou et al. (2014) to predict the slip–sinkage of Mars rovers while developing a simulation software. Ding et al. (2014b) carried out comprehensive research on using variant sinkage exponent of terrains to reflect the dynamic wheel sinkage. The angles θ1, θ2, and θm are functions of wheel sinkage z1 and coefficients c1, c2, and c3:

θ1 ¼ arccos½ðr  z1 Þ=r;

ð6Þ

θ 2  c3 θ 1 ;

ð7Þ

θm ¼ ðc1 þ c2 sÞθ1 :

ð8Þ

The shear stress is a function of the normal stress:

τðθÞ ¼ ½c þ σ ðθÞ tan φ  f1  exp½  jðθÞ=Kg ðh Z hmin Þ;

ð9Þ

where c is the cohesion of the soil, φ is the internal friction angle, and K is the shearing deformation modulus. The condition of hZ hmin should be satisfied to ensure that shearing occurs inside the soil, using the shearing stress model of (9). The wheel should typically be designed with a sufficient lug height to generate the maximum tractive force. The empirical equation for calculating the minimum lug height hmin is (Ding, 2009a) 0

″

hmin ¼ ð1  λÞhmin þ λhmin 0



hmin ¼ r cos ″

ð0 r λ r 1; λ  0:5Þ; 

2π 2π π þ sin tan X c þ nL nL nL

hmin ¼ r tan ðπ =nL Þ tan X c :



 1 ;

ð10aÞ ð10bÞ ð10cÞ

The coefficient λ ranges from 0 to 1 to reflect the friction condition between the lugs and the soil. If the lugs or the soils are lubricious, making no friction between the lugs and the soil, (10b) should be used, and λ ¼0. If the friction coefficient between the lugs and the soil is sufficiently large, (10c) could be used and λ ¼1. If the friction coefficient is not sufficiently large, λ can take an approximate value of 0.5. If the entrance angle is θ1 and the wheel sinkage is small, the shearing soil displacement can be considered as starting from 0 angle θ1 calculated by the maximum radius r þh. If the wheel

123

sinkage is large, the next lug cannot contact the soil immediately after the former one completely enters the soil. The interaction between the lugs and soil is quite complex. The radius Rj for cal0 culating θ1 is deduced for approximating the effect on the shearing displacement of the soil caused by the wheel lugs. Transitional slip ratios sj1 and sj2 are adopted because the wheel sinkage is related to the slip ratio. If the slip ratio is large enough, the influence of the wheel lugs on the starting angle of the soil deformation can be ignored. The empirical values of the parameters sj1 and sj2 are approximately 0.15 and 0.5, respectively. An improved equation for calculating the soil deformation is (Ding et al., 2014a, 2009c) 0

0

jðθÞ ¼ r s ½ðθ1  θÞ  ð1  sÞð sin θ1  sin θÞ; 

ð11Þ



θ01 ¼ acos ðr  zÞ=Rj ;

ð12Þ

8 r þh > < Rj ¼ r þ hðsj2 sÞ=ðsj2  sj1 Þ > :r

ð0 r s r sj1 Þ ðsj1 o s o sj2 Þ :

ð13Þ

ðsj2 r s r 1Þ

The concentrated forces and moment are calculated using (14). The improved models for calculating the normal and shearing stresses are used, and the equivalent shearing radius rs is used instead of r to consider the shearing process. Eqs. (5) and (11) contribute the most to the improvement of the model, because they reflect the wheels' slip–sinkage phenomena and the lug effect quite well. More details can be found in (Ding et al., 2014a). At quasi-static state with low velocity, the normal force FN, drawbar pull FDP, and resistance moment MDR can be predicted using the following integral model (M1): (Z θm FN ¼ b ½r σ 2 ðθÞ cos θ þ r s τ2 ðθÞ sin θdθ θ2

Z θ1

þ

θm

(Z F DP ¼ b þ

) ½r σ 1 ðθÞ cos θ þ r s τ1 ðθÞ sin θdθ θm

θ2

Z θ1 θm

ð14aÞ

¼ f DP ;

ð14bÞ

½r s τ2 ðθÞ cos θ  r σ 2 ðθÞ sin θdθ )

½r s τ1 ðθÞ cos θ  r σ 1 ðθÞ sin θdθ

"Z M DR ¼ r 2s b

¼ W;

θm θ2

τ2 ðθÞdθ þ

Z θ1 θm

#

τ1 ðθÞdθ ¼ T D :

ð14cÞ

3. Experimental study with single-wheel testbed A systematic experimental study on wheel–terrain interaction mechanics was conducted using a single-wheel testbed, and the soil parameters were tested by performing standard pressure– sinkage experiments and shearing experiments. 3.1. Single-wheel experiments The wheel–terrain interaction testbed (Fig. 2) developed by the State Key Laboratory of Robotics and System in the Harbin Institute of Technology, China, was used to study the sinkage and interaction mechanics of different wheels. In order to create various slip ratios, the testbed contains a driving motor to drive the wheel traversing forward and a carriage motor with a conveyance belt to imitate the influence of the rover on the wheel. The wheel sinkage was measured using a linear potentiometer displacement sensor; the interaction forces/torques were measured by a six-

124

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

Fig. 2. Wheel–terrain interaction testbed and experimental wheel (Ding et al., 2010b, 2013).

dimension F/T sensor over the wheel; and a torque sensor was installed on the driving axle to measure the driving torque directly. Six types of cylindrical metal wheels with different radii, widths, and lugs were used (Table 1). The wheels were comparable to most of that of current planetary rovers in terms of both their dimensions and mechanical properties. The minimum heights of the wheel lugs estimated by (10a)-(10c) for Whij (i¼ 1, 2; j¼ 1, 2) ″ 0 were hmin ¼ 14:81 mm and hmin ¼ 9:87 mm, and those for Wh3j ″ 0 (j¼ 1, 2) were hmin ¼ 14:73 mm and hmin ¼ 9:18 mm, respectively. The wheels with a lug height of 15 mm could satisfy the condition of h Zhmin , whereas the wheels with a lug height of 10 mm were marginal. If there was enough friction between the lugs and the soil, a 10-mm lug height was sufficient. However, if λ ¼0.5, the minimum lug height hmin for wheels Wh1j and Wh2j was 12.34 mm and that for Wh3j was 11.96 mm, which are both larger than 10 mm. The lunar soil simulant HIT-LSS1 used in this study was made from soft sand after removing the impurities, sieving, ventilating, and drying, and it was provided by the Shanghai Academy of Spaceflight Technology. The forward velocity of wheels has little influence on the interaction mechanics of planetary exploration rovers (Ding et al., 2013). Thus, it was set to a small value of 10 mm/s to obtain a sufficient amount of original data. The experimental slip ratios were 0,0.1, 0.2, 0.3, 0.4, and 0.6 by setting different velocities of the driving motor and carriage motor, and the slip ratios of 0.05 and 0.5 were set for some experimental conditions. The slip ratios were calculated using the radius r, and the values were amended with the shearing radius rs. The values of λs when calculating rs were 0.75 for wheels with 15-mm-high lugs and 0.65 for wheels with 10-mm-high lugs, which are close to the values estimated by (2). The vertical wheel load was approximately 80 N. All these experimental conditions were comparable to those of planetary rovers. Hundreds of raw data points were obtained from each test. The measured data fluctuated periodically in correspondence with the entrance and exit of the wheel lugs. The experimental results showed that the mean values of several tests were almost the same, regardless of the data fluctuations. The steady state data were used to calculate the mean values of FDP, FN, MDR and z, after filtering. The entrance angle θ1 was then calculated according to the wheel sinkage. The curves of FDP, FN, MDR and z(θ1) versus the slip ratio were obtained for the six types of experimental wheels. 3.2. Soil parameter tests using standard methods The standard methods developed for conventional terrestrial rovers (Bekker, 1969; Wong, 2010) were used to test the soil

Table 1 Parameters of experimental wheels. Test no.

Wheel code

r (mm)

b (mm)

h (mm)

nL

T1 T2 T3 T4 T5 T6

Wh11 Wh12 Wh21 Wh22 Wh31 Wh32

135 135 135 135 157.35 157.35

165 165 110 110 165 165

15 10 15 10 15 10

24 24 24 24 30 30

parameters, which were used to validate the parameter identification approaches. The mechanical properties of the terrain were divided into the bearing property in the normal direction and shearing property in the tangential direction. The bearing property was characterized by the pressure–sinkage relationship, whereas the latter was characterized by the shear stress–displacement relationship. The single-wheel testbed was used to test the soil parameters by replacing the driving wheels with experimental plates. The parameters kc, kφ, and n were estimated by measuring the sinkage in pressure–sinkage experiments using plates with different radii (35 mm and 50 mm) under vertical loads varying from 20 N to 300 N, as shown in Fig. 3(a). Shearing experiments were carried out under vertical loads varying from 30 N to 180 N to measure the maximum steering moment, based on which the shearing property parameters were estimated, including c and φ, as shown in Fig. 3(b). The shearing deformation modulus K was not estimated as the precision is low. The maximum radius of the shearing plate is 150 mm and the radius of the circle surrounded by the inner side of the shearing lugs is 115 mm. Table 2 lists the measured physical and mechanical parameters of the experimental soil, including its density ρs.

4. Sensitivity analysis of mechanics model In the wheel–terrain interaction mechanics model, the wheel parameters PW ¼{r, b, h} could be measured directly from a wheel, and the motion state parameters PM ¼{s, z (θ1)} could be estimated using the data acquired from the sensors. The remaining parameters PS ¼ {kc, kφ, n0, n1, c, φ, K, c1, c2, c3, λs, sj1, sj2} related to the terrain properties were divided into four groups, as listed in Table 3. Because the parameters λs, sj1, and sj2 could be determined empirically, they are called empirical parameters and are represented by PIV. Apart from λs, sj1, and sj2, there are 10 unknown model parameters that still need to be identified: parameters PI ¼{kc, kφ, n0, n1}, reflecting the bearing performance

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

125

Fig. 3. Experiments of measuring soil property parameters (Ding, 2009a), (a) Pressure–sinkage experiment (b) Shearing experiment. Table 2 Parameters of lunar soil simulant HIT-LSS1. kc (kPa/mn  1)

kφ (kPa/mn)

n

c (kPa)

φ (deg)

ρs (kg/m3)

Value

15.6

2407.4

1.10

0.25

31.9

1610

θm

θ1 Table 3 Interaction mechanics model parameters related to soil properties. Category PI

PII

Bearing performance parameters

Shearing performance parameters

PIII Contact angle coefficients

kc , k n0 , n1

σ (θ )

c, , k

τ (θ )

θ2 −( FDPσ 1 + FDPσ 2 )

FDPτ 1 + FDPτ 2

Symbol Remark kc kφ n0 , n1

Cohesive modulus of soil (kPa/mn  1) Frictional modulus of soil (kPa/mn) Coefficients for calculating the sinkage exponent N, N ¼ n0 þ n1s

c φ K

Cohesion of soil (kPa) Internal friction angle of soil (deg) Shearing deformation modulus of soil (m)

c1, c2

Coefficients for predicting maximum stress angle θm, θm ¼ (c1 þc2)θ1 Coefficient for predicting leaving angle θ2, θ2 ¼ c3θ1

c3

PIV Empirical parameters

FNτ 1 + FNτ 2

FNσ 1 + FNσ 2

M DRτ 1 + M DRτ 2

Parameter

λs

sj1, sj2

Coefficient for calculating equivalent shearing radius rs, λs A[0,1], calculated with Eq. (2). Transitional slip ratios, approximate values are: sj1 A [0.1, 0.2], sj2 A [0.4, 0.6].

of the soil; PII ¼{c, φ, K}, reflecting the shearing performance of the soil; and PIII ¼ {c1, c2, c3}, which are the contact angle coefficients. To understand the influence of different groups of parameters on the interaction mechanics, including the normal force, drawbar pull, and driving resistance moment, the coupling degree and parameter sensitivity should be analyzed. 4.1. Coupling degree analysis Eqs. (14a)–(14c) are three coupled functions, each of which relates to all the parameters in Table 3. The parameters influence the interaction mechanics through the normal and shearing stresses. Thus, the roles of the stresses on the interaction mechanics should be analyzed.

Fig. 4. Diagram of coupled wheel–terrain interaction mechanics model.

Rθ Rθ Let F Nσ1 ¼ rb θm1 σ 1 ðθÞ cos θ dθ; F Nσ2 ¼ rb θ2m σ 2 ðθÞ cos θ dθ, R θ1 R θm F Nτ1 ¼ r s b θm τ1 ðθÞ sin θ dθ; F Nτ2 ¼ r s b θ2 τ 2 ðθÞ sin θ dθ Rθ Rθ F DPσ1 ¼ rb θm1 σ 1 ðθÞ sin θ dθ; F DPσ2 ¼ rb θ2m σ 2 ðθÞ sin θ dθ, R θ1 Rθ F DPτ1 ¼ r s b θm τ1 ðθÞ cos θ dθ; F DPτ2 ¼ r s b θ2m τ2 ðθÞ cos θ dθ, M DRτ1 R R θ θ ¼ r 2s b θm1 τ1 ðθÞdθ; M DRτ2 ¼ r 2s b θ2m τ2 ðθÞdθ, then 8 > < F N ¼ F Nσ1 þ F Nσ2 þ F Nτ1 þ F Nτ2 F DP ¼ F DPτ1 þ F DPτ2  ðF DPσ1 þ F DPσ2 Þ : ð15Þ > :M ¼M DR DRτ1 þ M DRτ2 The coupled influence of the parameters and stresses on the integral mechanics is shown in Fig. 4. It is clear that the normal force, drawbar pull, and resistance moment are functions of all the unknown parameters. Eq. (15) shows highly coupled integral equations, and the number of parameters is much larger than the number of equations. As a result, it is infeasible to identify all the unknown parameters directly from the equations using the measured data. In Fig. 4, F Nσ ¼ F Nσ1 þ F Nσ2 , representing the component of the normal force generated by the radial normal stress; F Nτ ¼ F Nτ1 þ F Nτ2 , representing the component caused by the tangential shearing stress; F DPτ ¼ F DPτ1 þ F DPτ2 , representing the tractive force generated by the shearing stress; and F DPσ ¼ F DPσ1 þ F DPσ2 , representing the resistance force caused by the normal stress. The driving resistance moment M DR ¼ M DRτ1 þ M DRτ2 is only a function of the shearing stress. Wong (2008) provided the parameters for 21 different types of terrain, including dry sand, sandy loam, clay, snow, etc., where the internal friction angle φ is from 6° to 35°. The angle φ of the lunar soil is relatively large, ranging from 25° to 50° and usually smaller than 45°. The internal cohesion c is generally very small compared

126

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

Fig. 5. F Nτ =F N and F DPσ vs. F DPτ for wheel Wh31.

Table 4 Normal Force Component Caused by Shearing Stress (Using Parameters of Group G1, s ¼ 0.6). Wheel

Wh11

Wh12

Wh21

Wh22

Wh31

Wh32

Wh32

Wh32

FN (N) FNτ (N) FNτ/W(%)

80 13.2 16.6

80 12.3 15.4

80 139 17.3

80 13.1 16.4

80 12.0 15.0

80 11.4 14.4

35 4.14 11.8

150 24.7 16.4

with the normal stress and shearing stress for a frictional terrain with low cohesion such as dry sand. According to (9), the maximum shearing stress is c þ σtanφ, which is usually lower than the normal stress, and if the slip ratio is not large, the shearing stress will decrease because of the term ð1  e  j=K Þ. Figures for the normal stress and corresponding shearing stress distribution in Iagnemma et al. (2004) and Ray (2009) also verify that the shearing stress is obviously lower than the normal stress. The normal force component F Nτ is the integration of τsinθ, whereas F Nσ is the integration of σcosθ. The angle θ should not be too large so as to avoid getting a wheel stuck, particularly wheels moving on deformable terrain. As a result, the force component F Nτ is much smaller than F Nσ . Fig. 5 shows F Nτ =F N of Wh31 using the parameters of groups G1, G2, G3, G11, G14, and G20, which are described in Section IV.B. The ratio of F Nτ =F N increases with an increase in the slip ratio, whereas that for Wh31 is less than 24% if the slip ratio is not larger than 0.6, even if the parameters vary considerably. Table 4 indicates that F Nτ =F N is less than 18% for all the experimental wheels under different vertical loads, using the parameters of group G1. The force component F Nτ could be neglected to eliminate the coupling effect caused by the shearing stress and the related parameters c, φ, and K. The simplification error could be compensated by the identified parameters. Fig. 5(b) and (c) shows the components of F DPσ and F DPτ , respectively. The value of F DPσ is lower than that of F DPτ , and the difference is the effective drawbar pull. F DPσ increases linearly with an increase in the slip ratio because the wheel sinkage increases rapidly with a slip ratio increase, whereas the tractive force F DPτ increases exponentially, primarily as a result of the term ð1  e  j=K Þ when calculating the shearing stress. F DPτ is only influenced by the shearing performance parameters PII but F DPσ is primarily influenced by parameters PI and PIII. 4.2. Parameter sensitivity analysis According to Fig. 5, given the entrance angle θ1 (or wheel sinkage z) and the slip ratio s, the interaction mechanics could be

Fig. 6. Combined sinkage modulus versus b for lunar soil.

predicted. Actually, the vertical load W of a wheel is mainly determined by the rover's weight and the geometry of terrain surface that the rover is traversing on; the slip ratio s of a wheel is independent as it can be controlled by adjusting the angular velocity of the wheel. Therefore, W and s will be considered as given values, and the sensitivity of z, FDP, and MDR to the unknown parameters will be analyzed. Because of the highly coupled characteristics and the nonlinear nature of the equations, it is hard to analyze the parameter sensitivity using an analytical method. Thus, a numerical analysis method is primarily used. 4.2.1. Sensitivity of bearing performance parameters There are four parameters that influence the bearing performance of the soil in group PI, i.e., kc, kφ, n0, and n1. The parameters kc and kφ cannot be identified separately by a rover because the widths of its wheels are usually the same. If we define the combined sinkage modulus Ks ¼kc/bþ kθ, which is called lumped pressure–sinkage coefficient in Hutangkabodee et al. (2006), then the relationships between Ks, kc, and kφ should be discussed. Parameter kc reflects the influence of the wheel width on the sinkage modulus. For frictional soil with a small internal cohesion, such as planetary soil, the value of kc is also very low compared with kφ. If the smaller value between wheel width and contact patch length changes from 0.1 m to 0.5 m, the ratio of (kc/b)/Ks decreases from 6.08% to 1.28% for the HIT-LSS1 soil. The typical values of kc and kφ for lunar soil are 1.4 kPa and 820 kPa/m (n ¼1), respectively, and the variation in Ks versus b, the smaller value between wheel width and contact patch length, is shown in Fig. 6. If the wheel width is larger than 0.1 m, the ratio of (kc/b)/ Ks is less than 1.68%. According to the above analysis, identifying

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

parameter Ks using a wheeled rover is sufficient for a real-time application. The influence of parameter kc on the bearing performance of the soil is relatively small for frictional soil with little cohesion. If Ks is obtained, kφ can be considered to have approximately the same value as Ks, and kc can be simplified to zero. Ks will be considered as the sinkage modulus instead of kc and kφ. The experimental results for Wh31 were predicted using the measured soil parameters, and the unknown parameters were determined using a trial-and-error method; these parameters are listed in Table 5. The parameters of Ks, n0, and n1 can then be changed to investigate their influence on the interaction mechanics (Table 6 the other parameters are the same as those in Table 5). The parameters are used to predict the interaction mechanics (Fig. 7). Fig. 8 shows the stress distributions predicted with different parameters. The maximum relative errors compared with the experimental results are calculated using (16) as presented in Table 6. Err ¼ ðmax j F m  F p j Þ=maxF m  100%

ð16Þ

where Fm is the measured data of the interaction mechanics, and Fp is the predicted results. According to the results in Table 6, Figs. 7 and 8, the influence of the bearing performance parameters of PI on the wheel sinkage is much larger than their influence on the drawbar pull and driving resistance moment. The influence of the parameters in PI on MDR is ignorable, although the distribution of shearing stress τ is obviously influenced. Because the wheel sinkage z is much lower than 1 m, according to the property of exponential functions, the normal stress is a decreasing function of the sinkage exponent parameters n0 and n1. The vertical load is mainly sustained by the component of F Nσ as the integration of the normal stress, which is a decreasing function of the sinkage exponent parameters. As a result, the wheel sinkage z caused by a constant vertical load W is an increasing function of parameters n0 and n1. In contrast, an increase in Ks could decrease the wheel sinkage, indicating that the bearing performance of the soil is enhanced. An increase in wheel sinkage can cause a slight increase in the shearing stress τ, MDR, and F DPτ by increasing the shearing displacement of the soil. The component F DPσ also increases to a relatively larger extent than F DPτ . The drawbar pull, that equals F DPτ  F DPσ , i.e., the difference between soil traction and compaction resistance, decreases with an increase in the wheel sinkage, which could be influenced by the bearing performance parameters. Although the parameters of G1, G8, and G9 include quite different Ks, by modulating the parameters of n0 and n1, we could

Table 5 Parameters for wheel Wh31 (Group G1). Ks (kPa/mN)

n0

n1

c (kPa)

φ (deg)

K (mm)

c1

c2

c3

2500

0.76

1.30

0.25

31.9

25

0.5

 0.3

0

127

obtain almost the same predicted results for the interaction mechanics with high accuracy compared with the experimental results. This implies that the bearing performance of the soil could be characterized using only two independent parameters. The parameter Ks can be empirically estimated, and the accuracy of the interaction models could be ensured by modulating the parameters n0 and n1. 4.2.2. Sensitivity of shearing performance parameters Different groups of shearing performance parameters, as listed in Table 7, were used to analyze their influence on the wheel– terrain interaction mechanics. Fig. 9 shows the results of MDR and FDP predicted using these parameters, and the maximum errors are presented in Table 7. Based on these figures and errors, the following conclusions can be drawn. (1) The parameters of PII have no direct influence on the normal stress, and the influence of the shearing stress on force FN is small. Thus, they have little influence on the wheel sinkage z. (2) The driving resistance moment MDR and drawbar pull FDP are increasing functions of soil parameters c and φ, and decreasing functions of parameter K. This conclusion can also be drawn from analyzing the influence of these parameters on the shearing stress according to (9). (3) The influence of internal cohesion c on MDR and FDP is obvious when the slip ratio is large, whereas parameter K mainly influences MDR and FDP when the slip ratio is small, and it determines their increasing tendency. (4) Different groups of shearing performance parameters (groups 1 and 16) could obtain almost the same prediction results. Therefore, the three parameters c, φ, and K are not independent. 4.2.3. Sensitivity of contact angle coefficients Little research has been done on the influence of the contact angle coefficients, and many researchers assume that c1 ¼ 0.5, c2 ¼0, and c3 ¼0. The parameters of groups G17 to G22, as listed in Table 8, were used to analyze their influence on z, FDP, and MDR. From the errors shown in Table 8 and the curves in Fig. 10, the influence on MDR is very small and could be neglected. Except for group G20, which changed c2 greatly from 0.3 to  0.3, the influence on the wheel sinkage is also not obvious, especially when the results are compared with the influence seen in the bearing performance experiment. The influence of parameters c1, c2, and c3 on the drawbar pull FDP is relatively large compared with their influence on the others. However, compared with the influence of parameters PI and PII on the interaction mechanics, the influences of parameters c1, c2, and c3 are not very significant. Because no standard method yet exists for measuring the contact angle coefficients, choosing the typical values from the literature or setting c1 ¼0.5, c2 ¼0, and c3 ¼ 0 for simplification is acceptable if model accuracy is not strictly required, especially when the error could be compensated by the other identified parameters. However, if the soil parameters obtained from classical terramechanics experiments are used, the influence of the typical values for the contact angle coefficients on the interaction mechanics, especially the drawbar pull, should be noted.

Fig. 7. Influence of bearing performance parameters on z, FDP, and MDR.

128

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

Fig. 8. Influence of bearing performance parameters on stresses (s¼ 0.2). Table 6 Bearing performance parameters and their influence on prediction error of interaction mechanics. Group

Ks (kPa/mN)

n0

n1

z (%)

FDP (%)

MDR (%)

G1 G2 G3 G4 G5 G6 G7 G8 G9

2500 2500 2500 2500 2500 5000 800 5000 800

0.76 0.6 0.9 0.76 0.76 0.76 0.76 0.88 0.60

1.30 1.30 1.30 1.10 1.50 1.30 1.30 1.39 1.08

0.93 27.79 29.33 21.09 25.05 28.09 74.38 1.67 1.53

8.20 10.56 12.34 7.83 8.57 8.91 19.99 8.59 8.78

4.55 6.11 9.19 4.55 4.55 5.56 11.89 4.90 4.96

Table 7 Shearing performance parameters and their influence on prediction error of interaction mechanics. Group

c (kPa)

φ (deg)

K (m)

z (%)

FDP (%)

MDR (%)

G1 G10 G11 G12 G13 G14 G15 G16

0.25 0 2.5 0.25 0.25 0.25 0.25 1.00

31.9 31.9 31.9 42 25 31.9 31.9 25

0.012 0.012 0.012 0.012 0.012 0.020 0.005 0.010

0.93 1.33 4.22 2.43 2.53 1.85 1.20 0.95

8.20 15.02 97.11 60.70 40.76 36.85 49.96 10.40

4.55 8.02 53.53 33.03 21.94 20.49 29.54 4.40

Based on the theoretical equations, the influence tendency of parameters PIII could also be analyzed. If the sinkage z is given, the maximum stress angle becomes larger with increases in parameters c1 and c2, as shown in Fig. 11. Thus, the maximum normal stress decreases, and the generated normal force decreases. In contrast, if the vertical load is given, the wheel sinkage z is an increasing function of parameters c1 and c2. Because the integration effect of the shearing stress implied by MDR is not obvious, the drawbar pull will decrease with an increase in wheel sinkage because F Nσ is larger. The influence of parameter c3 on the calculation of the leaving angle is easy to imagine. With an increase in c3, the wheel sinkage decreases, and the drawbar pull increases. It is found that modulating two of the three interaction mechanics parameters is sufficient to ensure the fidelity of the model. Because the soil in contact with the rear part of a wheel cannot generate a large normal stress as a result of the flow of the soil, parameter c3 can be set to zero. 4.2.4. Concluding remarks The influences of the soil parameters on the wheel–terrain interaction are summarized in Table 9 based on the above analysis. The following conclusions could be drawn:

(1) The pressure–sinkage parameters in PI have very good relevancy with the wheel sinkage z. On one hand, PI has a very large influence on the wheel sinkage and a small influence on FDP, whereas the influence on MDR is ignorable. On the other hand, the wheel sinkage is mainly influenced by parameter PI, whereas the influence of PII is ignorable, and the influence of PIII is small. Two of the three parameters, Ks, n0, and n1, are independent. Thus, Ks can be empirically estimated to identify n0 and n1. (2) The resistance moment MDR has good relevancy with the shear stress–displacement parameters in PII because it is hardly influenced by the other parameters. Meanwhile, parameter PII also has considerable influence on the drawbar pull because of the term F DPτ . (3) The contact angle parameters do not have a significant influence on the interaction mechanics. Thus, because there is no standard method to measure these parameters, it is reasonable to use their typical values. However, they do influence the interaction mechanics, and their influence on the drawbar pull FDP is relatively large compared with their influence on z (FN) and MDR. These parameters should be identified if a high accuracy model is required. If the parameters in PII are known, FDP is mainly influenced by the parameters in PIII. The value of parameter c3 could be given, such as zero, and the other two parameters, c1and c2, can be identified. According to the above analysis, the parameter identification process is obtained: identifying Ks, n0, and n1 using the normal force and wheel sinkage data; identifying parameters c, φ, and K using the driving resistance moment; givinga typical value to c3 (zero is given in this study) to identify c1 and c2 using the drawbar pull force and the identified c, φ, and K.

5. Parameter identification based on integral-form model (IM1) 5.1. Parameter identification method: IM1 The parameters PI ¼{Ks, n0, n1}, PII ¼{c, φ, K}, and PIII ¼{c1, c2} could be identified based on the measured data for z, MDR, and FDP, respectively. The other parameters must be given or known when identifying one group of parameters. The accuracy of the given parameters could influence the parameter identification accuracy because of the coupling effect. A method to solve this deadlock problem is the most challenging issue. If the influence of the parameters in PIII is neglected by using F N  F Nσ , and setting typical values for P~ III (c1 ¼0.5, c2 ¼ 0, c3 ¼0), the parameters of group PI could be identified according to the measured FN and z. Then, P~ III and the identified P^ I could be used to

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

129

Fig. 9. Influence of shear performance parameters on FDP and MR.

Table 8 Contact angle coefficients and their influence on prediction error of interaction mechanics. Group

c1

c2

c3

z (%)

FDP (%)

MDR (%)

G1 G17 G18 G19 G20 G21 G22

0.5 0.4 0.6 0.5 0.5 0.5 0.5

 0.3  0.3  0.3 0 0.3  0.3  0.3

0 0 0 0 0 0.2  0.2

0.93 10.34 9.58 19.28 58.30 11.70 7.22

8.20 16.62 14.87 23.92 55.28 24.35 19.06

4.55 3.88 5.11 4.55 4.55 5.28 4.05

identify the shearing property parameters of PII. The parameters in PIII could in turn be identified using the identified P^ I and P^ II . However, this step-by-step identification process is not strict. Because a simplified FN is used, and the identified P^ III is not the same as the P~ III used for estimating the parameters of the other two groups, the identified P^ I and P^ II are not accurate, which further influences the accuracy of P^ III . It is found that an iterative parameter identification process can solve this problem, i.e., the identified parameters in P ^ III are used to obtain the parameters of

^P 0 ; P^ 0 , and P^ 0 . If

P^ 0  P^ III

o δ (δ is the error tolerance), the paraI II III III 0

0

0

meters P^ I ; P^ II , and P^ III are the results; otherwise, the iteration process is continued. The Matlab function lsqcurvefit(), which was developed to solve nonlinear data-fitting problems in a least-squares sense, is applied to realize the parameter identification. It can find parameter x to be identified by solving the problem min‖Fðx; sÞ  y‖22 ¼ min x

x

m X

½Fðx; si Þ  yi 2 ;

the measured soil parameters kc and kφ. The identified results of Ks are quite close to the initial value for different wheels, with the parameters n0 and n1 ranging 0.776–0.913 and 0.824–1.266, respectively. This further verifies that the degree of freedom of the bearing performance parameters is 2, and Ks is redundant. c1 and c2 range from 0.299 to 0.465 and  0.258 to 0.017, respectively. These parameters vary for different wheels. A Whi1 wheel with 15-mm-high lugs has a smaller n0, larger n1, smaller c, larger φ, and larger K than a Whi1 wheel with 10mm-high lugs. The shearing performance parameters for different wheels are close; c and φ can be considered to be soil parameter constants for different wheels, and K can be considered to be a wheel–terrain interaction parameter, which changes for different wheels. In the bearing performance parameters, Ks can be considered to be a soil parameter, whereas n0 and n1 can be treated as wheel–terrain interaction parameters that vary with different wheels. Parameters c1 and c2 are easily influenced by the dimensions of different wheels because the integral forces and moment are not sensitive to them. Thus, they are also treated as wheel– terrain interaction parameters. The maximum data fitting error for FDP is 7.56%, whereas those for FN and MDR are less than 5%. Fig. 12 shows the convergence process for the identified parameters. With the convergence of parameters c1 and c2, the other parameters also converge to get reasonable results. Several iterations should be performed before satisfying the convergence conditions, which takes up tens of seconds for the calculation (with a laptop PC@2 GHz). This approach could be used for high-accuracy soil parameter identification, but the long calculation period restricts its application to real-time situations.

ð17Þ

i¼1

given input data s (slip ratio) and the observed output y (measured FN, FDP, or MDR), where s and y are vectors of length m (the number of experiments with different slip ratios), and F(x, s) is a function of the interaction mechanics, i.e., FN(  ), FDP(  ), or MDR(  ). The other parameters included in F(  ), such as the entrance angle θ1 and the other known parameters, could be input as global variables. The Matlab programming is shown in Appendix (Table 17).

6. Closed-form decoupled analytical model of wheel–terrain interaction 6.1. Linearization of stress distribution models Let σm denote the maximum normal stress, maximum shearing stress, then

τm denote the

σ m ¼ K s rN ð cos θm  cos θ1 ÞN ; 5.2. Experimental verification If we let sj1 ¼ 0.1, sj2 ¼0.5, and δ ¼0.002, the vertical load is equal to the average measured normal force, which is approximately 80 N. The identified parameters are listed in Table 10. The internal cohesion of the soil, c, ranges 0.097–0.268 kPa; the internal friction angle φ ranges 29.2–31.1°. Both these parameters are close to the measured values of c and φ, 0.251 kPa and 31.9°, respectively. The shearing deformation modules are 10.1–13.8 mm, with small variations for different wheels. The results are comparable to the soil parameters found in the literature such as Wong (2008). The initial value of Ks is set to 2500 kPa/mN according to

ð18Þ

τm ¼ Eðc þ σ m tan φÞ:

ð19Þ 0 1

0 1

θm Þ  ð1  sÞð sin θ sin θm Þ=Kg. The where E ¼ 1  expf  r s ½ðθ normal stress and shearing stress can be simplified using the linearized method (Shibly et al., 2005): ( σ 1 ðθÞ  σ L1 ðθÞ ¼ σ m ðθ1  θÞ=ðθ1  θm Þ ðθm r θ r θ1 Þ ; ð20Þ σ 2 ðθÞ  σ L2 ðθÞ ¼ σ m ðθ  θ2 Þ=ðθm  θ2 Þ ðθ2 r θ o θm Þ (

τ1 ðθÞ  τL1 ðθÞ ¼ τm ðθ1  θÞ=ðθ1  θm Þ ðθm r θ r θ1 Þ : τ2 ðθÞ  τL2 ðθÞ ¼ τm ðθ  θ2 Þ=ðθm  θ2 Þ ðθ2 r θ o θm Þ

ð21Þ

130

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

Fig. 10. Influence of contact angle coefficients on z and FDP. Table 9 Influence of soil parameters on interaction mechanics. Parameter PI{Ks, n0, n1} Parameter PII{c, φ, K} Parameter PIII {c1, c2, c3} z (FN) Very large FDP Small MDR Very Small

Very small Large Large

Small Relatively large Very small

τ1 ðθÞ cos θ  τ01 ðθÞ ¼ τ0m ðθ1  θÞ=ðθ1  θm Þ τ2 ðθÞ cos θ  τ02 ðθÞ ¼ τ0m ðθ  θ2 Þ=ðθm  θ2 Þ :

ð23Þ

6.2. Derivation of closed-form analytical model 6.2.1. Closed-form model-CM1 By substituting (20) and (21) for (14) and integrating the equations, one arrives at (Ding et al., 2009c) 3 2 3 2 FN A B   6 F 7 6 B A 7 X : ð24Þ 5 4 DP 5 ¼ 4 Y M DR 0 rs C cos θm  cos θ2  cos θ 1 þ cos θθm  ; X ¼ rb m ; B ¼ θm  θ 2 θm 1 sin θm  sin θ1 θ1  θ 2 , Y ¼ r s b m and C ¼ 2 . A, B, and C are θ 1  θm

σ

Z θ1 θm

θm

θ2

#

τL1 ðθÞ sin θ dθ ¼ rbðθ1  θ2 Þσ m cos θm =2

þr s bτm B  rbðθ1  θ2 Þσ m cos θm =2 ¼ DX;

The simplification errors of (22) and (23) are close to those of (20) and (21), respectively, and data analysis shows that their differences are less than 1% for most cases.

where A ¼

θ2

þ

In Shibly et al. (2005), the use of a linearized method for soils with a sinkage exponent between 0.5 and 1.6 is validated. Because the entrance angle θ1 is not very large, and the normal stress and shearing stress are zero or very small, the product of cos θ and the stress could also be linearized. If we let σ ‘m ¼ σ m cos θm and τ‘m ¼ τm cos θm , then ( σ 1 ðθÞ cos θ  σ 01 ðθÞ ¼ σ 0m ðθ1  θÞ=ðθ1  θm Þ ð22Þ σ 2 ðθÞ cos θ  σ 02 ðθÞ ¼ σ 0m ðθ  θ2 Þ=ðθm  θ2 Þ ; (

is simpler than model CM1, is obtained: "Z # "Z Z θ1 θm θm F N  rb σ 02 ðθÞdθ þ σ 01 ðθÞdθ þ r s b τL2 ðθÞ sin θdθ

sin θm  sin θ 2 θm  θ 2

þ τ functions of entrance angle θ1 and the parameters of PIII; X is related to the parameters of PI and PIII, whereas Y is related to all the parameters. It is clear that FN, FDP, and MDR are functions of all the soil parameters. AX and BY are components of the normal force generated by the normal stress and shearing stress, respectively; whereas AY and BX are components of the drawbar pull generated by the shearing stress and normal stress, respectively. 6.2.2. Closed-form model-CM2 By using the stress simplification equations, (20)–(23), in the integral wheel–terrain interaction model and neglecting the normal force caused by the shearing stress, closed-form model II-CM2, which

F DP  r s b rb

Z θm θ2 Z θ1 θm

τ02 ðθÞdθ þ rs b

Z θ1 θm

τ01 ðθÞdθ  rb

Z θm θ2

ð25aÞ

σ L2 ðθÞ sin θ dθ

σ L1 ðθÞ sin θdθ ¼ rs bðθ1  θ2 Þτm cos θm =2

rbσ m B ¼ DY  BX;

ð25bÞ

M DR ¼ r s CY:

ð25cÞ

where D ¼ ðθ1  θ2 Þ cos θm =2. D has the dimension of radian; rD has the dimension of length (m); and rbD has the dimension of area (m2). DX ¼ rbDσ m can be considered to be the volume of a triangular prism with a base length of rðθ1  θ2 Þ, width of b, and height of σ m cos θm ; and it can also be considered to be the normal stress σ m cos θm =2 acting on a flat surface with an area of rðθ1  θ2 Þ U b. Accordingly, parameters A, B, and C also have the dimension of radian, and the meaning of AX, BY, AY, BX, and CY can also be analyzed using a method similar to that used for analyzing DX. The compact form of (25) is 2 3 2 3 FN D 0   6 F 7 6 B D 7 X : ð26Þ 4 DP 5 ¼ 4 5 Y M DR 0 rs C 6.2.3. Decoupling of closed-form analytical models Eqs. (24) and (26) are closed-form analytical equations as functions of X and Y. However, most of the equations still contain all three groups of parameters: PI, PII, and PIII. To realize parameter identification, decoupled equations as functions of a limited number of parameters should be deduced. 6.2.3.1. Decoupled closed-form model based on CM1–M2. According to (24), one obtains X ¼ ðF N  BYÞ=A ¼ ðAY  F DP Þ=B;

ð27Þ

Y ¼ M DR =r s C:

ð28Þ

By substituting (27) and (28) into (24), the mechanics models become (Ding et al., 2009c) F DP ¼

! B2 BF N ¼ Aþ Y A A

! A2 þ B2 M DR B  FN A AC rs

¼ F DP ðc1 ; c2 ; c3 ; θ1 ; s; M DR ; F N Þ;

ð29Þ

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

131

Fig. 11. Influence on stress distribution of contact angle coefficients (s¼ 0.2).

F N ¼ AX þ BY ¼ AX þ BM DR =ðr s CÞ ¼ rbAσ m þBM DR =ðr s CÞ ¼ F N ðK s ; n0 ; n1 ; P III ; θ1 ; s; M DR Þ:

ð30Þ

According to (30),

σ m ¼ F N =ðrbAÞ  BMDR =ðrrs bACÞ:

ð31Þ

By substituting (19), (31), and Y ¼ r s bτm , into (24), it is possible to obtain  r 2 CE bc þ F N tan φ=ðrAÞ ¼ M DR ðc; φ; K; P III ; θ1 ; s; F N Þ: ð32Þ M DR ¼ s 1 þ r s BE tan φ=ðrAÞ

Table 10 Parameter identification results of method IM1. Wheel

Wh11

Wh12

Wh21

Wh22

Wh31

Wh32

c1  c2 Ks* n0 n1 c (kPa) φ (deg) K (mm)

0.388 0.258 2499.0 0.823 1.257 0.103 30.9 13.1

0.299 0.083 2499.3 0.913 1.013 0.218 29.6 11.3

0.443 0.234 2499.3 0.873 1.059 0.212 31.1 13.8

0.465 0.204 2499.6 0.933 0.824 0.236 29.2 10.6

0.400 0.253 2499.0 0.776 1.266 0.097 30.5 11.0

0.339 0.017 2499.3 0.849 1.063 0.268 29.3 10.1

*

6.2.3.2. Decoupled closed-form model based on CM2-M3 . F N ¼ rbðθ1  θ2 Þσ m cos θm =2 ¼

F 0N ðK s ; n0 ; n1 ; P III ;

θ1 ; sÞ:

ð33Þ

According to (33) and (25c), the maximum normal stress

σmand shear stress τm are approximately: σ m  2F N =½rbðθ1  θ2 Þ cos θm :

ð34Þ

τm ¼ Y=ðrs bÞ ¼ 2MDR =½r 2s ðθ1  θ2 Þ:

ð35Þ

By substituting (34) and (35) into (25b) and (25c), the drawbar pull and resistance moment are obtained: F DP ¼ M DR cos θm =r s  2BF N =½ðθ1  θ2 Þ cos θm  ¼ F 0DP ðc1 ; c2 ; c3 ; θ1 ; s; M DR ; F N Þ:

ð36Þ

M DR ¼ Er 2s ½bcðθ1  θ2 Þ=2 þ F N tan φ=ðr cos θm Þ ¼ M 0DR ðc; φ; K; P III ; θ1 ; s; F N Þ:

ð37Þ Eqs. (29), (30), and (32) are decoupled equations for simplified model M2, whereas (33), (36), and (37) are the decoupled equations of simplified model M3. FDP is a function of the parameters in PIII; MDR is a function of the parameters in PII and PIII; and FN is a function of the parameters in PI and PIII. Model II is simpler than Model I because of the simplification method for the stress distribution and the omission of the term F Nτ . Model M3 is easier to understand and costs less computation time. However, its effectiveness should be evaluated.

7. Parameter identification based on closed-form decoupled analytical model The closed-form decoupled analytical models are used to realize real-time parameter identification for the planetary rovers. The

The unit of Ks is kPa/mN.

identification method that uses simplified model M2 is called IM2, and the method that uses simplified model M3 is called IM3.

7.1. Parameter identification methods: IM2 and IM3 The drawbar pull is only a function of the unknown parameters of PIII ¼ {c1, c2}, and the influence of the other parameters are included in the measured FDP and FN. Therefore, these parameters could be identified using the drawbar pull and the other observed data directly. If c1 and c2 are known, the parameters PI ¼{Ks, n0, n1} and PII ¼ {c, φ, k} could be identified using the measured normal force and drawbar pull, respectively. 7.1.1. Identification of c1 and c2 According to (29) and (36), FDP is a function of FN, MDR, θ1, s, and the unknown parameters are c1 and c2, which can be identified if the other parameters have been measured, as shown in Fig. 13. FN,FDP, MDR, θ1, and s for a wheeled rover can be measured or estimated (Iagnemma et al., 2004). 7.1.2. Identifying Ks, n0, and n1 By substituting (5) and (18) into (30) and (33), respectively, it is possible to obtain  B M DR ; F N ¼ rbA K s r n0 þ n1 s ð cos θm  cos θ1 Þn0 þ n1 s þ C rs

ð38Þ

 F N ¼ rbðθ1  θ2 Þ K s r n0 þ n1 s ð cos θm  cos θ1 Þn0 þ n1 s cos θm =2:

ð39Þ

According to (38), three parameters, Ks, n0, and n1, can be identified based on the measured θ1, FN, MDR, and s, and the identified c1 and c2. If (39) is used, it is unnecessary to measure the moment of the resistance. The parameter identification process is shown in Fig. 14.

132

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

Fig. 12. Convergence process for identified parameters.

Fig. 15. Identifying c, φ, and K using methods IM2 and IM3. Fig. 13. Diagram for identifying c1 and c2 with methods IM2 and IM3. Table 11 Parameter identification results using constant sinkage exponent. Wheel code N

Ks (Kpa/m ) n

Fig. 14. Identifying Ks, n0, and n1 using methods IM2 and IM3.

7.1.3. Identifying c, φ, and K By substituting (19) into (32) and (37), we can obtain 0

M DR ¼

0

r 2s C½bc þ F N tan φ=ðrAÞ  ð1  expf r s ½ðθ1  θm Þ  ð1  sÞð sin θ1  sin θ m Þ=KgÞ ; 0 0 1 þ ½r s B tan φ=ðrAÞ  ð1 expf  r s ½ðθ1  θm Þ  ð1 sÞð sin θ1  sin θm Þ=KgÞ

ð40Þ M DR ¼ r 2s ½bcðθ1  θ2 Þ=2 þ F N tan φ=ðr cos θm Þ

¼ f1  expð r s ½ðθ‘1  θm Þ  ð1  sÞð sin θ‘1  sin θm Þ=KÞg: ð41Þ

Based on (40) and (41), the parameters c, φ, and K can be identified by using the measured θ1, FN, MDR, and s, as shown in Fig. 15. 7.2. Parameter identification results with method IM2 The least square method was used to identify the model parameters with the experimental data using the lsqcurvefit() function of Matlab. The relative data fitting error values of FDP, FN, and MDR for all the wheels are less than 5%. The computation time for fitting FDP

Wh11

Wh12

Wh21

Wh22

Wh31

Wh32

1.32  0.541

1.62  0.506

2.26  0.513

2.13  0.539

1.55  0.502

1.35  0.535

and FN is less than 50 ms using a 2-GHz laptop PC, but it takes hundreds of milliseconds to fit MDR because of the complexity of (40). The largest error for the identified φ is only 3.1°. The error in c is also very small compared to the shear stress. Parameter K has an acceptable range (9.7–13.1 mm). Parameter Ks is sensitive to the initial value, and the identified result is quite close to the initial value of 2500 kPa/mN, indicating that it is sufficient to use only n0 and n1 for fitting the vertical load. For dry sand, Ks could be fixed to the value of estimated soil frictional modulus kφ, because kc is small enough to be negligible (Wong, 2010). For example, for lunar soil, Ks can be set to kφ ¼820 kPa/mN by ignoring kc ¼1.4 Kpa/ mN þ 1. Ks can also be set to the value of stiffness modulus which is an intrinsic terrain parameter with determined unit of Pa/m that has dominant role in governing the bearing performance (Ding et al., 2014b). The bearing performance of soil could be quite well characterized using n0 and n1. n0 þn1s for the experimental lunar soil stimulant is between 0.73 (Wh31, s¼ 0) and 2.10 (Wh11, s¼1), from which the wheel sinkage into the soil can be estimated. c1 is between 0.37 and 0.53, and c2 is between  0.38 and  0.04. Their values exhibit a wide range, but the other identified parameters are not obviously influenced by them. This also implies that the wheel–terrain interaction mechanics are not sensitive to c1 and c2. If we let c1 ¼0.5 and c2 ¼0, the drawbar pull can be fitted

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

133

Table 15 Data fitting error and calculation time of IM3. Wheel

Wh11 Wh12 Wh21 Wh22 Wh31 Wh32

Fig. 16. Data fitting result of FN for Wh11.

Table 12 Parameter identification results For Rj ¼ r. Wheel

Wh11

Wh12

Wh21

Wh22

Wh31

Wh32

c (kPa) φ (deg) K (mm) Relative error of MR (%)

3.071 15.7 2.1 7.72

0.448 26.2 5.5 7.74

2.972 19.7 2.5 7.65

0.634 25.9 4.8 8.8

3.215 15.8 1.3 10.58

1.847 20.0 2.9 8.64

Table 13 Parameter identification results for Rj ¼ R. Wheel

Wh11

Wh12

Wh21

Wh22

Wh31

Wh32

c (kPa) φ (deg) K (mm) Relative error of MR (%)

0.281 29.9 12.3 15.21

0.214 28.2 11.3 13.19

0.263 30.0 12.1 14.94

1.035 26.0 9.6 13.95

0.731 27.3 9.8 14.68

0.227 28.0 9.5 15.77

Table 14 Identified soil parameters of IM3. Wheel

Wh11

Wh12

Wh21

Wh22

Wh31

Wh32

c1  c2 Ks * n0 n1 c (kPa) φ (deg) k(mm)

0.458 0.416 2499.0 0.772 1.293 0.183 27.7 11.2

0.353 0.265 2499.2 0.866 1.054 0.344 25.5 10.1

0.471 0.361 2499.2 0.828 1.097 0.601 25.9 10.6

0.386 0.159 2499.6 0.914 0.824 0.740 24.2 9.0

0.460 0.395 2498.7 0.727 1.302 0.500 26.1 9.2

0.365 0.204 2499.2 0.805 1.124 0.193 26.0 8.6

*

The unit of Ks is kPa/mN.

with a maximum relative error in FDP of 8.57% for Wh11. The data fitting errors for FN and MDR for all the wheels are not significantly influenced, but the maximum relative errors in FDP for all the wheels are approximately 10%, whereas Wh12 has the largest error of 18.7%. Moreover, the results for the other parameters are only slightly influenced. If high accuracy is desired in predicting FDP or identifying the soil parameters, c1and c2 should be identified, rather than given typical values. If a constant sinkage exponent is used instead of N in (5), the identification results for Wh11 are Ks ¼14.2 kPa/mn and n¼ 0 (lowest limit). Extending the limits of n, the best data fitting results for Ks are from 1.32 kPa/mn (Wh11) to 2.66 kPa/mn (Wh21), and n shows results from  0.54 (Wh11) to  0.50 (Wh31), as

Relative error (%)

Calculation time (ms)

FDP

FN

MDR

FDP

FN

MDR

2.67 3.14 3.59 2.04 3.86 3.67

1.98 1.37 0.40 0.73 4.08 0.85

3.34 1.20 0.67 1.22 3.85 2.10

14.788 14.137 22.247 13.065 17.436 21.916

20.415 21.664 21.730 18.332 23.583 17.842

24.019 51.576 22.867 64.362 22.255 29.209

shown in Table 11. The fitting errors for FN of using constant sinkage exponent n are quite large, as Fig. 16 shows. If Eqs. (11)–(13) are not used, i.e., the influence of the lugs on the starting angle of the shearing deformation of the soil is ignored, the identified results using Rj ¼r are shown in Table 12: c ranges from 448 Pa (Wh12) to 3071 Pa (Wh11); φ, from 15.7° (Wh11) to 26.2° (Wh12); and k, from 1.3 mm (Wh31) to 5.5 mm (Wh12). These results, however, are untenable, and the fitting error of MDR reaches 10.58% (Wh31). If Rj ¼R is used, the data fitting error for MDR is too large, with values of approximately 15%, as shown in Table 13. The effectiveness of (11)–(13), which reflects the wheel lug effects, is thus verified. The process of predicting the interaction mechanics based on the given vertical load (W¼80 N) can be found in Ding et al. (2014a). The maximum relative error values for predicting z, FDP, and MDR using the identified parameters are 1.89%, 9.09%, and 4.28%, respectively. The calculation time for z is less than 30 ms, and that for FDP and MDR is about 30 μs. These are feasible values for real-time application. 7.3. Parameter identification results with method IM3 Tables 14 and 15 present the results for the parameter identification and calculation time, respectively, using method IM3. The identified internal friction angle φ is smaller than the measured results, whereas the internal cohesion is larger. This is caused by the effect of the simplification error, especially the omission of the coupled term, on the normal force by the shearing stress. The identified parameters have errors to compensate for the model simplification error. However, the maximum relative error of the data fitting is less than 5%. If the identified parameters are used to predict the interaction mechanics with model M3, the accuracy is acceptable. The advantages of this approach are that it is the simplest and the calculation time is the shortest. 7.4. Comparison of different parameter identification methods Figs. 17 and 18 compare the identified results, data fitting errors, and calculation times. All the identification methods have high goodness-of-fit values, implying that all three models can be used with the corresponding identified parameters. Both IM1 and IM2 could be used to estimate the soil parameters c and φ with high accuracy. IM1 has high accuracy in identifying the soil parameters, but it takes a long time. The accuracy of the parameters identified with IM3 is the lowest, but it has the simplest form and costs only tens of ms for the calculation. Model M3 with the corresponding parameters can be used in a situation that requires very high speed calculations. The decoupled analytical model M2 and the identification method IM2 that was developed based on it are a compromise, because the identified parameters c and φ also have high accuracy, and the calculation speed is fast.

134

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

Fig. 17. Comparison of identified parameters.

Fig. 18. Comparison of data fitting errors and calculating time.

8. Experimental validation with a four-wheeled rover If the contact mechanics model is applied to a multi- wheeled robot moving over very rough terrain, the loads distributed on the different wheels differ. The load-effect of the interaction models should be considered. The decoupled closed-form model M2 is improved to (42) to reflect the load-effect: F N ¼ rbAσ m þ BM DR =ðr s CÞ;

ð42aÞ

F DP ¼ ½M DR ðA2 þB2 Þ=ðr s ACÞ  BF N =A ð1 þ cP1 þ cP2 sÞ½1 þ cDP ðW  F N Þ=W; M DR ¼

r 2s CE½bc þ ½1 þ cM ðW  F N Þ=WF N tan φ=ðrAÞ : 1 þr s BE tan φ=ðrAÞ

Fig. 19. Photos of El-Dorado II moving on flat terrain.

ð42bÞ ð42cÞ

where W is the average vertical load of the robot on the wheels; cDP and cM are used to reflect the influences of the vertical load on the normal force and resistance moment, respectively; and the term ð1 þ cP1 þ cP2 sÞ is adopted to reflect the influence of the slip ratio as the contact angle coefficients are set to constant values: c1 ¼ 0.5, c2 ¼0, c3 ¼0. To apply the model to a four-wheeled El-dorado II robot (Yoshida, 2009) that was developed by Tohoku University in Japan, experiments were performed to obtain the basic input data. Fig. 19 shows the robot moving on flat terrain made from Toyoura soft sand, with a normal force generated by the weight to generate different slip ratios. The resistance forces were set to 0, 10, 20, 30, 35, 40, 50, 60, 70 N. In addition, the robot was controlled to climb up slopes with angles of 0°, 3°, 6°, 9°, 12°, and 15°.

The unknown parameters were identified using the in-situ wheel–soil contact data, and the results are presented in Table 16. The soil parameters were kept almost the same, and the coefficients were changed to reflect the difference between climbing up slopes and moving on flat terrain, as well as the experimental errors. The identified results of soil parameter Ks was close to the initial value provided by Ishigami (2008), and the identified results of c and φ were close to the values provided in the same reference. The identified parameters were then used for a dynamics simulation of the robot's motion. Comparisons of the simulation and experimental results are shown in Figs. 20 and 21. In addition to the motion of the robot, which was indicated by the slip ratio, the drawbar pull, moment of resistance, normal force, and wheel sinkage could also be predicted with high fidelity. In operational missions, different angular velocities can be set to the wheels to modulate their slip ratios. The actual slip ratios of wheels can be measured by visual odometry and encoders. The wheel sinkage can be measured by cameras when the rover stops

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

135

Table 16 Parameter identification results for interaction model of Eldorado II robot and Toyoura sand. Parameters

cP1

cP2

cDP

cM

Ks (kPa/mN)

n0

n1

c (kPa)

φ (deg)

K (mm)

Moving on flat terrain Climbing slope Reference (Ishigami, 2008)

 0.276  0.379 –

0.633 0.616 –

 0.304  0.448 –

0.354 0.214 –

1796 1796 1790

0.63 0.66 1.00

0.72 0.72 0

0.025 0.025 0.0

35.8 35.8 38.0

10.5 10.5 –

Fig. 20. Simulation and experimental results while moving on flat terrain.

Fig. 21. Simulation and experimental results while climbing up slopes.

moving. The driving torque can be estimated from the motor currents, while the drawbar pull and normal force can be calculated through quasi-static analysis. Alternatively, the contact responses can be measured directly by installing a F/T sensor at the wheel.

9. Conclusions This study provided parameter identification approaches for the improved terramechanics model with high fidelity; a model that could be used to estimate soil parameters with high accuracy was also provided. Concluding remarks on the work reported in this paper are provided as follows: (1) Because the influence of the contact angle parameters is not significant, it is reasonable to use their typical values. Their minor influence could be only considered if a high accuracy model is required, and they can be identified by the drawbar pull. The leaving angle coefficient c3 can be set to zero for simplification. (2) The pressure–sinkage parameters can be identified according to the wheel sinkage z because they have very high relevancy with z. Ks can be empirically estimated according to the soil property, and n0 and n1 can be identified as wheel–terrain interaction parameters. (3) The shear stress–displacement parameters can be identified using the driving resistance moment MDR because of their close relevancy. The parameters c and φ identified with the

integral model are very close to the soil parameters, and K can be considered to be a wheel–terrain interaction parameter. (4) Either the integral model or the decoupled analytical models could be selected according to different mission requirements because their accuracy can be ensured by the identified parameters based on in-situ experimental data. The iterative identification method IM1 has higher accuracy in identifying the soil parameters. Model M3 has the simplest form and IM3 requires the least computation. The decoupled analytical model M2 and the associated identification method IM2 represent a compromise. The accuracy of parameter identification methods and models has been verified using experimental data obtained from both a single-wheel testbed and a four-wheeled rover.

Acknowledgment This study was supported in part by the National Natural Science Foundation of China (Grant no. 61370033/51275106), National Basic Research Program of China (Grant no. 2013CB035502), Foundation of Chinese State Key Laboratory of Robotics and Systems (Grant no. SKLRS201401A01/SKLRS-2014-MS-09), the Fundamental Research Funds for the Central Universities (Grant no. HIT.BRETIII.201411), Harbin Talent Programme for Distinguished Young Scholars (No. 2014RFYXJ001), Postdoctoral Youth Talent Foundation of

136

L. Ding et al. / Planetary and Space Science 119 (2015) 121–136

Heilongjiang Province, China (Grant no. LBH-TZ0403), and the “111 Project” (Grant no. B07018).

Appendix

Table 17 Parameter identification programming with integral-form model (Method IM1). 1. 2.

Give wheel parameters; empirical parameters: λs (0.75 for wheels Whi1, and 0.65 for wheels Whi2), λj1 ¼ 0.1 and λj2 ¼ 0.5; δ ¼0.002. Load measured parameters: ðiÞ ðiÞ m fsðiÞ ; zðiÞ ; F ðiÞ N ; F DP ; M DR gi ¼ 1 , where m is the number of experimental slip ratio set up for a certain wheel.

3.

4.

5. 6. 7. 8. 9.

ðiÞ m m Calculate entrance angle fθðiÞ 1 gi ¼ 1 (Eq. (6)), radius fRj gi ¼ 1 (Eq. (14)), angle

fθ10ðiÞ gm i ¼ 1 (Eq. (12)), and r s ¼ r þ λs h (Eq. (1)). Give: x0_sink¼ [2.5  106, 0.5, 0.5], lb_sink ¼ [0.1  106,0,0], ub_sink ¼ [10  106, 2, 2]; x0_shear¼ [200, 30, 0.01], lb_shear¼ [0,0,0.002], ub_shear ¼ [5000,50,0.050]; x0_angcoef ¼[0.5, 0], lb_angcoef ¼[  1, 1], ub_angcoef ¼[1,1]. Give c1(k) ¼0.5, c2(k) ¼ 0, c3 ¼ 0. Set k ¼ 1. Do if (k ¼ ¼ 1) ½ K^s ðkÞ n^0 ðkÞ n^1 ðkÞ  ¼lsqcurvefit(@Fun_Fn_sigma ðiÞ m (s; θ1 ; c1 ðkÞ; c2 ðkÞ; c3 ; K s ; n0 ; n1 ), x0_sink,fsðiÞ gm i ¼ 1 ,fF N gi ¼ 1 , lb_ sink,

10. 11.

a m ub_sink,fθðiÞ 1 gi ¼ 1 ) else ½ K^s ðkÞ n^0 ðkÞ n^1 ðkÞ  ¼lsqcurvefit(@Fun_Fn_Int

(s; θ1 ; c1 ðkÞ; c2 ðkÞ; c3 ; K s ; n0 ; n1 ; cðk  1Þ; φðk  1Þ; Kðk  1Þ), 12. 13.

ðiÞ m ðiÞ m b x0_sink,fsðiÞ gm i ¼ 1 ,fF N gi ¼ 1 , lb_sink, ub_sink, fθ 1 gi ¼ 1 ). end ^ ^ ½ c^ ðkÞ φðkÞ KðkÞ  ¼ lsqcurvefit(@Func_Mdr_Intðs; θ1 ; c1 ðkÞ; c2 ðkÞ; c3 ; K s ðkÞ;

n0 ðkÞ; n1 ðkÞ; c; φ; KÞ, c m fθðiÞ 1 gi ¼ 1 ). 14.

½ c^ 1 ðk þ 1Þ

ðiÞ m x0_shear,fsðiÞ gm i ¼ 1 ,fF N gi ¼ 1 , lb_ shear, ub_ shear,

c^ 2 ðk þ 1Þ  ¼lsqcurvefit

(@Func_Fdp_Intðs; θ1 ; K s ðkÞ; n0 ðkÞ; n1 ðkÞ; cðkÞ; ðiÞ m φðkÞ; KðkÞ; c1 ; c2 ; c3 Þ, x0_shear,fsðiÞ gm i ¼ 1 ,fF N gi ¼ 1 , lb_shear, ub_shear, d m fθðiÞ 1 gi ¼ 1 ). 15. if (c^ 1 ðk þ 1Þ  c^ 1 ðkÞ o δ&&c^ 2 ðk þ 1Þ c^ 2 ðkÞ o δ) 16. break; 17. else 18. k ¼k þ 1; 19. End 20. End 21. Print parameters in PI:K^ s ¼ K s ðiÞ; n^ 0 ¼ n^ 0 ðkÞ; n^ 1 ¼ n^ 1 ðkÞ; ^ 22. Print parameters in PII: c^ ¼ c^ ðkÞ; φ^ ¼ φðkÞ; ^ K^ ¼ KðkÞ;

23. Print parameters in PIII:c^ 1 ¼ c^ 1 ðk þ 1Þ; c^ 1 ¼ c^ 2 ðk þ 1Þ. 24. End

ðiÞ m a In lsqcurvefit(@Fun_Fn_sigma, x0_sink,fsðiÞ gm i ¼ 1 ,fF N gi ¼ 1 , lb_ sink, ub_sink), F(  ) ¼ Func_Fn_sigma(  ), which is a function F Nσ ¼ F Nσ1 þ F Nσ2 according to (21); x0_sinkis the initial point of the pressure–sinkage parameters to identify; lb_ sink is the vector of the lower bounds; ub_sink is the vector of the upper bounds; fsðiÞ gm i¼1 ðiÞ m m is the input xdata, and fF ðiÞ N gi ¼ 1 is the observed ydata. fθ1 gi ¼ 1 and the given values of c1(k), c2(k), and c3 are input as global variables. ðiÞ m b In lsqcurvefit(@ Fun_Fn_Int, x0_sink, fsðiÞ gm i ¼ 1 , fF N gi ¼ 1 , lb_ sink, ub_sink, m fθðiÞ g ), the difference with *1 is that Fun_Fn_Int(  ) is the integral function of i ¼ 1 1 (15), and there are three additional global parameters: c(k  1), φ(k  1), and K (k  1). c F(  ) ¼ Func_Mdr_Int(  ) is the integral function of (15). θ1, c1(k), c2(k), c3, Ks(k), n0(k), and n1(k) are the global input parameters. d F(  ) ¼ Func_Fdp_Int(  ) is the integral function of (15). θ1, c3, Ks(k), n0(k), n1(k), c(k),φ(k), and K(k) are the global input parameters.

References Arvidson, R.E., Anderson, R.C., Haldemann, A.F.C., et al., 2003. Physical properties and localization investigations associated with the 2003 mars exploration rovers. J. Geophys. 108 (E12), 8070. Arvidson, R.E., Anderson, R.C., Bartlett, P., et al., 2004a. Localization and physical properties experiments conducted by Spirit at Gusev Crater. Science 305, 821–824. Arvidson, R.E., Anderson, R.C., Bartlett, P., et al., 2004b. Localization and physical property experiments conducted by Opportunity at Meridiani Planum. Science 306, 1730–1733. Bekker, M.G., 1969. Introduction to Terrain-rover Systems. The University of Michigan Press, Ann Arbor, MI, USA. Carrier, W.D., Olhoeft, G.R., Mendell, W., 1991. Physical properties of the lunar surface. In: Heiken, G.H., Vaniman, D.T., French, B.M. (Eds.), Lunar Source Book. The Cambridge University Press, New York. Cross, M., Ellery, A., Qadi, A., 2013. Estimating terrain parameters for a rigid wheeled rover using neural networks. J. Terramech. 50, 165–174. Ding, L., 2009a. Wheel–terrain Interaction Terramechanics For Lunar/planetary Exploration Rovers: Modeling and Application. Doctoral thesis of Harbin Institute of Technology, Harbin, China (in Chinese). Ding, L., Gao, H.B., Deng, Z.Q., et al., 2009b. Slip ratio for lugged wheel of lunar rover in deformable soil: definition and estimation. In: Proceedings of IEEE/RSJ International Conference of Intelligent Robots and Systems. St. Louis, MO, USA, October 11–15, pp. 3343–3348. Ding, L., Yoshida, K., Nagatani, K., et al., 2009c. Parameter identification for planetary soil based on a decoupled analytical wheel–terrain interaction terramechanics model. In: Proceedings of IEEE/RSJ International Conference on Intelligent Robots and Systems. St. Louis, MO, USA, October 11–15, pp. 4122–4127. Ding, L., Gao, H.B., Deng, Z.Q., Tao, J.G., 2010a. Wheel slip–sinkage and its prediction model of lunar rover. J. Central South Univ. Technol. 17, 129–135. Ding, L., Gao, H.B., Deng, Z.Q., Nagatani, K., Yoshida, K., 2010b. Experimental study and analysis on driving wheels' performance for planetary exploration rovers moving in deformable soil. J. Terramech. 48, 27–45. Ding, L., Deng, Z.Q., Gao, H.B., et al., 2013. Experimental study and analysis of the wheels' steering mechanics for planetary exploration WMRs moving on deformable terrain. Int. J. Robot. Res. 32, 712–743. Ding, L., Deng, Z.Q., Gao, H.B., et al., 2014a. Interaction mechanics model of driving wheels for planetary WMRs moving on deformable soil considering multiple effects. J. Field Robot. . http://dx.doi.org/10.1002/rob.21533 Ding, L., Gao, H.B., Deng, Z.Q., et al., 2014b. New perspective on characterizing pressure–sinkage relationship of terrains for estimating interaction mechanics. J. Terramech. 52, 57–76. ESA, 2013. Robotic exploration of Mars: ExoMars. 〈http://www.esa.int/SPECIALS/ ExoMars/SEM10VLPQ5F_0.html〉. Hutangkabodee, S., Zweiri, Y.H., Seneviratne, L.D., Althoefer, K., 2006. Performance prediction of a wheeled rover on unknown terrain using identified soil parameters. In: Proceedings of 2006 IEEE International Conference on Robotics and Automation. Orlando, FL, USA, pp. 3356–3361. Iagnemma, K., Kang, S., Shibly, H., Dubowsky, S., 2004. Online terrain parameter estimation for wheeled mobile robots with application to planetary rovers. IEEE Trans. Robot. 20, 921–927. Ishigami, G., 2008. Terramechanics-based Analysis and Control for Lunar/planetary Exploration Robots. Doctoral Thesis of Tohoku University Sendai, Japan. Japan Aerospace Exploration Agency, 2013. Moon lander SELENE-2. 〈http://www. jspec.jaxa.jp/e/activity/selene2.html〉. Lutz, D., 2013. New Mars rover's Mechanics to be used to Study Martian Soil Properties. 〈http://news.wustl.edu/news/pages/23139.aspx〉. NASA JPL, 2013. Mars Exploration Rover. 〈http://marsrovers.jpl.nasa.gov/home/ index.html〉. NASA, 2013. Mars Science Laboratory, Curiosity: Could Mars Have Once Harbored Life. 〈http://www.nasa.gov/mission_ pages/msl/index.html〉. Perko, H.A., Nelson, J.D., Green, J.R., 2004. Green mars soil mechanics investigation. NASA Technical Rep. No. NRA-00-01-MDAP-059, Colorado State University Fort Collins CO, USA. Ray, L.E., 2009. Estimation of terrain forces and parameters for rigid-wheeled rovers. IEEE Trans. Robot. 25, 717–726. T.P. Setterfield A. Ellery. Terrain response estimation using an instrumented rockerbogie mobility system IEEE Trans. Robot. 29 2013 172 188. Shibly, H., Iagnemma, K., Dubowsky, S., 2005. An equivalent soil mechanics formulation for rigid wheels in deformable terrain, with application to planetary exploration rovers. J. Terramech. 42, 1–13. The Rover Team, 1997. Characterization of the Martian surface deposits by the Mars pathfinder rover, sojourner. Science 278, 1765–1768. Wong, J.Y., 2008. Theory of Ground Rovers, fourth edition. John Wiley and Sons, Inc. Wong, J.Y., 2010. Terramechanics and Off-road Rover Engineering, 2nd edition. Elsevier. Yoshida, K., 2009. Achievements in space robotics. IEEE Robot. Autom. Mag. 16, 20–28. Zhou, F., Arvidson, R.E., Bennett, K., et al., 2014. Simulations of Mars rover traverses. J. Field Robot. 31 (1), 141–160.