#### Abstract

In the design of planetary sampling devices, calculating the reaction force acting on the sampling devices is crucial. According to related research, the influence zone caused by sampling plays an important role in calculating the reaction force. A new method for estimating the range of the influence zone based on 3D DEM simulation is discussed in this paper. Taking lunar soil as an example, first, via validation of physical experiments, the DEM lunar soil simulant was proven to have mechanical properties similar to those of real lunar soil. Second, stress was selected as an indicator to identify the influence zone by computing the match percentage via a comparison between classical soil mechanics and DEM simulation. Using the proposed calculation method, it can be observed that the trend of change of influence zone at different sampling moments showed similar to the change of reaction force. Calculation of the influence zone can be used to analyze the reaction force of different gravity environments, sampling device structures, and motions.

#### 1. Introduction

In the design of planetary sampling devices, the reaction force acting on the sampling devices is indispensable for calculating the minimum driving force of the sampling equipment as a whole because the driving force directly affects the selection of the power support equipment, power consumption, external shape, and internal structure of the sampling rover. Even the mass of the entire device is partially dependent on the driving force because weight can be calculated according to the device’s shape, structure, and material. Hence, calculation of the reaction force acting on the sampling devices plays a highly important role in the design of sampling devices.

According to related research, the influence zone caused by sampling is crucial for the calculation of the reaction force. In classical soil mechanics, such as the Prandtl theory, the reaction forces are always deduced from the assumed influence zone, such as zones I/II/III shown in Figure 1(a) [1]. In physical experiments, the influence zone can also be easily observed (Figure 1(b)).

**(a)**

**(b)**

Although the influence zone can be identified from the hypothesis in classical soil mechanics theory, obvious disadvantages still exist, and the change in the mechanical properties due to the influence is ignored, taking the example of lunar soil, which is easily influenced by sampling. When the influence is small (as shown in Figure 2 when digging angle smaller than 100°), the calculation results from classical theory (red dashed line) agreed with the experimental results (blue solid line). In contrast, when obvious influence exists (as shown in Figure 2 when digging angle greater than 100°), the theoretical results are much larger than the experimental results, and the maximum error can increase such that the measured data are only 20% of the theoretical results [2].

To avoid potential error, research on the influence zone is necessary, also because it is the key reason leading to the change of mechanical properties of lunar soil. Traditionally, only two methods (in situ tests and laboratory tests) can be used to detect the change in mechanical properties during sampling in traditional physical experiments, but the equipment size for in situ tests is too large, and the results are dependent on an empirical formula instead of accurate mathematical deduction. For laboratory tests, the mechanical properties of soil immediately change once it has been removed from its initial position. Therefore, neither of these methods can obtain accurate data. With the development of the discrete element modeling (DEM for short) simulation method [3], which can reach accuracies similar to those of the results from physical tests [4], DEM has become a powerful supplement to physical experiments. In DEM simulation, it is easy to observe the interaction effect among particles in the sampling process, and hence, it can be used to monitor changes in mechanical properties. Few researchers had reported an accurate calculation of the influence zone. Only Li et al. had tried to predict the influence zone to calculate the reaction force during the sampling process. However, the influence zone in this research was predicted from DEM simulation with two dimensions, and the core criterion to identify the influence zone was the displacement for each particle [5].

To achieve this goal, the first step is to determine the influence zone caused by sampling in DEM simulation, and a new method for estimating the range of the influence zone based on 3D DEM simulation and the potential application of the influence zone is discussed in this paper.

#### 2. Generation of the DEM Lunar Soil Simulant

##### 2.1. Random Generation

The first step in generating the lunar soil simulant is the random generation in which spheres are used to fill a fixed domain, and the diameter of these spheres is uniformly distributed within a specified range. In this simulation, the diameter of the radius of particle distributions from 0.25 to 0.5 cm is used. The particle size distribution curves and their control sizes are shown in Figure 3. This step is designed to limit the position and particle size of nonspherical particles instead of conducting a simulation.

It should be noticed that this size distribution does not coincide with the real lunar soil size distribution, because there will be an extreme increase of the computing time and capability for simulation if using real size distribution which is too small. Considering about the target of simulation is to figure out the mechanical performance for lunar soil, if the simulated macroperformance is consistent with the physical experiments (which will be described in Section 2.4), we still believe that this particle distribution was acceptable.

##### 2.2. Particle Replacement

The second step is to replace the spherical particles generated from random generation with nonspherical particles. Considering the computational efficacy, the rigid spherical cluster is selected, and no overlap occurs among the spherical elements to avoid the potential problems of particle mass and density, as shown in Figure 4, 6-element cluster were selected here. The particle shape in this work is defined from the main features of real lunar soil particles, and specifically, an ellipsoidal approximation method is applied in which the long, intermediate, and short axes of a fitted ellipsoid are defined as , , and , respectively [5, 6]. Based on the measurements of real lunar soil particles from the Apollo 16 mission, the ratios of the three axes can be obtained and are listed in Table 1 [6, 7].

##### 2.3. Sedimentation

As shown in Figure 5, after replacement, obvious voids are present among particles, and hence, the final step in establishing the lunar soil simulant is to sediment these particles in the lunar gravity environment (cm/s^{2}). During sedimentation, the contact force among particles gradually reaches equilibrium status and generates the initial earth pressure. To ensure accuracy, the following simulation or calculation must be conducted after equilibrium status is attained to prevent the external forces acting on every particle from changing. Hence, observing whether the forces acting on each particle are stable is an effective way to assess whether the sedimentation is completed. After statistics of the force of all particles during the whole process of sedimentation, it is found that the force trend of each particle is roughly the same, and the difference is only a slight fluctuation in the magnitude of the force, taking the particle of No. 300 as an example, as shown in Figure 5. At the beginning of sedimentation, a large void is present around this particle, and the external forces are small. After a period of downward motion, the particle starts to contact other particles, and the forces begin to change dramatically. As the simulation continues, the forces become steady to complete the sedimentation and reach the equilibrium status.

##### 2.4. Validation

An actual lunar soil simulant FJS-1 developed in Japan was used to validate whether the simplified model proposed in this work can obtain similar mechanical properties. First, as listed in Table 2, the microscopic parameters were used in DEM to generate the 6-element particles. The static mechanical properties were validated as listed in Table 3, and it can be concluded that the properties of the DEM simulant are similar to those of FJS-1, except for particle size (which is limited by the calculation capacity of the workstation) and cohesion (which is ignored in the DEM in this work).

Considering the validation of the dynamic mechanical behavior, direct shear tests and penetration experiments were carried out. The DEM simulations of direct shear tests were conducted with the number of 2000 particles placed in the domain with the size of mm, and the experimental data of direct shear tests of FJS-1 were collected [9]. As shown in Figure 6(a), the comparison showed that the DEM simulation results were in good agreement with the measured result. And the penetration experiment was performed using a flat plat, of which the size of the cross-section is mm, and the diameter of the container is 100 mm (Figure 6(b)). A total of 8000 particles was placed in the domain with a size of mm, as shown in Figure 6(c), and the boundary size was slightly smaller to ensure sufficient penetration depth (because the number of particles is limited by the calculation capacity of the workstation). The validation result (Figures 6(a) and 6(d)) indicates that this DEM simulant consisting of 6-element particles can obtain dynamic mechanical properties similar to those of FJS-1.

**(a)**

**(b)**

**(c)**

**(d)**

#### 3. Calculation of the Influence Zone

##### 3.1. Definition of the Influence Zone

In soil mechanics, stress is an essential component of research on soil behavior, and therefore, stress was selected to estimate the range of influence zone. In the DEM simulation, stress is calculated from the contact forces among particles in a fixed domain, known as a stress monitor, as shown in Equation (1) [3]. where is the total number of contacts in the volume , is the force vector for contact , and is the branch vector for contact . For the nonspherical particles in this simulation, as shown in Figure 7, the branch vector must be considered the vector between the centers of the particles instead of the centers of the elements [3].

**(a)**

**(b)**

According to Equation (1), the stress in the DEM can be considered the average value of the contact force, and hence, the number of contacts, which only depends on the monitor size and particle size, is the key to obtaining accurate results. As shown in Figure 8, a stress monitor can be defined as a range, in which all the contact forces of every particle in this range will be introduced into Equation (1) for stress calculation. And the calculation result can be regarded as the stress value at the centroid position of this range. Several experiments were conducted with increasing size of stress monitors with the same central coordinate, which was the centroid of all the DEM particles after sedimentation. Similar results were achieved, as shown in Figure 9, unless the stress monitor reaches the minimum suitable size (nearly equal to 4 cm), below which the calculated stress is not stable.

To simplify the following calculation, principal stresses are introduced from classical soil mechanics. As shown in Figure 9, the global coordinate system contains 6 types of stress for a random position in soil, including 3 types of normal stress () and 3 types of shear stress (). By rotating the coordinate system, the 6 types of stress can be replaced by 3 types of normal stress, referred to as principal stress (), where the maximum principal stress is selected as representative of the entire stress state at this random position.

With the accurate calculation of the stress value and its monitor size, as shown in Figure 9, the stress of a series of positions around the sampling devices can be calculated. To determine the change in stress result due to the sampling process, the stresses of the same position were also calculated without sampling. According to the two stress calculations, no mutation occurs for the stress change value in the DEM simulation, and the change value of stress was directly related to the calculation position. Therefore, a threshold value should be defined to estimate the range of the influence zone, as shown in Equation (3). If the stress change value of any position is larger than the threshold value, this position should be located in the influence zone.

where is the change value of the maximum principal stress at the same position, is the maximum principal stress value during sampling, is the maximum principal stress value without sampling, and is the threshold value of stress change.

##### 3.2. Calculation of the Threshold Stress Change Value

In the classical Prandtl theory, the range of influence zone does not change if the penetration depth, device width, and length remain the same, as shown by the area filled by blue circles in Figure 10. According to the definition in Equation (3), the value of the threshold stress is the key parameter for identifying the most suitable influence zone in which the different threshold stress change value leads to a different influence zone in DEM simulation (areas indicated by red circles shown in Figure 10).

According to Figure 8, the size of the influence zone is fixed in the Prandtl theory and changes with the magnitude of the threshold stress in the DEM simulation. Hence, the most suitable threshold stress results in the influence zone of the DEM simulation that is closest to that from the theoretical calculation. To determine the most suitable threshold value, a parameter known as the match percentage was defined to describe the difference between the range in the DEM simulation and classical soil mechanics, as shown in Equation (4). where is the match percentage, indicates the volume of the influence zone in the DEM simulation, indicates the volume of influence zone in Prandtl theory, and indicates the volume of the influence zone in both thr DEM simulation and Prandtl theory.

Considering the calculating process, the volume of the influence zone can be easily calculated from thr Prandtl theory [1]. And the volume can be calculated by counting the center position and radius of all selected particles and finding the outermost layer of the particle group. It should be noticed that the particles selected here indicated that the particles meet the requirement of Equation (3). With the knowledge of volume and volume , it is easy to compute volume .

Via Equation (4), when the DEM influence zone perfectly matches the theoretical range, where , the match percentage . A match percentage closer to 1 indicates that the size of the DEM influence zone is closer to the size of the theoretical influence zone. Therefore, the most suitable value for the threshold stress can be obtained from the maximum match percentage .

We take an example from the DEM simulation in which the penetration depth is 35 mm, the device width is 5 mm, and the length is 30 mm. As shown in Figure 11(a), the match percentage changes with the change in the threshold stress . When the threshold stress is quite small, the DEM influence zone is much larger than the theoretical range (), which leads to a rather small match percentage. With the increase in the threshold value, the DEM range decreases dramatically, leading to a significant increase in the match percentage. When the threshold stress reaches the most suitable value (approximately 150 Pa), the peak value of the match percentage is also obtained. Subsequently, with the continual growth of the threshold stress , the DEM range continues to decrease and becomes smaller than the theoretical range (), and the match percentage gradually decreases and approaches zero.

**(a)**

**(b)**

#### 4. Potential Application of the Influence Zone

##### 4.1. Analysis of Penetration in a Specific Moment

The calculation of the influence zone can help understand the stress distribution for any specific moment in the sampling process. As shown in Figure 12, for a penetration device such as a thin plate (the entire influence zone can be observed in Figure 11(b)), the influence zone distribution is nonuniform around the penetration devices. First, for a penetration motion, the stresses in the influence zone are the total compressive stress because the magnitude of the stress is below zero. Second, it can be concluded that the closer the location to the penetration device, the greater the change in stress, and the maximum stress appears at the position directly below the penetration device. Third, the area near the short side of the penetration device is more susceptible to influence, as suggested by the observation that the cross-section along the length of the influence zone is larger than that along the width.

##### 4.2. Analysis of the Penetration Process

With the proposed calculation method, the influence zone in the entire sampling process can also be observed. As shown in Figure 13, in the simulation of penetration, with increasing penetration depth, an increasing trend of the influence zone range is noted, and the magnitudes of the changed stress are also increased. This changing trend can be viewed as the direct reason for the increasing reaction force during the penetration process and can also be used to support the future calculation of the changing trend for the mechanical properties of soil.

##### 4.3. Other Potential Applications

In addition to the application mentioned above, the influence zone calculation method proposed in this work can also be used in potential research as follows: (i)First and most importantly, we must determine a method for calculating the trend of changing mechanical properties of soil in the area of the influence zone because the influence zone is calculated from the stress, which is closely related to the mechanical properties of soil, such as bulk density, compressive strength, and shear strength. Therefore, it is easy to identify the changing pattern of these properties from the known changed stress(ii)Second, we can compare the influence zone generated from various structures of sampling devices. This paper primarily focuses on the simple structure of thin-plate penetration devices because a similar classical soil mechanics method can be used to calculate and validate the accuracy of the influence zone. Indeed, the structures of sampling devices used around the world are varied, and thus, it is important and meaningful to study the difference in the influence zones for different sampling structures(iii)Third, we can compare the influence zone generated from various sampling motions. For planetary sampling, multiple types of sampling motions are available, especially for surface sampling, which is free from space limitations compared with deep sampling. The different sampling motion generates a different influence zone, which leads to a different reaction force and sampling efficiency(iv)Fourth, we can compare the influence zone generated from different gravity environments. The different gravity environment is another important factor that affects different soil mechanical behavior and reaction force calculation. More importantly, the determination of this difference is expected to help researchers to fix the experimental data that are conducted on earth to predict the result with higher accuracy under different gravity environments(v)Finally, we can determine the sampling reaction force calculated from the known influence zone. After completing the work mentioned above, the next crucial goal is to deduce the reaction force during the sampling process, which is also the main target for studying the influence zone. With the accurate calculation of the reaction force, the cost of experiments for a new sampling structure can be reduced, and this information can also support the study of the internal relationship between the sampling devices and samples

#### 5. Conclusion

Reaction force calculation plays an important role in the design of planetary sampling devices. Because soil (especially soil from extraterrestrial planets, such as lunar soil) can be easily influenced during the sampling process, this influence leads to a change in the mechanical properties of soil, and hence, the influence zone is the key to calculating the reaction force during sampling. To obtain an accurate influence zone, accurate DEM simulation is the primary goal. First, DEM lunar soil was generated after three basic steps of random generation, particle replacement, and sedimentation. Second, according to static and dynamic validation, DEM lunar soil was proven to have mechanical properties similar to real lunar soil. Stress, which is crucial for soil mechanical behavior, was selected as the criterion to identify the influence zone. By calculating the match percentage in comparison with the theoretical influence zone of the classical Prandtl theory, the threshold stress used to determine the DEM influence zone was determined. The application of the influence zone is wide, e.g., the stress at a specific moment or during the entire sampling process can be observed. Furthermore, with an accurate influence zone, the change in the mechanical properties of soil in the influence zone and the accurate reaction force from various sampling structures, sampling motions, and different gravity environments can also be calculated in future work.

#### Data Availability

All the data needed to evaluate the conclusions are present in the paper. Additional data related to this paper may be requested from the corresponding author.

#### Conflicts of Interest

No conflict of interest exists in the submission of this paper.

#### Authors’ Contributions

The manuscript has been approved by all authors for publication. All authors listed have approved the manuscript that is enclosed.

#### Acknowledgments

The authors gratefully acknowledge Professor Takashi Matsushima and Miss Midori Morikawa of the University of Tsukuba for access to DEM software, equipment, and experiments during this research. This project was financially supported by the National Natural Science Foundation of China (No. 42072344 and No. 11502034) and the China Scholarship Council (No. 201608510014). Chengdu University of Technology, the State Key Laboratory of Geohazard Prevention and Geoenvironment Protection, and the Institute of Exploration Technology (CAGS) are also acknowledged for their theoretical support.