# Interfacial mechanical behavior of epoxy-quartz: MD nanoindentation and nanoscratching study

Pengchang Wei, Zhen-Yu Yin, Pierre-Yves Hicher, Wangqi Xu

► To cite this version:

Pengchang Wei, Zhen-Yu Yin, Pierre-Yves Hicher, Wangqi Xu. Interfacial mechanical behavior of epoxy-quartz: MD nanoindentation and nanoscratching study. Acta Geotechnica, 2024, 20 (4), pp.1599-1620. 10.1007/s11440-024-02503-9. hal-05095012

HAL Id: hal-05095012

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

Submitted on 3 Jun 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
# Interfacial mechanical behavior of epoxy-quartz: MD nanoindentation and nanoscratching study

Pengchang Wei<fnref n="1" /> · Zhen-Yu Yin<fnref n="1" /> · Pierre-Yves Hicher<fnref n="2" /> · Wangqi Xu<fnref n="1" />

Received: 18 June 2024 / Accepted: 5 December 2024
© The Author(s) 2024

## Abstract

Fiber-reinforced polymer (FRP) is widely used in various engineering fields due to its several outstanding properties. In geotechnical engineering, the interactions between FRP and soil play an essential role. In this paper, molecular dynamics (MD) simulation method has been performed to study the interfacial mechanical behavior of epoxy-quartz interface as a subsystem of FRP-soil structure. Uniaxial traction on bulk epoxy was conducted to verify the accuracy of the model. The nanoindentation and nanoscratching mechanisms of epoxy-quartz interface were analyzed, considering the effect of loading rate, sliding velocity, and indentation depth. Abrasion models have been proposed based on the relationship between forces and displacements. Simulation results indicated that the indenter force and the indentation hardness of epoxy substrate increased with the loading rate during nanoindentation, and the relationship between indenter force and indentation depths could be expressed by a power law. The forces along three directions increased with the sliding velocity or indentation depths during nanoscratching, the sliding force and the sliding distance following an exponential function. The numerical simulations demonstrated that the surface wear of the epoxy substrate had the shape of a groove in nanoindentation and a fan-shaped distribution during the nanoscratching process.

**Keywords** FRP-soil interface · Loading rate · Molecular dynamics · Nanoindentation · Nanoscratching · Sand

## 1 Introduction

Fiber-reinforced polymer (FRP) composites are composed of fibers within a polymer matrix, where fibers are usually made of glass, carbon, or aramid, and polymer matrix are epoxy resin, polyester resin, or polypropylene, etc. Due to

several outstanding properties, such as high-strength, low weight, corrosion resistance, electromagnetic transparency, and ease of handling, FRP composites have been widely used in various fields, such as marine engineering, building structures, aerospace, and civil and hydraulic engineering. [2, 19] In civil engineering, FRP composites are effectively applied to protect, reinforce, and repair structures. For example, compared to traditional concrete and steel piles with poor corrosion resistance, FRP piles are constantly used in corrosive environments (e.g., marine geotechnical engineering) due to their corrosion and chemical resistance properties [41], such as concrete-filled FRP tube (CFFT) pile [1] and FRP seawater sea-sand concrete (SSC) pile [7]. Thus, it could be an alternative construction material for piles in hostile environments [35]. Moreover, the interface friction coefficient between FRP composite and sand was higher than that between conventional mild steel and sand [35]. However, it should be noted that the stiffness and surface hardness of FRP composites are significantly lower than those of concrete and steel structures, which can affect the pile-soil interaction characteristics for FRP piles. Furthermore, with the promotion and application of FRP piles

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

Zhen-Yu Yin
zhenyu.yin@polyu.edu.hk

Pengchang Wei
pengchang.wei@polyu.edu.hk

Pierre-Yves Hicher
pierre-yves.hicher@ec-nantes.fr

Wangqi Xu
wangqi.xu@connect.polyu.hk

<footnote n="1">¹ Department of Civil and Environmental Engineering, Research-Centre-for Nature-based. Urban Infrastructure Solutions, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong, China</footnote>

<footnote n="2">² Ecole Centrale de Nantes, UMR CNRS 6183, Research Institute of Civil Engineering and Mechanics (GeM), Nantes, France</footnote>
in foundations, the interfacial characteristics between FRP composites and soil significantly impact the integrity and stability of the structure-soil system [41, 42, 46]. Therefore, the interfacial properties between FRP and soil are worth investigating.

The mechanical behavior between FRP and soil has been mainly studied through experiments [13, 28, 35] and numerical analysis [1] at the macroscale. Based on direct shear and interface shear tests, Frost et al. [13] reported that the interface shear behavior between FRP composites and granular materials depended on the relative roughness (surface roughness/particle mean size), the normal stress level, the initial density of the soil mass, and the angularity of the particles. Pando et al. [28] indicated that the interface friction of FRP composite piles depended on the values of relative roughness parameters through interface shear tests. A similar shear test with a modified direct shear apparatus was also used to study the interfacial friction between FRP composites and granular materials from Shaia et al. [35]. Moreover, more studies [4, 5, 51] focused on the mechanical properties of FRP and soil/sand mixes, rather than on their interface. The previous studies help us understand the mechanical behavior of FRP composite and other materials at the macroscale, but the main mechanisms at the microscale remain still unclear. Moreover, the macroscale properties of materials are significantly affected by their structural features and properties at the micro-/nano-scale where the deterioration, failure initiation, and micro-mechanical behavior of the constituents control the macroscopic behavior [15, 56]. Obtaining the mechanical properties of materials from a microscale study is the first step toward understanding the mechanical properties at a larger scale [24].

The molecular dynamics (MD) simulation method is an excellent tool for investigating the interfacial mechanical behavior of materials at the microscale, revealing their fundamental mechanisms [37, 43, 50, 54, 55]. Epoxy resin, an important FRP component, has been commonly studied in MD simulations [46, 50]. Since quartz/silica is a dominant component of the soil fraction [41], the atomistic model of epoxy-quartz/silica interface has often been established to simulate the FRP-soil interactions at the microscale [41, 46]. So far, the adhesion [50], interfacial friction [41], and creep [46] behavior of epoxy-silica interface have been recently investigated using MD simulation, but the interfacial studies between epoxy and silica at the microscale remain relatively scarce.

Nanoindentation and nanoscratching methods are effective tools for measuring the microscale elastic properties and deformation behavior of materials at a length scale of nanometers (10⁻⁹ m). The nanoindentation method has several advantages compared with other conventional laboratory tests, such as low price, less-time consuming,

non-destructive features, higher accuracy, and small sample demand. [22, 49] Previous studies of nanoindentation in MD simulations and experiments focused on the pure epoxy resin [12, 14], the mixtures epoxy composites and other modified materials (such as nanoclay, fumed silica, single-walled carbon nanotube, and graphene) [9, 20, 40], and quartz [18]. These results have shown that their nanoindentation mechanical properties were significantly affected by loading rate, temperature, content of modified material, size of the indenter tip radius, penetration depth, etc. To our knowledge, the interface mechanical behavior between epoxy and quartz has never been investigated through MD nanoindentation and nanoscratching.

In this study, to explore the interfacial properties of FRP-soil at the microscale, MD nanoindentation and nanoscratching simulations were performed to investigate the interfacial mechanical behavior between epoxy and quartz. A cross-linked epoxy resin was set as the substrate and quartz as the indenter. Uniaxial traction of the bulk epoxy was first conducted to obtain its tensile stress–strain curve, density, Young’s modulus, and Poisson’s ratio. These mechanical parameters of the bulk epoxy were compared with previous experimental and simulation results to verify the accuracy of the model. Finally, MD nanoindentation and nanoscratching simulations for epoxy-quartz system were conducted, considering the effect of loading rate, sliding velocity, and indentation depth, to explore the corresponding mechanical mechanisms. An abrasion model could be suggested depending on the indenter force Fz and indentation hardness H for nanoindentation, as well as a relationship between the sliding force Fy and the sliding distance d for nanoscratching.

# 2 Simulation procedure

## 2.1 System setup

As shown in Fig. 1a, the pyramid type quartz model with a simulation box of 35.99 Å × 33.71 Å × 34.05 Å is established via VMD software, including 1380 atoms. Figure 1b displays an epoxy resin with a simulation box of 72.4612 Å × 72.4612 Å × 55.8246 Å, which is composed of 500 diglycidyl ether of bisphenol A (DGEBA) as an organic compound commonly used as a constituent of epoxy and 250 m-phenylenediamine (mPDA) as a curing agent, where the cross-linking reaction takes place between the epoxy group and an amino group (C-O bonds broke and C-N bonds formed) at 600 K. Thereafter, the conversion percentage of 84% for the epoxy resin was obtained, which was the ratio of the number of reacted carbon atoms to the number of epoxy groups [46].

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Epoxy-quartz model: a pyramid type quartz model, b epoxy resin [46], c epoxy resin-quartz system for nanoindentation and nanoscratching simulation, where the quartz model is set as an indenter, and the epoxy as a substrate(Note that d1, d2, l1, l2, φ, θ are corresponding length and angle in the simulation box of quartz,where φ and θ are set as 27.87°and 26.35°, respectively. hc, hz, hoare indentation depth (i.e., contact depth between the epoxy substrate and the quartz indenter), displacement of quartz indenter along the z-direction, and initial distance between epoxy and quartz, respectively.)</figcaption></figure>

Figure 1c shows the epoxy resin-quartz system with 115,380 atoms composed of a pyramid type quartz model (Fig. 1a) as an indenter and a 2a × 2b × 1c supercell of epoxy resin (Fig. 1b) as a substrate. The epoxy substrate in the epoxy-quartz system was divided into three layers, a Newtonian layer, a thermostat layer, and a boundary layer, where (1) the Newtonian layer with micro-canonical (NVE) ensemble was free to move, and this Newtonian layer interacted with the quartz indenter during nanoin- dentation and nanoscratching processes. (2) The thermostat layer was set as a temperature-controlled layer to keep the whole system at a targeted temperature, by Langevin thermostat [34] and NVE ensemble for this atom group. (3) The boundary layer was fixed in its position during the whole simulation to maintain the stability of the whole system. Moreover, the quartz indenter was set as a rigid body (i.e., a single entity) during the whole simulation.

## 2.2 Force fields

### 2.2.1 PCFF force field used for epoxy

The Polymer Consistent Force Field (PCFF) force field has been successfully used to describe the interaction of atoms or molecules for organic compounds, especially polymers [6, 30, 38], and it has been successfully employed for epoxy system in previous works [38, 46]. Thus, the PCFF force field was also used for the epoxy resin system in this study, where the potential parameters of several atoms are

shown in Table S1-S5 (see Supplementary materials), and its formula was expressed as follows:

$$ E_{\text{total}} = \sum_{b} \left[ k_{2}(b-b_{0})^{2} + k_{3}(b-b_{0})^{2} + k_{4}(b-b_{0})^{2} \right] \\ + \sum_{\theta} \left[ H_{2}(\theta-\theta_{0})^{2} + H_{3}(\theta-\theta_{0})^{2} + H_{4}(\theta-\theta_{0})^{2} \right] \\ + \sum_{\phi} \left\{ V_{1}[1-\cos(\phi-\phi_{1}^{0})] \right. \\ + V_{2}[1-\cos(\phi-\phi_{2}^{0})] + V_{3}[1-\cos(\phi-\phi_{3}^{0})] \left\} \\ + \sum_{\chi} k_{\chi} \chi^{2} + \sum_{b} \sum_{b'} (b-b_{0})(b'-b_{0}') \\ + \sum_{\theta} \sum_{\theta'} (\theta-\theta_{0})(\theta'-\theta_{0}') \\ + \sum_{b} \sum_{\theta} (b-b_{0})(\theta-\theta_{0}) \\ + \sum_{b} \sum_{\phi} (b-b_{0})(V_{1}\cos\phi + V_{2}\cos2\phi + V_{3}\cos3\phi) \\ + \sum_{b'} \sum_{\phi} (b'-b_{0}')(V_{1}\cos\phi + V_{2}\cos2\phi + V_{3}\cos3\phi) \\ + \sum_{\theta} \sum_{\phi} (\theta-\theta_{0})(V_{1}\cos\phi + V_{2}\cos2\phi + V_{3}\cos3\phi) \\ + \sum_{\phi} \sum_{\theta} \sum_{\theta'} k_{\phi\theta\theta'}\cos\phi(\theta-\theta_{0})(\theta'-\theta_{0}') \\ + \frac{e^{2}}{4\pi\epsilon_{0}} \sum_{ij} \frac{q_{i}q_{j}}{r_{ij}} + \sum_{ij} \epsilon_{ij} \left[ 2\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{9} - 3\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6} \right] \quad (1) $$

where Terms 1–3 represent the bond stretching b, the angle bending (θ), and the torsion angle rotation (φ) energy,
respectively. The out-of-plane (χ) energy is defined in Term 4. Terms 5–11 are the cross-interaction terms, among which the bond-bond, bond-angle, and bond-torsion angle are the most frequently used terms. Term 12 is the Coulombic interaction between the atomic charges, where qᵢ and qⱼ are the charges of atoms i and j, ε₀ is the dielectric constant, and rᵢⱼ is the distance between them. Term 13 represents the van der Waals interactions with a 9–6 Lennard–Jones (LJ) potential function, where σᵢⱼ and εᵢⱼ are the size and energy parameters between atoms i and j, respectively. Moreover, σᵢⱼ and εᵢⱼ can be obtained by Mixing Lorentz-Berthelot’s law [11], as shown in Eqs. (2) and (3).

$$ \sigma_{ij} = \frac{\sigma_i + \sigma_j}{2} \qquad (2) $$

$$ \varepsilon_{ij} = \sqrt{\varepsilon_i \varepsilon_j} \qquad (3) $$

## 2.2.2 CVFF force field used for quartz

The Consistent Valence Force Field (CVFF) force field [17] has been successfully used for inorganic materials, such as quartz [10, 47], where their potential parameters are displayed in Table S6-S9, and its formula is given by the following equation:

$$ E_{\text{total}} = \sum_{b} D_{b} \left[ 1 - e^{-\alpha(b-b_{0})} \right] + \sum_{\theta} H_{\theta}(\theta - \theta_{0})^{2} \\ + \sum_{\phi} H_{\phi}[1 + \cos(n\phi)] + \sum_{\chi} k_{\chi}\chi^{2} \\ + \sum_{b} \sum_{b'} F_{bb'}(b-b_{0})(b'-b'_{0}) + \sum_{\theta} \sum_{\theta'} F_{\theta\theta'}(\theta - \theta_{0})(\theta' - \theta'_{0}) + \sum_{b} \sum_{\theta} F_{b\theta}(b-b_{0})(\theta - \theta_{0}) \\ + \sum_{\phi} F_{\phi\theta\theta'}\cos\phi(\theta - \theta_{0})(\theta' - \theta'_{0}) \\ + \sum_{\chi} \sum_{\chi'} F_{\chi\chi'}\chi\chi' + \frac{e^{2}}{4\pi\varepsilon_{0}} \sum_{ij} \frac{q_{i}q_{j}}{r_{ij}} \\ + \sum_{ij} 4\varepsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right] \qquad (4) $$

where Terms 1–4 represent the energy of deformation of bond lengths, bond angles, torsion angles, and out-of-plane interactions, respectively. Terms 5–9 describe the couplings between deformations of internal coordinates. For example, Term 5 corresponds to the coupling between the stretching of adjacent bonds. Terms 10–11 describe the non-bonded interactions. Term 10 is the Coulombic representation of electrostatic interactions, and Term 11 represents the van der Waals interactions with a 12–6 LJ potential function.

## 2.2.3 Hybrid force field used for the interaction between epoxy and quartz

The hybrid force field was used to describe the epoxy-quartz system, where PCFF and CVFF were for epoxy and quartz, respectively. Moreover, the standard 12–6 LJ potential was chosen for the interaction between epoxy and quartz. However, based on the above Eqs. (1) and (4), it should be noted that the terms corresponding to van der Waals interactions in PCFF and CVFF force fields are different, PCFF force field depends on 9–6 LJ potential, whereas CVFF force field depends on 12–6 LJ potential.

To determine the non-bonded interaction between epoxy and quartz due to different LJ potentials, based on the 9–6 LJ potential parameters of each atom in epoxy system, a fitting method was proposed to obtain their corresponding 12–6 LJ potential parameters. As shown in Fig. 2 and Table 1, to obtain the best fitting result, the atomic radius and the position of the first peak of RDF of atom pairs were first considered, determining the range of the starting point position. Thereafter, based on the total deviation, the effects of each atom pair at various distance intervals were compared (see Fig. 3), so that the position of the final selected fitting starting point could be determined. Finally, the corresponding new 12–6 LJ potential parameters along suitable distance intervals were obtained, as shown in Table 1.

Then, the non-bonded interaction between epoxy and quartz could be calculated by the same 12–6 LJ potential (i.e., term 10 in Eq. 4), where σᵢⱼ and εᵢⱼ between atom i and atom j were calculated by Mixing Lorentz-Berthelot’s law [11] (see Eqs. (2) and (3)). This fitting method could also be used in similar studies for materials with different LJ potentials. Furthermore, the pair interaction potential parameters between epoxy and quartz systems are shown in Table 2.

## 2.3 MD nanoindentation simulations

All MD simulations were performed using LAMMPS code [29]. Three-dimensional periodic boundary conditions and the velocity-Verlet algorithm for integrating motion equations were applied in all simulations. The cut-off radius for the van der Waals forces was 10 Å, and that for long-range electrostatic interactions using the PPPM method was 8.5 Å. The time step was 1.0 fs, and the temperature was 300 K in all simulations.

As shown in Fig. 4a, MD nanoindentation simulation for epoxy-quartz system was divided into two stages: the relaxation process and the loading process. (1) The relaxation process was performed with Langevin thermostat at 300 K and NVE ensemble (the layered setting method see Sect. 2.1) for 2 ns, the quartz indenter kept still during this
<figure><img src="image_5.png" /><figcaption>Fig. 2 Flowchart of the fitting method to obtain new potential parameters</figcaption></figure>

<table><caption>**Table 1** Based on the 9–6 LJ potential parameters of each atom in epoxy system, its 12–6 LJ potential parameters were obtained through the fitting method</caption><thead><tr><th rowspan="2">Atom type of epoxy resin</th><th colspan="2">9–6 LJ potential</th><th colspan="2">12–6 LJ potential</th><th rowspan="2">Atomic radius (Å)</th><th rowspan="2">The position of the first peak of RDF (Å)</th><th rowspan="2">The range of starting point positions for fitting (Å)</th><th rowspan="2">The position of the final selected fitting starting point (Å)</th></tr><tr><th>ε (kcal/mol)</th><th>σ (Å)</th><th>ε (kcal/mol)</th><th>σ (Å)</th></tr></thead><tbody><tr><td>c3</td><td>0.054</td><td>4.01</td><td>0.056</td><td>3.479</td><td>0.7</td><td>4.05</td><td>1.4 ~ 4.05</td><td>3.3</td></tr><tr><td>hc</td><td>0.02</td><td>2.995</td><td>0.020</td><td>2.589</td><td>0.25</td><td>2.45</td><td>0.5 ~ 2.45</td><td>2.4</td></tr><tr><td>c</td><td>0.054</td><td>4.01</td><td>0.056</td><td>3.479</td><td>0.7</td><td>5.65</td><td>1.4 ~ 4.05</td><td>3.3</td></tr><tr><td>cp</td><td>0.064</td><td>4.01</td><td>0.066</td><td>3.479</td><td>0.7</td><td>6.05</td><td>1.4 ~ 6.05</td><td>3.3</td></tr><tr><td>oc</td><td>0.24</td><td>3.535</td><td>0.246</td><td>3.065</td><td>0.6</td><td>4.85</td><td>1.2 ~ 4.85</td><td>2.9</td></tr><tr><td>c2</td><td>0.054</td><td>4.01</td><td>0.056</td><td>3.479</td><td>0.7</td><td>4.35</td><td>1.4 ~ 4.05</td><td>3.3</td></tr><tr><td>c1</td><td>0.054</td><td>4.01</td><td>0.056</td><td>3.479</td><td>0.7</td><td>4.45</td><td>1.4 ~ 4.05</td><td>3.3</td></tr><tr><td>oh</td><td>0.24</td><td>3.535</td><td>0.246</td><td>3.065</td><td>0.6</td><td>2.85</td><td>1.2 ~ 4.85</td><td>2.9</td></tr><tr><td>c3h</td><td>0.054</td><td>4.01</td><td>0.056</td><td>3.479</td><td>0.7</td><td>3.75</td><td>1.4 ~ 4.05</td><td>3.3</td></tr><tr><td>o3e</td><td>0.24</td><td>3.535</td><td>0.246</td><td>3.065</td><td>0.6</td><td>4.35</td><td>1.2 ~ 4.85</td><td>2.9</td></tr><tr><td>ho</td><td>0.013</td><td>1.098</td><td>0.013</td><td>0.949</td><td>0.25</td><td>2.35</td><td>0.5 ~ 2.35</td><td>0.9</td></tr><tr><td>nb</td><td>0.065</td><td>4.07</td><td>0.065</td><td>3.524</td><td>0.65</td><td>3.55</td><td>1.3 ~ 3.55</td><td>3.3</td></tr><tr><td>hn</td><td>0.013</td><td>1.098</td><td>0.013</td><td>0.949</td><td>0.25</td><td>2.45</td><td>0.5 ~ 2.35</td><td>0.9</td></tr></tbody></table>

process. Thereafter, (2) there are usually two loading modes for the loading process of nanoindentation: displacement-controlled mode and load-controlled mode. The displacement-controlled mode was used during the loading process in this work, where the displacement of the quartz indenter was controlled by the constant displacement rate (i.e., loading rate), which was 0.1 ~ 2.0 Å/1000 steps (i.e., 0.1 ~ 2.0 Å/ps, based on the timestep of 1.0 fs). The quartz indenter was initially located 20 Å ($h_o$) above the epoxy substrate surface (see Fig. 1c), and then gradually approached and indented into the epoxy substrate during the loading process. The final displacement ($h_z$) of the quartz indenter along z-direction was set at 45 Å, and the corresponding final indentation depth ($h_c$) was 25 Å.

Moreover, the indentation depth (i.e., contact depth, $h_c$) was calculated by the following equation:

$$h_c = h_z - h_0 = \Delta h \times n_{\text{step}}^{\text{loading}} - h_0 \quad (5)$$

where $h_c$, $h_z$, $h_o$ are the indentation depth (corresponding to the depth of indentation into the epoxy resin substrate), the displacement of the quartz indenter along z-direction, and the initial distance between epoxy and quartz ($h_o = 20$ Å), respectively. $n_{\text{step}}^{\text{loading}}$ is the number of loading steps, where the timestep is 1 fs. $\Delta h$ is the displacement rate (i.e., loading rate) of the quartz indenter, which is equal to 0.1 ~ 2.0 Å/1000 steps (i.e., 0.1 ~ 2.0 Å/ps = 10 ~ 200 m/s). Note that in the actual nanoindentation experiments, the loading rate is in the order of µm/s (10⁻⁶ ~ 10⁻⁹ m/s) instead of m/s [23]. However, in MD nanoindentation simulations, the loading rate is in the order of m/s (0.5 ~ 400 m/s for various materials [22]) and is
<figure><img src="image_6.png" /><figcaption>Fig. 3 Evolution of van der Waals potential energy of C–C atom pairs with its distance based on different distance intervals: a 1.6 ~ 12.5 Å; b 2.0 ~ 12.5 Å; c 2.6 ~ 12.5 Å; d 2.8 ~ 12.5 Å; e 3.0 ~ 12.5 Å; f 3.2 ~ 12.5 Å; g 3.3 ~ 12.5 Å; h 3.4 ~ 12.5 Å. (The blue line refers to the original 9–6 LJ potential, and the orange line is the 12–6 LJ potential obtained from the fitting method. Note that the van der Waals potential energy of other atom pairs (such as H–H, O–O, N–N, etc.) were also analyzed, their fitting curves are not shown here to reduce the length of this article.)</figcaption></figure>

usually used due to the limitation of computational
capacity [16].

The indentation hardness (H) represents the ability of
the tested material to withstand contact loads [22], and is
defined as the ratio of the maximum indentation force
(P<sub>max</sub>) to the projected contact area (A<sub>c</sub>) of the indenter as
follows:

$$
d_1 = 2h_c \times \tan\varphi \qquad (6)
$$

$$
d_2 = 2h_c \times \tan\theta \qquad (7)
$$

$$
A_c = d_2 \times d_1 = 4h_c^2 \times \tan\varphi \times \tan\theta \quad (8)
$$

$$
H = \frac{P_{\max}}{A_c} \tag{9}
$$

where A_c is the projected contact area of the indenter at the
indentation depth of h_c. H and P_max are the indentation
hardness of epoxy substrate and maximum indentation
force at maximum indentation depth, respectively. More-
over, other physical quantities are introduced in the caption
of Fig. 1.

The reduced modulus (alias indentation modulus, E<sub>r</sub>) is
applicable to various axisymmetric shapes of the indenter
and elastic-plastic contact, where the relationship between
<table><caption>**Table 2** Epoxy-quartz pair interactions potential parameters</caption><thead><tr><th rowspan="2">Model</th><th>lj/cut (12–6 LJ potential)</th><th rowspan="2">Symbol</th><th rowspan="2">ε<br>(kcal/mol)</th><th rowspan="2">σ (Å)</th></tr><tr><th>Atom species</th></tr></thead><tbody><tr><td rowspan="12">Epoxy</td><td>Carbon in the methyl group</td><td>c3</td><td>0.056</td><td>3.479</td></tr><tr><td>Hydrogen, bonded to carbon</td><td>hc</td><td>0.02</td><td>2.589</td></tr><tr><td>Generic carbon</td><td>c</td><td>0.056</td><td>3.479</td></tr><tr><td>Aromatic carbon</td><td>cp</td><td>0.066</td><td>3.479</td></tr><tr><td>Oxygen in ether</td><td>oc</td><td>0.246</td><td>3.065</td></tr><tr><td>Carbon bonded to 2 H's</td><td>c2</td><td>0.056</td><td>3.479</td></tr><tr><td>Carbon bonded to 1 H</td><td>c1</td><td>0.056</td><td>3.479</td></tr><tr><td>Hydroxyl oxygen</td><td>oh</td><td>0.246</td><td>3.065</td></tr><tr><td>Carbon in 3-membered ring<br>with hydrogens</td><td>c3h</td><td>0.056</td><td>3.479</td></tr><tr><td>Oxygen in 3-membered ring</td><td>o3e</td><td>0.246</td><td>3.065</td></tr><tr><td>Hydroxyl hydrogen</td><td>ho</td><td>0.013</td><td>0.949</td></tr><tr><td>Nitrogen in aromatic amines</td><td>nb</td><td>0.065</td><td>3.524</td></tr><tr><td>Hydrogen bonded to nitrogen</td><td>hn</td><td>0.013</td><td>0.949</td></tr><tr><td rowspan="5">Quartz</td><td>Silicon in silicate</td><td>sz</td><td>0.040</td><td>4.053</td></tr><tr><td>Hydroxyl oxygen</td><td>oh</td><td>0.228</td><td>2.860</td></tr><tr><td>Hydroxyl hydrogen</td><td>ho</td><td>0.000</td><td>0.000</td></tr><tr><td>Oxygen</td><td>oz</td><td>0.228</td><td>2.860</td></tr></tbody></table>

Young's modulus (E) of the substrate epoxy and $E_r$ is expressed as follows [23]:

$$ \frac{1}{E_r} = \frac{1 - v^2}{E} + \frac{1 - v_i^2}{E_i} \quad (10) $$

where E and v are Young's modulus and Poisson's ratio of the substrate material (epoxy in this work), and $E_i$, $v_i$ are Young's modulus and Poisson's ratio of the indenter material (quartz in this work). The quartz indenter being set as a rigid body, $E_i$ and $v_i$ in Eq. (10) are equal to $\infty$ and 0,

respectively, and Eq. (10) could be further expressed as follows:

$$ \frac{1}{E_r} = \frac{1 - v^2}{E} \quad (11) $$

To quantify the plastic deformation behavior of atoms in the whole system, the shear strain (i.e., von mises strain, $\eta_i^{\text{Mises}}$) and volumetric strain ($\eta_i^{\text{Vol}}$) of atoms were calculated using OVITO from the atomic Green-Lagrangian strain tensor ($\eta$) and the atomic deformation gradient tensor ($J$) for each atom, as shown in Eqs. (12–14) [36]. Note that this atomic-level deformation gradient and the strain tensor at each particle could be calculated from the relative motion of its neighbors.

$$ \eta = \frac{J^T J - I}{2} \quad (12) $$

$$ \eta_i^{\text{Mises}} = \sqrt{\eta_{xy}^2 + \eta_{xz}^2 + \eta_{yz}^2 + \frac{(\eta_{xx} - \eta_{yy})^2 + (\eta_{xx} - \eta_{zz})^2 + (\eta_{yy} - \eta_{zz})^2}{6}} \quad (13) $$

$$ \eta_i^{\text{Vol}} = \frac{\eta_{xx} + \eta_{yy} + \eta_{zz}}{3} \quad (14) $$

where $\eta$ and $J$ are the atomic Green-Lagrangian strain tensor and atomic deformation gradient tensor for each particle, respectively. $\eta_i^{\text{Mises}}$ and $\eta_i^{\text{Vol}}$ are the shear strain and volumetric strain of the atom. $\eta_{xx}$, $\eta_{yy}$, $\eta_{zz}$, $\eta_{xy}$, $\eta_{xz}$, and $\eta_{yz}$ are the six components of strain for each atom.

The displacement magnitude of atoms during nanoindentation and nanoscratching processes is computed by subtracting its reference (or initial) position from its current position, as shown in Figures 6d, 8a, and 14a.

<figure><img src="image_7.png" /><figcaption>**Fig. 4** Epoxy-quartz system for a nanoindentation and b nanoscratching simulation</figcaption></figure>
## 2.4 MD nanoscratching simulations

As shown in Fig. 4b, MD nanoscratching simulation was also divided into two stages: the holding process and the sliding process. (1) after MD nanoindentation simulation, the holding process was conducted for 100 ps, where the quartz indenter was held in the target indentation depth without displacement. Five target indentation depths ($h_c$: 5, 10, 15, 20, and 25 Å) were considered. Thereafter, (2) during the sliding process, the quartz indenter was slid along y-direction in displacement-controlled mode at a given velocity. Five sliding velocities (0.1, 0.5, 1.0, 1.5, and 2.0 Å/ps) were applied, and the final sliding distance of the quartz indenter along y-direction was set at 30 Å.

# 3 Results and discussion

## 3.1 Uniaxial traction on bulk epoxy

The bulk epoxy system with 114,000 atoms (see Fig. 5d) was further relaxed with NPT ensemble for 1.0 ns, with the temperature and pressure equal to 300 K and 1.0 atm, respectively. Thereafter, the tensile mechanical properties of epoxy system along three directions were evaluated by applying uniaxial traction at 300 K, under a strain-controlled loading approach [52] at a tensile strain rate of $5 \times 10^{-7}/\text{fs}$.

Figure 5a–c shows the evolution of the tensile stress, total energy, and density of the bulk epoxy system with a tensile strain along three directions. The calculation method of the tensile stress, tensile strain, and Poisson's ratio (v) in MD simulation using LAMMPS can be found in Eqs. S3–S5 (see Supplementary materials).

As shown in Fig. 5a, compared to the mechanical behavior of crystal materials in MD simulations, such as clay minerals [44, 48, 52], the tensile stress–strain curve of the bulk epoxy presents a high discreteness without an obvious yield failure point, which agrees well with the MD simulation results from Xu et al. [46]. This can be explained by the fact that the epoxy resin, as an amorphous material, did not present obvious pores during the whole tensile process in the simulation box (Fig. 5d). The density of the bulk epoxy system decreased from 1.13 to around $1.10 \text{ kg/m}^3$ at a tensile strain equal to 0.3, demonstrating that the difference in density ($\sim 0.03 \text{ kg/m}^3$) was very small, as shown in Fig. 5c. Moreover, Fig. 5b shows that the total energy of the whole system increased with the tensile strain.

Tables 3 and 4 give Young's modulus and Poisson's ratio of the epoxy resin, which was calculated by linear regression of the tensile stress–strain curve for a tensile strain of 0 ~ 0.04 (Fig. 5a). As shown in Tables 3 and 4, the overall Young's modulus ($\bar{E}$) of 2.075 GPa and the average Poisson's ratio ($\bar{\nu}$) of 0.377 are both within the range of previous experimental and simulation results, indicating that the model built in this work appears well

<figure><img src="image_8.png" /><figcaption>**Fig. 5** a–c Evolution of tensile stress, total energy, and density of the bulk epoxy system with tensile strain along three directions. **d** Epoxy structure and its displacement magnitude under different tensile strains along x-direction</figcaption></figure>
<table><caption>**Table 3** Young's modulus of epoxy resin in this work and similar materials in previous studies at room temperature</caption><thead><tr><th>Methods</th><th>Young's modulus of epoxy (GPa)</th><th>E<sub>x</sub></th><th>E<sub>y</sub></th><th>E<sub>z</sub></th><th>Overall (Ē)</th><th>Reduced modulus E<sub>r</sub></th></tr></thead><tbody><tr><td rowspan="4">MD simulations</td><td>Epoxy resin (containing 2000 DGEBA + 1000 mPDA) with a conversion percentage of 84% through uniaxial tension using MD simulation in this work</td><td>2.203</td><td>1.906</td><td>2.105</td><td>2.075</td><td>2.419</td></tr><tr><td>Epoxy resin (containing 200 DGEBA + 100 mPDA) with a conversion percentage of 86% through uniaxial tension using MD simulation [46]</td><td>2.61 ± 0.28</td><td></td><td></td><td></td><td></td></tr><tr><td>Epoxy resin composed of EPON862 and curing agent DETDA with conversion degrees of 0% ~ 86% using MD simulation [21]</td><td>1.0 ~ 3.0</td><td></td><td></td><td></td><td></td></tr><tr><td>Epoxy resin (containing EPON862 + DETDA) [39]</td><td>2.56 (experimental), 1.61 (MD simulation)</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="7">Various experiments in previous works</td><td>DGEBA epoxy networks (containing DGEBA/D230, DGEBA/TETA, DGEBA/TEPA) through dynamic mechanical analyses experiments [3]</td><td>2.21 ~ 2.72</td><td></td><td></td><td></td><td></td></tr><tr><td>A DGEBA-based epoxy resin through tensile test [51]</td><td>3.5 (unmodified epoxy resin), 3.85 (epoxy + 8.5% nanofiller)</td><td></td><td></td><td></td><td></td></tr><tr><td>Matrix epoxy Epon 862 [25]</td><td>2.83</td><td></td><td></td><td></td><td></td></tr><tr><td>Epoxy resin composite by ultrasonic velocity Measurement [27]</td><td>3.441 ~ 4.601 (DGEBA/AGE), 3.441 ~ 6.293 (DGEBA/GMA)</td><td></td><td></td><td></td><td></td></tr><tr><td>Epoxy resin through uniaxial tension and nanoindentation test [40]</td><td>3.3 ± 0.12 (uniaxial tension), 3.9 ± 0.12 (nanoindentation test)</td><td></td><td></td><td></td><td></td></tr><tr><td>Epoxy resin with 0.78 wt% Polyaniline added by 0 ~ 1.0 wt% ZnO nanoparticles using nanoindentation test [31]</td><td>3.76 (0 wt%), 3.49 (0.5 wt%), 3.47 (1.0 wt%)</td><td></td><td></td><td></td><td></td></tr><tr><td>Epoxy composite with 0 ~ 5% CNF using nanoindentation test [32]</td><td>4.0 ~ 7.1</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>E<sub>r</sub> of epoxy resin at different loading rates of 0.25 ~ 6.0 mN/s using nanoindentation test [14]</td><td>4.78 ~ 5.35</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>E<sub>r</sub> of nanoclay/epoxy composite at 40 °C [20]</td><td>1.5 ~ 7.5</td><td></td><td></td><td></td><td></td></tr></tbody></table>

The overall Young's modulus (Ē) is the square average value of Young's modulus along three directions (i.e., E<sub>x</sub>, E<sub>y</sub>, E<sub>z</sub>). The reduced modulus E<sub>r</sub> was calculated through Eq. (10)

DGEBA Diglycidyl Ether of Bisphenol A, mPDA m-phenylenediamine, DETDA Diethyltoluenediamine, EPON862 Bisphenol F diglycidyl ether, TETA triethylenetetramine, TEPA tetraethylenepentamine, D230 Jeffamine D230, AGE Allyl glycidyl ether, GMA glycidyl methacrylate, CNF carbon nanofiber

suited for modeling the mechanical behavior of epoxy resin. Note that the data in Tables 3 and 4 were not the same, because the epoxy resin composition was diverse with different matrixes and curing agents, and various test conditions and methods were used in the experiments. In MD simulation, the sizes of the epoxy model, the conversion percentage of cross-linking, the equilibrium process, etc., were also somewhat different.

Moreover, as shown in Fig. 5, Tables 3, and 4, the tensile stress, density, Young's modulus, and Poisson's ratio along three directions for the bulk epoxy system are relatively close, indicating that the epoxy resin can be considered as an isotropic material [8, 53].

## 3.2 MD nanoindentation behavior

Figure 6(a, b) shows the evolution of the indentation force (F<sub>z</sub>) with the displacement of the quartz indenter (h<sub>z</sub>) and the indentation depth (h<sub>c</sub>) at different loading rates. The indentation force gradually increased with h<sub>z</sub> and the loading rate. As shown in Fig. 6b, the relationship between

F<sub>z</sub> and h<sub>c</sub> follows a power law, which has been proven to be suitable for most cases in nanoindentation via MD and experiments [23, 26]:

$$F_z = \alpha h_c^m \quad (15)$$

where F<sub>z</sub> and h<sub>c</sub> are the indentation force of the quartz indenter and the indentation depth, respectively. α and m are constants.

During MD nanoindentation, the epoxy-quartz structure and the surface mesh of the epoxy at different h<sub>z</sub> are shown in Figs. 6b and 7, respectively. The quartz indenter was pressed into the interior of the epoxy substrate, squeezing the atoms on the epoxy substrate surface and forming significant pits. The surface wear of the epoxy substrate during MD indentation was an obvious “groove.” Figure 8 shows that the displacement magnitude, shear strain, and volumetric strain of atoms within the whole system at a loading rate of 1.0 Å/ps constantly increased with the increase of h<sub>z</sub>, and the surface wear (i.e., the deformation area) in the whole system gradually expanded when h<sub>z</sub> increased. The higher h<sub>z</sub>, the more atoms of the epoxy
<table><caption>**Table 4** Poisson's ratio of epoxy resin in this work and similar materials in previous studies at room temperature</caption><thead><tr><th>Poisson's ratio of epoxy</th><th>v<sub>xy</sub></th><th>v<sub>xz</sub></th><th>v<sub>yx</sub></th><th>v<sub>yz</sub></th><th>v<sub>zx</sub></th><th>v<sub>zy</sub></th><th>v&#773;</th></tr></thead><tbody><tr><td>Uniaxial tension using MD simulation in this work</td><td>0.351</td><td>0.407</td><td>0.338</td><td>0.420</td><td>0.368</td><td>0.378</td><td>0.377</td></tr><tr><td>Epoxy resin composed of EPON862 and curing agent DETDA (conversion degrees of 0% ~ 86%) [21]</td><td colspan="7">0.27 ~ 0.475</td></tr><tr><td>Epoxy resin (containing EPON862 + DETDA) [39]</td><td colspan="7">0.353 (experimental), 0.448 (MD simulation)</td></tr><tr><td>Matrix epoxy Epon 862 [25]</td><td colspan="7">0.3</td></tr><tr><td>Epoxy resin composite by ultrasonic velocity Measurement [27]</td><td colspan="7">0.3461 ~ 0.3911 (DGEBA/AGE), 0.341 ~ 0.346 (DGEBA/GMA)</td></tr><tr><td>Araldite LY5052 epoxy resin by uniaxial relaxation tests [33]</td><td colspan="7">0.32 ~ 0.44</td></tr></tbody></table>

The average Poisson's ratio ($\bar{\nu}$) is the geometric average value of the six components of the matrix

substrate were displaced by the quartz indenter, and the greater the resistance the indenter had to overcome. This is why the indentation force gradually increased with $h_z$.

The indentation force had a significant dispersion at a loading rate of 0.1 Å/ps (see Fig. 6(a, b)), and the fitting coefficient $R^2$ was 0.78, which was significantly lower than the ones obtained for the other loading rates (see Fig. 6b). Thus, the loading rate of 0.1 Å/ps was not considered to calculate the indentation hardness of the epoxy substrate. As shown in Table 5, the indentation hardness of the epoxy substrate obtained by Eq. (9) increased with the increase in the loading rate from 0.25 to 0.75 Å/ps, and then tended to stabilize in the interval 0.75 ~ 2.0 Å/ps. These results agree well with some previous nanoindentation experimental works [12, 14], which reported that the hardness of epoxy resin increased with the loading rate.

As shown in Table 3, the value of the reduced modulus $E_r$ of 2.419 GPa is consistent with previous experimental results [20]. However, as shown in Table 5, the indentation hardness of the epoxy resin in MD simulations appears to be one order of magnitude larger than the one obtained experimentally, which is fully discussed in Sect. 3.4.1.

To understand the interaction between quartz and epoxy during nanoindentation, Fig. 9a shows the pairwise interaction force between quartz and epoxy groups with $h_z$ at a loading rate of 1.0 Å/ps. At $h_z$ of 15 ~ 20 Å, the pairwise interaction between quartz and epoxy groups originated from their non-bonded interactions, including coulombic and van der Waals interactions. Furthermore, the pairwise interaction force along z-direction increased significantly with rising $h_z$ at 20 ~ 45 Å, while the pairwise interaction force along x- and y-directions fluctuated around 0 nN, and the fluctuation increased with increasing $h_z$, especially at 35 ~ 45 Å.

Figure 10a-c shows the evolution of various energy terms with $h_z$. The potential energy, kinetic energy, and total energy of the whole system gradually increased with $h_z$. During the process of the quartz indenter entering into the epoxy substrate, some atoms within the epoxy substrate could move and rearrange, increasing various energy terms

(Fig. 10b, c). At the same time, this MD nanoindentation process required overcoming the attraction between atoms or molecules, thus requiring energy consumption. Therefore, an absorbing energy process is performed during MD nanoindentation.

Figures 10a, d, 6b, 9, and Table 5 illustrate that the loading rate could affect the energy of the whole system and the mechanical parameters of nanoindentation: the higher the loading rate, the lower all the energy terms, the higher the indented force, the indentation hardness of the epoxy, and the average overall pairwise interaction force between quartz and epoxy groups.

To explore the underlying mechanisms of the above phenomenon, Fig. 11 presents the MD nanoindentation scheme: the initial model (Fig. 1c) would gradually deform with the increase of $h_z$ when the quartz indenter constantly moves into the epoxy substrate (see Figs. 7 and 8). Then, to try to reach a new equilibrium state, the movement and rearrangement of atoms in the whole system occurred, but at a high loading rate, there was not enough time to complete this equilibrium process at a given timestep, leading to a smaller movement and displacement of atoms in the epoxy corresponding to a non-equilibrium state (see Fig. 6d). Thus, at a higher loading rate, the distance of the surface connecting atoms between quartz and epoxy was smaller, causing $F_z$ and $H$ to be higher. Moreover, at a higher loading rate, more atoms in the whole system could not undergo sufficient movement and rearrangement in time, so the change of potential and kinetic energy of the whole system was relatively smaller. Thus, the higher the loading rate, the lower the total energy of the whole system.

### 3.3 MD nanoscratching behavior

Several MD nanoscratching simulations along y-direction (i.e., sliding direction) were performed with five indentation depths ($h_c$: 5, 10, 15, 20, and 25 Å) and five sliding velocities (0.1, 0.5, 1.0, 1.5, and 2.0 Å/ps).
<figure><img src="image_9.png" /><figcaption>**Fig. 6** a, b Evolution of the indentation force ($F_z$) of the whole system with the displacement of the indenter ($h_z$) and the indentation depth ($h_c$) at different loading rates. c Epoxy-quartz structure at different $h_z$ and $h_c$. d Displacement magnitude of atoms within the whole system at different loading rates and $h_c = 25$ Å</figcaption></figure>

To explore the effect of indentation depths on the nanoscratching behavior, Fig. 12(a-e) presents the evolution of the forces on the quartz indenter along three directions with the sliding distance for different $h_c$ and a sliding velocity of 1.0 Å/ps. As shown in Fig. 12a, for the case of $h_c = 5$ Å, the forces along the three directions fluctuated around a stable value, and this fluctuation was relatively small. However, as shown in Fig. 12b-e, the fluctuation of the three forces gradually increased when $h_c$

increased. With the growth of $h_c$, the contact area between the quartz indenter and the epoxy gradually rose, indicating that the quartz indenter needed to overcome more interaction between atoms, resulting in a higher resistance during nanoscratching process. Thus, the higher the $h_c$, the higher the shear strain of atoms and the surface wear of the epoxy substrate (Fig. 13), and the higher the total energy of the whole system (Fig. 12f).
<figure><img src="image_10.png" /><figcaption>Fig. 7 Surface mesh of the epoxy system under different displacements of the indenter ($h_z$) at a loading rate of 1.0 Å/ps (the indenter is hidden for better visual clarity)</figcaption></figure>

As shown in Fig. 12d–f, the three forces and the total energy of the whole system gradually increased with the sliding distance at $h_c$ equal to 20 and 25 Å. Figure 14 displays the deformation pattern of the epoxy system under different sliding distances at $h_c$ equal to 20 Å and the sliding velocity equal to 1.0 Å/ps. The surface wear gradually expanded with the increase of the sliding distance, showing higher displacement, shear and volumetric strain for more atoms within the epoxy. As shown in Fig. 14a, with the increase in the sliding distance, more atoms were packed in front of the indenter during the nanoscratching process, causing more resistance to slide. Meanwhile, there was a fan-shaped distribution of the displacement magnitude and surface mesh of atoms (see Fig. 14a, b). Thus, for $h_c$ of 20 and 25 Å, the higher the sliding distance, the higher the sliding force $F_y$.

Figure 12f shows that the total energy of the whole system increases with rising the sliding distance, illustrating that MD nanoscratching is also an energy absorption process, and the higher the sliding distance, the higher the total energy.

To study the effect of the sliding velocity on nanoscratching behavior, Fig. 15a–e shows the evolution of the forces along three directions with the sliding distance at five sliding velocities from 0.1 to 2.0 Å/ps at $h_c$ of 20 Å. It can be concluded that the higher the sliding velocity, the higher the forces along the three directions and their fluctuation, especially the force along y-direction. Moreover, as shown in Fig. 16, based on the shear strain and displacement magnitude of atoms, the sliding velocity has little effect on the surface wear of the epoxy substrate during nanoscratching.

During the nanoscratching process, the forces along the sliding direction ($F_y$) are the main forces. Figure 15f shows the evolution of $F_y$ with the sliding distance at 0.1 ~ 2.0 Å/ps and $h_c$ of 20 Å. A relationship between the

sliding force and the sliding distance could be established in the form of an exponential function (see Eq. 16), but the coefficient of correlation (R²) of 0.42 ~ 0.63 is not very high due to the high dispersivity of $F_y$ values.

$$F_y = ae^{(-d/b)} + c \quad (16)$$

where $F_y$ and $d$ are the sliding force and sliding distance, respectively. $a$, $b$, and $c$ are constants.

Figure 17a–d shows the evolution of the average forces along three directions with sliding velocity at different $h_c$, where the forces are averaged along the sliding distance between 15 and 30 Å (see Figs. 12a–e and 15a–e). As shown in Fig. 17, the average forces along y- and z-directions are significantly higher than that along x-direction, since the indenter was sliding along y-direction, causing the influence of wear in x-direction to be relatively small. Moreover, the higher the sliding velocity or indentation depths, the higher the force along the three directions. Based on the above discussion, the MD nanoscratching mechanism diagram could be obtained, as shown in Fig. 18.

## 3.4 Discussion

According to the above research results, the nanoindentation and nanoscratching mechanisms of epoxy-quartz interface are deeply discussed, considering the effect of loading rate, sliding velocity, and indentation depth. Moreover, there are some limitations and prospects of this study, which are discussed as follows:

### 3.4.1 Indentation hardness of epoxy in MD simulation and experiments

As shown in Table 5, although the indentation hardness of epoxy resin in MD simulation is one order of magnitude
<figure><img src="image_11.png" /><figcaption>**Fig. 8** a Displacement magnitude of atoms for the whole system at different displacements of the indenter ($h_z$) under two different slices: one slice along x-direction from 0 to 74 Å, and another one from 64 to 74 Å. b Shear strain and c volumetric strain of atoms at different depths under the slice along x-direction from 64 to 74 Å. (The above figures refer to the case of a loading rate of 1.0 Å/ps. The indenter is hidden for better visual clarity.)</figcaption></figure>

greater than that obtained from previous nanoindentation experiments, it does not mean that the MD simulation results are unreasonable.

On the one hand, this is due to the difference of nanoindenter materials, such as rigid quartz as nanoindenter in this MD simulations and diamond (its chemical formula: C) as nanoindenter in conventional nanoindentation experiments [12]. Since diamond is very hard, but the quartz indenter is set as a rigid body in this MD simulation, the influence of the nanoindenter hardness could be eliminated. However, the molecular interface interactions between diamond nanoindenter and epoxy are significantly different from those between quartz nanoindenter and epoxy, so their friction coefficients and adhesion effect are

also different, which will affect the mechanical response during the indentation process. Meanwhile, the defects of epoxy resin in experiment are also an important reason, which could cause the decrease of hardness. In MD simulations, the defects in epoxy resins are usually reduced due to size effects.

On the other hand, this may be due to the type of applied force field, which is a more critical influencing factor because the potential function parameters of materials significantly affect calculating the maximum indentation force ($P_{max}$). In this work, PCFF and CVFF forcefield are performed for epoxy and quartz, respectively. This PCFF-CVFF force field combination has been successfully used to study the friction and creep behavior of epoxy-quartz
<table><caption>**Table 5** The indentation hardness (H) of the epoxy resin obtained by Eq. (9) in MD nanoindentation simulation at different loading rates and in previous experimental works</caption><thead><tr><th rowspan="2">Methods</th><th colspan="2">Conditions and other notes</th><th rowspan="2">H (GPa)</th></tr><tr><th>Loading rate (Å/ps)</th><th>P<sub>max</sub> (nN)</th></tr></thead><tbody><tr><td rowspan="6">MD nanoindentation simulation in this work, where nanoindenter is rigid quartz</td><td>0.25</td><td>9.757</td><td>1.492</td></tr><tr><td>0.5</td><td>10.594</td><td>1.620</td></tr><tr><td>0.75</td><td>13.052</td><td>1.997</td></tr><tr><td>1.0</td><td>13.115</td><td>2.006</td></tr><tr><td>1.5</td><td>14.813</td><td>2.266</td></tr><tr><td>2.0</td><td>14.215</td><td>2.174</td></tr><tr><td rowspan="10">Nanoindentation experiments in previous works, where nanoindenter is usually diamond</td><td>Epoxy-reinforced nanocomposites with fumed silica of 0 ~ 3% [40]</td><td></td><td>0.256 ~ 0.318</td></tr><tr><td>RTM6 Epoxy resin with strain rate of 0.01 ~ 0.25 and 25 °C [12]</td><td></td><td>0.23 ~ 0.35</td></tr><tr><td>SWNT of 0% ~ 1.0% reinforced epoxy composites [9]</td><td></td><td>0.23 ~ 0.31</td></tr><tr><td>Epoxy resin at different loading rate of 0.25 ~ 6.0 mN/s [14]</td><td></td><td>0.26 ~ 0.276</td></tr><tr><td>Epoxy resin with 0.78 wt% Polyaniline added by 0 ~ 1.0 wt% ZnO nanoparticles [31]</td><td></td><td>0.18 (0 wt%) ~ 0.206 (1.0 wt%)</td></tr><tr><td>Epoxy composite with 0 ~ 5% CNF [32]</td><td></td><td>0.26 ~ 0.67</td></tr></tbody></table>

RTM resin transfer molding, SWNT single-walled carbon nanotube, CNF: carbon nanofiber

systems [46]. Moreover, other force field combinations have been used to investigate the physio-chemical properties of epoxy-quartz systems, such as CVFF-CVFF [41, 50] and ReaxFF [45]. So far, no other MD nanoindentation simulation of epoxy has been performed, which prevents us from comparing our simulations' results.

### 3.4.2 Scale difference between MD nanoindentation simulation and experiments

The macroscopic results often reflect the overall mechanical properties of the epoxy resin and quartz mixture, but the microscopic nanoindentation and nanoscratching studies can better observe the interaction and deformation characteristics of the epoxy resin-quartz interface, reflecting their local deformation behavior, such as the strength, toughness, and indentation response of the interface.

It is noted that there is a scale difference (i.e., time and space scales) between MD nanoindentation simulation and experiments [23]. For example, in MD nanoindentation simulations, its time and space scales are ~ ps to ~ ns, and ~ nm for nanoindenter, respectively. In nanoindentation experiments, its time and space scales are ~ s, and ~ ten-hundreds nm for the standard nanoindenter, respectively. Moreover, due to the computational limitations, the loading rate of nanoindenters is ~ Å/ps (approximate ~ m/s) for MD simulations, but ~ µm/s (10⁻⁶ ~ 10⁻⁹ m/s) for nanoindentation experiments.

### 3.4.3 Elastic-plastic deformation analysis at macro-micro scales

#### 1. Elastic-plastic deformation of epoxy resin in MD nanoindentation

The MD nanoindentation process of epoxy resin can usually be analyzed by the relationship between indentation depth and load, which is mainly divided into two stages: elastic and plastic deformation stages (Fig. 6b).

* Elastic deformation stage: When the indenter initially contacts the epoxy resin and begins to apply pressure, exhibiting elastic deformation. The load-depth curve in this stage usually shows a linear relationship, which follows Hooke's law.

* Plastic deformation stage: As the indentation depth increases, the epoxy resin begins to show plastic deformation characteristics. Meanwhile, the load-depth curve deviates from linearity and the hysteresis occurs. Plastic deformation involves irreversible processes, such as sliding and deformation of molecular chains, and breaking or reorganization of crosslinking points.

#### 2. Elasticity of intermolecular interaction potential and permanent deformation of materials in MD

The interaction potentials in MD simulations (such as Lennard-Jones potential, bond length and bond angle potential) behave elastically within a small deformation range. Because these potential energy functions are parabolic near the equilibrium position,
<figure><img src="image_12.png" /><figcaption>**Fig. 9** a Evolution of pairwise interaction force between quartz and epoxy groups with $h_z$ at a loading rate of 1.0 Å/ps, where the overall pairwise interaction force is the square average of the forces in the three directions. The right sketch map shows the quartz-epoxy system containing the two groups and the pairwise interaction energy/force between them. b Relationship between the average overall pairwise interaction force between the two groups and the loading rate at two different intervals of $h_z$, i.e., 40 ~ 45 Å and 20 ~ 45 Å</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10** a Relationship between total energy of the whole system and $h_z$ at different loading rates. b, c Evolution of energy increment of energy terms with $h_z$ at a loading rate of 1.0 Å/ps. d Energy difference ($\Delta E$) between the cases of $h_z = 45$ Å and $h_z = 0$ Å at different loading rates, including potential energy (pe), kinetic energy (ke), and total energy (total)</figcaption></figure>

following Hooke's law.

Although those interaction potentials itself is elastic, the microstructure inside the material will undergo

plastic deformation under large deformation or high stress, resulting in permanent deformation. These irreversible processes are beyond the scope of simple
<figure><img src="image_14.png" /><figcaption>**Fig. 11** MD nanoindentation mechanism diagram. *d*<sub>A-B</sub>, *F*<sub>z</sub>, and *H* are the distance of the surface connecting atoms between quartz and epoxy, the indenter force, and the indentation hardness of the epoxy substrate, respectively</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12** Evolution of *a–e* forces along three directions, *f* total energy, *g* epoxy-quartz structure and shear strain of atoms with sliding distance at different indentation depths (*h*<sub>c</sub>) and the same sliding velocities of 1.0 Å/ps</figcaption></figure>

elastic potential description, and the microstructure and
dynamic behavior of the material need to be
considered.

3. The connection and differences between elastic–plastic theory at micro- and macro-scales

The classical elastic–plastic theory at the macro-
scale is based on continuum mechanics, where the
stress–strain relationship is used to describe the elastic
and plastic behavior of materials. The connection and

differences between elastic–plastic theory at micro-
and macro-scales are as follows:

*   **Connection:** In MD simulation, the elastic behavior can be quantified by the second-order derivative of the potential energy function (i.e., elastic modulus), which corresponds to the concept of elastic modulus in classical theory at the macroscale. Furthermore, the plastic behavior is manifested as irreversible changes in the microstructure of the
<figure><img src="image_16.png" /><figcaption>**Fig. 13** Epoxy-quartz structure and shear strain of atoms with sliding distance at different indentation depths ($h_c$) and the same sliding velocities of 1.0 Å/ps</figcaption></figure>

material during MD simulations, such as molecular chain slip (amorphous materials) and dislocation movement (crystalline materials). These microscopic mechanisms correspond to the yield and hardening phenomena in classical theory.

* Differences: MD simulations could provide the ability to observe and understand material behavior from the atomic insights, and reveal the microscopic mechanisms that cannot be directly described by the classical theory, such as molecular chain slip and dislocation formation. On the other hand, the classical theory is applicable to the macroscopic scale and relies on empirical formulas and assumptions.

### 3.4.4 Future works

To sum up, this study is an exploratory work, and there is still more research work worth carrying out in the future.

1. To determine the most appropriate force field combination, an analysis comparing multiple force fields may elucidate which force fields are closer to the experimental data, improving the reliability of the research results. Various force field combinations could be worth investigating, such as CVFF, Tersoff, BKS, ClayFF, Morse, ReaxFF for quartz nanoindenter, and CVFF, PCFF, ReaxFF for epoxy.

2. There is a different adhesion effect between various quartz nanoindenters and epoxy substrate due to different molecular interfacial interactions [50]. Thus, the diamond and virtual nanoindenters are necessary to be investigated and compared to the quartz nanoindenter in MD nanoindentation simulation, where the diamond nanoindenter is usually used in nanoindentation experiments, and there is without adhesion effect between virtual nanoindenter and epoxy substrate. Furthermore, the adhesion effect between nanoindenter

and epoxy substrate is further investigated during the loading, holding, and unloading stage in nanoindentation process.

3. There are several nanoindenter shapes in nanoindentation experiments, such as Vickers and spherical indenters. Moreover, the shape and size of soil particles are usually diverse and complex in practical geotechnical engineering, so various shapes and sizes of nanoindenters are necessary to simulate the actual contact conditions between soil particles and FRP. Furthermore, the angle of the pyramid-shaped quartz nanoindenter is also worth studying.

4. To conduct the MD simulations that are a twin of the experiments, the coarse-graining MD simulation method for the epoxy-quartz system is necessary to establish the large MD model size and large time scale, which is worth investigating.

## 4 Conclusions

The molecular dynamics (MD) simulation method was used to investigate the interfacial mechanical behavior of epoxy and quartz, considering the effect of loading rate on their nanoindentation behavior, and the effect of sliding velocity and indentation depths on their nanoscratching behavior. The main conclusions of this study are as follows:

1. The fitting method developed in this study is used to obtain the consistent pair interaction potential parameters between epoxy and quartz, when different force fields with different LJ potentials (e.g., 12–6 and 9–6 LJ potentials) are performed for epoxy and quartz, respectively. Moreover, the fitting method could also be applied to various materials using different LJ potentials in MD simulations.
<figure><img src="image_17.png" /><figcaption>**Fig. 14** a Displacement magnitude, b surface mesh, c shear strain, and d volumetric strain of atoms within the whole system at different sliding distances under one slice along x-direction from 64 to 74 Å for the case of indentation depths ($h_c$) of 20 Å and sliding velocity of 1.0 Å/ps (the quartz indenter is hidden for better visual clarity)</figcaption></figure>

2. The overall Young's modulus of 2.075 GPa and average Poisson's ratio of 0.377 of the bulk epoxy were obtained through uniaxial traction, and these values are consistent with previous simulation and experimental works, validating the accuracy of the epoxy model.

3. Values of the reduced modulus $E_r = 2.419$ GPa and the indentation hardness $H = 1.492 \sim 2.174$ GPa were obtained. The value of $H$ was one order of magnitude larger than the one obtained by experimental nanoindentation, which might be due to the difference in nanoindenter materials, the defects of epoxy resin,

and the type of applied force field. Thus, this work suggests that various force field combinations for epoxy-quartz systems, and various nanoindenter materials for comparison are worth investigating in the future.

4. From the simulation results, abrasion models could be derived. The relationship between the indentation force $F_z$ and the indentation depth $h_c$ followed a power law for the nanoindentation process. The relationship between the sliding force $F_y$ and the sliding distance $d$ followed an exponential function for the nano-scratching process. The surface wear of the epoxy
<figure><img src="image_18.png" /><figcaption>Fig. 15 Evolution of a–e forces along three directions and f sliding force ($F_y$) with sliding distance at different sliding velocities of 0.1 ~ 2.0 Å/ps and $h_c$ of 20 Å</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16 Evolution of a shear strain and b displacement magnitude of atoms with sliding distance at different sliding velocities</figcaption></figure>

substrate was an obvious “groove” in nanoindentation and a “fan-shaped distribution” during the nano-scratching process.

5. MD nanoindentation and nanoscratching mechanisms between epoxy substrate and quartz indenter were revealed. During the nanoindentation process, the higher the loading rate, the smaller the distance of surface connecting atoms between quartz and epoxy, and the higher the indentation force $F_z$ and hardness $H$ of the epoxy substrate. Furthermore, for the

nanoscratching process, the higher the sliding velocity or the higher the indentation depths were, the higher the forces along the three directions were.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s11440-024-02503-9.

**Acknowledgements** We are very grateful to Dr. Lu-Yu WANG for his invaluable discussion and suggestions in this work. This research was financially supported by the Research Grants Council (RGC) of
<figure><img src="image_20.png" /><figcaption>**Fig. 17** Evolution of the average forces along three directions with sliding velocity at different indentation depths ($h_c$) for sliding distance between 15 and 30 Å (see Fig. 12a-e and Fig. 15a-e)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 18** MD nanoscratching mechanism diagram. ($d_{A-B}$: the distance of surface connecting atoms between quartz and epoxy; A: quartz indenter; B: epoxy substrate; v: sliding velocity; $F_y$: sliding force.)</figcaption></figure>

Hong Kong Special Administrative Region Government (HKSARG) of China (Grant No.: N_PolyU534/20, 15217220), Research Centre for Nature-based Urban Infrastructure Solutions (Project no. P0053045) of The Hong Kong Polytechnic University, and. the Project of RCRE (grantNo.:1-BBEM) of The Hong Kong Polytechnic University.

**Author's contribution** PW: Conceptualization, Methodology, Investigation, Data Curation, Writing—Review and Editing, Writing—

Original Draft, Visualization; ZYY: Conceptualization, Writing-Review and Editing, Supervision, Project administration; PYH: Writing-Review and Editing, Supervision; WX: Data Curation, Methodology.

**Funding** Open access funding provided by The Hong Kong Polytechnic University.
**Data availability** No datasets were generated or analyzed during the current study.

Declarations

**Conflict of interest** The authors declare no competing interests.

**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/.

References

1. Abyaneh MJ, El Naggar H, and Sadeghian P (2020) Numerical modeling of the lateral behavior of concrete-filled frp tube piles in sand. Int J Geomech 20(8). https://doi.org/10.1061/(ASCE)GM.1943-5622.0001725

2. Aliasghar-Mamaghani M, Khaloo A (2021) Effective flexural stiffness of beams reinforced with frp bars in reinforced concrete moment frames. J Compos Constr 25(1). https://doi.org/10.1061/(ASCE)CC.1943-5614.0001101

3. Amaral CR, Rodriguez RJS, Garcia FG, Junior LPB, Carvalho EA (2014) Impact of aliphatic amine comonomers on dgeba epoxy network properties. Polym Eng Sci 54(9):2132-2138. https://doi.org/10.1002/pen.23751

4. Anagnostopoulos CA, Kandiliotis P, Lola M, Karavatos S (2014) Improving properties of sand using epoxy resin and electrokinetics. Geotech Geol Eng 32(4):859-872. https://doi.org/10.1007/s10706-014-9763-6

5. Anagnostopoulos CA, Papaliangas TT (2012) Experimental investigation of epoxy resin and sand mixes. J Geotech Geoenviron Eng 138(7):841-849. https://doi.org/10.1061/(ASCE)GT.1943-5606.0000648

6. Chang K, Yoshioka T, Kanezashi M, Tsuru T, Tung K (2010) A molecular dynamics simulation of a homogeneous organic-inorganic hybrid silica membrane. Chem Commun 46(48):9140. https://doi.org/10.1039/c0cc02531c

7. Chen Z, Chen W, Yin J, and Malik N (2021) Shaft friction characteristics of two frp seawater sea-sand concrete piles in a rock socket with or without debris. Int J Geomech 21(7). https://doi.org/10.1061/(ASCE)GM.1943-5622.0002059

8. Dong K, Zhang J, Cao M, Wang M, Gu B, Sun B (2016) A Mesoscale study of thermal expansion behaviors of epoxy resin and carbon fiber/epoxy unidirectional composites based on periodic temperature and displacement boundary conditions. Polym Test 55:44-60. https://doi.org/10.1016/j.polymertesting.2016.08.009

9. Dutta AK, Penumadu D, Files B (2004) Nanoindentation testing for evaluating modulus and hardness of single-walled carbon nanotube-reinforced epoxy composites. J Mater Res 19(1):158-164. https://doi.org/10.1557/jmr.2004.19.1.158

10. Emami FS, Puddu V, Berry RJ, Varshney V, Patwardhan SV, Perry CC, Heinz H (2014) Force field and a surface model

database for silica to simulate interfacial properties in atomic resolution. Chem Mater 26(8):2647-2658. https://doi.org/10.1021/cm500365c

11. Frenkel D, Smit B (2001) Understanding molecular simulation: from algorithms to applications. Academic Press. https://doi.org/10.1016/B978-012267351-1/50012-2

12. Frontini P, Lotfian S, Monclús MA, Molina-Aldareguia JM (2015) High temperature nanoindentation response of Rtm6 epoxy resin at different strain rates. Exp Mech 55(5):851-862. https://doi.org/10.1007/s11340-015-9985-4

13. Frost JD, Han J (1999) Behavior of interfaces between fiber-reinforced polymers and sands. J Geotech Geoenviron 125(8):633-640. https://doi.org/10.1061/(ASCE)1090-0241(1999)125:8(633)

14. Fu K, Chang Y, Tang Y, Zheng B (2014) Effect of loading rate on the creep behaviour of epoxy resin insulators by nanoindentation. J Mater Sci Mater Electron 25(8):3552-3558. https://doi.org/10.1007/s10854-014-2055-3

15. Gao F, Stead D, Elmo D (2016) Numerical simulation of microstructure of brittle rock using a grain-breakable distinct element grain-based model. Comput Geotech 78:203-217. https://doi.org/10.1016/j.compgeo.2016.05.019

16. Goel S, Joshi SS, Abdelal G, Agrawal A (2014) Molecular dynamics simulation of nanoindentation of Fe3C and Fe4C. Mater Sci Eng A 597:331-341. https://doi.org/10.1016/j.msea.2013.12.091

17. Hagler A T, Lifson S, and Dauber P (1979) Consistent force field studies of intermolecular forces in hydrogen-bonded crystals. A benchmark for the objective comparison of alternative force fields. J Am Chem Soc 101(18):5122-5130. https://doi.org/10.1021/ja00512a002

18. Hu H, Zhang W, Xing Y, and Li X (2017) Molecular dynamics study on the effect of single grain boundary on the mechanical behavior of alpha-quartz. Int J Rock Mech Min Sci (Oxford, England : 1997) 93:94-100. https://doi.org/10.1016/j.ijrmms.2017.01.009

19. Khaloo A, Tabatabaeian M, Khaloo H (2020) The axial and lateral behavior of low strength concrete confined by Gfrp wraps: an experimental investigation. Structures 27:747-766. https://doi.org/10.1016/j.istruc.2020.06.008

20. Lam CK, Lau KT (2006) Localized elastic modulus distribution of nanoclay/epoxy composites by using nanoindentation. Compos Struct 75(1-4):553-558. https://doi.org/10.1016/j.compstruct.2006.04.045

21. Li C, Strachan A (2011) Molecular dynamics predictions of thermal and mechanical properties of thermoset polymer Epon862/detda. Polymer 52(13):2920-2928. https://doi.org/10.1016/j.polymer.2011.04.041

22. Ma X, Kang X, Cao J (2023) Origin of the elastic anisotropy of silica particles: insights from first-principles calculations and nanoindentation molecular dynamic simulations. Comput Geotech 159:105489. https://doi.org/10.1016/j.compgeo.2023.105489

23. Ma Z, Gamage RP, Zhang C (2021) Mechanical properties of A-quartz using nanoindentation tests and molecular dynamics simulations. Int J Rock Mech Min 147:104878. https://doi.org/10.1016/j.ijrmms.2021.104878

24. Ma Z, Pathegama Gamage R, Zhang C (2021) Effects of temperature and grain size on the mechanical properties of polycrystalline quartz. Comp Mater Sci 188:110138. https://doi.org/10.1016/j.commatsci.2020.110138

25. Muller A, Jang J S, Suhr J, and Gibson R F (2010) Influence of particle diameter on coefficient of thermal expansion of sio2/epoxy particulate composites. In Process SAMPE fall technology conference, Salt Lake City.

26. Oliver WC, Pharr GM (1992) An improved technique for determining hardness and elastic modulus using load and displacement
sensing indentation experiments. J Mater Res 7(6):1564-1583.
https://doi.org/10.1557/JMR.1992.1564

27. Oral I, Guzel H, Ahmetli G (2013) Determining the mechanical properties of epoxy resin (Dgeba) composites by ultrasonic velocity measurement. J Appl Polym Sci 127(3):1667-1675.
https://doi.org/10.1002/app.37534

28. Pando M, Filz G, Dove J, Hoppe E (2002) Interface shear tests on frp composite piles. international deep foundations congress, pp 1486-1500. https://doi.org/10.1061/40601(256)10

29. Plimpton S (1995) Fast parallier algorithms for short-range molecular dynamics. J COMPUT PHYS 117:1-19. https://doi.org/10.1006/jcph.1995.1039

30. Rahman R, Haque A (2013) Molecular dynamic simulation of graphene reinforced nanocomposites for evaluating elastic constants. Proc Eng 56:789-794. https://doi.org/10.1016/j.proeng.2013.03.197

31. Samad UA, Alam MA, Chafidz A, Al-Zahrani SM, Alharthi NH (2018) Enhancing mechanical properties of epoxy/polyaniline coating with addition of Zno nanoparticles: nanoindentation characterization. PROG ORG COAT 119:109-115. https://doi.org/10.1016/j.porgcoat.2018.02.018

32. Sánchez M, Rams J, Campo M, Jiménez-Suárez A, Ureña A (2011) Characterization of carbon nanofiber/epoxy nanocomposites by the nanoindentation technique. Compos B Eng 42(4):638-644. https://doi.org/10.1016/j.compositesb.2011.02.017

33. Saseendran S, Wysocki M, Varna J (2017) Cure-state dependent viscoelastic Poisson'S ratio of Ly5052 epoxy resin. Adv Manuf Polym Comp Sci 3(3):92-100. https://doi.org/10.1080/20550340.2017.1348002

34. Schneider T, Stoll E (1978) Molecular-dynamics study of a three-dimensional one-component model for distortive phase transitions. Phys Rev B 17(3):1302-1322. https://doi.org/10.1103/PhysRevB.17.1302

35. Shaia HA, Al-Asadi AK, Ramadan SH (2018) Evaluation of the interface friction between fiber-reinforced polymers and granular materials using modified shear apparatus. Int J Civil Eng Technol 9(13):1010-1016

36. Shimizu F, Ogata S, Li J (2007) Theory of shear banding in metallic glasses and molecular dynamics calculations. Mater Trans 48(11):2923-2927. https://doi.org/10.2320/matertrans. mj20076910.2320/matertrans.mj200769

37. Song X, Wang MC (2019) Molecular dynamics modeling of a partially saturated clay-water system at finite temperature. Int J Numer Anal Methods Geomech 43(13):2129-2146. https://doi.org/10.1002/nag.2944

38. Sun Y, Chen L, Cui L, Zhang Y, Du X (2018) Molecular dynamics simulation of cross-linked epoxy resin and its interaction energy with graphene under two typical force fields. Comp Mater Sci 143:240-247. https://doi.org/10.1016/j.commatsci.2017.11.007

39. Tack J (2006) Thermodynamic and mechanical properties of epon 862 with curing agent detda by molecular simulation. Texas A&M University

40. Tzetzis D, Mansour G, Tsiafis I, Pavlidou E (2013) Nanoindentation measurements of fumed silica epoxy reinforced nanocomposites. J Reinf Plast Comp 32(3):160-173. https://doi.org/10.1177/0731684412463978

41. Wan J, Zaoui A, Li X, Zheng Y (2024) Molecular dynamics simulations of the interface friction behavior between fiber-reinforced polymer pile and sand. Tribol Int 192:109288. https://doi.org/10.1016/j.triboint.2024.109288

42. Wang P, Yin Z, Zhou W, Chen W (2022) Micro-mechanical analysis of soil-structure interface behavior under constant

normal stiffness condition with dem. Acta Geotech 17(7):2711-2733. https://doi.org/10.1007/s11440-021-01374-8

43. Wei P, Yin Z, Yao C, Ren Z, Zheng Y, Ma W (2024) Atomistic origin of montmorillonite clay subjected to freeze-thaw hysteresis. J Rock Mech Geotech Eng. https://doi.org/10.1016/j.jrmge.2024.11.004

44. Wei P, Zheng Y, Zaoui A, Yin Z (2023) Atomistic study on thermo-mechanical behavior and structural anisotropy of montmorillonite under triaxial tension and compression. Appl Clay Sci 233:106817. https://doi.org/10.1016/j.clay.2023.106817

45. Xie J, Chen K, Yan M, Guo J, Xie Q, Lü F (2022) Effect of temperature and water penetration on the interfacial bond between epoxy resin and glass fiber: a molecular dynamics study. J Mol Liq 350:118424. https://doi.org/10.1016/j.molliq.2021.118424

46. Xu WQ, Yin ZY, Zheng YY (2023) Frp-soil interfacial mechanical properties with molecular dynamics simulations: insights into friction and creep behavior. Num Anal Meth Geo-mech. https://doi.org/10.1002/nag.3607

47. Xu W, Yin Z, Zheng Y (2023) Investigating silica interface rate-dependent friction behavior under dry and lubricated conditions with molecular dynamics. Acta Geotech. https://doi.org/10.1007/s11440-022-01792-2

48. Yang H, He M, Lu C, Gong W (2019) Deformation and failure processes of kaolinite under tension: insights from molecular dynamics simulations. Sci China Phys Mech Astron 62(6):64612. https://doi.org/10.1007/s11433-018-9316-3

49. Yang W, Ayoub G, Salehinia I, Mansoor B, Zbib H (2018) The effect of layer thickness ratio on the plastic deformation mechanisms of nanoindented Ti/Tin nanolayered composite. Comp Mater Sci 154:488-498. https://doi.org/10.1016/j.commatsci.2018.08.021

50. Yaphary YL, Yu Z, Lam RHW, Hui D, Lau D (2017) Molecular dynamics simulations on adhesion of epoxy-silica interface in salt environment. Compos B Eng 131:165-172. https://doi.org/10.1016/j.compositesb.2017.07.038

51. Zappalorto M, Pontefisso A, Fabrizi A, Quaresimin M (2015) Mechanical behaviour of epoxy/silica nanocomposites: experiments and modelling. Compos A Appl Sci Manuf 72:58-64. https://doi.org/10.1016/j.compositesa.2015.01.027

52. Zhang L, Zheng Y, Wei P, Diao Q, Yin Z (2021) Nanoscale mechanical behavior of kaolinite under uniaxial strain conditions. Appl Clay Sci 201:105961. https://doi.org/10.1016/j.clay.2020.105961

53. Zhang X, Xiong W (2009) Effect of bonding process on the properties of isotropic epoxy resin-bonded Nd-Fe-B magnets. Rare Met 28(3):248-252. https://doi.org/10.1007/s12598-009-0049-8

54. Zhang Z, Song X (2022) Nonequilibrium molecular dynamics (nemd) modeling of nanoscale hydrodynamics of clay-water system at elevated temperature. Int J Numer Anal Met 46(5):889-909. https://doi.org/10.1002/nag.3325

55. Zhang Z, Song X (2023) Nanoscale crack propagation in clay with water adsorption through reactive Md modeling. Int J Numer Anal Met 47(7):1103-1133. https://doi.org/10.1002/nag.3507

56. Zhu W, Hughes JJ, Bicanic N, Pearce CJ (2007) Nanoindentation mapping of mechanical properties of cement paste and natural rocks. Mater Charact 58(11-12):1189-1198. https://doi.org/10.1016/j.matchar.2007.05.018

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