# Using indentation to identify the parameters of a numerical third-body model: Application to a rubber-like material

Kévin Daigne, Guilhem Mollon, Sylvie Descartes, Nicolas Fillot, Paul Baral, Romain Jeanneret-Dit-Grosjean, Frederic Biesse, Antoine Perriot, Gaëtan Maurel

► To cite this version:

Kévin Daigne, Guilhem Mollon, Sylvie Descartes, Nicolas Fillot, Paul Baral, et al.. Using indentation to identify the parameters of a numerical third-body model: Application to a rubber-like material. Tribology International, 2025, 212 (110950), 10.1016/j.triboint.2025.110950. emse-05190306

HAL Id: emse-05190306

https://hal-emse.ccsd.cnrs.fr/emse-05190306v1

Submitted on 28 Jul 2025

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
due to adhesion forces, which can be significant. For example, in [9], indentations were performed on polystyrene and polyethylene. It was shown that adhesion is higher with a diamond tip than with a tungsten tip, resulting in a significant error if not accounted for (a hardness 2 to 3 times higher for the diamond tip). This must be taken into account during experimental investigations and will be discussed later.

A common way of extracting information from an indentation test is to find properties by means of an analytical solution. There are several approaches, but most of them are based on the Oliver–Pharr method [10]. One of its main advantages is that it requires no measurements other than those made by the indenter. This consists in establishing a link between the initial slope at unloading, which can be extracted from the measurements, and the reduced Young modulus. For a soft material (compared to the Young modulus of the tip), this corresponds to Eq. (1).

$$S = \left. \frac{dP}{dh} \right|_{h=h_{max}} = \frac{2E}{1-v^2} \sqrt{\frac{A}{\pi}} \quad (1)$$

with S the contact stiffness, P the load, h the indentation depth, $h_{max}$ the maximum indentation depth, E the sample Young modulus, v the sample Poisson ratio and A the contact area. In cases where no analytical solution exists, reverse analysis can be carried out using numerical models.

The majority of the numerical models consider the material as continuous and use the finite element method (FEM). Although many works deal with elastic materials, some have focused on visco-elastic materials [11,12]. There are two main reasons to model the material as a packing of several bodies (i.e. using discrete modeling). The first reason is that it provides an access to fracture mechanisms for brittle materials, with the ability to follow the crack path [13,14]. The second reason is that the medium is granular, like powders [15]. In the current work, a discrete model will be used due to the final objective, which is a friction model to study mineral diffusion.

Indentations using a discrete model above the molecular scale are mainly performed with the discrete element method (DEM). However, the main limitation of DEM is well known: it is limited to small deformations at the grain scale [16]. Consequently, during loading, equilibrium will be satisfied by the displacement of the bodies, and a limited amount of elastic energy will be recovered during unloading. This implies that these models are more accurate for indentations mainly governed by irreversible processes.

In the present study, soft bodies will be used. This type of model combines the advantages of a continuous elastic medium with those of a discrete one. The originality of the work lies in developing an identification procedure compatible with elasto-plastic third-body models. This approach will therefore differ from conventional indentation DEM-based models.

# 1. Methods

## 1.1. Framework

The final objective of this study is to allow the modeling of mineral diffusion in a sliding interface (i.e. tire-road contact) [17,18]. Tire tread surface is modeled as a packing of separable bodies to allow mineral diffusion within it. Tire tread Young modulus is relatively low compared to the stress expected in a sliding interface. Consequently, large deformations must be computed, which is not possible under the DEM framework. In addition, the combined use of rigid (i.e. minerals) and highly deformable (i.e. rubber) bodies results in significant geometric distortions, which cannot be treated by FEM. The interactions between these very different bodies will require a robust contact detection algorithm. Finally, all these bodies must interact through consistent contact models and be subjected to tribological stresses. For all these reasons, a multibody meshfree approach implemented in a code called MELODY is used [19]. A plane-strain kinematics is assumed, which will

<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of the friction model.</figcaption></figure>

be discussed later. An illustration of the model (without the minerals as only the tire tread will be characterized in this study) can be found in Fig. 1.

Tire tread surface is modeled as a collection of soft bodies. To model large-scale phenomena (e.g. stick-slip events), the bulk of the tire tread is modeled as a single soft body (which is not fully displayed in the illustration for the sake of clarity). All subsequent modeling decisions are based on this model, with the objective of characterizing the tire tread material with the closest configuration. The tire tread physical behavior is described by a constitutive model and a contact model. In these models, three main parameters are involved (E the Young modulus, $C_N$ the normal cohesion coefficient and $C_T$ the tangential cohesion coefficient), as shown in Fig. 2 and described below.

A meshless method (similar to FEM) is employed to compute the deformation of each body. A quasi-incompressible neo-Hookean hyper-elastic model is used. The corresponding strain energy density function (W) is written in Eq. (2) with $C_1$ a material parameter and $I_1$ the first invariant of the right Cauchy–Green strain tensor [20]. This constitutive model well describes the behavior of rubber-like materials, even with large deformations [21,22]. A density of 940 kg m⁻³ is assumed. If a tensile test is performed, as shown in Fig. 2(a), the Young modulus (E) corresponds to the slope of the linear region of the stress–strain curve. Although the neo-Hookean model is non-linear and valid for large deformations, E will be used in the following analysis for the sake of clarity (as it is more commonly employed).

$$W = C_1(I_1 - 3) \quad (2)$$

Cohesion coefficient is used to compute the contact forces between the bodies. More details about the contact algorithm can be found in [23]. Contact between two bodies is detected if their contours overlap. The contact model defines the value of the resulting force as a function of several parameters. In the current study, normal ($C_N$) and tangential ($C_T$) cohesion coefficients are defined. By convention, these coefficients are positive. The normal cohesion coefficient is a tensile cut-off stress. It means that if the normal force per unit area ($\sigma_N$) on the contour exceeds this value, the two bodies are no longer considered in contact (i.e. are free to move relatively). An example is shown in Fig. 2(b) and the system of equations is written in Eq. (3) for a given node, with $k_N$ the normal contact stiffness, $\delta_N$ the normal interparticle distance of the node with the contour of another body (negative if interpenetrated and otherwise positive), $\delta_N^c$ the critical normal interparticle distance, and $A_C$ the contact area. This model can be viewed as a linear elastic link between a node and the contour of another body, which breaks if $\delta_N \ge \delta_N^c$ (i.e. $-\sigma_N \ge C_N$). This is therefore similar to what can be used in DEM with rigid grains, with
<figure><img src="image_3.png" /><figcaption>Fig. 2. Elementary simulations showing the three main parameters describing material behavior (E, C<sub>N</sub> and C<sub>T</sub>); the gray body represents the rigid frame in contact with a deformable body that is blue in the initial state and red in a more advanced state; with σ the tensile stress, ε the tensile strain, E the Young modulus, σ<sub>N</sub> the normal force per unit area exerted by the frame, C<sub>N</sub> the normal cohesion coefficient between bodies, t the time, t<sub>max</sub> the simulation duration, σ<sub>T</sub> the tangential force per unit area exerted by the frame and C<sub>T</sub> the tangential cohesion coefficient between bodies.</figcaption></figure>

the difference that here δ<sub>N</sub> represents a purely numerical parameter and will be minimized as much as possible while keeping computation time reasonable. The tangential cohesion coefficient is a resistive stress. It means that two bodies can move tangentially, without breaking the contact, if the tangential force per unit area (σ<sub>T</sub>) on the contour exceeds this value. An example is shown in Fig. 2(c) and the system of equations is written in Eq. (4) for a given node, with k<sub>T</sub> the tangential contact stiffness, δ<sub>T</sub> the tangential interparticle distance of the node with the contour of another body, δ<sup>c</sup><sub>T</sub> the critical tangential interparticle distance, and A<sub>C</sub> the contact area.

$$
(3) \quad \sigma_N = \begin{cases} -\frac{k_N \delta_N}{A_C}, & \text{if } \delta_N < \delta_N^c \\ 0, & \text{otherwise} \end{cases} \qquad (3a) \\
\delta_N^c = \frac{C_N A_C}{k_N} \qquad (3b)
$$

$$
(4) \quad \sigma_T = \begin{cases} -\frac{k_T \delta_T}{A_C}, & \text{if } |\delta_T| < \delta_T^c \text{ & } \delta_N < \delta_N^c \\ -C_T \operatorname{sign}(\delta_T), & \text{if } |\delta_T| \ge \delta_T^c \text{ & } \delta_N < \delta_N^c \\ 0, & \text{otherwise} \end{cases} \quad (4a) \\
\delta_N^c = \frac{C_N A_C}{k_N} \quad (4b) \\
\delta_T^c = \frac{C_T A_C}{k_T} \quad (4c)
$$

To reduce the number of parameters under study, normal and tan-
gential cohesion coefficients are assumed to be equal and are merged
into a single parameter, hereafter referred to as cohesion coefficient
(C). Note that each body contour is discretized as a succession of
nodes and segments (which can be compared to the mesh used in
FEM). The cohesion coefficient corresponds to the load per unit length
divided by the length of the element in contact. For the sake of brevity,
the cohesion coefficient will hereafter be referred to as cohesion. In
summary, two parameters will be identified: the Young modulus (E)
and the cohesion (C).

The characteristic time scale of the friction model (~10⁻⁷ s) is well
below the relaxation time of the rubber-like material (~10² s) [24].
Moreover, as a reminder, the cross-linking network is considered dam-
aged and the indentations will be performed on an uncured mate-
rial. Therefore, the viscous component will be considerable, further

increasing the relaxation time. In addition, tribological stresses are
significant in a sliding interface. All these points suggest a highly
elasto-plastic flow. Consequently, the effect of viscosity is neglected
in the indentation model, in order to identify only the elasto-plastic
component. However, such strain and loading rates cannot be reached
experimentally with indentation. Therefore, viscosity will unavoidably
be an experimental parameter. The model does not aim to capture
this viscosity dependency. The aim is to reproduce a given experi-
mental configuration, with assumptions that seem consistent at the
scale of the friction model. Consequently, the numerical indentations
will exhibit an elasto-plastic behavior, while the experimental one a
visco-elasto-plastic behavior.

The aim of reverse analysis is to compare a set of numerical and
experimental data. Thus, assuming a strictly monotonic dependence, it
is necessary to extract two experimental data in order to identify E and
C. In Fig. 3 is displayed an overview of the procedure according to the
experimental behavior of the material. The details of each stage will be
discussed in the corresponding sections.

For an elasto-plastic material, residual depth (h<sub>r</sub>) and hardness (H)
can be used. Residual depth corresponds to the depth at which contact
is lost between the tip and the surface. Hardness corresponds to the
ratio of the maximum load (P<sub>max</sub>) divided by the maximum contact area
(A<sub>max</sub>). For this type of material, reverse analysis is therefore performed
based on {h<sub>r</sub>, H} to identify {E, C}. The indentation model assumes
an elasto-plastic behavior and this method will be presented (with
numerical data only) to provide inputs for further studies. However,
the final procedure, which requires comparison with experimental data
on the viscous rubber-like material, will follow a different method.

For a visco-elasto-plastic material, h<sub>r</sub> has a limited physical meaning and the use of contact stiffness (S) will be preferred. For this type of material, an analytical solution (Eq. (1)) based on {S} will be used to determine {E} and the reverse analysis will be performed based on {H} to identify {C}. As rubber has a visco-elasto-plastic behavior, this method will be employed in this study. Note that for these indentations, a continuous stiffness measurement (CSM) mode is used, which consists in adding small oscillations to the main indentation. The analytical solution is then applied to each oscillation. This allows the Young modulus plateau to be reached and obtain the instantaneous Young modulus, which is more relevant for describing a third body, as will be discussed later.

1.2. Indentation model

An indentation model is created to perform a reverse analysis for
the material of the tire tread (i.e. rubber-like material), which is part
<figure><img src="image_4.png" /><figcaption>Fig. 3. Overview of the procedure according to the experimental material behavior with *P* the load, *h* the indentation depth, *h*<sub>r</sub> the residual depth, *P*<sub>max</sub> the maximum load, *h*<sub>max</sub> the maximum indentation depth, *H* the hardness, *A*<sub>max</sub> the maximum contact area, CSM for continuous stiffness measurement, *S* the contact stiffness, *ν* the Poisson ratio and *A* the contact area.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Illustration of the indentation model; *h*<sub>max</sub> corresponds to the maximum indentation depth.</figcaption></figure>

of the third body. The objective is to identify *E* and *C* while keeping certain properties (constitutive model, plane-strain kinematics, etc.) of the friction model. The physical meaning of *C* strongly depends on the configuration (friction or indentation). In the context of indentation, it corresponds to the material ability to flow. This will be discussed further, but it may encompass both viscous and plastic effects. Indentations are performed under quasi-static conditions. The velocity of the tip (hereinafter referred to as the displacement rate) during loading and unloading is equal and constant (with a smooth transition between the two). The model assumes an elasto-plastic behavior. The maximum indentation depth (*h*<sub>max</sub>) is equal to 1 µm. An illustration of the model can be found in Fig. 4.

Near the indenter tip, the rubber-like material is modeled as a packing of 1463 circular soft bodies. The diameters follow a Gaussian distribution with a mean value (*D*<sub>M</sub>) of 0.15 *h*<sub>max</sub> and a standard deviation of 0.03 *h*<sub>max</sub>. This distribution is close to the one used in the friction model. Moreover, the *h*<sub>max</sub> value was chosen to perform the indentations on a significant number of bodies while keeping computation times affordable. The discrete domain is bounded by a

continuous soft body. The prescribed displacement is zero on its lower boundary.

The discrete bodies and the continuous one have the same local mechanical properties (different properties emerge if the whole domain is considered). These properties, which are described in the constitutive model (i.e. neo-Hookean) and contact model (i.e. cohesive), are identical to those used in the friction model. The tip is considered rigid. Contact between the tip and the rubber-like material is assumed to be frictionless. Indeed, friction for a high tip angle, such as the one used in the model, remains low [25,26]. Tip adhesion is also neglected to ensure the unicity of the solution. This point will be discussed in the experimental section.

Experimental indentations are performed using a Berkovich tip. To simulate this type of tip, the majority of 2D models use axisymmetric simulations (whereas the present method can only model plane-strain kinematics). This consists in assuming a conical tip with a half apex angle of 70.3°. This gives the same contact area (A) as a 3D model for a given indentation depth (h). Several studies have shown that a 2D axisymmetric model faithfully reproduces some of the results (e.g. load-depth curves) obtained by a 3D model [27,28]. Note that the contact area is the area projected onto the plane orthogonal to the indentation axis. However, the previous technique cannot be used with the current model (2D with plane-strain assumption). Indeed, in this case, the contact area varies linearly (A ∝ h), whereas for a Berkovich tip (3D or 2D with axisymmetric assumption), it varies quadratically (A ∝ h²). Nevertheless, these studies show that for a given configuration, the results depend mainly on the contact area. For this reason, the reverse analysis will be carried out on quantities that do not depend on this area (e.g. load per unit area). For an elasto-plastic material, residual depth and hardness can be used. As a reminder, residual depth (*h*<sub>r</sub>) corresponds to the depth at which contact is lost between the tip and the surface. Hardness (*H*) corresponds to the ratio of the maximum load (*P*<sub>max</sub>) divided by the maximum contact area (*A*<sub>max</sub>). For the sake of comparison, the half apex angle will be equal to 70.3°.

In order to replicate the properties of the third body described previously, the indentation model relies on two main assumptions. The first is that plasticity is modeled by the relative displacements of the
<figure><img src="image_6.png" /><figcaption>Fig. 5. Main preliminary stages.</figcaption></figure>

bodies, and not in the bulk of each of them. The second assumption concerns viscosity. Indeed, as explained above, it is assumed that for the friction model, the viscosity of the material can be neglected. The same assumption is therefore made for the indentation model for the sake of consistency, although experimentally such strain or loading rates can- not be reached. For all these reasons, the discrete domain represents an elasto-plastic region where relative internal displacements are allowed, whereas the continuous domain is purely elastic.

A preliminary study was carried out to determine the size of these domains. The continuous domain has a semicircle shape with a radius of 100 $h_{max}$. This radius allows to consider the domain as infinitely large (the error on maximum load being 1.6% for an elastic case). In addition, a smooth transition between the discrete and continuous domains must be guaranteed. To this end, a sinusoidal surface roughness at the interface of these domains is defined, with an amplitude of 1 $D_M$ and a period of $2D_M$. This ensures a good force transmission between the two domains while using the same contact model as the one used between the discrete bodies (and therefore preserving the unicity of the solution). Before each indentation, a preliminary procedure is carried out, as shown in Fig. 5.

The bodies are initially free of contact and cohesionless (Fig. 5(a)). The discrete bodies are compacted by a top plate inside the discrete domain area, bounded by the continuous domain (Fig. 5(b)). The cohesion is then activated and the top plate is removed (Fig. 5(c)).

Before performing a reverse analysis, the objective is to demonstrate the unicity of the solution by defining a parametric space {$E, C$}. As discussed later, cohesion can be reduced to a normalized cohesion equal to $\frac{C}{E}$. The parametric space is composed of 12 simulations, with 3 Young moduli (1, 5 and 10 MPa) and 4 cohesions (0.05, 0.10, 0.20 and 0.40 E). Based on these simulations, a procedure for the reverse analysis is proposed in Section 3. An alternative procedure, using the same method, is also proposed for finding E and C without a pre-existing study and with a few iterations.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Surface topography measured by interferometry.</figcaption></figure>

## 1.3. Experimental set-up

Indentations are performed on a Nanoindenter XP (KLA Corporation) equipped with a dynamic contact module head (DCM V1), allowing accurate measurements at low displacements and loads. A Berkovich diamond tip is used (three-sided pyramid with two reference angles: 65.3° and 76.9°). The tip defect is taken into account using the area function given in [29]. The sample is an uncured rubber-like material. The composition is as follows (in phr<fnref n="1" />): 100 isoprene rubber (IR), 50 N347, 3 6PPD, 1.5 stearic acid, 3 ZnO, 1.3 CBS and 1.3 sulfur. Note that sulfur is added even if no vulcanization is performed. This choice is made to be as close as possible to a damaged cross-linking network. The sample is 2.8 mm thick and fixed to a metallic support using a heat-softening glue (CRYSTALBOND 555).

The current material is very adhesive and soft, which leads to several difficulties. Firstly, the adhesive force is close to the withdrawal force capacity of the indenter. Secondly, the low Young modulus combined with stress relaxation makes surface detection difficult with a potential gap between the detected and the actual surface. If this is not taken into account in the reverse analysis, it will lead to a significant error. To overcome these difficulties, the surface is coated with a 15 nm thick layer of gold, using a Safematic CCU-010 compact coating unit. This thickness is a compromise to minimize adhesion effect without significantly affecting the results. Note that the final reverse analysis will be performed on hardness (i.e. before the unloading for which adhesion is most significant). Therefore, although adhesion may have an effect on the contact area, it will be considered negligible thanks to the coating.

Several characterizations of the topography were carried out by interferometry using a 3D non-contact optical surface profiler (Zygo, ZeGage TM Pro). An example of one of these surfaces is shown in Fig. 6. It should be noted that in some cases, the coating did not seem to adhere properly to the surface (a wavy pattern was observed). These samples were withdrawn from the study. For surfaces of 174 by 174 µm after gold coating, the arithmetic mean height is on average equal to 0.63 µm. For an indentation depth of 1 µm, the contact area is approximately $10^3$ times smaller than the area over which the roughness was measured. Since roughness decreases with the size of the measured surface (except in the presence of isolated irregularities), roughness can reasonably be considered negligible at the indentation scale. In addition, several tests will be carried out at various locations on the sample to account for potential local heterogeneities.

First, conventional indentations are performed. A constant and identical displacement rate is applied during loading and unloading (without holding time). The maximum indentation depth ($h_{max}$) is equal to 1 µm. This configuration is chosen for the sake of comparison with

<footnote n="1">¹ parts per hundred of rubber</footnote>
the numerical model. Indentations are performed for 9 displacement rates (V), respectively equal to: 10, 20, 30, 40, 50, 60, 100, 200 and 500 nm s⁻¹. The objective is to determine the sensitivity of the material to strain rate (ϵ). The strain rate is approximated by $\frac{V}{h_{max}}$. In the post-processing, only the stiffness of the tip holder is corrected. Moreover, thermal drift measurements (using a holding time at reduced load) could not be performed, as the material would continue to deform due to its viscosity. Consequently, experiments are carried out overnight to reduce the potential sources of drift, and no correction of the thermal drift is made. In the reverse analysis procedure, these results will be used to compare numerical and experimental hardness and identify cohesion.

Then, indentations using continuous stiffness measurement (CSM) are performed. This consists in adding small oscillations to a main indentation. These tests make it possible to carry out a large number of indentations without complete unloading, and to extract the sensitivity of material properties according to the loading frequency (expressed as a strain rate). This is equivalent to performing successive indentations for an increasing depth, allowing Young modulus to be derived analytically as a function of depth (using Eq. (1)). The load is controlled to ensure a constant strain rate during indentation (displacement has an exponential form). Indentations will be carried out for two different strain rates with a value of 0.05 and 0.25 s⁻¹ respectively. It is important to note that for this type of CSM tests, the strain rate that is relevant for the modulus is the one at the oscillation scale. Those mentioned above correspond to the scale of the main indentation (calculated in the same way as conventional indentations). Indeed, it will be shown that the derived Young modulus corresponds to an instantaneous modulus ($E_0$) and is independent of these strain rates. The maximum indentation depth is set to 2 µm (for the CSM tests only) and oscillations with an amplitude of 2 nm at a frequency of 75 Hz are added. In the reverse analysis procedure, these results will be used to determine the Young modulus (by assuming $E \approx E_0$).

## 2. Numerical results

### 2.1. Kinematics

Ensuring a consistent kinematics (i.e. velocity field) is mandatory for a reverse analysis, otherwise an arbitrary set of inputs might erroneously lead to the expected target output. For irrelevant kinematics, there is no guarantee that if the configuration changes, the set of inputs (i.e. *E* and *C*) will be in agreement with the material properties. Kinematics can be studied by examining, after indentation, the residual displacement field $U_r$ (which corresponds to the norm between the initial and the final state) and by tracking centroid positions. The results are presented in Fig. 7 (the 0.40 *E* cohesion is not displayed as it is very similar to the 0.20 *E* cohesion).

As discussed earlier, plasticity is described as irreversible displacements of matter, which corresponds to the bodies in this model. Consequently, $\frac{U_r}{h_{max}}$ can be considered as an analogous quantity to plasticity. Simulations range from elasto-plastic to elastic indentations. Indeed, for the lowest cohesion, $U_r$ is locally higher than 0.70 $h_{max}$, whereas for the highest cohesion it is closer to 0.05 $h_{max}$ in the whole domain. For a cohesion of 0.05 *E*, the highest $U_r$ values are located close to the tip and the bodies form a bump on the surface. This bump decreases with increasing cohesion. A small asymmetry can be observed with $U_r$ higher on the right side. This is attributed to the size and position of the bodies, which are slightly unevenly distributed. For cohesions of 0.05, 0.10, 0.20 and 0.40 *E*, the maximum $U_r$ is respectively equal to 0.77, 0.57, 0.15 and 0.11 $h_{max}$. However, for cohesions of 0.20 and 0.40 *E*, $U_r$ is only above 0.05 $h_{max}$ for a small volume. It therefore appears negligible. These results are consistent with other studies using continuous models for elasto-plastic indentations. For example, [30] used FEM with elasto-plastic solids. This study shows that the bump, often referred to as the pile-up phenomenon, is proportional to the ratio

of yield strength divided by the Young modulus, which is closely related to $\frac{C}{E}$. For an elasto-plastic case, *C* can therefore be considered as a yield strength.

To get a better view of material flow, several centroid positions during indentation are represented by lines in Fig. 7. The color gradient corresponds to the time (*t*) normalized by the time at which the tip has returned to its initial position (*t*<sub>max</sub>). For a cohesion of 0.05 *E*, a displacement in the normal ($\bar{n}$) and tangential ($\bar{t}$) directions can be observed. For bodies below the tip apex, normal displacement increases during loading (*t* < 0.5 *t*<sub>max</sub>) and decreases during unloading (*t* > 0.5 *t*<sub>max</sub>). For bodies more distant from the apex and close to the surface, the normal displacement increases continuously. Regarding tangential displacement, it increases continuously (or is null) during loading and is constant during unloading. This shows that tangential flow is mainly plastic (there is no tangential elastic recovery) whereas normal flow is elastic or elasto-plastic. To support this analysis, it is possible to look at the cohesion of 0.20 *E*. This indentation is mainly elastic. Indeed, centroid positions during unloading overlap those during loading. In this simulation, tangential flow is negligible compared to normal flow. Finally, the cohesion of 0.10 *E* is an intermediate case in which tangential flow is non-zero, but lower than the cohesion of 0.05 *E*. It can also be noted that the displacement of deeper bodies is lower for lower cohesion. This is explained by a decrease in load for an increasing degree of plasticity and will be discussed later.

Several studies have investigated the dependence of the velocity field on several parameters (material, tip geometry, etc.). For example, in [31] is measured how the velocity field is modified as a function of the apex angle for a wedge tip. Indentations are made on copper and the velocity field is acquired by image correlation using a high-speed camera. The results show that below the tip apex, material flow is mainly normal. On the sides of the tip, the tangential component increases compared to the normal one. These results are consistent with the kinematics found in the current model.

Bodies located close to the tip apex have been subjected to significant stresses which do not result in relative displacements. Consequently, during unloading, these bodies will recover their initial strains and positions. If plasticity were modeled within the bulk of each body, $U_r$ under the tip apex would be higher [32,33]. However, the effect of this plastic deformation only concerns a small volume and will be considered negligible.

### 2.2. Load-depth curves

For several cohesion and Young modulus values, the corresponding load-depth curves are displayed in Fig. 8. Normal load (*P*) is normalized by *E* and the maximum contact area ($A_{max} = 5.58h_{max}L$) with *L* the length in the transverse direction (i.e. one unit of length). Indentation depth (*h*) is normalized by $h_{max}$. If the ratio $\frac{C}{E}$ is kept constant, the load-depth curves overlap. Firstly, this shows that as $A_{max}$ is constant for all simulations, *P* is linearly proportional to *E*. Secondly, if $\frac{C}{E}$ varies, then *P* is modified. These properties will be discussed in detail in the reverse analysis section. Furthermore, there is no significant change in curvature at the very beginning of each indentation. This shows that even though the contact is applied to a few bodies, the results remain consistent with higher *h*. For cohesions of 0.05 (Fig. 8(a)) and 0.10 (Fig. 8(b)) *E*, results are noisier during the second half of loading. This is attributed to local reorganization reinforced by low cohesion values.

During the loading, the slope of the load-depth curves changes as a function of several parameters. Firstly, plastic flow leads to a decrease in the slope. As shown in Fig. 7, material tends to flow tangentially for elasto-plastic indentations. As a result, the volume of material under the tip also decreases, and therefore the load. This can be assessed by looking at the maximum load. For cohesions of 0.05, 0.10, 0.20 and 0.40 *E*, the normalized maximum load is on average equal to 0.022, 0.035, 0.045 and 0.047 respectively. The second parameter affecting the slope is the contact area. For a wedge tip, the contact area increases
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Kinematics according to normalized cohesion ($E = 5$ MPa); the background corresponds to the residual displacement field ($U_r$) normalized by the maximum indentation depth ($h_{max}$); $U_r$ is equal to the centroid displacement norm between the initial and final states; the lines correspond to various centroid positions during indentation and the color to time ($t$) normalized by the time at which the tip has returned to its initial position ($t_{max}$).</figcaption></figure>

linearly during loading ($A \propto h$). Therefore, in the absence of plasticity, the slope should increase with $h$. For example, the increase in slope can be observed for a cohesion of 0.40 E (Fig. 8(d)), which is mainly elastic. Unloading, like loading, depends on plasticity and contact area. However, as discussed in the kinematics section, tangential flow only occurs during loading. A key quantity during the unloading is the

residual depth. This depth provides an evaluation of the degree of plasticity.

From a cohesion of 0.20 E, the results start to become less sensitive to normalized cohesion. This gives an upper limit to the normalized cohesion that can be identified by the model. Furthermore, for these cohesions, the unloading and loading stages are not perfectly overlapping.
<figure><img src="image_9.png" /><figcaption>Fig. 8. Numerical load-depth curves according to the Young modulus (1, 5 and 10 MPa) and for several normalized cohesions (0.05, 0.10, 0.20 and 0.40).</figcaption></figure>

Additional simulations show that this gap can be reduced by decreasing the displacement rate. Inertial forces must also be negligible compared to cohesive forces, otherwise non-physical decohesions close to the surface will occur. This gives a lower limit to the normalized cohesion, which is mainly limited by computation time. However, lower values are beyond the scope of this study. From the point of view of a friction model, lower cohesions would be well below the typical stress level and would therefore be negligible.

### 3. Reverse analysis

#### 3.1. Elasto-plastic behavior

This section refers to numerical data only, and may correspond to any material with a usual elasto-plastic behavior (i.e. not suitable for rubber). The results depend on two main parameters: *E* and *C*. As it could be seen on the indentation curves (Fig. 8), these two parameters have an effect on the load and the residual depth. However, as previously mentioned, it is not possible to directly compare numerical and experimental load due to plane-strain kinematics. For this reason, it was decided to only use *H* (hardness) and *h*<sub>r</sub> (residual depth). As the experimental real contact area is unknown, for the sake of comparison, *H* is assumed to be equal to *F*<sub>max</sub>/<i>A</i><sub>max</sub>. Note that using results at *h*<sub>max</sub> reduces the inertial effects, as the velocity gradually decreases at this moment. In addition, it maximizes the number of involved bodies which reduces the dispersion (less dependent on one body position). The values of *H* and *h*<sub>r</sub> as a function of *E* and *C* are displayed in Fig. 9. Data highlighted in red with the superscript “exp” for experimental, refer to data obtained experimentally. Data highlighted in blue refer to numerical data only.

In Fig. 9(a), dependency of normalized hardness on normalized cohesion is displayed. Hardness is normalized by the Young modulus. It can be seen that $\frac{H}{E}$ increases as $\frac{C}{E}$ increases, until it reaches an asymptote. This trend is consistent with what is obtained in finite element simulations, where a yield strength is defined instead of cohesion [34,35]. This reinforces the physical validity of the discrete

model and supports the interpretation of the cohesion parameter as analogous to a yield strength. Furthermore, one of the main challenges of reverse analysis is to guarantee the unicity of the solution. For values of *C* greater than 0.20 *E*, hardness becomes independent of cohesion, implying that the unicity criterion is no longer satisfied.

In Fig. 9(b), dependency of normalized hardness on normalized residual displacement is displayed. As *H* is almost linearly proportional to *E* (data points overlap closely for a given normalized cohesion in Fig. 8), a change in $\frac{H}{E}$ describes how hardness is influenced by plasticity, characterized by $\frac{h_r}{h_{max}}$. The ratio $\frac{H}{E}$ decreases with increasing $\frac{h_r}{h_{max}}$. The slope of the regression is equal to -0.043. For the lowest normalized cohesion compared to the highest one, hardness is reduced by half. The ratio $\frac{h_r^{exp}}{h_{max}^{exp}}$ can be obtained from the experimental data to find the corresponding $\frac{H}{E}$. Then, using the experimental hardness $H^{exp}$, *E* can be derived.

In Fig. 9(c), dependency of normalized cohesion on normalized residual displacement is displayed. Residual depth increases as cohesion decreases. For cohesions of 0.05, 0.10, 0.20 and 0.40 *E*, $h_r$ is on average respectively equal to 0.64, 0.42, 0.13 and 0.11 $h_{max}$. For the same normalized cohesion, lower Young moduli tend to give lower $h_r$. One way of minimizing this effect would be to take a depth before loss of contact (e.g. 90% of unloading). Indeed, $h_r$ is very sensitive to local reorganization near the surface. In addition, the displacement rate could also be reduced. As mentioned above, the ratio $\frac{h_r^{exp}}{h_{max}^{exp}}$ can be obtained from the experimental data, which allows to find the corresponding $\frac{C}{E}$. Since *E* is already known, *C* can be derived.

Note that the method is presented for a pre-existing numerical campaign exploring the parametric space. Without these results, an iterative way to proceed is described in Algorithm 1.

With this method, a close result for *E* is expected after 2–3 iterations, as $\frac{H}{E}$ as a function of $\frac{h_r}{h_{max}}$ should be linear. A few more iterations will be needed to obtain *C*. The use of ratios is not mandatory and a more classical minimization method can be used on $\{|\frac{h_r}{h_{max}} - \left(\frac{h_r}{h_{max}}\right)^{exp}|, |H - H^{exp}|\}$. However, particular attention must
<figure><img src="image_10.png" /></figure>

(a) Normalized hardness as a function of normalized cohesion; the dotted line connects the mean value for each normalized cohesion

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

(b) Normalized hardness as a function of normalized residual depth; the dotted line corresponds to a linear regression performed on all data points

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

(c) Normalized cohesion as a function of normalized residual depth; the dotted line connects the mean value for each normalized cohesion

<caption>Fig. 9. Method for identifying *E* and *C* from reverse analysis; markers correspond to numerical data only; numerical parameters are highlighted in blue, and those derived from measurements in red.</caption>

**Algorithm 1 Reverse analysis procedure**

Require: {$H^{exp}, (\frac{h_r}{h_{max}})^{exp}$}

1: Assume an initial Young modulus $E^0$

2: Take the initial parameters as follows:

$$k-1: \{E^{k-1} = E^0, C^{k-1} = 0.05 E^0\}$$

$$k: \{E^k = E^0, C^k = 0.40 E^0\}$$

with *k* the iteration number

3: Complete steps 4, 5 and 6 for k-1

4: Perform a simulation for {$E^k, C^k$}

5: Extract the following numerical ratios:

$$\{(\frac{h_r}{h_{max}})^k, (\frac{H}{E})^k, (\frac{C}{E})^k\}$$

6: Calculate $dh$:

$$dh^k = (\frac{h_r}{h_{max}})^k - (\frac{h_r}{h_{max}})^{exp}$$

7: Calculate $f'$: $(f')^k = \frac{dh^k-dh^{k-1}}{(\frac{H}{E})^k-(\frac{H}{E})^{k-1}}$

8: Calculate the expected normalized hardness:

$$(\frac{H}{E})^* = (\frac{H}{E})^k - \frac{dh^k}{(f')^k}$$

9: Update the Young modulus:

$$E^{k+1} = H^{exp} (\frac{H}{E})^*^{-1}$$

10: Calculate $g'$:

$$g'^k = \frac{dh^k-dh^{k-1}}{(\frac{C}{E})^k-(\frac{C}{E})^{k-1}}$$

11: Calculate the expected normalized cohesion:

$$(\frac{C}{E})^* = (\frac{C}{E})^k - \frac{dh^k}{(g')^k}$$

12: Update the cohesion:

$$C^{k+1} = E^{k+1} (\frac{C}{E})^*$$

13: Repeat steps 4 to 12 until the criterion (to be defined) is reached

be paid to the relevance of cohesion values (which must be positive and below 0.2–0.4 $E$).

## 3.2. Experimental results

Rubber-like materials have visco-elastic properties. Consequently, the mechanical response during indentation depends on the displacement rate. As a reminder, to evaluate this effect, indentations are performed for 9 displacement rates ($V$), respectively equal to: 10, 20, 30, 40, 50, 60, 100, 200 and 500 nm s⁻¹. The strain rate (ε̇) is approximated by $\frac{V}{h_{max}}$ with $V$ the displacement rate and $h_{max}$ equal to 1 µm. Results are presented in Fig. 10.

Indentation curves are displayed in Fig. 10(a). For the sake of clarity, only the mean load-depth curve has been plotted for each displacement rate. During loading, the more the strain rate increases, the more the load increases. For a depth below 50 nm (0.05 $h_{max}$), a change in curvature can be observed. This may correspond to gold coating (15 nm), rubber compounds (carbon black, zinc oxide, etc.) or tip defects. However, it remains at a small scale compared to $h_{max}$. The slope is almost constant during loading, except for a strain rate of 0.50 s⁻¹, with an inflection point at 0.27 $h_{max}$. Depth should decrease when the tip reaches 1 $h_{max}$. However, the indenter was less and less able to fully satisfy this condition for increasing displacement rates. During unloading, the load evolves until there is no longer a contact between the tip and the surface. The residual depth depends on several parameters. More importantly, residual depth is a contribution from plasticity and viscosity. Indeed, if the tip withdraws too fast, the viscous relaxation of the deformed surface cannot match the pace, and contact is quickly lost. The relevance of using $h_r$ to measure the visco-plastic effect will be discussed thereafter.

The dependency of hardness on strain rate can be observed in Fig. 10(b). As the real contact area is unknown, $H$ is assumed to be equal to
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Experimental results for several strain rates (0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.10, 0.20 and 0.50 s⁻¹).</figcaption></figure>

$P_{max}$ with $A_{max} = 24.56 h(P_{max})^2$. This area corresponds to an approximation in which adhesion, friction and plastic effects are neglected. The resulting hardness values may thus be slightly overestimated [36]. Note that the area is calculated for the *h* corresponding to $P_{max}$ (which may slightly differ from $h_{max}$). The more $\dot{\epsilon}$ increases, the more *H* increases. The standard deviation of hardness remains low, with a mean value of 0.04 MPa. A power law can be derived, giving an exponent (*m*) of 0.17. This exponent reflects the sensitivity of the material to strain rate. In contrast, for PC and PMMA, the *m* exponent is respectively equal to 0.05 and 0.10 [37]. The current material is therefore highly dependent on the strain rate.

It should be noted that the strain rates obtained here (≤ 0.50 s⁻¹) are much lower than those obtained in the friction model (which are closer to 10⁷ s⁻¹). The possibility of extrapolating the values obtained here over several orders of magnitude remains unknown, and would deserve further investigations. There are a few characterization techniques that can attain significant strain rates. For example, the Hopkinson technique allows strain rates of the order of 10⁴ s⁻¹ to be reached [38] but could hardly be applied on such a soft material.

In practice, the Young modulus of a visco-elastic material is time-dependent ($E(t)$ with *t* the time). When this type of material is subjected to a given strain, it exhibits an instantaneous response driven by the instantaneous Young modulus ($E_0$), followed by an evolution as a function of time. This instantaneous modulus therefore corresponds to $E(t \to 0)$. In the numerical models, *E* is constant and therefore not time-dependent. For the friction model, the characteristic time scale is

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Results according to indentation depth from CSM tests.</figcaption></figure>

close to 10⁻⁷ s and it will be assumed that $E \approx E_0$. This assumption will therefore also be made for the indentation model, which will be important during the reverse analysis procedure.

The Oliver–Pharr method mentioned above (Eq. (1)) can be used to find $E_0$. However, applying this method to indentation without holding time at $h_{max}$ and with a relatively slow unloading rate, leads to inaccurate results [39]. For this reason, it cannot be used to identify $E_0$ with the previous measurements and additional CSM tests have been carried out. One of the main advantages is that the derived Young modulus depends on the strain rate driven by the oscillations. This makes it possible to obtain much higher strain rates and reach a plateau for *E* (i.e. the instantaneous Young modulus) [40].

As a reminder, the maximum indentation depth is set to 2 µm, the strain rate to 0.05 and 0.25 s⁻¹, and oscillations with an amplitude of 2 nm at a frequency of 75 Hz are added. The results are presented in Fig. 11.

First, as the displacement rate is a bit different from the previous measurements (with constant displacement rates), the load in Fig. 11(a) differs slightly from the one in Fig. 10(a). For a strain rate of 0.05 s⁻¹, the load is higher (+29%) at 1 µm but is not significantly changed at 0.5 µm. In addition, the load increases with the strain rate, as previously described.

Regarding $E_0$ in Fig. 11(b), it can be seen that, unlike the load–displacement curve, $E_0$ does not depend on the strain rate. This is
because the modulus is derived from the nanoscillations added to the main indentation. The following comments therefore apply to both strain rates. The depth-dependent modulus at shallow depths is well documented in the literature. Even in the case of a homogeneous material, gradients appear due to several factors such as plasticity, tip defect, or roughness [41]. However, particularly in polymers, size effects can also originate from elastic mechanisms. In such cases, generalized elasticity theories, such as couple stress elasticity, offer a suitable framework, as they account for strain gradient effects through the introduction of an internal length scale [42]. It can be seen that beyond 0.5 µm, the variation in modulus is relatively small. Since the typical length scale of the friction model is close to 1 µm, gradient effects will therefore be considered non-dominant. This also indicates that the influence of the gold coating is negligible. Moreover, the influence of the coating is generally considered minimal when the indentation depth exceeds ten times the film thickness, with the resulting error on *E* typically limited to a few percent [43,44]. In our case, the reference depth of 1 µm corresponds to 66 times the coating thickness, further supporting this assumption. For a depth of 1 µm, *E*₀ is equal to 18 MPa and the standard deviation to 2 MPa. This is the value of Young modulus that will be considered and it will be discussed later.

### 3.3. Visco-elasto-plastic behavior

This section refers to numerical data involving an elasto-plastic behavior, but which must be compared with experimental data involving a visco-elasto-plastic behavior. To fully identify *E* and *C*, it is necessary to extract *H* and *h*<sub>r</sub> from the measurements. However, for visco-elastic materials, *h*<sub>r</sub> corresponds more to a contact loss depth, than a real measure of visco-plastic properties. Indeed, for a strain rate of 0.50 s<sup>-1</sup>, *h*<sub>r</sub> is close to 1 *h*<sub>max</sub> (cf. Fig. 10(a)). Therefore, if a reverse analysis were performed, it could lead to *H* being zero, *C* being zero and *E* being undefined (cf. Fig. 9), which does not seem accurate. However, it should be noted that the previous method remains perfectly suitable for elasto-plastic indentations. In this context, where the ultimate goal is to model a material that is unable to relax internal stresses over time, viscosity will be treated as plasticity, as previously discussed. Therefore, during indentation, viscous and plastic effects will both be encompassed by the parameter *C*. This parameter describes how the material behavior deviates from a purely elastic response and is only meaningful for times shorter than the relaxation time.

First, only the dependency of hardness on the visco-elasto-plastic effect will be used. Reverse analysis is performed using the results presented in Fig. 9(a) and Fig. 10(b). The power law in Fig. 10(b) gives $H^{\text{exp}}$ at a chosen $\dot{\epsilon}$, from which $\frac{H}{E}$ can be computed using a specified value of $E$. Finally, by performing a linear interpolation between the mean normalized hardness values (dotted line in Fig. 9(a)), the corresponding $\frac{C}{E}$ can be derived. This procedure is then repeated for several values of $\dot{\epsilon}$ and $E$. It should be noted that the results are extrapolated to higher Young moduli, as the data points in Fig. 9(a) overlap closely. The results can be found in Fig. 12.

The upper boundary of the surface delimits cohesions below 0.05 *E*, which are outside the scope of this study. The lower boundary delimits the hardness that cannot be obtained for a given Young modulus (*C* having no effect beyond 0.20–0.40 *E*). The non-unicity of the solution lies on this boundary. For a given $\dot{\epsilon}$, the more the Young modulus increases, the more the normalized cohesion decreases. Indeed, to keep $\dot{\epsilon}$ (i.e. *H*) constant, the plastic contribution must reduce *H* as it is increased by the elastic contribution. For a given *E*, the more $\dot{\epsilon}$ increases, the more the normalized cohesion increases. This is due to a decrease in plasticity to increase *H* without a change in *E*.

One solution to fully identifying *E* and *C* is to determine *E* using a different method that does not involve plasticity (e.g. dynamic mechanical analysis). In this study, it is proposed to assume that Young modulus is equal to the instantaneous one measured by CSM. Indeed, for a sliding interface, the shear rate easily exceeds $10^7$ s<sup>-1</sup>, suggesting

<figure><img src="image_15.png" /><figcaption>Fig. 12. Young modulus and cohesion as a function of strain rate.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. Results from the reverse analysis for $E = E_0 = 18$ MPa.</figcaption></figure>

that Young modulus is close to the instantaneous Young modulus. The distinction between plasticity and time-dependent modulus is not easy. It is highly dependent on the chosen time scale. For this reason, both plasticity and time-dependent modulus, for a given strain rate, will be included in *C*. By taking $E_0$ equal to 18 MPa, the constitutive relationship displayed in Fig. 13 is obtained. It corresponds to a horizontal slice of the surface in Fig. 12. For the tested strain rates, cohesion should therefore range between 0.050 and 0.105 *E*. For a strain rate of 0.4 s<sup>-1</sup>, a Young modulus of 18 MPa and a cohesion of 1.75 MPa is then obtained.

### 4. Conclusion

The final objective of this study is to identify the material parameters within a friction model studying a rubber-like material. To this end, a method based on a reverse analysis of an indentation test was proposed. To be as close as possible to the configuration of the friction model with the indentation one, an elasto-plastic behavior is assumed. The material being a viscous rubber-like material, for a given strain rate, experimental viscous indentations will therefore be modeled as elasto-plastic indentations. The material is modeled as a collection of discrete and deformable bodies. Its mechanical behavior is described by two main parameters: the Young modulus (*E*) of each particle and the cohesion coefficient (*C*) between these particles. The objective is to determine the values of *E* and *C*.

A reverse analysis must satisfy several conditions. A first condition is that the kinematics obtained numerically must be consistent with what is expected experimentally. It has been confirmed by comparing
the velocity field with data available in the literature. This allows to reinforce the validity of the numerical parameters if the configuration changes.

Several data can be extracted from the model, the main ones being hardness (H) and residual depth ($h_r$). The use of these data allows simulations using plane-strain kinematics to be compared with three-dimensional experiments. The numerical model is able to simulate a wide range of behaviors, from elastic to elasto-plastic indentations. The corresponding load-depth curves could be obtained. The results show that H depends on E and C. Indeed, if C decreases, the tangential flow increases. As a result, the volume of material under the tip is reduced and therefore H. In addition, there is a close link between this tangential flow and $h_r$. Finally, it has been shown that H is linearly proportional to E.

A second condition which must be satisfied for the reverse analysis concerns the unicity of the solution. It has been shown that E and C depend on H and $h_r$ in a non-asymptotic and strictly monotonic way (for C less than 0.20-0.40 E). Consequently, the unicity of the solution is satisfied for the range of tested C.

Several conventional indentations were experimentally performed on the material. Measurements show that H depends on the strain rate ($\dot{\epsilon}$) according to a power law. These measurements were supplemented by continuous stiffness measurement (CSM) to obtain the instantaneous Young modulus ($E_0$). Assuming that $E \approx E_0$ and based on the reverse analysis of H, the values of E and C as a function of $\dot{\epsilon}$ were determined. These results will be introduced in the friction model of Fig. 1 in future work, in order to investigate tribological phenomena in the tire-road contact.

## CRediT authorship contribution statement

**Kévin Daigne:** Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Guilhem Mollon:** Writing – review & editing, Supervision, Software, Methodology. **Sylvie Descartes:** Writing – review & editing, Supervision, Methodology. **Nicolas Fillot:** Writing – review & editing, Supervision, Methodology, Funding acquisition. **Paul Baral:** Writing – review & editing, Methodology, Investigation, Data curation. **Romain Jeanneret-Dit-Grosjean:** Writing – review & editing, Supervision. **Frederic Biesse:** Writing – review & editing, Supervision. **Antoine Perriot:** Writing – review & editing, Supervision. **Gaëtan Maurel:** Writing – review & editing, Supervision.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

The authors sincerely thank the ANRT and the Michelin Chair for their financial support for this research, as well as Bérangère Lesaint and the CLYM platform (FED 4092, www.clym.fr) for their technical support.

## Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.triboint.2025.110950.

## Data availability

Data will be made available on request.

## References

[1] Kreider ML, Panko JM, McAtee BL, Sweet LI, Finley BL. Physical and chemical characterization of tire-related particles: Comparison of particles generated using different methodologies. *Sci Total Environ* 2010;408(3):652-9.

[2] Adachi K, Tainosho Y. Characterization of heavy metal particles embedded in tire dust. *Environ Int* 2004;30(8):1009-17.

[3] Mollon G. Solid flow regimes within dry sliding contacts. *Tribol Lett* 2019;67(4):120.

[4] Coran AY. Chemistry of the vulcanization and protection of elastomers: A review of the achievements. *J Appl Polym Sci* 2003;87(1):24-30, *eprint*: https://onlinelibrary.wiley.com/doi/pdf/10.1002/app.11659.

[5] Smith RW, Veith AG. Electron microscopical examination of worn tire treads and tread Debris. *Rubber Chem Technol* 1982;55(2):469-82.

[6] Chen Z, Diebels S, Peter NJ, Schneider AS. Identification of finite viscoelasticity and adhesion effects in nanoindentation of a soft polymer by inverse method. *Comput Mater Sci* 2013;72:127-39.

[7] Klapperich C, Komvopoulos K, Pruitt L. Nanomechanical properties of polymers determined from nanoindentation experiments. *J Tribol* 2000;123(3):624-31.

[8] Oyen ML, Cook RF. Load-displacement behavior during sharp indentation of viscous-elastic-plastic materials. *J Mater Res* 2003;18(1):139-50, Publisher: Cambridge University Press.

[9] Grunlan JC, Xia X, Rowenhorst D, Gerberich WW. Preparation and evaluation of tungsten tips relative to diamond for nanoindentation of soft materials. *Rev Sci Instrum* 2001;72(6):2804, Publisher: American Institute of Technology.

[10] Oliver W, Pharr G. An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments. *J Mater Res* 1992;7(06):1564-83, Publisher: Cambridge University Press (CUP).

[11] Liu K, VanLandingham MR, Ovaert TC. Mechanical characterization of soft viscoelastic gels via indentation and optimization-based inverse finite element analysis. *J Mech Behav Biomed Mater* 2009;2(4):355-63.

[12] Ngan AHW, Wang HT, Tang B, Sze KY. Correcting power-law viscoelastic effects in elastic modulus measurement using depth-sensing indentation. *Int J Solids Struct* 2005;42(5):1831-46.

[13] Huang H, Detournay E. Discrete element modeling of tool-rock interaction II: rock indentation. *Int J Numer Anal Methods Geomech* 2013;37(13):1930-47, *eprint*: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nag.2114.

[14] André D, Jebahi M, Iordanoff I, Charles J-L, Néauport J. Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter. *Comput Methods Appl Mech Engrg* 2013;265:136-47.

[15] Pasha M, Hare C, Hassanpour A, Ghadiri M. Analysis of ball indentation on cohesive powder beds using distinct element modelling. *Powder Technol* 2013;233:80-90.

[16] Mora S, Vu T-L, Barés J, Nezamabadi S. Highly deformed grain: from the Hertz contact limitation to a new strain field description in 2D. *EPJ Web Conf* 2017;140:05011, Publisher: EDP Sciences.

[17] Daigne K, Mollon G, Descartes S, Fillot N, Jeanneret-Dit-Grosjean R, Biesse F, et al. The tire-road contact: A mechanical mixing seen as a shear-induced diffusive process. *Wear* 2024;546-547:205339.

[18] Daigne K. Shear-induced diffusion of road minerals within the tire tread: a solid flow modeling using a soft multibody approach [Ph.D. thesis], Insa Lyon; 2024.

[19] Mollon G. A multibody meshfree strategy for the simulation of highly deformable granular materials. *Internat J Numer Methods Engrg* 2016;108(12):1477-97, *eprint*: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.5258.

[20] Comellas E, Bellomo FJ, Oller S. A generalized finite-strain damage model for quasi-incompressible hyperelasticity using hybrid formulation. *Internat J Numer Methods Engrg* 2016;105(10):781-800, *eprint*: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.5118.

[21] Aloui S, El Yaagoubi M. Determining the compression-equivalent deformation of SBR-based rubber material measured in tensile mode using the finite element method. *Appl Mech 2021*;2(1):195-208, Number: 1 Publisher: Multidisciplinary Digital Publishing Institute.

[22] Fujikawa M, Maeda N, Yamabe J, Kodama Y, Koishi M. Determining stress-strain in rubber with in-plane biaxial tensile tester. *Exp Mech* 2014;54(9):1639-49.

[23] Mollon G. A unified numerical framework for rigid and compliant granular materials. *Comput Part Mech* 2018;5(4):517-27.

[24] Khajehsaeid H. Application of fractional time derivatives in modeling the finite deformation viscoelastic behavior of carbon-black filled NR and SBR. *Polym Test* 2018;68:110-5.

[25] Warren AW, Guo YB. Machined surface properties determined by nanoindentation: Experimental and FEA studies on the effects of surface integrity and tip geometry. *Surf Coat Technol* 2006;201(1):423-33.

[26] Tabor D. The hardness of metals. In: Oxford classic texts in the physical sciences, Oxford, New York: Oxford University Press; 2000.

[27] Lichinchi M, Lenardi C, Haupt J, Vitali R. Simulation of Berkovich nanoindentation experiments on thin films using finite element method. *Thin Solid Films* 1998;312(1):240-8.

[28] Shim S, Oliver WC, Pharr GM. A comparison of 3D finite element simulations for Berkovich and conical indentation of fused silica. *Int J Surf Sci Eng* 2007. Publisher: Inderscience Publishers.
[29] Oliver WC, Pharr GM. Measurement of hardness and elastic modulus by instru-
mented indentation: Advances in understanding and refinements to methodology.
J Mater Res 2004;19(1):3–20, Publisher: Cambridge University Press.

[30] Taljat B, Pharr GM. Development of pile-up during spherical indentation of elastic-plastic solids. Int J Solids Struct 2004;41(14):3891–904.

[31] Udupa A, Sundaram N, Sugihara T, Chandrasekar S. Direct in situ obser-
vation of deformation modes in wedge indentation of metals. Mater Trans
2019;60(8):1442–9.

[32] Dean J, Clyne B. Extraction of plasticity parameters from a single test using a spherical indenter and FEM modelling. Mech Mater 2016;105.

[33] Sakharova NA, Fernandes JV, Antunes JM, Oliveira MC. Comparison between Berkovich, Vickers and conical indentation tests: A three-dimensional numerical simulation study. Int J Solids Struct 2009;46(5):1095–104.

[34] Vandeperre L, Clegg WJ. The correlation between hardness and yield strength of hard materials. Mater Sci Forum 2005;492–493:555–60, Publisher: Trans Tech Publications Ltd.

[35] Cheng Y-T, Cheng C-M. Relationships between hardness, elastic modulus, and the work of indentation. Appl Phys Lett 1998;73(5):614–6.

[36] Choi Y, Lee H-S, Kwon D. Analysis of sharp-tip-indentation load-depth curve for contact area determination taking into account pile-up and sink-in effects. J Mater Res 2004;19(11):3307–15.

[37] Kermouche G, Loubet JL, Bergheau JM. A new index to estimate the strain rate sensitivity of glassy polymers using conical/pyramidal indentation. Phil Mag 2006;86(33–35):5667–77, Publisher: Taylor & Francis eprint: https://doi.org/10.1080/14786430600778682.

[38] Sobczyk K, Chmielewski R, Kruszka L, Rekucki R. Strength characterization of soils' properties at high strain rates using the Hopkinson technique—A review of experimental testing. Materials 2022;15(1):274, Number: 1 Publisher: Multidisciplinary Digital Publishing Institute.

[39] Cheng Y-T, Ni W, Cheng C-M. Determining the instantaneous modulus of viscoelastic solids using instrumented indentation measurements. J Mater Res 2005;20(11):3061–71.

[40] Cifuentes SC, Frutos E, Benavente R, González-Carrasco JL, Lorenzo V. Strain rate effect on semi-crystalline PLLA mechanical properties measured by instrumented indentation tests. Eur Polym J 2014;59:239–46.

[41] Xia Y, Bigerelle M, Marteau J, Mazeran P-E, Bouvier S, Iost A. Effect of surface roughness in the determination of the mechanical properties of material using nanoindentation test. Scanning 2014;36(1):134–49.

[42] Peng C, Zeng F. Modeling the indentation size effects of polymers, based on couple stress elasticity and shear transformation plasticity. Arch Appl Mech 2022;92(12):3661–81.

[43] Wald MJ, Considine JM, Turner KT. Measuring the elastic modulus of soft thin films on substrates. In: Proulx T, editor. In: Experimental and applied mechanics, vol. 6, New York, NY: Springer; 2011, p. 741–7.

[44] Jung Y-G, Lawn BR, Martyniuk M, Huang H, Hu XZ. Evaluation of elastic modulus and hardness of thin films by nanoindentation. J Mater Res 2004;19(10):3076–80.