Grain size effects on nanocutting behaviour modelling based on molecular dynamics simulations

Grain size is one of the most critical factors affecting the mechanical and thermal properties of metallic materials. In this study the effect of the grain size of a workpiece made of pure Aluminium on the nanocutting process has been investigated via means of Molecular Dynamics simulations. The polycrystalline workpiece has been generated starting from a Face-Centred Cubic block of Aluminium atoms which was melted and subsequently quenched under various cooling rates in order to control the average grain size. The case of a monocrystalline workpiece has been considered as well. The generated workpieces were ground by a diamond abrasive. Simulations have been repeated in order to eliminate any statistical errors. The obtained results suggest that the average grain size of the workpiece significantly influences almost every aspect related to the nanocutting process. More specifically, it has been found that the cutting forces increase and the friction coefficient decreases with the grain size. Very small grain sizes lead to lower thermal conductivity and consequently high temperature at the cutting region. Finally, it has been shown that the high residual stresses at the grain boundaries can be relieved as the tool passes on top of the workpiece; this phenomenon resembles heat treatment. In summary, the nanocutting behaviour of polycrystalline materials depends on the average grain size and significantly differs from the case of monocrystalline materials. This should be taken into account in future numerical models of nanocutting processes.


Introduction
Since the early 1950s, there has been clear scientific evidence that the average grain size plays a significant role in the mechanical properties of nanocrystalline materials. According to the empirical Hall and Petch relationship [1,2] the yield strength is inversely proportional to the grain size whereas this trend appears to change for nanocrystalline materials with a grain diameter below a critical value [3]. Although the effect of the grain size on the mechanical properties of bulk and nanocrystalline materials has been extensively investigated by both experimental and theoretical/numerical studies [4,5], there are only a few experimental studies discussing the effects of the workpiece grain size on the cutting characteristics presenting occasionally conflicting results. For instance, in [6] the authors have found that smaller grain size increases the cutting forces while in [7] the opposite conclusion was drawn. This discrepancy is attributed to the different grain sizes and materials used in the aforementioned studies. Thus it can be concluded that experimental results on the topic under investigation are quite sensitive to the process/material parameters adopted.
The rapid evolution of computational power over the past decades has allowed for the investigation of phenomena taking place in material removal processes over a wide range of length scales using numerical modelling tools. Some of the most common numerical modelling methods used for modelling material removal processes are Kinematic Models, Finite Element Modelling (FEM) and Molecular Dynamics (MD) simulations [8]. Although Kinematic Models and FEM can be applied in both micro-and macroscale, they are based on constitutive relations and assumptions such as a moving heat source [9] or pre-estimates of the cutting forces and the friction coefficient [10]. In contrast to the aforementioned modelling techniques, MD simulations are based on interatomic potentials derived from experimental data or quantum mechanical calculations [11] while no further assumptions need to be made. As a consequence, the high computational cost of MD is offset by its high resolution and thus has been used for modelling a wide variety of phenomena including nanoflows [12], fluid [13,14] and material properties [15]. Since the current investigation is focused on nanoscale phenomena occurring in nanocutting, MD has been selected as the most suitable modelling technique among the available alternatives.
Moreover, MD is suitable for modelling the nanocutting process of polycrystalline materials as no assumptions need to be made on the dependence of material properties on the average grain size.
There is a vast body of literature on modelling nanocutting processes via means of MD simulations. One of the most extensively topics studied is the effect of the tool geometry on grinding process characteristics [16][17][18]. These studies were mainly concerned with the effect of the tool geometry (rake angle, clearance angle and tool edge radius) on the cutting forces, cutting specific energy and subsurface temperature field induced. The high resolution of MD simulations has also been exploited to investigate subsurface damage evolution during nanocutting. Li et al. [19] quantified subsurface damage by evaluating the population of the Intrinsic Stacking Fault (ISF) and amorphous atoms over the course of their simulations. They concluded that low values of the depth of cut and the tool radius as well as high grinding speeds reduce the extent of subsurface damage. In another study, Wang et al. [20] adopted the Centro-Symmetry Parameter (CSP) and the spherical harmonic method (SHM) in order to characterise subsurface defects and the local atomic structure at the ground workpiece region. The authors visualised and analysed the structural evolution of the workpiece during the cutting process by monitoring the population of atoms being part of subsurface defect structures (dislocation/vacancy/atomic cluster defects, stacking faults and stair-rod dislocations). MD simulations have also been employed to investigate the effects of many other nanocutting process parameters such as grinding speed [21], fractally rough tool-workpiece interface [22] and wheel binder stiffness [23] on both thermal and mechanical behaviour of the components, i.e. workpiece and abrasive grit.
A shared characteristic between the aforementioned MD studies is that the workpiece is considered as a monocrystalline block of atoms placed on the sites of a crystal lattice. However, this is an assumption that does not hold in reality as metals exhibit a polycrystalline grain structure while only a few MD studies have considered the scenario of a polycrystalline workpiece. One of the first MD investigations focusing on the effects of polycrystalline grain structure on the nanocutting mechanism was performed by Venkatachalam et al. [24] in 2015. The authors attempted to describe the effects of grain size, grain boundaries and crystallographic orientation on the flow stress of the polycrystalline workpiece and the cutting forces via means of both experiments and MD simulations. Their 2-dimensional MD model consisted of a workpiece containing only 2 grains with different crystallographic orientation and it was observed that the cutting forces as well as deformation were affected by the lattice orientation. In another study, Liu et al. [25] considered the effect of polycrystalline copper on the formation mechanism of grain boundary steps. Their copper workpiece consisted initially of 2 grains with different orientations. The authors observed that the lattice dislocations generated at the interface between the cutting tool and the workpiece were piled up at the grain boundary between the initial 2 grains; this subsequently lead to the formation of sub-grains at this area in order for the system energy to be minimised. Moreover, it was observed that the cutting forces reached a peak as the tool passed on top of the grain boundary and this was attributed to hardening. Shi et al. [26] investigated the effect of the grain size of a copper workpiece and some machining parameters (depth of cut, grinding speed and tool rake angle) on the cutting forces and the stress/dislocation distribution over the workpiece. In their model grains were represented as hexahedral patterns while the polycrystalline workpiece geometry was prepared by adopting the Voronoi site-rotation and cut package of AtomEye [27]. Liu et al. [28] employed the same model in order to build a polycrystalline silicon carbide workpiece and investigated the material removal mechanisms and the dependence of the cutting forces and temperature on the grinding speed and depth of cut. However, the aforementioned studies employ patterned grains which are not encountered in nature often; this is an assumption which might affect the accuracy of the obtained results.
In the current investigation the effects of the grain size of a polycrystalline Aluminium (Al) workpiece on the nanocutting process characteristics will be presented and discussed. In contrast to the previously mentioned investigations, in which the polycrystalline structure was generated using the Voronoi site-rotation and cut package of AtomEye, in this research the polycrystalline geometry was built by melting and subsequently quenching a block of Aluminium atoms in order to obtain a realistic grain distribution over the simulation domain. Workpiece geometries with various average grain sizes were constructed by tuning the cooling rate while the case of a monocrystalline workpiece was considered as well. The MD simulations performed were also repeated in favour of statistical accuracy and consistency. The obtained results are focused on the cutting forces, friction coefficient, subsurface temperature over the simulation domain which appear to be highly dependent on the average grain size. Last but not least, it is shown that the average grain size affects the distribution of stress and temperature at the workpiece subsurface. More specifically, it has been shown that the grain size plays a significant role in the morphology of the shear zone as well as the heat dissipation over the workpiece.

Methodology
In the present investigation the workpiece topology was generated starting from a Face-Centred Cubic (FCC) block of Aluminium which was melted and subsequently quenched under various cooling rates in order to obtain different grain sizes. The obtained Al structure was then used as a workpiece in the grinding simulations. In this section, the steps employed for setting up the simulations will be presented in the following two subsections (2.1 & 2.2). In the first one the methodology followed for generating the workpiece grain structure will be presented while the latter one will be focused on the setup of the nanocutting simulations.

Melting and solidification
The starting point of this research was a rectangular block of 480,000 Al atoms with dimensions 80 × 6 × 16 nm 3 arranged in the sites of a FCC lattice with a lattice constant equal to 4.046 Å. The Finnis-Sinclair (FS) potential [29,30], which has been widely used in solidification MD studies [31,32], was employed to model the interactions between the Al atoms. The potential energy of a system of N atoms (E tot ) is estimated as follows: where V ij (r ij ) is the potential energy between the ij pair of atoms as a function of their interatomic distance r ij , A the binding energy and ρ i the local electronic charge in the vicinity of an atom i which can be calculated as follows: In Eq. (2) φ ij is the atomic charge density, β a parameter used to introduce a maximum in the atomic charge density φ ij and d the cutoff distance.
The methodology followed to generate various workpiece grain structures is similar to the one presented in [33]. The simulations were performed using the isothermal-isobaric ensemble (NPT) and a Nose-Hoover thermostat and barostat were used for controlling the temperature and pressure of the simulation domain respectively. In the beginning, the Al block was equilibrated at 273 K and then isothermally annealed to 1173 K in order to eliminate all FCC structures. The melt was subsequently equilibrated and quenched under 3 different cooling rates (0.5, 2 and 8 K/ps) in order to obtain various grain structures. For each cooling rate the simulations were repeated 3 times under different initial conditions and 3 workpiece topologies were obtained (9 in total); this was to enhance the statistical accuracy of this investigation as it will be discussed in the Results section. As shown in previous investigations, higher cooling rates lead to lower average grain size [33,34]. The temperature and pressure damping parameters during quenching were equal to 100 and 1000 timesteps (0.2 and 2 ps) respectively. Finally, the simulation domain was equilibrated for another 20,000 timesteps until temperature and pressure were stabilised. The melting and solidification simulation parameters are summarised in Table 1.
In Fig. 1 the workpiece topologies corresponding to 3 different cooling rates are illustrated. It can be observed that the cooling rate significantly affects the workpiece grain structure. Low cooling rates lead to large but few grains while on the contrary, the higher the cooling rate the more the grains (or the smaller the average grain size). It has to be noted that due to the periodic boundary conditions applied, some grains extend through the periodic boundaries; this is the reason why some red grains appear to be smaller in size compared to some green ones. In the case of c R = 0.5 K/ps the obtained grains are fewer and larger compared to the other two. The obtained grain structures were used as an input for the grinding simulations as it will be further elaborated in the following subsection. The main advantage of the method employed to generate the atom configuration compared to Voronoi-site rotation [35] is that the final grain structure (average grain size, grain size distribution and grain orientation) is dependent on the cooling rate and resembles natural patterns to a greater extent.
Finally, it has to be noted the various grain structures generated via the current approach are characterised by intrinsic defects including lattice dislocations, crystal twins and stacking faults. In order to provide an overview of the crystal structure corresponding to each cooling rate the population of Face Centred Cubic (FCC), boundary or amorphous (AMO) and Hexagonal Close Packed (HCP) atoms, the latter being indicative of the initial defect content of the workpiece, was plotted in Fig. 2 as a function of the cooling rate. The population of each atom type was estimated at the final timestep of the equilibration process following solidification. As illustrated in Fig. 2 the boundary atoms increase with the cooling rate; this is because higher cooling rates lead to additional grains and consequently increased area of the grain boundaries. The opposite observation can be made for the FCC atoms for the exact same reason; increased grain boundaries imply reduced volume of the FCC phase. Finally, as far as HCP atoms are concerned, their population appears to increase for cooling rates lower than 4 K/ps and subsequently decrease. This is due to the superposition of two phenomena: (a) higher cooling rates promote the formation of lattice defects [36] and (b) higher cooling rates hinder the formation of the FCC and consequently the HCP phase [34]. As advocated by Li et al. [37] the presence of defects is expected to significantly affect the strength-ductility trade-off of the workpiece.

Nanocutting
The nanocutting MD model was consisted of one of the solidified workpiece topologies and a cutting tool made of diamond as illustrated in Fig. 3. The cutting tool was assumed to have a negative rake angle which is quite common in nanocutting simulations using MD [38] while the selected rake angle (-45 ) is much below the critical value for chip formation (≈-60 ) [39]. Finally, the width of the cutting tool (z-direction) was selected so as to make sure that the material pile-up in the zdirection does not interfere with the periodic boundary conditions imposed on the z-axis, as this would alter the system dynamics (Fig. 7).   The interaction between the workpiece atoms (Al-Al) was modelled using the Finnis-Sinclair (FS), potential as described in the previous subsection. Besides solidification, FS potentials have been used in the past for simulating nanocutting processes and proven to be capable of accurately estimating the stacking fault energy [40,41]. The Morse potential [42] was used to model the pair interaction energy (E) between the diamond (C) and workpiece (Al) atoms: The parameters of the Morse potential, D 0 and α were set equal to 0.28 eV and 2.78 Å −1 respectively according to [43]. Moreover, the equilibrium distance was set equal to 2.2 Å according to [44] while the cutoff was considered equal to 6 Å. Although the interaction between the C atoms of the tool has been occasionally neglected due to the superior hardness of diamond compared to other metals [45], in this work it was described by the Tersoff potential [46] in order to enhance accuracy. The Tersoff potential has been frequently used for modelling the structural properties of diamond [23,28] and is given by: where f c (r ij ) a cutoff function, f R (r ij ) and f A (r ij ) the repulsive and attractive components of the potential respectively and b ij a monotonically decreasing function of the coordination of atoms i and j. The simulation setup contained 3 atom types as illustrated in Fig. 3, namely (a) rigid, (b) thermostat and (c) Newtonian atoms. The rigid atoms of the workpiece were used for holding the workpiece still during the cutting process whereas the rigid atoms of the cutting tool were assigned a velocity equal to 100 m/s in the x-direction following an initial equilibration; these atoms were allowed to move in the x-direction but their relative position was constant and they were not time integrated. A Langevin thermostat was applied to the thermostat atoms in order to maintain the temperature of this region at 300 K. During the initial equilibration process, which took place over the first 10,000 timesteps of the nanocutting simulation, the temperature of all the simulation atoms reached the desired temperature (300 K) due to conduction. The thickness of both the rigid and thermostat layers was equal to 8 Å. In contrast to the two previous atom types, the Newtonian atoms were time integrated via the microcanonical ensemble (NVE) and their trajectories were dictated by the Newton's second law of motion. In contrast to the solidification simulation, periodic boundary conditions were applied solely on the z-direction while fixed boundary conditions were applied to the other two directions. The depth of cut (d c ) was set at 2 nm and the timestep was equal to 2 ps. The nanocutting simulation parameters are listed in Table 2. The MD simulations of this study were run on 16 cores using the LAMMPS Molecular Dynamics simulator [47] and the running time was about 47 h per nanocutting simulation. The running time of the solidification simulations ranged between 6 and 40 h depending on the cooling rate applied. The Ovito open source software [48] was used for visualisation and partially for post-processing.

Results & discussion
The obtained results will be presented in three sections. In the first one (3.1) the effects of the grain structure on the cutting forces and friction coefficient will be discussed. Subsequently the spotlight will be shifted towards the dependence of subsurface temperature (3.2) and stresses (3.3) on the average grain size both visually and quantitively. The average grain size was calculated in terms of grain diameter (d g ) and by assuming that grains are spherical according to the following relationship: where N ave the average number of Al atoms contained by the system grains and V a the per atom volume of workpiece Al atoms which was calculated using Voronoi tessellation. The effect of the cooling rate on the average grain size is presented in the errorbar of Fig. 4. The range of the errorbars is due to the 3 repetitions performed for each case in order  to enhance statistical accuracy. As expected, the higher the cooling rate the smaller the average grain size; the results obtained are consistent with previous studies [33,49]. It has to be noted that the calculated average grain size values seem to be small when compared to Fig. 1. This is due of the fact there are numerous very small grains lying at the interstices of the larger ones and consequently downgrade this value. In the following subsections, for ease of understanding and to avoid the use of ranges, the term average grain size d g,ave will refer to the following quantity: where i corresponds to the repetition index (max = 3). The average grain size values considered in this study will be 4.49, 2.58 and 1.98 nm, corresponding to cooling rates of 0.5, 2 and 8 K/ps respectively as shown in Fig. 4. The case of a monocrystalline workpiece will be considered as well.
Besides the average grain size, the cooling rate also affects the grain size distribution across the workpiece as illustrated in Fig. 5. Moreover, the properties of a polycrystal are not dependent solely on the average grain size but also to the grain size distribution. By comparing Fig. 5(a), (b) and (c) it can be observed that faster quenching leads to fewer coarse grains and additional finer ones as expected. Moreover, it can be seen that for a cooling rate of 0.5 K/ps the largest grain encountered corresponds to d g ≈15.28 nm while in the case of c R = 8 K/ps the largest grain encountered is less than half in size (d g ≈6.85 nm).

Cutting forces and friction coefficient
The evolution of the cutting forces over time is illustrated in Fig. 6. It can be observed that for all the graphs there is a linear increase of both the cutting forces in the time interval between 0 and 300 ps; this corresponds to the period when the cutting tool penetrates the workpiece. Following this timespan, the tangential (F x ) forces appear to increase over the course of the simulation; this is due to the material pilling-up at the front of the cutting tool as shown in Fig. 7. From Fig. 6 it can be observed that the cutting forces exhibit large fluctuations. These are attributed to 2 main reasons, namely (a) the stick slip phenomenon [25] and (b) the grain structure in the vicinity of the workpiece; as advocated by Liu et al. [25] peaks at the cutting forces can be observed when the tool passes on top of a grain boundary.
Another observation that can be extracted from Fig. 6 is that both tangential (F x ) and normal (F y ) forces seem to decrease on average in the time interval between 300 and 800 ps as the average grain size decreases; this is more pronounced for F y . Moreover, it can be observed that in the monocrystalline workpiece case, F x is significantly higher than F y while this trend appears to become inverted as the average grain size decreases. This is expected to affect the value of the friction coefficient η = F x /F y and will be further discussed in the following lines.
Although Fig. 6 provides significant information on the evolution of the cutting forces over time for a single set of simulations the quite large fluctuation of the forces does not ease the extraction of safe conclusions. In order to obtain a more systematic and cumulative overview of the effect of the average grain size on the cutting forces, the mean values of F x and F y were obtained over the timespan [300 ps, 800 ps]. The calculated means were averaged over the 3 sets of simulations in a manner similar to the calculation of the average grain size (Eq. (6)) and represented in the bar charts of Fig. 8 and Fig. 9.
A conclusion that can be drawn by observing Fig. 8 and Fig. 9 is that smaller grains lead to the reduction of the cutting forces. This is due to grain boundary sliding which can be observed for very small grain sizes and deteriorate the material strength. Grain boundary sliding has been proven to play a major role in the deformation of nanocrystalline materials with fine grains even at room temperature [50]. Moreover, according to [51] cutting forces are proportional to the shear strength of materials which, for FCC metals, is about 30 times lower than their shear modulus [52]. Our results are in accordance with the ones presented in [26]. However, it has to be noted that the grain size seems to have a more dominant effect on the tangential forces rather on the normal ones. This observation is associated with the thermal softening of the material at the vicinity of the tool and will be further elaborated in the following section.
As discussed above, the grain structure does not affect both F x and F y in the same manner. This makes the friction coefficient η dependent on the grain structure as illustrated in Fig. 10. In the cases of a monocrystalline workpiece and d g,ave = 4.49 nm the friction coefficient maintains a value below unity which resembles real material removal  applications. For further reduced values of d g,ave the friction coefficient becomes greater than 1 which however is not uncommon in nanocutting simulations and depends on various parameters including the material under consideration and the interface area between the tool and the workpiece [53].

Subsurface temperature
In order to study the evolution of the subsurface temperature the ensemble average temperature T of the workpiece atoms contained by a spherical region with a diameter of 1.5 nm and attached to the diamond tool was measured as follows: In Eq. (7) KE is the kinetic energy of the N the number of atoms contained in the spherical region, k B = 1.38•10 -23 kg s −2 K −1 the Boltzmann constant and dim the number of cartesian coordinates participating in the calculation of temperature. In this study only the zcomponent of the atoms' velocity (v zi ) has been used for the calculation of temperature; the other 2 components (v xi and v yi ) were excluded as they are associated to the translational motion of the cutting tool and would lead to miscalculation of temperature. Since only one component of the velocity vector was used dim was set equal to 1. In order to reduce noise the temperature of the workpiece atoms at the tool vicinity was time-averaged over time windows of 20,000 timesteps and the average temperature over the course of each simulation was calculated. Subsequently, these values were averaged for the 3 sets of simulations performed and plotted in Fig. 11.
As illustrated in Fig. 11 the workpiece temperature at the cutting region is dependent on the grain size: the lower the grain size the higher the temperature. It has to be noted that the subsurface temperature does not significantly increase for grains with d g,ave < 4.49 nm although there is a steep increase when moving from the monocrystalline case to d g,ave = 4.49 nm. Therefore it can be concluded that thermal softening is the reason behind the sharp decrease of the F x forces (Fig. 8) which appear to fall inversely to the rise of subsurface temperature at the cutting region. However, the reason behind the rise of temperature in the cutting area is still unclear.
In order to shed light to the aforementioned observation the per workpiece-atom temperature was calculated according to Eq. (7) and the 3-dimensional workpiece temperature profiles were obtained at identical timesteps for the 4 different cases under examination: (a) monocrystalline workpiece, (b) d g,ave = 4.49 nm, (c) d g,ave = 2.58 nm, and (d) d g,ave = 1.98 nm as illustrated in Fig. 12. In favour of clarity, a cross section at the midpoint of the workpiece width (z-direction) was used. It can be observed that the temperature in the regions 1(a-d) (T 1a-d ) is lower than in the corresponding regions 2(a-d) (T 2a-d ).This is a very interesting phenomenon and has been attributed to the effect of the grain size on the heat transfer mechanism. According to the kinetic theory the thermal conductivity (κ) of metals is proportional to the phonon mean free path (l) [54]: where C is the specific heat per unit volume and v RMS the root-mean-    squared velocity of phonons. The presence of grain boundaries leads to phonon scattering and consequently reduces the mean free path. As a result the thermal conductivity of metals reduces for smaller grain sizes [55]. The observations that T 1a < T 1b < T 1c < T 1d and T 2a > T 2b > T 2c > T 2d are attributed to the slower heat dissipation emerging from the lower thermal conductivity values for finer grain structures. Moreover, by observing Fig. 11 it can be concluded that thermal conductivity is negligibly affected when the average grain size becomes lower than 2.58 nm. This is the first MD investigation on nanocutting capturing this nanoscale phenomenon. Finally, the obstacles imposed on heat transfer imposed by the presence of grain boundaries lead to the higher temperature in the cutting area as observed in Fig. 12.

Stresses & defects
The principal and shear stress distribution during nanocutting has been illustrated for the monocrystalline case in Fig. 13 in terms of absolute values at a specific timestep. The cutting tool atoms have been made transparent and a cross section at the middle of the z-axis has been used to enhance clarity. With regard to the principal stresses, σ xx obtains high values at the front of the cutting tool due to the compression of this area. For the exact same reasons, σ yy obtains negative values at the front and below the cutting tool. However, in the case of σ zz the high stresses observed in the vicinity of the tool are a result of (a) the dislocation of the atoms to the sides of the tool and (b) the periodic boundary conditions applied on this direction which lead to compression. Moreover, it has to be noted σ zz is high only in the very close vicinity of the tool and lower in magnitude compared to the other two principal stresses; this is due to the fact that the magnitude of z-forces is much smaller compared to the other 2-directions [22]. With regard to the shear stresses, σ xy appears to be much lower compared to the principal ones as expected.
In the last part of this investigation the effect of the grain size on the distribution of Von Mises (VM) stresses over the simulation domain will be discussed. Equivalent Von Mises stress given by Eq. (9) was used to describe the stress distribution on the workpiece. In Eq. (9) σ xx, σ yy, σ zz the principal and σ xy, σ yz , σ zx the shear stresses respectively. In order to obtain 3-dimensional Von Mises stress profiles the per-atom stress tensor components were calculated and averaged over time windows of 20,000 timesteps in order to minimise noise.
In Fig. 14 the obtained VM stress profiles for the monocrystalline and polycrystalline workpiece morphologies have been present in parallel to the results of the Common Neighbour Analysis (CNA) at identical timesteps. With respect to the CNA results, green, white and red atoms represent Face-Centred Cubic (FCC), amorphous (AMO) and Hexagonal Closed Packed (HCP) structures respectively. The morphology of the shear zone of Fig. 14(a) is similar to FEA results for negative rake angle tools [56]. It has to be noted that the CNA of the corresponding timestep shows that there are almost no planar defects (or stacking faults) in the vicinity of the tool. However, the landscape changes completely when planar defects are present in the cutting area as shown in Fig. 14(b). In this case the obtained stress distribution is not uniform but VM stresses are concentrated in the vicinity of the dislocations/planar defects as expected. In the case of polycrystalline workpiece morphology the shear zone is even less uniform while it can be observed that there is high stress concentration at the grain boundaries and planar defects as shown in Fig. 14(c) and Fig. 14(d). It has to be noted that the stress concentration at the grain boundaries is not attributed only to cutting forces but also to the residual stresses emerging from solidification. The aforementioned observations are consistent with [26] and highlight the advantages of MD simulations over FEA regarding modelling nanoscale phenomena. FEA is not capable of capturing the planar defect generation or grain boundaries and thus the dependence of the shear zone   profile on the planar defects or grain boundaries is ignored.
The effects of stress concentration at the grain boundaries have been analysed with the help of the Dislocations Analysis (DXA) toolkit [57] of OVITO and presented in Fig. 15. A cross section at the midpoint of the workpiece width (z-axis) was used in favour of clarity while a transparent surface mesh was employed in order to make dislocations distinguishable. In Fig. 15(a) there are only a few dislocations observed at the bulk of the monocrystalline workpiece and thus, it can be concluded that the cutting mechanism for the monocrystalline workpiece is lattice fracture. On the other hand, in Fig. 15(b) the number of dislocations at the final timestep is significantly higher compared to the monocrystalline workpiece. Moreover, the vast majority of dislocations lie at the grain boundaries. This observation indicates that grain boundary sliding is one of the cutting mechanisms, is in accordance with previous investigations [26,28] and supports the claim that the inverse Hall-Petch relationship becomes dominant for very small grain sizes.
Although Fig. 15 provides a nice illustration of the dislocation distribution over the workpiece, it is of significant importance to quantify these observations. In Fig. 16 the population of various dislocation types segments at the uncut region of the workpiece has been plotted for the various workpiece grain topologies. It has to be noted that in the monocrystalline case the cutting direction (or crystal orientation) might affect defect formation as advocated by Guo et al. [58]. As shown in Fig. 16(a), which corresponds to the monocrystalline workpiece case, there are no dislocations before nanocutting, as expected. However, after the nanocutting process has been completed, Shockley and Stairrod dislocations are prevalent; this is in agreement with [59].
However, the dislocation distribution changes when grains are present in the workpiece. As shown in Fig. 16(b-d) other, perfect and Shockley dislocations are now prevalent. It should be highlighted that the dislocations before nanocutting are attributed to the solidification process. Be comparing Fig. 16(b-d) it can be observed the higher the cooling rate the more the lattice dislocations (before nanocutting), as expected. Moreover, in the case of d g,ave = 4.49 nm, it can be seen that almost all dislocation types segments (with the exception of Frank dislocations) in the uncut region increase after nanocutting. This trend is reversed for lower grain sizes. More specifically, the population of the various dislocation types appears to decrease after nanocutting. This phenomenon will be discussed and explained in the following paragraphs.
In the last part of this study the evolution of the residual stresses during the nanocutting process was investigated. For this purpose three regions (rectangular bins) were considered as shown in Fig. 17. The length of R1, R2 and R3 was 20 nm (x-direction) while the total length of the workpiece was 80 nm. Moreover, they extended all along the z-axis. In order to analyse the evolution of the residual stresses during nanocutting, the Von Mises stresses at the grain boundaries were summed for each one of the regions. Therefore, these regions were extended in the ydirection below the level of the cutting line in order to exclude the amorphous atoms of the chip/cutting region and provide misleading VM stress values. For the same reason the considered regions (R1-3) contained only Newtonian atoms; thermostat or rigid atoms would possibly provide a false estimation of the residual stresses. The current analysis excluded the monocrystalline workpiece scenario as there are no grain boundary atoms within the volumes covered by R1, R2 and R3. The obtained results are summarised in Fig. 18. The values of residual stresses at the first timestep are attributed to quenching while the values at the last timesteps correspond to the remaining residual stresses after the action of the cutting tool has been removed. A common observation between the graphs of Fig. 18 is that in all cases the total VM stress is higher for finer grain structures. This is because smaller grains increase the area of grain boundaries. Moreover, this is in accordance with experimental results suggesting that high cooling rates lead to high levels of residual stresses [60].
In Fig. 18, the value of the total VM stress at the initial averaging point of the curves corresponding to different regions but identical grain sizes varies. This is because σ VM is dependent on the local atomic structure which is not uniform across the x-direction of the workpiece (or cutting direction). However, a common trend between all curves of Fig. 18 can be observed: the total VM stress (a) starts from a specific value which corresponds to the residual stresses left after quenching, (b) VM stress drops as the cutting tool approaches the location of each region and (c) VM stress asymptotically converges towards a constant value. It can be observed that the drop of VM stresses occurs at a later stage in time in R2 and R3 compared to R1; this is because additional time is required for the cutting tool to reach those regions. The sudden drop of the VM stresses when the tool approaches a specific region is attributed to the thermal softening phenomenon which is more dominant for smaller grain sizes. This is due to the reduction of thermal conductivity and slower heat dissipation for finer grains, as discussed in the previous sections. It is evident that the final value of VM stress is lower than the initial one (t = 0) for d g,ave values of 1.98 and 2.58 nm but almost equal to the initial one for d g,ave = 4.49 nm. This phenomenon can be understood if we make a parallelisation of the nanocutting process with heat treatment. As the cutting tool passes on top of each one of the region under consideration, a temperature gradient is applied to the workpiece atoms in the vicinity of the cutting tool; this leads to the relaxation of the pre-existing residual stresses (which emerged following the quenching and equilibration processes). However, this phenomenon is not so pronounced for large grain sizes for 2 main reasons: (a) grain boundaries lie away from the cutting region and only a small part of them can be heat-treated and relieved and (b) the temperature at the cutting region is lower compared to finer grain structures (see Fig. 11and  Fig. 12). This phenomenon is also justified by the reduction of dislocation segments after nanocutting for low grain sizes as illustrated in    16. According to the authors' best knowledge this is the first time that this phenomenon has been captured via MD simulations.

Conclusions
In reality metallic materials are polycrystalline; this is a fact that has been often disregarded in numerical modelling of nanocutting processes. In this investigation a novel methodology for modelling nanocutting in polycrystalline Aluminium using Molecular Dynamics simulations has been presented. In order to generate the polycrystalline workpiece structure a FCC block of Aluminium atoms was melted, quenched and subsequently equilibrated. During the quenching process various cooling rates were implemented in order obtain various grain sizes ranging from 4.49 nm and 1.98 nm. The obtained topology was used as in nanocutting simulations in order to investigate the effect of the average grain size on the process characteristics, namely (a) cutting forces, (b) friction factor, (c) subsurface temperature, (d) defects and (e) residual stresses. The obtained results suggest that almost every one of the aforementioned characteristics are significantly affected by the average grain size. The main conclusions drawn from this study are briefly summarised below: • Both tangential (F x ) and normal (F y ) cutting forces drop for lower grain sizes within the examined range. This is attributed to grain boundary sliding.
• The friction coefficient decreases with grain size. This is again due to grain boundary sliding as well as thermal softening of the workpiece.
• The subsurface temperature is dependent on the grain size. More specifically, small grain lead to reduced phonon mean free path which in turn lowers thermal conductivity. This leads to slower heat dissipation and higher temperature in the cutting region for fine grain structures. This behaviour has been captured for the first time using MD simulations.
• The high residual stresses at the grain boundaries can be relieved for finer grains lying in the vicinity of the cutting region. This phenomenon resembles heat treatment and has been effectively captured via means of MD simulations for the first time according to the authors' best knowledge.

Declaration of Competing Interest
The authors declared that there is no conflict of interest.