The Influence of the Grain Size Effect on the Mechanical Properties of Metallic Tungsten during Nanoindentation

1  Introduction

Semiconductor technology is extensively applied across critical fields such as artificial intelligence and electronic communication. The continuous downscaling of transistor dimensions has been driven by Moore’s Law, proposed in 1965 [1], and further enabled by Dennard scaling, introduced in 1974 [2]. Chau et al. successfully demonstrated CMOS transistors with physical gate lengths of 80 nm using high-k/metal-gate electrodes employing n+ and p+ work functions [3]. To meet the requirements of nanoscale device fabrication, advanced computational lithography techniques such as Inverse Lithography Technology (ILT) have been widely adopted in semiconductor manufacturing [4,5]. Reliable electrical interconnection is essential for the integration of semiconductor devices [6]. Tungsten (W), due to its relatively low bulk resistivity and high diffusion activation energy, has emerged as a promising alternative to copper for interconnect applications [7]. In Very-Large-Scale Integration (VLSI) technology, chemical vapor deposition of tungsten (W-CVD) is considered an attractive solution for filling and interconnects [8]. Atomic layer deposition (ALD) of tungsten thin films have broad applications in semiconductor device fabrication, including tungsten plug processes and tungsten metal lines such as gate electrodes, bit lines in Dynamic Random-Access Memory (DRAM), and gate electrodes in 3D NAND flash memory, owing to its low resistivity, high melting point, high mechanical hardness, and chemical inertness [9]. Recently, in-situ W-doped ALD has been employed to fabricate In2O3 semiconductor layers, achieving a 24.4-fold performance enhancement in Compute-In-Memory (CIM) hardware compared to state-of-the-art resistive synapses [10]. These advancements underscore the critical role of tungsten in semiconductor interconnect technology. Therefore, investigating nanoscale processing techniques for tungsten is essential for advancing next-generation interconnect solutions in nanoscale semiconductor devices.

Nanoscale simulations provide a valuable tool for characterizing the mechanical properties of nanostructured materials, particularly in understanding plastic deformation mechanisms and identifying dislocation sources [11–13]. Molecular dynamics (MD) simulations, in particular, offer critical insights into the microscopic deformation behavior of tungsten during nanoscale processing [14]. At present, molecular dynamics simulations of metallic tungsten have covered a wide range of topics, with a particular focus on its behavior during cutting processes. Guo et al. investigated crystal structure evolution and dislocation dynamics during cutting, revealing that dislocation nucleation within the workpiece significantly influences subsurface defect formation [14]. Wang et al. developed an MD simulation framework to study the surface formation mechanism in ultrahigh-speed cutting of single-crystal tungsten, incorporating ultrasonic elliptical vibration cutting [15]. Their work provided atomic-level insights into the surface/subsurface formation mechanisms during ultrasonic-assisted eutectic processes and examined the influence of ultrasonic amplitude and frequency on deformation behavior, as well as the stress-induced plasticity and microstructural evolution in tungsten [15,16]. Dong et al. proposed a novel method for determining the minimum chip thickness based on atomic trajectory recognition and relative coordinate analysis, and identified key factors influencing chip formation during tungsten removal [16]. Additionally, Wang et al. conducted a comprehensive study on the effect of tool geometry on the removal behavior of metallic tungsten, concluding that larger rake and clearance angles, as well as smaller tool edge radii, are beneficial for improving post-removal surface integrity [17]. Beyond cutting, several studies have also focused on the indentation behavior of tungsten. For example, Valencia et al. employed MD simulations to investigate the indentation response of nanoporous tungsten, highlighting the roles of twin nucleation and dislocation motion in influencing its plasticity and mechanical performance, thereby offering insights into future applications of this novel material [18]. Moreover, Lin et al. combined crystal defect theory, MD simulations, and machine learning to study the impact of helium bubbles on the indentation behavior of pure tungsten, providing a deeper understanding of helium-induced damage mechanisms [19].

Grain refinement is widely recognized as an effective strategy for enhancing the mechanical properties of materials [20], and it can provide crucial support for the nanoscale processing of tungsten. Shi and Meng found that grain refinement leads to an increase in dislocation density during scratching, where most dislocations are short and confined by grain boundaries [21]. In nanocrystalline polycrystalline substrates, Hall-Petch breakdown was observed in their work, with the dominant deformation mechanism shifting from dislocation motion to lattice rotation and grain boundary sliding. Cao et al. investigated the effects of grain size on atomic-scale plastic flow mechanisms, defect evolution and distribution, average cutting forces, and stress variations during the nanocutting of polycrystalline γ-TiAl alloys [22]. Their work addressed the lack of research on atomic plasticity in polycrystalline γ-TiAl during nanoscale machining. Most existing studies have focused on the effects of grain size on the mechanical properties of bulk and nanocrystalline materials. However, the influence of workpiece grain size on cutting forces remains controversial and inconsistent in the literature [23]. During the cutting process, especially as the machining scale decreases, the grain size, crystal orientation, and mechanical properties of polycrystalline materials play increasingly important roles in determining the overall mechanical response. Moreover, variations in dislocation density and strain-rate hardening behavior may also contribute to the grain-size dependence of the Rehbinder effect. It is worth noting that previous studies on the grain size effect have predominantly focused on materials with face-centered cubic (FCC) crystal structures, while significantly fewer investigations have been conducted on body-centered cubic (BCC) systems. Furthermore, existing research on tungsten has largely concentrated on nanoscale cutting processes. As a result, there is a lack of comprehensive evidence supporting the influence of grain size on the mechanical behavior of tungsten under nanoindentation, leaving a critical gap in understanding grain size effects in BCC-structured materials such as tungsten. Unlike FCC metals, BCC tungsten exhibits temperature-dependent screw dislocation mobility and asymmetric slip systems [24]. Recent nanoindentation studies reveal unique 12<111> dislocation nucleation in tungsten [19], but grain size effects remain unexplored.

In this work, the grain size effect on tungsten, a material with a typical body-centered cubic (BCC) crystal structure, is systematically investigated. Six groups of tungsten crystals with different grain sizes are examined to evaluate their behavior in terms of surface pile-up, atomic displacement, elastic recovery, loading force, hardness, stress and strain response, as well as dislocation nucleation and evolution. The contribution of this work is a systematic, grain-size-resolved mapping of nanoindentation responses in tungsten, indicating a phenomenological transition around ~7.3 nm within the tested range and its mechanistic signatures across multiple metrics. This study contributes to a deeper understanding of the plastic deformation mechanism and dislocation behavior of tungsten, providing important guidance for the practical application and material processing of metallic tungsten.

2  Computation Settings

In this simulation, the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [25] program is used to achieve an MD simulation of nanoindentation of crystal W. The overall schematic diagram of the constructed MD model is shown in Fig. 1. The crystal orientation of the workpiece is randomly generated, with an initial size of 24 × 24 × 24 nm3, containing 871,143 atoms. Because the grains are randomly oriented, no single crystallographic direction is aligned with the indentation axis; the indenter moves along the global -Z direction. The spherical diamond indenter was treated as a rigid body, and elastic deformation of the indenter was not allowed. To ensure precise depth measurement, the spherical indenter with a radius of 4 nm was initially positioned such that its lowest point was 1 nm above the highest atom of the tungsten surface. This predefined gap established a consistent reference plane (Z = 0) for all simulation. The tool then moves at a constant speed of 1 nm/s in the -Z direction, indicated by the arrow in Fig. 1. The indentation depth h is defined as the vertical displacement of the indenter center from this initial position, ensuring all depth values are measured from the same reference plane. The indentation depth h is set with a gradient change within the range of 1 to 5.3 nm to explore the grain size effect. Although the maximum depth (5.3 nm) exceeds the indenter radius, it remains within the nanoscale and is selected to ensure clear plastic deformation and dislocation activity in MD simulations; therefore, it is still considered nanoindentation in the present context. The workpiece is divided into a Boundary layer, a Thermostat layer, and a Newtonian layer. The 1 nm area at the bottom of the workpiece is set as a fixed layer to suppress the overall displacement, and the 1 nm area above it is set as a constant temperature layer to maintain the stability of the system. The 1 nm thermostat layer was controlled by a Langevin thermostat with a damping constant of 1.0 to dissipate the heat generated by plastic deformation. The polycrystalline tungsten models with varying grain sizes (4.8–7.9 nm) were generated via Voronoi tessellation algorithm [26]. A single representative Voronoi realization was generated for each grain size, as the present study focuses on systematic grain-size trends rather than statistical variability. Each grain was assigned a random crystallographic orientation to eliminate preferential texture effects, following the methodology established for BCC metals [27]. The grain boundary structure was optimized through energy minimization with conjugate gradient method, ensuring equilibrium at 293 K prior to indentation.

Figure 1: MD model of W nanoindentation. The arrow indicates the indentation direction (-Z). The bottom fixed boundary layer is shown in deep blue, the thermostat layer in light blue, and the Newtonian layer in orange.

In the simulation process, the system is relaxed via energy minimization using the conjugate gradient method. Subsequently, it is balanced for 30 ps at a temperature of 293 K under the canonical ensemble (NVT) to ensure that the system reaches the thermodynamic equilibrium state. The nanoindentation process adopts the microcanonical ensemble (NVE). The model employs mixed boundary conditions, with the Y direction specified as a periodic boundary condition and the X and Z directions designated as free boundary conditions. The accuracy of MD simulation determines its reliability [28]. To ensure the accuracy of MD simulation, an appropriate potential function needs to be selected. The interactions between different atomic pairs were assigned as follows:

C-C atomic interactions (within diamond indenter) are modeled by the Tersoff potential [29], suitable for covalent materials like diamond. The Tersoff potential function takes into account covalent bonds and bond angles, which can be used to simulate the interactions between atoms, as shown in Eq. (1).

{E=12∑i∑i≠jVijVij=fc(rij)[fR(rij)+bijfA(rij)](1)

where rij represents the distance between particle i and particle j, Vij represents the bond energy between atom i and atom j, fR is the repulsive pair potential, fA is the attractive pair potential, and fc is the smooth cutoff function.

W-W atomic interactions are described by the EAM/FS potential [30], which accurately reproduces the vacancy formation energy and mechanical properties of BCC tungsten. The EAM/FS potential is a modified and more general form of the traditional Embedded Atom Method (EAM) potential. Compared with the original EAM potential, it offers improved accuracy in reproducing experimental vacancy formation energies and does not require externally applied pressure to balance the Cauchy pressure. As shown in Eq. (2).

Ei=Fα(∑j≠iραβ(rij))+12∑j≠i∅αβ(rij)(2)

where Fα is the embedding energy, a function of the atomic electron density. ραβ and ∅αβ are a pair potential interaction, α and β are the element types of atom i and atom j, respectively.

W-C cross-interactions are governed by the Morse potential [31], effectively capturing the pairwise interactions between dissimilar elements. The specific Morse parameters for W-C were used without additional fitting for the indenter-workpiece interface. The Morse potential energy function (PEF) is highly suitable for describing the effective pair interaction forces in metals, as shown in Eq. (3)

E=D[e−2α(rij−r0)−2e−α(rij−r0)](3)

where the constant D=0.9001 eV represents the dissociation energy, α=26.364nm−1 is a constant with dimensions of reciprocal length, and r0=0.16277 nm is the equilibrium distance between particle i and particle j [32,33].

In conclusion, the interatomic interactions are described by Tersoff, EAM/FS, and Morse potential functions, while potential-model sensitivity is not evaluated here. The detailed parameters during the simulation process can be seen in Table 1. Furthermore, the results calculated using the LAMMPS program cannot be directly observed; therefore, the OVITO software [34] is selected for visualizing the simulation results. However, a dedicated timestep convergence test and a computational cost analysis were not included in this study.

3  Results and Discussion

3.1 Surface Morphology and Atomic Displacement

Fig. 2 shows the machined surface and pile-up morphology of the W matrix with a grain size of 7.9 nm, where h represents indentation depth, and t represents time. For ease of observation, the tool model is removed in Fig. 2a–f, and atoms in the matrix are colored according to the height in the Z direction. It can be seen from the figure that when the indentation depth reaches h = 3.5 nm, the surface pile-up phenomenon begins to occur, and with the increase of the indentation depth, the number of atoms within the pile-up regions gradually increases. When the loading process is completed and the indentation depth reaches h = 5.3 nm, the number of pile-up atoms reaches 699. Unloading begins when the indentation depth reaches h = 5.3 nm. Contrary to expectations, the number of pile-up atoms increases rather than decreases during the initial stage of unloading. After unloading for 12 ps (at this time, the total indentation time is 72 ps), the number of pile-up atoms reached the maximum value of 2675. Subsequently, as the unloading process continued, the number of pile-up atoms showed a non-monotonic variation, but overall, it still showed a downward trend, as shown in Fig. 2i. It can be inferred that the fluctuation of the number of pile-up atoms during the unloading process is due to the combined effect of elastic recovery and the interaction force between the tool and the workpiece. When the unloading is completed, that is, t = 120 ps, the number of pile-up atoms is 1939, which is lower than the maximum value.

Figure 2: At a grain size of 7.9 nm, the machined surface and pile-up morphology of the W matrix. (a–c) Surface pile-up morphology during loading (0–60 ps); (d–f) Surface pile-up morphology during unloading (60–120 ps); (g) Schematic diagram of loading/unloading; (h) Unloading cross-section view; (i) Number of pile-up atoms during the process.

The simultaneous pile-up peak at t = 72 ps is attributed to delayed elastic recovery during unloading. When the indenter retracts, grain boundary atoms experience reduced constraints, enabling dislocation rearrangement that temporarily increases surface pile-up. This transient phenomenon aligns with the stress relaxation mechanism in ultrasonic-assisted processing [15], where delayed recovery occurs under dynamic unloading conditions.

At different grain sizes, the surface pile-up morphologies exhibit slight variations, as shown in Fig. 3. In addition, the area of the surface pile-up also differs. Specifically, except for the case of a grain size of 7.9 nm, the surface pile-up area tends to increase overall as the grain size decreases.

Figure 3: The surface pile-up morphology of the W matrix at t = 72 ps (when the number of pile-up atoms reaches its maximum) for different grain sizes. (a) Surface pile-up morphology of grain size at 7.9 nm; (b) 7.3 nm; (c) 6.9 nm; (d) 6.6 nm; (e) 6.4 nm; (f) 4.8 nm.

The pile-up size is defined as the height of the atomic surface pile-up, as illustrated in Fig. 4a. As shown in Fig. 4b, the pile-up sizes remain below 1.2 nm. In general, the pile-up size increases with increasing grain size. The grain size effect has a significant influence on the pile-up size. In addition, Fig. 4c shows the variation trend of the number of pile-up atoms over time at different grain sizes during the nanoindentation process. The trends of the number of pile-up atoms at other grain sizes exhibit similarity to those observed at 7.9 nm. Among these six groups of data, at a grain size of 6.4 nm, the number of pile-up atoms is the largest, and the time when pile-up occurred is also the earliest (at t = 36 ps). At a grain size of 7.3 nm, the number of pile-up atoms is the smallest, and the time of pile-up atom occurrence is the latest (at t = 48 ps). As shown in Fig. 4c, the pile-up area fraction is represented by the number of surface pile-up atoms, and its quantitative variation with grain size. It can be observed from Fig. 4c that for the four groups of data with grain sizes of 7.3, 6.9, 6.6, and 6.4 nm, as the grain size decreases, the number of pile-up atoms gradually increases, and the time of pile-up atoms also gradually advances. These four groups of data have a relatively obvious grain size effect.

Figure 4: (a) Schematic diagram of the pile-up size; (b) Curve of the pile-up size at different grain sizes; (c) Curve of the number of surface pile-up atoms. For clarity, the data from the first 30 ps (when the number of atoms is zero) is omitted.

In conclusion, within the tested range, the trends differ for grain sizes smaller than ~7.3 nm vs. larger than ~7.3 nm. For grain sizes larger than ~7.3 nm, further increases lead to opposite trends in the number of pile-up atoms and the pile-up area.

It can be seen from Fig. 4c that the number of pile-up atoms decreases after unloading is completed, which is related to elastic recovery. This can be observed in Fig. 5a–f. After unloading is completed (t = 120 ps), the width and depth of the indentation are reduced compared with those after loading is completed (t = 60 ps). As shown in Fig. 5g, define the widths of the workpiece after loading and unloading as w0 and w1, respectively, and the depths of the workpiece after loading and unloading as h0 and h1, respectively. The elastic recovery ratio is defined as the ratio of the change in the width or depth direction before and after unloading. Specifically, the elastic recovery ratio in the width direction is shown as Eq. (4), and the elastic recovery ratio in the depth direction is shown as Eq. (5).

wer=(w0−w1)/w0(4)

her=(h0−h1)/h0(5)

Figure 5: Elastic recovery of W before and after unloading at different grain sizes. (a–f) Elastic recovery morphology diagrams, where the upper one is the image after loading is completed (t = 60 ps), and the lower one is the image after unloading is completed (t = 120 ps); (g) Characteristic dimensions of depth and width before and after unloading; (h) Elastic recovery ratio curve.

These two elastic recovery ratios at different grain sizes are shown in Fig. 5h. The results show that when the grain size is lower than the phenomenological transition size observed in this study, which is 7.3 nm, with the increase in the grain size, the elastic recovery ratio in the depth direction shows a significant increasing trend in general, whereas the elastic recovery ratio in the width direction fluctuates within the range of 0%–2%. The relationship between grain size and elastic recovery ratio is related to the elastic, inelastic, and plastic deformations during the nanoindentation process [35]. Specifically, the increase of the elastic recovery ratio in the depth direction is related to the number of grain boundaries. The atoms at the grain boundaries are disordered in arrangement and have a lower density than those in the grains. Compared with the atoms in the grains, the atoms at the grain boundaries have higher potential energy and are more likely to move under load; they are more likely to undergo plastic deformation. Therefore, during the nanoindentation process, the atoms at the grain boundaries store less elastic energy than those in the grains. The above conclusions are consistent with the findings reported by Li et al. [36]. When the grain size is below the phenomenological transition size, the increase in grain size leads to a decrease in the number of grain boundaries; meanwhile, the weak interfacial connection of the grain boundaries reduces the interatomic repulsive and attractive forces [35], resulting in a stronger elastic recovery ratio. Furthermore, since the nanoindentation process mainly occurs in the depth direction, the elasticity in this direction is greater [36], and the elastic recovery ratio is more obvious after unloading.

For the W crystal with different grain sizes, the overall atomic displacement after loading remains largely similar, as shown in Fig. 6a1–f1. Similar to the case at a grain size of 7.9 nm, surface atoms exhibit atomic displacement under various grain sizes, whereas atoms located directly beneath the spherical indenter exhibit relatively smaller atomic displacements.

Figure 6: Atomic displacement vector diagrams of the central cross-section of the W matrix at different grain sizes. (a1–f1) Top views of the atomic displacement of the W matrix; (a2–f2) Views in the Y direction at an indentation depth of 5.3 nm.

To better observe atomic displacements beneath the indenter, a central cross-sectional view is generated, as illustrated in Fig. 6a2–f2. In most cases, the atomic displacement beneath the indenter exhibits an asymmetric conical morphology. However, at grain sizes of 7.3 and 4.8 nm, the atomic displacement pattern is nearly spherical. By comparing the crystal structures at different grain sizes, it is observed that, unlike other grain sizes, no grain boundary intersection occurs in the region beneath the indenter in these two specific cases. Therefore, it is inferred that when grain boundary intersections exist beneath the indenter, atoms tend to move along the grain boundaries, resulting in atomic displacements that predominantly occur at the grain boundaries [37]. In summary, the differences in atomic displacement primarily depend on the distribution of grain boundaries within the crystal, whereas their correlation with grain size remains relatively weak.

3.2 Force and Hardness

In nanocrystalline materials, pile-up can be highly localized, which may alter the true contact area. To ensure consistent comparison across grain sizes, we use the same contact-area definition for all cases; therefore, the reported hardness values should be interpreted as comparative trends rather than absolute values. At different grain sizes, the force of W shows an approximate changing trend. With the increase of indentation depth, the normal force Fz increases significantly. Approximately at the end of the loading process (t = 60 ps, h = 5.3 nm), the normal force Fz reaches its maximum value. Then, during the progress of the unloading process, the normal force Fz drops sharply, as shown in Fig. 7a–f. However, the forces Fx and Fy in the X and Y directions remain nearly constant with increasing indentation depth, exhibiting only minor fluctuations around zero. For the normal force Fz, the force–indentation depth curve exhibits a growth trend that approximately follows a power-law function, and it can be fitted through the Hertz elastic model, whose formula is shown in Eq. (6).

F=34Er×R12×h32(6)

where the R represents the indenter radius, and h represents the indentation depth.

Figure 7: The relationship between triaxial forces and indentation depth under different grain sizes during the crystal W processing. (a–f) The variation curves of Fx, Fy, and normal force Fz with indentation depth; (g) The variation curves of the normal force Fz of different grains during the nanoindentation process; (h) The average value of the normal force Fz during the loading process under different grain sizes.

The reduced elastic modulus is shown as Eq. (7).

Er=Es/(1−vs2)(7)

where Es and vs are Young’s modulus and Poisson’s ratio, respectively.

During the loading process, the force-indentation depth curve reveals that as the indentation depth increases, the normal force Fz exhibits non-monotonic behavior with occasional drops, due to plastic rearrangement and strain hardening [38]. Specifically, at indentation depths of approximately 1.5, 2.5, and 4.25 nm, noticeable drop points appear, which are attributed to misalignment nucleation occurring beneath the indenter. It can be reasonably inferred that dislocation nucleation begins at an indentation depth of approximately 1.5 nm and then propagates within the crystal as the indenter further penetrates the matrix. It is precisely because of the dislocation nucleation on the indentation surface that the force steadily increases with occasional large or small decreases [19]. This phenomenon can also be seen from the hardness-indentation depth curve, as shown in Fig. 7.

However, at different grain sizes, the variation trend of the normal force Fz is slightly different. Although in most cases, the force decreases near indentation depths h = 1.5 nm, h = 2.5 nm, and h = 4.25 nm, the extent of the decrease is different at different grain sizes. In addition, the decrease is less pronounced at grain sizes of 6.6 and 4.8 nm than at other sizes. However, for the workpiece with a grain size of 6.9 nm, a relatively pronounced drop occurs at an indentation depth of approximately 5 nm, as shown in Fig. 7c. From this phenomenon, it can be speculated that the grain size affects the dislocation nucleation behavior.

In addition, at the grain size of 7.3 nm, the force reaches the maximum value of approximately 2651 nN when the indentation depth h = 4.91 nm. Although occasional fluctuations occur, as the grain size decreases, the maximum value of the corresponding force shows an overall downward trend and reaches the minimum value of 2235 nN at the grain size of 4.8 nm. At this time, the indentation depth is h = 4.84 nm. It can be observed from Fig. 7g that the force–indentation depth curve for the grain size of 7.3 nm generally lies at the top, whereas the curve for the grain size of 4.8 nm is predominantly located at the bottom. Approximately at h = 5.3 nm, the difference in force between the 7.3 nm grain size and the 4.8 nm grain size reaches a maximum of approximately 393 nN. These results suggest that the force tends to decrease with decreasing grain size, although the effect of grain size appears to be relatively minor. This can be seen more clearly from the bar chart of average force in Fig. 7h. At the grain sizes of 7.3, 6.9, 6.6, 6.4, and 4.8 nm, as the grain size decreases, the average normal force also decreases, although the reduction is minimal. The magnitude of the indentation force is closely related to the mode of plastic deformation [39]. Under the action of the external load, the sliding and rotation of the grain boundaries release stress successively, when the grain size is smaller and the number of grain boundaries increases, more stress is released, facilitating dislocation motion and resulting in increased plastic deformation [36,38,39], which leads to the lowest force at 4.8 nm. For the unloading process, at the grain sizes of 7.3, 6.9, 6.6, 6.4, and 4.8 nm, the grain size effect is still followed, with the reduction of the grain size, the average normal force decreases.

Hardness is defined as the resistance of a material to localized plastic deformation [36], which significantly influences its wear performance [40]. It can be calculated by Eq. (8).

H=F/A(8)

where F is the load applied during the nanoindentation process, and A is the contact area of the load applied between the indenter and the crystal W.

For spherical indenters, the contact area can be defined by Eq. (9).

A=π×[R2−(R−h)2](9)

where R is the indenter radius and h is the indentation depth.

It can be observed from Fig. 8a–f that the hardness generally exhibits an upward trend during the loading process, increasing approximately according to a power-law relationship. This may be attributed to the existence of the fixed dislocation network that hinders the mobility of dislocations and thereby leads to strain hardening [38]. As mentioned above, due to plastic rearrangement and strain hardening, the hardness exhibits slight fluctuations with increasing indentation depth, which can be confirmed from Fig. 8g. Moreover, due to dislocation nucleation, three decline points also appeared on most of the hardness-indentation depth curves, and the timing of their occurrence corresponds to the decline points on the force-indentation depth curves. Similarly, at the grain size of 6.8 nm, there is an additional decline point. The hardness-indentation depth curve illustrates a steady increase, possibly because in body-centered cubic nanocrystals, a continuously increasing stress is required to initiate dislocation nucleation, whereas grain boundary plastic deformation tends to be more pronounced [18].

Figure 8: The relationship between hardness and indentation depth under different grain sizes during the W processing. (a–f) The variation curve of hardness with indentation depth; (g) The variation curves of hardness of different grains during the indentation process; (h) The average hardness under different grain sizes.

Similar to the conclusion obtained from the force-indentation depth curve, the value of hardness decreases with the reduction of grain size. However, the grain size effect on hardness is relatively minor, which can be confirmed from Fig. 8h. In Fig. 8h, at the grain sizes of 7.3, 6.9, 6.6, 6.4, and 4.8 nm, the average hardness decreases with the reduction of the grain size, but the decrease is very small. It can be observed in Fig. 8g that approximately at h = 5.3 nm, the hardness difference of the 7.3 nm grain size reaches the maximum, approximately 11 GPa.

Consistent with the conclusions obtained from atomic pile-up and elastic recovery, it can also be inferred from the variation trends of force and hardness that after exceeding the phenomenological transition size, the variation trends of force and hardness undergo a reversal [36], that is, they will decrease as the grain size increases. Since grain orientations are randomly assigned, the observed trends represent qualitative correlations rather than precise quantitative relationships. Notably, the grain-size dependence of the normal force and hardness (Figs. 7 and 8) is consistent with the dislocation-density evolution shown in Figs. 9 and 10: grain sizes that exhibit lower dislocation density and a lower 1/2<111> proportion also show lower hardness/force at the same indentation depth.

Figure 9: The dislocation density curves of different indentation depths during the processing of the W matrix at different grain sizes. (a) The dislocation density curve of grain size at 7.9 nm; (b) 7.3 nm; (c) 6.9nm; (d) 6.6 nm; (e) 6.4 nm; (f) 4.8 nm.

Figure 10: Bar chart of the proportion of dislocation density at different indentation depths during the processing of W matrix at different grain sizes. (a) The proportion of dislocation density of grain size at 7.9 nm; (b) 7.3 nm; (c) 6.9 nm; (d) 6.6 nm; (e) 6.4 nm; (f) 4.8 nm.

3.3 Stress and Shear Strain

In this section, we first show the depth-dependent evolution of stress and shear strain for a representative large grain size (7.9 nm), and then compare different grain sizes at the maximum indentation depth (5.3 nm) to highlight grain-size effects.

Fig. 11 demonstrates the local von Mises stress during the loading and unloading process. It can be observed that the local von Mises stress increases with increasing indentation depth and propagates downward, exhibiting an asymmetric distribution. Fig. 11e,f further illustrates that, during the unloading process, a small amount of residual stress remains and does not fully return to the initial state, which is typically concentrated near the grain boundaries. In addition, low levels of local von Mises stress are observed at the grain boundaries.

Figure 11: The distribution of the local von Mises stress within the central cross-section in the Y direction of the matrix for the crystal W with a grain size of 7.9 nm at different indentation depths. (a–d) The distribution of the local von Mises stress during the loading process; (e,f) The distribution of the local von Mises stress during the unloading process.

Fig. 12 shows that the high-stress region at 5.3 nm depth is confined near the indenter and does not extend to the fixed bottom layer, indicating limited substrate influence. When the grain size is smaller than the phenomenological transition size of 7.3 nm, with the decrease of the grain size, the local von Mises stress also decreases accordingly. Conversely, as the grain size decreases, the area of the region affected by higher stress increases accordingly, as indicated by the green area in Fig. 12. At the maximum indentation depth, the high-stress region remains confined near the indenter and does not reach the fixed bottom layer, suggesting that substrate constraint does not dominate the measured forces. Although the indenter diameter is 8 nm, the 24 nm × 24 nm × 24 nm workpiece provides sufficient clearance; the high-stress region at the maximum indentation depth remains localized beneath the indenter, suggesting limited boundary influence on the trends reported here. In addition, the green region in Fig. 12 indicates that the locations of higher stress are consistent with the distribution of grain boundaries. The occurrence of the above phenomena can be attributed to the higher energy of atoms located at grain boundaries compared to those within grain interiors. This energy difference results in elevated stress levels at the grain boundaries, which in turn facilitate stress transmission along these interfaces. When the grain size increases, the number of atoms situated at the grain boundaries decreases, reducing the obstacles to stress transfer and allowing stress to propagate more easily through the boundaries and into the grain interiors. In contrast, when the grain size is small, the relative fraction of grain boundaries increases, leading to a larger total interfacial area between grains. This amplifies stress concentration at the grain boundaries. The above discussions of the results are consistent with the conclusions reported by Cao et al. [39].

Figure 12: The distribution of the average local von Mises stress within the central cross-section in the Y direction of the matrix at the indentation depth of 5.3 nm for different grain sizes. (a) The average local von Mises stress of grain size at 7.9 nm; (b) 7.3 nm; (c) 6.9 nm; (d) 6.6 nm; (e) 6.4 nm; (f) 4.8 nm.

A phenomenological transition around ~7.3 nm is observed because several independent metrics show a consistent turning point at this grain size. As the grain size increases from 6.4–7.3 nm, the number of pile-up atoms and pile-up area decrease while elastic recovery and hardness increase; beyond 7.3 nm, these trends weaken or reverse. The stress/strain localization near the indenter also shows a corresponding change. The convergence of these indicators supports 7.3 nm as a transition point in deformation behavior.

As shown in Fig. 13a–d, atoms located beneath the indenter are subjected to higher shear strain, which intensifies with increasing indentation depth and extends both downward and laterally around the indenter. For the convenience of observation, atoms in the matrix are colored according to the shear strain magnitude. In other words, the region of high strain is centered around the indenter, exhibits near-uniform distributions [39]. Specifically, the higher strain primarily occurs in the Type 6 region and gradually spreads toward the boundary between Type 6 and Type 7 as the indentation depth increases. Fig. 13e,f depicts the shear strain distribution after unloading. It can be observed that the shear strain does not exhibit a significant reduction after unloading, which provides evidence of the presence of plastic deformation.

Figure 13: The distribution of shear strain within the central cross-section in the Y direction of the matrix for crystal W with a grain size of 7.9 nm. (a–d) The distribution during the loading process; (e,f) The distribution during the unloading process; (g,h) Two grain types in the cross-sectional view at an indentation depth of 5.3 nm.

It can be seen from Fig. 14 that the degree of shear strain presents an asymmetric characteristic due to the existence of different types of grains in the crystal. For the crystal W at grain sizes of 6.9, 6.6, 6.4, and 4.8 nm, as the grain size decreases, the shear strain shows an overall decreasing trend. However, a reduction in grain size leads to a relative expansion of the shear-strain-affected region, as highlighted by the green areas in Fig. 14. The results demonstrate that when grain boundaries are densely distributed, changes are more readily propagated into the matrix via the boundaries [36].

Figure 14: The distribution of shear strain within the central cross-section in the Y direction of the matrix at the indentation depth of 5.3 nm for different grain sizes. (a) The distribution of shear strain of grain size at 7.9 nm; (b) 7.3 nm; (c) 6.9 nm; (d) 6.6 nm; (e) 6.4 nm; (f) 4.8 nm.

3.4 Dislocation and Dislocation Density

Due to the extrusion of the indenter, the atomic energy beneath the machined surface increases and accumulates. Upon reaching a sufficient energy level, energy is released and dislocation nucleation is initiated [17], and the deformation takes place under the condition of dislocation accumulation [38].

Dislocation density is an important indicator of the dislocation response during nanoindentation, and the residual dislocation density has a significant effect on the mechanical properties of the material [39]. Therefore, it is an important parameter for evaluating the mechanical properties and plastic deformation of materials. Although the grains are randomly oriented, a qualitative correlation is observed: higher indentation force and hardness generally coincide with denser dislocation networks and more pronounced closed-loop dislocations, whereas lower mechanical responses correspond to more scattered dislocation structures. Quantitative dislocation metrics are provided in Figs. 9 and 10. These plots quantify nucleation/evolution trends across grain sizes. Fig. 9 provides the evolution of total dislocation density with indentation depth, which qualitatively tracks the grain-size trend of hardness/force in Figs. 7 and 8. Different types of dislocation densities exhibit similar trends with increasing indentation depth at the same grain size; however, these trends differ markedly across different grain sizes. Specifically, at a grain size of 7.9 nm, when the indentation depth reaches approximately 0.5 nm, the dislocation density begins to increase sharply, which is consistent with the beginning of dislocation nucleation and corresponds to the significant plastic deformation caused by the spherical indenter [19]. The dislocation density starts to decline at an indentation depth of around 1.5 nm and drops sharply near 2.25 nm. Subsequently, it tends to stabilize near the indentation depth h = 2.5 nm. At the grain size of 6.9 nm, the dislocation density begins to decrease sharply when the indentation depth reaches approximately 0.25 nm, then increases significantly near h = 1 nm. Subsequently, the dislocation densities of different types exhibit fluctuations with closely aligned rates as the indentation depth increases. The irregular fluctuations described above are attributed to dislocation nucleation-induced plastic deformation occurring in the matrix under high stress; meanwhile, elastic recovery takes place, releasing stored elastic energy, which leads to either the regeneration or annihilation of dislocations [39]. As a result, the dislocation density varies. There is also mutual transformation among dislocations. It can be observed from Fig. 9c that when the indentation depth reaches approximately 4 nm at a grain size of 6.9 nm, other types of dislocations gradually transform into 12<111>, <110>, and <100> dislocations. Similarly, in Fig. 9d, it can also be clearly seen that 12<111> dislocations are mutually transformed with other types of dislocations at a grain size of 6.6 nm. In brief, dislocation evolution occurs in intermittent surges rather than as a monotonic process [19], which is attributed to plastic deformation and elastic recovery. In addition, dislocations already exist before the indentation loading begins, and they do not disappear completely at the end of the unloading process. However, by the end of unloading, the dislocation density returns to its initial state. Fig. 10 shows the grain-size dependence of dislocation-type proportions, consistent with the relative hardness/force ranking across grain sizes. It indicates that plastic deformation is mainly attributed to the formation of 12<111> dislocations, whereas the contributions of <110> and <100> dislocations are very small. This conclusion can also be observed through the dislocation maps, which are illustrated in Fig. 15. As the nanoindentation process proceeds, transformations and reactions occur between other types of dislocations and 12<111>, <110>, as well as <100> dislocations. Thus, their proportions fluctuate with the progress of the indentation process.

Figure 15: The dislocation maps of the W matrix at an indentation depth of 2.3 nm for different grain sizes. For the convenience of observation, only the dislocations are retained. (a) The dislocation map of grain size at 7.9 nm; (b) 7.3 nm; (c) 6.9 nm; (d) 6.6 nm; (e) 6.4 nm; (f) 4.8 nm.

The dislocation structures in Figs. 15 and 16 were obtained using the Dislocation Extraction Algorithm (DXA) method in OVITO. By comparing the six groups of dislocation images with different grain sizes in Fig. 15, it can be observed that different types of dislocations are interconnected to form nodes and junction networks [38]. Among them, green represents 12<111>, blue represents <110>, purple-red represents <100>, and red represents dislocation types other than the above three. Some planes are connected through dislocations, whereas some planes are relatively independent. There are various angles between each plane. With decreasing grain size, the number of dislocations declines, the interconnected dislocation networks shrink in area, and the spacing between individual dislocations increases, leading to a more scattered distribution. At different grain sizes, the dislocation planes formed by the interconnection of dislocation curves are mostly concentrated near the -XOZ plane.

Figure 16: (a) The dislocation maps of the W matrix at a grain size of 7.9 nm; (b–d) Magnified local images of different areas at an indentation depth of 2.3 nm and a grain size of 7.9 nm; (e) Dislocation maps of the W matrix at a grain size of 7.3 nm and an indentation depth of 2.3 nm; (f) Enlarged local view of (e).

For better observation, a crystal with a grain size of 7.9 nm is selected, and the highlighted region in Fig. 15a is locally magnified as shown in Fig. 16a. In Fig. 16c, the dislocation [111] and the adjacent dislocation 12[−111] connect two nodes, linking the two dislocation planes at a specific angle, which are plane 1 and plane 2 in Fig. 16a. The dislocation [111] and 12[−111] form a dislocation loop, indicated by the white dotted line in Fig. 16d. Within this dislocation loop, segments parallel to the Burgers vector 12<111> coexist with non-parallel segments. According to the definition of different types of dislocation loops, this dislocation loop is a mixed dislocation loop. Meanwhile, the dislocation [111] also forms another mixed dislocation loop with 12[111], indicated by the gray dotted line part in Fig. 16d. Furthermore, interconnected dislocation looping [38] is observed on Plane 1, where the dislocations interact, intersect, and link together, as shown in Fig. 16b.

Fig. 16e demonstrates a view of the -XOZ plane in a crystal with a grain size of 7.3 nm. It can be observed that different dislocations form planes at various angles to each other; however, unlike the case at a grain size of 7.9 nm, these planes are not interconnected by certain dislocation lines. On the contrary, planes 3 and 4 are nearly parallel. Plane 5 is almost entirely composed of the same dislocation line, namely dislocation 12[1−1−1]. By magnifying the lower half locally, it can be observed that the dislocation lines advance layer by layer along the crystal plane and eventually form a stepped, closed-loop dislocation loop structure, known as closed stair-rod dislocation looping [38]. For the W crystals with grain sizes of 6.9, 6.6, 6.4, and 4.8 nm, the dislocation plane and dislocation loop structures are similar to those at 7.9 nm.

4  Conclusion

This work clarifies the grain-size effect on nanoindentation responses of tungsten and the associated dislocation evolution by MD simulation. The main conclusions are:

(1)   A phenomenological transition around ~7.3 nm is observed within the tested grain-size range. Below this value, increasing grain size raises elastic recovery, load, hardness, and local stress/strain, while reducing pile-up atoms/area; above this value the trends weaken or reverse.

(2)   Indentation produces asymmetric pile-up and stress/strain localization governed mainly by grain-boundary configuration, with grain size modulating the magnitude.

(3)   The load–depth curve follows a Hertz-type power-law trend with intermittent drops due to dislocation nucleation; similar behavior appears in hardness–depth curves.

(4)   Dislocation density evolves dynamically through nucleation, interaction, and transformation, with 12<111>. dislocations dominating plastic deformation.

A quantitative comparison with experimental nanoindentation data, as well as a systematic assessment of temperature and loading-rate effects, would further strengthen the validation of these results and will be pursued in future work. Extending the same analysis framework to other BCC metals will help assess the transferability of the grain-size trends reported here.

Acknowledgement: Not applicable.

Funding Statement: This work was supported by the National Natural Science Foundation of China (52405585) and the Postdoctoral Fellowship Program of CPSF under Grant Number GZC20242223.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Duo Li, Huan Liu; methodology, Shuhao Kang; software, Yang Shen; formal analysis, Duo Li; investigation, Yukun Liu, Xin Wu; resources, Huan Liu; data curation, Shuhao Kang; writing—original draft preparation, Duo Li, Shuhao Kang, Yuhu Liu; writing—review and editing, Shujun Huang, Xin Wu, Huan Liu; visualization, Ruihan Li; project administration, Shujun Huang; funding acquisition, Huan Liu. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the corresponding author upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: Huan Liu reports financial support was provided by the National Natural Science Foundation of China. Huan Liu reports financial support was provided by the Postdoctoral Fellowship Program of CPSF. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

1. Mack CA. Fifty years of Moore’s law. IEEE Trans Semicond Manuf. 2011;24(2):202–7. doi:10.1109/TSM.2010.2096437. [Google Scholar] [CrossRef]

2. Dennard RH, Gaensslen FH, Yu HN, Rideout VL, Bassous E, LeBlanc AR. Design of ion-implanted MOSFET’s with very small physical dimensions. IEEE J Solid State Circuits. 1974;9(5):256–68. doi:10.1109/JSSC.1974.1050511. [Google Scholar] [CrossRef]

3. Chau R, Datta S, Doczy M, Doyle B, Kavalieros J, Metz M. High-κ/metal-gate stack and its MOSFET characteristics. IEEE Electron Device Lett. 2004;25(6):408–10. doi:10.1109/LED.2004.828570. [Google Scholar] [CrossRef]

4. Cecil T, Peng D, Abrams D, Osher SJ, Yablonovitch E. Advances in inverse lithography. ACS Photonics. 2022;10(4):910–8. doi:10.1021/acsphotonics.2c01026. [Google Scholar] [CrossRef]

5. Kaestner M, Krivoshapkina Y, Rangelow IW. Next generation lithography—the rise of unconventional methods? Front Nanosci. 2016;11:479–95. doi:10.1016/b978-0-08-100354-1.00014-4. [Google Scholar] [CrossRef]

6. Baglin JEE, Harrison HB, Tandon JL, Williams JS. Chapter 11—contacts and interconnections on semiconductors. In: Williams JS, Poate JM, editors. Ion implantation and beam processing. Cambridge, MA, USA: Academic Press; 1984. p. 357–409. [Google Scholar]

7. Choi D, Barmak K. On the potential of tungsten as next-generation semiconductor interconnects. Electron Mater Lett. 2017;13(5):449–56. doi:10.1007/s13391-017-1610-5. [Google Scholar] [CrossRef]

8. Mouroux A, Zhang SL, Petersson CS, Palmans R, Maex K, Ahlgren T, et al. An experimental study of chemical vapour deposition of tungsten on Ti/TiN adhesion bilayers: mechanical properties. Surf Coat Technol. 1998;99(3):274–80. doi:10.1016/S0257-8972(97)00567-7. [Google Scholar] [CrossRef]

9. Lee M, Hidayat R, Nandi DK, Kim TH, Kim Y, Kim S, et al. Atomic layer deposition of tungsten and tungsten-based compounds using WCl5 and various reactants selected by density functional theory. Appl Surf Sci. 2021;563:150373. doi:10.1016/j.apsusc.2021.150373. [Google Scholar] [CrossRef]

10. Lee J, Zhang C, Shon M, Read J, Deng S, Phadke O, et al. BEOL-compatible non-volatile capacitive synapse with ALD W-doped In2O3 semiconductor layer. In: 2024 IEEE International Electron Devices Meeting (IEDM); 2024 Dec 7–11; San Francisco, CA, USA. p. 1–4. doi:10.1109/IEDM50854.2024.10873481. [Google Scholar] [CrossRef]

11. Bringa EM, Caro A, Wang Y, Victoria M, McNaney JM, Remington BA, et al. Ultrahigh strength in nanocrystalline materials under shock loading. Science. 2005;309(5742):1838–41. doi:10.1126/science.1116723. [Google Scholar] [PubMed] [CrossRef]

12. Hasnaoui A, Van Swygenhoven H, Derlet PM. Cooperative processes during plastic deformation in nanocrystalline FCC metals: a molecular dynamics simulation. Phys Rev B. 2002;66(18):184112. doi:10.1103/physrevb.66.184112. [Google Scholar] [CrossRef]

13. Zepeda-Ruiz LA, Stukowski A, Oppelstrup T, Bertin N, Barton NR, Freitas R, et al. Atomistic insights into metal hardening. Nat Mater. 2021;20(3):315–20. doi:10.1038/s41563-020-00815-1. [Google Scholar] [PubMed] [CrossRef]

14. Guo X, Gou Y, Dong Z, Yuan S, Li M, Du W, et al. Study on subsurface layer of nano-cutting single crystal tungsten in different crystal orientations. Appl Surf Sci. 2020;526(7088):146608. doi:10.1016/j.apsusc.2020.146608. [Google Scholar] [CrossRef]

15. Wang H, Dong Z, Wang C, Kang R, Guo X, Bao Y. Surface/subsurface formation mechanism of tungsten during ultrasonic elliptical vibration cutting. Int J Mech Sci. 2024;262:108725. doi:10.1016/j.ijmecsci.2023.108725. [Google Scholar] [CrossRef]

16. Dong Z, Wang H, Qi Y, Guo X, Kang R, Bao Y. Effects of minimum uncut chip thickness on tungsten nano-cutting mechanism. Int J Mech Sci. 2023;237(5):107790. doi:10.1016/j.ijmecsci.2022.107790. [Google Scholar] [CrossRef]

17. Wang H, Dong Z, Yuan S, Guo X, Kang R, Bao Y. Effects of tool geometry on tungsten removal behavior during nano-cutting. Int J Mech Sci. 2022;225:107384. doi:10.1016/j.ijmecsci.2022.107384. [Google Scholar] [CrossRef]

18. Valencia FJ, Ortega R, González RI, Bringa EM, Kiwi M, Ruestes CJ. Nanoindentation of nanoporous tungsten: a molecular dynamics approach. Comput Mater Sci. 2022;209(1):111336. doi:10.1016/j.commatsci.2022.111336. [Google Scholar] [CrossRef]

19. Lin PD, Nie JF, Cui SG, Cui WD, He L, Xiao GY, et al. Exploring the impact of pre-existing helium bubbles on nanoindentation in tungsten through molecular dynamics simulation. J Mater Res Technol. 2024;31(1):2708–22. doi:10.1016/j.jmrt.2024.06.222. [Google Scholar] [CrossRef]

20. Wang BB, Xie GM, Wu LH, Xue P, Ni DR, Xiao BL, et al. Grain size effect on tensile deformation behaviors of pure aluminum. Mater Sci Eng A. 2021;820:141504. doi:10.1016/j.msea.2021.141504. [Google Scholar] [CrossRef]

21. Shi Z, Meng Y. Effects of indentation depth and grain size on scratching behavior of nanograin FCC Fe polycrystalline substrate. Tribol Int. 2024;193:109464. doi:10.1016/j.triboint.2024.109464. [Google Scholar] [CrossRef]

22. Cao H, Guo Z, Feng R, Li H, Fu R, Zhou Y, et al. Numerical simulation of nano-cutting behaviors for polycrystalline γ-TiAl alloy: the effect of grain sizes. J Manuf Process. 2023;102(5):169–81. doi:10.1016/j.jmapro.2023.07.047. [Google Scholar] [CrossRef]

23. Papanikolaou M, Salonitis K. Grain size effects on nanocutting behaviour modelling based on molecular dynamics simulations. Appl Surf Sci. 2021;540(9):148291. doi:10.1016/j.apsusc.2020.148291. [Google Scholar] [CrossRef]

24. Hachet G, Caillard D, Ventelon L, Clouet E. Mobility of screw dislocation in BCC tungsten at high temperature in presence of carbon. Acta Mater. 2022;222(11):117440. doi:10.1016/j.actamat.2021.117440. [Google Scholar] [CrossRef]

25. Thompson AP, Aktulga HM, Berger R, Bolintineanu DS, Brown WM, Crozier PS, et al. LAMMPS—a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comput Phys Commun. 2022;271(4):108171. doi:10.1016/j.cpc.2021.108171. [Google Scholar] [CrossRef]

26. Hirel P. Atomsk: a tool for manipulating and converting atomic data files. Comput Phys Commun. 2015;197:212–9. doi:10.1016/j.cpc.2015.07.012. [Google Scholar] [CrossRef]

27. Nakamura K, Kumagai T, Ohnuma T. Atomistic simulation of shear deformation at bcc-Fe grain boundary and precipitation strengthening by Cr23C6. Mater Today Commun. 2022;33:104711. doi:10.1016/j.mtcomm.2022.104711. [Google Scholar] [CrossRef]

28. Du X, Zhao H, Zhang L, Yang Y, Xu H, Fu H, et al. Molecular dynamics investigations of mechanical behaviours in monocrystalline silicon due to nanoindentation at cryogenic temperatures and room temperature. Sci Rep. 2015;5(1):16275. doi:10.1038/srep16275. [Google Scholar] [PubMed] [CrossRef]

29. Tersoff J. Modeling solid-state chemistry: interatomic potentials for multicomponent systems. Phys Rev B. 1989;39(8):5566–8. doi:10.1103/physrevb.39.5566. [Google Scholar] [PubMed] [CrossRef]

30. Farkas D, Caro A. Model interatomic potentials for Fe-Ni–Cr-Co–Al high-entropy alloys. J Mater Res. 2020;35(22):3031–40. doi:10.1557/jmr.2020.294. [Google Scholar] [CrossRef]

31. Singh JK, Adhikari J, Kwak SK. Vapor-liquid phase coexistence curves for Morse fluids. Fluid Phase Equilib. 2006;248(1):1–6. doi:10.1016/j.fluid.2006.07.010. [Google Scholar] [CrossRef]

32. Li D, Zhang Z, Long C, Li R, Kang S, Li Z, et al. Molecular dynamics of the nanogrinding process for single-crystal tungsten: investigation of the regulatory role of grinding speed. Int J Refract Met Hard Mater. 2026;134:107491. doi:10.1016/j.ijrmhm.2025.107491. [Google Scholar] [CrossRef]

33. Wang H, Kang R, Bao Y, Wang K, Guo X, Dong Z. Microstructure evolution mechanism of tungsten induced by ultrasonic elliptical vibration cutting at atomic/nano scale. Int J Mech Sci. 2023;253:108397. doi:10.1016/j.ijmecsci.2023.108397. [Google Scholar] [CrossRef]

34. Stukowski A. Visualization and analysis of atomistic simulation data with OVITO—the open visualization tool. Model Simul Mater Sci Eng. 2010;18(1):015012. doi:10.1088/0965-0393/18/1/015012. [Google Scholar] [CrossRef]

35. Li J, Guo J, Luo H, Fang Q, Wu H, Zhang L, et al. Study of nanoindentation mechanical response of nanocrystalline structures using molecular dynamics simulations. Appl Surf Sci. 2016;364:190–200. doi:10.1016/j.apsusc.2015.12.145. [Google Scholar] [CrossRef]

36. Li J, Lu B, Zhang Y, Zhou H, Hu G, Xia R. Nanoindentation response of nanocrystalline copper via molecular dynamics: grain-size effect. Mater Chem Phys. 2020;241:122391. doi:10.1016/j.matchemphys.2019.122391. [Google Scholar] [CrossRef]

37. Zhou J, Jiao Z, Zhang J, Zhong Z. Nanoindentation of single-crystal and polycrystalline yttria-stabilized zirconia: a comparative study by experiments and molecular dynamics simulations. J Alloys Compd. 2021;878(5566):160336. doi:10.1016/j.jallcom.2021.160336. [Google Scholar] [CrossRef]

38. Shinde AB, Owhal A, Sharma A, Ranjan P, Roy T, Balasubramaniam R. Comparative analysis of mechanical properties for mono and poly-crystalline copper under nanoindentation-Insights from molecular dynamics simulations. Mater Chem Phys. 2022;277(80):125559. doi:10.1016/j.matchemphys.2021.125559. [Google Scholar] [CrossRef]

39. Cao H, Yu Z, Zhou B, Li H, Guo Z, Wang J, et al. Molecular dynamics simulation of the effect of supersonic fine particle bombardment on the mechanical behaviour of γ-TiAl alloy: the effect of grain size. Vacuum. 2023;217:112498. doi:10.1016/j.vacuum.2023.112498. [Google Scholar] [CrossRef]

40. Fan J, Li J, Huang Z, Wen PH, Bailey CG. Grain size effects on indentation-induced plastic deformation and amorphization process of polycrystalline silicon. Comput Mater Sci. 2018;144:113–9. doi:10.1016/j.commatsci.2017.12.017. [Google Scholar] [CrossRef]