## Abstract

Optimization of the wire sawing process of silicon ingots requires a profound understanding of the dynamic interaction of wire, slurry and silicon material. In this paper, the influence of wire velocity and applied wire stress on the process is investigated using dissipative particle dynamics and discrete element simulations for modelling the fluid and the grains in the abrasive suspension. In our simulations, different contact regimes occur depending on grain shape and a stress balance within the system. We observed semi-contact for high wire stress and low wire velocity and non-contact for low stress and high velocities in agreement with predictions from elasto-hydrodynamic modelling. Our simulations suggest the usage of sharp grains, since in this case, stress localization on the base of the sawing groove occurs even in the non-contact regime. These insights are likely to provide a scientific base for the optimization of sawing rates and reduction of kerf loss.

## 1. Introduction

Wire sawing is an efficient process for wafer slicing from a silicon ingot. It is widely employed both in the solar and in the semiconductor industries. Multi-wire saws use an assembly of several hundred parallel running wires to simultaneously separate the wafers. The material removed from the sawing groove gives non-recyclable Si dust and is referred to as kerf loss. A wafer thickness of about 220 μm and a kerf loss of about 180 μm can be achieved by using a typical wire diameter of 140 μm. Yet, there is still an economic demand for a decrease in wafer thickness and kerf loss while meeting high standards for wafer quality, e.g. high strength and low roughness. A systematic optimization under these conditions requires a profound and detailed understanding of the wire sawing process.

The abrasive slurry used in the sawing process typically consists of polyethylene glycol (PEG) with suspended silicon carbide (SiC) grains. The sawing wire drags the slurry through a groove in the Si ingot, where the SiC grains cause the abrasion. Experimental studies with industrial wire saws by Möller [1] showed increasing sawing rates if the applied stress on the wire or the wire velocity is increased. Sawing defects in multi-crystalline Si wafers were attributed to embedded SiC particles within the ingot material by Du *et al.* [2]. However, microscopic details of the sawing process are not fully understood, partly because experimental *in situ* investigations are difficult to conduct. This strongly motivates numerical analyses of the process. Few simulation studies have been carried out so far. Bhagavat *et al.* [3] and Möller [4] developed elasto-hydrodynamic models that predict contact regimes and pressure distributions depending on wire stress and velocity. Bhagavat & Kao [5] proposed a model to determine the optimal direction of wire approach with respect to the crystallographic orientation on the Si ingot.

In this study, numerical investigations of the dynamics of wire and slurry in the sawing groove are carried out. The influence of the applied wire stress and the wire velocity on wire kinematics and slurry dynamics is analysed. For this purpose, particle-based simulations are performed in order to take interactions of individual SiC grains within the hydrodynamic environment into account.

## 2. Numerical method

The system consisting of Si ingot, sawing wire and PEG/SiC slurry was modelled using a hybrid approach, i.e. the fluid was described by dissipative particle dynamics (DPD) and the motion of the SiC grains was treated by a discrete element approach. Since both methods are based on Newton’s equations for particles, fluid and rigid body dynamics can be easily coupled in the simulation.

All components of the system were composed of spherical particles as basic building blocks. The particle movement is governed by Newton’s equations of motion,
2.1where *m*_{i}, **r**_{i} and **v**_{i} are the mass, position and velocity of particle *i* and **F**_{ij} is the total force exerted from particle *j* on particle *i*. A velocity Verlet scheme [6] with time step Δ*t* is used for explicit time integration. The actual force laws are constructed to yield physically reasonable behaviour of the components, e.g. hydrodynamic flow of the PEG or repulsion between distinct SiC grains. Details of the modelling approach are given in the following sections.

### (a) Modelling of polyethylene glycol

An extension of the DPD method [7] is used to describe the PEG fluid. Basic DPD particles can be interpreted as single molecules or lumps of fluid molecules depending on the applied level of coarse graining. Their interaction is given by a density-dependent force [8]. The total pairwise forces are decomposed into conservative forces
2.2dissipative forces
2.3and stochastic forces
2.4where *r*_{ij}=|**r**_{ij}|=|**r**_{i}−**r**_{j}| and . *K* is a compressibility constant, *n* the average fluid number density, *n*_{i} the local number density at **r**_{i}, *γ*_{D} a dissipative constant and *σ*_{R} a random force constant. *w*(*r*_{ij}) is a monotonically decreasing weight function that vanishes for particle distances larger than the interaction range *r*_{c},
2.5*ζ*_{ij} is a Gaussian white noise with vanishing mean value and a variance of unity. Equation (2.2) models the fluid’s thermodynamics and allows for the formation of free surfaces. Equation (2.3) describes viscous damping of the fluid, and equation (2.4) is motivated by the influence of suppressed atomic degrees of freedom. A fluctuation-dissipation theorem [9] is satisfied via *σ*^{2}_{R}=2*γ*_{D}*k*_{B}*T*, where *k*_{B} is the Boltzmann constant and *T* is the temperature.

### (b) Modelling of silicon carbide grains

Each SiC grain consists of an ensemble of rigidly connected particles [10]. For this purpose, a rigid body motion solver [11] is used. The basic particles of the SiC grains interact with the DPD particles of the PEG via equations (2.2)–(2.4). Thereby, hydrodynamic coupling between the PEG model and the grains is ensured. Linear and angular momenta of the fluid are transferred to the grains causing drag and rotation. Collisions between SiC grains are modelled using a discrete element approach. Basic particles of two distinct SiC grains interact via Hertzian repulsion [12],
2.6to satisfy the excluded volume of each particle; *κ* controls the particle stiffness. An overview of the simulation parameters is given in table 1.

A microscopic image of typical SiC F400 grains is displayed in figure 1. Figure 2 shows the corresponding model grains. The relative grain size distributions used for the simulations resemble the experimentally measured one. Two types of model grains are compared, namely spheres and tetrahedrons. Although no explicit frictional force is used, torque is transmitted between grains owing to their morphological ability to interlock.

### (c) Modelling of wire and ingot

Both the sawing wire and the Si ingot are represented by rigidly connected particles. The elastic interaction between the particles of the wire or ingot and those of the SiC grains is defined by equation (2.6). At the interfaces between the wire or ingot, respectively, and the PEG fluid no-slip boundary conditions are ensured by using an adhesive DPD wall model [13]. No fracture mechanism is included in the model. Yet, stresses exerted on the Si ingot owing to the SiC grains can be extracted from the simulations. Figure 3*a* shows a schematic drawing of the wire and the ingot. The set-up is periodic along the *y*-axis thereby resembling an infinitely long wire and cutting groove.

## 3. Model calibration

DPD simulations in intrinsic units have been carried out. However, it is possible to compare numerical and experimental results by using two non-dimensional hydrodynamic numbers which characterize the system.

The Reynolds number,
3.1is the ratio of inertial to viscous forces. In a typical wire sawing process, a wire velocity *v*_{y}=15 m s^{−1} is applied. The SiC grain size distribution (compare figures 1 and 2) yields an approximate distance *d*=20 μm between the wire and the ingot. The slurry with 25 vol% SiC grains has a density *ρ*=1.7 g cm^{−3} and a viscosity *η*=200 mPa s. Thus, the Reynolds number of the system is *Re*=2.6. The estimation for the distance between the wire and the ingot in the simulations is *d*=4 *r*_{c}. The model slurry with 25 vol% SiC has a density *ρ*=8.4. Lees–Edwards boundary conditions [6] are used to measure the viscosity *η*=7. The requirement to meet the order of magnitude of the experimental Reynolds number then defines the wire velocity in the simulations. Two values are used in this study: *v*_{y}=0.16 yields *Re*=0.8 and *v*_{y}=0.64 yields *Re*=3.1.

The ratio of stress *p*_{z}, which is applied vertically on the wire, to the kinetic energy per volume (i.e. the dynamic pressure) of the slurry,
3.2renders the conversion of the applied stress between simulation and experiment possible. The stress *p*_{z} is defined by the quotient of the force in the *z*-direction on the wire and the cross-sectional area of the wire in the *x*–*y*-plane. A typical experimental value is *p*_{z}=1 MPa which yields *C*=5.2. In the simulations, *p*_{z}=1.7 and *p*_{z}=17 are used which yield in combination with the applied velocities values for *C* in the range of 1–160.

## 4. Wire sawing simulations

### (a) Set-up

Figure 3*c* shows a frontal snapshot of the sawing simulations. The wire diameter is 20 *r*_{c}, and the groove width is chosen to be 28 *r*_{c} according to a typical value of the kerf loss for the used grain size distribution. The slurry is located initially in the cutting groove while the wire is positioned above the slurry surface. Then, the dynamic simulation starts and the wire is forced into the slurry owing to an applied stress *p*_{z} while it is dragged with constant velocity *v*_{y}. The simulations are stopped after a stationary wire floating height is reached. Stress distributions on the ingot are sampled in the stationary regime.

### (b) Contact regimes

Möller [4] distinguishes for the steady-state sawing process between a non-contact and a semi-contact regime on the basis of elasto-hydrodynamic modelling depending on the wire stress and velocity. This classification is based on the distance between the wire and the ingot surface with respect to the SiC grain size. For high stresses and low velocities a so-called semi-contact regime is predicted and for low stresses and high velocites a non-contact regime is predicted. Both regimes are observed in the present simulations. When applying a low stress (*p*_{z}=1.7) and a high wire velocity (*v*_{y}=0.64), several layers of grains are located between the wire and the ingot, i.e. the system is in the non-contact regime (figure 4*a*,*c*). For a high wire stress (*p*_{z}=17) at a low velocity (*v*_{y}=0.16), semi-contact is observed (figure 4*b*,*d*): grains within a single layer are in contact with both the wire and the ingot. At low stress and low velocity, an intermediate regime exists where only the biggest grains are in contact with the wire and the ingot. High stress and high velocity yield full contact between the wire and the ingot, as all grains are pushed aside. Measured distances between wire and ingot surface in the *z*-direction for all simulations are listed in tables 2 and 3. The non-contact regime is observed for a stress ratio value of *C*=1, while semi-contact or even full contact is observed for *C* larger than 10.

### (c) Stress on ingot surface

The stress exerted from the SiC grains on the surface of the ingot groove is tracked during the simulations. The stress is averaged over time and along the *y*-coordinate because the system is translationally invariant in that direction. Figure 5 shows distributions of the stress normal to the ingot surface along the tangential coordinate of the groove in the non-contact and semi-contact regimes. The stress distribution is strongly localized around the centre of the groove in the case of semi-contact for both spherical and tetrahedral grains. In the non-contact regime, the stress is evenly distributed along the surface for spherical grains while it shows localization around the groove centre for tetrahedral grains. Note that the distributions are normalized by the applied wire stress *p*_{z}.

## 5. Conclusion

The non-contact and the semi-contact wire sawing regimes could be observed in the particle-based simulations. Furthermore, the occurrence of semi-contact (non-contact) for high wire stress and low wire velocity (low stress and high velocity) is in agreement with predictions from elasto-hydrodynamic modelling [4]. The ratio between applied wire stress and dynamic pressure of the slurry is identified as an indicator for the contact regime. For a ratio of unity only non-contact was observed, while ratios of larger than 10 yielded semi-contact or full contact. This finding provides a direct link to experimental investigations, because the stress ratio can be easily calculated based on process parameters.

It is found that in the semi-contact regime a multiple of the applied normal stress on the wire is transmitted directly to the base of the cutting groove. In the non-contact case for spherical grains, the stress exerted on the ingot is evenly distributed along the base and the sides of the groove and is equal to the applied stress in magnitude. In the non-contact case for tetrahedral grains, the stress distribution shows slight localization at the base of the groove. We assume that the abrasion is rather localized if the stress distribution is localized, which would imply reduced kerf loss. Running the process in the semi-contact regime and usage of sharp grains might therefore be beneficial.

In the future, the modelling approach should be extended by including phase transformations [14] and fracture mechanisms [15] of the silicon ingot. Further investigations should focus on the influence of the grain size distribution in order to tailor abrasive slurries for minimal kerf loss and high sawing efficiency.

## Acknowledgements

The German Ministry for the Environment, Nature Conservation and Nuclear Safety (BMU) is gratefully acknowledged for financially supporting this work within the project Kerf loss (ref. no. 0327601E). The Deutsche Solar AG, PV Silicon AG and Ersol Wafers ASi Industries GmbH are also gratefully acknowledged for financial support of this work.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society