Numerical Simulations of Regolith Sampling Processes

TL;DR

This paper improves numerical simulations of regolith sampling in microgravity using SPH, analyzing how material properties and environmental factors affect sampling efficiency for asteroid missions.

Contribution

It introduces enhanced SPH-based modeling of regolith sampling, incorporating elastic-plastic soil behavior and varying environmental conditions.

Findings

Simulation results align well with Earth experiments.

Material cohesion and friction influence sampling rates.

Surface gravity and brush rotation affect collection efficiency.

Abstract

We present recent improvements in the simulation of regolith sampling processes in microgravity using the numerical particle method smooth particle hydrodynamics (SPH). We use an elastic-plastic soil constitutive model for large deformation and failure flows for dynamical behaviour of regolith. In the context of projected small body (asteroid or small moons) sample return missions, we investigate the efficiency and feasibility of a particular material sampling method: Brushes sweep material from the asteroid's surface into a collecting tray. We analyze the influence of different material parameters of regolith such as cohesion and angle of internal friction on the sampling rate. Furthermore, we study the sampling process in two environments by varying the surface gravity (Earth's and Phobos') and we apply different rotation rates for the brushes. We find good agreement of our sampling simulations on Earth with experiments and provide estimations for the influence of the material properties on the collecting rate.

Figures

Tables

Table 1. Table 1: Material and setup parameters for the sampling simulations. The rotation speed of the motor is positive in the collecting phase and negative for the clearing phase by definition.

Gravity $g$ $\left[m/s^2\right]$	Motor [RPM]	Angle of friction [deg]	Cohesion [kPa]	Run Identifier
9.81 (Earth)	100	42	0	E0
9.81	100	42	0.5	E1
9.81	300	42	0	E2
9.81	500	42	0	E3
$5.7\text{×}{10}^{-3}$ (Phobos)	100	42	0	P0
$5.7\text{×}{10}^{-3}$	100	42	0.5	P1
$5.7\text{×}{10}^{-3}$	100	20	1	P2
$5.7\text{×}{10}^{-3}$	100	42	1	P3
$5.7\text{×}{10}^{-3}$	100	60	1	P4
$5.7\text{×}{10}^{-3}$	500	42	0	P5
$5.7\text{×}{10}^{-3}$	500	42	1	P6
$5.7\text{×}{10}^{-3}$	-100	42	0	P7
$5.7\text{×}{10}^{-3}$	-100	42	0.1	P8

Equations

$${\dot{\epsilon }}^{\alpha \beta }=\frac{1}{2}\left(\frac{\partial v^{\alpha }}{\partial x^{\beta }}+\frac{\partial v^{\beta }}{\partial x^{\alpha }}\right),$$

$${\dot{\epsilon }}^{\alpha \beta }={\dot{\epsilon }}_{\mathrm{e}}^{\alpha \beta }+{\dot{\epsilon }}_{\mathrm{p}}^{\alpha \beta }.$$

$${\dot{\epsilon }}_{\mathrm{e}}^{\alpha \beta }=\frac{1}{2\mu }{\dot{S}}^{\alpha \beta }+\frac{1-2\nu }{3E}{\dot{\sigma }}^{\gamma \gamma }{\delta }^{\alpha \beta }.$$

$${\dot{\epsilon }}_{\mathrm{p}}^{\alpha \beta }=\dot{\lambda }\frac{\partial g}{\partial {\sigma }^{\alpha \beta }}.$$

$${\dot{\epsilon }}^{\alpha \beta }=\frac{1}{2\mu }{\dot{S}}^{\alpha \beta }+\frac{1-2\nu }{3E}{\dot{\sigma }}^{\gamma \gamma }+\dot{\lambda }\frac{\partial g}{\partial {\sigma }^{\alpha \beta }},$$

$${\dot{\sigma }}^{\alpha \beta }=$$

$$-\lambda \left(\left(K-\frac{2\mu }{3}\right)\frac{\partial g}{\partial {\sigma }^{\phi \vartheta }}{\delta }^{\phi \vartheta }{\delta }^{\alpha \beta }+2\mu \frac{\partial g}{\partial {\sigma }^{\alpha \beta }}\right).$$

$$K=\frac{E}{3-6\nu }.$$

$$\dot{\lambda }=\frac{2\mu {\dot{\epsilon }}^{\alpha \beta }\frac{\partial f}{\partial {\sigma }^{\alpha \beta }}+\left(K-\frac{2\mu }{3}\right){\dot{\epsilon }}^{\gamma \gamma }\frac{\partial f}{\partial {\sigma }^{\alpha \beta }}{\delta }^{\alpha \beta }}{2\mu \frac{\partial f}{\partial {\sigma }^{\phi \vartheta }}\frac{\partial g}{\partial {\sigma }^{\phi \vartheta }}+\left(K-\frac{2\mu }{3}\right)\frac{\partial f}{\partial {\sigma }^{\phi \vartheta }}{\delta }^{\phi \vartheta }\frac{\partial g}{\partial {\sigma }^{\phi \vartheta }}{\delta }^{\phi \vartheta }}.$$

$$\frac{\mathrm{d}{v}^{\alpha }}{\mathrm{d}t}=\frac{1}{\varrho }\frac{\partial {\sigma }^{\alpha \beta }}{\partial x^{\beta }},$$

$$\frac{\mathrm{d}\varrho }{\mathrm{d}t}=-\varrho \frac{\partial v^{\alpha }}{\partial x^{\alpha }}.$$

$$f(I_1,J_2)=\sqrt{J_2}+{\alpha }_{\varphi }{I}_1-k_c=0.$$

$$I_1=\mathrm{tr}(\sigma )={\sigma }^{\gamma \gamma }\text{mbox}and{J}_2=\frac{1}{2}{S}^{\alpha \beta }{S}_{\alpha \beta }.$$

$${\alpha }_{\varphi }=\frac{\tan \varphi }{\sqrt{9+12{\tan }^2\varphi }}\text{mbox}and{k}_c=\frac{3c}{\sqrt{9+12{\tan }^2\varphi }}.$$

$$g=\sqrt{J_2}+3I_1\sin \psi ,$$

$$\frac{\partial g}{\partial {\sigma }^{\alpha \beta }}$$

$$=\frac{\partial g}{\partial I_1}{\delta }^{\alpha \beta }+\frac{1}{2\sqrt{J_2}}\frac{\partial g}{\partial \sqrt{J_2}}{S}^{\alpha \beta }.$$

$$\frac{\partial g}{\partial {\sigma }^{\alpha \beta }}=3\sin \psi {\delta }^{\alpha \beta }+\frac{1}{2\sqrt{J_2}}{S}^{\alpha \beta }$$

$$\frac{\partial f}{\partial {\sigma }^{\alpha \beta }}={\alpha }_{\varphi }{\delta }^{\alpha \beta }+\frac{1}{2\sqrt{J_2}}{S}^{\alpha \beta }.$$

$${\dot{\sigma }}^{\alpha \beta }=$$

$$-\dot{\lambda }(\left(K-\frac{2\mu }{3}\right)9\sin \psi {\delta }^{\alpha \beta }$$

$$+2\mu \left(3\sin \psi {\delta }^{\alpha \beta }+\frac{1}{2\sqrt{J_2}}{S}^{\alpha \beta }\right))$$

$$=$$

$$-\dot{\lambda }\left(9K\sin \psi {\delta }^{\alpha \beta }+\frac{\mu }{\sqrt{J_2}}{S}^{\alpha \beta }\right).$$

$$\dot{\lambda }$$

$$g=\sqrt{J_2}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\sigma =\dot{\sigma }-\sigma \cdot R^T-R\cdot \sigma ,$$

$$R^{\alpha \beta }=\frac{1}{2}\left(\frac{\partial v^{\alpha }}{\partial x^{\beta }}-\frac{\partial v^{\beta }}{\partial x^{\alpha }}\right).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{x}^{\alpha }=v^{\alpha },$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{v}^{\alpha }=\frac{1}{\varrho }\frac{\partial {\sigma }^{\alpha \beta }}{\partial x^{\beta }}+k_{\mathrm{ext}}^{\alpha },$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\varrho =-\varrho \frac{\partial v^{\alpha }}{\partial x^{\alpha }},$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\sigma }^{\alpha \beta }={\sigma }^{\alpha \gamma }{R}^{\beta \gamma }+{\sigma }^{\gamma \beta }{R}^{\alpha \gamma }+2\mu {\dot{\epsilon }}^{\alpha \beta }$$

$$+\left(K-\frac{2\mu }{3}\right)\mathrm{tr}(\dot{\epsilon }){\delta }^{\alpha \beta }-\dot{\lambda }\frac{\mu }{\sqrt{J_2}}{S}^{\alpha \beta },$$

$$\displaystyle\dot{\lambda}=\left\{\begin{array}[]{cl}\frac{3K}{\mu}\alpha_{\phi}\mathrm{tr}(\dot{\varepsilon})+\frac{1}{\sqrt{J_{2}}}S^{\alpha\beta}\dot{\varepsilon}_{\alpha\beta}&\mbox{if }\sqrt{J_{2}}+\alpha_{\phi}I_{1}-k_{c}<0,\\ 0&\mbox{else.}\end{array}\right.$$

$$p(\varrho )=\frac{K_0}{n}{\left(\frac{\varrho }{{\varrho }_0}-1\right)}^n.$$

$$\frac{\mathrm{d}{\varrho }_a}{\mathrm{d}t}={\varrho }_a\sum _b\frac{m_b}{{\varrho }_b}(v_a^{\alpha }-v_b^{\alpha })\frac{\partial W_{ab}}{\partial x_{\alpha }},$$

Full text

††journal: Planetary and Space Science

1 Introduction

Regolith is a layer of loose, heterogeneous material. It includes dust, soil, brittle broken rock and can be found on terrestrial planets, moons and asteroids. ESA conducted a feasibility study for a European Phobos sample return mission, which was called Phootprint (Koschny et al., 2014; Barraclough et al., 2014), followed by a Phase A study called Phobos Sample Return (PhSR), which investigated a joint ESA/Roscosmos mission as well as an ESA standalone scenario. Phobos was chosen as a scientifically interesting destination due to its unknown formation process and as a unique opportunity to test some of the key components of a potential follow-up Mars sample return mission. Phobos has a semi-major axis of $9378\text{\,}\mathrm{k}\mathrm{m}$ and a sidereal orbital period of $0.3\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$ in bound rotation. Phobos is irregularly shaped and has dimensions of $27$ $\times$ $22$ $\times$ $18$ $\mathrm{km}$. Phobos’ internal structure is poorly constrained, Andert et al. (2010) provide values for the mean density of about $1876(20)\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$ and a porosity of 30 % $\pm$ 5 %. The surface of Phobos is covered with regolith and the thickness of the layer varies between $5\text{\,}\mathrm{m}$ and $100\text{\,}\mathrm{m}$ (Basilevsky et al., 2014). Various spacecraft missions have investigated the moons of Mars, see the publication by Duxbury et al. (2014) for a complete review of Russian and American missions. However, no spacecraft has landed on Phobos yet.

The PhSR study is a technical baseline for a sample return mission to identify the key technological requirements and address their development requirements for a sample return mission. One of the major milestones is the design of a innovative sampling tool installed on the lander that will allow to collect at least $100\text{\,}\mathrm{g}$ of Phobos’ surface material using a novel rotary brush mechanism. In addition to Phobos’ material, a collected sample will presumably contain pieces of material originated from Mars.

Computer simulations and computer aided testing become more and more important in the field of spacecraft design. By the help of modern computers, we can investigate processes under physical conditions that are not easily available for experiments on Earth, e.g. low-gravity fields. Recently, new models for the dynamic behaviour of regolith have been developed using different numerical schemes: Schwartz et al. (2014) use a soft-sphere discrete element method (SSDEM) in the N-body gravity tree code pkdgrav to investigate low-speed impact simulations into regolith in support of asteroid sampling mechanisms in the context of the JAXA’s Hayabusa and Hayabusa2 missions, and Bui et al. (2008) developed a model for granular flow using SPH in the context of solid mechanics combined with a Drucker-Prager yield criterion for plastic flow.

The study at hand is performed as part of ESA’s PhSR mission and investigates the efficiency of the rotary brush sampling mechanism by numerical modeling with the latter SPH approach. As part of ESA’s Mars Robotic Exploration technology development programme, a prototype of the rotary brush sampling tool was also tested under laboratory conditions on Earth and will be tested in microgravity on a parabolic flight (Ortega Juaristi et al., 2015). We developed a two-dimensional SPH model for one of the discussed sampling tool designs and performed several simulations with varying material properties of the simulated regolith. The main goal was to show the feasibility of the numerical method for the application of the sampling mechanism and to provide relative values for the sampling rate depending on the environmental properties. At first, we will present the physical model for regolith in the next section, followed by the numerical model and the implementation. In Section 4, we will show and discuss the simulations of the sampling process. We will give a conclusion in Section 5.

2 Physical model for regolith

The dynamical behaviour of regolith is usually described by a special constitutive model for soil. For the simulations carried out in this study, we applied the elastic-plastic model introduced by Bui et al. (2008). We will provide a short version of their derivation of the model in the following to point out the differences to our standard model for solid bodies. For the comprehensive treatise, we refer to their complete analysis in Bui et al. (2008). For an elastic-plastic material, the total strain rate tensor $\dot{\varepsilon}^{\alpha\beta}$, given by the derivation of the velocity $\mathbf{v}$ with respect to the coordinate system for small deformations

[TABLE]

can be written as the composition of a purely elastic and a totally plastic strain rate tensor

[TABLE]

The elastic strain rate tensor $\dot{\varepsilon}^{\alpha\beta}_{\mathrm{e}}$ is calculated by the three-dimensional version of Hooke’s law

[TABLE]

Here, $\dot{S}^{\alpha\beta}$ denotes the deviatoric stress rate tensor, $\dot{\sigma}^{\alpha\beta}$ is the stress rate tensor, $\mu$, $E$, and $\nu$ are the material dependent parameters shear modulus, Young’s modulus and Poison’s ratio, respectively, and the Einstein sum convention is used ($\dot{\sigma}^{\gamma\gamma}=\mathrm{tr}\left({\dot{\sigma}}\right)$).

The plastic strain rate tensor is given by the following relation of the rate of change of the plastic multiplier $\lambda$ and the plastic potential function $g$.

[TABLE]

If the plastic potential $g$ is equal to the yield function $f$ of the material, the flow rule is called associated, and otherwise it is called non-associated. The plastic multiplier $\lambda$ has to satisfy the conditions of the yield criterion: $\lambda=0$ for $f<0$ or $f=0$ and $\mathrm{d}f<0$, which corresponds to elastic or plastic unloading and $\lambda>0$ for $f=0$ and $\mathrm{d}f=0$, which corresponds to plastic loading. Following Bui et al. (2008), the two expressions for the elastic and the plastic strain rate tensors can be substituted into the equation for the strain rate tensor

[TABLE]

which finally yields an expression for the general stress-strain relationship for an elastic-plastic material

[TABLE]

Here, $K$ denotes the bulk modulus, which is given by the following expression

[TABLE]

In order to calculate the stress in the regolith, one needs an expression for the rate of change of the plastic multiplier $\dot{\lambda}$. The general formulation for $\dot{\lambda}$ is given by

[TABLE]

Hence, as soon as the yield function $f$ and the plastic potential function $g$ are known for a specific material, the rate of change of the plastic multiplier $\dot{\lambda}$ can be computed, and consequentially, the stress rate tensor by the use of equation (6). The stress rate tensor is then integrated and the acceleration due to the stress is determined by the equation for the conservation of momentum

[TABLE]

where $\varrho$ is the density of the regolith. The term $k_{\mathrm{ext}}$ corresponds to a specific volume force, e.g., an additional external gravity term.

The evolution of the mass density of the material is governed by the continuity equation

[TABLE]

In summary, the set of partial differential equations (6,10,9) is sufficient to describe the dynamical behaviour of a solid, plastic body. The specific plastic behaviour of the material is given by the yield function $f$ and the plastic potential function $g$. The plastic behaviour of granular media can be described by the model introduced by Drucker and Prager (1952). In the Drucker-Prager model, the yield function is given by the following relation between the first and second invariants of the stress tensor

[TABLE]

The invariants are given by the expressions

[TABLE]

The two material constants $\alpha_{\phi}$ and $k_{c}$ are called Drucker-Prager’s constant and are related to the Coulomb’s material constants cohesion $c$ and angle of internal friction $\phi$. The dependence between the four material parameters is different in plane strain and plane stress conditions. Throughout this study, we focus on plane strain conditions, where the relation is given by (see Bui et al. 2008)

[TABLE]

In addition to the yield function $f$, the plastic potential function has to be determined to specify the stress-strain relationship. In the simulations carried out in this project, we use the non-associated flow rule

[TABLE]

where $\psi$ denotes the dilatancy angle. We can write the derivative of the plastic potential function with respect to the stress tensor as

[TABLE]

In the case of a non-associated plastic flow rule given by equation (14), the derivative of the plastic potential function reads

[TABLE]

and the derivative of the yield criterion function is given by

[TABLE]

By the use of equations (8) and (16), we can finally write the complete stress-strain relationship as

[TABLE]

The rate of change of the plastic multiplier is given by

[TABLE]

or, if the stress state is not on the yield surface (e.g., $f<0$), $\dot{\lambda}=0$. In the simulations performed throughout of this work, we set the dilatancy angle $\psi$ to zero, which means there is no volume change of the material due to shearing deformations. For the non-associated flow rule eq. (14) follows the expression of the potential function

[TABLE]

In order to have a stress rate that is invariant under rotations of the reference frame, we apply the Jaumann stress rate (see, e.g., Gray et al. 2001 and references therein) and write the stress change with the help of the rotation rate tensor $R$ as

[TABLE]

where the components of $R$ are calculated according to

[TABLE]

Finally, we have the complete set of equations to describe the motion of regolith for a zero dilatancy angle

[TABLE]

In contrast to the regolith model, we need an additional equation of state (EoS) to model the materials of brushes and the collecting tray. For these objects, we use the Murnaghan EoS (see, e.g., Melosh 1996), which reads

[TABLE]

Here, $K_{0}$, $\varrho_{0}$ and $n$ are material parameters denoting the bulk modulus at zero pressure, some scaling parameter and the density at zero pressure, respectively.

3 Numerical method and implementation

We have implemented the set of equations from Section 2 in our existing code miluphCUDA (Schäfer et al., 2016) as a module to simulate regolith within the SPH framework.

Smooth particle hydrodynamics is a meshless Lagrangian particle method that was first introduced by Lucy (1977) and Gingold and Monaghan (1977) for the solution of hydrodynamic equations for compressible flows in astrophysical applications. The SPH method was later extended to model solid bodies. This work was pioneered by Libersky and Petschek (1991) with various improvements from the same source later on (Randles and Libersky, 1996; Libersky et al., 1997). The first astrophysical application of SPH with strength of material was by Benz and Asphaug (1994). Our implementation of the SPH equations in the code mainly follows their work. For a comprehensive introduction to the basic SPH idea and algorithm, we refer to the excellent publications by Benz (1990) and Monaghan (1992). In this section, we present the regolith module and SPH equations that we have implemented in our code for this project. For a complete description of the implementation, we refer to the code paper (Schäfer et al., 2016). The main difference to the usual approach to model solid bodies is the lack of an explicit equation of state for the regolith material. Instead, the time evolution of the whole stress tensor and not only the deviatoric stresses is calculated.

In the following, roman indices $a$, $b$ denote particle indices and the sum in a SPH equation runs over all of the interaction partners.

3.1 SPH equations

The conservation of mass is given by the continuity equation. The change of the density of particle $a$ is given by

[TABLE]

where $W_{ab}$ denotes the kernel function for the interaction of particle $a$ and particle $b$. Throughout our study, we applied the cubic B-spline kernel introduced by Monaghan and Lattanzio (1985). The kernel function with the normalisation for two dimensions is written as

[TABLE]

where $r$ is the distance between two interacting particles and $h$ denotes the smoothing length. The first derivative of the kernel is given by

[TABLE]

In order to calculate the acceleration for particle $a$, we use the SPH representation of the equation of motion

[TABLE]

We add some additional artificial viscosity terms to eq. (34) to prevent particles from unphysical mutual penetration and to dissipate kinetic energy in the shock front of shock waves. The artificial viscosity terms were introduced by Monaghan and Gingold (1983). They can be regarded as an additional artificial pressure term in the equations for the conservation of momentum. The additional pressure is calculated for each interaction pair $a$ and $b$ as follows:

[TABLE]

where $\alpha$ and $\beta$ are free parameters that determine the strength of the viscosity, and $\overline{\varrho}_{ab}$ and $\overline{c}_{\mathrm{s}_{ab}}$ are the averaged quantities for density and sound speed for the two interacting particles, $\overline{\varrho}_{ab}=(\varrho_{a}+\varrho_{b})/2$ and ${\overline{c}_{\mathrm{s}}}_{ab}=({c_{\mathrm{s}}}_{a}+{c_{\mathrm{s}}}_{b})/2$. The term $\nu_{ab}$ is an approximation for the divergence and is given by

[TABLE]

Here, $\overline{h}_{ab}^{2}$ is the average of the smoothing length of particle $a$ and $b$. The term $\varepsilon_{v}\overline{h}_{ab}^{2}$ in the denominator prevents divergence of $\nu_{ab}$ for particles with small separations. This artificial pressure is added to the acceleration of the particle

[TABLE]

with $\Pi^{\star}_{ab}=\Pi_{ab}$ for $(\mathbf{v_{a}}-\mathbf{v_{b}})\cdot(\mathbf{x_{a}}-\mathbf{x_{b}})<0$ and [math] elsewhere.

Since the model for the regolith does not depend on the internal energy and the Murnaghan EoS is applied for the solid material, we do not have to solve the energy equation.

In order to determine the evolution of the stress tensor, we calculate the strain rate tensor according to eq. (1). The SPH representation of the equation reads

[TABLE]

With the use of eq. (38) we can calculate the rate of change of the plastic multiplier given by eq. (8) and finally the evolution of the stress tensor.

We have used miluphCUDA’s default integrator, a Runge-Kutta 2nd order with adaptive time step, to integrate the evolution of the density, stress tensor, velocities and positions of the SPH particles. The relative accuracy of the integrator was set to $10^{-6}$.

3.2 Validation and test run

In order to validate our implementation of the model by Bui et al. (2008), we performed a comparison simulation of a two dimensional experiment of soil collapse. They used small aluminium bars of diameters $1\text{\,}\mathrm{m}\mathrm{m}1.5\text{\,}\mathrm{m}\mathrm{m}$, length $50\text{\,}\mathrm{m}\mathrm{m}$ and a density of $2650\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$ to model soil. The bars were initially arranged into a rectangular area of size 20\text{\,}\mathrm{cm}$\times$10\text{\,}\mathrm{cm}, surrounded by four flat solid walls on a flat surface. For their soil model, they estimated a friction angle of $19.8\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$, a constant Poisson’s ratio of $\nu=0.3$ and a bulk modulus $K=$700\text{,}\mathrm{k}\mathrm{P}\mathrm{a}$$. The experiment starts by removing one of the four solid walls, leading to the collapse of the soil.

The parameters for the artificial viscosity were $\alpha=0.3$, $\beta=0.7$, a fixed sound speed of $c_{s}=$600\text{,}\mathrm{m}\text{,}{\mathrm{s}}^{-1} and no-slip boundary conditions were applied. The initial number of particles was $2926$ ($77\times 38$), the initial particle separation was $\Delta=$2.582\text{\times}{10}^{-03}\text{\,}\mathrm{m} and the smoothing length was set to $2.5\times\Delta$. The initial and final locations of the SPH particles are shown in Figure 1. The plot additionally shows the outcome from the cliff collapse experiment by Bui et al. (2008) as red connected data points. The simulation successfully reproduces the shape and length of the collapsed cliff.

4 Numerical simulations of the sampling process

In this section, we present and discuss the simulations regarding the sampling process with the device prototype.

4.1 Investigated parameter space

We have chosen the material properties for the simulations to fit into the range given by Ortega Juaristi et al. (2015) for the regolith simulant that they used for the experiments. The material parameters for regolith that remain unchanged for all simulations are bulk modulus $3.33\text{\,}\mathrm{M}\mathrm{P}\mathrm{a}$, shear modulus $370\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}$, and bulk density $1.5\text{\,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$. We have varied the angle of internal friction and cohesion. An overview of all simulations performed for this study is provided in Table 1.

4.2 Initial setup, model for the brushes and the collecting tray

In order to raise the spatial resolution, we have developed a two dimensional model of the existing CAD data of the sampling tool (see Figure 2).

The initial setup and geometry of the sampling tool are shown in Figure 3. The collecting tray is located above two brushes. The brushes ought to shovel regolith from the surface into the tray. The collecting tray has a height of $63\text{\,}\mathrm{m}\mathrm{m}$ and a width of $82\text{\,}\mathrm{m}\mathrm{m}$. On the bottom of the tray is a opening with a diameter of $26\text{\,}\mathrm{m}\mathrm{m}$. The final sampling tool CAD design includes a closure for the opening at the bottom of the collecting tray.

The two brushes are modelled with 12 individual bristle, separated by $30\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$. The inner part of a brush is a cylinder with a radius of $15\text{\,}\mathrm{m}\mathrm{m}$, the bristles sit on top of this inner part and have a length of $20\text{\,}\mathrm{m}\mathrm{m}$ each. On the outer part of each brush is a circular shielding, covering an opening angle of $80\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$. The purpose of the shielding is to lower the contamination of parts of the landing device during the sampling process and reduce staining of cameras and sensor devices. The distance between the pivot point of the brushes is $95\text{\,}\mathrm{m}\mathrm{m}$.

In the collecting phase, the orange brush rotates counterclockwise and the yellow one clockwise. During the clearing phase, both brushes rotate in the opposite direction. We have applied the Murnaghan EoS for the brushes and the collecting tray including the shielding with the following material parameters: bulk modulus $24.2\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$, bulk density at zero pressure $2.7\text{\,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$, and Murnaghan parameter $n=6$.

Initially, the lowest part of the brush is $8.4\text{\,}\mathrm{mm}$ above the surface. Due to the spatial resolution given by the smoothing length, the brush touches the surface at a distance of $3\text{\,}\mathrm{mm}$. The brushes and the collecting tray move at constant speed of $v_{y}=$-5\text{,}\mathrm{cm}\text{,}{\mathrm{s}}^{-1}$$ during the first $0.2\text{\,}\mathrm{s}$ after simulation start and remain at their final location after this time until the end of the simulation. Additionally, the brushes start rotating from the beginning and their rotation speed is fixed. The total simulation time is $1\text{\,}\mathrm{s}2\text{\,}\mathrm{s}$. The number of particles is approximately $280\,000$ and the smoothing length is set to $3\text{\,}\mathrm{m}\mathrm{m}$ (for the red region of the regolith) in each run.

In order to analyze the efficiency of the collecting process and the dependency on the material properties of regolith, we monitor the number of particles that get captured in the collecting tray during the collection phase. This study does not aim to provide absolute values for collecting rates, we rather investigate the relative difference and the influence of material properties on the collection process. Hence, the number of collected particles in the tray is a sufficient indicator for the quality of the sampling rate for the purpose of this study.

In the remaining parts of this section, we present the results of 13 selected simulations of the sampling process on Earth and Phobos. We will discuss the influence of particular material parameters such as cohesion and angle of internal friction (see Table 1 for the parameters in each run) on the collection rate.

4.3 Sampling on Earth

We have started our investigations with the sampling process under Earth’s gravity conditions to compare with the experiments performed by Ortega Juaristi et al. (2015). They conducted a test campaign with an early breadboard design to narrow down bristle material, configuration, motor torque and rotation speed. In total, more than 300 collecting sequences were performed.

The minimum required rotation speed to collect material in Earth’s gravity field is investigated by performing four different simulations with three different rotation speeds. Figure 4 shows the particle distribution after $1\text{\,}\mathrm{s}$ simulation time for these runs. The colour code indicates the speed of the particles. The runs E0, E2 and E3 use cohesionless regolith with the same angle of internal friction of $42\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ for three different rotation speeds of the brushes, $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}300\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}500\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$. Run E1 has the same parameters as E0 but a cohesion of $500\text{\,}\mathrm{Pa}$. The particle distribution in between the bristles for runs E0 and E1 and the corresponding speed of the particles are alike. In the cohesionless case, the brushes have already left deeper traces on the surface. In both runs, the collecting tray remains empty, only single particles stray in. Runs E2 and E3 with the higher rotation speeds of $300\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ and $500\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ display higher velocities of the particles. Additionally, the impact on the surface morphology raises with the rotation speed: E3 with the highest rotation speed reveals a cavity almost $20\text{\,}\mathrm{cm}$ in diameter. While some parts of the material from the cavity were transported to the collecting tray, some parts were also pushed sidewards. In both runs with the higher rotation speed, the sampling process works effectively and material from the surface gets lifted into the tray.

To analyze the collecting process quantitavely, the number of particles in the collecting tray during the simulation time is tracked and plotted in Figure 5. A rotation speed of the brushes of $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ is not enough to transfer material into the collecting tray during the first second (run E0). Adding a little bit of cohesion to the regolith enhances the sticking between the particles and leads to a slightly higher collecting rate in this case (run E1).

Using a speed of $300\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ (run E2), material can be collected. However, material that reaches the tray after $0.25\text{\,}\mathrm{s}$ can also fall back through the opening and the actual collecting phase starts at $0.5\text{\,}\mathrm{s}$, as soon as the number of particles in the tray rises monotonically. At a speed of $500\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$, the material transport from the surface to the tray becomes fast enough to prevent regolith in the tray from falling back down, and the number of collected particles increases continuously.

In agreement with the experiments, we find that a rotation speed significantly lower than $500\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ does not lead to mass transfer into the tray. Ortega Juaristi et al. (2015) applied rotation rates between $500\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}1000\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ in their study in order to achieve sampling.

4.4 Sampling on Phobos

We use the same material properties and change the gravitational acceleration to a mean value on Phobos’ surface. The gravitational acceleration on its surface varies between $4\text{\times}{10}^{-3}\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}^{2}8\text{\times}{10}^{-3}\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}^{2}$, and we use a mean value of $5.7\text{\times}{10}^{-3}\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}^{2}$.

The positions of the particles after $1\text{\,}\mathrm{s}$ can be seen in the plots of Figures 6 and 7. The colour code indicates the speed of the particles. Runs P0-3 (Figure 6) have the same brush rotation speed of $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$, the material properties of the regolith were varied: The runs P0, P1 and P3 include an angle of internal friction of $42\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ and a cohesion of $0\text{\,}\mathrm{Pa}$, $500\text{\,}\mathrm{Pa}$ and $1\text{\,}\mathrm{kPa}$, respectively, and P2 a lower angle of internal friction of $20\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ and a cohesion of $1\text{\,}\mathrm{kPa}$. The runs with $42\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ and cohesion show some small clumping of the particles in the tray (P1 and P3). In the run with $20\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ and $1\text{\,}\mathrm{kPa}$, the opening of the collecting tray gets blocked by a smaller pile, and material that entered the tray sticks for the most part on the ceiling.

The rotation speed of the brushes was varied in runs P4-7 (Figure 7): Runs P5 and P6 have a rotation speed of $500\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$, P4 of $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ and P7 of $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ in the opposite direction. Runs P5-7 include an angle of internal friction of $42\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ and P4 of $60\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$. The runs P5 and P7 are without cohesion and runs P4 and P6 include a cohesion of $1\text{\,}\mathrm{kPa}$. Note the different velocity scaling compared to Figure 6. Both runs with the highest rotation rate (P5 and P6) show the largest impact on the surface. In the case with a cohesion of $1\text{\,}\mathrm{kPa}$ (P6), the particles fill the space in the collecting tray and in between the bristles of both brushes. On both sides, larger piles of regolith get pushed away. The same simulation setup without any cohesion (P5) shows a less dramatic picture with slower velocities. The cavity below the brushes is deeper than in the other runs and less particles stick in between the bristles.

To analyze the collecting process in more detail, the time evolution of the amount of sampled material in the collecting tray on Phobos’ surface is shown in Figure 8. Higher rotation speeds yield higher sampling rates. In the case of $500\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ and a cohesion of $1\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}$ (blue curve, run P6), the collecting tray is filled after $0.5\text{\,}\mathrm{s}$ after simulation start, upper right panel), and subsequently the amount of regolith in the tray decreases.

4.4.1 Varying angle of internal friction and cohesion

The results for the three simulations with varying angle of internal friction on Phobos are shown in Figure 9. The internal friction leads to higher shear in the regolith, the force between the particles raises, which eventually supports the sampling process, since the brushes break off larger clumps at first contact with the surface. In the case of a small internal friction (P2), the collecting tray is much less filled compared to the simulation with $42\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ (P3). In this particular simulation, we observe that one small clump of regolith sticks at the left end of the opening and blocks it partially (see Figure 6).

A further increase of the internal friction leads to an even higher collection rate (P4). However, the value $60\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ is higher than the expected value for regolith on Phobos. The most probable range of lunar soil friction angle is about $30\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}50\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ (Mitchell et al., 1972) and Phobos’ regolith is presumably in a range comparable to lunar and martian soil.

The influence of the cohesion on the collecting rate is indistinct. On the one hand, cohesion supports the effect of clumping and larger chunks of regolith are transported by the brushes, while on the other hand cohesion also hinders the brushes to lift the chunks from the surface. We see a lower collection rate in the Phobos runs with $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ rotation speed for cohesion of $500\text{\,}\mathrm{P}\mathrm{a}$ and $1\text{\,}\mathrm{kPa}$ compared to the cohesionless case for the first second of the simulation, see curves P0, P1 and P3 in Figure 8. The run with $1\text{\,}\mathrm{kPa}$ cohesion, however, yields a slightly larger collection rate than the run with $500\text{\,}\mathrm{Pa}$.

Using a five times higher rotation speed of $500\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ brings particles to the tray already after $0.19\text{\,}\mathrm{s}$. Interestingly, a higher cohesion in combination with a higher rotation speed at $500\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ yields to a faster collection rate (run P6) compared to no cohesion (run P5), whereas for a lower rotation speed of $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$, cohesion lowers the collection rate (runs P0, P1, P3). This is caused by the fact that higher cohesion leads to larger clumps of regolith which requires more energy to be lifted by the brushes as can be seen quantitatively in Figure 6 (panels P0, P1, P3 with cohesion $0\text{\,}\mathrm{Pa}$, $500\text{\,}\mathrm{Pa}$, and $1\text{\,}\mathrm{kPa}$, respectively).

In the case of the highest cohesion and highest rotation speed (run P6), the tray reaches its maximum charging level already after $0.51\text{\,}\mathrm{s}$. Subsequently, material from the tray falls back through the opening until eventually the pressure in the tray and the pressure in the space between the brushes balance. Since the material stream from the surface is not unlimited as the brushes stop their vertical motion after $0.2\text{\,}\mathrm{s}$, the number of particles in the tray slowly decays until the end of the simulation.

4.4.2 Clearing phase

In the clearing phase, the brushes rotate the other way. In this manner, potential obstacles on the surface of Phobos may be cleared and material from below the top layer can be collected in the subsequent collecting phase. We compare the collecting rate for a rotation speed of $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ to the rotation speed of $-100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ in Figure 10, where the time evolution of the number of particles in the collecting tray for the clearing phase simulation is shown. In the first case, the first material reaches the tray after about $0.35\text{\,}\mathrm{s}$, independent from the cohesion. In the latter case, in which no material at all should be collected, we still find small number of particles that reach the tray. These particles are transported by sticking between the bristles of each brush. Their path to the collecting tray is longer than for the particles that are lifted between the two brushes and they arrive later at the opening of the tray. The first particles enter the tray after $0.53\text{\,}\mathrm{s}$.

4.5 Comparison Phobos vs. Earth

The time evolution of the collected number of particles in direct comparison between Earth and Phobos is shown in Figure 11.

In contrast to the sampling on Earth, regolith can be transferred to the collecting tray on Phobos even for the lowest investigated rotation speed of the brushes ($100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$). At this speed, the first material enters the tray after $0.35\text{\,}\mathrm{s}$ on Phobos (run P0), compared to $0.75\text{\,}\mathrm{s}$ on Earth (run E0). However, the sampling rate on Earth is significantly lower, and after $1\text{\,}\mathrm{s}$ there are less than 20 particles in the tray. At this time, there are already 400 collected particles in run P0.

5 Conclusions

The simulations give some promising insight into the possible application of SPH to regolith modelling. The basic functionality of the sample process could be investigated for different material properties of regolith. Since our model cannot provide absolute data for the collecting phase, we compare the relative values between the individual simulations. We have simulated the collection phase of the brush sampling tool for varying cohesion and angles of internal friction on Earth and on Phobos.

For non-cohesive soil, Earth’s gravity is high enough to prevent a high collecting rate of dust material at brushes speed lower than $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}300\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$. If the soil is slightly cohesive (about $500\text{\,}\mathrm{P}\mathrm{a}$ cohesion), the situation improves. For Phobos’ gravity and lunar soil properties ($42\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$, $500\text{\,}\mathrm{P}\mathrm{a}$), the collecting rate is more than ten times higher at $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ than for Earth’s gravity.

We found one case, in which the collecting tray hole was blocked by particles ($20\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$, $1\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}$, $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$, Phobos g, run P2) and the collecting rate was slowed down in the following. This might come from the 2D geometrical setup and may not happen in 3D simulations. However, to obtain a reasonable high spatial resolution, we chose a 2D geometrical setup. The speed of the dispersed dust particles is mainly determined by the rotation rate of the brushes and depends only slightly on the material properties angle of internal friction and cohesion.

The escape velocity of Phobos ($11.39\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$) was only reached in one simulation ($42\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$, $1\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}$, $500\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$, run P6) after $0.6\text{\,}\mathrm{s}$. However, this resulted from some escaping particles that were pressed through the slit between the collecting tray and the shield of the brushes. This slit is not in the CAD data and only in the SPH setup.

The shielding around the brushes is very effective: We do not see any material flying directly upwards to the spacecraft. Qualitatively, the simulations do not show high dust contamination of the spacecraft from the sampling process. The maximum settling time of the dispersed dust in our simulations can be estimated to be between $40\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{u}\mathrm{t}\mathrm{e}\mathrm{s}80\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{u}\mathrm{t}\mathrm{e}\mathrm{s}$ (for Phobos g varying from $4\text{\times}{10}^{-3}\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}^{2}8\text{\times}{10}^{-3}\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}^{2}$).

During the clearing (brushes turning the other way round), we find some particles that get eventually in the collecting tray. These collected particles are transported within one single brush only and enter the tray after $0.53\text{\,}\mathrm{s}$. Cohesion enhances this effect since it supports the sticking of material within the bristles of a single brush. Therefore, it is recommended to keep the tray closed during the clearing phase, and additionally a process to clean out the material in the bristles of the brushes after the cleaning phase may be required to avoid material from the clearing phase to be sampled.

Our analysis is based on a series of simulations that do not cover the entire parameter space. It is desirable to broaden the parameters even more and additionally vary the density of the regolith or add a layered structure to the surface model. Based on our experience with this study, we suggest following further steps to improve the sampling model simulations

Subsequent simulations of the collecting process need higher spatial resolution and should be three-dimensional. Especially, the latter requirement is obvious if one wants to calculate absolute values for the collection rate, since the opening at the bottom of the collecting tray cannot be modelled correctly in two dimensions. Recently, a study about the Mars Exploration Rover wheel performance was published by Johnson et al. (2015). They successfully applied the three-dimensional discrete element method (DEM) to investigate wheel slip conditions for various contact friction coefficients. They underline the need for three-dimensional simulations in their study to accomplish good agreement with experimental data.

• 2.

Simulations with different angles between sampling brushes and surface: The simulations for the sampling process presented in this work feature brushes, where the axes of the brushes are oriented parallel to Phobos’ surface and move only up and down. Subsequent simulations may be more realistic if the sampling tool moves with more degrees of freedom. Furthermore, we have simulated only one path of the brush axes: from the starting point going with constant speed in vertical direction for $0.2\text{\,}\mathrm{s}$ and halting. A proper model of the brushes may include a feedback reaction of the brushes depending on the sampling rate. The prototype by Ortega Juaristi et al. (2015) had the additional possibility to manually adjust the angle between the surface of the regolith and the brushes with three different values. This could be added quite simply to our initial setup for a future study.

• 3.

In this study, we have applied the soil model presented by Bui et al. (2008). In a further study, we will compare this model to another model for granular media using continuous mechanics, where the soil can be treated as a viscous material (Ulrich et al., 2013).

In summary, in this study we have successfully applied the soil model by Bui et al. (2008) for regolith and built a simulation testbed for modeling the sampling process considered for PhSR. Although absolute values for the sampling rate cannot be provided, the dependencies of the collecting process on various material properties such as the cohesion and angle of internal friction can be studied in detail.

Moreover, in agreement with the laboratory experiments with the sampling prototype by Ortega Juaristi et al. (2015), we find that a rotation speed of $100\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ for the brushes is not sufficient to transfer regolith from Earth’s surface into the tray.

In this study we cannot provide absolute numbers for sampling rates, and give estimations about the influence of the mechanical properties of the regolith on the sampling rate. Using an improved three dimensional model and higher spatial resolution, future numerical models will allow for more precise conclusions.

Acknowledgements

This material is based on the work performed on behalf of the industrial Phobos Sample Return Phase A study team of Airbus Defense and Space under ESA contract 4000114302/15/NL/PA (technical officer: Thomas Voirin). We thank Gerhard Wurm and an anonymous reviewer whose comments and suggestions helped to improve and clarify this manuscript.

All plots in this publication have been made by the use of the matplotlib package by Hunter (2007).

TM appreciates support by the FWF Austrian Science Fund project S 11603-N16. CS wants to thank Daniel Thun for helpful discussions during the course of this project.