Simulation of the ductile machining mode of silicon

Diamond wire sawing has been developed to reduce the cutting loss when cutting silicon wafers from ingots. The surface of silicon solar cells must be flawless in order to achieve the highest possible efficiency. However, the surface is damaged during sawing. The extent of the damage depends primarily on the material removal mode. Under certain conditions, the generally brittle material can be machined in ductile mode, whereby considerably fewer cracks occur in the surface than with brittle material removal. In the presented paper, a numerical model is developed in order to support the optimisation of the machining process regarding the transition between ductile and brittle material removal. The simulations are performed with an GPU-accelerated in-house developed code using mesh-free methods which easily handle large deformations while classic methods like FEM would require intensive remeshing. The Johnson-Cook flow stress model is implemented and used to evaluate the applicability of a model for ductile material behaviour in the transition zone between ductile and brittle removal mode. The simulation results are compared with results obtained from single grain scratch experiments using a real, non-idealised grain geometry as present in the diamond wire sawing process.


Introduction
Hard and brittle materials, such as silicon, are difficult to machine as they exhibit high hardness and withstand high temperatures but show low resistance in shear and tension. In cutting operations, the material removal can change from ductile to brittle when increasing the depth of cut. Brittle cutting allows for higher material removal rates but leads to surface damage which is unwanted in the manufacturing process of silicon wafers and has to be removed by subsequent etching processes. Therefore, it is desired to generate the wafer surface in the ductile Inspire AG, Zürich, Switzerland mode only. Figure 1 shows resulting surface morphologies when cutting with different cut depths where at a low cut depth, no surface damage occurs but at higher cut depths, larger surface damage and breakouts are visible. Silicon wafers are nowadays most efficiently produced by means of diamond wire sawing (DWS) [1], where geometrically undefined diamond grains fixed on a thin steel wire are indented into and moved across the workpiece material. Models used to describe the material removal process are typically built upon idealised geometries, leading to limitations when transferring results to the real process. Numerical simulations can help not only to understand phenomena associated with material removal processes but also allow to optimise the cutting process parameters such that a high material removal rate with acceptable surface properties can be achieved. Key problems are the modelling of the material behaviour as well as the numerical method to resolve the continuum. This paper presents a simulation of silicon cutting with a single diamond grain. Its geometry was determined before cutting using optical microscopy. The resulting grain geometry was 3D meshed and used to drive a process simulation with numerical methods. Here, meshfree methods, in particular the smoothed particle hydrodynamics Fig. 1 Diamond scribing of monocrystalline silicon: surface morphologies resulting from different depth of cut h cu : ductile removal h cu = 0.123μm (left), transition from ductile to brittle removal h cu = 0.722μm (middle), and brittle removal h cu = 1.225μm (right) [2] (SPH), were applied as they offer advantages in terms of computation times compared to classical mesh-based approaches: Upon large deformations, the continuum needs to be remeshed regularly, while meshfree methods automatically adapt to changes in the continuum. The goal of the study is to evaluate the applicability of the Johnson-Cook flow stress model to the cutting process of silicon in the transition zone between ductile and brittle material removal mode. Theoretical aspects of the SPH and cutting of silicon are introduced hereafter.

SPH
SPH was introduced 1977 in astrophysics for the calculation of a smoothed density from point clouds [3]. A simple derivation of the method is based on the partition of unity [4], where a field value at a spatial location x can be determined as: The Dirac-delta function δ(x) in Eq. 1 has two important properties: and Replacing the Dirac-delta function δ with a smoothing function, the so-called Kernel, W (x − x, h), e.g. the Gaußfunction, with h being a smoothing length, the behaviour of the Dirac-delta can be reproduced for the limit: Inserting (4)- (3) gives an approximation of the function value f(x) at x: This can be integrated within a discrete neighborhood using a Riemann-sum: with the point index i at which the function value is to be approximated by its neighbor points j , x ij is the spatial distance between point i and j and V j being an integration weight of point j . Computation of the function's derivative leads to where only the derivative of the Kernel W (x ij , h) is required. In this way, derivatives of values given at point cloud locations can be computed without the requirement of a functional description or a mesh-based relation between these points (particles). With this meshfree approximation derivatives in continuum mechanics equation can be computed by sums of discrete values in the particles neighborhood. Meshfree techniques were adopted in early 1990s to structural simulations [5] and for numerical cutting simulations first by [6]. At IWF of ETH Zürich the software mfree iwf was developed in the past years for SPH based cutting simulations [7]. The software is capable of performing CPU as well as GPU-enhanced computations of metal cutting simulations. It is establishing the stateof-the-art in SPH simulations as it facilitates the most recent correctors and stabilisation measures for mechanical [8] and thermal simulations [9]. The software successfully demonstrated computationally highly efficient metal cutting simulations [10], identification of friction parameters [11], and its ability to adaptive refinements [12]. Furthermore, single grain grinding simulations of engineered grinding tools [13] were successfully demonstrated.
In the presented investigation, the software is adopted to simulate the cutting of silicon with a single diamond grain assuming isotropic behaviour of silicon with the Young's modulus of Silicon in the 100 directions [14]. A flow stress model according to Johnson and Cook [15] came to application. The model is commonly used to describe metal plasticity within machining simulations and is given as: (8) with A, B, C, m, and n being material parameters, ε pl the current plastic strain,ε pl the current plastic strain rate, and T the current temperature. T f is the melting temperature, T ref is the reference temperature, andε 0 pl the reference plastic strain rate. The first two terms describe hardening due to plastic strain and plastic strain rate, respectively. The third term controls thermal softening upon increasing temperature. In this investigation, material parameters for silicon according to [16,17] were used.

Machining of silicon
Important contributions to the understanding of cutting of hard and brittle materials and the damage associated were first introduced by Hamilton and Goodman [18] and Lawn [19] who presented analytical models of the pressure field under a sliding indenter. The mechanics and mechanisms of plastic deformation and fracture of brittle materials upon indentation have been studied by Lawn et al. They introduced the median-radial crack system [20], as well as the lateral crack system [21]. Bifano et al. introduced the possibility of machining brittle materials in ductile cutting regime [22] and presented a model describing the critical depth of cut h cu,crit which presents a threshold above which brittle fracture is introduced. Based on these concepts, the ductile machining of silicon has been extensively studied taking into account the effects of various aspects. Larger cutting edge radii lead to higher normal and cutting forces [23], increase the critical depth of cut [24], and lead to horseshoe-like cracks in the grooves rather than radialchevron cracks [25]. The critical depth of cut is very sensitive to the crystallographic plane and direction and varies between 0.112μm in (001)[010] and 1.270μm in (110)[111] [26]. High cutting speed, as present in DWS of silicon, leads to amorphisation of silicon which is accompanied by a volumetric expansion and consequently lower residual cutting depth [27]. Furthermore, phase transitions happen due to high contact pressure [28][29][30], which affect the mechanical properties of the material and therefore the material removal behaviour. Recently, modelling the elastic stress field has gained attention again to predict damage [31][32][33] with good success in predicting the direction of cracks, their frequency, and their length.
With exception of Bifano's work, all studies used idealised grain shapes to study the respective phenomena, presenting a limitation to their applicability towards grinding processes. In order to exploit the freedom of geometry presented by the SPH method, the experiments conducted in this study use a non-idealised grain shape.

Experimental methodology and results
The cutting tests were performed with one single grain from a diamond wire used for the diamond wire sawing process of silicon wafers. The tests were conducted for different depths of cut h cu and the process forces F N and F C as well as the resulting surface properties were recorded.
A summary is presented in the following sub-chapters, along with a discussion of experimental results.

Optical measurement of the grain geometry
The grain chosen for the scratch test is shown in Fig. 2. The grain is embedded on a commercial diamond wire; the metallic matrix covering it has been carefully removed with sand paper. The wire is bent over a pin in such a way that the selected grain is isolated on the highest point.
Prior to the scratch test, this diamond grain was optically measured with an Alicona microscope and the surface data was retrieved as a point cloud. This point cloud was triangulated in MATLAB and a 3D mesh with tetrahedrons was created from the top part of the diamond as the cutting takes place in this region only. The surface mesh is shown in Fig. 3, the geometry part above the black plane indicates the relevant part for the simulations conducted.
Due to the geometric properties of fractured diamond grains used on electroplated diamond wires, only some of the parameters typically used to describe the tool geometry in micro machining can be applied: The projected crosssection of the grain is approximately trapezoidal with a slight tilt of the flat as shown in Fig. 2. The radii at the grain flanks are approximately 3μm on both sides, while the left side in direction of cut (right side in Fig. 2) is only engaged at large depth of cut of approximately 3μm. The cutting edge radius varies between 1.4 and 6.7 μm over the cutting edge and is smallest at the sides and larger in the centre of the grain. Rounding of the edge due to wear leads to an increase with increasing scratch length. The width-to-height ratio r of the grain is evaluated for each simulated depth of cut h cu separately: r h cu =0.2μm ≈ 27, r h cu =0.5μm ≈ 17, r h cu =1μm ≈ 13, and r h cu =1.5μm ≈ 11.

Cutting Experiments
The cutting tests were performed on a Fehlmann Picomax Versa 825 5-axis milling machine. The set-up is shown in Fig. 4.
A polished silicon surface is attached to and rotated with the spindle with angular velocity ω. The diamond grain on the wire is isolated by bending it over a pin and mounted onto a three-component force dynamometer (Kistler 9256C, eigenfrequency f ≈ 6kH z) fixed to the stationary machine table. The force measurement set-up is completed by a Kistler Type 5080A charge amplifier with active drift compensation and 1000-Hz low-pass filter and a National Instruments NI 9222 A/D converter, sampling each channel at 211.64kS/s. A pyrometer fibre was furthermore installed and connected to a Fire-3 two-colour fibre optic pyrometer for the measurement of flash temperatures. The diamond grain is kept stationary while the rotating spindle is fed vertically into the grain with normal feed speed v f and horizontally across the grain with lateral feed speed v p in such a way that the vertical displacement is equal to 3μm and the horizontal displacement is equal to 1000μm over 20 spindle rotations. The angular velocity ω is chosen in such a way that the resulting cutting speed v c equals 10m/s. The motion results in spiral-shaped contact path; due to a small but relevant perpendicularity error of the Si surface with The topography of the scratches is optically evaluated at their centre, where the engaged depth is largest. The residual surface aligned with the resulting topography profile line is shown in the bottom of Fig. 5. The residual profile is generated by plotting the median profile height for all data points that lie on a line in cutting direction at a given radial distance x. The median value has proven to be closer to the profile depth of a manually chosen representative section of each scratch than the arithmetic mean.
The residual scratch depth h res for each pass is identified as the vertical distance from a reference height determined in the unscratched area of the workpiece, shown as a dashed line "assumed surface" in Fig. 5. Cutting force F c and normal force F n are identified from the peaks in the respective signals and also aligned with the surface topography. All experimentally determined values are compiled in Table 1.

Discussion of experimental results
The determination whether the removal mode is mainly ductile, brittle, or in the transition is subjective and based on the optical appearance of the scratch. Appearance of cracks and pittings in the bottom of scratch is used as a transition criterion; however, also scratches classified as ductile show some fracture on the side.
Cracks on the side of the scratch appear first on the side with the larger engagement (left side in the bottom of Fig. 2). Their appearance is furthermore dominant in the regions of the grain where the cutting edge radius is small in agreement with [24].
It is noteworthy that scratches with large proportions of material removed in ductile cutting mode are observed with residual depths much larger than the critical depth of cut proposed by Bifano et al. and determined to be in the vicinity of 60nm for silicon [34], that is, in spite of elastic recovery which happens after the cut and renders the engaged depth deeper than the residual depth. The observation is, however, in agreement with findings published in [26], where the critical depth of cut for ductile removal of silicon with a cone shaped indenter with a cutting edge radius of 8μm was determined to lay between 0.122 − 1.270μm depending on the crystallographic plane and orientation. It is furthermore observed, that the scratches 16-18 show a ductile bottom of the scratch with comparably few breakouts on the edges while the scratches before and especially the scratches after them show significant brittle fracture also on the bottom of the scratch. The residual scratch depth decreases over the last scratches in spite of the increasing in-feed (evident in the increasing normal force), which might be attributed to increased amorphisation and associated volumetric expansion [27] and increased elastic recovery.

Fig. 6 Workpiece dimensions and boundary conditions
It was further attempted to quantify the elastic recovery and the undeformed chip thickness from the experimentally determined residual depths by applying models found in literature. For ductile removal, a reasonable undeformed chip thicknesses that can be used in this context could not be determined, which is mainly due to the large deviations between the idealised indenter geometries used to develop the model and the real grain geometry used in this study. The model has to be developed further which goes beyond the scope of this study. For this reason, comparison between experimental and stimulative results is based on residual scratch depth only. The models evaluated are only described briefly and can be referred to in the references mentioned.  For ductile material removal, indentation theory proposed by Oliver and Pharr [35] and a model developed and validated for single grain scratching of silicon by Huang et al. [36] is applied to consider elastic recovery, with the assumption that the grain tip is spherical. For a spherical indenter, the actively engaged tip radius can be determined from the residual scratch width and the contact depth. With a known tip radius and the assumption that the residual scratch width is equal to the width of the tool tip in engagement [35], it is possible to solve the equations for the undeformed chip thickness .
In case of brittle material removal, a model developed by Marshall et al. [21] can be applied to accommodate spalling of material due to brittle fracture. Since both the scratch width and depth are affected by brittle removal, consideration of forces and material parameters is required to determine the undeformed chip thickness. The determination requires the assumption that the depth of the material removed equals the depth of the plastic zone [21] Fig. 9 Temperature dependent thermal conductivity [39] and approximation used in this investigation Fig. 10 Temperature dependent specific heat [39] and approximation used in this investigation and that material on the side of the groove is removed due to the propagation of lateral cracks to the surface [37]. The plastic zone extends equally in half width and depth and equals the residual, measured depth of the scratch, while the length of the lateral cracks equals half of the measured width of the residual scratch. The engaged depth equals the undeformed chip thickness [21]. Details of the application of this approach to real grain geometries can be found in [38].

Numerical model
The modelled workpiece geometry has a size of 40μm by 40μm by 6μm (length, width, and height) and is discretised with 259 200 (120×120×18) particles. The workpiece is restrained at the bottom, left, and right faces and a temperature boundary condition of T = 300K is applied. Geometry and boundary conditions are depicted in Fig. 6.
The diamond grain is modelled as a rigid body with a mesh consisting of 4426 tetrahedron elements. The discretised workpiece and meshed diamond grain are shown in Fig. 7 together with the directions of the acting process forces. Heat conduction is considered in the workpiece using the Particle Strength Exchange (PSE) approximation of the Laplacian. Tool heat transfer was not considered.
Coulomb friction is used with a friction coefficient assumed to be μ = 0.3. The fraction of plastic work that is converted into heat (Taylor-Quinney coefficient) in silicon cutting is unknown and is assumed here to be η T Q = 90% for all simulations. Heat generation due to friction is considered with η W P f ric = 8.5% of the dissipated frictional energy being transferred into the workpiece (WP). The ratio stems from the assumption of a proportional split of the frictional heat transferred into tool and workpiece according to their respective thermal diffusivities α W P and α Dia . This Temperature-dependent properties for the Young's modulus (Fig. 8), the thermal conductivity (Fig. 9), and the specific heat ( Fig. 10) are used. All material parameters used in the simulations are summarised in Table 2.
The simulations are performed for cut depths h cu of 200nm, 500nm, 1000nm, and 1500nm and a cutting speed of v c = 10m/s as in the conducted cutting experiments. An overview of all simulated process parameters is provided with Table 3, along with the resulting forces. All simulations are carried out with GPU acceleration on a NVIDIA Quadro GP100 graphics card. The computation time of each scratch simulation was in the order of 5.75 h for ∼ 270 000 time steps. 1 With material parameters for diamond from [7]:

Surface topography and elastic recovery
For the generation of the residual profile and estimation of the elastic recovery, the SPH particle positions at the simulation end are triangulated with Open3D [40] and the depth coordinates are mapped from the particles to the mesh. Figure 11 shows the residual depth after the cut. The numerical simulation predicts for the lowest depths of cut h cu = 200nm and h cu = 500nm no chip or debris formation, while for cut depths of h cu = 1000nm and h cu = 1500nm, a chip develops in front of the tool as well as smaller breakouts and debris sideways of the tool are created. These breakouts are much smaller than those observed in the experiments.
The first two simulated scratches with h cu = 200nm and h cu = 500nm result in very shallow groves of h res = 10nm and h res = 100nm respectively which cannot be matched up with experimental results as the lowest experimentally determined residual depth is 550nm and therefore much larger. The 3rd and 4th simulated scratches produce groves of a residual depth of h res = 800nm and h res = 1300nm respectively, which lays in the vicinity of the experimental scratches number 3 to 6, for which an enlarged view is provided in Fig. 12.
The general appearance of the simulated scratches agrees with that of the experimentally cut ones: the left flank in direction of cut is steeper with "break-outs" which are significantly deeper than those on the right flank. The bottom of the grooves gets rougher with increasing residual depths in both cases.
For the lowest simulated depth of cut of h cu = 200nm, the imprint has a maximum depth of around 10nm which implies an elastic recovery depth h rec of 190nm. At h cu = 500nm, the elastic recovery is around 400nm and for h cu = 1000nm and h cu = 1500nm about 200nm. In comparison with experimental results obtained from grinding of silicon with grains of similar average size (≈ 25.3μm) and cutting speed of v c = 44m/s [36], the simulated recovery depth is higher for small depth of cut and lower for high depth of cut. Recovery depth of h rec,h cu =200nm ≈ 100nm, h rec,h cu =500nm ≈ 200nm, and h rec,h cu =800nm ≈ 300nm were observed in the grinding experiments.

Process forces
In comparison to the experimental values from Table 1, the predicted process forces are generally higher. The errors in the prediction become smaller towards higher depth of cut h cu . A comparison of the experimentally and numerically determined forces is provided with Fig. 13 where forces are plotted over the residual depth h res .
While the numerical cutting force is in an acceptable agreement to the experimental values, the simulated normal force tends to be higher. The reason for the over-estimation of the normal forces may lie with the very low fracture toughness of silicon of < 1 MPam 1/2 [16] which is not considered with the modelling approach here. A small flaw in the workpiece initiates crack growth and brittle fracture with low energy intake. The cracks will not lead to breakouts as long as the hydrostatic stress is compressive; however, the residual hydrostatic stresses in the wake of the tool are partially in tension (refer to chapter Section 4.4) and indicate that the breakouts may likely occur. The ductile deformation and shearing of the material which The cutting direction is from bottom to top. Note that different colour scales are used. Dark-blue-coloured regions are related to material which is still in the initial zero-plane of the surface, while material accumulation and debris extending above the initial surface are coloured in grey are simulated require significantly more energy leading to higher forces since fracture is not modelled.
The simulated cutting and normal force characteristics during the cut are shown for all simulations in Fig. 14. At increased depth of cut, the process forces are slightly varying due to formation of chips and small debris as well as torn out material.

Hydrostatic pressure
The hydrostatic pressure is investigated below the cutting tool where positive values are in compression. The pressure distributions are displayed in Fig. 15 with the tool shown transparent to visualise the hydrostatic pressure distribution below the tool as well. A video of this simulation is uploaded to https://youtu.be/nMTaoP0XIj4.
The increasing contact area with increasing depth of cut can be seen. The region of maximum hydrostatic  pressure shifts towards the right side (in cutting direction) of the grain. At the lowest simulated depth of cut h cu = 200nm, the maximum hydrostatic pressure is in the order 4GPa, at h cu = 500nm increasing to 13 GPa while with further increasing cut depth the area of highly compressed zones increased steadily and surpassed in some spots even 60GPa. The levels of hydrostatic pressures found in the numerical prediction give rise to the assumption that phase transformation in the workpiece will occur, as for example described in [28] and simulated and validated in [27]. At the lowest cut depth, very small hydrostatic pressures in compression remain in the wake of the tool, while at h cu = 500nm they are in the order of 3 − 4GP a. At both larger depths of cut of h cu = 1000nm and h cu = 1500nm, the residual hydrostatic pressures are lower in the cut surface compared to h cu = 500nm but exceed 5 GPa locally (green to red in Fig. 15). On the other hand, they turn partially into tension (dark-blue in Fig. 15). This observation is in general agreement with a SPH simulation presented in [41] where cutting of KDP crystal was modelled by means of constitutive model based on the relative length of crack. Figure 16 shows the stress distributions for the four depths of cut . The equivalent stress reaches up to 5.5 GPa under the tool at h cu = 200nm and increases locally up to 10-11 GPa for the three higher depths of cut of h cu = 500nm, h cu = 1000nm, and h cu = 1500nm. A video of the simulations is uploaded to https://youtu.be/QZw88fNM1CA.

Stress distribution
In the wake of the tool, residual stresses form which, when overlayed with tensile hydrostatic pressures, potentially lead to opening of cracks and breakouts. The cracks

Predicted temperatures
The predicted peak temperatures in the workpiece strongly depend on the cut depth h cu and increase from around T = 310K (h cu = 200nm) up to the melting temperature T = 1700K for h cu = 1000nm and h cu = 1500nm, see Table 3 and Fig. 17. The temperatures at the higher depths of cut are in line with observations for experiments conducted in the same manner as the one presented in this study [38] where the flash temperature was measured to lie in the range around 1500K; however, the experimental results show that the temperature decreases with increasing depth of cut and brittle behaviour. At lower cut depths of h cu = 200nm and h cu = 500nm, the temperatures are therefore well below the experimental values. For this reason, the accuracy of the temperature prediction is questionable. A video of the simulations is uploaded to https://youtu.be/ Vw4zPtsh3Do.

Conclusion
Single grain cutting in silicon in ductile mode was simulated with the SPH and compared with experimental findings applying the same non-idealised grain geometry. As expected, the results indicate that the literature data for the Johnson-Cook flow stress model parameters have limited validity for depths of cut that lay in and above the ductile-to-brittle transition.
The process forces are predicted higher than experimentally observed, where the cutting forces only slightly exceed the experimental values which indicates that the assumed While at the lowest simulated cutting depth compressive hydrostatic pressure reached around 4GP a, it then increased to values of 13GP a at a cutting depths of h cu = 500nm which is sufficient for phase transformations in the . The cut direction is from right to left. Note that different temperature scales are used silicon. Simulations with increased depth of cut showed strongly increased hydrostatic pressures under the cutting tool in excess of 60GP a locally. These high values at larger cutting depths have to be treated with care as the normal force predictions deviate from the experimental values and therefore are subject to changes upon improved constitutive modelling.
Predicting the transition from ductile to brittle cutting mode is not possible with the Johnson-Cook flow stress model only. Instead, at least extensions, e.g. to a fracture criterion and subsequent crack propagation analysis is required. In addition, more elaborated constitutive models which consider the anisotropic nature of silicon as well as phase transformations under high compressive stresses are expected to enhance the quality of the prediction of the forces and the critical depth of cut for the ductile-to-brittle transition.
Since the temperature predictions are not satisfactory, especially for lower cut depths, the fraction of plastic work dissipated into heat (Taylor-Quinney coefficient) and the friction coefficient between silicon and diamond need to be investigated experimentally. This, together with the modelling of tool heat transfer is expected to further improve the results.
On the experimental side of the study, it was shown that, with a real grain geometry, predominantly ductile scratches with only few breakouts on the sides and bottom of the grooves are generated at depths of cut that are much larger than typically observed in cutting experiments.
Acknowledgements The authors hereby would like to thank the SNF for partial funding of this work from the projects: -

Competing interests
The authors declare no competing interests.
Availability of data and material (data transparency) All data generated or analysed during this study are included in this published article.
Code availability (software application or custom code) The software tool mfree iwf used in this investigation is open source and can be obtained from here.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.