Available online at www.sciencedirect.

com

ScienceDirect
Advances in Space Research 70 (2022) 3231–3248
www.elsevier.com/locate/asr

Investigating particle-particle electrostatic effects on charged lunar

dust transport via discrete element modeling
Hao Wang a,⇑, James R. Phillips III b, Adrienne R. Dove c, Tarek A. Elgohary a
a
University of Central Florida, Department of Mechanical and Aerospace Engineering, 12760 Pegasus Drive, Orlando, FL 32816, USA
b
NASA Kennedy Space Center, Electrostatics and Surface Physics Laboratory, NASA KSC Mailcode UB-G, Kennedy Space Center, FL 32899, USA
c
University of Central Florida, Department of Physics, 4111 Libra Drive, Physical Sciences Bldg. 430, Orlando, FL 32816, USA

Received 19 March 2022; received in revised form 3 August 2022; accepted 30 August 2022
Available online 6 September 2022

Abstract

NASA surface exploration missions have always seen negative effects of dust including the Apollo missions. The astronaut-witnessed
unusual behavior of the dust particles that surround the vehicle after engine cutoff has the potential to have more of an influence on
surface systems dust loading than the high velocity lunar rocket plume ejecta in the landing process. The levitation and transport of
the fine components of regolith on lunar surface has been linked to electrostatic effects and electric field, but so far there is no accurate
model considering the inter-particle electrostatic interactions, especially when the particles are charged by rocket plume or other mechan-
ical interactions due to exploration activities. This study is proposed to investigate the dynamics of charged lunar regolith with a discrete
element method (DEM) approach focusing on the inter-particle interactions and contact charge transfer. The grain dynamics is coupled
with mechanical and electrical particle interactions, and both short- and long-range interactions between spherical particles are incorpo-
rated. A tribo-charging model based on instantaneous collisions between particles is adopted and validated by comparing the simulation
results to existing experimental data. Sensitivity analysis is conducted to quantify the effects of initial charge, tribo-charging, and E-field
on transport of lunar dust based on JSC-1 simulants with a radius of 50 lm. DEM simulations are also conducted in a near realistic lunar
environment with the estimations of initial conditions that shows the difference in position and velocity distributions between charged
particles and uncharged particles. The results indicate that the charged dust particles have higher dispersion of position and velocity by
several orders of magnitude due to electrostatic effects. This provides a potential explanation for the phenomena of the approximately 30
s dust lofting following Apollo Lunar Module landing.
Ó 2022 COSPAR. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/
by-nc-nd/4.0/).

Keywords: Lunar dust; Discrete element method; Planetary science; Aerospace

1. Introduction have subsided took much longer than expected for ballistic
trajectories, which is best noted in the Apollo 14 landing
Dust problems raise significant concerns in planetary video, in which visible dust takes nearly thirty seconds to
surface exploration for NASA and other agencies. This dissipate (Metzger et al., 2011). This astronaut-witnessed
starts with landings. For instance, in the ‘‘clearing” stage phenomenon has the potential to have more of an influence
of Apollo landings, settling of the grains after plume effects on surface systems dust loading than the high velocity
lunar rocket plume ejecta that is thought of as the primary
product of the landing process. A prevailing hypothesis
⇑ Corresponding author.
attributes this unusual behavior to the accumulated charge
E-mail addresses: wanghao@knights.ucf.edu (H. Wang), james.r.
phillips.iii@nasa.gov (J.R. Phillips III), adove@ucf.edu (A.R. Dove), on the particles due to the plume and due to interactions
elgohary@ucf.edu (T.A. Elgohary).

https://doi.org/10.1016/j.asr.2022.08.080
0273-1177/Ó 2022 COSPAR. Published by Elsevier B.V.
This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

with the surrounding natural environment (Metzger et al., focus on dust charging behavior, especially to understand
2011; Wang et al., 2021). natural lofting processes. Starting from the electrostatic
The lunar surface is collisionally evolved as a result of effect on dust grains, Hartzell and Scheeres (2011) evalu-
meteoroid bombardment, thermal, and space weathering ated the electric field strength required for launching dust
processes. The lunar regolith consists of particles in a particles from the surface, for which it must overcome
broad size distribution, including dust, which is usually the cohesive force between grains. However, significant
defined as particles at the level of 100 lm and below uncertainty exists in the method of launching small dust
(Colwell et al., 2007). The lunar surface is not a static envi- particles, and no particle–particle electrostatic interactions
ronment - grains are loosed and ejected from the surface are taken into account. Yeo et al. (2021) investigated the
due to natural processes as well as from human activities. dynamics of charged dust above the surfaces of airless bod-
Thus, lunar dust behavior must be adequately character- ies across a wide range of dust grain size and of the surface
ized and addressed for safe exploration and In-Situ gravity to identify the lofting heights using initial condi-
Resource Utilization (ISRU) (Taylor et al., 2005; Heiken tions from recent laboratory results. However, these initial
et al., 1991). conditions are derived from the experiment of simulant
In addition to the mechanical complexities of the lunar grains charging by UV light or electron beam, and the cal-
regolith, the plasma and electrical environment at the lunar culation of the charge and the trajectories are based on the
surface presents challenges to operating on the Moon. UV and plasma environment without inter-particle interac-
While the surface of the Moon is exposed to the vacuum tions. An electrostatic lofting model based on adhesion was
of space, it is also immersed in the Solar Wind plasma. This proposed in Wang et al. (2020) for accelerated separation
results in surface charges, potential differences, and electric process and ballistic projection process with lunar surface
fields in the near-surface region, with strong variations over charging, but the interaction force between the lofted par-
the course of the lunar day and night (Colwell et al., 2007). ticles is neglected.
The charging currents on the lunar surface come from four A few DEM models incorporated inter-particle interac-
main sources: photoemission of electrons, plasma elec- tions with a different focus. In Hou et al. (2016), the adhe-
trons, plasma ions, and secondary electrons. There remain sion mechanism, especially Van der Waals force among
significant uncertainties in the charging processes and their inter-particle interactions, is mainly discussed. The inter-
natural variations (Stubbs et al., 2007). In this environ- particle electrostatic interactions were also included in
ment, dust particles accumulate charge and are strongly modeling and simulation of electrostatic lunar dust collec-
affected by electrostatic forces that compete with the grav- tor to evaluate the collection efficiency and determine the
itational forces above the lunar surface, leading to unusual external E-field required to collect dust particles (Afshar-
behavior including levitation and transport across the sur- Mohajer et al., 2011; Afshar-Mohajer et al., 2012). These
face (Colwell et al., 2007). Additionally, grains in motion studies investigated the dust particles collection on the
will be affected by contact forces, including mechanical charged conducting plates with provided E-field, neither
and electrostatic forces that arise from effects such as tri- in the natural environment nor under the influence of
bocharging. The levitation and transport of the fine com- rocket plume or other exploration activities.
ponents of regolith on the lunar surface has been linked In this paper, a discrete element method (DEM)
to electrostatic effects and electric field, the same as on approach is developed to explore charged particle behavior
Martian and/or asteroid surfaces (Colwell et al., 2007; on the lunar surface, focusing on the inter-particle interac-
Colwell et al., 2005; Hughes et al., 2008; Hartzell et al., tions of the macroscopic grains, especially when lofted by
2021; Patel and Hartzell, 2021). plume or other mechanical interactions due to exploration
In situ observations and samples returned from the activities. So far, no generalized granular dynamics model
lunar surface provide a picture of lunar regolith mechanical has been built that also fully takes into account electro-
behavior, but observations of these electrostatic properties static interactions.
are sparse. A combination of experimental work and mod- DEM techniques have been employed to simulate rego-
eling is needed to build a full picture of regolith grain lith behavior in a number of relevant environments, includ-
behavior in relevant environments, and this is essential ing the lunar surface, the 3-D axisymmetric triaxial
for continued lunar (and other dusty planetary surface) compression of lunar regolith simulant by Hasan and
exploration. Alshibli (2010), and with Martian rovers by Knuth et al.
Experiments have been undertaken to understand grain (2012). In another highly relevant study, DEM techniques
charging with relevant materials. For instance, Sternovsky were also used to investigate the erosion flux of lunar soil
et al. (2002) explored contact charging of lunar and Mar- from rocket exhaust plume in Berger and Hrenya (2017).
tian simulants, providing useful information on relevant More broadly, DEMs have been widely used in molecular
charging behaviors. Tribochargeability of lunar dust simu- dynamics (Alder and Wainwright, 1957; Alder and
lants with various properties was studied experimentally, Wainwright, 1959; Alder and Wainwright, 1960) and have
such as in Afshar-Mohajer et al. (2014) and in Forward been proven very useful in granular dynamics, with the
et al. (2009). However, no theoretical models or simula- advent of modern computational power monitoring each
tions are produced from these works. A number of models particle’s dynamics and allowing for more realistic force
3232
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

models (Kloss et al., 2012). DEM was adopted for the sim- (2008). In this method, velocities in the upcoming time step
ulation of non-smooth granular dynamics at the scale of t þ Dt=2 and positions at t þ Dt are calculated as in Eq. 2.
particle rearrangements in Radjai and Richefeu (2009),
x_ tþDt=2 ¼ x_ tDt=2 þ x
€ðxt ; x_ tDt=2 ÞDt
and was implemented to study the particle rearrangement ð2Þ
during powder compaction in Martin et al. (2003). A num- xtþDt ¼ xt þ x_ tþDt=2 Dt
ber of other works have taken various approaches to solv-
ing dust dynamics problems. Regarding the transport of For the present work, the open-source software
dust particles, Marshall et al. (2011) employed fracture LIGGGHTS is adopted and modified. The public version
mechanics to explain lofting particulate clouds of lunar of the Large-scale Atomic/Molecular Massively Parallel
dust in low ballistic trajectories rather than individual sub- Simulator (LAMMPS) Improved for General Granular
micron dust particles in high trajectories, and integrated and Granular Heat Transfer Simulations (LIGGGHTS)
models including a hybrid continuum–kinetic solver, a cou- software package provides a DEM simulation framework
pled two-phase flow model, and a model for inelastic for granular interactions. LIGGGHTS is supported by
grain–grain collisions were adopted in Morris et al. (2016). CFDEM Project and is written in C++ code with parallel
DEM techniques are highly appropriate for simulating computing (Kloss et al., 2012).
complex lunar regolith interactions, and thus plume LIGGGHTS is often implemented to solve industrial
dynamics, but it is also essential to include electrostatic manufacturing problem, such as in application for charging
effects to accurately model particle trajectories. This work system of ironmaking blast furnace (Wei et al., 2017), in
adds electrostatic effects to the well-used LIGGGHTS soft- analysis of the compaction of a cohesive material
ware (Kloss et al., 2012) to model charged particles in lunar (Ramı́rez-Aragón et al., 2018), and in simulation of magne-
environment to explain the phenomenon of the ejecta that torheological fluids (Leps and Hartzell, 2021). In this work,
surround a vehicle after engine cutoff by numerical simula- the module of short-range forces is enhanced to start the
tion of particle dynamics with inter-grain interactions. exploration of charge transfer, and more importantly,
The paper is organized as follows: description of the long-range electrostatic and gravitational forces that are
DEM techniques, improvements to LIGGGHTS, and val- essential to exploring the behavior of charged dust grains
idation of the tribo-charging model are introduced in Sec- are implemented. Both standard and screened Coulomb
tion 2. Section 3 discusses the sensitivity analysis and the potentials are integrated into the LIGGGHTS framework
DEM simulations of charged lunar dust particles with con- to provide a basis for long-range electrostatic interactions.
tact forces, Coulomb force, and charge transfer. A discus- A gravitational potential is implemented to enable inter-
sion of results and numerical errors are presented in grain gravitational interactions, but the inter-particle gravi-
Section 4, with concluding remarks in Section 5. tational force is relatively small and can be neglected in cur-
rent simulation compared to the electrostatic force which is
the main focus of the present research. The properties such
2. Discrete element method
as work function and electrical conductivity are incorpo-
rated to the library of available material characteristics to
The discrete element method (DEM) is a particle-based
investigate the charge transfer between grains. There are
numerical model computing the individual behaviour of
numerous challenges to incorporate realistic interactions
each particle where the interaction of the particles is mon-
between complex, realistic particles. Currently, grains are
itored contact by contact and the motion of the particles is
modeled to behave as if the entirety of the charge acts at
modeled particle by particle (Cundall and Strack, 1979).
the center of mass, such as conductors with spherical symme-
DEM provides a high-fidelity solution for each particle
try and insulators with homogeneously distributed charge.
dynamics; however, it is also known to be a computation-
ally intensive method (Lommen et al., 2018).
The charged particles dynamics is modeled as a multi- 2.1. Long-range interactions: Coulomb force
body dynamical system. In DEM multi-body dynamics
can be generally described as in Eq. 1. Electrostatic interactions between charged particles can
be modeled using a basic Coulomb force, which will be
m x ¼ RF ð1Þ
repulsive or attractive based on the sign of the charge.
where m is the mass, x is the position vector, and F is the When the particles are in a natural plasma environment,
€ is assumed constant over the time step Dt. The sys-
force. x or when there are additional interactions due to the
tem can also handle externally applied forces via F. charged particles themselves, there can be an additional
The resulting system of second order differential equa- ‘‘screening” of these forces, typically represented by a fac-
tions is solved numerically mostly via low-order finite dif- tor related to Debye length, which is a scale of how far a
ference methods (Matuttis and Chen, 2014). The charge carrier’s electrostatic effect persists. Based on the
Velocity-Verlet method is adopted to perform the numeri- local charge concentration, a target charged particle will
cal integration which is a widely used integration method experience a screened Coulomb force due to the presence
of second order as initially proposed by Cundall and of other charged particles in close vicinity because in reality
Strack (1979) and summarized in Kruggel-Emden et al. the long-range electrostatic effects are suppressed by the
3233
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

flow of particles and the effective Coulomb interaction

between particles is reduced to a finite range, as discussed
in both Hogue et al. (2008, 2016).
The Coulomb force, F c , is implemented between a pair
of particles i and j as in Eq. 3.
qi qj j 1
Fc;ij ¼ ð þ Þejrij ^rij rij < rc ð3Þ
4pe0 er rij r2ij

where qi and qj are the charges of the two interacting par-

ticles, r is the separation between their centers of mass, j is
the inverse of the Debye length, e0 is the permittivity of free
space, and er is the relative permittivity of the space
between the particles. To improve the computation time
of the simulation, the force is assumed to be negligible if
the particles are separated by a distance equal to or larger
than a cutoff distance, rc . While the particles have a phys-
ical size in the model, the charge is treated as if it is homo-
geneously distributed around the particle center of mass.
Particle and plasma characteristics er and j are user-
specified constants, while the charges on the particles, q,
are defined via initial conditions and then allowed to vary
in the simulation. The special case where j is zero repre-
sents a completely Coulombic interaction with no screen-
ing. Use of the j values allows refinement for more
realistic interaction values. The performance of initial force
and time to contact between two particles with varying j Fig. 1. Hertzian contact model.
has been shown in Phillips et al. (2021). To simplify the
problem, j is set to zero in this paper.
(Kloss et al., 2012) is used to calculate the contact force
between two granular spherical particles in this Hertzian
2.2. Contact interactions: Contact forces granular model, when the distance r between two particles
of radii Ri and Rj is less than their contact distance
When charged or uncharged particles come into contact d ¼ Ri þ Rj. There is no contact force between the particles
they will experience a charge transfer that is based on vari- when r > d.
ables including contact type, contact area, material proper-
ties, and particle charge state. Interactions may include F ¼ ðk n dnij  cn vnij Þ þ ðk t dtij  ct vtij Þ ð4Þ
brief contacts (bouncing), sticking, rolling, or contacts of The first term of the RHS of Eq. (4) is the normal force
multiple particles, with the area of contact, length of time between the two particles including a spring force and a
in contact, and number of contacts potentially influencing damping force while the second term is the tangential force
the distribution of charge. In order to study these effects, including a shear force and a damping force. The limit of
conduction and tribo-charging are integrated with the maximum Hertzian tangential force is governed by
Hertz contact model in LIGGGHTS to produce a new con- F t 6 xl F n , where F t is the tangential force, F n is the normal
tact model that accounts for electrostatics. Conduction and force, and xl is the coefficient of static friction.
triboelectric charging added to LIGGGHTS have been The coefficients in Eq. 4 are calculated from the material
introduced in Phillips et al. (2021). The present work will properties as:
focus on the inter-particle triboelectric charging process pffiffiffiffiffiffiffiffiffi qffiffi pffiffiffiffiffiffiffiffiffiffi
considering the insulated properties of the lunar dust. k n ¼ 43 Y  R dn ; cn ¼ 2 56b S n m P 0
The particle–particle contact forces govern how particles pffiffiffiffiffiffiffiffiffi qffiffi pffiffiffiffiffiffiffiffiffiffi
resist interpenetration and rebound off each other in colli- k t ¼ 8G R dn ; ct ¼ 2 56b S t m P 0
sions. The Hertzian model is a contact model with high fide- pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi
S n ¼ 2Y  R dn ; S t ¼ 8G R dn
lity in LIGGGHTS, which treats the particles as an elastic
material with a defined Young’s modulus and Poisson ratio b ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi
lnðeÞ
2
ffi ð5Þ
ln ðeÞþp2
(Leps and Hartzell, 2021). Fig. 1 shows the particle to parti- ð1m22 Þ
¼ ð1m Þ
2

cle contact and particle to wall interaction in the Hertzian

1
Y Y1
þ Y2
contact model. 1
G
¼ 2ð2m1YÞð1þm 1Þ
þ 2ð2m2YÞð1þm 2Þ

Given Young’s Moduli Y i ; Y j , Poisson ratio mi ; mj , and 1 2

coefficient of restitution e, Eq. 4 from LIGGGHTS

1
R
¼ R11 þ R12 ; 1
m
¼ m11 þ m12
3234
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

where dnij is normal overlap distance of the i-th and j-th charge transfer. The charge transferred by triboelectric
particles, dtij is tangential overlap distance, vnij is normal charging, Dq is calculated as in Eq. 6.
component of the relative velocity, vtij is tangential compo-
e0 Aij;max  
nent of the relative velocity, k n is elastic constant for nor- Dqij ¼ ui  uj ð6Þ
zqe
mal contact, k t is elastic constant for tangential contact,
cn is viscoelastic damping constant for normal contact, ct where Dq is the charge transferred between the particle i
is viscoelastic damping constant for tangential contact, and the particle j during a single collision, e0 is the permit-
and R is effective radius. tivity of free space, A is the contact area, z is the cutoff dis-
In the next section, the tribo-charging model in Eq. 7 tance for charge transfer, qe is the charge of an electron, u
from Laurentie et al. (2013), is combined with the Hertzian is the work function of the particle, and q is the charge on
contact model in Eq. 4, to quantify the charge transfer the particle. In order to determine the contact area, A, at
between particles. To validate the new combined model, that time step, the Hertzian contact model described above
simulations from the new model is validated against is used and the overlap distance between two particles in
numerical and experimental results shown in Laurentie the normal direction at the current time step is compared
et al. (2013). with the same variable at the previous two time steps.
When this distance is maximized, the contact area is the
maximum, representing the moment of maximum impact.
2.3. Contact interactions: Triboelectric charging In Laurentie et al. (2013), the charge transferred by tri-
boelectric charging, Dq is calculated as in Eq. 7.
In electrostatics, triboelectric charging is important as  
well as challenging with no exact model existing so far. e0 DAij;n dij
Dqij;n ¼ ui  uj  Eij;n zqe ð7Þ
Tribo-charging is a type of contact electrification that zqe kdij k
mostly takes place between different materials but can also
happen between particles of the same material with differ- where Dq is the charge transferred between the particle i
ent sizes (Carter and Hartzell, 2017). The charge amount and the particle j during time step n to time step n þ 1; e0
could depend on material, size, surface roughness, temper- is the permittivity of free space, DA is the contact area
ature, and other properties. (as defined in Eq. 8) variation between particles i and j dur-
In Hogue et al. (2008), the trajectories of triboelectri- ing time step n to time step n þ 1; z is the cutoff distance for
cally charged particles with electrostatic forces were mod- charge transfer, qe is the charge of an electron, u is the
elled in a commercial DEM software. The model was work function of the particle, Eij;n is the electric field
developed depending on a key parameter a, which is the located at the point of impact between the two particles
time constant of the charge generation. a was used as the at time step n; d ij is the position vector between the two
main driver for the charge transfer process and was deter- particles, and q is the charge on the particle.
mined by the experimental results, while it was specific to
Ri Rj
the materials used in the simulation and experiment. Aij;n ¼ pR dnij ; R ¼ ð8Þ
Ri þ Rj
Mukherjee et al. (2016) introduced a tribo-charging model
depending on contact area and work function difference to where R is the effective radius and dnij is the normal over-
predict tribocharging in Hopper-chute flow with DEM-
lap between the particle i and the particle j.
based simulation and experiment. The charge transfer via
The charge transferred is only taken into account when
tribocharging only counts once per collision at the maxi-
DA is positive.
mum impact time step and does not take the electric field
The electric field at the center of a particle is derived by
into account. The model in Laurentie et al. (2013) also
Coulomb’s law and the principle of superposition
depends on contact area and work function difference,
(Laurentie et al., 2013),
but it includes the total electric field and is evaluated at
each time step. In addition, the experiment of granular 1 N dik
Ei ¼ Eext þ R q ð9Þ
insulating material particles was conducted to validate 4pe0 k¼1;k–i k kdik k3
the DEM numerical simulation model using in-house code
in Laurentie et al. (2013). It is worth noting that the latest where Eext represents an external electric field which should
tribo-charging model in Mukherjee et al. (2021) developed be added in virtue of the principle of superposition to the
from the model in Mukherjee et al. (2016) also includes electric field induced by each discrete charge qk , and N is
electric field. the total number of particles in the domain.
It is assumed as above that the net charge on any parti- The second term of the RHS of Eq. 9 represents the
cle is distributed homogeneously at each time step. In total E-field induced by each discrete charged particles to
Mukherjee et al. (2016), the charge transferred via tribo- include the influence of other charged particles nearby
electric charging is assumed to be a function of the differ- and is evaluated in the tribo-charging model during contact
ence between the work functions of the contacting as indicated. The electric field Eij is computed at the con-
particles, the contact area, and the cutoff distance for tact point between the colliding particles i and j, which is
3235
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

linearly interpolated from the values of Ei and Ej located at model in virtue of the principle of superposition to the elec-
each particle center. tric field due to charged particles nearby.
Although the model in Eq. 6 from Mukherjee et al. The computation time for a 90 s simulation can be quite
(2016) is good for computational efficiency as the charge prohibitive, especially for ‘‘hard” particles with high
transfer is evaluated less frequently, the varying electric Young’s Modulus and high coefficient of restitution since
field, that is not included, could be critical in the lunar dust it would be difficult to track the short contacts during
environment. Thus, to improve the simulation accuracy, one time step. To overcome that limitation, it is a common
the model in Eq. 7 from Laurentie et al. (2013) is adopted practice to adopt a lower Young’s Modulus based on the
where the electric field at the point of contact between the previous investigations of Young’s Modulus, such as in
two particles are incorporated and the charge transfer are Lommen et al. (2014, 2009). In addition, the details of
evaluated at each time step. Meanwhile, the experimental some parameters including coefficient of restitution and
results from the same paper are used to validate the simu- coefficient of friction are not stated in Laurentie et al.
lation results conducted in LIGGGHTS. (2013). Different from the spherical particles used in
LIGGGHTS simulation, Laurentie et al. (2013) used cylin-
drical particles. The height of cylinder is not specified in
2.3.1. Validation of tribo-charging model Laurentie et al. (2013) but mentioned in their preliminary
To validate the tribo-charging model in Eq. 7, a simula- work Laurentie et al. (2010). The densities are kept the
tion of a vibrating box containing 20 g of Polyamide (PA) same in this simulation as in the reference, but the differ-
and 20 g of Polycarbonate (PC) is implemented in ence in shape leads to difference mass per particle and sur-
LIGGGHTS to reproduce the results in Laurentie et al. face area.
(2013). The results from LIGGGHTS simulation are then However, theoretically, changing Young’s Modulus has
validated with both numerical simulation and experimental direct influence on the contact dynamics and results in dif-
results published in Laurentie et al. (2013). In this simula- ferent contact force according to Eq. 5. According to Eq. 7
tion, the neutral mixture is introduced in the vibrating box and Eq. 8, increasing the radius and contact area also
made of Aluminum (Al) connected to the ground and ver- affects the amount of charge transferred via tribo-
tically vibrated with an amplitude of 1 mm at a frequency charging. Larger surface area generally contributes to more
of 50 Hz during 90 s before being separated in a free-fall charge transfer. Additionally, varying the coefficient of
electrostatic separator and measured in a Faraday cage. restitution and the coefficient of friction can change the
A vibrating box is set up in LIGGGHTS, shown in number of collisions and the contact time. As a result,
Fig. 2. The particle–surface charge transfer interactions the charge transfer is affected.
are included and the wall is treated as if it is grounded. Thus, two simulations (Sim 1 and Sim 2) with different
In this numerical experiment, the vibrating box contains parameters including Young’s Modulus, coefficients of
particles made of 2 materials, Polyamide (PA) and Polycar- restitution, coefficients of friction and surface areas within
bonate (PC), with work functions of 3.9 eV and 4.2 eV, a reasonable range are implemented to reduce the potential
respectively (Laurentie et al., 2013). The system is simu- errors. Parameters used in LIGGGHTS simulation and in
lated for 90 s real world time with a time step size of 1 ls Laurentie et al. (2013) including Young’s Modulus, radius,
as 0.1 to 10 ls is suggested in Laurentie et al. (2013). The shape, etc. are summarized in Table 1. In Sim 1, ‘‘softer”
cutoff distance of tribo-charging is set to 500 nm as in particles with lower Young’s Modulus and relatively low
Laurentie et al. (2013). According to Laurentie et al. coefficient of restitution and coefficient of friction are
(2013), the charge limits for PA and PC have been set to adopted. The radius of the spherical PA and PC particles
70 pC per particle and 190 pC per particle, respectively. are the same as the radius of cylindrical PA and PC parti-
The cutoff for electrostatic force is set to 2 cm. An external cles from Laurentie et al. (2013). The same total mass is
E-field of 100000 V/m is also added to the charge transfer used in Sim 1 as in Laurentie et al. (2013). In Sim 2, real-
istic Young’s Modulus is used together with higher coeffi-
cient of restitution and coefficient of friction. To get the
same surface area as the cylinder used in Laurentie et al.
(2013), the radii of PA and PC particles are increased from
1 mm and 1.25 mm to 1.41 mm and 1.63 mm, respectively.
The number of PA and PC particles in Sim 2 is the same as
in Laurentie et al. (2013).
The outcome of the two LIGGGHTS simulations are
plotted along with the simulation and experimental results
from Laurentie et al. (2013) in Fig. 3. However, the unit for
the average charge per particle on the plot from Laurentie
et al. (2013) was corrected from nC to pC. Considering the
context and the preliminary results from Laurentie et al.
Fig. 2. Vibrating box in LIGGGHTS.
3236
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

r1 = 1, r2 = 1.25, h = 3
Not mentioned, 0.5–0.9
Laurentie et al. (2013)

n1 = 1877, n2 = 1133
PA = 2.5, PC = 2.3,

Not mentioned

Cylinder
Al = 69

PA = 0.1, PC = 0.31, Al = 1.05, PA&PC = 0.9, PA&Al = 0.9,

PA = 0.7, PC = 0.65, Al = 0.1, PA&PC = 0.9, PA&Al = 0.9,
PA = 2.5, PC = 2.3, Al = 2.5

n1 = 1877, n2 = 1133

r1 = 1.41, r2 = 1.63
PC&Al = 0.9

PC&Al = 0.9

Sphere
Sim 2

Fig. 3. Comparison of the mean charge per particle measured in

LIGGGHTS simulation and in Laurentie et al. (2013).
PA = 0.1, PC = 0.31, Al = 1.05, PA&PC = 0.21, PA&Al = 0.21,
PA = 0.7, PC = 0.65, Al = 0.1, PA&PC = 0.5, PA&Al = 0.5,

(2010), the unit used for the charge is pC instead of nC as

Comparison of parameters used in LIGGGHTS simulation and in Laurentie et al. (2013).

shown on the plot here.

In the 90 s simulation, the particles experience a signif-
PA = 0.25, PC = 0.25, Al = 0.25

icant amount of collisions. Charge transfer by tribo-

charging happens during those contacts until the charge
n1 = 4225, n2 = 2038

amount on individual particle reaches the charge limit.

r1 = 1, r2 = 1.25
PC&Al = 0.21
PC&Al = 0.5

PC has a higher work function than PA which means it

Sphere
Sim 1

is more difficult to remove an electron from PC than from

PA. As a result, PA charges positively and PC charges neg-
atively eventually. From the comparison, it is not hard to
notice that the results from the LIGGGHTS simulation
matched well with the results reported by Laurentie
et al. (2013). There are some discrepancies in the first 15
s, but when it comes closer to the end, the LIGGGHTS
simulation matched better to both numerical and experi-
mental results presented in Laurentie et al. (2013). The dis-
crepancies may result from the differences in the
properties, such as shape and mass per particle, but they
are within a reasonable range. This comparison validates
Coefficient of Restitution
Young’s Modulus(GPa)

the tribo-charging model in LIGGGHTS.

Coefficient of Friction

Dimension(mm)
No. of Particles
Properties

Shape

3. Sensitivity analyses of charged lunar dust dispersion and

velocity
Table 1

A sensitivity analysis is conducted to quantify the

effects of initial charge, tribo-charging, and E-field on dis-
3237
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

persion and velocity of lunar dust. The chemical composi-

tion of lunar dust varies across the lunar surface. One of
the most widely distributed lunar soil simulant is JSC-1
(Sternovsky et al., 2002; Colwell et al., 2007). Defined as
a industry standard simulant for education and research,
JSC-1 is a low-titanium mare-like soil with a high percent-
age of glass comparable to Apollo 14 lunar sample 14163
(Colwell et al., 2007).
In this simulation, two types of particles based on JSC-1
and lunar samples are used consisting of 50 % of each type.
The particle properties are summarized in Table 2 (Colwell
et al., 2007; Morgan et al., 2018; Brisset et al., 2018).
According to Eq. 7, the work function is the driving factor
of tribo-charging; however, it is difficult to be measured
accurately. Thus, an estimation of the work function is
adopted where the two types of particles have work func-
tions of 5 eV and 5.3 eV, respectively (Trigwell et al., 2009).
Fig. 4 introduces the setup of the sensitivity analysis,
where 100 particles will fall from rest for 1 s under lunar
gravity in each run. The starting height is assumed to be
0 m. The particles start from rest on a fictitious plane at
the height of h0 ¼ 0 m with the dimension of 1 mm by
1 mm. The particles are then released at the start of the
simulation. After the particles are released, there is no wall
in the space.
Lunar gravity 1.62 ms2 is applied downward in the Z- Fig. 4. Sensitivity analysis simulation setup.
direction as in Fig. 4. Each simulation is run for 1 s with a
time step size of 1 ls. Fig. 4b illustrates the charged parti-
cles in motion, where the color scale represents the amount rest. For each simulation, the average position is recorded
of charge on each particle. with the standard deviations of the final X, Y, Z positions
to represent the dispersion of the particles. The average
velocity in the three directions with the standard deviations
3.1. Initial charge
of final velocity is also recorded to investigate the motion
of particles. The average percentage difference between
According to Sternovsky et al. (2002), a JSC-1 particle
final charge and initial charge on each particle due to
with a radius of 50 lm can get a charge on the order of
tribo-charging is recorded as well. The standard deviations
0.016 pC - 0.16 pC from contact charging. Thus, in order
of final position and the percentage difference of 10 cases
to conduct a sensitivity analysis of the initial charge, the
with varying initial charge are summarized in Table 3,
initial charge is randomly assigned in uniform distributions
and the standard deviations of final velocity in the X-,
around these values. For each case, the seed used for the
Y-, Z-direction are shown in Table 4.
random number generator is the same to keep the simula-
tion reproducible and the results comparable.
A case with no charge, no initial velocity and no E-field 3.2. Tribo-charging
is first run as a baseline in Case 1. Tribo-charging is
included in all 10 cases except the baseline case. In all cases, To evaluate the effect of tribo-charging on the dispersion
no external E-field is applied and all particles start from of final position and final velocity of charged particles,
Case 1 to 10 are re-run without tribo-charging.
Table 5 summarizes the comparison between the stan-
Table 2
dard deviations of the cases without tribo-charging and
Particle properties used in lunar environment simulations,
their counterpart of the cases with tribo-charging. A posi-
Properties P1 P2 P1&P2
tive value represents the standard deviation for the case
Young’s Modulus(MPa) 65 65 without tribo-charging is bigger than the one for the case
Poisson’s Ratio 0.33 0.33
with tribo-charging, and vice versa. The comparison
Density(kg/m3) 2900 2900
Work Function(eV) 5.0 5.3 between the standard deviations of the final velocity in
Radius(lm) 50 50 the three directions of the cases without tribo-charging
Percentage by Number(%) 0.5 0.5 and their counterpart of the cases with tribo-charging is
Coefficient of Restitution 0.45 0.45 0.45 recorded in Table 6. Similarly, a positive value represents
Coefficient of Friction 0.84 0.84 0.84
the standard deviation of velocity for the case without
3238
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

Table 3
Table of standard deviation of final position and charge transfer varying initial charge for Case 1–10.
Case Range of Initial Charge STD of Final X (mm) STD of Final Y (mm) STD of Final Z (mm) Average Percent of Charge Change
(pC) (%)
1 0 0.3071 0.3166 0.1551 0
2 [-6.4E-5,6.4E-5] 0.3049 0.3151 0.1469 4.56
3 [-6.4E-4,6.4E-4] 0.3076 0.3159 0.1499 0.51
4 [-6.4E-3,6.4E-3] 0.4870 0.5612 0.2955 0.31
5 [-3.2E-2,3.2E-2] 1.632 2.234 1.294 0.38
6 [-6.4E-2,6.4E-2] 3.121 4.368 2.560 0.59
7 [-3.2E-1,3.2E-1] 14.10 21.13 8.282 0.39
8 [-6.4E-1,6.4E-1] 28.13 39.83 15.40 0.26
9 [-3.2,3.2] 140.4 189.5 73.65 0.32
10 [-6.4,6.4] 271.4 402.3 136.7 0.62

Table 4
Table of standard deviation of final velocity varying initial charge for Case 1–10.
Case Range of Initial Charge (pC) STD of Final Velocity in X (mm/s) STD of Final Velocity in Y (mm/s) STD of Final Velocity in Z (mm/s)
1 0 0.05520 0.05659 0.1359
2 [-6.4E-5,6.4E-5] 0.04986 0.05930 0.1290
3 [-6.4E-4,6.4E-4] 0.05645 0.07419 0.1332
4 [-6.4E-3,6.4E-3] 0.3346 0.4477 0.2994
5 [-3.2E-2,3.2E-2] 1.558 2.192 1.307
6 [-6.4E-2,6.4E-2] 3.090 4.354 2.665
7 [-3.2E-1,3.2E-1] 14.30 21.42 8.885
8 [-6.4E-1,6.4E-1] 28.39 40.45 16.63
9 [-3.2,3.2] 147.9 190.4 79.45
10 [-6.4,6.4] 288.2 415.3 153.4

Table 5
Table of standard deviation of final position without tribo-charging for Case 1–10.
Case STD of X Difference to Tribo-charging STD of Y Difference to Tribo-charging STD of Z Difference to Tribo-charging
(mm) Cases (%) (mm) Cases (%) (mm) Cases (%)
1 0.3071 0 0.3166 0 0.1550 0
2 0.3071 0.66 0.3166 0.63 0.1551 5.44
3 0.3096 0.65 0.3166 0.32 0.1581 5.33
4 0.5102 4.72 0.5226 6.77 0.3420 15.54
5 1.7742 8.84 2.064 7.60 1.462 13.00
6 3.586 14.92 3.818 12.59 2.911 13.70
7 16.91 19.91 16.45 22.15 9.751 17.74
8 32.36 15.03 32.14 19.30 18.02 16.98
9 152.6 8.69 157.3 16.99 73.68 0.046
10 308.1 13.52 314.2 21.90 136.7 21.95

tribo-charging is bigger than the one for the case with the standard deviations of final velocity in the X-, Y-,
tribo-charging, and vice versa. Z-direction with and without tribo-charging are shown in
The average positions of particles in the X-, Y-, Fig. 8.
Z-direction with error bars are plotted on semi-log scale When the initial charge is very small in the first few
as shown in Fig. 5, while the standard deviations of final cases, the standard deviations of the final positions vary
X, Y, Z positions with and without tribo-charging are plot- a little. Starting from [-6.4E-4,6.4E-4] pC, the standard
ted on a log–log scale versus the initial charge in Fig. 6. The deviations are almost in linear relationship with the initial
average velocities of particles in the X-, Y-, Z-direction charge range on log–log scale for all three components
with error bars are plotted on semi-log scale in Fig. 7, while according to Fig. 6. Thus, more charge leads to wider final

3239
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

Table 6
Table of standard deviation of final velocity without tribo-charging for Case 1–10.
Case STD of Final Difference to Tribo- STD of Final Velocity Difference to Tribo- STD of Final Velocity Difference to Tribo-
Velocity in X (mm/s) charging Cases (%) in Y (mm/s) charging Cases (%) in Z (mm/s) charging Cases (%)
1 0.05520 0 0.05659 0 0.1359 0
2 0.05526 10.82 0.05669 4.38 0.1359 5.43
3 0.06589 16.84 0.07590 2.29 0.1399 5.26
4 0.3629 8.36 0.4024 10.27 0.3722 24.41
5 1.760 12.81 2.039 6.99 1.511 15.61
6 3.573 15.63 3.803 12.66 2.992 12.25
7 17.38 21.50 16.92 21.02 10.10 13.71
8 33.12 16.66 33.42 17.38 19.20 15.43
9 162.0 9.53 164.5 13.60 81.90 3.08
10 318.6 10.55 323.4 22.13 175.8 14.60

Fig. 5. Distribution of charged particles with/without tribo-charging for Fig. 6. Dispersion of charged particles with/without tribo-charging for
case 1–10. case 1–10.

3240
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

Fig. 7. Distribution of final velocity of charged particles with/without Fig. 8. Dispersion of final velocity of charged particles with/without tribo-
tribo-charging for case 1–10. charging for case 1–10.

distribution of dust particles. Regarding the average charge Thus, more charge leads to higher dispersion and wider
change per particle, tribo-charging does not play a big role distribution of final velocity of dust particles.
in the 1 s experiment. However, the difference in standard
deviations of final position between cases with and without 3.3. E-field
tribo-charging becomes larger with increased initial charge.
The difference can be around 22% for one direction as According to the estimation of the electric field on the
shown in Table 5. The standard deviations of final velocity lunar surface in Stubbs et al. (2007) and Colwell et al.
have similar trends with the standard deviations of final (2007), a range of 1 V/m to 10 V/m is reasonable, although
position that are approximately in linear relationship with it could be as high as 1000 V/cm over the major portion of
the initial charge range on log–log scale according to Fig. 7. the sunlit area in the terminator region according to
Meanwhile, the difference in standard deviation of final Colwell et al. (2007), De and Criswell (1977). The lunar
velocity between cases with and without tribo-charging surface on the dayside typically charges positive while on
can be as large as 22% for one direction as shown in the nightside the surface usually charges negative (Stubbs
Table 6. et al., 2007). Thus, varying external E-field in the
3241
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

Z-direction is applied to conduct the sensitivity analysis. The histograms of final X, Y, Z positions of the case ini-
The standard deviations and percentage difference of 9 tially charged compared to the case without initial charge
cases are recorded in Table 7, while the standard deviations are shown in Fig. 10 and the histograms of final velocity
of final velocity in the X-, Y-, Z-direction are recorded in in the X-, Y-, Z-direction are shown in Fig. 11.
Table 8. The initial charge is randomly assigned in uniform The 3-D plot directly shows the comparison of charged
distribution of 0.064 pC - 0.064 pC for all cases from particles and uncharged particles. For the uncharged case,
Case 11 to 18. Case 6 from Tables 3 and 4 is the starting the initial and final locations are seen as two dots (in blue
point with the same initial charge as in Case 11 to 18 and and black, respectively) on the plot since all the particles
no E-field applied. are close to each other. For the case with assigned initial
According to Tables 7 and 8, E-field shows less impact charge, the initial locations are shown in blue which are
while the inter-particle electrostatic effect is dominant. at the same location as the uncharged case, while the red
dots represent the final locations of the spread out parti-
4. DEM simulations of charged lunar dust cles. According to the plot, the particles with initial charge
shown in red loft around the group of particles without ini-
4.1. Lofted particles starting from rest tial charge shown in black (where the purple arrow is point-
ing to) at the final moment.
In this section, a scenario is presented where all the afore- As shown in Fig. 10c, charged particles have higher dis-
mentioned forces are combined to produce a near realistic persion of final Z position than uncharged particles, the
simulation in the lunar environment. In the sensitivity anal- same as in the histograms of final X, Y positions shown
ysis, the effects of initial charge, tribo-charging, and E-field in Fig. 10a and Fig. 10b, respectively. Without initial
are isolated. The DEM simulations in this section are con- charge, the final X positions stay within a range of
ducted under the influence of all of them with more particles 3.9246 mm to 2.1773 mm, the final Y positions stay
for a longer simulation time to more clearly show the electro- within a range of 5.1289 mm to 2.0299 mm, and the final
static effects on charged dust particles transport. Z positions stay within a range of 3243.5 mm to
In this simulation, 500 particles are simulated. The 3236.1 mm. When initially charged, the final X positions
properties of the particles are the same as in Table 2. The are in a range of 268.62 mm to 489.81 mm, the final Y
contact model in Eq. 5 combined with Coulomb force in positions are in a range of 250.54 mm to 424.64 mm,
Eq. 3 and the triboelectric charge transfer model in Eq. 7 and the final Z positions are in a range of 3559.9 mm
is applied in the simulation of lunar dust particle transport. to 3099.1 mm. The particles are more widely distributed
The values of initial charge and E-field are picked based on with initial charge. This result is consistent with the 3-D
the sensitivity analysis. Lunar gravity 1.62 ms2 is applied plot. According to Fig. 11, charged particles have wider
downward in the Z-direction while a E-field of 4 V/m is distribution of final velocity in the X-, Y-, Z-direction than
applied in the opposite direction. In this case, the initial uncharged particles. Without initial charge, the final veloc-
charge is assigned randomly in the uniform distribution ities in the X-direction stay within a range of 1.4873 mm/
of 0.64 pC - 0.64 pC as in Case 8. The starting height is s to 0.6137 mm/s, the final velocities in the Y-direction stay
assumed to be 0 m. The simulation starts at the height of within a range of 2.0901 mm/s to 0.5399 mm/s, and the
h0 ¼ 0 m where all the particles lying or piling up on a final velocities in the Z-direction stay within a range of
plane with the dimension of 2 mm by 2 mm. After being 3241.8 mm/s to 3238.1 mm/s. When initially charged,
released from rest, the particles are falling instantly. There the final velocities in the X-direction are in a range of
is no wall in the space once the simulation started. The sim- 134.24 mm/s to 244.87 mm/s, the final velocities in the
ulation runs for 2 s real time with a time step size of 1 ls. Y-direction are in a range of 178.18 mm/s to
The 3-D plot of the initial and final locations of particles 212.11 mm/s, and the final velocities in the Z-direction
with initial charge and of particles without charge is shown are in a range of 3400.9 mm/s to 3168.8 mm/s. The
in Fig. 9.

Table 7
Table of standard deviation of final position and charge transfer varying E-field for Case 6, Case 11–18.
Case E-field (V/m) STD of Final X (mm) STD of Final Y (mm) STD of Final Z (mm) Average Percent of Charge Change (%)
6 0 3.121 4.368 2.560 0.59
11 4 3.204 4.375 2.563 0.56
12 4 3.762 4.167 2.628 0.49
13 10 3.239 4.362 2.725 0.60
14 10 3.457 4.117 2.642 0.47
15 50 3.299 4.349 2.686 0.62
16 50 3.453 4.084 2.854 0.63
17 100 3.431 4.141 2.737 0.31
18 100 3.263 4.096 2.794 0.61

3242
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

Table 8
Table of standard deviation of final velocity varying E-field for Case 6, Case 11–18.
Case E-field (V/m) STD of Final Velocity in X (mm/s) STD of Final Velocity in Y (mm/s) STD of Final Velocity in Z (mm/s)
6 0 3.090 4.354 2.665
11 4 3.151 4.429 2.646
12 4 3.801 4.188 2.762
13 10 3.250 4.355 2.825
14 10 3.480 4.119 2.759
15 50 3.294 4.364 2.925
16 50 3.451 4.101 2.964
17 100 3.359 4.162 2.934
18 100 3.257 4.098 3.037

Fig. 9. Initial and final locations of charged and uncharged particles

starting from rest.

charged particles can gain high velocity compared to the

uncharged particles.

4.2. Lofted particles in motion

In addition to the basic scenario that particles start fall-

ing from rest, the dust particles can gain initial velocity due
to exploration activities, such as when kicked up by the
wheels of the Lunar Roving Vehicle (LRV). In this section,
initial velocity is added according to Hsu and Horányi
(2012) to represent the motion of dust clouds kicked up
by the wheels of the lunar rover following ballistic
trajectories.
In this simulation, 500 particles are simulated. The
properties of the particles are the same as in Table 2. Lunar
gravity 1.62 ms2 is applied downward in the Z-direction
while a E-field of 4 V/m is applied in the opposite direction.
The initial charge is assigned randomly in the uniform dis-
tribution of 0.64 pC - 0.64 pC. Based on the calculation
in Hsu and Horányi (2012), the rover is driving in the
X-direction, while the Z-direction points upward perpen-
dicular to the X-direction. The simulation starts at the
height of h0 ¼ 0 m, where the particles are kicked up by
the lunar rover. After being assigned initial velocity, the
particles are released instantly. For initial velocity, the X Fig. 10. Histograms of final positions with/without initial charge for
particles starting from rest.
component is assigned in the normal distribution with a
mean of 2.99 m/s and a standard deviation of 0.05 m/s,

3243
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

Fig. 12. Initial and final locations of charged and uncharged particles with
initial velocity.

particles initially charged compared to the particles without

initial charge are shown in Fig. 13, and the histograms of
final velocity in the X-, Y-, Z-direction are shown in
Fig. 14.
The 3-D plot directly shows the comparison of charged
particles and uncharged particles. The initial locations are
shown in blue as a dot. For the uncharged case, particles
followed ballistic trajectories and the final locations are
seen as dust clouds in black on the plot. For the case with
assigned initial charge, while the particles have the same
Fig. 11. Histograms of final velocities with/without initial charge for initial locations as the uncharged case, their final locations
particles starting from rest.
are represented by the red dots. With initial velocity
assigned, the particles spread out at the final moment.
However, the particles with initial charge shown in red
the Y component is assigned randomly in the uniform dis- have higher dispersion than the group of particles without
tribution of 0.05 m/s - 0.05 m/s, and the Z component is initial charge shown in black according to Fig. 12.
also assigned randomly in the normal distribution with a As shown in Fig. 13, charged particles have higher dis-
mean of 2.99 m/s and a standard deviation of 0.05 m/s. persion of final position than uncharged particles, although
The simulation runs for 4 s real time with a time step size it is not as evident as in the case where initial velocity is not
of 1 ls. Most of the particles would hit the ground after involved. Without initial charge, the final X positions stay
about 3.7 s; however, to better show the dispersion of within a range of 11.7236 m to 12.2187 m, the final Y posi-
charged dust particles without interacting with the ground, tions stay within a range of 0.2211 m to 0.1927 m, and the
particles are allowed to keep falling below the starting final Z positions stay within a range of 1.4693 m to
height. 0.4969 m. When initially charged, the final X positions
The 3-D plot of the initial and final locations of charged are in a range of 11.2771 m to 12.7147 m, the final Y posi-
and uncharged particles with initial velocity is shown in tions are in a range of 0.5337 m to 0.7080 m, and the final
Fig. 12. The histograms of final X, Y, Z positions of the Z positions are in a range of 1.8217 m to 0.2118 m. The
3244
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

Fig. 13. Histograms of final positions with/without initial charge for Fig. 14. Histograms of final velocities with/without initial charge for
particles starting with initial velocity. particles starting with initial velocity.

particles are more widely distributed with initial charge. 3.2400 m/s. The difference between cases with and
This result is consistent with the 3-D plot. According to without initial charge is not significant in the Z-direction
Fig. 14, charged particles have wider distribution of final since there is lunar gravity in the Z-direction affecting the
velocity in the X-, Y-, Z-direction than uncharged particles. motion of dust particles.
Without initial charge, the final velocities in the X-direction
stay within a range of 2.9311 m/s to 3.0545 m/s, the final 5. Discussion
velocities in the Y direction stay within a range of
0.0550 m/s to 0.0479 m/s, and the final velocities in the The simulations show that the dust particles will have
Z-direction stay within a range of 3.6073 m/s to higher dispersion and wider distribution of final velocity
3.3643 m/s. When initially charged, the final velocities when they are charged due to electrostatic effects. Accord-
in the X-direction are in a range of 2.8190 m/s to ing to the sensitivity analysis, the standard deviations of
3.1786 m/s, the final velocities in the Y-direction are in a final position increase with increasing initial charge. The
range of 0.1333 m/s to 0.1769 m/s, and the final velocities standard deviations are almost in linear relationship with
in the Z-direction are in a range of 3.6977 m/s to the range of initial charge on log–log scale for all X, Y,
3245
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

Z components except when the initial charge is very small. In the DEM simulation of 500 spherical particles for 2 s,
The simulations also show that more charge leads to wider the 3-D plot and histograms, as shown in Fig. 9–11, indi-
distribution of final velocity of dust particles. The standard cate that the charged particles are hovering around the
deviations of final velocity in all three directions are uncharged particles while the uncharged particles stay close
approximately in linear relationship with the initial charge to each other falling under lunar gravity. Dust particles
range on log–log scale and are increasing with the increased spread out and are moving at different speeds after they
initial charge. are initially charged. The charged particles have wider dis-
The dust particles can spread out widely and gain high tributions of final position and of final velocity compared
velocity when they are charged. Although the average posi- to the uncharged particles that have almost same distribu-
tion and velocity do not significantly change with varying tions as they had initially. When initial velocity is added to
initial charge, the standard deviations vary significantly. represent the scenario that dust particles are kicked up by
When the range of initial charge increases from 0 to [-6.4, the wheels of lunar rover in motion, the charged particles
6.4] pC with tribo-charging turned on, the standard devia- still have higher dispersion in final position and in final
tion of X position can increase from 0.3049 mm to velocity compared to the uncharged particles, while the dif-
271.4 mm, the standard deviation of Y position can ference between charged and uncharged particles is not as
increase from 0.3151 mm to 402.3 mm, and the standard significant as in the case where particles started from rest,
deviation of Z position can increase from 0.1469 mm to according to Fig. 13 and 14. The difference in the
136.7 mm. The standard deviation of velocity in the X- Z-direction is smaller than the difference in the X- and
direction increases from 0.04986 mm/s to 288.2 mm/s, the Y-direction because of the lunar gravity in the Z-direction.
standard deviation of velocity in the Y-direction increases In the sensitivity analysis and the lunar simulations,
from 0.05659 mm/s to 415.3 mm/s, and the standard devi- charge, size and E-field work coordinately. The three
ation of velocity in the Z-direction increases from long-range forces on individual particle are gravity, the
0.1290 mm/s to 153.4 mm/s. It is worth noting that the X electric force exerted by the environment E-field, and
and Y components are in uniform distributions of a inter-particle electrostatic forces. This research uses JSC-1
1 mm range initially and almost all the particles start at particles, of which the mean radius is about 50 lm, while
the same height, thus the standard deviations of final posi- the range of initial charge is relatively small and derived
tion are considerable which can be up to 402.3 mm for one from a contact charging only experiment (Sternovsky
direction with an increase of three orders of magnitude, et al., 2002). For an individual particle in the simulation,
according to Table 5. Although all the particles start from the mass is 1.5184 lg, the radius is 50 lm, the acceleration
rest, the standard deviations of final velocity can be as high due to gravity on the surface of the Moon is approximately
as 415.3 mm/s for one direction, according to Table 6. The 1.62 ms2, the E-field is 4 V/m, and the assigned charge is
difference in standard deviation of final velocity between 0.64 pC. Thus, the gravity G ¼ m  g is about 2.46  109
cases with and without tribo-charging can be as large as N, and the force exerted by the environment E-field
22% for one direction as shown in Table 6. Although the F ei ¼ E  qi is 2.56  1012 N. When in contact with a par-
simplified Johnson-Kendall-Roberts (SJKR) model for ticle of the same amount of charge, the Coulombic force
cohesion is available in LIGGGHTS, cohesion is not k q  q
between the two particles F eij ¼ di 2 j is about
included in the simulations to isolate the electrostatic ij

effects on the charged dust particles, which is the focus of 3.68  107 N. Compared to the gravity and inter-
this study. particle electric force, the force exerted by the lunar E-
The results indicate that the E-field does not play a sig- field is not dominant.
nificant role in particle dispersion and in final velocity dis- In the real world, dust particles could get much more
tribution while the inter-particle electrostatic effect is charge in the plasma environment due to rocket plume
dominant. There is not too much change in the standard when perturbed by lunar landers. In the simulation, the ini-
deviations of final position and of final velocity with vary- tial velocity is estimated from the calculation of ballistic
ing E-field. When the E-field varies from 100 V/m to trajectory when dust particles are kicked up by the lunar
100 V/m, the standard deviation of X position stays within rover. The charged particles may have a different range
a range from 3.121 mm to 3.762 mm, the standard devia- of initial velocity due to the impingement of rocket exhaust
tion of Y position stays within a range from 4.117 mm to and other mechanical interactions on the lunar surface, and
4.375 mm, and the standard deviation of Z position stays thus the effects of electrostatics on the dispersion of
within a range from 2.560 mm to 2.854 mm, according to charged dust particles may differ from the present case.
Table 7. Similarly, the standard deviation of velocity in Spherical particles are used in the simulation to simplify
the X-direction is within a range from 3.090 mm/s to the problem, while most of the real lunar regolith is highly
3.801 mm/s, the standard deviation of velocity in the angular for which the specific surface area is nearly eight
Y-direction is within a range from 4.098 mm/s to times a sample of spheres with the same size distribution
4.429 mm/s, and the standard deviation of velocity in the (Colwell et al., 2007). The increased surface area with com-
Z-direction is within a range from 2.646 mm/s to plex shapes could lead to more charge on individual parti-
3.037 mm/s, according to Table 8. cle. According to Eq. 3, more charge initially assigned leads
3246
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

to larger electrostatic force and thus larger acceleration on slower than uncharged particles due to electrostatic effects.
individual particle. As a result, the motion of charged par- When initial velocity is involved to represent that dust par-
ticles increases and more collisions between particles could ticles are kicked up, the charged particles still have higher
happen which results in more tribo-charging. Thus, electro- dispersion in final position and wider distribution of final
static interactions could have a more significant effect on velocity compared to the uncharged particles, although
real lunar dust. the difference is not as significant as in the case where par-
The simulations start with initial charge and initial posi- ticles started from rest. This result provides a potential
tions randomly assigned in uniform distributions for the explanation for the phenomena of the approximately 30 s
sensitivity analysis. Each set of initial conditions may lead dust hovering after engine cut off following Apollo Lunar
to slightly different behavior, which explains why the stan- Module landing. Dust particles can have a high dispersion
dard deviation of X component might be different from the of position and also gain high velocity when charged by
standard deviation of Y component. Many processes might rocket plume or other mechanical interactions. The disper-
be contributing to charging of the dust particles, including sion in position and the high velocity of the charged dust
dust-dust collisions and exposure to the plasma. Both pos- particles due to electrostatic effects become more apparent
itive and negative charges are expected to be generated in with increased initial charge, and tribo-charging plays a
dust-dust collisions (Horányi and Szalay, 2015). In the pre- more important role in the dynamics of charged particles
sent study, each particle has the same probability to get a with a large amount of charge. Future investigations will
positive or negative initial charge to generalize the model. employ aspherical geometries and non-uniform charge dis-
In reality, the type of charge depends on the surroundings, tributions on particles as well as consider other factors such
such as plume charging, dayside or nightside, and tribo- as screened electrostatic force and cohesion in the dust par-
electric chargeability. If most of the particles are charged ticle model.
positively or negatively, the effect of dispersion may be
more apparent because of the repulsion between particles Declaration of Competing Interest
with the same charge.
Finally, the simulations are only run for up to 5 s. The The authors declare that they have no known competing
effect could be more noticeable in a longer time, for which financial interests or personal relationships that could have
more computational power is required. However, for appeared to influence the work reported in this paper.
longer time simulations, the accumulated errors from com-
putation should be considered, especially with particles on Acknowledgments
the micrometer scale.
This work has been supported by NASA STMD
6. Conclusion Center Innovation Fund (CIF) (Grant
#80NSSC20P0391), The Center for Lunar and Asteroid
With limited knowledge of the real lunar environment, Surface Science (CLASS) (Grant #80NSSC19M0214),
this work aims to advance our understanding of the elec- and UCF seed funding.
trostatic effects on the dynamics of charged dust particles
in the lunar environment. A DEM approach is developed References
to investigate the dynamics of charged dust particles on
the lunar surface, focusing on their inter-particle interac- Afshar-Mohajer, N., Damit, B., Wu, C.-Y., Sorloaica-Hickman, N., 2011.
tions and contact charge transfer. This model of dust par- Efficiency evaluation of an electrostatic lunar dust collector. In: 41st
International Conference on Environmental Systems, p. 5201.
ticle dynamics includes both short- and long-range Afshar-Mohajer, N., Wu, C.-Y., Sorloaica-Hickman, N., 2012. Efficiency
interactions between spherical particles. A tribo-charging determination of an electrostatic lunar dust collector by discrete
model based on instantaneous collisions between particles element method. J. Appl. Phys. 112 (2), 023305.
is adopted and validated by comparing the simulation Afshar-Mohajer, N., Wu, C.-Y., Sorloaica-Hickman, N., 2014. Electro-
static collection of tribocharged lunar dust simulants. Adv. Powder
results to existing data. In the simulations, JSC-1 simulant
Technol. 25 (6), 1800–1807.
particles with a radius of 50 lm are adopted to represent Alder, B.J., Wainwright, T.E., 1957. Phase transition for a hard sphere
the lunar dust behavior under the estimated initial condi- system. J. Chem. Phys. 27 (5), 1208–1209.
tions. Sensitivity analysis is conducted to evaluate the Alder, B.J., Wainwright, T.E., 1959. Studies in molecular dynamics. i.
effects of initial charge, tribo-charging, and E-field on the general method. J. Chem. Phys. 31 (2), 459–466.
dispersion and velocity of lunar dust. DEM simulations Alder, B.J., Wainwright, T.E., 1960. Studies in molecular dynamics. ii.
behavior of a small number of elastic spheres. J. Chem. Phys. 33 (5),
of 500 lunar dust particles are then implemented with 1439–1451.
aforementioned forces and charge transfer to more clearly Berger, K.J., Hrenya, C.M., 2017. Predicting regolith erosion during a
show the effects of electrostatics on charged particles in a lunar landing: Role of continuous size distribution. J. Aerospace Eng.
near realistic simulation in the lunar environment. Accord- 30 (5), 04017027.
Bierwisch, C., Kraft, T., Riedel, H., Moseler, M., 2009. Three-dimensional
ing to the simulations, the charged dust particles can be
discrete element models for the granular statics and dynamics of
hovering around while the uncharged particles stay close powders in cavity filling. J. Mech. Phys. Solids 57 (1), 10–31.
under the lunar gravity. The charged particles can dissipate
3247
H. Wang et al. Advances in Space Research 70 (2022) 3231–3248

Brisset, J., Colwell, J., Dove, A., Abukhalil, S., Cox, C., Mohammed, N., Lommen, S., Schott, D., Lodewijks, G., 2014. Dem speedup: Stiffness
2018. Regolith behavior under asteroid-level gravity conditions: low- effects on behavior of bulk material. Particuology 12, 107–112.
velocity impact experiments. Progress Earth Planet. Sci. 5 (1), 1–21. Marshall, J., Richard, D., Davis, S., 2011. Electrical stress and strain in
Carter, D., Hartzell, C., 2017. Extension of discrete tribocharging models lunar regolith simulants. Planet. Space Sci. 59 (14), 1744–1748.
to continuous size distributions. Phys. Rev. E 95 (1), 012901. Martin, C., Bouvard, D., Shima, S., 2003. Study of particle rearrangement
Colwell, J., Batiste, S., Horányi, M., Robertson, S., Sture, S., 2007. Lunar during powder compaction by the discrete element method. J. Mech.
surface: Dust dynamics and regolith mechanics. Rev. Geophys. 45 (2). Phys. Solids 51 (4), 667–693.
Colwell, J.E., Gulbis, A.A., Horányi, M., Robertson, S., 2005. Dust Matuttis, H.-G., Chen, J., 2014. Understanding the discrete element
transport in photoelectron layers and the formation of dust ponds on method: simulation of non-spherical particles for granular and multi-
eros. Icarus 175 (1), 159–169. body systems. John Wiley & Sons.
Cundall, P.A., Strack, O.D., 1979. A discrete numerical model for Metzger, P.T., Smith, J., Lane, J.E., 2011. Phenomenology of soil erosion
granular assemblies. Geotechnique 29 (1), 47–65. due to rocket exhaust on the moon and the mauna kea lunar test site. J.
De, B.R., Criswell, D.R., 1977. Intense localized photoelectric charging in Geophys. Res.: Planets 116 (E6).
the lunar sunset terminator region, 1. development of potentials and Morgan, P., Grott, M., Knapmeyer-Endrun, B., Golombek, M., Delage,
fields. J. Geophys. Res. 82 (7), 999–1004. P., Lognonné, P., Piqueux, S., Daubar, I., Murdoch, N., Charalam-
Forward, K.M., Lacks, D.J., Sankaran, R.M., 2009. Triboelectric bous, C., et al., 2018. A pre-landing assessment of regolith properties
charging of lunar regolith simulant. J. Geophys. Res.: Space Phys. at the insight landing site. Space Sci. Rev. 214 (6), 1–47.
114 (A10). Morris, A.B., Goldstein, D.B., Varghese, P.L., Trafton, L.M., 2016.
Hartzell, C., Zimmerman, M., Hergenrother, C., 2021. Simulations of Lunar dust transport resulting from single-and four-engine plume
electrostatic lofting and subsequent particle motion at bennu. Bull. impingement. AIAA J. 54 (4), 1339–1349.
AAS 53 (7), URL: https://baas.aas.org/pub/2021n7i104p02. Mukherjee, R., Gupta, V., Naik, S., Sarkar, S., Sharma, V., Peri, P.,
Hartzell, C.M., Scheeres, D.J., 2011. The role of cohesive forces in particle Chaudhuri, B., 2016. Effects of particle size on the triboelectrification
launching on the moon and asteroids. Planet. Space Sci. 59 (14), 1758– phenomenon in pharmaceutical excipients: Experiments and multi-
1768. scale modeling. Asian J. Pharm. Sci. 11 (5), 603–617.
Hasan, A., Alshibli, K.A., 2010. Discrete element modeling of strength Mukherjee, R., Sansare, S., Nagarajan, V., Chaudhuri, B., 2021. Discrete
properties of johnson space center (jsc-1a) lunar regolith simulant. J. element modeling (dem) based investigation of tribocharging in the
Aerospace Eng. 23 (3), 157–165. pharmaceutical powders during hopper discharge. Int. J. Pharm. 596,
Heiken, G.H., Vaniman, D.T., French, B.M., 1991. Lunar Sourcebook: a 120284.
user’s guide to the Moon. Cambridge University Press. Patel, A., Hartzell, C., 2021. A model to predict the size of regolith clumps
Hogue, M.D., Calle, C.I., Weitzman, P.S., Curry, D.R., 2008. Calculating on planetary bodies. Planet. Sci. J. 2 (5), 196.
the trajectories of triboelectrically charged particles using discrete Phillips III, J.R., Wang, H., Hillegass, A., Esparza, A., Dove, A.R.,
element modeling (dem). J. Electrostat. 66 (1–2), 32–38. Elgohary, T.A., 2021. Implementation of charged particle behavior in
Horányi, M., Szalay, J., 2015. Dust charge measurements by the lunar discrete element method (dem) simulations. In: Earth and Space 2021.
dust experiment. In: 2015 IEEE Aerospace Conference. IEEE. pp. 1–5. pp. 188–199.
Hou, X., Zhang, K., Jiang, Y., Xue, P., Cao, P., Jiang, J., 2016. Study on Radjai, F., Richefeu, V., 2009. Contact dynamics as a nonsmooth discrete
electrostatic adhesion mechanism of lunar dust based on dem. In: 2016 element method. Mech. Mater. 41 (6), 715–728.
IEEE International Conference on Robotics and Biomimetics Ramı́rez-Aragón, C., Ordieres-Meré, J., Alba-Elı́as, F., González-Marcos,
(ROBIO). IEEE. pp. 191–196. A., 2018. Comparison of cohesive models in edem and liggghts for
Hsu, H.-W., Horányi, M., 2012. Ballistic motion of dust particles in the simulating powder compaction. Materials 11 (11), 2341.
lunar roving vehicle dust trails. Am. J. Phys. 80 (5), 452–456. Sternovsky, Z., Robertson, S., Sickafoose, A., Colwell, J., Horányi, M.,
Hughes, A.L., Colwell, J.E., DeWolfe, A.W., 2008. Electrostatic dust 2002. Contact charging of lunar and martian dust simulants. Journal
transport on eros: 3-d simulations of pond formation. Icarus 195 (2), of Geophysical Research: Planets 107 (E11), 15–1–15–8.
630–648. Stubbs, T.J., Halekas, J.S., Farrell, W.M., Vondrak, R.R., 2007. Lunar
Kloss, C., Goniva, C., Hager, A., Amberger, S., Pirker, S., 2012. Models, surface charging: A global perspective using lunar prospector data.
algorithms and validation for opensource dem and cfd–dem. Progress Dust Planet. Syst. 643, 181–184.
in Computational Fluid Dynamics, An Int. J. 12 (2–3), 140–152. Taylor, L., Schmitt, H., Carrier, W., Nakagawa, M., 2005. Lunar dust
Knuth, M.A., Johnson, J., Hopkins, M., Sullivan, R., Moore, J., 2012. problem: From liability to asset. In: 1st Space Exploration Conference:
Discrete element modeling of a mars exploration rover wheel in Continuing the Voyage of Discovery. p. 2510.
granular material. J. Terrramech. 49 (1), 27–36. Trigwell, S., Captain, J.G., Arens, E.E., Quinn, J.W., Calle, C.I., 2009.
Kruggel-Emden, H., Sturm, M., Wirtz, S., Scherer, V., 2008. Selection of The use of tribocharging in the electrostatic beneficiation of lunar
an appropriate time integration scheme for the discrete element simulant. IEEE Trans. Ind. Appl. 45 (3), 1060–1067.
method (dem). Comput. Chem. Eng. 32 (10), 2263–2279. Wang, X., Blewett, D.T., Kramer, G., Barker, D., Hartzell, C., Han, D.,
Laurentie, J., Traoré, P., Dascalescu, L., 2013. Discrete element modeling Horányi, M., Garrick-Bethell, I., Hsu, H.-W., Wang, J., et al., 2021.
of triboelectric charging of insulating materials in vibrated granular Electrostatic dust transport effects on shaping the surface properties of
beds. J. Electrostat. 71 (6), 951–957. the moon and airless bodies across the solar system. Bull. Am. Astron.
Laurentie, J.-C., Traoré, P., Dragan, C., Dascalescu, L., 2010. Numerical Soc. 53 (4), 104.
modeling of triboelectric charging of granular materials in vibrated Wang, Z., Tian, D., Ma, Z., Zhao, C., Liu, Y., Liu, Y., Ding, Y., Shen, Z.,
beds. In: 2010 IEEE Industry Applications Society Annual Meeting. 2020. An electrostatic lofting model based on adhesion separation
https://doi.org/10.1109/IAS.2010.5615669. characteristics for lunar dust. Planet. Space Sci. 184, 104883.
Leps, T., Hartzell, C., 2021. High fidelity, discrete element method Wei, H., Zhao, Y., Zhang, J., Saxén, H., Yu, Y., 2017. Liggghts and edem
simulation of magnetorheological fluids using accurate particle size application on charging system of ironmaking blast furnace. Adv.
distributions in liggghts extended with mutual dipole method. Mater. Powder Technol. 28 (10), 2482–2487.
Res. Express 8 (8), 085701. Yeo, L.H., Wang, X., Deca, J., Hsu, H.-W., Horányi, M., 2021. Dynamics
Lommen, S., Lodewijks, G., Schott, D.L., 2018. Co-simulation framework of electrostatically lofted dust on airless planetary bodies. Icarus 366,
of discrete element method and multibody dynamics models. Eng. 114519.
Comput.

3248