Sensitivity Analysis of the DEM Model Numerical Parameters on the Value of the Angle of Repose of Lunar Regolith Analogs

The growing interest in planetary in situ resources utilization (ISRU) activities opens new frontiers to the design of space mining equipment. The development of space technologies is generally very costly and time-consuming because standard design approaches used in terrestrial applications are not fully applicable to extraterrestrial scenarios. It has been observed that while pressure, radiation, and temperature conditions existing in space are relatively simple to reproduce on Earth, low-gravity conditions are very difficult to recreate (Wilkinson and DeGennaro, 2007; Just et al., 2020). On the other hand, the reduced gravity of celestial bodies such as Mars and Moon would surely affect the in situ excavation process due to the peculiar interaction between soil particles and machine elements. It is one of the reasons why the design of tools for space applications is very challenging.

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.

REPOSE ANGLE TEST DEM MODEL

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.

2.1.

Model setup

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.

Geometrical setup for the DEM repose angle test in three views

Subsequent steps of the repose angle test

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.

Table 1.

Particles’ material parameters

Parameter	Value
Density (bulk density) (kg/m³)	2900 (1740)
Young's modulus (MPa)	45
Poisson's ratio (-)	0.5

2.2.

Calculation of AR

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.

Example of the 3D DEM obtained pile (A), averaged particles pile (B), and the result of the angle of repose evaluation (C)

2.3.

DEM contact force models

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.

2.3.1.

Normal forces

Hysteretic linear spring model

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). (1) $F_{n}^{t} = min (K_{nl} \cdot s_{n}^{t}, F_{n}^{t - Δ t} + K_{nu} \cdot Δ s_{n})$ F_n^t = {\rm{\min}}\left( {{K_{nl}} \cdot s_n^t ,F_n^{t - \Delta t} + {K_{nu}} \cdot \Delta {s_n}} \right) (2) $F_{n}^{t} = max (λ \cdot K_{nl} \cdot s_{n}^{t}, F_{n}^{t - Δ t} + K_{nu} \cdot Δ s_{n})$ F_n^t = {\rm{\max}}\left( {\lambda \cdot {K_{nl}} \cdot s_n^t ,F_n^{t - \Delta t} + {K_{nu}} \cdot \Delta {s_n}} \right)

In (1) and (2), $F_{n}^{t}$ F_n^t denotes the normal force at time t, Δt the simulation time step, $s_{n}^{t}$ s_n^t the NCO, Δs_n the change of NCO, and λ is (usually) a small empirical constant (in this study, its value is set to 0.001). K_nl and K_nu are the loading and unloading contact stiffnesses, respectively, defined as: (3) $K_{nl} = E \cdot L$ {K_{nl}} = E \cdot L (4) $K_{nu} = \frac{K_{nl}}{ε^{2}}$ {K_{nu}} = {{{K_{nl}}} \over {{\varepsilon ^2}}} where ɛ is the coefficient of restitution, L the particle size, and E the Young’s modulus.

Linear spring dashpot model

In this model, normal forces are given by (5) $F_{n} = K_{nl} \cdot s_{n} + C_{n} \cdot {\dot{s}}_{n}$ {F_n} = {K_{nl}} \cdot {s_n} + {C_n} \cdot {\dot s_n} with C_n denoting the damping coefficient calculated by (6) $C_{n} = 2 \cdot η \cdot \sqrt{m^{*} \cdot K_{nl}}$ {C_n} = 2 \cdot \eta \cdot \sqrt {{m^*} \cdot {K_{nl}}}

In (6), m^* is the effective mass of the contact and η is the damping ratio obtained by solving (7)–(9) by means of the numerical procedure described in Schwager and Pöschel (2007). (7) $ε = e^{[- \frac{η}{\sqrt{1 - η^{2}}} \cdot (π - arctan \frac{2 \cdot η \cdot \sqrt{1 - η^{2}}}{1 - 2 \cdot η^{2}})]} for 0 \leq η \leq \frac{1}{\sqrt{2}}$ \varepsilon = {e^{\left[ { - {\eta \over {\sqrt {1 - {\eta ^2}} }}\cdot\left( {\pi - {\rm{\;arctan\;}}{{2 \cdot \eta \cdot \sqrt {1 - {\eta ^2}} } \over {1 - 2 \cdot {\eta ^2}}}} \right)} \right]}}\;\;\;{\rm{for}}\;0 \le \eta \le {1 \over {\sqrt 2 }} (8) $ε = e^{[- \frac{η}{\sqrt{1 - η^{2}}} \cdot (arctan \frac{2 \cdot η \cdot \sqrt{1 - η^{2}}}{2 \cdot η^{2} - 1})]} for \frac{1}{\sqrt{2}} \leq η \leq 1$ \varepsilon = {e^{\left[ { - {\eta \over {\sqrt {1 - {\eta ^2}} }}\cdot\left( {{\rm{\;arctan\;}}{{2 \cdot \eta \cdot \sqrt {1 - {\eta ^2}} } \over {2 \cdot {\eta ^2} - 1}}} \right)} \right]}}\;\;\;{\rm{for}}\;{1 \over {\sqrt 2 }} \le \eta \le 1 (9) $ε = e^{[- \frac{η}{\sqrt{η^{2} - 1}} \cdot (\ln \frac{η + \sqrt{η^{2} - 1}}{η - \sqrt{η^{2} - 1}})]} for η \geq 1$ \varepsilon = {e^{\left[ { - {\eta \over {\sqrt {{\eta ^2} - 1} }}\cdot\left( {ln{{\eta + \sqrt {{\eta ^2} - 1} } \over {\eta - \sqrt {{\eta ^2} - 1} }}} \right)} \right]}}\;\;\;\;\;{\rm{for}}\;\eta \ge 1

Hertzian spring dashpot model

This model takes the form (10) $F_{n} = K_{H} \cdot s_{n}^{\frac{3}{2}} + C_{H} \cdot s_{n}^{\frac{1}{4}} \cdot {\dot{s}}_{n}$ {F_n} = {K_H} \cdot s_n^{{3 \over 2}} + {C_H} \cdot s_n^{{1 \over 4}} \cdot {\dot s_n} where K_H is the Hertzian contact stiffness (11) $K_{H} = \frac{4}{3} \cdot E^{*} \cdot \sqrt{R^{*}}$ {K_H} = {4 \over 3} \cdot {E^*} \cdot \sqrt {{R^*}} and C_H is the Hertzian damping coefficient given by (12) $C_{H} = 2 \cdot η_{H} \cdot \sqrt{m^{*} \cdot K_{H}}$ {C_H} = 2 \cdot {\eta _H} \cdot \sqrt {{m^*} \cdot {K_H}}

The parameter R^* denotes the effective equivalent radius, while E^* and η_H are additional coefficients calculated from (13) $\frac{1}{E^{*}} = \frac{1 - ν_{1}^{2}}{E_{1}} + \frac{1 - ν_{2}^{2}}{E_{2}}$ {1 \over {{E^*}}} = {{1 - \nu _1^2} \over {{E_1}}} + {{1 - \nu _2^2} \over {{E_2}}} (14) $η_{H} = \frac{\sqrt{5}}{2} \cdot η$ {\eta _H} = {{\sqrt 5 } \over 2} \cdot \eta where ν is the Poisson’s ratio.

2.3.2.

Tangential forces

Mindlin–Deresiewicz model

In this model, tangential forces are expressed by (15) $F_{T} = - f \cdot F_{n} \cdot (1 - {(1 - \frac{min (|s_{T}|, s_{T, max})}{s_{T, max}})}^{\frac{3}{2}}) \cdot \frac{s_{T}}{|s_{T}|} + η_{T} \cdot \sqrt{\frac{6 \cdot f \cdot m^{*} \cdot F_{n}}{s_{T, max}}} \cdot {(1 - \frac{min (|s_{T}|, s_{T, max})}{s_{T, max}})}^{\frac{1}{4}} \cdot {\dot{s}}_{T}$ {F_T} = - f \cdot {F_n} \cdot \left( {1 - {{\left( {1 - {{{\rm{\min}}\left( {\left| {{s_T}} \right|,{s_{T,{\rm{max\;}}}}} \right)} \over {{s_{T,{\rm{max\;}}}}}}} \right)}^{{3 \over 2}}}} \right) \cdot {{{s_T}} \over {\left| {{s_T}} \right|}} + {\eta _T} \cdot \sqrt {{{6 \cdot f \cdot {m^*} \cdot {F_n}} \over {{s_{T,{\rm{max\;}}}}}}} \cdot {\left( {1 - {{{\rm{\min}}\left( {\left| {{s_T}} \right|,{s_{T,{\rm{max\;}}}}} \right)} \over {{s_{T,{\rm{max\;}}}}}}} \right)^{{1 \over 4}}} \cdot {\dot s_T} where f is the friction coefficient, s_T is the tangential displacement, and η_T is the tangential damping ratio calculated as (16) $η_{T} = \frac{ln ε}{\sqrt{{ln}^{2} ε + π^{2}}}$ {\eta _T} = {{{\rm{\;ln\;}}\varepsilon } \over {\sqrt {{\rm{l}}{{\rm{n}}^2}\varepsilon + {\pi ^2}} }}

The maximum relative tangential displacement s_T_,max is defined by (17) $s_{T, max} = f \cdot {(\frac{1 - ν_{1}}{2 - ν_{1}} + \frac{1 - ν_{2}}{2 - ν_{2}})}^{- 1} \cdot s_{n}$ {s_{T,{\rm{max\;}}}} = f \cdot {\left( {{{1 - {\nu _1}} \over {2 - {\nu _1}}} + {{1 - {\nu _2}} \over {2 - {\nu _2}}}} \right)^{ - 1}} \cdot {s_n}

This model works only with the Hertzian spring dashpot model.

Linear spring Coulomb limit model

In this model, the tangential force is elastic-frictional. Its expression is (18) $F_{T}^{t} = min (|F_{T, e}^{t}|, f F_{n}^{t}) \cdot \frac{F_{T, e}^{t}}{|F_{T, e}^{t}|}$ F_T^t = {\rm{\min}}\left( {\left| {F_{T,e}^t} \right|,f\;F_n^t} \right) \cdot {{F_{T,e}^t} \over {\left| {F_{T,e}^t} \right|}} where the elastic tangential force is described by (19) $F_{T, e}^{t} = F_{T}^{t - Δ t} K_{T} \cdot Δ s_{T}$ F_{T,e}^t = F_T^{t - \Delta t}{K_T} \cdot \Delta {s_T} with K_T denoting the tangential stiffness, which, in turn, depends on K_nl and the tangential stiffnes ratio r_K according to the equation (20) $K_{T} = r_{K} \cdot K_{nl}$ {K_T} = {r_K} \cdot {K_{nl}}

2.3.3.

Rolling resistance model

All simulations that we ran used a linear spring rolling limit model. The resistant moment is calculated by (21) $M_{r}^{t} = min (|M_{r, e}^{t}|, M_{r, lim}) \cdot \frac{M_{t, e}^{t}}{|M_{r, e}^{t}|}$ M_r^t = {\rm{\min}}\left( {\left| {M_{r,e}^r} \right|,{M_{r,\lim }}} \right) \cdot {{M_{r,e}^t} \over {\left| {M_{r,e}^t} \right|}} where $M_{t, e}^{t}$ M_{r,e}^t denotes the pure elastic rolling resistance moment (22) $M_{r, e}^{t} = M_{r}^{t - Δ T} - K_{r} \cdot ω_{rel} \cdot Δ t$ M_{r,e}^t = M_r^{t - \Delta T} - {K_r} \cdot {\omega _{{\rm{rel}}}} \cdot \Delta t and M_r_,lim is a factor that limits the maximum value of the rolling resistance moment by (23) $M_{r, lim} = f_{r} \cdot R_{r} \cdot F_{n}$ {M_{r,\lim }} = {f_r} \cdot {R_r} \cdot {F_n}

In (22), ω_rel is the relative angular velocity of the grain, R_r the particle rolling radius, f_r the rolling resistance coefficient, and K_r the rolling stiffness given by (24) $K_{r} = R_{r}^{2} \cdot K_{T}$ {K_r} = R_r^2 \cdot {K_T}

SENSITIVITY ANALYSIS

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 (f), the tangential stiffness ratio (r_K), the restitution coefficient (ɛ), and the rolling resistance coefficient (f_r). The friction coefficient and the tangential stiffness are the main parameters for tangential forces. The restitution coefficient is responsible for the nature of the collision. When its value is 1, the contact is perfectly elastic; when it is 0, the contact is perfectly rigid. The rolling resistance affects the particle moment and permits to partially include in the model the torque caused by nonspherical grains while using spherical particles. Furthermore, the addition of adhesive forces to the model (the model of the adhesive force is not described in Section 2) brings two other relevant parameters. They are the adhesive distance (δ_adh) and the force fraction (f_adh). The adhesive distance defines the distance of the application of adhesive forces, whereas the force fraction is a dimensionless parameter that determines the fraction of the adhesive force that will be used in the calculations. The literature review shows that including adhesive forces in the model is particularly important for low-gravity conditions (Jiang et al., 2017).

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.

3.1.

Analysis of the grain size

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.

Table 2.

Materials’ interaction parameters used in the grain size and contact model analyses

Parameter	Particle–particle	Particle–boundary
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.

Angle of repose versus particle size (left) and computational time versus particle size (right) diagrams from the simulation results

Morris plot of the repose angle. The diagram shows that the rolling resistance is the most important factor with a significant linear correlation with the repose angle.

Table 3.

AR for different DEM contact models

Tangential forces models	Normal forces models
Tangential forces models	Hertzian spring dashpot	Linear spring dashpot	Hysteretic linear spring
Mindlin–Deresiewicz	42.58	-	-
Linear spring Coulomb limit	43.00	45.03	41.40

The AR values are expressed in degrees.

3.2.

Influence of contact models

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).

Table 4.

Computation time for different DEM contact models

Tangential forces models	Normal forces models
Tangential forces models	Hertzian spring dashpot	Linear spring dashpot	Hysteretic linear spring
Mindlin–Deresiewicz	117.73	-	-
Linear spring Coulomb limit	122.93	181.41	237.09

The CPU time is in min.

Table 5.

Ranges of variability of the contact model parameters used in the sensitivity analysis

	Parameter	Lower limit	Upper limit
p1	Friction coefficient (f)	0.2	0.5
p2	Adhesive distance (δ_adh)	0.0005	0.0015
p3	Force fraction (f_adh)	0.1	0.4
p4	Restitution coefficient (ɛ)	0.1	0.4
p5	Rolling resistance (f_r)	0.2	0.8
p6	Tangential stiffness ratio (r_K)	0.2	1.0

3.3.

Influence of contact model parameters

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.

3.4.

Morris method

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, μ, and the elementary effect standard deviation, σ. Large μ values indicate a global influence of the parameter. The larger the σ value, the larger the nonlinearity of the parameter or its interaction with another parameter(s). In this work, a Morris setup with six-level hyper-grid and seven trajectories was employed. It resulted in 49 different sample points. The parameter configurations and simulation results for all sample points are given in Table 7 of the Appendix. The two resultant Morris metrics for the repose angle and the simulation time, normalized to their maximum values, are presented in the plots of Figs 5 and 6, respectively.

Morris plot of the simulation time. The diagram shows that the adhesive distance is the predominant factor on the simulation time.

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.

3.5.

Stepwise regression

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 (25) $AR = 30.1 + 22.1 x_{1} + 7.2 x_{2} + 5.9 x_{3} + 5.4 x_{5} + 13 x_{1} x_{5} - 17.3 x_{1}^{2}$ {\rm{AR}} = 30.1 + 22.1{x_1} + 7.2{x_2} + 5.9{x_3} + 5.4{x_5} + 13{x_1}{x_5} - 17.3x_1^2

The regression model of the simulation time T is (26) $T = 6400 - 973 x_{1} + 3163 x_{2} - 637 x_{1} x_{2} + 591 x_{1}^{2} + 3092 x_{2}^{2}$ T = 6400 - 973{x_1} + 3163{x_2} - 637{x_1}{x_2} + 591x_1^2 + 3092x_2^2

The regression variables in (25) and (26) read as follows:

x₁: coded variable representing the friction coefficient parameter

x₂: coded variable representing the adhesive distance

x₃: coded variable representing the friction force

x₅: coded variable representing the rolling resistance

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.

CONCLUSIONS

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.

eISSN:: 2083-6104
Język:: Angielski

Częstotliwość wydawania:: 4 razy w roku
Dziedziny czasopisma:: Geosciences, other

Kanał RSS czasopisma

Sensitivity Analysis of the DEM Model Numerical Parameters on the Value of the Angle of Repose of Lunar Regolith Analogs

Data publikacji: 29 gru 2023

Zakres stron: 188 - 202

Otrzymano: 16 gru 2022

Przyjęty: 01 cze 2023

DOI: https://doi.org/10.2478/arsa-2023-0022

Słowa kluczowe
DEM, regolith, space mining, repose angle test, sensitivity analysis

© 2023 Przemysław Młynarczyk et al., published by Sciendo

This work is licensed under the Creative Commons Attribution 4.0 International License.

Figure 1.

Figure 2.

Figure 3.

Figure 4.

Figure 5.

Figure 6.

Sensitivity Analysis of the DEM Model Numerical Parameters on the Value of the Angle of Repose of Lunar Regolith Analogs

Data publikacji: 29 gru 2023

Zakres stron: 188 - 202

Otrzymano: 16 gru 2022

Przyjęty: 01 cze 2023

DOI: https://doi.org/10.2478/arsa-2023-0022

Słowa kluczoweDEM, regolith, space mining, repose angle test, sensitivity analysis

© 2023 Przemysław Młynarczyk et al., published by Sciendo

This work is licensed under the Creative Commons Attribution 4.0 International License.

Figure 1.

Figure 2.

Figure 3.

Figure 4.

Figure 5.

Figure 6.

Słowa kluczowe
DEM, regolith, space mining, repose angle test, sensitivity analysis