- Email: [email protected]
Mechanism of ring crack initiation in Hertz indentation of monocrystalline silicon analyzed by controlled molecular dynamics
CIRP Annals - Manufacturing Technology 59 (2010) 559–562
Contents lists available at ScienceDirect
CIRP Annals - Manufacturing Technology jou rnal homep age : ht t p: // ees .e lse vi er. com/ci rp/ def a ult . asp
Mechanism of ring crack initiation in Hertz indentation of monocrystalline silicon analyzed by controlled molecular dynamics T. Inamura (2)*, Y. Shishikura, N. Takezawa Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya 466-8555, Japan
A R T I C L E I N F O
A B S T R A C T
Keywords: Simulation Fracture analysis Molecular dynamics
The mechanism of ring crack initiation in the Hertz indentation of monocrystalline silicon with no preexisting defect has been analyzed by controlled molecular dynamics. It has been found that microvoids that develop into a ring crack can be generated outside the outer periphery of the contact area between the silicon and an indenter, such that the static stress and dynamic stress associated with acoustic waves locally transform the monocrystal structure to a polycrystal one, and then the static stress causes cross slips at grain boundaries to cause microvoids. ß 2010 CIRP.
1. Introduction
these studies indicate that there can be defect initiation, it has been concluded that more basic and in-depth studies on the process are necessary using more clearly defined stress fields than those of cutting. In this study, we focus on the mechanism of ring crack initiation in the well-defined stress field of Hertz indentation. The material to be indented is semi-infinite monocrystalline silicon with no preexisting defect. We use the controlled MD to observe and analyze the phenomena.
Brittle materials such as monocrystalline silicon exhibit socalled brittle/ductile transition phenomena in machining. Thus, the machining of these materials must be carried out by choosing conditions for ductile-mode machining in order to produce a highquality surface but at the cost of low productivity. Because brittlemode machining is accompanied by the brittle fracture of workpiece materials and brittle fracture has been analyzed on the basis of fracture mechanics, the brittle/ductile transition phenomena have also been studied in the framework of fracture mechanics [1]. Fracture mechanics discusses the condition under which a microcrack develops into a macrocrack assuming that the microcrack or a defect already exists. Because intrinsic defects exist in many brittle materials, fracture mechanics have been applied to brittle materials without seriously considering the above assumption. However, it is also true that carefully prepared monocrystalline silicon with almost no preexisting defects causes fracture in machining. This indicates that there is a mechanism by which defects can be initiated in the machining itself. Therefore, if the mechanism of defect initiation is clarified and combined with the theory of existing fracture mechanics, we will be able to understand the fracture phenomena more deeply and thus be able to control them more efficiently. On the basis of the above idea, the authors have carried out studies on the mechanism of defect initiation in cutting [2,3]. In these studies, new methods of studying the process based on molecular dynamics (MD) (see [4–6] for MD) have been devised because fracture occurs only in macroscale materials and defects of atomic size maybe initiated in picoseconds. Although the results of
2. Experimentally observed phenomena in Hertz indentation Facts contradicting the simple idea that fracture occurs when the tensile stress is larger than the theoretical tensile strength of a material are as follows [7]. A ring crack appears not on the periphery of the contact area at which the tensile stress is maximum but at an area outside the periphery where the tensile stress is about one-tenth of the theoretical tensile strength; the experimental data do not follow the Hertz contact theory, according to which the critical load that causes a ring crack should be proportional to the square of the radius of the indenter. The mechanisms of these phenomena have been explained by assuming preexistence of microcracks in materials. In these explanations, additional assumptions for the preexisting microcracks have usually been introduced with regard to their length and spatial distributions. Although the statistical nature of the preexisting defects accounts for the scattering of experimental data, it has been claimed that the observed data are less scattered than that predicted from the explanation [8]. 3. Analysis method
* Corresponding author at: Department of Mechanical Engineering, Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya 466-8555, Japan. 0007-8506/$ – see front matter ß 2010 CIRP. doi:10.1016/j.cirp.2010.03.139
The controlled MD [3,9] is a method that satisfies the condition described in Section 1. In this method, an MD model is further
560
T. Inamura et al. / CIRP Annals - Manufacturing Technology 59 (2010) 559–562
Fig. 1. Entire view of the indentation model and the MD model embedded in it.
coarse-grained and is regarded as being composed of a specified number of atomic clusters. Then, the motions of those clusters, which are the mean motions of their respective member atoms, are controlled so as to coincide with the displacements of an analytical solution at the center of each cluster. The important point here is that the free relative motion of each atom is still possible so that dislocations and/or slips can occur, as in an actual material. Because analytical solutions can be derived with macroscopic boundary conditions such as an infinite or a semi-infinite boundary, the MD model in these cases can be regarded as being embedded in the analytical solution fields. Then, the results of the controlled MD can be regarded as locally magnified images of the analytical solutions. The analytical displacement field of the Hertz contact theory is shown in Fig. 1 together with the MD model used in this study. Because the analytical field is axisymmetric, we take its cross-section and embed the MD model in it with a periodic boundary condition in the direction normal to the cross-section in Fig. 1. The analytical displacement field for the tip radius of an indenter (100 mm) has been derived on the basis of the method described in [10] as a function of the indentation load P, which increases with time t as P = ct3, c = 9.62 1044 kg m s 5. This change in the load P is set so that the indentation depth reaches 7 mm after 104 steps for a 10 15 s time step of the MD simulation. Note that 7 mm is sufficiently deep to cause a ring crack. We use the elastic Hertz contact theory to control the MD because microscopic plastic phenomena, such as defect initiation, occur in a macroscopically elastic field. The MD model is composed of 6 clusters in the r and z directions and 2 clusters in the direction normal to the figure, where the size of a cluster is 6a 6a 6a (a: lattice constant of silicon). This MD model is embedded on the surface at a distance of 36.8 mm from the center so that it passes through the area of ring crack initiation and reaches the outer periphery of the contact area after 104 steps. The value of Young’s modulus and Poisson’s ratio used in the analytical solution are 786 MPa and 0.2 for the diamond indenter, and 170 MPa and 0.2 for the silicon specimen, respectively. The interatomic potential used in the MD is that proposed by Tersoff [11]. 4. Results and discussion The deformation of the MD model during indentation is shown in Fig. 2 by taking a slice with a thickness 2a. It is shown in the figure that the model is inclined along the slope of the specimen surface outside the indenter and is being stretched in the r direction with a shear deformation in the r and z directions. It can be seen in Fig. 2 that the model in (a) is deformed only elastically but the model in (b) contains a crystal structure different from that of the diamond structure around the center of the model. It is also shown in (b) that amorphous areas are
Fig. 2. Deformation and phase transformation of the embedded MD model during indentation. (a) 7400 step, (b) 9000 step, (c) 8700 step, and (d) mechanism of amorphization.
formed at the right diagonal corners of the model. Then, the above results should be analyzed on the basis of the following facts. It is known theoretically as well as experimentally that the diamond structure of silicon transforms to a b-tin structure under a hydrostatic pressure of about 10 GPa but never transforms to an amorphous structure under any hydrostatic pressure [12]. The latter fact implies that the amorphous structure is more unstable than the b-tin structure because the b-tin structure has a lower average atomic potential energy than the amorphous structure. On the other hand, it is also known experimentally that an amorphous layer is formed on the silicon surface during indentation and/or cutting [13]. Thus, it is concluded that the phase transformation to the amorphous state in the indentation and/or cutting must be caused by the introduction of a shear stress. The latter fact indicates that the shear stress can cause a transformation from the diamond structure to an amorphous state with an increase in the average atomic potential energy. Then, a b-tin structure may appear as an interim structure during the phase transformation to the amorphous state. On the basis of the above facts, the results in Fig. 2(b) indicate that the new crystal structure that has appeared around the center of the model is b-tin with a cubic structure. In addition, it is found that the amorphous areas at the corners of the model in Fig. 2(b) have been created as the result of cross slips that occurred in the btin structure. This has been found by observing the deformation of
T. Inamura et al. / CIRP Annals - Manufacturing Technology 59 (2010) 559–562
Fig. 3. Dynamic change in density distribution due to elastic waves emitted, transmitting and reflected in the model. Waves interfere with each other. (a) 7800 step, (b) 8200 step, (c) 8600 step, and (d) 9000 step.
a red ribbon in Fig. 2(c). The process illustrated in Fig. 2(d) shows that a slip occurs in the right diagonal direction ([1 1 1] direction) owing to compression and shearing in the z direction, leaving the shrunk part of the red ribbon as indicated by a circle in Fig. 2(c). Then, another slip occurs in the left diagonal direction (also the [1 1 1] direction) owing to tension in the r direction, which bends both ends of the shrunk part of the red ribbon as shown in Fig. 2(c). Although the amorphous areas are formed during this process, the results in Fig. 2 show that defects such as microcracks cannot be created through the phase transformations. Note that the results obtained above are of high accuracy because the parameters in the interatomic potential used are those adjusted so as to reproduce the actual atomic potential energy of both the diamond and the btin structures. Thus, the results indicate that there is still another factor that must be introduced into our simulation to reproduce the ring crack initiation. It is reported that acoustic emission (AE) is observed in monocrystalline silicon when a phase change to the amorphous state takes place [14]. This fact can also be confirmed in Fig. 3, which shows the dynamic change in the density distribution due to the AE in the model. The result in Fig. 3 has been obtained by the following procedure: because the clusters in the model are controlled by the analytical solution such that they move in a quasi-static manner, only acoustic waves whose wavelengths are shorter than the cluster size and thus can travel without disturbing the cluster motions can exist in the model [9]. Then, on the basis of the wave speed and cluster size, we can calculate the time in which the wave with the longest possible wavelength travels its half wavelength. We then calculate the difference in the atomic density distributions between the start and the end of this time interval. By this procedure, we can extract the density distribution due to the acoustic wave with only the longest wavelength by doubly magnifying its amplitude, and we can approximately filter out the density distributions due to the other wave components because they almost cancel each other. It is shown in Fig. 3 that the density distribution due to the acoustic wave thus extracted changes dynamically, indicating complicated interference among waves emitted from several locations in the model and between the original waves and those reflected at the four boundaries of the model. Although the reflection occurs only at the top free surface in actual indentation, the complicated nature of the density distribution due to
561
interference will remain unchanged even in this case. On the other hand, the acoustic waves in actual indentation will propagate outward beyond the model and inward from outside the model, repeatedly interfering with each other and with waves reflected at the top free surface along the traveling path. We thus take these phenomena into account in our simulation. The most accurate and direct method of introducing the above effects will be to use an MD model sufficiently large to include the emission of waves as well as their transmission. However, because this would require a high computational cost and a long calculation time, we use an approximate method that uses the model in Fig. 2. Here, we first consider the fact that the outer periphery of the contact area between the indenter and the silicon approaches our model, moving along the surface of the silicon. We then focus on the effect of surface waves that travel along this path. Because the left-outside area of the model approaches the contact area earlier than the model itself and thus will cause a phase transformation prior to inside the model, waves emitted in the area will propagate, with some attenuation, into the model after a time determined by the wave speed and the distance from the wave source. Here, because the outer periphery of the contact area is moving toward our model, the area of the wave source will be the area of the model after an appropriate time. This means that the waves emitted in the model can be superimposed, with some attenuation, on the event that occurred in the same model at an appropriate time before. Thus, in this study, we consider three models with different time steps that are aligned side by side and approaching together the contact area. Then, we superimpose the waves extracted from the middle and left-end models on the events in the right-end model after multiplying the extracted waves by an attenuation factor and maintaining a time difference of 2000 and 4000 steps, respectively. Here, the time difference of 2000 steps is determined so that the waves travel across the model. The method of wave extraction is that described in Fig. 3, and the amplitudes of the waves extracted from the middle and left-end models have been multiplied by the reference values of 0.3 and 0.09, respectively, before superimposing them on the atomic motion in the right-end model. This is because an accurate value for the attenuation factor cannot be determined easily. The superposition has been carried out in the following manner: the extracted dynamic density distribution is evaluated at the center of each cluster and the values are transferred to the corresponding clusters in the rightend model by giving their member atoms the displacement required to cause the corresponding volume changes. These procedures are carried out at each step of the MD simulation during the period from 4000 to 6000 steps of the right-end model. The above method will contain errors in the sense that the transferred dynamic density change in the right-end model includes the effect of fictitious waves reflected at the boundaries in the original models and does not include that of actual waves reflected at the top free surface along the actual traveling path. The method also contains an error in the sense that the indentation depth of 4 mm, at which the left-outside area of the right-end model causes a phase transformation, is different from the indentation depth of 5 mm when the corresponding event occurs in the middle model. However, it is considered that the errors in both cases are not critical because the density distribution is in any case of random nature and because the microscopic states at the local area close to the outer periphery of the contact surface do not differ significantly between the indentation depths of 4 and 5 mm. Then, it is found that the distance to the nearest-neighbor atoms in the right-end model is distributed between 2.3 and 2.5 A˚ during the period without the effect of transmitted waves, whereas the corresponding value is spread between 2.1 and 2.6 A˚ when the effect of transmitted waves is introduced. The result of the simulation including the effect of transmitted waves is shown in Fig. 4 using the same slice as that in Fig. 2. It can be seen in Fig. 4 that only a simple shear deformation occurs in (a) but a polycrystal structure that is composed of a mixture of diamond and b-tin structures appears in (b). The reason why the
562
T. Inamura et al. / CIRP Annals - Manufacturing Technology 59 (2010) 559–562
forms a cross slip, as shown in (4). Here, because the cross slip occurs at a grain boundary, the atoms at the cross point cannot move freely, and thus a sufficiently high tensile stress to create a void appears locally. Then, the model collapses because the microvoids create a stress field of the typical mode-I crack front and trigger the formation of a ring crack. The reason why a ring crack occurs not at the outer periphery of the contact area but at an area outside of it is because the condition mentioned above is satisfied there first. It can also be concluded that the positions of ring cracks thus created will scatter even in the case of no preexisting defect because the dynamic stress distribution has a random nature to a certain extent owing to the complex interference between emitted and reflected waves. Finally, we point out that the time from the formation of a local polycrystal structure to the initiation of microvoids is only about 2 ps, and thus the process can be studied only by the method used in this study at present. However, we should endeavor to observe and confirm the process experimentally in the future. 5. Conclusions (1) Microvoids that trigger the formation of a ring crack can be initiated even in monocrystalline silicon with no preexisting defect. The void initiation occurs through the combined effect of dynamic stress and static stress. (2) The dynamic stress associated with acoustic waves originates from an amorphous area that is created by the static stress around the outer periphery of the contact area between an indenter and silicon. (3) The dynamic stress superimposed on the static stress field of Hertz theory locally transforms a monocrystal structure into a polycrystal structure outside the amorphous area, and then the static stress causes cross slips at the resulting grain boundaries in the area to cause microvoids. References
Fig. 4. Phase transformation to polycrystal structure followed by void initiation caused by cross slips. (a) 4300 step, (b) 5200 step, (c) 5300 step, (d) 5800 step, (e) 5290 step, and (f) mechanism of void creation.
b-tin structure appears locally and earlier than in Fig. 2 can be attributed to the fact that the dynamic compressive/tensile stresses associated with acoustic waves act in addition to the static stress, such that the combined shear stress causes a local phase transformation. Then, as shown in Fig. 4(c), a microvoid is created around a grain boundary. Other microvoids then appear at several grain boundaries, and these microvoids gradually coalesce leaving amorphous areas around them, as shown in Fig. 4(d). The ribbon mark in Fig. 4(e) indicates how the silicon has been deformed around the void in Fig. 4(c). It is shown in the figure that the ribbon undulates violently and is almost torn off, indicating that a cross slip has occurred there. Then, the deformation of the ribbon is observed during the void creation, and the mechanism of which is shown schematically in Fig. 4(f). In Fig. 4(f), a polycrystal structure is formed in (1), as indicated by different colors, and a local slip is about to occur in the direction of 45 degrees owing to shear and compression in the z direction. In (2), the local slip described above has occurred along a grain boundary. Then in (3), another slip is about to occur in the direction of 45 degrees owing to the tensile stress in the r direction. The slip described in (3)
[1] Ueda K, Sugita T, Hiraga H (1991) A J-integral Approach to Material Removal Mechanisms in Microcutting of Ceramics. Annals of the CIRP 40(1):61–64. [2] Inamura T, Shimada S, Takezawa N, Ikawa N (1999) Crack Initiation in Machining Monocrystalline Silicon. Annals of the CIRP 48(1):81–84. [3] Inamura T, Shishikura Y, Takezawa N (2009) Digital Microscope Observation of the Initial Stage of Cutting Monocrystalline Silicon. Annals of the CIRP 58(1):69–72. [4] Rentsch R, Inasaki I (1994) Molecular Dynamics Simulation for Abrasive Process. Annals of the CIRP 43(1):327–330. [5] Komanduri R, Chandrasekaran N, Raff LM (1999) Orientation Effects in Nanometric Cutting of Single Crystal Materials: An MD Approach. Annals of the CIRP 48(1):67–72. [6] Shimada S, Ikawa N, Tanaka H, Ohmori G, Uchikoshi J, Yoshinaga H (1993) Feasibility Study on Ultimate Accuracy in Microcutting Using Molecular Dynamics Simulation. Annals of the CIRP 42(1):91–94. [7] Wilks E, Wilks J (1991) Properties and Application of Diamond. ButterworthHeinemann. 192–211. [8] Harrison J, Wilks J (1978) The Hertz Indentation Test and Auerbach’s Law. Journal of Physics D 11(1):73–81. [9] Inamura T, Takezawa N, Shibuya K, Yamada K (2007) Dynamic Phenomena at Mode-I Crack Front in Silicon Simulated by Extended Molecular Dynamics. Annals of the CIRP 56(1):561–564. [10] Timoshenko SP, Goodier JN (1970) Theory of Elasticity. third edition. McGrawHill. 398–414. [11] Tersoff J (1989) Modeling Solid-State Chemistry: Interatomic Potentials for Multicomponent Systems. Physical Review B 39(8):5566–5568. [12] Hu JZ, Merkle LD, Menou CS, Spain IL (1986) Crystal Data for High-Pressure Phases of Silicon. Physical Review B 34(7):4679–4684. [13] Clarke DR, Kroll MC, Kirchner PD, Cook RF (1988) Amorphization and Conductivity of Silicon and Germanium Induced by Indentation. Physical Review Letters 60(21):2156–2159. [14] Shimada S, Inamura T, Takezawa N, Ohmori H, Sata T, Ikawa N (1995) BrittleDuctile Transition Phenomena in Microindentation and Micromachining. Annals of the CIRP 44(1):523–526.
Contents lists available at ScienceDirect
CIRP Annals - Manufacturing Technology jou rnal homep age : ht t p: // ees .e lse vi er. com/ci rp/ def a ult . asp
Mechanism of ring crack initiation in Hertz indentation of monocrystalline silicon analyzed by controlled molecular dynamics T. Inamura (2)*, Y. Shishikura, N. Takezawa Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya 466-8555, Japan
A R T I C L E I N F O
A B S T R A C T
Keywords: Simulation Fracture analysis Molecular dynamics
The mechanism of ring crack initiation in the Hertz indentation of monocrystalline silicon with no preexisting defect has been analyzed by controlled molecular dynamics. It has been found that microvoids that develop into a ring crack can be generated outside the outer periphery of the contact area between the silicon and an indenter, such that the static stress and dynamic stress associated with acoustic waves locally transform the monocrystal structure to a polycrystal one, and then the static stress causes cross slips at grain boundaries to cause microvoids. ß 2010 CIRP.
1. Introduction
these studies indicate that there can be defect initiation, it has been concluded that more basic and in-depth studies on the process are necessary using more clearly defined stress fields than those of cutting. In this study, we focus on the mechanism of ring crack initiation in the well-defined stress field of Hertz indentation. The material to be indented is semi-infinite monocrystalline silicon with no preexisting defect. We use the controlled MD to observe and analyze the phenomena.
Brittle materials such as monocrystalline silicon exhibit socalled brittle/ductile transition phenomena in machining. Thus, the machining of these materials must be carried out by choosing conditions for ductile-mode machining in order to produce a highquality surface but at the cost of low productivity. Because brittlemode machining is accompanied by the brittle fracture of workpiece materials and brittle fracture has been analyzed on the basis of fracture mechanics, the brittle/ductile transition phenomena have also been studied in the framework of fracture mechanics [1]. Fracture mechanics discusses the condition under which a microcrack develops into a macrocrack assuming that the microcrack or a defect already exists. Because intrinsic defects exist in many brittle materials, fracture mechanics have been applied to brittle materials without seriously considering the above assumption. However, it is also true that carefully prepared monocrystalline silicon with almost no preexisting defects causes fracture in machining. This indicates that there is a mechanism by which defects can be initiated in the machining itself. Therefore, if the mechanism of defect initiation is clarified and combined with the theory of existing fracture mechanics, we will be able to understand the fracture phenomena more deeply and thus be able to control them more efficiently. On the basis of the above idea, the authors have carried out studies on the mechanism of defect initiation in cutting [2,3]. In these studies, new methods of studying the process based on molecular dynamics (MD) (see [4–6] for MD) have been devised because fracture occurs only in macroscale materials and defects of atomic size maybe initiated in picoseconds. Although the results of
2. Experimentally observed phenomena in Hertz indentation Facts contradicting the simple idea that fracture occurs when the tensile stress is larger than the theoretical tensile strength of a material are as follows [7]. A ring crack appears not on the periphery of the contact area at which the tensile stress is maximum but at an area outside the periphery where the tensile stress is about one-tenth of the theoretical tensile strength; the experimental data do not follow the Hertz contact theory, according to which the critical load that causes a ring crack should be proportional to the square of the radius of the indenter. The mechanisms of these phenomena have been explained by assuming preexistence of microcracks in materials. In these explanations, additional assumptions for the preexisting microcracks have usually been introduced with regard to their length and spatial distributions. Although the statistical nature of the preexisting defects accounts for the scattering of experimental data, it has been claimed that the observed data are less scattered than that predicted from the explanation [8]. 3. Analysis method
* Corresponding author at: Department of Mechanical Engineering, Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya 466-8555, Japan. 0007-8506/$ – see front matter ß 2010 CIRP. doi:10.1016/j.cirp.2010.03.139
The controlled MD [3,9] is a method that satisfies the condition described in Section 1. In this method, an MD model is further
560
T. Inamura et al. / CIRP Annals - Manufacturing Technology 59 (2010) 559–562
Fig. 1. Entire view of the indentation model and the MD model embedded in it.
coarse-grained and is regarded as being composed of a specified number of atomic clusters. Then, the motions of those clusters, which are the mean motions of their respective member atoms, are controlled so as to coincide with the displacements of an analytical solution at the center of each cluster. The important point here is that the free relative motion of each atom is still possible so that dislocations and/or slips can occur, as in an actual material. Because analytical solutions can be derived with macroscopic boundary conditions such as an infinite or a semi-infinite boundary, the MD model in these cases can be regarded as being embedded in the analytical solution fields. Then, the results of the controlled MD can be regarded as locally magnified images of the analytical solutions. The analytical displacement field of the Hertz contact theory is shown in Fig. 1 together with the MD model used in this study. Because the analytical field is axisymmetric, we take its cross-section and embed the MD model in it with a periodic boundary condition in the direction normal to the cross-section in Fig. 1. The analytical displacement field for the tip radius of an indenter (100 mm) has been derived on the basis of the method described in [10] as a function of the indentation load P, which increases with time t as P = ct3, c = 9.62 1044 kg m s 5. This change in the load P is set so that the indentation depth reaches 7 mm after 104 steps for a 10 15 s time step of the MD simulation. Note that 7 mm is sufficiently deep to cause a ring crack. We use the elastic Hertz contact theory to control the MD because microscopic plastic phenomena, such as defect initiation, occur in a macroscopically elastic field. The MD model is composed of 6 clusters in the r and z directions and 2 clusters in the direction normal to the figure, where the size of a cluster is 6a 6a 6a (a: lattice constant of silicon). This MD model is embedded on the surface at a distance of 36.8 mm from the center so that it passes through the area of ring crack initiation and reaches the outer periphery of the contact area after 104 steps. The value of Young’s modulus and Poisson’s ratio used in the analytical solution are 786 MPa and 0.2 for the diamond indenter, and 170 MPa and 0.2 for the silicon specimen, respectively. The interatomic potential used in the MD is that proposed by Tersoff [11]. 4. Results and discussion The deformation of the MD model during indentation is shown in Fig. 2 by taking a slice with a thickness 2a. It is shown in the figure that the model is inclined along the slope of the specimen surface outside the indenter and is being stretched in the r direction with a shear deformation in the r and z directions. It can be seen in Fig. 2 that the model in (a) is deformed only elastically but the model in (b) contains a crystal structure different from that of the diamond structure around the center of the model. It is also shown in (b) that amorphous areas are
Fig. 2. Deformation and phase transformation of the embedded MD model during indentation. (a) 7400 step, (b) 9000 step, (c) 8700 step, and (d) mechanism of amorphization.
formed at the right diagonal corners of the model. Then, the above results should be analyzed on the basis of the following facts. It is known theoretically as well as experimentally that the diamond structure of silicon transforms to a b-tin structure under a hydrostatic pressure of about 10 GPa but never transforms to an amorphous structure under any hydrostatic pressure [12]. The latter fact implies that the amorphous structure is more unstable than the b-tin structure because the b-tin structure has a lower average atomic potential energy than the amorphous structure. On the other hand, it is also known experimentally that an amorphous layer is formed on the silicon surface during indentation and/or cutting [13]. Thus, it is concluded that the phase transformation to the amorphous state in the indentation and/or cutting must be caused by the introduction of a shear stress. The latter fact indicates that the shear stress can cause a transformation from the diamond structure to an amorphous state with an increase in the average atomic potential energy. Then, a b-tin structure may appear as an interim structure during the phase transformation to the amorphous state. On the basis of the above facts, the results in Fig. 2(b) indicate that the new crystal structure that has appeared around the center of the model is b-tin with a cubic structure. In addition, it is found that the amorphous areas at the corners of the model in Fig. 2(b) have been created as the result of cross slips that occurred in the btin structure. This has been found by observing the deformation of
T. Inamura et al. / CIRP Annals - Manufacturing Technology 59 (2010) 559–562
Fig. 3. Dynamic change in density distribution due to elastic waves emitted, transmitting and reflected in the model. Waves interfere with each other. (a) 7800 step, (b) 8200 step, (c) 8600 step, and (d) 9000 step.
a red ribbon in Fig. 2(c). The process illustrated in Fig. 2(d) shows that a slip occurs in the right diagonal direction ([1 1 1] direction) owing to compression and shearing in the z direction, leaving the shrunk part of the red ribbon as indicated by a circle in Fig. 2(c). Then, another slip occurs in the left diagonal direction (also the [1 1 1] direction) owing to tension in the r direction, which bends both ends of the shrunk part of the red ribbon as shown in Fig. 2(c). Although the amorphous areas are formed during this process, the results in Fig. 2 show that defects such as microcracks cannot be created through the phase transformations. Note that the results obtained above are of high accuracy because the parameters in the interatomic potential used are those adjusted so as to reproduce the actual atomic potential energy of both the diamond and the btin structures. Thus, the results indicate that there is still another factor that must be introduced into our simulation to reproduce the ring crack initiation. It is reported that acoustic emission (AE) is observed in monocrystalline silicon when a phase change to the amorphous state takes place [14]. This fact can also be confirmed in Fig. 3, which shows the dynamic change in the density distribution due to the AE in the model. The result in Fig. 3 has been obtained by the following procedure: because the clusters in the model are controlled by the analytical solution such that they move in a quasi-static manner, only acoustic waves whose wavelengths are shorter than the cluster size and thus can travel without disturbing the cluster motions can exist in the model [9]. Then, on the basis of the wave speed and cluster size, we can calculate the time in which the wave with the longest possible wavelength travels its half wavelength. We then calculate the difference in the atomic density distributions between the start and the end of this time interval. By this procedure, we can extract the density distribution due to the acoustic wave with only the longest wavelength by doubly magnifying its amplitude, and we can approximately filter out the density distributions due to the other wave components because they almost cancel each other. It is shown in Fig. 3 that the density distribution due to the acoustic wave thus extracted changes dynamically, indicating complicated interference among waves emitted from several locations in the model and between the original waves and those reflected at the four boundaries of the model. Although the reflection occurs only at the top free surface in actual indentation, the complicated nature of the density distribution due to
561
interference will remain unchanged even in this case. On the other hand, the acoustic waves in actual indentation will propagate outward beyond the model and inward from outside the model, repeatedly interfering with each other and with waves reflected at the top free surface along the traveling path. We thus take these phenomena into account in our simulation. The most accurate and direct method of introducing the above effects will be to use an MD model sufficiently large to include the emission of waves as well as their transmission. However, because this would require a high computational cost and a long calculation time, we use an approximate method that uses the model in Fig. 2. Here, we first consider the fact that the outer periphery of the contact area between the indenter and the silicon approaches our model, moving along the surface of the silicon. We then focus on the effect of surface waves that travel along this path. Because the left-outside area of the model approaches the contact area earlier than the model itself and thus will cause a phase transformation prior to inside the model, waves emitted in the area will propagate, with some attenuation, into the model after a time determined by the wave speed and the distance from the wave source. Here, because the outer periphery of the contact area is moving toward our model, the area of the wave source will be the area of the model after an appropriate time. This means that the waves emitted in the model can be superimposed, with some attenuation, on the event that occurred in the same model at an appropriate time before. Thus, in this study, we consider three models with different time steps that are aligned side by side and approaching together the contact area. Then, we superimpose the waves extracted from the middle and left-end models on the events in the right-end model after multiplying the extracted waves by an attenuation factor and maintaining a time difference of 2000 and 4000 steps, respectively. Here, the time difference of 2000 steps is determined so that the waves travel across the model. The method of wave extraction is that described in Fig. 3, and the amplitudes of the waves extracted from the middle and left-end models have been multiplied by the reference values of 0.3 and 0.09, respectively, before superimposing them on the atomic motion in the right-end model. This is because an accurate value for the attenuation factor cannot be determined easily. The superposition has been carried out in the following manner: the extracted dynamic density distribution is evaluated at the center of each cluster and the values are transferred to the corresponding clusters in the rightend model by giving their member atoms the displacement required to cause the corresponding volume changes. These procedures are carried out at each step of the MD simulation during the period from 4000 to 6000 steps of the right-end model. The above method will contain errors in the sense that the transferred dynamic density change in the right-end model includes the effect of fictitious waves reflected at the boundaries in the original models and does not include that of actual waves reflected at the top free surface along the actual traveling path. The method also contains an error in the sense that the indentation depth of 4 mm, at which the left-outside area of the right-end model causes a phase transformation, is different from the indentation depth of 5 mm when the corresponding event occurs in the middle model. However, it is considered that the errors in both cases are not critical because the density distribution is in any case of random nature and because the microscopic states at the local area close to the outer periphery of the contact surface do not differ significantly between the indentation depths of 4 and 5 mm. Then, it is found that the distance to the nearest-neighbor atoms in the right-end model is distributed between 2.3 and 2.5 A˚ during the period without the effect of transmitted waves, whereas the corresponding value is spread between 2.1 and 2.6 A˚ when the effect of transmitted waves is introduced. The result of the simulation including the effect of transmitted waves is shown in Fig. 4 using the same slice as that in Fig. 2. It can be seen in Fig. 4 that only a simple shear deformation occurs in (a) but a polycrystal structure that is composed of a mixture of diamond and b-tin structures appears in (b). The reason why the
562
T. Inamura et al. / CIRP Annals - Manufacturing Technology 59 (2010) 559–562
forms a cross slip, as shown in (4). Here, because the cross slip occurs at a grain boundary, the atoms at the cross point cannot move freely, and thus a sufficiently high tensile stress to create a void appears locally. Then, the model collapses because the microvoids create a stress field of the typical mode-I crack front and trigger the formation of a ring crack. The reason why a ring crack occurs not at the outer periphery of the contact area but at an area outside of it is because the condition mentioned above is satisfied there first. It can also be concluded that the positions of ring cracks thus created will scatter even in the case of no preexisting defect because the dynamic stress distribution has a random nature to a certain extent owing to the complex interference between emitted and reflected waves. Finally, we point out that the time from the formation of a local polycrystal structure to the initiation of microvoids is only about 2 ps, and thus the process can be studied only by the method used in this study at present. However, we should endeavor to observe and confirm the process experimentally in the future. 5. Conclusions (1) Microvoids that trigger the formation of a ring crack can be initiated even in monocrystalline silicon with no preexisting defect. The void initiation occurs through the combined effect of dynamic stress and static stress. (2) The dynamic stress associated with acoustic waves originates from an amorphous area that is created by the static stress around the outer periphery of the contact area between an indenter and silicon. (3) The dynamic stress superimposed on the static stress field of Hertz theory locally transforms a monocrystal structure into a polycrystal structure outside the amorphous area, and then the static stress causes cross slips at the resulting grain boundaries in the area to cause microvoids. References
Fig. 4. Phase transformation to polycrystal structure followed by void initiation caused by cross slips. (a) 4300 step, (b) 5200 step, (c) 5300 step, (d) 5800 step, (e) 5290 step, and (f) mechanism of void creation.
b-tin structure appears locally and earlier than in Fig. 2 can be attributed to the fact that the dynamic compressive/tensile stresses associated with acoustic waves act in addition to the static stress, such that the combined shear stress causes a local phase transformation. Then, as shown in Fig. 4(c), a microvoid is created around a grain boundary. Other microvoids then appear at several grain boundaries, and these microvoids gradually coalesce leaving amorphous areas around them, as shown in Fig. 4(d). The ribbon mark in Fig. 4(e) indicates how the silicon has been deformed around the void in Fig. 4(c). It is shown in the figure that the ribbon undulates violently and is almost torn off, indicating that a cross slip has occurred there. Then, the deformation of the ribbon is observed during the void creation, and the mechanism of which is shown schematically in Fig. 4(f). In Fig. 4(f), a polycrystal structure is formed in (1), as indicated by different colors, and a local slip is about to occur in the direction of 45 degrees owing to shear and compression in the z direction. In (2), the local slip described above has occurred along a grain boundary. Then in (3), another slip is about to occur in the direction of 45 degrees owing to the tensile stress in the r direction. The slip described in (3)
[1] Ueda K, Sugita T, Hiraga H (1991) A J-integral Approach to Material Removal Mechanisms in Microcutting of Ceramics. Annals of the CIRP 40(1):61–64. [2] Inamura T, Shimada S, Takezawa N, Ikawa N (1999) Crack Initiation in Machining Monocrystalline Silicon. Annals of the CIRP 48(1):81–84. [3] Inamura T, Shishikura Y, Takezawa N (2009) Digital Microscope Observation of the Initial Stage of Cutting Monocrystalline Silicon. Annals of the CIRP 58(1):69–72. [4] Rentsch R, Inasaki I (1994) Molecular Dynamics Simulation for Abrasive Process. Annals of the CIRP 43(1):327–330. [5] Komanduri R, Chandrasekaran N, Raff LM (1999) Orientation Effects in Nanometric Cutting of Single Crystal Materials: An MD Approach. Annals of the CIRP 48(1):67–72. [6] Shimada S, Ikawa N, Tanaka H, Ohmori G, Uchikoshi J, Yoshinaga H (1993) Feasibility Study on Ultimate Accuracy in Microcutting Using Molecular Dynamics Simulation. Annals of the CIRP 42(1):91–94. [7] Wilks E, Wilks J (1991) Properties and Application of Diamond. ButterworthHeinemann. 192–211. [8] Harrison J, Wilks J (1978) The Hertz Indentation Test and Auerbach’s Law. Journal of Physics D 11(1):73–81. [9] Inamura T, Takezawa N, Shibuya K, Yamada K (2007) Dynamic Phenomena at Mode-I Crack Front in Silicon Simulated by Extended Molecular Dynamics. Annals of the CIRP 56(1):561–564. [10] Timoshenko SP, Goodier JN (1970) Theory of Elasticity. third edition. McGrawHill. 398–414. [11] Tersoff J (1989) Modeling Solid-State Chemistry: Interatomic Potentials for Multicomponent Systems. Physical Review B 39(8):5566–5568. [12] Hu JZ, Merkle LD, Menou CS, Spain IL (1986) Crystal Data for High-Pressure Phases of Silicon. Physical Review B 34(7):4679–4684. [13] Clarke DR, Kroll MC, Kirchner PD, Cook RF (1988) Amorphization and Conductivity of Silicon and Germanium Induced by Indentation. Physical Review Letters 60(21):2156–2159. [14] Shimada S, Inamura T, Takezawa N, Ohmori H, Sata T, Ikawa N (1995) BrittleDuctile Transition Phenomena in Microindentation and Micromachining. Annals of the CIRP 44(1):523–526.