The growing interest in planetary
The discrete element method (DEM) was originally developed in the mining industry to examine rock mechanics problems. Growing in popularity, DEM has been applied to many other fields of science and engineering (Barreto and Leak, 2021). The possibility of easily modeling an environment with different gravity levels has also attracted the attention of scientists working on terramechanical problems in space. For example, an alterable DEM model of lunar soil was presented by Liu et al. (2020). The authors validated the model by comparing triaxial test simulations with those obtained from physical experiments. In Jiang et al. (2017), DEM simulations of soil cutting tests performed under different gravity conditions were examined. DEM modeling of a lunar soil simulant was also proposed in Xi et al. (2021). In Li et al. (2021), analog method experiments, which have a long history in tectonic modeling, were used to calibrate DEM models by comparing the granular structures formed during experiments with those yielded by DEM simulations.
In this paper, we present a study conducted within an ISRU project whose one of the tasks is the analysis by DEM of the mechanical behavior of the lunar soil under conditions of reduced gravity. Our modeling activities commenced with the creation of a DEM analog of a repose angle test built in-house. Comparing the results of the physical model with its digital twin allows us to calibrate and validate the model and use it for prediction. However, the calibration phase is not trivial. DEM models are complex and have numerous parameters to tune. The terramechanical phenomena they describe are also complex, highly nonlinear, and not easily repeatable. Thus, the calibration phase can be preceded and supported by a sensitivity analysis, which gives insight into the model behavior when changing its parameters. In this context, past research includes the assessment of the influence of the size and shape of DEM grain parameters (Elekes and Parteli, 2021), the friction coefficient (Li et al., 2021), the damping coefficient (Jiang et al., 2017; Xi et al., 2021), and the friction and restitution coefficients (Zhu et al., 2022; Wang et al., 2021) on the angle of repose (AR) in a repose angle test model. The research conducted so far is extensive, but not exhaustive in any way. In this article, we present a new sensitivity analysis that investigates the influence of a large set of DEM parameters on both AR and the simulation time. We adopt the Morris method (Saltelli et al., 2004) and the stepwise regression (Myers et al. 2016) for the sensitivity assessment. Both analyses represent the first step of a broader ongoing model calibration activity.
The paper is structured as follows. Section 2 describes the DEM model of the repose angle test. Section 3 illustrates the results of the sensitivity analysis. Section 4 provides a short summary and highlights future research.
A repose angle experiment is a simple test which allows the researcher to assess some fundamental mechanical properties of the soil related to the inter-particle friction. The test consists of dropping soil on a horizontal plate from a funnel located above it. The soil piles up on the plate with a cone-like shape. AR is the angle that the side of the soil pile forms with the plate. The three-dimensional DEM model we prepared simulates a repose angle test. The model is made in ANSYS/Rocky software. Below, we present details of this model.
The geometrical model is composed of three rigid elements: a circular plate, a funnel, and a cap, as presented in Figure 1. These elements have properties compliant with the ISO-4324-1977 standard. The main steps of a DEM repose angle simulation are (refer to Figure 2) as follows. The bulk of soil is generated inside the funnel volume (1). The particles settle down and fill the bottom of the funnel volume, while the cap is closed (2). The funnel cap is removed, and the particles start to flow down and pile up on the plate (3). The simulation ends when a minimum kinetic energy criterion is met (4). The AR is then calculated via post-processing simulation results.
The bulk mechanical properties we adopted in our model relate to the JSC-1A (Johnson Space Center Number One) regolith analog (Liu et al., 2020) and are listed in Table 1.
Particles’ material parameters
Density (bulk density) (kg/m3) | 2900 (1740) |
Young's modulus (MPa) | 45 |
Poisson's ratio (-) | 0.5 |
AR is obtained from DEM results by means of a Python script that uses the OpenCV library developed by Bradski et al. (2000). The script fits a line to the pile contour at multiple cross sections, calculates the slopes of the fitted lines, and averages them. In this process, the upper and lower regions of the contours are cropped. The height of these regions is at least 10 times the particle size, as proposed by Chen et al. (2019). It excludes the influence of factors related to the impact of particles at the top, as well as the friction of the table. Profile fitting is based on the M-estimator technique with the use of a Welsch function as the distance type. An example of application of the method is presented in Figure 3.
Simulation results strongly depend on the contact model that is employed. ANSYS/Rocky includes different contact models for normal and tangential forces and rolling resistant torque. We examined in our analysis three different models of normal forces, two models of tangential forces, and a model of rolling resistance. A synthetic description of these models is provided below.
In this model, the normal forces are calculated in two different ways, depending on the value of the normal contact overlap (NCO). NCO is defined as the difference between the normal overlap values of successive simulation time steps. When the NCO is positive, the normal force is calculated from (1), otherwise from (2).
In (1) and (2),
In this model, normal forces are given by
In (6),
This model takes the form
The parameter
In this model, tangential forces are expressed by
The maximum relative tangential displacement
This model works only with the Hertzian spring dashpot model.
In this model, the tangential force is elastic-frictional. Its expression is
All simulations that we ran used a linear spring rolling limit model. The resistant moment is calculated by
In (22),
The model described in Section 2 entails many parameters. Past research and initial model examination evidenced a subset of parameters that are expected to largely affect the simulation results. This subset includes the friction coefficient (
The shape and size of the particles are parameters that also play a big role in the model behavior (see, e.g., Elekes and Parteli, 2021). The lunar regolith analog is a heterogeneous material with grains of various shapes and sizes. However, the implementation of variable particle shapes and sizes results in models that are cumbersome from a computational point of view. Since the final goal of this research is the creation of a large model that simulates the entire excavation process, the possibility of using only spherical particles is of interest. This research focuses on particles with a spherical shape and uniform size. Further research will examine whether this assumption is admissible or it leads to an unacceptably poor description of the physical behavior of the soil dynamics. Notwithstanding the uniformity of the particles, the impact of their size was analyzed. We initially performed experiments with variable particle size, while keeping all other parameters fixed. In general, it is expected that reducing the size can improve the ability of the model to replicate the physics of the terramechanical phenomena, but, on the downside, it increases the computational time. Thus, the best simulation time/simulation accuracy trade-off is to be chosen.
This preliminary study involved particles with six different diameters ranging from 0.75 to 2 mm. The other contact model parameters considered to be highly influential were set to the values listed in Table 2.
Materials’ interaction parameters used in the grain size and contact model analyses
Static friction (−) | 0.7 | 0.7 |
Dynamic friction (−) | 0.7 | 0.7 |
Tangential stiffness ratio (−) | 1 | 1 |
Adhesive distance (m) | 0.0001 | 0.0001 |
Force fraction (−) | 0.2 | 0 |
Restitution coefficient (−) | 0.2 | 0.3 |
The results of the repose angle test simulations for all particle sizes are presented in Table 3 and Figure 4. They show that AR tends to decrease with particle size. Observing AR versus the particle size diagram of Figure 4, we note that the AR values at particle sizes of 1.75 and 0.75 mm do not follow the trend because they are affected by the “avalanching effect”. This effect occurs when AR of the pile is exceeding because of the addition of a large amount of granular matter to the top of the pile. Then the soil becomes unstable and, losing its shear strength, undergoes a rapid flow or sliding. If the simulation were stopped at the very moment of its occurrence or right after it would yield AR that are artificially greater (in the first case) or smaller (in the second case). This effect, poorly documented in the literature, is under further investigation. The computational time versus particle size diagram shows a typical inversely proportional relationship that is highlighted in the plot with a red dashed line. If we assume that the model accuracy always increases with the particle size reduction, then the constraint for the size is given only by the computational time. Looking at the time versus size plot of Figure 4, it seems that a particle size of 1.5 mm can provide a good trade-off between these two parameters; therefore, it has been chosen in the next simulations.
AR for different DEM contact models
Mindlin–Deresiewicz | 42.58 | - | - |
Linear spring Coulomb limit | 43.00 | 45.03 | 41.40 |
The AR values are expressed in degrees.
The configurations of the contact models described in Section 2 that we examined are as follows:
Hertzian spring dashpot | Mindlin–Deresiewicz; Hertzian spring dashpot | Linear spring Coulomb limit; Linear spring dashpot | Linear spring Coulomb limit; and Hysteretic linear spring | Linear spring Coulomb limit.
Other combinations are not possible because the Mindlin–Deresiewicz tangential model works only with the Hertzian spring dashpot normal model. The simulation results we obtained are given in Tabs 4 and 5. The goal of this analysis is to find an adequate contact model that ensures, at the same time, feasible computational time. The criterion we used for assessing the model adequacy is the comparison of its AR prediction with a documented value for JSC-1A. This value oscillates between 37° and 41° (Calle and Buhler, 2020; Kobaka et al., 2023). The simulations we performed show that the utilization of the linear spring dashpot yielded an AR that is sensibly larger than the reference value. Differently, both the hysteric linear spring and the Hertzian model provide similar results, although the former requires twice the computational time of the latter. Thus, we could conclude that the Hertzian spring dashpot–Linear spring Coulomb limit contact model combination should be used in the future. It should be noted that the criterion we used to assess the model fidelity has evident flaws in the fact that it was employing a model with non-tuned parameters and it was run of a single model configuration. Nevertheless, the validity of our choice of contact models is demonstrated in previous research, as, for instance, in Knuth et al. (2012) and Sanchez and Scheeres (2020).
Computation time for different DEM contact models
Mindlin–Deresiewicz | 117.73 | - | - |
Linear spring Coulomb limit | 122.93 | 181.41 | 237.09 |
The CPU time is in min.
Ranges of variability of the contact model parameters used in the sensitivity analysis
p1 | Friction coefficient ( |
0.2 | 0.5 |
p2 | Adhesive distance ( |
0.0005 | 0.0015 |
p3 | Force fraction ( |
0.1 | 0.4 |
p4 | Restitution coefficient ( |
0.1 | 0.4 |
p5 | Rolling resistance ( |
0.2 | 0.8 |
p6 | Tangential stiffness ratio ( |
0.2 | 1.0 |
We examined the influence of contact model parameters on AR through two numerical tools: the Morris parameter screening method and a linear regression built using a stepwise method. In these analyses, the parameters are assumed to vary within the ranges defined in Table 6.
The Morris method is a sample-based procedure that performs a global sensitivity analysis. This method is reported to yield qualitative results with a reduced number of simulations. In the Morris method, the samples are randomly picked from a regular hyper-grid according to a one-variable-at-the-time scheme. This sampling procedure defines a set of trajectories that try to uniformly cover the whole parameter space. The simulation responses are processed, and two summary sensitivity metrics are calculated for every model parameter (also called factor): the elementary effect mean,
Figure 5 shows that all parameters, except the tangential stiffness ratio, sensibly affect the repose angle, but the rolling resistance and the friction coefficient stand out for their level of importance. However, the rolling resistance appears to be more linearly correlated with AR than the friction coefficient. The simulation time appears to be highly and linearly correlated with the adhesive distance.
Stepwise regression is a smart way to build regression models without prior definition of the shape of the approximation model. The method progressively inserts and removes regressors (from a defined set of candidates) into the approximation until no further improvement is observed in the regression. This sequential procedure usually yields a response surface that well captures the response behavior with a limited number of regressors, therefore eliminating the overfitting problem. In this study, we used the 49 sample points of the Morris analysis to build the regression model. The candidate regressors were linear, quadratic, and two-way interaction terms. The regression model of AR obtained with the stepwise regression takes the form
The regression model of the simulation time
The regression variables in (25) and (26) read as follows:
A coded variable is the variable obtained by rescaling the physical parameter range, to a [0, 1] range. Using coded variables permits to directly read the influence of a regressor on the response by looking at its associated coefficient. Relative comparison of the regressor coefficients shows good agreement with the Morris sensitivity analysis. The accuracy of the regression models can be assessed by the coefficient of multiple determination. The calculated coefficients are 0.86 and 0.99 for the AR and T models, respectively. The coefficient of multiple determination has a maximum value of 1, meaning that the regression models we obtained fit well the training data.
In this article, we presented various aspects of DEM modeling for repose angle test simulations. At the beginning, we examined the influence of the particle size on the value of AR and the calculation time. Next, we performed computer simulations using different contact models to gain knowledge of their behavior in the tested scenario. Finally, we performed a sensitivity analysis to understand how six fundamental contact parameters may affect the simulation results. The analysis pointed out parameters that are influential in the simulation. This insight will be helpful in the DEM calibration phase. The analysis also revealed that repose angle test results may be affected by the avalanching effect, a phenomenon that is rarely documented in DEM calibration literature. This problem is currently under investigation.