# Analysis of impact deformation of elastic-perfectly plastic particles

Saba Saifoori, Saeid Nezamabadi, Mojtaba Ghadiri

► To cite this version:

Saba Saifoori, Saeid Nezamabadi, Mojtaba Ghadiri. Analysis of impact deformation of elastic-perfectly plastic particles. Computational Particle Mechanics, 2024, 10.1007/s40571-024-00742-x. hal-04799795

HAL Id: hal-04799795

https://hal.science/hal-04799795v1

Submitted on 23 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analysis of impact deformation of elastic-perfectly plastic particles

Saba Saifoori<fnref n="1" /> · Saeid Nezamabadi<fnref n="1" />,<fnref n="2" /> · Mojtaba Ghadiri<fnref n="1" />

Received: 11 December 2023 / Revised: 8 March 2024 / Accepted: 11 March 2024
© The Author(s) 2024

## Abstract

Material Point Method is used to study the impact deformation of elastic-perfectly plastic spherical particles. A wide range of material properties, i.e. density, Young's modulus and yield strength, are considered. The method is particularly suitable for simulating extensive deformation. The focus of the analysis is on linking the coefficient of restitution and the percentage of the incident kinetic energy dissipated by plastic deformation, $W_p/W_i \times 100$, to the material properties and impact conditions. Dimensionless groups which unify the data for the full range of material properties have been identified for this purpose. The results show that when the particle deforms extensively, $W_p/W_i \times 100$ and the equivalent plastic strain, are only dependent on the particle yield strength and the incident kinetic energy, as intuitively expected. On the other hand, when the deformation is small, Young's modulus of the particle also affects both $W_p/W_i \times 100$ and the equivalent plastic strain. Moreover, coefficient of restitution is insensitive to Young's modulus of the material. Dimensionless correlations are then suggested for prediction of the coefficient of restitution, the equivalent plastic strain and $W_p/W_i \times 100$. Finally, it is shown that the extent to which the particle flattens due to impact can be predicted using its yield strength and initial kinetic energy.

**Keywords** Grains · Powder · Particle impact · Cold spraying · Large plastic deformation · Material point method

## 1 Introduction

Particle impact is a common occurrence in numerous applications that involve handling and processing of powders. Depending on the impact details, it can have various implications for a process, e.g. it can affect the flow behaviour of powders due to kinetic energy dissipation [1] or influence the particle-particle and particle-substrate bonding mechanism, and consequently the quality of the final film in coating processes such as cold spraying and aerosol deposition [2, 3]. Thus, investigating the impact phenomenon is important for understanding and improving the efficiency of such processes.

Considering normal impact of an elastic-perfectly plastic spherical particle with another body, a fraction of the initial kinetic energy of the impact is stored in the contacting bodies as recoverable elastic strain energy. The remaining fraction is primarily dissipated by propagation of elastic

waves and, if the initial kinetic energy is sufficiently high to induce yielding, plastic deformation. In light of this, inspection of energy loss during impact has long been a common approach to study the impact phenomenon. While elastic wave propagation is inherent in any impact regardless of the impact velocity, it has been shown experimentally and analytically that energy losses due to this mechanism are less than 3–4% of the initial kinetic energy [4–7]. However, at relatively high impact velocities, plastic deformation becomes the dominant mechanism for energy dissipation. There are various experimental [8–19] and analytical [20–27] studies that have investigated the energy dissipation during elastic-plastic impact of spheres. However, experimental investigation of particle impact is precarious as the event takes place in an extremely short span of time. Additionally, the available theoretical models are developed based on simplifying assumptions and need to be validated by experiments. Therefore, numerical simulations provide a great means for analysis of the phenomena taking place throughout impact.

Finite Element Method (FEM) and Discrete Element Method (DEM) are the most popular numerical methods used to date for simulation of particle impact. Specifically, DEM has been broadly explored over the past three decades

<footnote>Saeid Nezamabadi
saeid.nezamabadi@umontpellier.fr</footnote>

<footnote n="1">¹ School of Chemical and Process Engineering, University of Leeds, Leeds LS2 9JT, UK</footnote>

<footnote n="2">² LMGC, CNRS, University of Montpellier, Montpellier, France</footnote>
for a wide range of particle-related problems. The accuracy
of any DEM simulation depends on selecting an appro-
priate contact law as input, which can accurately account
for deformation. Analytical contact laws have widely been
adopted for DEM simulations. There are numerous models
implemented in DEM that describe the normal elastic-plas-
tic contact with [28–32] and without [21, 33, 34] adhesion.
However, such models are based on simplifying assumptions
and are currently only adequate for modelling problems that
involve small strains, i.e. not accounting for extensive defor-
mation.

As for FEM, there are several studies that model contact
between an elastic–plastic sphere and a rigid surface [35–41].
However, these studies examine the aforementioned case in
a quasi-static analysis. In fact, to the best of our knowledge,
Wu et al.’s [42] work is the only study where the specific
case of contact between an elastic–plastic sphere with a rigid
surface is modelled in the dynamic framework, though there
are several FEM studies that consider impact of an elastic
sphere with an elastic–plastic substrate [22, 43, 44]. Wu
et al. [42] have carried out two-dimensional FEM simula-
tions of impact between an elastic-perfectly plastic sphere
and a rigid wall. They conduct the simulations for differ-
ent impact velocities and mechanical properties, to analyse
the rebound behaviour of the sphere. This is done by inves-
tigating the coefficient of restitution, which is the ratio of
rebound velocity to impact velocity, representing the energy
loss during impact. Their results suggest that there is a crit-
ical velocity above which the sphere undergoes finite plastic
deformation. At lower velocities, the coefficient of restitution
is only dependent on the ratio of impact velocity to the yield
velocity (Vᵢ/Vᵧ). However, at velocities higher than the crit-
ical value, the coefficient of restitution is not only dependent
on (Vᵢ/Vᵧ), but also on the ratio of the effective Young’s
modulus to yield strength (E*/Y) of the sphere. Based on
these findings, Wu et al. propose equations for coefficient
of restitution for both small deformation and finite plastic
deformation regimes.

Nevertheless, it should be noted that although FEM is
a well-established method, the problem of mesh distortion
and element entanglement during large deformations affects
the accuracy of the calculations. On the other hand, in
recent years, a new method known as the Material Point
Method (MPM) has been developed by Sulsky et al. [45]
that overcomes this drawback. MPM discretises the mate-
rial domain by a finite set of Lagrangian material points
(integration points in FEM). Additionally, the space occu-
pied by the body is discretised by an Eulerian background
mesh. These material points are then tracked during the
deformation process and each of them is assigned with a
position and carries the state variables. The algorithm used
in MPM is as follows: initially, information is mapped from
the material points to the mesh nodes. Then, the solution

to the momentum equations is calculated on the nodes. In
the end, the nodal solution is mapped back to the mate-
rial points to update their position and state variables [46].
The moving integration points naturally bypass the mesh
distortion problems typical of FEM and provide a more
convenient solution for problems concerning extensive defor-
mation. There are few instances of MPM being used to
study impact of materials that undergo plastic deformation
[47–49]. The work of Li et al. [48] is of particular interest
as it implements MPM to model impact of elastic-perfectly
plastic disks of different mechanical properties on a rigid
target. Their results show that the normalised contact law
(force normalised by the contact force at yield and displace-
ment normalised by radius of the particle) depends on and
can be determined from E*/Y and Vᵢ/Vᵧ. Moreover, when
the coefficient of restitution is expressed in terms of Vᵢ/Vγ,
three distinct zones of deformation behaviour are identified:
small deformation, full plasticity and large deformation. Li
et al. suggest that the coefficient of restitution is only depen-
dent on Vᵢ/Vγ in the first two zones and they express this
dependency by formulating their own analytical expressions.
For the third zone, their numerical results are in perfect
agreement with Wu et al. [42] equation for coefficient of
restitution.

Considering the drawbacks of the other methods and the
fact that MPM has seldom been used for simulation of par-
ticle impact, it is timely to carry out a comprehensive study
of the impact phenomenon by MPM, especially when large
deformation is concerned. Following the approach of Wu
et al. [42] and Li et al. [48], we analyse the effect of a
much wider range of material properties, especially parti-
cle density, using three-dimensional MPM simulations to
model high-velocity impact of single elastic-perfectly plas-
tic particles on a rigid target. The analysis is focussed on
the amount of energy dissipated through plastic deforma-
tion, and coefficient of restitution. Typically, the classical
updated Lagrangian formulation of MPM (ULMPM) defines
the reference configuration based on the configuration of the
previous time step. This can lead to cell-crossing instability
as the material points might not lie at an optimal position
inside the background mesh elements. Moreover, the refer-
ence configuration is updated at each time step, which makes
ULMPM relatively computationally expensive. Hence, the
current study utilises the somewhat novel total Lagrangian
formulation of MPM (TLMPM) to compensate for these
drawbacks. In TLMPM, the reference configuration is fixed
and the material points are always associated with their ini-
tial positions. This provides an efficient approach in terms of
numerical stability and computational expenses.
# 2 Methodology

## 2.1 Explicit material point method

An implicit formulation of MPM coupled with Contact Dynamics (CD) method for treatment of frictional contact can be found in [50–53]. In this work, an explicit total Lagrangian formulation of MPM [46, 54] is coupled with a new implicit contact algorithm to model the particle.

To describe a continuum body in its initial configuration, a domain $\Omega_0$ is considered in $\mathbb{R}^D$, D being the domain dimension, with an external boundary $\partial \Omega_0$. The body is subjected to prescribed displacements and forces on its independent complementary parts of the boundary, i.e. the Dirichlet boundaries, $\partial \Omega_0^u$ and the Neumann boundaries, $\partial \Omega_0^f$. The conservation of linear momentum of the body is expressed by Eq. (1) below:

$$ \nabla \cdot \Pi(X, t) + b(X, t) = \rho(X, t)a(X, t) \quad \text{in } \Omega_0, \quad (1) $$

where $\Pi(X,t)$, $b(X,t)$, $\rho(X,t)$ and $a(X,t)$ are respectively the first Piola–Kirchhoff stress tensor, the body force, the density and the acceleration of a point at position $X$ in its initial configuration at time $t$. The boundary conditions are described by:

$$ \left\{ \begin{array}{ll} u(X, t) = \hat{u}(t) & \text{on } \partial\Omega_0^u, \\ \Pi(X, t).n = f(t) & \text{on } \partial\Omega_0^f, \end{array} \right. \quad (2) $$

where $u(X,t)$ and $\hat{u}(t)$ are the displacement and the prescribed displacement fields, respectively. The terms $f(t)$ and $n$ respectively denote the prescribed load and the outward unit normal vector to $\partial \Omega_0$.

In the MPM, the body is divided into $N_p$ material points and mass is automatically conserved as each point is appointed with a fixed amount of mass. In order to solve Eq. (1), it is important to note that MPM is identical to FEM with the difference that the integration points (material points) are mobile. So, a weak form of Eq. (1), as described by Nezamabadi et al. [55] is first calculated using the relevant boundary conditions and is then discretized into Finite elements. Finally, MPM discretises the FEM integrals using a Dirac delta function and the resulting equation is shown by Eq. (3):

$$ M a_{node}(t) = f_{int}(t) + f_{ext}(t) \quad (3) $$

where $a_{node}$ is the nodal acceleration and

$$ M = \sum_{e=1}^{N_e} \sum_{p=1}^{N_p^e} m_p N_{p0} \quad \text{Lumped mass matrix,} $$

$$ f_{int}(t) = \sum_{e=1}^{N_e} \sum_{p=1}^{N_p^e} G_{p0} \Pi_p(t) V_{p0} \quad \text{Internal force vector,} $$

Sum of body forces and surface tractions, $f_S f_{ext}(t) =$

$$ \sum_{e=1}^{N_e} \sum_{p=1}^{N_p^e} N_{p0} b_p(t) + f_S(t). $$

where $N_e$ and $N_p^e$ denote the number of the background mesh elements and that of the material points for an element, respectively. Also, $m_p$ and $V_{p0}$ are respectively the mass and volume of the material point in the reference (initial) configuration. $N_{p0}$ is the shape function matrix (interpolation matrix) for a material point $p$ and its function is to relate the quantities associated with the material points to the variables associated with the nodes of the element to which the material point belongs, at the reference configuration. $G_{p0}$ is the gradient of $N_{p0}$. It should be noted that $f_S$ also includes the nodal contact forces, $f_C$, between several bodies. In order to determine the nodal velocities $v_{node}$ from the material point velocities $v_p$, a weighted squares approach is used by solving the equation below:

$$ P_{node}(t) = M v_{node}(t) = \sum_{e=1}^{N_e} \sum_{p=1}^{N_p^e} m_p N_{p0} v_p(t) \quad (4) $$

where $P_{node}$ is the nodal momentum.

Equation (4) is solved using an explicit time integration procedure. Details of the algorithm used to update the material points position and state variables and the Contact Dynamics method to calculate the nodal contact forces can be found in detail in [56] and are not repeated here. An in-house MPM code is used for the current study.

## 2.2 Verification

In order to verify the method, the case of small deformation for impact of an elastic spherical particle is investigated through comparison with the analytical contact model provided by Johnson [20]. This is done by simulating normal impact of an elastic sphere at 1 m/s with a rigid wall placed 1.5 mm below it. The time step is set as 1 ns for the simulation. The sphere is discretised into 74,227 material points and 132,651 eight-node elements along with the cubic-spline shape function are used. The properties of the sphere are as follows: Radius, $R = 12.5$ mm, Young's modulus, $E = 4.9$ MPa, Poisson's ratio, $\nu = 0.25$ and density, $\rho = 1404$ kg/m$^3$.

According to Johnson [20], Hertz equation can be used as the contact law for impact of an elastic sphere undergoing small deformation, as shown by Eq. (5):

$$ F = \frac{4}{3} R^{1/2} E^* \delta^{3/2} \quad (5) $$

where $F$ is the contact force, and $\delta$ is the displacement of the centre of the particle. $E^*$ is the effective Young's modulus
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Variation of contact force and displacement with time for impact of a sphere undergoing small elastic deformation; comparison between analytical [20] (dashed lines) and MPM results (discrete symbols)</figcaption></figure>

of the sphere defined as $E^* = E/(1-\nu^2)$. The relationship between $\delta$ and time $t$ is given by:

$$t = \frac{\delta^*}{V_i} \int \frac{d(\delta/\delta^*)}{[1 - (\delta/\delta^*)^{5/2}]^{1/2}} \quad (6)$$

where:

$$\delta^* = \left( \frac{15mV_i^2}{16R^{1/2}E^*} \right)^{2/5} \quad (7)$$

in which $V_i$ is the impact velocity, $\delta^*$ is the maximum displacement of the centre of the particle during impact and $m$ is the mass of the particle. Deresiewicz [57] has evaluated Eq. (6) numerically, providing the values for $\delta/\delta^*$ and $2t/T$, with $T$ being the total contact time given by Johnson [20] according to Eq. (8):

$$T = 2.87 \left( \frac{m^2}{RE^{*2}V_i} \right)^{1/5} \quad (8)$$

By calculating the maximum contact force, $F^*$, from Eqs. (5) and (7), the change in contact force and displacement with time can be obtained in a dimensionless form. This is shown in Fig. 1 with the MPM results overlaid for comparison. It should be noted that the numerical total contact time is considered from the first instance the particle makes contact with the wall, to the time where contact is last detected. Considering Fig. 1, there is a good agreement between the analytical and the numerical results obtained from the MPM simulation, confirming the accuracy of the method.

## 2.3 Case studies and simulation setup

As mentioned before, the current work aims to study the impact deformation behaviour of elastic-perfectly plastic particles having a wide range of material properties. Thus, four different densities of 1000, 2000, 4000 and 8000 kg/m³ are considered. As shown by Wu et al. [42] and Li et al. [48], the ratio of the effective Young's modulus to yield strength of the material ($E^*/Y$) plays an important role in the deformation behaviour of the particle. In the current work, three different values of 1, 10 and 100 GPa for Young's modulus, and eight different values of 20, 40, 80, 160, 320, 640, 1280 and 2560 for the ratio of Young's modulus to yield strength ($E/Y$) are investigated. This is done for the four different values of density, leading to a total of ninety six cases being studied. The values of yield strength corresponding to each case are shown in Table 1.

For all of the cases, particle radius is set as 250 µm and Poisson's ratio as 0.35. For each case, the particle is normally impacting a wall placed 30 µm below it with an impact velocity of 50 m/s. After conducting a sensitivity analysis, a time step of 10 ps is selected and the particle is discretised into 74,227 material points. The mesh size is adjusted in a way that the ratio of the element size to the distance between the material points in each dimension is 1.5 for cases with a Young's modulus of 100 GPa and 1.05 for cases with a Young's modulus of 1 and 10 GPa. Damping is not applied.
<table><caption>**Table 1** The values of yield strength, *Y*, corresponding to each of the studied cases</caption><thead><tr><th rowspan="2"><i>E</i>/<i>Y</i> (-)</th><th colspan="3"><i>E</i> (GPa)</th></tr><tr><th>1</th><th>10</th><th>100</th></tr><tr><th></th><th colspan="3"><i>Y</i> (MPa)</th></tr></thead><tbody><tr><td>20</td><td>50</td><td>500</td><td>5000</td></tr><tr><td>40</td><td>25</td><td>250</td><td>2500</td></tr><tr><td>80</td><td>12.50</td><td>125</td><td>1250</td></tr><tr><td>160</td><td>6.25</td><td>62.50</td><td>625</td></tr><tr><td>320</td><td>3.13</td><td>31.25</td><td>312.50</td></tr><tr><td>640</td><td>1.56</td><td>15.63</td><td>156.25</td></tr><tr><td>1280</td><td>0.78</td><td>7.81</td><td>78.13</td></tr><tr><td>2560</td><td>0.39</td><td>3.91</td><td>39.06</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Visualisation of the modelled particle after rebound for different cases: **a** $\rho = 8000 \text{ kg/m}^3$, $E = 100 \text{ GPa}$, $E/Y = 160$; **b** $\rho = 8000 \text{ kg/m}^3$, $E = 1 \text{ GPa}$, $E/Y = 40$; **c** $\rho = 2000 \text{ kg/m}^3$, $E = 1 \text{ GPa}$, $E/Y = 320$ and **d** $\rho = 1000 \text{ kg/m}^3$, $E = 1 \text{ GPa}$, $E/Y = 1280$. The orientation is chosen randomly to provide the best view of the extent of deformation</figcaption></figure>

# 3 Results and discussion

## 3.1 Extent of deformation

In order to display the diversity in the modelled extents of deformation, visualisations of the particle after rebound are shown in Fig. 2 for four different cases, as an example. Considering Fig. 2d, it is clear that the method allows for modelling very large deformation. However, it should be noted that for a number of the cases with a very low yield stress, as highlighted in the Appendix, the deformation is so extreme that no results are obtained.

## 3.2 Coefficient of restitution

The coefficient of restitution, $e$, is calculated using Eq. (9):

$$e = \frac{V_r}{V_i} \qquad (9)$$

where $V_r$ and $V_i$ are the absolute values of the rebound and impact velocity of the particle in the impact direction, respectively. In order to find $V_r$, the velocity of the particle in the impact direction is first calculated for each time step using the following equation:

$$\sum_{p=1}^{N_p} m_p v_p^z / \sum_{p=1}^{N_p} m_p \qquad (10)$$

where $v_p^z$ is the velocity of material point $p$ in the impact direction (z axis here). Subsequently, the first instant at which the calculated velocity reaches a constant value after contact is lost is considered the rebound instant and the corresponding velocity is taken as the rebound velocity.

The coefficient of restitution obtained from Eq. (9) is plotted as a function of $V_i/V_y$ for all the cases and presented in Fig. 3a, where $V_y$ is the velocity required for the onset of plastic deformation, given by Johnson [20], as shown in Eq. (11). It should be noted that due to the high number of the case studies, no legends are displayed on any of the forthcoming graphs, and the reader is referred to Table 2 for the designation of the symbols.

$$V_y = \left( \frac{26Y}{\rho} \right)^{1/2} \left( \frac{Y}{E^*} \right)^2 \qquad (11)$$

Considering Fig. 3a, when the coefficient of restitution is plotted as a function of $V_i/V_y$, a family of curves start to form where all the data points corresponding to a certain value of $E/Y$ group together (in line with the findings of Wu et al. [42]). As a general trend, there is decrease in the coefficient of restitution with increase in $V_i/V_y$. Since $V_i$ is the same for all the studied cases in the current work and the change in $V_i/V_y$ is simply due to the change in $V_y$, this trend is expected; decrease in $Y$ and increase in $\rho$, i.e. decrease in $V_y$ translates to an earlier onset for plastic deformation, which in turn leads to a diminish in the recovered elastic energy. The trend is rapid for smaller values of $E/Y$ (which generally correspond to small deformation) and slows down with further increase in $E/Y$. It should be noted that for the cases where the coefficient of restitution is less than 0.1, the trend becomes less orderly due to numerical errors, as the extreme deformation results in a lower number of material points per mesh cell.

According to Wu et al. [42], when $(V_i/V_y)/(E^*/Y)^2 \ge 0.008$, the particle is considered to have undergone large
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Coefficient of restitution, *e*, as a function of **a** $V_i/V_y$ and **b** dimensionless group $(V_i/V_y)/(E^*/Y)$. The dashed line is the model of Wu et al. [42], shown by Eq. (12). For designation of the symbols, see Table 2. Note that the data points in each data set correspond to different densities</figcaption></figure>

deformation and the coefficient of restitution can be determined from Eq. (12).

$$e = 0.62 \left( \frac{V_i / V_y}{E^* / Y} \right)^{-1/2} \qquad (12)$$

Consequently, the coefficient of restitution for the currently studied cases that fit this criterion (all cases except for three) is plotted as a function of $(V_i/V_y)/(E^*/Y)$ and shown in Fig. 3b, along with Wu et al. model superimposed for comparison. As seen in Fig. 3b, even though displaying *e* based on $(V_i/V_y)/(E^*/Y)$ unifies the data to a certain extent, Wu

<table><caption>**Table 2** Symbol reference for all the graphs. NA denotes "not applicable"</caption><thead><tr><th colspan="4">Legend</th></tr><tr><th colspan="4">E (GPa)</th></tr><tr><th>E/Y (-)</th><th>1</th><th>10</th><th>100</th></tr></thead><tbody><tr><td>20</td><td></td><td></td><td></td></tr><tr><td>40</td><td></td><td></td><td></td></tr><tr><td>80</td><td></td><td></td><td></td></tr><tr><td>160</td><td></td><td></td><td></td></tr><tr><td>320</td><td></td><td></td><td></td></tr><tr><td>640</td><td></td><td></td><td></td></tr><tr><td>1280</td><td></td><td></td><td></td></tr><tr><td>2560</td><td>NA</td><td></td><td></td></tr></tbody></table>

et al.'s model (dashed line) appears to deviate from the MPM simulation results for extensive deformations. This is due to the fact that only six different cases are considered in Wu et al. study (two and six different values for Young's modulus and yield strength, respectively and one value for density), whereas the current work studies a much wider range of material properties. Moreover, the criterion suggested by Wu et al. to mark the boundary between small and large deformation is not deemed reasonable based on the results of the current work: there are multiple cases that fit the criterion yet their plastic strain in the impact direction is no more than 1% of the original particle diameter.

Whilst using the dimensionless group $(V_i / V_y)/(E^*/Y)$ does not fully unify the data for such a wide range of material properties, expressing the coefficient of restitution by the group $[1-(2Y\epsilon_p/\rho V_i^2)]^{0.5}$ leads to better unification for the full range, as shown in Fig. 4. Here, $\epsilon_p$ is the equivalent plastic strain at the instant of rebound, calculated from Eq. (13) where $\epsilon^p_{ij}$ is the deviatoric plastic strain.

$$\epsilon_p = \sqrt{\frac{2}{3} \epsilon_{ij}^p \epsilon_{ij}^p} \qquad (13)$$

With regards to Fig. 4, coefficient of restitution correlates well with the group $2Y\epsilon_p/\rho V_i^2$ which describes the ratio of plastic work resulting from impact to the initial kinetic energy of an elastic-perfectly plastic particle. Moreover, presenting the data in this manner reveals that coefficient of restitution is not sensitive to Young's modulus for elastic-perfectly plastic particles and the resultant decrease of *e* from unity is uniquely due to plastic deformation. Interestingly, the fitted line intercepts the x axis at about $[1-(2Y\epsilon_p/\rho V_i^2)]^{0.5} = 0.2$, suggesting that the particle will not rebound if more than 96% of the initial kinetic energy is dissipated by plastic deformation.
<figure><img src="image_6.png" /><figcaption>Fig. 4 Coefficient of restitution, $e$, as a function of the dimensionless group $[1-(2Y\epsilon_p/\rho V_i^2)]^{0.5}$. For designation of the symbols, see Table 2</figcaption></figure>

## 3.3 Plastic deformation

In order to probe the link between material properties and the plastic deformation induced by particle impact, the equivalent plastic strain, $\epsilon_p$, is plotted against $\rho V_i^2/Y$, as displayed in Fig. 5a. Considering the graph, $\epsilon_p$ increases with increase in $\rho V_i^2/Y$ which is expected as the increase in the latter is due to the decrease in the yield strength of the material or increase in the incident kinetic energy. Moreover, the data points remarkably unify and follow a common trajectory up to a certain point, after which they start to deviate from the common course. The unification is observed for data points corresponding to cases with relatively higher values of $\rho V_i^2/Y$. This suggests that for high-energy impact of highly deformable materials (large deformation), the extent to which an elastic-perfectly plastic particle undergoes plastic strain is only affected by the yield strength and the incident kinetic energy of the particle, with no sensitivity to Young's modulus. Interestingly, when $\epsilon_p$ is expressed by the group $(\rho V_i^2 / Y)(E/Y)^{0.2}$, as shown in Fig. 5b, the data points corresponding to cases with relatively lower values of $\rho V_i^2 / Y$ unify only up to a certain value for the dimensionless group (approximately 0.1). This implies that for impact of elastic-perfectly plastic particles, if the yield strength of the material is relatively high or the incident kinetic energy is relatively low (small deformation), the plastic strain caused by impact is not only dependent on yield strength and the initial kinetic energy, but also on Young's modulus of the material.

Considering the above discussion, a threshold value of 0.1 is chosen for the dimensionless group $(\rho V_i^2 / Y)(E/Y)^{0.2}$ to distinguish between these two regimes of behaviour: for cases with $(\rho V_i^2 / Y)(E/Y)^{0.2} < 0.1$, the equivalent plastic

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5 Equivalent plastic strain, $\epsilon_p$, as a function of a $\rho V_i^2/Y$ and b dimensionless group $(\rho V_i^2 / Y)(E/Y)^{0.2}$ for different cases. For designation of the symbols, see Table 2</figcaption></figure>

strain is sensitive to Young's modulus and can be obtained from the equation of the fitted line in Fig. 6a. As for cases with $(\rho V_i^2 / Y)(E/Y)^{0.2} > 0.1$, plastic strain becomes independent of Young's modulus, and the fitted line in Fig. 6b can be used for estimating the equivalent plastic strain. A summary is shown in Eq. (14).

$$ \epsilon_p = \begin{cases} 0.32 \left[ \left( \frac{\rho V_i^2}{Y} \right) \left( \frac{E}{Y} \right)^{0.2} \right]^{1.25} & \text{if } \left( \frac{\rho V_i^2}{Y} \right) \left( \frac{E}{Y} \right)^{0.2} < 0.1 \\ 0.47 \left( \frac{\rho V_i^2}{Y} \right)^{1.03} & \text{if } \left( \frac{\rho V_i^2}{Y} \right) \left( \frac{E}{Y} \right)^{0.2} > 0.1 \end{cases} \quad (14) $$
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 6** Equivalent plastic strain, $\epsilon_p$, as a function of **a** dimensionless group $(\rho V_i^2/Y)(E/Y)^{0.2}$ for cases with $(\rho V_i^2/Y)(E/Y)^{0.2} < 0.1$ and **b** $\rho V_i^2/Y$ for cases with $(\rho V_i^2/Y)(E/Y)^{0.2} > 0.1$. For designation of the symbols, see Table 2</figcaption></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 7** Percentage of the initial kinetic energy dissipated by plastic deformation, $W_p/W_i \times 100$, as a function of **a** $V_i/V_y$ and **b** dimensionless group $(E/Y)^{1.5}[\rho(V_i^2-V_y^2)/Y]$ for different cases. For designation of the symbols, see Table 2</figcaption></figure>

So as to further inspect the effect of material properties on plastic deformation of the impacting particle, the percentage of the initial kinetic energy dissipated by plastic deformation, $W_p/W_i \times 100$, is plotted as a function of $V_i/V_y$ and shown in Fig. 7a. $W_p$ is the plastic work/energy, calculated, in simple terms, by deducting the recovered elastic energy from the total deformation energy for each time step. It should be noted that as damping is not considered in the system, after rebound, the recovered elastic energy is not dissipated. This results in fluctuations in the value of the elastic energy and

consequently $W_p$, especially for soft particles. Thus, the values shown in Fig. 7a for $W_p/W_i \times 100$ are average values, taken after rebound.

With regards to Fig. 7a, similar to what was observed for $e$, when $W_p/W_i \times 100$ is expressed by $V_i/V_y$, a family of curves start to form based on $E/Y$, but the unification of data is poor. In general, there is increase in the percentage of the energy dissipated by plastic deformation with
increase in $V_i/V_y$ due to the decrease in $V_y$, i.e., the material getting more deformable. This trend is rapid and slows down as $V_i/V_y$ is increased further. For the purpose of data unification, $W_p/W_i \times 100$ is subsequently plotted against the dimensionless group $(E/Y)^{1.5}[\rho(V_i^2-V_y^2)/Y]$, as shown in Fig. 7b. In doing so, two regimes of behaviour can be recognised: for cases corresponding to small deformation, the percentage of the energy dissipated by plastic deformation is significantly affected by the combination of material properties expressed by the dimensionless group, where a rapid increase in $W_p/W_i \times 100$ is observed with increase in $(E/Y)^{1.5}[\rho(V_i^2-V_y^2)/Y]$. For this regime, a somewhat good unification is achieved up to the value of 40 for the dimensionless group, as shown in Fig. 8a. However, for the cases that undergo large deformation, $W_p/W_i \times 100$ is less sensitive to $(E/Y)^{1.5}[\rho(V_i^2-V_y^2)/Y]$, suggesting that for extensive deformation, the amount of the induced plastic work as percentage of the initial kinetic energy cannot be fully described by this combination of material properties. It was previously observed that for the cases undergoing large deformation, the equivalent plastic strain is independent of Young's modulus. Thus, as shown in Fig. 8b, $W_p/W_i \times 100$ is alternatively expressed by $\rho(V_i^2-V_y^2)/Y$, for the cases which fit the criterion $(E/Y)^{1.5}[\rho(V_i^2-V_y^2)/Y] > 40$, i.e. the data points that deviate from the common trajectory in Fig. 7b. Noting Fig. 8b, a better unification is achieved when the term $(E/Y)$ is removed from the dimensionless group. This implies that for highly deformable materials, Young's modulus has no influence on the percentage of the incident kinetic energy dissipated by plastic deformation of the particle during impact. However, it should be noted that the dependency of $W_p/W_i \times 100$ on $\rho(V_i^2-V_y^2)/Y$ for large deformation is still much weaker compared to that of $W_p/W_i \times 100$ on $(E/Y)^{1.5}[\rho(V_i^2-V_y^2)/Y]$ for small deformation. This suggests that for extensively deformable materials (or very high-energy impacts), the deformation of the particle is potentially dominated by the rate of plastic flow. Nevertheless, the equations of the fitted lines in Fig. 8a and b can be used to get an approximate prediction for the percentage of the initial kinetic energy loss due to plastic deformation based on the material properties, as summarised in Eq. (15):

$$ \frac{W_p}{W_i} \times 100 = \begin{cases} 10.40 \ln \left[ \left( \frac{E}{Y} \right)^{1.5} \frac{\rho(V_i^2 - V_y^2)}{Y} \right] + 46.60 & \text{if } \left( \frac{E}{Y} \right)^{1.5} \left[ \frac{\rho(V_i^2 - V_y^2)}{Y} \right] < 40 \\ 2.82 \ln \left[ \frac{\rho(V_i^2 - V_y^2)}{Y} \right] + 93.52 & \text{if } \left( \frac{E}{Y} \right)^{1.5} \left[ \frac{\rho(V_i^2 - V_y^2)}{Y} \right] > 40 \end{cases} \quad (15) $$

Another perspective from which plastic deformation can be studied is via investigating the extent to which the particle flattens due to impact. To this end, $H_{deformed}/H_0$, which describes the “flattening extent” of the particle, is plotted as a function of $\rho V_i^2 / Y$, as shown in Fig. 9. $H_{deformed}$ and $H_0$

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 8 Percentage of the initial kinetic energy dissipated by plastic deformation, $W_p/W_i \times 100$, as a function of a dimensionless group $(E/Y)^{1.5}[\rho(V_i^2-V_y^2)/Y]$ for cases with $(E/Y)^{1.5}[\rho(V_i^2-V_y^2)/Y] < 40$, and b $\rho(V_i^2-V_y^2)/Y$ for cases with $(E/Y)^{1.5}[\rho(V_i^2-V_y^2)/Y] > 40$. For designation of the symbols, see Table 2</figcaption></figure>

denote the length of the imaginary centre line that connects the top and bottom of the particle parallel to the impact direction, before and after impact, respectively (refer to Fig. 9). Consequently, decrease in $H_{deformed}/H_0$ implies an increase in the flattening extent. It should be mentioned that $H_{deformed}$ is measured after rebound as an average value, due to the fluctuations caused by the elastic waves in the particle.

Considering Fig. 9, a master curve is obtained when $H_{deformed}/H_0$ is expressed by $\rho V_i^2 / Y$. Also, $H_{deformed}/H_0$ is intuitively decreasing with increase in $\rho V_i^2 / Y$, i.e.
<figure><img src="image_15.png" /><figcaption>**Fig. 9** $H_{deformed}/H_0$, as a function of $\rho V_i^2/Y$ for different cases. For designation of the symbols, see Table 2</figcaption></figure>

increase in the incident kinetic energy or decrease in the yield strength of the particle. The trend is initially slow when the material is less deformable, or the impact energy is low. However, the flattening extent increases at a faster rate for more deformable materials, or high-energy impacts. Figure 9 suggests that based on the value of yield strength and the incident kinetic energy, the flattening extent of the particle can be predicted. This is beneficial for high-velocity impact deposition techniques like cold spraying, where the particles experience extensive deformation, i.e. depending on the density and yield strength of the material, the impact velocity required to achieve a desired extent of flattening for the particles can be decided.

## 4 Conclusions

A numerical study of the impact phenomenon for elastic-perfectly plastic spherical particles has been carried out using MPM. The method is first verified through comparison with

the analytical contact model of Johnson [20] for impact of an elastic sphere with a rigid wall. Subsequently, MPM is used to study the impact of an elastic-perfectly plastic particle with a rigid wall, covering a wide range of material properties and extents of deformation. The results show that for elastic-perfectly plastic particles with the studied range of material properties, the coefficient of restitution is not affected by Young's modulus of the particle, but rather by the plastic flow. Additionally, it is found that for large deformation, the percentage of the initial kinetic energy loss due to plastic deformation, as well as the equivalent plastic strain are only affected by the particle yield strength and the impact energy. However, in the case of small deformation, the amount of the stored elastic energy also plays a role in the initial energy dissipation due to plastic deformation and the equivalent plastic strain. It is also shown that the extent to which the particle flattens due to impact correlates well with the incident kinetic energy and yield strength of the material. In the end, MPM proves to be a valuable tool for study of large deformation, compensating for the shortcomings of the other available methods.

**Acknowledgements** The support of Engineering and Physical Sciences Research Council (EPSRC) under grant EP/S029036/1 and its Principal Investigator, Professor Steven Milne is gratefully acknowledged. This work was undertaken on ARC4, part of the High Performance Computing facilities at the University of Leeds, UK. The authors also acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-20-CE08-0011 (PaMaCo project).

## Declarations

**Conflict of interest** On behalf of all authors, the corresponding author states that there is no conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## Appendix

See Table 3.
<table><caption>**Table 3** The studied cases are colour coded based on the simulation success: red signifies the cases for which no result is obtained and green signifies the cases for which the simulation is successful</caption><thead><tr><th rowspan="3">E/Y (-)</th><th colspan="3">1000</th><th colspan="3">2000</th><th colspan="3">4000</th><th colspan="3">8000</th></tr><tr><th colspan="9">E (GPa)</th></tr><tr><th>1</th><th>10</th><th>100</th><th>1</th><th>10</th><th>100</th><th>1</th><th>10</th><th>100</th><th>1</th><th>10</th><th>100</th></tr></thead><tbody><tr><td>20</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>40</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>80</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>160</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>320</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>640</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1280</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>2560</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

## References

1. Wu C-Y, Dihoru L, Cocks ACF (2003) The flow of powder into simple and stepped dies. Powder Technol 134:24–39. https://doi.org/10.1016/S0032-5910(03)00130-X

2. Assadi H, Kreye H, Gärtner F, Klassen T (2016) Cold spraying—a materials perspective. Acta Mater 116:382–407. https://doi.org/10.1016/j.actamat.2016.06.034

3. Akedo J (2020) Room temperature impact consolidation and application to ceramic coatings: aerosol deposition method. J Ceram Soc Japan 128:101–116. https://doi.org/10.2109/jcersj2.19196

4. Tillett JPA (1954) A study of the impact of spheres on plates. Proc Phys Soc Sect B 67:677–688. https://doi.org/10.1088/0370-1301/67/9/304

5. Hunter SC (1957) Energy absorbed by elastic waves during impact. J Mech Phys Solids 5:162–171. https://doi.org/10.1016/0022-5096(57)90002-9

6. Hutchings IM (1979) Energy absorbed by elastic waves during plastic impact. J Phys D Appl Phys 12:1819–1824. https://doi.org/10.1088/0022-3727/12/11/010

7. Reed J (1985) Energy losses due to elastic wave propagation during an elastic impact. J Phys D Appl Phys 18:2329–2337. https://doi.org/10.1088/0022-3727/18/12/004

8. Tabor D (1948) A simple theory of static and dynamic hardness. Proc R Soc London Ser A Math Phys Sci 192:247–274. https://doi.org/10.1098/rspa.1948.0008

9. Goldsmith W, Lyman PT (1960) The penetration of hard-steel spheres into plane metal surfaces. J Appl Mech 27:717–725. https://doi.org/10.1115/1.3644088

10. Jebelisiniak F, Boettcher R, van Wachem B, Mueller P (2019) Impact of dominant elastic to elastic-plastic millimeter-sized metal spheres with glass plates. Powder Technol 356:208–221. https://doi.org/10.1016/j.powtec.2019.07.105

11. Melo KRB, de Pádua TF, Lopes GC (2021) A coefficient of restitution model for particle-surface collision of particles with a wide range of mechanical characteristics. Adv Powder Technol 32:4723–4733. https://doi.org/10.1016/j.apt.2021.10.023

12. Dahneke B (1975) Further measurements of the bouncing of small latex spheres. J Colloid Interface Sci 51:58–65. https://doi.org/10.1016/0021-9797(75)90083-1

13. Labous L, Rosato AD, Dave RN (1997) Measurements of collisional properties of spheres using high-speed video analysis. Phys Rev E 56:5717–5725. https://doi.org/10.1103/PhysRevE.56.5717

14. Kharaz AH, Gorham DA, Salman AD (1999) Accurate measurement of particle impact parameters. Meas Sci Technol 10:31–35. https://doi.org/10.1088/0957-0233/10/1/009

15. Weir G, Tallon S (2005) The coefficient of restitution for normal incident, low velocity particle impacts. Chem Eng Sci 60:3637–3647. https://doi.org/10.1016/j.ces.2005.01.040

16. Stevens AB, Hrenya CM (2005) Comparison of soft-sphere models to measurements of collision properties during normal impacts. Powder Technol 154:99–109. https://doi.org/10.1016/j.powtec.2005.04.033

17. King H, White R, Maxwell I, Menon N (2011) Inelastic impact of a sphere on a massive plane: nonmonotonic velocity-dependence of the restitution coefficient. EPL Europhysics Lett 93:14002. https://doi.org/10.1209/0295-5075/93/14002

18. Marinack MC, Musgrave RE, Higgs CF (2013) Experimental investigations on the coefficient of restitution of single particles. Tribol Trans 56:572–580. https://doi.org/10.1080/10402004.2012.748233

19. Patil D, Fred Higgs C (2018) Experimental investigations on the coefficient of restitution for sphere-thin plate elastoplastic impact. J Tribol 140:1–13. https://doi.org/10.1115/1.4037212

20. Johnson KL (1985) Contact Mechanics. Cambridge University Press, Cambridge

21. Thornton C (1997) Coefficient of restitution for collinear collisions of elastic-perfectly plastic spheres. J Appl Mech 64:383–386. https://doi.org/10.1115/1.2787319

22. Li L-Y, Wu C-Y, Thornton C (2001) A theoretical model for the contact of elastoplastic bodies. Proc Inst Mech Eng Part C J Mech Eng Sci 216:421–431. https://doi.org/10.1243/0954406021525214

23. Jackson RL, Green I, Marghitu DB (2010) Predicting the coefficient of restitution of impacting elastic-perfectly plastic spheres. Nonlinear Dyn 60:217–229. https://doi.org/10.1007/s11071-009-9591-z

24. Krijt S, Güttler C, Heißelmann D et al (2013) Energy dissipation in head-on collisions of spheres. J Phys D Appl Phys 46:435303. https://doi.org/10.1088/0022-3727/46/43/435303

25. Zhou Y (2013) Modeling of softsphere normal collisions with characteristic of coefficient of restitution dependent on impact velocity. Theor Appl Mech Lett 3:021003. https://doi.org/10.1063/2.1302103

26. Xie J, Dong M, Li S (2016) Dynamic impact model of plastic deformation between micro-particles and flat surfaces without adhesion. Aerosol Sci Technol 50:321–330. https://doi.org/10.1080/02786826.2016.1150957

27. Green I (2022) The prediction of the coefficient of restitution between impacting spheres and finite thickness plates undergoing elastoplastic deformations and wave propagation. Nonlinear Dyn 109:2443–2458. https://doi.org/10.1007/s11071-022-07522-3

28. Thornton C, Ning Z (1998) A theoretical model for the stick/bounce behaviour of adhesive, elastic-plastic spheres. Powder Technol 99:154–162. https://doi.org/10.1016/S0032-5910(98)00099-0

29. Tomas J (2007) Adhesion of ultrafine particles-A micromechanical approach. Chem Eng Sci. https://doi.org/10.1016/j.ces.2006.12.055
30. Luding S (2005) Shear flow modeling of cohesive and frictional fine powder. Powder Technol 158:45–50. https://doi.org/10.1016/j.powtec.2005.04.018

31. Luding S (2008) Cohesive, frictional powders: contact models for tension. Granul Matter. https://doi.org/10.1007/s10035-008-0099-x

32. Pasha M, Dogbe S, Hare C et al (2014) A linear model of elasto-plastic and adhesive contact deformation. Granul Matter 16:151–162. https://doi.org/10.1007/s10035-013-0476-y

33. Vu-Quoc L, Zhang X (1999) An elastoplastic contact force–displacement model in the normal direction: displacement–driven version. Proc R Soc London Ser A Math Phys Eng Sci 455:4013–4044. https://doi.org/10.1098/rspa.1999.0488

34. Walton OR, Braun RL (1986) Viscosity, granular-temperature, and stress calculations for shearing assemblies of inelastic, frictional disks. J Rheol 30:949–980. https://doi.org/10.1122/1.549893

35. Kogut L, Etsion I (2002) Elastic-plastic contact analysis of a sphere and a rigid flat. J Appl Mech 69:657–662. https://doi.org/10.1115/1.1490373

36. Quicksall JJ, Jackson RL, Green I (2004) Elasto-plastic hemispherical contact models for various mechanical properties. Proc Inst Mech Eng Part J J Eng Tribol 218:313–322. https://doi.org/10.1243/1350650041762604

37. Jackson RL, Green I (2005) A finite element study of elasto-plastic hemispherical contact against a rigid flat. J Tribol 127:343–354. https://doi.org/10.1115/1.1866166

38. Lin LP, Lin JF (2006) A new method for elastic-plastic contact analysis of a deformable sphere and a rigid flat. J Tribol 128:221–229. https://doi.org/10.1115/1.2164469

39. Brizmer V, Kligerman Y, Etsion I (2006) The effect of contact conditions and material properties on the elasticity terminus of a spherical contact. Int J Solids Struct 43:5736–5749. https://doi.org/10.1016/j.ijsolstr.2005.07.034

40. Shankar S, Mayuram MM (2008) A finite element based study on the elastic-plastic transition behavior in a hemisphere in contact with a rigid flat. J Tribol. https://doi.org/10.1115/1.2958081

41. Sahoo P, Chatterjee B (2010) A finite element study of elastic-plastic hemispherical contact behavior against a rigid flat under varying modulus of elasticity and sphere radius. Engineering 02:205–211. https://doi.org/10.4236/eng.2010.24030

42. Wu C, Li L, Thornton C (2003) Rebound behaviour of spheres for plastic impacts. Int J Impact Eng 28:929–946. https://doi.org/10.1016/S0734-743X(03)00014-9

43. Wu C-Y, Li L-Y, Thornton C (2005) Energy dissipation during normal impact of elastic and elastic–plastic spheres. Int J Impact Eng 32:593–604. https://doi.org/10.1016/j.ijimpeng.2005.08.007

44. Mukhopadhyay S, Das PK, Abani N (2023) A theoretical model to predict normal contact characteristics for elasto-plastic collisions. Granul Matter 25:20. https://doi.org/10.1007/s10035-023-01307-0

45. Sulsky D, Zhou SJ, Schreyer HL (1995) Application of a particle-in-cell method to solid mechanics. Comput Phys Commun. https://doi.org/10.1016/0010-4655(94)00170-7

46. de Vaucorbeil A, Nguyen VP, Sinaie S, Wu JY (2020) Material point method after 25 years: Theory, implementation, and applications. In: Advances in Applied Mechanics. pp 185–398

47. Sulsky D, Schreyer HL (1996) Axisymmetric form of the material point method with applications to upsetting and Taylor impact problems. Comput Methods Appl Mech Eng 139:409–429. https://doi.org/10.1016/S0045-7825(96)01091-2

48. Li F, Pan J, Sinka C (2009) Contact laws between solid particles. J Mech Phys Solids 57:1194–1208. https://doi.org/10.1016/j.jmps.2009.04.012

49. Ye Z, Zhang X, Zheng G, Jia G (2018) A material point method model and ballistic limit equation for hyper velocity impact of multi-layer fabric coated aluminum plate. Int J Mech Mater Des 14:511–526. https://doi.org/10.1007/s10999-017-9387-0

50. Nezamabadi S, Nguyen TH, Delenne J-Y, Radjai F (2017) Modeling soft granular materials. Granul Matter 19:8. https://doi.org/10.1007/s10035-016-0689-y

51. Nezamabadi S, Radjai F, Averseng J, Delenne JY (2015) Implicit frictional-contact model for soft particle systems. J Mech Phys Solids. https://doi.org/10.1016/j.jmps.2015.06.007

52. Nezamabadi S, Frank X, Delenne J-Y et al (2019) Parallel implicit contact algorithm for soft particle systems. Comput Phys Commun 237:17–25. https://doi.org/10.1016/j.cpc.2018.10.030

53. Nezamabadi S, Ghadiri M, Delenne JY, Radjai F (2022) Modelling the compaction of plastic particle packings. Comput Part Mech. https://doi.org/10.1007/s40571-021-00391-4

54. de Vaucorbeil A, Nguyen VP (2021) Modelling contacts with a total Lagrangian material point method. Comput Methods Appl Mech Eng 373:113503. https://doi.org/10.1016/j.cma.2020.113503

55. Nezamabadi S, Radjai F, Averseng J, Delenne J-Y (2015) Implicit frictional-contact model for soft particle systems. J Mech Phys Solids 83:72–87. https://doi.org/10.1016/j.jmps.2015.06.007

56. Nezamabadi S, Radjai F (2023) Explicit total Lagrangian material point method with implicit frictional-contact model for soft granular materials

57. Deresiewicz H (1968) A note on Hertz’s theory of impact. Acta Mech 6:110–112. https://doi.org/10.1007/BF01177810

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.