- Email: [email protected]
Trafficability analysis of lunar mare terrain by means of the discrete element method for wheeled rover locomotion
Available online at www.sciencedirect.com
Journal of Terramechanics Journal of Terramechanics 47 (2010) 161–172 www.elsevier.com/locate/jterra
Trafficability analysis of lunar mare terrain by means of the discrete element method for wheeled rover locomotion Wen Li, Yong Huang *, Yi Cui, Sujun Dong, Jun Wang School of Aeronautical Science and Engineering, Beijing University of Aeronautics and Astronautics, Beijing 100191, China Received 15 February 2009; received in revised form 14 August 2009; accepted 11 September 2009
Abstract An irregularly shaped particulate system for simulation of lunar regolith is developed using discrete element modeling based on the fractal characteristics, particle shape, and size distribution of returned Apollo-14 samples. The model parameters are determined by dimensional analysis and biaxial test simulation with an improved boundary condition. Under terrestrial conditions, the trafficability of lunar mare terrain is estimated in terms of wheel–terrain interaction by experiment and simulation in order to validate the applicability of the wheel–terrain model employed here. The results show that the discrete element method combined with the wheel–terrain model is sufficiently accurate for mare terrain trafficability analysis without consideration of lunar environmental effects. To predict the trafficability of in situ lunar mare terrain, the non-contact forces attributed to the lunar surface environment are discussed and the initial mechanical model of discrete elements is modified by introduction of lunar gravitational force as well as electrostatic force. In the modified model, wheel–terrain interaction is analyzed under the same travel conditions as that of the experiment. The result shows the trafficability of the in situ lunar mare terrain is worse than that obtained by experiment and simulation with the initial model according to the value of horizontal force at any slip ratio. However, the wheel requires less drive torque on the moon than that on the earth. An explanation for these phenomena may be that lunar subsurface regolith particles are arranged in a looser manner under local environmental effects that effectively decrease the bearing and shearing strength of regolith. Ó 2009 ISTVS. Published by Elsevier Ltd. All rights reserved.
1. Introduction China is planning to send an unmanned lunar rover vehicle to the surface of the moon to make observations and measure the physical and mechanical properties of in situ regolith material during the second phase of the Chang’E Lunar Exploration Project. In consideration of mass and power consumption constraints, the lunar rover vehicle was generally designed with wheeled locomotion in previous and current lunar exploration missions. From prior experience reported by the Apollo and Luna missions, the lunar mare is a favorite landing site because of its relatively flat topographical condition and absence *
Corresponding author. Address: Mailbox 505, School of Aeronautical Science and Engineering, Beijing University of Aeronautics and Astronautics, 37 Xueyuan Road, Haidian District, Beijing 100191, China. Tel.: +86 10 8233 8498; fax: +86 10 8231 6654. E-mail address: [email protected] (Y. Huang).
of rocks. The Apollo-14 mission report [1] demonstrated that the lunar mare surface was generally covered with a loose mantle of material generated by micro-meteoritic impacts, and the population of rocks in the Fra Mauro area was much less than 0.5% of the total area. As a result, the soft-terrain trafficability of the lunar mare has become the foremost factor limiting the traction performance of a lunar rover during lunar surface exploration missions. Based on terramechanics, soft-terrain trafficability is defined as the capability of a soil-like material to support the traverse of vehicular traffic from one place to another [2]. It is related to the mechanical properties of the terrain on the macroscopic scale and the dynamic behavior of soillike particles on a granular scale. Current scientists have developed various methods to meet the requirements for rover design and mobility evaluation as well as lunar soft-terrain trafficability analysis. To predict the lunar rover’s mobility performance, Carrie [3] listed the recom-
0022-4898/$36.00 Ó 2009 ISTVS. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.jterra.2009.09.002
162
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
mended lunar soil trafficability parameters used for Bekker’s empirical equation. Tao et al. [4] and Patel et al. [5] analyzed the interaction between a rigid wheel and soft-terrain based on Terzaghi’s passive soil pressure theory and Bekker’s formula, respectively. Experimental validation was conducted using a specially designed and manufactured apparatus. Melzer [6] experimentally measured the Boeing LRV wheel’s performance by a dynamometer system and lunar soil simulant. Lizuka [7] studied the performance of a lugged wheel on sandy terrain using a model experiment in the laboratory. Kuroda et al. [8] evaluated the mobility performance of a planetary rover with a similar model experiment under low-gravity flight conditions. Khot et al. [9] simulated the dynamic wheel–soil interaction via the distinct element method and verified the simulation results with laboratory experiments. Nakashima et al. [10] numerically analyzed the influence of wheel configuration on traction performance for a lunar microrover by means of the discrete element method. In comparison of the existing methods for trafficability estimation, we find that: (i) The continuum mechanics method, based on the soil constitutive law and failure pattern, is unable to cover the flow of soil particles; therefore, it is difficult to analyze the trajectories of soil particles and explain soil deformation and failure mechanisms on the granular scale. (ii) The laboratory experiment measurement method can directly acquire trafficability parameters with a highaccuracy transducer system; however, lunar environmental simulation is costly and difficult to realize. (iii) The discrete or distinct element method (DEM), which was originally proposed by Cundall and Strack [11] to model the behavior of granular materials, seems to be a promising method for numerical analysis of the trafficability of lunar mare terrain, without the limitations of the methods mentioned above. In such an approach, the forces acting on an individual particle can be calculated using a well-founded contact reaction model. The trajectories and velocity vectors of particles can be determined based on Newton’s second law. With the development of a contact detection algorithm and discrete particle modeling method, DEM has become an ideal tool for terramechanics-based research. So far, however, the soil particles are commonly simplified as circular discrete element models in terramechanicsbased studies to reduce computational time during DEM simulation. Although Asaf [12,13] has developed a nonspherical DEM model in his research, such a model is far from realistic. Herein, we aim to develop a strategy for producing an irregularly shaped particulate system based on scanning electron microscopy (SEM) images and grain size distributions of lunar regolith and simulant, as described in Section 2. In addition, the mechanical model (contact reaction
model) and its corresponding microparameters are determined based on the physical and mechanical properties of lunar regolith. In Section 3, a series of DEM simulations with respect to validated experiments are performed to evaluate the applicability of the DEM model specified here. In Section 4, we discuss the influence of the lunar environment on dynamic behavior of regolith particles, and modify the initial DEM model by introducing non-contact forces attributed to lunar environmental effects. With the modified model, we estimate the trafficability of lunar mare terrain and explain the microscopic mechanism of terrain deformation and failure through trajectories and networks of contact force of regolith particles beneath the wheel. 2. DEM model for lunar regolith particles Lunar mare terrain trafficability is governed by the mechanical properties of in situ regolith and is mainly affected by geometric features and dynamic behavior of subsurface soil particles. Herein, we construct the geometric model and mechanical model of DEM for lunar mare regolith with high fidelity. 2.1. Geometric model for lunar regolith particles Based on knowledge of soil mechanics, the microscopically geometric parameters that influence the macromechanical properties of soil are particle shape, grain size distribution, and initial voidage. 2.1.1. Particle shape model Considering the fact that lunar regolith particles are generally irregular with re-entrant surfaces [14] and the irregularity of particle shape influences the interlock interactions among particles, we present an imaged-based modeling method to simulate the truly geometric features of lunar regolith particles. The required geometric feature parameters are acquired from SEM photography of returned lunar samples and corresponding simulants by using ImageJ software [15] and fractal calculation. The results show that the fractal dimension with perimeter–area relation for lunar regolith is 1.098, and particles with different sizes have a self-similar fractal behavior. Consequently, we build our DEM model with four typical particle shapes shown in Fig. 1. As is shown in Fig. 1, any irregular-shaped particle can be modeled by overlapping elements with a center element O and several angular elements Bi (i = 1, 2, 3, 4, . . .). The number of angular elements is determined by the complexity of the geometric features for a realistic particle. All elements are parametrically modeled in the Cartesian coordinate system, and both positions and radii are variables. To reduce the contact detection time, every particle shape model is clumped and behaves as a rigid body wherein the mutual contact check of constituent elements is neglected during DEM simulation.
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
Fig. 1. Particle shape model of DEM for lunar regolith.
2.1.2. Grain size distribution model The grain size distribution (GSD) by wt.% of returned lunar samples from the Apollo missions was experimentally measured and published in the literature [16]. Fig. 2 depicts the upper and lower bounds of the GSD for returned Apollo-14 samples as well as the GSD of the lunar mare simulant used in our experimental validation by solid and dashed curves, respectively. Here, we develop an effective method to generate the particulate system of the DEM with the best-fit GSD of the lunar mare simulant shown in Fig. 2. The strategy of producing irregularly shaped particles with a specified GSD is described as follows: (1) Generate ball elements with radii R and position (x0, y0) according to the GSD specified here, where x0, y0, R are variables in the Cartesian coordinate system. (2) Set contact stiffness of ball elements with high value (1.0 108 Pa) and without friction and gravity.
Fig. 2. Grain size distribution curves.
163
164
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
(3) Balance the particulate system under inter-element contact reactions in a virtual container constructed by rigid wall elements with higher stiffness than that of ball elements. (4) Reduce the radii of ball elements to R0 while keeping initial locations constant, then add angular elements based on the particle shape model given in Fig. 1. (5) Set best-fit values to DEM model parameters, such as density, spring stiffness, damping coefficient, and Coulomb friction coefficient, and balance the particulate system. Finally, we obtain the most appropriate model for our DEM analysis. 2.2. Mechanical model of lunar regolith particles The dynamic behavior of any particle is governed by the forces and torques acting upon it. As described by Apollo mission reports, there are no liquid water and other adhesive materials in lunar regolith; the cohesion and internal friction angle of regolith are mainly attributed to the interlock interactions among irregularly shaped particles. In our DEM model, we initially simplify the mechanical model of ball elements as a linear spring–dashpot system with Coulomb friction and non-tensile links, as shown in Fig. 3. As shown in Fig. 3, the inter-element contact force F can be decomposed into the normal force Fn and tangential force Fs. The normal force is determined by the elastic part F ne due to the normal spring, and the inelastic part F nd due to the normal damper and non-tensile link, i.e. As shown in Fig. 3, the inter-element contact force F can be decomposed into the normal force Fn and tangential force Fs. The normal force is determined by the elastic part F ne due to the normal spring, and the inelastic part F nd due to the normal damper and non-tensile link, i.e., kF n k ¼ F ne þ F nd
ð1Þ
F ne F nd
ð2Þ
¼ k n urn ¼ k n ðRj þ Ri dÞ ¼ cn vrn ¼ cn ðu_ j u_ i Þ n
ð3Þ
where kn is normal spring stiffness, urn is the penetration of two elements at the interface, d is the distance between element centers, Ri and Rj are element radii, cn is the normal damping coefficient, vrn is normal relative velocity, u_ i and u_ j
are the translational velocities of elements, and n is the unit normal vector. The tangential force includes tangential spring force, tangential damping force, and Coulomb friction force. It can be expressed by 8 old 1 Dt > < F s ðks 1 2csDtÞvrs Dt ðkF new s k 6 ljF n jÞ new ðks þ2cs Þ ð4Þ Fs ¼ > new s : ljF n j F new ðkF k > ljF jÞ new n s kF k s
where kF new s k is the norm of the current tangential contact force, F old is the tangential force in the previous timestep, s ks is the tangential spring stiffness, cs is the tangential damping coefficient, vrs is the relative tangential velocity at the contact point, l is the Coulomb friction coefficient, and Dt is the timestep able to avoid disturbance by propagating from the particle to a further range. 2.3. Determination of DEM model parameters The DEM model parameters, such as spring stiffness, damping coefficient, and Coulomb friction coefficient, are difficult to determine by direct physical measurement. As a result, various indirect methods have been developed. Cundall and Strack obtained the model parameters by using Hertzian contact theory and Mindlin and Deresiewich’s approach [17]. Liao and Chang [18] developed a method for determining the bounds of the spring constant between two particles by assuming an elastic bulk material. Oida [19] and Asaf [12,13] used the trial-and-error method to determine model parameters in their DEM simulation. In previous research, we confirmed that the method of dimensional analysis combined with biaxial test simulation can obtain the best-fit parameters of the DEM model [20]; therefore, we determine the model parameters required in Eqs. (1)–(4) by an improved biaxial test model and similarity criteria. Based on the Buckingham G theorem and dimensional matrix analysis, we obtain four dimensionless parameters Gi (i = 1, . . . , 4) and the relationship between microscopic parameters of the DEM model and macro-mechanical property parameters. F ðP1 ; P2 ; P3 ; P4 Þ ¼ 0 ck s , E2
Fig. 3. Contact reaction model between two elements.
ð5Þ ck n , s2max
kn , ks
Here, P1 ¼ P2 ¼ P3 ¼ P4 = l, where c is the bulk density of soil, E is the elastic modulus, and smax is the shear strength of the soil. According to the triaxial compression test of the lunar simulant mentioned above, we construct an improved DEM model for biaxial test simulation shown in Fig. 4. In this model, two flexible bounds are built with overlapped ball elements to accurately simulate the deformation of a soil sample during axial loading. Assuming that the radius of constituent elements of flexible bounds is Rb in addition to the width W0 and height H0 of the soil sample, then the number of elements for each flexible bound can be calculated by
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
165
Table 1 Parameters for boundary condition of DEM. Soil sample dimensions, W0 H0 (mm) Confining pressure, r3 (kPa) Density of boundary elements, qb (kg/m3) Radius of boundary elements, Rb (mm) Normal stiffness of spring, kn (kPa) Tangential stiffness of spring, ks (kPa) Bond between elements, sb (N/m) Friction coefficient of element, lb Element number of each bound, n Stiffness of walls, k (kPa) Friction coefficient of walls, lw Velocity of bottom wall, vwy (mm/min)
39.1 80.0 26.0 100.0 1.0 1000.0 1000.0 1.0 1010 0.2 80 1.0 106 0.0 0.8
and the best-fit parameters of DEM model are obtained and listed in Table 2 based on the dimensional analysis result. For comparison, the axial strain vs. axial stress difference curves obtained from both experiment and simulation are plotted in Fig. 5. It can be observed from Fig. 5 that the DEM model and its best-fit parameters are able to accurately reflect the soil constitutive law for lunar regolith and the simulant, and the DEM biaxial simulation result agrees with the results of the laboratory experiment.
Fig. 4. DEM model for biaxial test simulation.
3. Trafficability analysis for lunar mare terrain n¼
H0 Rd
ð6Þ
To control the boundary condition and acquire the radial deformation of soil sample during DEM simulation, the constituent elements of both bounds are indexed with an ID number: for the left bound, the ID number varies from 1 to n, and for the right bound, the ID number varies from n + 1 to 2n. Furthermore, every element of the bounds is preloaded with an x-component force Fbx. Given that the triaxial test confining pressure is r3, the force imposed on every bound element can be written by F bx ¼
2Rb H 0 r3 n
ð7Þ
where the positive or negative sign on the force is consistent with the Cartesian coordinate system shown in Fig. 4. For the left bound’s constituent elements (ID = 1, . . . , n), Fbx > 0; for the right bound’s constituent elements (ID = n + 1, . . . , 2n), Fbx < 0. During the biaxial test simulation, two rigid wall elements are generated on the top and bottom sides of the soil sample. The top wall is fixed and functions as a load cell with which axial force Fwy can be acquired by summing up the contact reaction forces of the soil particle elements acting on the wall. The bottom wall is given by a constant translational velocity vwy on the y-axis in order to simulate the loading process. The parameters of the boundary condition for our biaxial test simulation are listed in Table 1. Using the lunar regolith DEM model (geometric model and mechanical model) and the boundary conditions of the biaxial test, a series of biaxial test simulations are performed
In terms of soft-terrain, trafficability is usually analyzed by the wheel–soil interaction under various slip ratios. In this section, we first perform a series of experiments under terrestrial conditions to validate the applicability of the wheel–terrain DEM model. Then, the DEM models for both wheel configuration and lunar mare terrain are developed. Third, a series of DEM simulations for wheel–terrain interaction are carried out with respect to experiments and to validate the DEM models. Finally, the influence of lunar environmental effects on trafficability is discussed. The non-contact forces, which are critical to the microdynamic behavior of lunar regolith particles, are introduced into the initial mechanical model of DEM for soil particles. With this modified model, lunar mare terrain trafficability is predicted with high fidelity. 3.1. Experimental validation 3.1.1. Experimental apparatus and conditions Experimental measurements of lunar mare terrain trafficability are performed with a specially designed and manTable 2 Best-fitted model parameters for lunar soil particles. Density of element, q (kg/m3) Radius of center element, R0 (mm) Friction coefficient, l Normal stiffness, kn (kPa) Tangential stiffness, ks (kPa) Damping coefficient ratio, c/ccrit Increment of timestep, Dt (s)
2770 0.25–0.625 0.18 5.0 104 5.0 104 0.7 8.5 106
166
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
sured during the experiment. Because the carriage is mounted on the framework of the soil bin by two horizontally fixed guiding axles, it can traverse smoothly along the horizontal axis shown in Fig. 6. In our experimental measurement of lunar mare terrain trafficability, a rigid wheel with thin lugs is adopted and the lunar mare terrain is reproduced by the lunar simulant in the soil bin. The mechanical property of simulated terrain in the soil bin is measured with the cutting ring sample method and simple shear test. The wheel configuration and travel conditions are summarized in Table 3.
Fig. 5. Axial strain vs. axial stress difference curves.
ufactured wheel–terrain testbed for the study of planetary rover locomotion performance. The testbed primarily consists of a soil bin and lunar simulants, a carriage assembly, wheel assembly, and the data acquisition/control/processing system shown in Fig. 6. The soil bin is constructed from a metallic frame and transparent sidewalls (tempered glass) to aid visibility of soil particle flow beneath the tested wheel. The soil bin is designed to be large enough to eliminate boundary effects during the experiment. The wheel assembly is mainly composed of the tested wheel, a Maxon DC motor with controller, and six-axial force/torque sensors, and it is able to move freely along the vertical axis of the carriage. During experimental measurements, the tested wheel is driven and controlled by the motor module while the forces and torques exerted by the ground on the tested wheel are acquired by the six-axial force/torque sensor. The carriage assembly is equipped with a potentiometer, photoelectric encoder, and the driven motor module, with which the sinkage and translational velocity of the wheel can be mea-
3.1.2. DEM model for wheel–soil interaction In the DEM model for wheel–soil interaction analysis, the tested wheel configuration is modeled by clumped ball elements. The travel conditions of the wheel with respect to the experiment are input to the DEM model. To quantitatively estimate the soft-terrain trafficability by DEM simulation, we propose a method used to build the relationship between the microscopic results of DEM and the macromechanical parameters of trafficability based on coordinate transformation and algebraic summation. Fig. 7 shows a contact pair between wheel element i and soil element j as they contact each other, and also gives the schematic of forces exerted on element i. As shown in Fig. 7, element i is acted on by normal force F N i and tangential force F Si ; which can be calculated using Eqs. (1)–(4). By means of coordinate transformation, we obtain the contact reaction forces of element i in the Cartesian coordinate system. ( S F xi ¼ F N i sin hi F i cos hi ð8Þ S F yi ¼ F N i cos hi þ F i sin hi where hi = Arc tan [(x xj)/(y yj)]. Supposing that the number of constituent elements of the wheel model is Nw, then the horizontal force, wheel load, and drive torque of wheel can be expressed as FH ¼
Nw X
F xi
i¼1
T ¼
Nw X
F Si
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxi x0 Þ2 þ ðy i y 0 Þ2
ð9Þ ð10Þ
i¼1
Table 3 Parameters for experiment.
Fig. 6. Experimental apparatus.
Soil bin dimensions, L W H (mm) Cohesion of test terrain, c (kPa) Internal friction angle of test terrain, u (°) Initial void ratio of lunar simulant, e0 Wheel diameter without lug, D (mm) Lug dimensions, LT LH (mm) Wheel wide, B (mm) Number of lug on wheel, NL Vertical load on wheel, W (kg) Slip ratio, i (%) Angular velocity of wheel, x (rpm) Gravitational acceleration, g (m/s2)
3.0 0.8 0.6 0.28 39.95 0.85 220.0 3.0 25.0 160.0 18 11.8 0–100 5.0 9.8
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
167
Fig. 7. Contact reaction force on individual wheel element.
In this study, the lunar mare terrain is reproduced by an irregularly shaped particulate system (described in Section 2) in a virtual container. The container is 1.0-m long and 0.2-m deep, and constructed with rigid walls. To improve computational efficiency, we filled the soil particles with various scales. Namely, the soil particle size is enlarged with the increase in container depth. The equivalent diameter of a soil particle in our DEM model varies from 0.5 mm to 2.0 mm. The process of DEM simulation is divided into four phases: (i) lunar mare terrain reproduction with respect to the experiment; (ii) construction of a wheel model above the terrain surface, without any contact between wheel and terrain; (iii) wheel sinkage under a vertical load imposed on its center; (iv) wheel–terrain interaction simulation with specified travel conditions. In the wheel–terrain DEM model, there are 166, 135 discrete elements used, and the increment of the timestep is set to 8.85 106 s for the sake of stability and equilibrium of simulation. 3.1.3. Experiment and simulation results The trafficability of lunar mare terrain is analyzed by experiments and simulations in terms of wheel–soil interaction under specified conditions listed in Table 3. Figs. 8 and
9 show terrain deformation with a lugged wheel for 0.0 and 0.8 slip ratios by experiment and simulation, respectively. It is obvious that the terrain surface forms clear ruts along the travel direction as the wheel moves without slip; however, at higher slip (i = 0.8), the upper layer particles of the terrain shift backward and heap up behind the wheel due to the lug digging and shearing. Soil dilatancy or shear deformation occurs at that moment as well, as shown in Fig. 8b. Corresponding DEM simulations are performed with the wheel–terrain DEM model. The trajectories of soil particles and terrain deformation at 0.0 and 0.8 slip ratios are analyzed qualitatively and demonstrated in Fig. 9. Compared with experimental results, DEM simulation shows the same tendency in terrain deformation, that is to say, the terrain surface forms clear and wavy ruts at zero slip while soil dilatancy occurs behind the wheel at the 0.8 slip ratio. Consequently, we draw the conclusion that the DEM model for wheel–soil interaction can qualitatively estimate the terrain deformation and trajectories of subsurface soil particles with sufficient accuracy. To quantitatively validate the applicability of our DEM model and its parameters, the horizontal force and drive torque acting on the wheel center are calculated by means
Fig. 8. Experimental result for terrain deformation (a: i = 0; b: i = 0.8).
168
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
Fig. 9. DEM simulation for terrain deformation (i = 0.0 and 0.8).
of Eqs. (8)–(10). The relationship between horizontal force and drive torque vs. slip ratio is included and is depicted in Fig. 10. As shown in Fig. 10, the distributions of horizontal force and drive torque vs. slip ratio can be approximately fitted with polynomial curves shown by solid and dashed curves for both experimental datum and DEM simulation results, respectively. Compared with the data points obtained from laboratory experiments, the DEM simulation results generally agree well with that of the experiment in spite of a slight difference in drive torque between experiment and simulation. As a result, we can draw the conclusion that DEM simulation with the wheel–soil DEM model proposed here is suitable for numerical analysis of the trafficability of lunar mare terrain under terrestrial conditions; the DEM simulation can quantitatively and qualitatively estimate terrain trafficability without considering the environmental effects of the moon. In Section 3.2, we will further discuss the effect of the lunar environment on lunar mare terrain trafficability and modify the initial DEM model for regolith particles. 3.2. Trafficability analysis for lunar mare terrain 3.2.1. Modification of DEM model The lunar surface environment is characterized by reduced gravity, ultra-high vacuum (1012 Torr), and a large temperature difference between day and night or the shaded region of the moon, as well as the intense and unim-
peded radiation from solar and high-energy cosmic and galactic rays in the absence of atmosphere shielding. The lunar environmental affects the microdynamic behavior of regolith particles on the moon subsurface, and generates various non-contact forces among regolith particles, such as Van der Waals forces, electrostatic force, and gravitational force. The Van der Waals force is a short-range molecular attraction. The force range falls off as 1/R6 with distance, that is to say, it is only important for extremely small particles (<1 lm) that are extremely close to very smooth surfaces (<1010 m) [21]. For regolith particles with any noticeable size(>1 lm), which is a significant portion of lunar regolith, particles will not experience this force since the range of force is beyond the Van der Waals range. Consequently, we neglect the influence of Van der Waals forces in our DEM analysis. Electrostatic force is a significant non-contact force on regolith particles that is responsible for the microdynamic behavior of regolith particles. With virtually no atmospheric shielding, subsurface regolith particles are radiated by solar electromagnetic radiation, solar wind plasma, and high-energy galactic cosmic rays, and consequently, part of them are charged. Considering that the lunar rover vehicle merely works during the lunar day specified by the Chinese Lunar Surface Exploration mission, we only discuss the influence of electrostatic force generated on the lunar lit side. On the directly illuminated lunar surface, subsurface regolith particles are mostly charged by photoelectric emis-
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
169
where md and qd are mass and density of the particle, respectively, and gm is the gravitational acceleration on the moon surface. In this study, we modify our particle model by introducing electrostatic force and gravitational force into the initial mechanical model of DEM. The required environmental parameters are summarized in Table 4 [22]. Based on Newton’s second law, the dynamic equation for the ith element of a soil particle can be written by 8 dvi P c g E > < mi dt ¼ j F ij þ F i þ F i ð13Þ P dx > M ij : I i dt i ¼ j
F cij ;
where M ij are contact force and torque acting on element i by element j, F gi ; F Ei are gravitational force and electrostatic force of element i, and mi, vi, xi, and Ii are mass, translational velocity, angular velocity, and moment of inertia of element i, respectively.
Fig. 10. Experimental and simulation results.
sions due to solar UV and X-ray radiation. Furthermore the high-energy photoelectrons emitted from the directly illuminated material surface create a large electrostatic field characterized by the Debye length near the lunar surface. In this case, the electrostatic force is produced on all of the charged particles. For a spherical particle with radius Rd, the electrostatic force acting on it can be expressed by Ec W E ð11Þ F E ¼ Q E 4pe0 Rd e where Q is the charge on particle, E is the surrounding electric field strength, e is the element charge, Ec is the energy of the solar spectrum (which is typically dominated by the Lyman-alpha emission), W is the work function of the particle, and e0 is the permittivity of free space. In addition to the electrostatic force, the gravitational force is also important to particle dynamic behavior and terrain trafficability. For a spherical particle under a lunar gravity field, the gravitational force acting on the particle can be written by 4p 3 R q gm ð12Þ F g ¼ md g m 3 d d
3.2.2. DEM simulation results A series of analyses are performed for wheel–terrain interaction at various slip ratios using the DEM model. The terrain deformation and trafficability parameters are estimated by quantitative and qualitative results, as shown in Figs. 11–15. Figs. 11 and 12 depict the terrain deformation and contact force distribution of particles beneath the wheel under static sinkage with the initial and modified mechanical DEM model. The DEM simulation results show that the maximal values of static sinkage for the initial and modified model are 37.8 mm and 26.0 mm respectively, and the network of contact forces beneath wheel is more loosely distributed for the modified mechanical model results than that obtained from the initial DEM model. This phenomenon can be explained: reducing the gravitational acceleration from 9.8 m/s2 to 1.63 m/s2 decreases the wheel load (unbalanced force in vertical direction), and results in a relatively small contact pressure produced at the wheel–terrain interface, as represented by the relatively thin contact force line distribution in the particulate system of the terrain. Although lunar gravity and repulsive electrostatic forces acting on regolith particles decrease the bearing strength of regolith due to the more loosely arranged soil structure, the major contributor responsible for wheel sinkage is wheel load; as a result, less sinkage is observed in the DEM simulation with the modified model.
Table 4 Lunar surface environmental parameters. Lunar gravity, gm (m/s2) Element charge, e (C) Angle from subsolar point, h (°) Surface electric field, E (V/m) Energy of the solar spectrum, Ec (eV) Work function of particle, W (eV) Permittivity of free space, e0 (F/m)
1.63 1.602 10–19 0–6 +9.9 10.2 (k = 121.6 nm) 5.8 8.854 1012
170
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
Fig. 11. Network of contact force (on the Earth).
Fig. 12. Network of contact force (on the moon).
Fig. 13. Trajectories of soil particles and dynamic sinkage of wheel at 0.5 slip ratio calculated by the modified model.
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
Fig. 14. Dynamic sinkage of wheel under the moon environment.
171
gradually is enhanced and reaches its maximal value at the 1.0 slip ratio. Furthermore, at low slip ratio, the dynamic sinkage of wheel tends to reach a stationary state, whereas at a high slip ratio, such as 1.0, the dynamic sinkage always increases and reaches its maximal value, since the wheel motor torque is less than the drive torque. Based on Eqs. (8)–(10) and the modified mechanical model, we calculated the horizontal force and drive torque at various slip ratios to quantitatively estimate the trafficability of lunar mare terrain, as shown in Fig. 15. Compared with the result obtained from the initial model shown in Fig. 10, it is noted that the lunar environment affects the value of horizontal force at any slip ratio, and it is generally smaller than that obtained by experiment and simulation with the initial model under the same traveling conditions. Consequently, we draw the conclusion that the trafficability of in situ lunar mare terrain is worse than what we predict by laboratory experiment, and the wheel requires less drive torque on the moon surface. In future work, it will be necessary to carefully design and estimate lunar rover locomotion while taking lunar environmental effects into account. 4. Conclusions
Fig. 15. Trafficability parameters of lunar mare terrain obtained by modified model.
Fig. 13 plots the terrain deformation and trajectories of soil particles with a wheel running at a 0.5 slip ratio. It is observed that soil particles beneath wheel are dug up and sheared by the lugs and then move backward and upward with the rotation of the wheel. Some of the regolith particles leave the terrain surface with a certain acceleration due to the wheel disturbance. After a period of time, the levitated particles fall back to the terrain surface under gravitational force and form a hump behind the wheel. Additionally, the obvious dynamic sinkage occurs at the 0.5 slip ratio as the wheel travels along the horizontal direction with a given translational velocity vx and angular velocity x, as shown in Fig. 13. Fig. 14 shows the dynamic sinkage of the wheel under a 11.8-kg mass load as it horizontally traverses from 0 to 200,000 timesteps at 0.0, 0.5, and 1.0 slip ratios. It is observed that the magnitude of dynamic sinkage varies over a small range with the wheel running at 0.0 slip; with the increase of slip ratio, the dynamic sinkage of the wheel
We have developed a strategy for producing an irregularly shaped particulate system by discrete element modeling based on physical and mechanical properties of lunar regolith. The model parameters are determined by means of dimensional analysis and biaxial test simulation with an improved boundary condition. A series of experiments and simulations were performed under terrestrial conditions to validate the applicability of the wheel–terrain model established here. Taking lunar environmental effects into account, we modified the initial mechanical model by introducing non-contact forces attributed to the lunar environment into the model. With the modified model, the wheel–terrain interaction is qualitatively and quantitatively analyzed, and the terrain deformation and dynamic sinkage of wheel as well as the trajectories of regolith particles are obtained. Using the relationship between microscopic DEM results and macroscopic traction performance parameters, the trafficability of the in situ lunar mare terrain is estimated. The results show that trafficability is worse on the moon than what we predict by laboratory experiment, and the wheels require less drive torque on the moon surface. Therefore, in future work, it is necessary to carefully design and estimate lunar rover locomotion while taking lunar environmental effects into account. Acknowledgements This study is funded by the China Academy of Space Technology and Research Fund for the Doctoral Program of Higher Education of China (No. 20070006012); we are
172
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
sincerely grateful to them. We also express our sincere gratitude to Prof. F. Gao and Dr. P. Sun for their assistance with experiment validation. References [1] Apollo 14 mission report. National Aeronautics and Space Administration, Manned Spacecraft Center, Houston, Texas; 1971 [MSC-04112]. [2] Li W, Gao F, Xu GY. Design and trafficability study of flexible wheel for planetary exploration rover. J Beijing Inst Technol 2007;16(3): 279–83. [3] Carrier WD, Olhoeft GR, Mendell W. Physical properties of the lunar surface. In: Heiken G, Vaniman D, French BM, editors. Lunar source book. Cambridge: Cambridge University Press; 1991. [4] Tao JG, Wang L, Deng ZQ, et al. Mechanical analysis and measurement of parameters of wheel–soil interaction for a lunar rover. J Phys: Conf Ser 2006;48:1222–6. [5] Patel N, Scott GP, Ellery A. Application of Bekker theory for planetary exploration through wheeled, tracked and legged vehicle locomotion. In: Space 2004 conference and exhibit, 28–30 September 2004, San Diego, California [AIAA 2004-6091]. [6] Melzer KJ. Performance of the Boeing LRV wheel in a lunar soil simulant-effcts of speed, wheel load, and soil [R]. Technical report M71-10; 1971. [7] Lizuka K, Sato Y, Kuroda Y, et al. Study on wheel of exploration robot on sandy terrain. In: Intelligent robots and systems, 2006 IEEE/ RSJ international conference. Beijing, China; 2006.p. 4272–77. [8] Kuroda Y, Teshima T, Sato Y, et al. Mobility performance evaluation of planetary rover with similarity model experiment. In: Robotics and automation, 2004. Proceedings. ICRA ’04. 2004 IEEE international conference, vol. 2. New Orleans, LA, USA; 2004. p. 2098–103. [9] Khot LR, Salokhe VM, Jayasuriya HPW, et al. Experimental validation of distinct element simulation for dynamic wheel–soil interaction. J Terramech 2007;44:429–37.
[10] Nakashima H, Fujii H, Oida A, et al. Parametric analysis of lugged wheel performance for a lunar microrover by means of DEM. J Terramech 2007;44:153–62. [11] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. J Geotech 1979;29(1):47–65. [12] Asaf Z, Rubinstein D, Shmulevich I. Evaluation of link-track performances using DEM. J. Terramech 2006;43:141–61. [13] Shmulevich I, Asaf Z, Rubinstein D. Interaction between soil and a wide cutting blade using the discrete element method. Soil Tillage Res 2007;97:37–50. [14] Gromov V. Physical and mechanical properties of lunar and planetary soils. Earth Moon Planets 1999;80:51–72. [15] Liu Y, Schnare D, Park JS, et al. Shape analyses of lunar dust particles for astronaut toxicological studies. Lunar Planet Sci XXXVIII; 2007 [1383.pdf]. [16] Carrier DW. Lunar soil grain size distribution. Moon 1973;6:250–63. [17] Cundall PA, Strack ODL. Modeling of microscopic mechanics in granular material. In: Jenkins, JT, Satake M, editors. Mechanics of granular materials: new models and constitutive relations. Amsterdam: Elsevier; 1983. p. 113–49. [18] Liao CL, Chang TP, Young DH, et al. Stress strain relationship for granular materials based on the hypothesis of best fit. Int J Solids Struct 1997;34(31):4087–100. [19] Oida A, Schwanghart H, Ohakubo S. Effect of tire lug cross section on tire performance simulated by distinct element method. In: Proceedings of the 13th international conference of the ISTVS, Munich; 1999. p. 345–52. [20] Li W, Gao F, Jia Y. Tractive performance analysis on radially deployable wheel configuration of lunar rover vehicle by discrete element method. Chinese J Mech Eng 2008;21(4):66–72. [21] Buhler CR, Calle CI, Clements JS, et al. Test method for in situ electrostatic characterization of lunar dust. In: Aerospace Conference, 2007. IEEE; 2007. p. 1–19. [22] Stubbs TJ, Vondrak RR, Farrell WM. A dynamic fountain model for lunar dust. Adv Space Res 2006;37:59–66.
Journal of Terramechanics Journal of Terramechanics 47 (2010) 161–172 www.elsevier.com/locate/jterra
Trafficability analysis of lunar mare terrain by means of the discrete element method for wheeled rover locomotion Wen Li, Yong Huang *, Yi Cui, Sujun Dong, Jun Wang School of Aeronautical Science and Engineering, Beijing University of Aeronautics and Astronautics, Beijing 100191, China Received 15 February 2009; received in revised form 14 August 2009; accepted 11 September 2009
Abstract An irregularly shaped particulate system for simulation of lunar regolith is developed using discrete element modeling based on the fractal characteristics, particle shape, and size distribution of returned Apollo-14 samples. The model parameters are determined by dimensional analysis and biaxial test simulation with an improved boundary condition. Under terrestrial conditions, the trafficability of lunar mare terrain is estimated in terms of wheel–terrain interaction by experiment and simulation in order to validate the applicability of the wheel–terrain model employed here. The results show that the discrete element method combined with the wheel–terrain model is sufficiently accurate for mare terrain trafficability analysis without consideration of lunar environmental effects. To predict the trafficability of in situ lunar mare terrain, the non-contact forces attributed to the lunar surface environment are discussed and the initial mechanical model of discrete elements is modified by introduction of lunar gravitational force as well as electrostatic force. In the modified model, wheel–terrain interaction is analyzed under the same travel conditions as that of the experiment. The result shows the trafficability of the in situ lunar mare terrain is worse than that obtained by experiment and simulation with the initial model according to the value of horizontal force at any slip ratio. However, the wheel requires less drive torque on the moon than that on the earth. An explanation for these phenomena may be that lunar subsurface regolith particles are arranged in a looser manner under local environmental effects that effectively decrease the bearing and shearing strength of regolith. Ó 2009 ISTVS. Published by Elsevier Ltd. All rights reserved.
1. Introduction China is planning to send an unmanned lunar rover vehicle to the surface of the moon to make observations and measure the physical and mechanical properties of in situ regolith material during the second phase of the Chang’E Lunar Exploration Project. In consideration of mass and power consumption constraints, the lunar rover vehicle was generally designed with wheeled locomotion in previous and current lunar exploration missions. From prior experience reported by the Apollo and Luna missions, the lunar mare is a favorite landing site because of its relatively flat topographical condition and absence *
Corresponding author. Address: Mailbox 505, School of Aeronautical Science and Engineering, Beijing University of Aeronautics and Astronautics, 37 Xueyuan Road, Haidian District, Beijing 100191, China. Tel.: +86 10 8233 8498; fax: +86 10 8231 6654. E-mail address: [email protected] (Y. Huang).
of rocks. The Apollo-14 mission report [1] demonstrated that the lunar mare surface was generally covered with a loose mantle of material generated by micro-meteoritic impacts, and the population of rocks in the Fra Mauro area was much less than 0.5% of the total area. As a result, the soft-terrain trafficability of the lunar mare has become the foremost factor limiting the traction performance of a lunar rover during lunar surface exploration missions. Based on terramechanics, soft-terrain trafficability is defined as the capability of a soil-like material to support the traverse of vehicular traffic from one place to another [2]. It is related to the mechanical properties of the terrain on the macroscopic scale and the dynamic behavior of soillike particles on a granular scale. Current scientists have developed various methods to meet the requirements for rover design and mobility evaluation as well as lunar soft-terrain trafficability analysis. To predict the lunar rover’s mobility performance, Carrie [3] listed the recom-
0022-4898/$36.00 Ó 2009 ISTVS. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.jterra.2009.09.002
162
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
mended lunar soil trafficability parameters used for Bekker’s empirical equation. Tao et al. [4] and Patel et al. [5] analyzed the interaction between a rigid wheel and soft-terrain based on Terzaghi’s passive soil pressure theory and Bekker’s formula, respectively. Experimental validation was conducted using a specially designed and manufactured apparatus. Melzer [6] experimentally measured the Boeing LRV wheel’s performance by a dynamometer system and lunar soil simulant. Lizuka [7] studied the performance of a lugged wheel on sandy terrain using a model experiment in the laboratory. Kuroda et al. [8] evaluated the mobility performance of a planetary rover with a similar model experiment under low-gravity flight conditions. Khot et al. [9] simulated the dynamic wheel–soil interaction via the distinct element method and verified the simulation results with laboratory experiments. Nakashima et al. [10] numerically analyzed the influence of wheel configuration on traction performance for a lunar microrover by means of the discrete element method. In comparison of the existing methods for trafficability estimation, we find that: (i) The continuum mechanics method, based on the soil constitutive law and failure pattern, is unable to cover the flow of soil particles; therefore, it is difficult to analyze the trajectories of soil particles and explain soil deformation and failure mechanisms on the granular scale. (ii) The laboratory experiment measurement method can directly acquire trafficability parameters with a highaccuracy transducer system; however, lunar environmental simulation is costly and difficult to realize. (iii) The discrete or distinct element method (DEM), which was originally proposed by Cundall and Strack [11] to model the behavior of granular materials, seems to be a promising method for numerical analysis of the trafficability of lunar mare terrain, without the limitations of the methods mentioned above. In such an approach, the forces acting on an individual particle can be calculated using a well-founded contact reaction model. The trajectories and velocity vectors of particles can be determined based on Newton’s second law. With the development of a contact detection algorithm and discrete particle modeling method, DEM has become an ideal tool for terramechanics-based research. So far, however, the soil particles are commonly simplified as circular discrete element models in terramechanicsbased studies to reduce computational time during DEM simulation. Although Asaf [12,13] has developed a nonspherical DEM model in his research, such a model is far from realistic. Herein, we aim to develop a strategy for producing an irregularly shaped particulate system based on scanning electron microscopy (SEM) images and grain size distributions of lunar regolith and simulant, as described in Section 2. In addition, the mechanical model (contact reaction
model) and its corresponding microparameters are determined based on the physical and mechanical properties of lunar regolith. In Section 3, a series of DEM simulations with respect to validated experiments are performed to evaluate the applicability of the DEM model specified here. In Section 4, we discuss the influence of the lunar environment on dynamic behavior of regolith particles, and modify the initial DEM model by introducing non-contact forces attributed to lunar environmental effects. With the modified model, we estimate the trafficability of lunar mare terrain and explain the microscopic mechanism of terrain deformation and failure through trajectories and networks of contact force of regolith particles beneath the wheel. 2. DEM model for lunar regolith particles Lunar mare terrain trafficability is governed by the mechanical properties of in situ regolith and is mainly affected by geometric features and dynamic behavior of subsurface soil particles. Herein, we construct the geometric model and mechanical model of DEM for lunar mare regolith with high fidelity. 2.1. Geometric model for lunar regolith particles Based on knowledge of soil mechanics, the microscopically geometric parameters that influence the macromechanical properties of soil are particle shape, grain size distribution, and initial voidage. 2.1.1. Particle shape model Considering the fact that lunar regolith particles are generally irregular with re-entrant surfaces [14] and the irregularity of particle shape influences the interlock interactions among particles, we present an imaged-based modeling method to simulate the truly geometric features of lunar regolith particles. The required geometric feature parameters are acquired from SEM photography of returned lunar samples and corresponding simulants by using ImageJ software [15] and fractal calculation. The results show that the fractal dimension with perimeter–area relation for lunar regolith is 1.098, and particles with different sizes have a self-similar fractal behavior. Consequently, we build our DEM model with four typical particle shapes shown in Fig. 1. As is shown in Fig. 1, any irregular-shaped particle can be modeled by overlapping elements with a center element O and several angular elements Bi (i = 1, 2, 3, 4, . . .). The number of angular elements is determined by the complexity of the geometric features for a realistic particle. All elements are parametrically modeled in the Cartesian coordinate system, and both positions and radii are variables. To reduce the contact detection time, every particle shape model is clumped and behaves as a rigid body wherein the mutual contact check of constituent elements is neglected during DEM simulation.
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
Fig. 1. Particle shape model of DEM for lunar regolith.
2.1.2. Grain size distribution model The grain size distribution (GSD) by wt.% of returned lunar samples from the Apollo missions was experimentally measured and published in the literature [16]. Fig. 2 depicts the upper and lower bounds of the GSD for returned Apollo-14 samples as well as the GSD of the lunar mare simulant used in our experimental validation by solid and dashed curves, respectively. Here, we develop an effective method to generate the particulate system of the DEM with the best-fit GSD of the lunar mare simulant shown in Fig. 2. The strategy of producing irregularly shaped particles with a specified GSD is described as follows: (1) Generate ball elements with radii R and position (x0, y0) according to the GSD specified here, where x0, y0, R are variables in the Cartesian coordinate system. (2) Set contact stiffness of ball elements with high value (1.0 108 Pa) and without friction and gravity.
Fig. 2. Grain size distribution curves.
163
164
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
(3) Balance the particulate system under inter-element contact reactions in a virtual container constructed by rigid wall elements with higher stiffness than that of ball elements. (4) Reduce the radii of ball elements to R0 while keeping initial locations constant, then add angular elements based on the particle shape model given in Fig. 1. (5) Set best-fit values to DEM model parameters, such as density, spring stiffness, damping coefficient, and Coulomb friction coefficient, and balance the particulate system. Finally, we obtain the most appropriate model for our DEM analysis. 2.2. Mechanical model of lunar regolith particles The dynamic behavior of any particle is governed by the forces and torques acting upon it. As described by Apollo mission reports, there are no liquid water and other adhesive materials in lunar regolith; the cohesion and internal friction angle of regolith are mainly attributed to the interlock interactions among irregularly shaped particles. In our DEM model, we initially simplify the mechanical model of ball elements as a linear spring–dashpot system with Coulomb friction and non-tensile links, as shown in Fig. 3. As shown in Fig. 3, the inter-element contact force F can be decomposed into the normal force Fn and tangential force Fs. The normal force is determined by the elastic part F ne due to the normal spring, and the inelastic part F nd due to the normal damper and non-tensile link, i.e. As shown in Fig. 3, the inter-element contact force F can be decomposed into the normal force Fn and tangential force Fs. The normal force is determined by the elastic part F ne due to the normal spring, and the inelastic part F nd due to the normal damper and non-tensile link, i.e., kF n k ¼ F ne þ F nd
ð1Þ
F ne F nd
ð2Þ
¼ k n urn ¼ k n ðRj þ Ri dÞ ¼ cn vrn ¼ cn ðu_ j u_ i Þ n
ð3Þ
where kn is normal spring stiffness, urn is the penetration of two elements at the interface, d is the distance between element centers, Ri and Rj are element radii, cn is the normal damping coefficient, vrn is normal relative velocity, u_ i and u_ j
are the translational velocities of elements, and n is the unit normal vector. The tangential force includes tangential spring force, tangential damping force, and Coulomb friction force. It can be expressed by 8 old 1 Dt > < F s ðks 1 2csDtÞvrs Dt ðkF new s k 6 ljF n jÞ new ðks þ2cs Þ ð4Þ Fs ¼ > new s : ljF n j F new ðkF k > ljF jÞ new n s kF k s
where kF new s k is the norm of the current tangential contact force, F old is the tangential force in the previous timestep, s ks is the tangential spring stiffness, cs is the tangential damping coefficient, vrs is the relative tangential velocity at the contact point, l is the Coulomb friction coefficient, and Dt is the timestep able to avoid disturbance by propagating from the particle to a further range. 2.3. Determination of DEM model parameters The DEM model parameters, such as spring stiffness, damping coefficient, and Coulomb friction coefficient, are difficult to determine by direct physical measurement. As a result, various indirect methods have been developed. Cundall and Strack obtained the model parameters by using Hertzian contact theory and Mindlin and Deresiewich’s approach [17]. Liao and Chang [18] developed a method for determining the bounds of the spring constant between two particles by assuming an elastic bulk material. Oida [19] and Asaf [12,13] used the trial-and-error method to determine model parameters in their DEM simulation. In previous research, we confirmed that the method of dimensional analysis combined with biaxial test simulation can obtain the best-fit parameters of the DEM model [20]; therefore, we determine the model parameters required in Eqs. (1)–(4) by an improved biaxial test model and similarity criteria. Based on the Buckingham G theorem and dimensional matrix analysis, we obtain four dimensionless parameters Gi (i = 1, . . . , 4) and the relationship between microscopic parameters of the DEM model and macro-mechanical property parameters. F ðP1 ; P2 ; P3 ; P4 Þ ¼ 0 ck s , E2
Fig. 3. Contact reaction model between two elements.
ð5Þ ck n , s2max
kn , ks
Here, P1 ¼ P2 ¼ P3 ¼ P4 = l, where c is the bulk density of soil, E is the elastic modulus, and smax is the shear strength of the soil. According to the triaxial compression test of the lunar simulant mentioned above, we construct an improved DEM model for biaxial test simulation shown in Fig. 4. In this model, two flexible bounds are built with overlapped ball elements to accurately simulate the deformation of a soil sample during axial loading. Assuming that the radius of constituent elements of flexible bounds is Rb in addition to the width W0 and height H0 of the soil sample, then the number of elements for each flexible bound can be calculated by
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
165
Table 1 Parameters for boundary condition of DEM. Soil sample dimensions, W0 H0 (mm) Confining pressure, r3 (kPa) Density of boundary elements, qb (kg/m3) Radius of boundary elements, Rb (mm) Normal stiffness of spring, kn (kPa) Tangential stiffness of spring, ks (kPa) Bond between elements, sb (N/m) Friction coefficient of element, lb Element number of each bound, n Stiffness of walls, k (kPa) Friction coefficient of walls, lw Velocity of bottom wall, vwy (mm/min)
39.1 80.0 26.0 100.0 1.0 1000.0 1000.0 1.0 1010 0.2 80 1.0 106 0.0 0.8
and the best-fit parameters of DEM model are obtained and listed in Table 2 based on the dimensional analysis result. For comparison, the axial strain vs. axial stress difference curves obtained from both experiment and simulation are plotted in Fig. 5. It can be observed from Fig. 5 that the DEM model and its best-fit parameters are able to accurately reflect the soil constitutive law for lunar regolith and the simulant, and the DEM biaxial simulation result agrees with the results of the laboratory experiment.
Fig. 4. DEM model for biaxial test simulation.
3. Trafficability analysis for lunar mare terrain n¼
H0 Rd
ð6Þ
To control the boundary condition and acquire the radial deformation of soil sample during DEM simulation, the constituent elements of both bounds are indexed with an ID number: for the left bound, the ID number varies from 1 to n, and for the right bound, the ID number varies from n + 1 to 2n. Furthermore, every element of the bounds is preloaded with an x-component force Fbx. Given that the triaxial test confining pressure is r3, the force imposed on every bound element can be written by F bx ¼
2Rb H 0 r3 n
ð7Þ
where the positive or negative sign on the force is consistent with the Cartesian coordinate system shown in Fig. 4. For the left bound’s constituent elements (ID = 1, . . . , n), Fbx > 0; for the right bound’s constituent elements (ID = n + 1, . . . , 2n), Fbx < 0. During the biaxial test simulation, two rigid wall elements are generated on the top and bottom sides of the soil sample. The top wall is fixed and functions as a load cell with which axial force Fwy can be acquired by summing up the contact reaction forces of the soil particle elements acting on the wall. The bottom wall is given by a constant translational velocity vwy on the y-axis in order to simulate the loading process. The parameters of the boundary condition for our biaxial test simulation are listed in Table 1. Using the lunar regolith DEM model (geometric model and mechanical model) and the boundary conditions of the biaxial test, a series of biaxial test simulations are performed
In terms of soft-terrain, trafficability is usually analyzed by the wheel–soil interaction under various slip ratios. In this section, we first perform a series of experiments under terrestrial conditions to validate the applicability of the wheel–terrain DEM model. Then, the DEM models for both wheel configuration and lunar mare terrain are developed. Third, a series of DEM simulations for wheel–terrain interaction are carried out with respect to experiments and to validate the DEM models. Finally, the influence of lunar environmental effects on trafficability is discussed. The non-contact forces, which are critical to the microdynamic behavior of lunar regolith particles, are introduced into the initial mechanical model of DEM for soil particles. With this modified model, lunar mare terrain trafficability is predicted with high fidelity. 3.1. Experimental validation 3.1.1. Experimental apparatus and conditions Experimental measurements of lunar mare terrain trafficability are performed with a specially designed and manTable 2 Best-fitted model parameters for lunar soil particles. Density of element, q (kg/m3) Radius of center element, R0 (mm) Friction coefficient, l Normal stiffness, kn (kPa) Tangential stiffness, ks (kPa) Damping coefficient ratio, c/ccrit Increment of timestep, Dt (s)
2770 0.25–0.625 0.18 5.0 104 5.0 104 0.7 8.5 106
166
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
sured during the experiment. Because the carriage is mounted on the framework of the soil bin by two horizontally fixed guiding axles, it can traverse smoothly along the horizontal axis shown in Fig. 6. In our experimental measurement of lunar mare terrain trafficability, a rigid wheel with thin lugs is adopted and the lunar mare terrain is reproduced by the lunar simulant in the soil bin. The mechanical property of simulated terrain in the soil bin is measured with the cutting ring sample method and simple shear test. The wheel configuration and travel conditions are summarized in Table 3.
Fig. 5. Axial strain vs. axial stress difference curves.
ufactured wheel–terrain testbed for the study of planetary rover locomotion performance. The testbed primarily consists of a soil bin and lunar simulants, a carriage assembly, wheel assembly, and the data acquisition/control/processing system shown in Fig. 6. The soil bin is constructed from a metallic frame and transparent sidewalls (tempered glass) to aid visibility of soil particle flow beneath the tested wheel. The soil bin is designed to be large enough to eliminate boundary effects during the experiment. The wheel assembly is mainly composed of the tested wheel, a Maxon DC motor with controller, and six-axial force/torque sensors, and it is able to move freely along the vertical axis of the carriage. During experimental measurements, the tested wheel is driven and controlled by the motor module while the forces and torques exerted by the ground on the tested wheel are acquired by the six-axial force/torque sensor. The carriage assembly is equipped with a potentiometer, photoelectric encoder, and the driven motor module, with which the sinkage and translational velocity of the wheel can be mea-
3.1.2. DEM model for wheel–soil interaction In the DEM model for wheel–soil interaction analysis, the tested wheel configuration is modeled by clumped ball elements. The travel conditions of the wheel with respect to the experiment are input to the DEM model. To quantitatively estimate the soft-terrain trafficability by DEM simulation, we propose a method used to build the relationship between the microscopic results of DEM and the macromechanical parameters of trafficability based on coordinate transformation and algebraic summation. Fig. 7 shows a contact pair between wheel element i and soil element j as they contact each other, and also gives the schematic of forces exerted on element i. As shown in Fig. 7, element i is acted on by normal force F N i and tangential force F Si ; which can be calculated using Eqs. (1)–(4). By means of coordinate transformation, we obtain the contact reaction forces of element i in the Cartesian coordinate system. ( S F xi ¼ F N i sin hi F i cos hi ð8Þ S F yi ¼ F N i cos hi þ F i sin hi where hi = Arc tan [(x xj)/(y yj)]. Supposing that the number of constituent elements of the wheel model is Nw, then the horizontal force, wheel load, and drive torque of wheel can be expressed as FH ¼
Nw X
F xi
i¼1
T ¼
Nw X
F Si
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxi x0 Þ2 þ ðy i y 0 Þ2
ð9Þ ð10Þ
i¼1
Table 3 Parameters for experiment.
Fig. 6. Experimental apparatus.
Soil bin dimensions, L W H (mm) Cohesion of test terrain, c (kPa) Internal friction angle of test terrain, u (°) Initial void ratio of lunar simulant, e0 Wheel diameter without lug, D (mm) Lug dimensions, LT LH (mm) Wheel wide, B (mm) Number of lug on wheel, NL Vertical load on wheel, W (kg) Slip ratio, i (%) Angular velocity of wheel, x (rpm) Gravitational acceleration, g (m/s2)
3.0 0.8 0.6 0.28 39.95 0.85 220.0 3.0 25.0 160.0 18 11.8 0–100 5.0 9.8
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
167
Fig. 7. Contact reaction force on individual wheel element.
In this study, the lunar mare terrain is reproduced by an irregularly shaped particulate system (described in Section 2) in a virtual container. The container is 1.0-m long and 0.2-m deep, and constructed with rigid walls. To improve computational efficiency, we filled the soil particles with various scales. Namely, the soil particle size is enlarged with the increase in container depth. The equivalent diameter of a soil particle in our DEM model varies from 0.5 mm to 2.0 mm. The process of DEM simulation is divided into four phases: (i) lunar mare terrain reproduction with respect to the experiment; (ii) construction of a wheel model above the terrain surface, without any contact between wheel and terrain; (iii) wheel sinkage under a vertical load imposed on its center; (iv) wheel–terrain interaction simulation with specified travel conditions. In the wheel–terrain DEM model, there are 166, 135 discrete elements used, and the increment of the timestep is set to 8.85 106 s for the sake of stability and equilibrium of simulation. 3.1.3. Experiment and simulation results The trafficability of lunar mare terrain is analyzed by experiments and simulations in terms of wheel–soil interaction under specified conditions listed in Table 3. Figs. 8 and
9 show terrain deformation with a lugged wheel for 0.0 and 0.8 slip ratios by experiment and simulation, respectively. It is obvious that the terrain surface forms clear ruts along the travel direction as the wheel moves without slip; however, at higher slip (i = 0.8), the upper layer particles of the terrain shift backward and heap up behind the wheel due to the lug digging and shearing. Soil dilatancy or shear deformation occurs at that moment as well, as shown in Fig. 8b. Corresponding DEM simulations are performed with the wheel–terrain DEM model. The trajectories of soil particles and terrain deformation at 0.0 and 0.8 slip ratios are analyzed qualitatively and demonstrated in Fig. 9. Compared with experimental results, DEM simulation shows the same tendency in terrain deformation, that is to say, the terrain surface forms clear and wavy ruts at zero slip while soil dilatancy occurs behind the wheel at the 0.8 slip ratio. Consequently, we draw the conclusion that the DEM model for wheel–soil interaction can qualitatively estimate the terrain deformation and trajectories of subsurface soil particles with sufficient accuracy. To quantitatively validate the applicability of our DEM model and its parameters, the horizontal force and drive torque acting on the wheel center are calculated by means
Fig. 8. Experimental result for terrain deformation (a: i = 0; b: i = 0.8).
168
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
Fig. 9. DEM simulation for terrain deformation (i = 0.0 and 0.8).
of Eqs. (8)–(10). The relationship between horizontal force and drive torque vs. slip ratio is included and is depicted in Fig. 10. As shown in Fig. 10, the distributions of horizontal force and drive torque vs. slip ratio can be approximately fitted with polynomial curves shown by solid and dashed curves for both experimental datum and DEM simulation results, respectively. Compared with the data points obtained from laboratory experiments, the DEM simulation results generally agree well with that of the experiment in spite of a slight difference in drive torque between experiment and simulation. As a result, we can draw the conclusion that DEM simulation with the wheel–soil DEM model proposed here is suitable for numerical analysis of the trafficability of lunar mare terrain under terrestrial conditions; the DEM simulation can quantitatively and qualitatively estimate terrain trafficability without considering the environmental effects of the moon. In Section 3.2, we will further discuss the effect of the lunar environment on lunar mare terrain trafficability and modify the initial DEM model for regolith particles. 3.2. Trafficability analysis for lunar mare terrain 3.2.1. Modification of DEM model The lunar surface environment is characterized by reduced gravity, ultra-high vacuum (1012 Torr), and a large temperature difference between day and night or the shaded region of the moon, as well as the intense and unim-
peded radiation from solar and high-energy cosmic and galactic rays in the absence of atmosphere shielding. The lunar environmental affects the microdynamic behavior of regolith particles on the moon subsurface, and generates various non-contact forces among regolith particles, such as Van der Waals forces, electrostatic force, and gravitational force. The Van der Waals force is a short-range molecular attraction. The force range falls off as 1/R6 with distance, that is to say, it is only important for extremely small particles (<1 lm) that are extremely close to very smooth surfaces (<1010 m) [21]. For regolith particles with any noticeable size(>1 lm), which is a significant portion of lunar regolith, particles will not experience this force since the range of force is beyond the Van der Waals range. Consequently, we neglect the influence of Van der Waals forces in our DEM analysis. Electrostatic force is a significant non-contact force on regolith particles that is responsible for the microdynamic behavior of regolith particles. With virtually no atmospheric shielding, subsurface regolith particles are radiated by solar electromagnetic radiation, solar wind plasma, and high-energy galactic cosmic rays, and consequently, part of them are charged. Considering that the lunar rover vehicle merely works during the lunar day specified by the Chinese Lunar Surface Exploration mission, we only discuss the influence of electrostatic force generated on the lunar lit side. On the directly illuminated lunar surface, subsurface regolith particles are mostly charged by photoelectric emis-
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
169
where md and qd are mass and density of the particle, respectively, and gm is the gravitational acceleration on the moon surface. In this study, we modify our particle model by introducing electrostatic force and gravitational force into the initial mechanical model of DEM. The required environmental parameters are summarized in Table 4 [22]. Based on Newton’s second law, the dynamic equation for the ith element of a soil particle can be written by 8 dvi P c g E > < mi dt ¼ j F ij þ F i þ F i ð13Þ P dx > M ij : I i dt i ¼ j
F cij ;
where M ij are contact force and torque acting on element i by element j, F gi ; F Ei are gravitational force and electrostatic force of element i, and mi, vi, xi, and Ii are mass, translational velocity, angular velocity, and moment of inertia of element i, respectively.
Fig. 10. Experimental and simulation results.
sions due to solar UV and X-ray radiation. Furthermore the high-energy photoelectrons emitted from the directly illuminated material surface create a large electrostatic field characterized by the Debye length near the lunar surface. In this case, the electrostatic force is produced on all of the charged particles. For a spherical particle with radius Rd, the electrostatic force acting on it can be expressed by Ec W E ð11Þ F E ¼ Q E 4pe0 Rd e where Q is the charge on particle, E is the surrounding electric field strength, e is the element charge, Ec is the energy of the solar spectrum (which is typically dominated by the Lyman-alpha emission), W is the work function of the particle, and e0 is the permittivity of free space. In addition to the electrostatic force, the gravitational force is also important to particle dynamic behavior and terrain trafficability. For a spherical particle under a lunar gravity field, the gravitational force acting on the particle can be written by 4p 3 R q gm ð12Þ F g ¼ md g m 3 d d
3.2.2. DEM simulation results A series of analyses are performed for wheel–terrain interaction at various slip ratios using the DEM model. The terrain deformation and trafficability parameters are estimated by quantitative and qualitative results, as shown in Figs. 11–15. Figs. 11 and 12 depict the terrain deformation and contact force distribution of particles beneath the wheel under static sinkage with the initial and modified mechanical DEM model. The DEM simulation results show that the maximal values of static sinkage for the initial and modified model are 37.8 mm and 26.0 mm respectively, and the network of contact forces beneath wheel is more loosely distributed for the modified mechanical model results than that obtained from the initial DEM model. This phenomenon can be explained: reducing the gravitational acceleration from 9.8 m/s2 to 1.63 m/s2 decreases the wheel load (unbalanced force in vertical direction), and results in a relatively small contact pressure produced at the wheel–terrain interface, as represented by the relatively thin contact force line distribution in the particulate system of the terrain. Although lunar gravity and repulsive electrostatic forces acting on regolith particles decrease the bearing strength of regolith due to the more loosely arranged soil structure, the major contributor responsible for wheel sinkage is wheel load; as a result, less sinkage is observed in the DEM simulation with the modified model.
Table 4 Lunar surface environmental parameters. Lunar gravity, gm (m/s2) Element charge, e (C) Angle from subsolar point, h (°) Surface electric field, E (V/m) Energy of the solar spectrum, Ec (eV) Work function of particle, W (eV) Permittivity of free space, e0 (F/m)
1.63 1.602 10–19 0–6 +9.9 10.2 (k = 121.6 nm) 5.8 8.854 1012
170
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
Fig. 11. Network of contact force (on the Earth).
Fig. 12. Network of contact force (on the moon).
Fig. 13. Trajectories of soil particles and dynamic sinkage of wheel at 0.5 slip ratio calculated by the modified model.
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
Fig. 14. Dynamic sinkage of wheel under the moon environment.
171
gradually is enhanced and reaches its maximal value at the 1.0 slip ratio. Furthermore, at low slip ratio, the dynamic sinkage of wheel tends to reach a stationary state, whereas at a high slip ratio, such as 1.0, the dynamic sinkage always increases and reaches its maximal value, since the wheel motor torque is less than the drive torque. Based on Eqs. (8)–(10) and the modified mechanical model, we calculated the horizontal force and drive torque at various slip ratios to quantitatively estimate the trafficability of lunar mare terrain, as shown in Fig. 15. Compared with the result obtained from the initial model shown in Fig. 10, it is noted that the lunar environment affects the value of horizontal force at any slip ratio, and it is generally smaller than that obtained by experiment and simulation with the initial model under the same traveling conditions. Consequently, we draw the conclusion that the trafficability of in situ lunar mare terrain is worse than what we predict by laboratory experiment, and the wheel requires less drive torque on the moon surface. In future work, it will be necessary to carefully design and estimate lunar rover locomotion while taking lunar environmental effects into account. 4. Conclusions
Fig. 15. Trafficability parameters of lunar mare terrain obtained by modified model.
Fig. 13 plots the terrain deformation and trajectories of soil particles with a wheel running at a 0.5 slip ratio. It is observed that soil particles beneath wheel are dug up and sheared by the lugs and then move backward and upward with the rotation of the wheel. Some of the regolith particles leave the terrain surface with a certain acceleration due to the wheel disturbance. After a period of time, the levitated particles fall back to the terrain surface under gravitational force and form a hump behind the wheel. Additionally, the obvious dynamic sinkage occurs at the 0.5 slip ratio as the wheel travels along the horizontal direction with a given translational velocity vx and angular velocity x, as shown in Fig. 13. Fig. 14 shows the dynamic sinkage of the wheel under a 11.8-kg mass load as it horizontally traverses from 0 to 200,000 timesteps at 0.0, 0.5, and 1.0 slip ratios. It is observed that the magnitude of dynamic sinkage varies over a small range with the wheel running at 0.0 slip; with the increase of slip ratio, the dynamic sinkage of the wheel
We have developed a strategy for producing an irregularly shaped particulate system by discrete element modeling based on physical and mechanical properties of lunar regolith. The model parameters are determined by means of dimensional analysis and biaxial test simulation with an improved boundary condition. A series of experiments and simulations were performed under terrestrial conditions to validate the applicability of the wheel–terrain model established here. Taking lunar environmental effects into account, we modified the initial mechanical model by introducing non-contact forces attributed to the lunar environment into the model. With the modified model, the wheel–terrain interaction is qualitatively and quantitatively analyzed, and the terrain deformation and dynamic sinkage of wheel as well as the trajectories of regolith particles are obtained. Using the relationship between microscopic DEM results and macroscopic traction performance parameters, the trafficability of the in situ lunar mare terrain is estimated. The results show that trafficability is worse on the moon than what we predict by laboratory experiment, and the wheels require less drive torque on the moon surface. Therefore, in future work, it is necessary to carefully design and estimate lunar rover locomotion while taking lunar environmental effects into account. Acknowledgements This study is funded by the China Academy of Space Technology and Research Fund for the Doctoral Program of Higher Education of China (No. 20070006012); we are
172
W. Li et al. / Journal of Terramechanics 47 (2010) 161–172
sincerely grateful to them. We also express our sincere gratitude to Prof. F. Gao and Dr. P. Sun for their assistance with experiment validation. References [1] Apollo 14 mission report. National Aeronautics and Space Administration, Manned Spacecraft Center, Houston, Texas; 1971 [MSC-04112]. [2] Li W, Gao F, Xu GY. Design and trafficability study of flexible wheel for planetary exploration rover. J Beijing Inst Technol 2007;16(3): 279–83. [3] Carrier WD, Olhoeft GR, Mendell W. Physical properties of the lunar surface. In: Heiken G, Vaniman D, French BM, editors. Lunar source book. Cambridge: Cambridge University Press; 1991. [4] Tao JG, Wang L, Deng ZQ, et al. Mechanical analysis and measurement of parameters of wheel–soil interaction for a lunar rover. J Phys: Conf Ser 2006;48:1222–6. [5] Patel N, Scott GP, Ellery A. Application of Bekker theory for planetary exploration through wheeled, tracked and legged vehicle locomotion. In: Space 2004 conference and exhibit, 28–30 September 2004, San Diego, California [AIAA 2004-6091]. [6] Melzer KJ. Performance of the Boeing LRV wheel in a lunar soil simulant-effcts of speed, wheel load, and soil [R]. Technical report M71-10; 1971. [7] Lizuka K, Sato Y, Kuroda Y, et al. Study on wheel of exploration robot on sandy terrain. In: Intelligent robots and systems, 2006 IEEE/ RSJ international conference. Beijing, China; 2006.p. 4272–77. [8] Kuroda Y, Teshima T, Sato Y, et al. Mobility performance evaluation of planetary rover with similarity model experiment. In: Robotics and automation, 2004. Proceedings. ICRA ’04. 2004 IEEE international conference, vol. 2. New Orleans, LA, USA; 2004. p. 2098–103. [9] Khot LR, Salokhe VM, Jayasuriya HPW, et al. Experimental validation of distinct element simulation for dynamic wheel–soil interaction. J Terramech 2007;44:429–37.
[10] Nakashima H, Fujii H, Oida A, et al. Parametric analysis of lugged wheel performance for a lunar microrover by means of DEM. J Terramech 2007;44:153–62. [11] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. J Geotech 1979;29(1):47–65. [12] Asaf Z, Rubinstein D, Shmulevich I. Evaluation of link-track performances using DEM. J. Terramech 2006;43:141–61. [13] Shmulevich I, Asaf Z, Rubinstein D. Interaction between soil and a wide cutting blade using the discrete element method. Soil Tillage Res 2007;97:37–50. [14] Gromov V. Physical and mechanical properties of lunar and planetary soils. Earth Moon Planets 1999;80:51–72. [15] Liu Y, Schnare D, Park JS, et al. Shape analyses of lunar dust particles for astronaut toxicological studies. Lunar Planet Sci XXXVIII; 2007 [1383.pdf]. [16] Carrier DW. Lunar soil grain size distribution. Moon 1973;6:250–63. [17] Cundall PA, Strack ODL. Modeling of microscopic mechanics in granular material. In: Jenkins, JT, Satake M, editors. Mechanics of granular materials: new models and constitutive relations. Amsterdam: Elsevier; 1983. p. 113–49. [18] Liao CL, Chang TP, Young DH, et al. Stress strain relationship for granular materials based on the hypothesis of best fit. Int J Solids Struct 1997;34(31):4087–100. [19] Oida A, Schwanghart H, Ohakubo S. Effect of tire lug cross section on tire performance simulated by distinct element method. In: Proceedings of the 13th international conference of the ISTVS, Munich; 1999. p. 345–52. [20] Li W, Gao F, Jia Y. Tractive performance analysis on radially deployable wheel configuration of lunar rover vehicle by discrete element method. Chinese J Mech Eng 2008;21(4):66–72. [21] Buhler CR, Calle CI, Clements JS, et al. Test method for in situ electrostatic characterization of lunar dust. In: Aerospace Conference, 2007. IEEE; 2007. p. 1–19. [22] Stubbs TJ, Vondrak RR, Farrell WM. A dynamic fountain model for lunar dust. Adv Space Res 2006;37:59–66.