- Email: [email protected]
Molecular dynamics simulation of amorphous silica under uniaxial tension: From bulk to nanowire
Journal of Non-Crystalline Solids 358 (2012) 3481–3487
Contents lists available at SciVerse ScienceDirect
Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/ locate/ jnoncrysol
Molecular dynamics simulation of amorphous silica under uniaxial tension: From bulk to nanowire Fenglin Yuan, Liping Huang ⁎ Department of Materials Science and Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180, United States
a r t i c l e
i n f o
Article history: Received 27 March 2012 Received in revised form 17 May 2012 Available online 14 July 2012 Keywords: Silica glass; Amorphous silica nanowire; Brittle-to-ductile transition
a b s t r a c t Molecular dynamics (MD) simulations were carried out to study bulk silica glass and amorphous silica nanowire under uniaxial tension. Periodic boundary conditions were employed to mimic infinite bulk samples. Cutting and casting methods were used to prepare nanowires. Our study shows that simulation parameters, such as system size, cooling rate, working temperature and strain rate, need to be carefully chosen in order to correctly reproduce the brittle fracture behavior of amorphous silica. The stiffness of silica glass is less sensitive to these parameters than the tensile strength and the failure strain. However, the sample density and the anomalous nonlinear elasticity of silica glass should be correctly taken into account to get an accurate estimate of its stiffness from MD simulations. Our study also shows that, with proper simulation parameters, amorphous silica nanowires down to 1 nm in radius still exhibits a brittle fracture behavior. Nanowires prepared by the cutting method have a lower stiffness and tensile strength but a higher failure strain than the cast ones, due to more surface defects generated during the cutting process at low temperatures. Defects-induced ductility could be an effective way to make less brittle nanostructures of amorphous silica. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Classical molecular dynamics (MD) simulation has been extensively used to study the deformation and fracture behavior of bulk silicate glass [1–12] and amorphous silica nanowire [13–15]. However, a straightforward comparison between MD simulations and experiments is still not always possible due to the limited time scale (~μs) and length scale (~μm) accessible by the current computation power [16]. Often an extrapolation relationship between simulations and experiments can be assumed; therefore computational experiments in MD simulations can be conducted to understand the structure and properties of materials at the atomic level. In doing so, cautions need to be taken to ensure samples are correctly prepared and test conditions are properly set up. For sample preparation, an accurate force field, a reasonable cooling rate and a large enough system size are needed to represent an amorphous system correctly. On the other hand, test conditions like working temperature and strain rate would inevitably influence the response of glass to a mechanical stimulus. However, to our best knowledge, a systematic study on different factors influencing mechanical tests of glass systems in MD is very scarce, even for the most widely studied amorphous silica. Thus, in this paper, we primarily focused on a few important factors that should be chosen properly for doing MD mechanical tests in both bulk silica glass and amorphous silica nanowire. We chose the widely used pair-wise BKS potential [17] for our MD simulations, so we can compare our results with previous studies based on it. ⁎ Corresponding author. E-mail address: [email protected] (L. Huang). 0022-3093/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2012.05.045
We also compared two different methods (cutting vs. casting) to prepare amorphous silica nanowire and studied the effect of defects on the stiffness, tensile strength, and failure strain of amorphous BKS silica nanowires. Our study will help explain some discrepancies in the literature, which are likely due to the difference in sample preparation and testing conditions. 2. Computational details 2.1. Force field and simulation techniques In order to overcome the short-range collision problem in the original BKS potential [17], a modified version proposed by Vollmayr et al. [18] was implemented in LAMMPS [19] package and used throughout this paper. The short-range interaction was truncated and shifted at 0.55 nm in order to obtain a density of ~2.29 g/cm 3 (close to experimental value of 2.20 g/cm3) of as-quenched silica glass at room temperature. For bulk systems, the Columbic interaction was calculated via the Ewald's summation technique [20] with a relative precision of 10−5 in force. For amorphous silica nanowires, Particle–Particle Particle-Mesh (PPPM) technique [21], with a relative precision of 10−4 in force, was adopted to speed up calculations. Velocity Verlet algorithm with a time step of 1.6 fs was used to integrate the Newton's equations of motion. Periodic boundary conditions (PBC) were imposed in all simulations of bulk silica glass. For nanowires, a vacuum region of 5 nm was added to each side to minimize the long-range Coulombic interaction between the nanowire in the central simulation box and its surrounding images, and the periodic boundary condition was only applied along the
3482
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
axial direction. Furthermore, Nose–Hoover thermostat [22,23] and barostat [24] were utilized to control the temperature and pressure of the sample when necessary.
2.2. Sample preparation Silica liquid was obtained by melting a 5184 atom system of α-cristobalite silica crystal at 7000 K. A NPT (constant number of atoms, pressure and temperature) quenching scheme from 7000 K to 300 K was adopted to quench silica liquid to bulk glass, with a nominal quenching rate of 1, 10, 50 or 100 K/ps. Unless otherwise specified, the default quenching rate was 10 K/ps in this study. A duplication method demonstrated in Fig. 1 was applied to study the system size effect. The as-quenched bulk silica glass with a cubic geometry was duplicated along x-, y- and z-axis for different times (z-axis is the default tensile axis). A (111) duplication means no duplication was done, or equivalently, the sample was the as-quenched silica glass. While a (113) duplication means the as-quenched sample was duplicated along the z-axis for two times without any duplication along the x- and y-axis. Structural relaxation for at least 0.1 ns at 300 K was performed after the duplication operation and before the system was subjected to uniaxial tension tests. This duplication method provided us a quick way to generate large systems without quenching every sample from the high temperature liquid and allowed us to evaluate the system size effect in mechanical tests systematically. For nanowires, as shown in Fig. 2, we either directly cut out a cylinder-shaped nanowire from the as-quenched bulk silica glass (we tentatively called it the “cutting” method) or quenched a silica liquid within a repulsive wall (the “casting” method). In the casting process, we prepared a bulk liquid using the periodic boundary conditions, and then a liquid of a cylinder-shaped shape was taken out and a repulsive wall was added to it to mimic a container in the experimental quenching process. Most studies on either amorphous or crystalline nanowires used the cutting method [14,25–27] or simply released the lateral boundary conditions of a bulk sample at room temperature [15]. Techniques similar to the casting method were used to prepare amorphous silica nanospheres [28] or metallic glass nanowires [29]. One important point worth noting is that the density of nanowires will increase and correspondingly the radius will shrink during the quenching process, so the effect of the
Fig. 2. Schematic of the cutting and casting methods for preparing amorphous silica nanowires. The gray layer indicates a repulsive wall only applied during the quenching process in the casting method.
repulsive wall gradually diminishes on cooling and has a negligible effect the structure and properties of as-quenched nanowires. 2.3. Uniaxial tension tests Uniaxial tension tests were conducted along the z-axis in both bulk silica glass and amorphous silica nanowire by rescaling the z-coordinate of each atom in a certain frequency defined by an effective strain rate. For bulk silica glass, ensemble effect (i.e., NVT vs. NPT uniaxial tension test), system size effect, cooling rate effect, strain rate effect and working temperature effect were studied systematically. We also studied the effect of surface defects on the stiffness, the tensile strength, and the failure strain of amorphous silica nanowires. 2.4. Calculation of fracture surface energy (FSE) and elastic strain energy (ESE) for bulk silica glass The fracture surface energy (FSE) was calculated by two steps: a) cutting the simulation box into two halves, and b) separating one half from the other without any structural relaxation. Total internal energy of the two halves would converge with respect to the separation distance. The fracture surface energy was calculated as the difference between the total internal energy of the original sample and the converged total internal energy of the two halves, i.e., FSE= E(two-halves) − E(original). The elastic strain energy (ESE) was approximated by ESE= 0.5 ⁎ E(ave)⁎ ε(elastic)2 ⁎ V0, where ESE is the total elastic strain energy before any plastic deformation, E(ave) is the nominal Young's modulus from a linear fitting of the stress–strain curve in the strain range 0–12%, V0 is the total volume of the sample before deformation, ε(elastic) is the elastic strain limit and taken as 12% for the BKS silica glass, consistent with the elastic strain range used by Pedone et al. [9]. 2.5. Calculation of Young's modulus of bulk silica glass To obtain the zero-strain Young's modulus of bulk silica glass, we either fitted the stress–strain curve to the Hook's law within the linear elasticity region (taken as 0–0.5% strain range) or to a third-order polynomial within the elastic limit (taken as 0–12% strain range). 2.6. Calculation of nanowire radius
Fig. 1. Schematic of the duplication method to generate large samples for uniaxial tension tests in bulk BKS silica glass. Note: 111 duplication with 5184 atoms is the as-quenched sample.
In order to obtain an accurate stress–strain curve for uniaxial tension tests of nanowires, a good way to estimate the average radius of nanowire is needed. Here we adopted an algorithm as used in experiments [15] to compute the average radius of nanowire in MD: a) a nanowire is oriented along the z-axis and divided into sections of 0.1 nm thick; b) for each section, the outmost atoms are identified and the radius is calculated by averaging the distances between the outmost atoms and
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
3483
the central axis; c) radii of all sections are averaged to give the nanowire radius. 3. Results 3.1. Uniaxial tension test in bulk silica glass Previous MD simulations of uniaxial tension tests extensively employed the NVT (constant number of atoms, volume and temperature) method, in which the lateral axes were fixed while the tensile axis was elongated, and reported a brittle fracture behavior in bulk silica glass based on BKS and CTBKS potentials [8]. However, in a recent paper by Pedone et al. [9] on uniaxial tension tests of silica glass using a pair-wise empirical potential developed by the same authors [30], a NPT (constant number of atoms, pressure and temperature) tension test or unconstraint tension test was used. In this case, the lateral axes were relaxed to keep zero lateral pressure while the tensile axis was elongated, and a less brittle fracture was reported than in the NVT tension test [9]. We observed similar behaviors in NPT and NVT tension tests of bulk BKS silica glass as shown in Fig. 3. Tested at 0.1 Tg (i.e., 300 K since Tg is ~3000 K measured in our MD simulations) with a strain rate of 10 9/s, the as-quenched 5184-atom silica glass didn't fracture even up to 50% strain in the NPT tension test, while the same system fractured rapidly after 30% strain in the NVT tension test. Such a discrepancy is attributed to the fact that geometrical constraints in lateral axes in the NVT tension test facilitate the brittle fracture by imposing a triaxial tensile stress state in the system. Therefore, the NPT ensemble mimics the experimental uniaxial tension test more realistically than the NVT ensemble, and was used for bulk silica glass in the rest of this study. Since the brittle fracture of bulk silica glass has been well known experimentally, the NPT uniaxial tension test in MD should be able to reproduce a similar behavior given the force field is accurate enough. Driven by such a conjecture, we investigated the system size effect in uniaxial tension tests of bulk silica glass by using the duplication method described in Section 2.2. Fig. 4 compares the stress–strain curves of 111- (as-quenched silica glass), 112- and 113-duplicated samples in the uniaxial NPT tension test at 0.1 Tg with a strain rate of ~109/s. It can be seen that the elastic region and tensile strength do not depend on the system size appreciably, but the failure strain changes dramatically when the system size increases. The as-quenched sample does not fracture even up to 50% strain, the 112-duplicated sample behaves in the same manner initially, and then fractures at ~36% strain. The 113-duplicated sample exhibits a very clean brittle fracture with a failure
Fig. 3. Comparison of stress–strain curves of NVT and NPT uniaxial tension tests in as-quenched BKS silica glass (111 duplication, 10 K/ps cooling rate, 0.1 Tg working temperature, 109 strain rate).
Fig. 4. Comparison of stress–strain curves of NPT uniaxial tension test in BKS silica glasses of different sizes (10 K/ps cooling rate, 0.1 Tg working temperature, 109 strain rate).
strain of 22%, similar to the results from MD simulations of NPT tension tests in silica glass by Pedone et al. [9], which were carried out on systems of 12,288 atoms at a strain rate of 109/s. In the rest of this paper, if a sample fractures (its stress drops to zero suddenly) below 25% strain, we call it “brittle”. To explain the above system size effect, two more MD simulations of 116 and 222-duplicated samples were performed and also shown in Fig. 4. The 222-duplicated sample behaves quite similarly to the as-quenched sample. In contrast, the 116-duplicated sample differs from the 113-duplicated sample only by a very small decrease in the fracture strain. Following the energy consideration of brittle fracture in classical mechanics, the fracture surface energy (FSE) and the elastic strain energy (ESE) were calculated as described in Section 2.4 and plotted in Fig. 5. For “ductile” 111-, 112- and 222-duplicated samples, the FSE is always higher than the ESE; while for “brittle” 113 and 116-duplicated samples, the FSE is far smaller than the ESE. Hence, a consistent explanation of the system size effect could be constructed by plotting the ESE/FSE ratio versus the tensile axis length (i.e., the length of the system along the tensile loading direction before any deformation occurs), as seen in Fig. 6(a). The ESE should overwhelm the FSE, i.e., the ESE/FSE ratio should exceed 1.0, in order to allow the system to undergo a brittle fracture. Otherwise, the system doesn't
Fig. 5. The fracture surface energy (FSE) and the elastic strain energy (ESE) as a function of the initial volume in different BKS silica glasses.
3484
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
~10 nm in Fig. 6(a), in order for the system to see the benefit of going through a brittle fracture. In other words, under such a condition, creating fracture surface causes the system less energy than being further elastically deformed. The key strategy to resolve such a system size issue is to make sure the tensile axis is larger than the critical length. The critical value can be also calculated from the ratio between the area density of FSE and the volume density of ESE, which is 36.608 eV/nm 2 and 3.5979 eV/nm3, from the slope in Fig. 6(b) and (c), respectively. This ratio gives a critical length of 10.2 nm, consistent with the value obtained in Fig. 6(a). Therefore, a critical tensile axis length of ~10 nm is needed for BKS silica glass to exhibit a brittle fracture in a NPT uniaxial tension test, given the lateral dimensions are big enough to ensure the periodic boundary conditions can mimic the “bulk” properties. Further increase in lateral dimensions is not necessary to reproduce the brittle fracture behavior in bulk BKS silica glass, but would inevitably increase the computational cost. After fully understanding the origin of the system size effect, we further studied other important factors, such as the cooling rate, the strain rate and the working temperature effect. All subsequent simulations in bulk silica glass utilized the 113-duplication sample to avoid any possible system size effect. To study the cooling rate effect, we compared four different samples quenched under a nominal cooling rate of 1 K/ps, 10 K/ps, 50 K/ps or 100 K/ps. Fig. 7 depicts the tensile stress– strain curves of these samples tested at 0.1 Tg with a strain rate of ~109/s. The sample quenched with the cooling rate of 100 K/ps, does not fracture up to 50% tensile strain and a brittle fracture is seen for samples quenched with 1 and 10 K/ps cooling rates, while the sample quenched under 50 K/ps behaves in between these two extremes. A glass formed under a higher cooling rate retains a liquid structure at a higher temperature, i.e., has a higher fictive temperature [31]. The higher fictive temperature state of silica glass is accompanied by a lower viscosity and a higher defect concentration [32]; both facilitate a plastic deformation instead of a brittle fracture under a uniaxial tension test. Thus, a nominal cooling rate of no more than 10 K/ps is required for reproducing the brittle fracture behavior in BKS bulk silica glass. In addition, the strain rate effect in uniaxial tension tests of BKS silica glass was investigated by conducting four different MD simulations at 0.1 Tg under strain rates of 2.48 ⁎ 10 8/s, 2.48 ⁎ 10 9/s, 2.48 ⁎ 10 10/s and 2.48⁎1011/s, spanning the ranges used by Pedone et al. [9] (108/s– 109/s) and Muralidharan et al. [8] (5⁎109/s–1011/s). The stress–strain curves of the four different simulations are shown in Fig. 8. As found by Muralidharan et al. [8] and Pedone et al. [9], the tensile strength and the failure strain increase with increasing strain rate, while the initial linear elastic region seems to be independent on the strain rate. An increase of
Fig. 6. (a) The ESE /FSE ratio versus the sample tensile axis length, (b) the FSE versus the sample cross section area, (c) the ESE versus the sample volume.
possess enough elastic strain energy to promote the brittle fracture when deformed beyond the elastic limit. The horizontal line across the ESE/FSE=1.0 marks the boundary between two regions in Fig. 6(a): in the ESE/FSEb 1.0 region, due to a small tensile axis length, a brittle fracture behavior wouldn't be observed in BKS silica even up to 50% tensile strain. Conversely, in the ESE/FSE>1.0 region, due to a large tensile axis length, a brittle fracture would be reproduced in BKS silica. In short, the geometry and size of the simulated system need to be carefully chosen to satisfy ESE/FSE>1.0, which indicates a critical tensile axis length of
Fig. 7. Stress–strain curves of NPT uniaxial tension test in BKS silica glasses quenched with different cooling rates (113 duplication, 0.1 Tg working temperature, 109 strain rate).
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
3485
3.2. Uniaxial tension test in amorphous silica nanowire What we learned in uniaxial tension tests of bulk systems helped us determine the sample preparation and test conditions for nanowires. For nanowires, it is not so critical to choose NPT or NVT tension test due to the large vacuum region, they can relax laterally in either case. Fig. 10 shows the uniaxial tension stress–strain curves for cast and cut BKS amorphous silica nanowires with a radius of 1 nm and different axial lengths. For either cast or cut nanowires, among the three samples, the elastic response is very similar. Nanowire with a small tensile axis length (e.g., as-quenched silica nanowire of 4 nm in length) shows an artificial ductility, but the silica nanowire from the 118-duplication has an axial length of 32 nm and shows a clear brittle fracture. Such a system size effect in silica nanowire is consistent with the results from bulk samples, and further proves the importance of choosing a proper system size for MD mechanical tests, especially the axial length. 4. Discussion Fig. 8. Stress–strain curves of NPT uniaxial tension test in BKS silica glass under different strain rates (113 duplication, 10 K/ps cooling rate, 0.1 Tg working temperature).
failure strain with increasing strain rate was also observed in two-point bending tests of silica glass at 77 K [33]. Fig. 8 shows that samples tested under slow strain rates display a brittle fracture, while samples tested under high strain rates don't fracture up to 50% tensile strain. A slower strain rate allows the structural rearrangement to catch up and thus produces a “close-to-equilibrium” tension test, while a higher strain would generate a “far-from-equilibrium” test condition. If the strain rate is too high, voids don't have enough time to form and coalesce to nucleate cracks, and it is hard for the sample to fracture in a brittle manner. A similar brittle fracture was observed by Pedone et al. [9] using slower strain rates. The brittle-to-ductile transition with increasing strain rate in Muralidharan et al.'s NVT tension tests [8] was not as obvious as in our NPT tests, because, as discussed above, the lateral constraints in the NVT ensemble facilitate a brittle fracture in general. Finally, Fig. 9 shows that a brittle-to-ductile transition can also be observed by increasing the working temperature in MD simulations of uniaxial tension tests under a strain rate of ~ 10 9/s, consistent with previous simulations [9] and experiments [34,35]. At 0.1 Tg, BKS silica glass is fully brittle. However, with increasing temperature, especially above 0.2 Tg, BKS silica glass doesn't fracture even up to 50% tensile strain, indicating a ductile behavior. Such a temperature-induced brittle-to-ductile transition usually originates from the increase of bond reformation rate with temperature, to be able to catch up with the bond breaking rate [36] during deformation.
Fig. 9. Stress–strain curves of NPT uniaxial tension test in BKS silica glass at different homologous temperatures (113 duplication, 10 K/ps cooling rate, 109 strain rate).
4.1. Bulk silica glass Through the above systematic study, we found the optimal range of parameters for simulating the uniaxial tension test of BKS bulk silica glass in MD to be: b10 K/ps for the cooling rate, >10 nm for the tensile axis length, b0.1 Tg for the working temperature and b 109/s for the strain rate.
Fig. 10. Stress–strain curves of uniaxial tension test in cast (a) and cut (b) BKS amorphous silica nanowires with a radius of ~1 nm (10 K/ps cooling rate, 0.1 Tg working temperature, 109 strain rate).
3486
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
After optimal conditions for mechanical tests were determined, we studied the elastic behavior of bulk silica glass. Young's modulus of bulk silica glass with the short-range cutoff of the BKS potential at 0.55 nm was calculated by: a) a third-order polynomial fit to the stress–strain curve in 0–12% strain range (shown in Fig. 11), and b) a linear fit to the stress–strain curve in 0–0.5% strain range. The first method gave a value of ~ 84 GPa and the second one gave ~ 85 GPa for the Young's modulus, higher than the experimental value of 72 GPa [37]. Our results are comparable to the Young's modulus of 75 GPa for bulk silica glass from MD simulations using the original BKS potential reported by Silva et al. [13]. However, using the same BKS potential, Muralidharan et al. obtained a value of ~ 100 GPa for the Young's modulus of bulk silica glass in their MD simulations [8]. A couple of reasons can help explain the different values obtained from different MD simulations based on the same BKS potential. As seen in Fig. 11, if the short range cutoff is increased to 1.0 nm, a higher density of ~ 2.38 g/cm 3 for as-quenched silica glass is obtained, compared with a density of 2.29 g/cm 3 using the cutoff of 0.55 nm. The zero-strain Young's modulus from the third-order polynomial fitting in the 0–12% strain range is 71.8 GPa. Therefore, the stiffness of the BKS silica glass is very sensitive to the short-range cutoff, which is usually not reported in papers. It is interesting to note from the above two simulations that a higher density leads to a lower Young's modulus in bulk silica glass. This is due to the anomalous nonlinear elasticity in silica glass, well known experimentally [38]. Fig. 12 compares the nonlinear elasticity of silica glass from our work, Pedone et al.'s MD simulations [9], as well as Gupta and Kurkjian's experiments [38]. Although the absolute value of Young's modulus is a little bit overestimated in our study, the overall nonlinear trend is similar to experiments [38] and Pedone's simulations [9]. Contrary to normal materials, the first-order nonlinear elasticity coefficient of silica glass is negative (i.e., a larger compressive strain would lead to a lower elastic modulus). Furthermore, the stiffness of silica glass is very sensitive to the strain state of the system. Taking BKS potential as an example, it can be seen from Fig. 12, between −5 to 5% strain, the Young's modulus changes from ~40 to above 100 GPa. If two different strain ranges were used to extract the stiffness, they would give very different values; even everything else is set the same in MD simulations. In our simulations, if 0–5% strain range is used to extract Young's modulus by a linear fit to the Hooke's law, the fitted Young's modulus is ~97 GPa, very close to the value of ~100 GPa from Muralidharan et al.'s MD results [8]. Thus, a very small linear strain range (e.g., 0–0.5% strain) is the key to an accurate estimate of zero-strain Young's modulus. A linear fit of the stress–strain curve in a larger tensile strain range will inevitably
Fig. 11. A third-order polynomial fit of the stress–strain curve in 0–12% strain range in uniaxial tension test of bulk BKS silica glass with a short range cutoff at 0.55 nm and 1.0 nm, respectively.
Fig. 12. A comparison of the nonlinear elasticity of silica glass between Pedone et al.’s MD simulations [9], Gupta and Kurkjian’s third-order polynomial fit to the experimental data of silica fibers [38] and this work.
overestimate it, due to the strong nonlinear elasticity in silica glass. It can be easily understood now from Fig. 12 that denser silica would have a lower stiffness. If Young's modulus were obtained from a stress– strain curve under compression, its value would be lower than that obtained under tension. This may partially explain the systematically lower Young's modulus of amorphous silica nanowire under compression in Lilian et al.'s MD simulations [14] than that in Silva et al.'s MD simulations [13] of silica nanowire under tension. Other reasons could be different potentials (Garofalini [39] vs. BKS [17]) and sample preparations (cutting vs. annealing) used in these studies. Therefore, in order to compare results from different simulations, even just for the elasticity using the sample potential, it is very important to be aware of the different cutoffs, which should be published as potential parameters. Otherwise, different cutoffs lead to different density, which in turn gives different stiffness. Another point should always be born in mind is that silica glass has the anomalous nonlinear elasticity; a comparison is only meaningful if the stiffness was obtained from the same strain state. 4.2. Amorphous silica nanowire Another important point worth noting in Fig. 10 is that silica nanowire still remains brittle even when the radius shrinks down to 1 nm. If the BKS potential is accurate enough to capture the deformation mechanism at the nanometer scale for silica glass, the above finding indicates that pristine silica nanowire probably never becomes ductile despite the large surface-to-volume ratio. In the meantime, ductile silica nanowires observed in experiments [15] are likely due to the defects creation under electron-beam radiation. A systematic study is undergoing to understand the radiation effect in the fracture behavior of BKS silica nanowire and will be published elsewhere. The large ductility observed in the companying MD simulations [15] of uniaxial tension test of a CHIK [40] (a modified version of BKS) silica nanowire (3.56 nm by 3.56 nm in cross section) may be due to the short axial length (7.12 nm) and surface defects created after relaxing the lateral boundary conditions at 300 K. These conjectures are supported by our MD results in Fig. 10(a) and (b). For either cast or cut nanowires, the tensile strength and the failure strain become higher with decreasing nanowire axial length. So a shorter nanowire tends to exhibit a larger apparent but artificial ductility. By comparing the corresponding curves in Fig. 10(a) and (b), it can be seen that the stiffness and the tensile strength are lower, but the failure strain is higher in cut nanowires than in cast ones. This is due to a higher concentration of defects (e.g., dangling bonds) on the surface of cut nanowires than on the cast ones as seen in Fig. 13. The dangling bonds from the cutting process at
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
3487
References [1] [2] [3] [4] [5]
[6] [7] [8] [9]
Fig. 13. Three-fold coordinated Si atoms (big spheres) on the surface of cast (left) and cut (right) nanowire with a radius of 1 nm. Note: O atoms are omitted for clarity; small spheres are four-fold coordinated Si atoms.
300 K don't have enough thermal energy to reform the bonds correctly. They facilitate the bond-switching events that impart flexibility to nanowire to accommodate the plastic deformation [15]. 5. Conclusions Our systematic study identified the optimal range of parameters for uniaxial tension test of BKS amorphous silica in MD simulations. To observe the brittle fracture behavior, b 10 K/ps for the cooling rate, >10 nm for the tensile axis length, b 0.1 Tg for the working temperature and ≲10 9/s for the strain rate should be used. Otherwise, different degree of artificial ductility may be observed. Moreover, due to the high surface-to-volume, it only makes sense to compare the mechanical properties of nanowires if their surface defect states are comparable. A higher concentration of defects on the surface would introduce a higher ductility, which can be exploited as a means to increase the ductility of silica nanowire, for example, by using energetic particle bombarding. Our study further points out that the density and the anomalous nonlinear elasticity of silica glass should be taken into account correctly to get an accurate estimate of its stiffness in MD simulations. Acknowledgment This work is supported by the National Science Foundation under Grant No. DMR-907076. We thank Drs. M. Tomozawa and Y. Shi at Rensselaer for stimulating discussions.
[10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40]
T. Soules, R. Busbey, J. Chem. Phys. 78 (1983) 6307–6316. R. Ochoa, J. Simmons, J. Non-Cryst. Solids 75 (1985) 413–418. J. Kieffer, C. Angell, J. Non-Cryst. Solids 106 (1988) 336–342. R. Ochoa, T. Swiler, J. Simmons, J. Non-Cryst. Solids 128 (1991) 57–68. L. Van Brutzel, C. Rountree, R. Kalia, A. Nakano, P. Vashishta, In: in: S. Komarneni, J. Parker, R. Vaia, G. Lu, J. Matsushita (Eds.), Nanophase and Nanocomposite Materials IV, vol. 703, Materials Research Society, Warrendale, 2002, pp. 117–122. K. Trachenko, M. Dove, Phys. Rev. B 67 (2003). L. Huang, J. Kieffer, Phys. Rev. B 69 (2004) 224203. K. Muralidharan, J.H. Simmons, P.A. Deymier, K. Runge, J. Non-Cryst. Solids 351 (2005) 1532–1542. A. Pedone, G. Malavasi, M.C. Menziani, U. Segre, A.N. Cormack, Chem. Mater. 20 (2008) 4356–4366. A. Pedone, J. Phys. Chem. C 113 (2009) 20773–20784. T. Taniguchi, S. Ito, J. Ceram. Soc. Jpn. 116 (2008) 885–889. A. Takada, J. Ceram. Soc. Jpn. 116 (2008) 880–884. E.C.C.M. Silva, L. Tong, S. Yip, K.J. Van Vliet, Small 2 (2006) 239–243. L. Davila, V. Leppert, E. Bringa, Scr. Mater. 60 (2009) 843–846. K. Zheng, C. Wang, Y.-Q. Cheng, Y. Yue, X. Han, Z. Zhang, Z. Shan, S.X. Mao, M. Ye, Y. Yin, E. Ma, Nat. Commun. 1 (2010). M. Rieth, W. Schommers, P.-G. de Gennes, N.P. Laureate, Handbook of Theoretical and Computational Nanotechnology, 1st ed. American Scientific Publishers, 2006. B. Vanbeest, G. Kramer, R. Vansanten, Phys. Rev. Lett. 64 (1990) 1955–1958. K. Vollmayr, W. Kob, K. Binder, Phys. Rev. B 54 (1996) 15808–15827. S. Plimpton, J. Comput. Phys. 117 (1995) 1–19. P.P. Ewald, Ann. Phys. 369 (1921) 253–287. R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles, Taylor & Francis, 1989. S. Nose, J. Chem. Phys. 81 (1984) 511–519. S. Nose, Mol. Phys. 52 (1984) 255–268. W. Shinoda, M. Shiga, M. Mikami, Phys. Rev. B 69 (2004) 134103. H. Liang, M. Upmanyu, H. Huang, Phys. Rev. B 71 (2005) 241403. T.-Y. Zhang, M. Luo, W.K. Chan, J. Appl. Phys. 103 (2008). H. Park, J. Zimmerman, Phys. Rev. B 72 (2005). A. Roder, W. Kob, K. Binder, J. Chem. Phys. 114 (2001) 7602–7614. Y. Shi, Appl. Phys. Lett. 96 (2010) 121909. A. Pedone, G. Malavasi, M.C. Menziani, A.N. Cormack, U. Segre, J. Phys. Chem. B 110 (2006) 11780–11795. C. Moynihan, A. Easteal, M. Debolt, J. Tucker, J. Am. Ceram. Soc. 59 (1976) 12–16. M. Tomozawa, C.-Y. Li, T.M. Gross, J. Non-Cryst. Solids 356 (2010) 1194–1197. C.R. Kurkjian, P.K. Gupta, R.K. Brow, Int. J. Appl. Glass Sci. 1 (2010) 27–37. E. Le Bourhis, T. Rouxel, J. Non-Cryst. Solids 316 (2003) 153–159. T. Rouxel, J. Sangleboeuf, J. Non-Cryst. Solids 271 (2000) 224–235. N. Mott, Philos. Mag. B-Phys. Condens. Matter Stat. Mech. Electron. Opt. Magn. Prop. 56 (1987) 257–262. H. Scholze, Glass: Nature, Structure, and Properties, Springer-Verlag, 1991. P.K. Gupta, C.R. Kurkjian, J. Non-Cryst. Solids 351 (2005) 2324–2328. E. Webb, S. Garofalini, J. Chem. Phys. 101 (1994) 10101–10106. A. Carre, J. Horbach, S. Ispas, W. Kob, Europhys. Lett. 82 (2008) 17001.
Contents lists available at SciVerse ScienceDirect
Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/ locate/ jnoncrysol
Molecular dynamics simulation of amorphous silica under uniaxial tension: From bulk to nanowire Fenglin Yuan, Liping Huang ⁎ Department of Materials Science and Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180, United States
a r t i c l e
i n f o
Article history: Received 27 March 2012 Received in revised form 17 May 2012 Available online 14 July 2012 Keywords: Silica glass; Amorphous silica nanowire; Brittle-to-ductile transition
a b s t r a c t Molecular dynamics (MD) simulations were carried out to study bulk silica glass and amorphous silica nanowire under uniaxial tension. Periodic boundary conditions were employed to mimic infinite bulk samples. Cutting and casting methods were used to prepare nanowires. Our study shows that simulation parameters, such as system size, cooling rate, working temperature and strain rate, need to be carefully chosen in order to correctly reproduce the brittle fracture behavior of amorphous silica. The stiffness of silica glass is less sensitive to these parameters than the tensile strength and the failure strain. However, the sample density and the anomalous nonlinear elasticity of silica glass should be correctly taken into account to get an accurate estimate of its stiffness from MD simulations. Our study also shows that, with proper simulation parameters, amorphous silica nanowires down to 1 nm in radius still exhibits a brittle fracture behavior. Nanowires prepared by the cutting method have a lower stiffness and tensile strength but a higher failure strain than the cast ones, due to more surface defects generated during the cutting process at low temperatures. Defects-induced ductility could be an effective way to make less brittle nanostructures of amorphous silica. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Classical molecular dynamics (MD) simulation has been extensively used to study the deformation and fracture behavior of bulk silicate glass [1–12] and amorphous silica nanowire [13–15]. However, a straightforward comparison between MD simulations and experiments is still not always possible due to the limited time scale (~μs) and length scale (~μm) accessible by the current computation power [16]. Often an extrapolation relationship between simulations and experiments can be assumed; therefore computational experiments in MD simulations can be conducted to understand the structure and properties of materials at the atomic level. In doing so, cautions need to be taken to ensure samples are correctly prepared and test conditions are properly set up. For sample preparation, an accurate force field, a reasonable cooling rate and a large enough system size are needed to represent an amorphous system correctly. On the other hand, test conditions like working temperature and strain rate would inevitably influence the response of glass to a mechanical stimulus. However, to our best knowledge, a systematic study on different factors influencing mechanical tests of glass systems in MD is very scarce, even for the most widely studied amorphous silica. Thus, in this paper, we primarily focused on a few important factors that should be chosen properly for doing MD mechanical tests in both bulk silica glass and amorphous silica nanowire. We chose the widely used pair-wise BKS potential [17] for our MD simulations, so we can compare our results with previous studies based on it. ⁎ Corresponding author. E-mail address: [email protected] (L. Huang). 0022-3093/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2012.05.045
We also compared two different methods (cutting vs. casting) to prepare amorphous silica nanowire and studied the effect of defects on the stiffness, tensile strength, and failure strain of amorphous BKS silica nanowires. Our study will help explain some discrepancies in the literature, which are likely due to the difference in sample preparation and testing conditions. 2. Computational details 2.1. Force field and simulation techniques In order to overcome the short-range collision problem in the original BKS potential [17], a modified version proposed by Vollmayr et al. [18] was implemented in LAMMPS [19] package and used throughout this paper. The short-range interaction was truncated and shifted at 0.55 nm in order to obtain a density of ~2.29 g/cm 3 (close to experimental value of 2.20 g/cm3) of as-quenched silica glass at room temperature. For bulk systems, the Columbic interaction was calculated via the Ewald's summation technique [20] with a relative precision of 10−5 in force. For amorphous silica nanowires, Particle–Particle Particle-Mesh (PPPM) technique [21], with a relative precision of 10−4 in force, was adopted to speed up calculations. Velocity Verlet algorithm with a time step of 1.6 fs was used to integrate the Newton's equations of motion. Periodic boundary conditions (PBC) were imposed in all simulations of bulk silica glass. For nanowires, a vacuum region of 5 nm was added to each side to minimize the long-range Coulombic interaction between the nanowire in the central simulation box and its surrounding images, and the periodic boundary condition was only applied along the
3482
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
axial direction. Furthermore, Nose–Hoover thermostat [22,23] and barostat [24] were utilized to control the temperature and pressure of the sample when necessary.
2.2. Sample preparation Silica liquid was obtained by melting a 5184 atom system of α-cristobalite silica crystal at 7000 K. A NPT (constant number of atoms, pressure and temperature) quenching scheme from 7000 K to 300 K was adopted to quench silica liquid to bulk glass, with a nominal quenching rate of 1, 10, 50 or 100 K/ps. Unless otherwise specified, the default quenching rate was 10 K/ps in this study. A duplication method demonstrated in Fig. 1 was applied to study the system size effect. The as-quenched bulk silica glass with a cubic geometry was duplicated along x-, y- and z-axis for different times (z-axis is the default tensile axis). A (111) duplication means no duplication was done, or equivalently, the sample was the as-quenched silica glass. While a (113) duplication means the as-quenched sample was duplicated along the z-axis for two times without any duplication along the x- and y-axis. Structural relaxation for at least 0.1 ns at 300 K was performed after the duplication operation and before the system was subjected to uniaxial tension tests. This duplication method provided us a quick way to generate large systems without quenching every sample from the high temperature liquid and allowed us to evaluate the system size effect in mechanical tests systematically. For nanowires, as shown in Fig. 2, we either directly cut out a cylinder-shaped nanowire from the as-quenched bulk silica glass (we tentatively called it the “cutting” method) or quenched a silica liquid within a repulsive wall (the “casting” method). In the casting process, we prepared a bulk liquid using the periodic boundary conditions, and then a liquid of a cylinder-shaped shape was taken out and a repulsive wall was added to it to mimic a container in the experimental quenching process. Most studies on either amorphous or crystalline nanowires used the cutting method [14,25–27] or simply released the lateral boundary conditions of a bulk sample at room temperature [15]. Techniques similar to the casting method were used to prepare amorphous silica nanospheres [28] or metallic glass nanowires [29]. One important point worth noting is that the density of nanowires will increase and correspondingly the radius will shrink during the quenching process, so the effect of the
Fig. 2. Schematic of the cutting and casting methods for preparing amorphous silica nanowires. The gray layer indicates a repulsive wall only applied during the quenching process in the casting method.
repulsive wall gradually diminishes on cooling and has a negligible effect the structure and properties of as-quenched nanowires. 2.3. Uniaxial tension tests Uniaxial tension tests were conducted along the z-axis in both bulk silica glass and amorphous silica nanowire by rescaling the z-coordinate of each atom in a certain frequency defined by an effective strain rate. For bulk silica glass, ensemble effect (i.e., NVT vs. NPT uniaxial tension test), system size effect, cooling rate effect, strain rate effect and working temperature effect were studied systematically. We also studied the effect of surface defects on the stiffness, the tensile strength, and the failure strain of amorphous silica nanowires. 2.4. Calculation of fracture surface energy (FSE) and elastic strain energy (ESE) for bulk silica glass The fracture surface energy (FSE) was calculated by two steps: a) cutting the simulation box into two halves, and b) separating one half from the other without any structural relaxation. Total internal energy of the two halves would converge with respect to the separation distance. The fracture surface energy was calculated as the difference between the total internal energy of the original sample and the converged total internal energy of the two halves, i.e., FSE= E(two-halves) − E(original). The elastic strain energy (ESE) was approximated by ESE= 0.5 ⁎ E(ave)⁎ ε(elastic)2 ⁎ V0, where ESE is the total elastic strain energy before any plastic deformation, E(ave) is the nominal Young's modulus from a linear fitting of the stress–strain curve in the strain range 0–12%, V0 is the total volume of the sample before deformation, ε(elastic) is the elastic strain limit and taken as 12% for the BKS silica glass, consistent with the elastic strain range used by Pedone et al. [9]. 2.5. Calculation of Young's modulus of bulk silica glass To obtain the zero-strain Young's modulus of bulk silica glass, we either fitted the stress–strain curve to the Hook's law within the linear elasticity region (taken as 0–0.5% strain range) or to a third-order polynomial within the elastic limit (taken as 0–12% strain range). 2.6. Calculation of nanowire radius
Fig. 1. Schematic of the duplication method to generate large samples for uniaxial tension tests in bulk BKS silica glass. Note: 111 duplication with 5184 atoms is the as-quenched sample.
In order to obtain an accurate stress–strain curve for uniaxial tension tests of nanowires, a good way to estimate the average radius of nanowire is needed. Here we adopted an algorithm as used in experiments [15] to compute the average radius of nanowire in MD: a) a nanowire is oriented along the z-axis and divided into sections of 0.1 nm thick; b) for each section, the outmost atoms are identified and the radius is calculated by averaging the distances between the outmost atoms and
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
3483
the central axis; c) radii of all sections are averaged to give the nanowire radius. 3. Results 3.1. Uniaxial tension test in bulk silica glass Previous MD simulations of uniaxial tension tests extensively employed the NVT (constant number of atoms, volume and temperature) method, in which the lateral axes were fixed while the tensile axis was elongated, and reported a brittle fracture behavior in bulk silica glass based on BKS and CTBKS potentials [8]. However, in a recent paper by Pedone et al. [9] on uniaxial tension tests of silica glass using a pair-wise empirical potential developed by the same authors [30], a NPT (constant number of atoms, pressure and temperature) tension test or unconstraint tension test was used. In this case, the lateral axes were relaxed to keep zero lateral pressure while the tensile axis was elongated, and a less brittle fracture was reported than in the NVT tension test [9]. We observed similar behaviors in NPT and NVT tension tests of bulk BKS silica glass as shown in Fig. 3. Tested at 0.1 Tg (i.e., 300 K since Tg is ~3000 K measured in our MD simulations) with a strain rate of 10 9/s, the as-quenched 5184-atom silica glass didn't fracture even up to 50% strain in the NPT tension test, while the same system fractured rapidly after 30% strain in the NVT tension test. Such a discrepancy is attributed to the fact that geometrical constraints in lateral axes in the NVT tension test facilitate the brittle fracture by imposing a triaxial tensile stress state in the system. Therefore, the NPT ensemble mimics the experimental uniaxial tension test more realistically than the NVT ensemble, and was used for bulk silica glass in the rest of this study. Since the brittle fracture of bulk silica glass has been well known experimentally, the NPT uniaxial tension test in MD should be able to reproduce a similar behavior given the force field is accurate enough. Driven by such a conjecture, we investigated the system size effect in uniaxial tension tests of bulk silica glass by using the duplication method described in Section 2.2. Fig. 4 compares the stress–strain curves of 111- (as-quenched silica glass), 112- and 113-duplicated samples in the uniaxial NPT tension test at 0.1 Tg with a strain rate of ~109/s. It can be seen that the elastic region and tensile strength do not depend on the system size appreciably, but the failure strain changes dramatically when the system size increases. The as-quenched sample does not fracture even up to 50% strain, the 112-duplicated sample behaves in the same manner initially, and then fractures at ~36% strain. The 113-duplicated sample exhibits a very clean brittle fracture with a failure
Fig. 3. Comparison of stress–strain curves of NVT and NPT uniaxial tension tests in as-quenched BKS silica glass (111 duplication, 10 K/ps cooling rate, 0.1 Tg working temperature, 109 strain rate).
Fig. 4. Comparison of stress–strain curves of NPT uniaxial tension test in BKS silica glasses of different sizes (10 K/ps cooling rate, 0.1 Tg working temperature, 109 strain rate).
strain of 22%, similar to the results from MD simulations of NPT tension tests in silica glass by Pedone et al. [9], which were carried out on systems of 12,288 atoms at a strain rate of 109/s. In the rest of this paper, if a sample fractures (its stress drops to zero suddenly) below 25% strain, we call it “brittle”. To explain the above system size effect, two more MD simulations of 116 and 222-duplicated samples were performed and also shown in Fig. 4. The 222-duplicated sample behaves quite similarly to the as-quenched sample. In contrast, the 116-duplicated sample differs from the 113-duplicated sample only by a very small decrease in the fracture strain. Following the energy consideration of brittle fracture in classical mechanics, the fracture surface energy (FSE) and the elastic strain energy (ESE) were calculated as described in Section 2.4 and plotted in Fig. 5. For “ductile” 111-, 112- and 222-duplicated samples, the FSE is always higher than the ESE; while for “brittle” 113 and 116-duplicated samples, the FSE is far smaller than the ESE. Hence, a consistent explanation of the system size effect could be constructed by plotting the ESE/FSE ratio versus the tensile axis length (i.e., the length of the system along the tensile loading direction before any deformation occurs), as seen in Fig. 6(a). The ESE should overwhelm the FSE, i.e., the ESE/FSE ratio should exceed 1.0, in order to allow the system to undergo a brittle fracture. Otherwise, the system doesn't
Fig. 5. The fracture surface energy (FSE) and the elastic strain energy (ESE) as a function of the initial volume in different BKS silica glasses.
3484
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
~10 nm in Fig. 6(a), in order for the system to see the benefit of going through a brittle fracture. In other words, under such a condition, creating fracture surface causes the system less energy than being further elastically deformed. The key strategy to resolve such a system size issue is to make sure the tensile axis is larger than the critical length. The critical value can be also calculated from the ratio between the area density of FSE and the volume density of ESE, which is 36.608 eV/nm 2 and 3.5979 eV/nm3, from the slope in Fig. 6(b) and (c), respectively. This ratio gives a critical length of 10.2 nm, consistent with the value obtained in Fig. 6(a). Therefore, a critical tensile axis length of ~10 nm is needed for BKS silica glass to exhibit a brittle fracture in a NPT uniaxial tension test, given the lateral dimensions are big enough to ensure the periodic boundary conditions can mimic the “bulk” properties. Further increase in lateral dimensions is not necessary to reproduce the brittle fracture behavior in bulk BKS silica glass, but would inevitably increase the computational cost. After fully understanding the origin of the system size effect, we further studied other important factors, such as the cooling rate, the strain rate and the working temperature effect. All subsequent simulations in bulk silica glass utilized the 113-duplication sample to avoid any possible system size effect. To study the cooling rate effect, we compared four different samples quenched under a nominal cooling rate of 1 K/ps, 10 K/ps, 50 K/ps or 100 K/ps. Fig. 7 depicts the tensile stress– strain curves of these samples tested at 0.1 Tg with a strain rate of ~109/s. The sample quenched with the cooling rate of 100 K/ps, does not fracture up to 50% tensile strain and a brittle fracture is seen for samples quenched with 1 and 10 K/ps cooling rates, while the sample quenched under 50 K/ps behaves in between these two extremes. A glass formed under a higher cooling rate retains a liquid structure at a higher temperature, i.e., has a higher fictive temperature [31]. The higher fictive temperature state of silica glass is accompanied by a lower viscosity and a higher defect concentration [32]; both facilitate a plastic deformation instead of a brittle fracture under a uniaxial tension test. Thus, a nominal cooling rate of no more than 10 K/ps is required for reproducing the brittle fracture behavior in BKS bulk silica glass. In addition, the strain rate effect in uniaxial tension tests of BKS silica glass was investigated by conducting four different MD simulations at 0.1 Tg under strain rates of 2.48 ⁎ 10 8/s, 2.48 ⁎ 10 9/s, 2.48 ⁎ 10 10/s and 2.48⁎1011/s, spanning the ranges used by Pedone et al. [9] (108/s– 109/s) and Muralidharan et al. [8] (5⁎109/s–1011/s). The stress–strain curves of the four different simulations are shown in Fig. 8. As found by Muralidharan et al. [8] and Pedone et al. [9], the tensile strength and the failure strain increase with increasing strain rate, while the initial linear elastic region seems to be independent on the strain rate. An increase of
Fig. 6. (a) The ESE /FSE ratio versus the sample tensile axis length, (b) the FSE versus the sample cross section area, (c) the ESE versus the sample volume.
possess enough elastic strain energy to promote the brittle fracture when deformed beyond the elastic limit. The horizontal line across the ESE/FSE=1.0 marks the boundary between two regions in Fig. 6(a): in the ESE/FSEb 1.0 region, due to a small tensile axis length, a brittle fracture behavior wouldn't be observed in BKS silica even up to 50% tensile strain. Conversely, in the ESE/FSE>1.0 region, due to a large tensile axis length, a brittle fracture would be reproduced in BKS silica. In short, the geometry and size of the simulated system need to be carefully chosen to satisfy ESE/FSE>1.0, which indicates a critical tensile axis length of
Fig. 7. Stress–strain curves of NPT uniaxial tension test in BKS silica glasses quenched with different cooling rates (113 duplication, 0.1 Tg working temperature, 109 strain rate).
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
3485
3.2. Uniaxial tension test in amorphous silica nanowire What we learned in uniaxial tension tests of bulk systems helped us determine the sample preparation and test conditions for nanowires. For nanowires, it is not so critical to choose NPT or NVT tension test due to the large vacuum region, they can relax laterally in either case. Fig. 10 shows the uniaxial tension stress–strain curves for cast and cut BKS amorphous silica nanowires with a radius of 1 nm and different axial lengths. For either cast or cut nanowires, among the three samples, the elastic response is very similar. Nanowire with a small tensile axis length (e.g., as-quenched silica nanowire of 4 nm in length) shows an artificial ductility, but the silica nanowire from the 118-duplication has an axial length of 32 nm and shows a clear brittle fracture. Such a system size effect in silica nanowire is consistent with the results from bulk samples, and further proves the importance of choosing a proper system size for MD mechanical tests, especially the axial length. 4. Discussion Fig. 8. Stress–strain curves of NPT uniaxial tension test in BKS silica glass under different strain rates (113 duplication, 10 K/ps cooling rate, 0.1 Tg working temperature).
failure strain with increasing strain rate was also observed in two-point bending tests of silica glass at 77 K [33]. Fig. 8 shows that samples tested under slow strain rates display a brittle fracture, while samples tested under high strain rates don't fracture up to 50% tensile strain. A slower strain rate allows the structural rearrangement to catch up and thus produces a “close-to-equilibrium” tension test, while a higher strain would generate a “far-from-equilibrium” test condition. If the strain rate is too high, voids don't have enough time to form and coalesce to nucleate cracks, and it is hard for the sample to fracture in a brittle manner. A similar brittle fracture was observed by Pedone et al. [9] using slower strain rates. The brittle-to-ductile transition with increasing strain rate in Muralidharan et al.'s NVT tension tests [8] was not as obvious as in our NPT tests, because, as discussed above, the lateral constraints in the NVT ensemble facilitate a brittle fracture in general. Finally, Fig. 9 shows that a brittle-to-ductile transition can also be observed by increasing the working temperature in MD simulations of uniaxial tension tests under a strain rate of ~ 10 9/s, consistent with previous simulations [9] and experiments [34,35]. At 0.1 Tg, BKS silica glass is fully brittle. However, with increasing temperature, especially above 0.2 Tg, BKS silica glass doesn't fracture even up to 50% tensile strain, indicating a ductile behavior. Such a temperature-induced brittle-to-ductile transition usually originates from the increase of bond reformation rate with temperature, to be able to catch up with the bond breaking rate [36] during deformation.
Fig. 9. Stress–strain curves of NPT uniaxial tension test in BKS silica glass at different homologous temperatures (113 duplication, 10 K/ps cooling rate, 109 strain rate).
4.1. Bulk silica glass Through the above systematic study, we found the optimal range of parameters for simulating the uniaxial tension test of BKS bulk silica glass in MD to be: b10 K/ps for the cooling rate, >10 nm for the tensile axis length, b0.1 Tg for the working temperature and b 109/s for the strain rate.
Fig. 10. Stress–strain curves of uniaxial tension test in cast (a) and cut (b) BKS amorphous silica nanowires with a radius of ~1 nm (10 K/ps cooling rate, 0.1 Tg working temperature, 109 strain rate).
3486
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
After optimal conditions for mechanical tests were determined, we studied the elastic behavior of bulk silica glass. Young's modulus of bulk silica glass with the short-range cutoff of the BKS potential at 0.55 nm was calculated by: a) a third-order polynomial fit to the stress–strain curve in 0–12% strain range (shown in Fig. 11), and b) a linear fit to the stress–strain curve in 0–0.5% strain range. The first method gave a value of ~ 84 GPa and the second one gave ~ 85 GPa for the Young's modulus, higher than the experimental value of 72 GPa [37]. Our results are comparable to the Young's modulus of 75 GPa for bulk silica glass from MD simulations using the original BKS potential reported by Silva et al. [13]. However, using the same BKS potential, Muralidharan et al. obtained a value of ~ 100 GPa for the Young's modulus of bulk silica glass in their MD simulations [8]. A couple of reasons can help explain the different values obtained from different MD simulations based on the same BKS potential. As seen in Fig. 11, if the short range cutoff is increased to 1.0 nm, a higher density of ~ 2.38 g/cm 3 for as-quenched silica glass is obtained, compared with a density of 2.29 g/cm 3 using the cutoff of 0.55 nm. The zero-strain Young's modulus from the third-order polynomial fitting in the 0–12% strain range is 71.8 GPa. Therefore, the stiffness of the BKS silica glass is very sensitive to the short-range cutoff, which is usually not reported in papers. It is interesting to note from the above two simulations that a higher density leads to a lower Young's modulus in bulk silica glass. This is due to the anomalous nonlinear elasticity in silica glass, well known experimentally [38]. Fig. 12 compares the nonlinear elasticity of silica glass from our work, Pedone et al.'s MD simulations [9], as well as Gupta and Kurkjian's experiments [38]. Although the absolute value of Young's modulus is a little bit overestimated in our study, the overall nonlinear trend is similar to experiments [38] and Pedone's simulations [9]. Contrary to normal materials, the first-order nonlinear elasticity coefficient of silica glass is negative (i.e., a larger compressive strain would lead to a lower elastic modulus). Furthermore, the stiffness of silica glass is very sensitive to the strain state of the system. Taking BKS potential as an example, it can be seen from Fig. 12, between −5 to 5% strain, the Young's modulus changes from ~40 to above 100 GPa. If two different strain ranges were used to extract the stiffness, they would give very different values; even everything else is set the same in MD simulations. In our simulations, if 0–5% strain range is used to extract Young's modulus by a linear fit to the Hooke's law, the fitted Young's modulus is ~97 GPa, very close to the value of ~100 GPa from Muralidharan et al.'s MD results [8]. Thus, a very small linear strain range (e.g., 0–0.5% strain) is the key to an accurate estimate of zero-strain Young's modulus. A linear fit of the stress–strain curve in a larger tensile strain range will inevitably
Fig. 11. A third-order polynomial fit of the stress–strain curve in 0–12% strain range in uniaxial tension test of bulk BKS silica glass with a short range cutoff at 0.55 nm and 1.0 nm, respectively.
Fig. 12. A comparison of the nonlinear elasticity of silica glass between Pedone et al.’s MD simulations [9], Gupta and Kurkjian’s third-order polynomial fit to the experimental data of silica fibers [38] and this work.
overestimate it, due to the strong nonlinear elasticity in silica glass. It can be easily understood now from Fig. 12 that denser silica would have a lower stiffness. If Young's modulus were obtained from a stress– strain curve under compression, its value would be lower than that obtained under tension. This may partially explain the systematically lower Young's modulus of amorphous silica nanowire under compression in Lilian et al.'s MD simulations [14] than that in Silva et al.'s MD simulations [13] of silica nanowire under tension. Other reasons could be different potentials (Garofalini [39] vs. BKS [17]) and sample preparations (cutting vs. annealing) used in these studies. Therefore, in order to compare results from different simulations, even just for the elasticity using the sample potential, it is very important to be aware of the different cutoffs, which should be published as potential parameters. Otherwise, different cutoffs lead to different density, which in turn gives different stiffness. Another point should always be born in mind is that silica glass has the anomalous nonlinear elasticity; a comparison is only meaningful if the stiffness was obtained from the same strain state. 4.2. Amorphous silica nanowire Another important point worth noting in Fig. 10 is that silica nanowire still remains brittle even when the radius shrinks down to 1 nm. If the BKS potential is accurate enough to capture the deformation mechanism at the nanometer scale for silica glass, the above finding indicates that pristine silica nanowire probably never becomes ductile despite the large surface-to-volume ratio. In the meantime, ductile silica nanowires observed in experiments [15] are likely due to the defects creation under electron-beam radiation. A systematic study is undergoing to understand the radiation effect in the fracture behavior of BKS silica nanowire and will be published elsewhere. The large ductility observed in the companying MD simulations [15] of uniaxial tension test of a CHIK [40] (a modified version of BKS) silica nanowire (3.56 nm by 3.56 nm in cross section) may be due to the short axial length (7.12 nm) and surface defects created after relaxing the lateral boundary conditions at 300 K. These conjectures are supported by our MD results in Fig. 10(a) and (b). For either cast or cut nanowires, the tensile strength and the failure strain become higher with decreasing nanowire axial length. So a shorter nanowire tends to exhibit a larger apparent but artificial ductility. By comparing the corresponding curves in Fig. 10(a) and (b), it can be seen that the stiffness and the tensile strength are lower, but the failure strain is higher in cut nanowires than in cast ones. This is due to a higher concentration of defects (e.g., dangling bonds) on the surface of cut nanowires than on the cast ones as seen in Fig. 13. The dangling bonds from the cutting process at
F. Yuan, L. Huang / Journal of Non-Crystalline Solids 358 (2012) 3481–3487
3487
References [1] [2] [3] [4] [5]
[6] [7] [8] [9]
Fig. 13. Three-fold coordinated Si atoms (big spheres) on the surface of cast (left) and cut (right) nanowire with a radius of 1 nm. Note: O atoms are omitted for clarity; small spheres are four-fold coordinated Si atoms.
300 K don't have enough thermal energy to reform the bonds correctly. They facilitate the bond-switching events that impart flexibility to nanowire to accommodate the plastic deformation [15]. 5. Conclusions Our systematic study identified the optimal range of parameters for uniaxial tension test of BKS amorphous silica in MD simulations. To observe the brittle fracture behavior, b 10 K/ps for the cooling rate, >10 nm for the tensile axis length, b 0.1 Tg for the working temperature and ≲10 9/s for the strain rate should be used. Otherwise, different degree of artificial ductility may be observed. Moreover, due to the high surface-to-volume, it only makes sense to compare the mechanical properties of nanowires if their surface defect states are comparable. A higher concentration of defects on the surface would introduce a higher ductility, which can be exploited as a means to increase the ductility of silica nanowire, for example, by using energetic particle bombarding. Our study further points out that the density and the anomalous nonlinear elasticity of silica glass should be taken into account correctly to get an accurate estimate of its stiffness in MD simulations. Acknowledgment This work is supported by the National Science Foundation under Grant No. DMR-907076. We thank Drs. M. Tomozawa and Y. Shi at Rensselaer for stimulating discussions.
[10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40]
T. Soules, R. Busbey, J. Chem. Phys. 78 (1983) 6307–6316. R. Ochoa, J. Simmons, J. Non-Cryst. Solids 75 (1985) 413–418. J. Kieffer, C. Angell, J. Non-Cryst. Solids 106 (1988) 336–342. R. Ochoa, T. Swiler, J. Simmons, J. Non-Cryst. Solids 128 (1991) 57–68. L. Van Brutzel, C. Rountree, R. Kalia, A. Nakano, P. Vashishta, In: in: S. Komarneni, J. Parker, R. Vaia, G. Lu, J. Matsushita (Eds.), Nanophase and Nanocomposite Materials IV, vol. 703, Materials Research Society, Warrendale, 2002, pp. 117–122. K. Trachenko, M. Dove, Phys. Rev. B 67 (2003). L. Huang, J. Kieffer, Phys. Rev. B 69 (2004) 224203. K. Muralidharan, J.H. Simmons, P.A. Deymier, K. Runge, J. Non-Cryst. Solids 351 (2005) 1532–1542. A. Pedone, G. Malavasi, M.C. Menziani, U. Segre, A.N. Cormack, Chem. Mater. 20 (2008) 4356–4366. A. Pedone, J. Phys. Chem. C 113 (2009) 20773–20784. T. Taniguchi, S. Ito, J. Ceram. Soc. Jpn. 116 (2008) 885–889. A. Takada, J. Ceram. Soc. Jpn. 116 (2008) 880–884. E.C.C.M. Silva, L. Tong, S. Yip, K.J. Van Vliet, Small 2 (2006) 239–243. L. Davila, V. Leppert, E. Bringa, Scr. Mater. 60 (2009) 843–846. K. Zheng, C. Wang, Y.-Q. Cheng, Y. Yue, X. Han, Z. Zhang, Z. Shan, S.X. Mao, M. Ye, Y. Yin, E. Ma, Nat. Commun. 1 (2010). M. Rieth, W. Schommers, P.-G. de Gennes, N.P. Laureate, Handbook of Theoretical and Computational Nanotechnology, 1st ed. American Scientific Publishers, 2006. B. Vanbeest, G. Kramer, R. Vansanten, Phys. Rev. Lett. 64 (1990) 1955–1958. K. Vollmayr, W. Kob, K. Binder, Phys. Rev. B 54 (1996) 15808–15827. S. Plimpton, J. Comput. Phys. 117 (1995) 1–19. P.P. Ewald, Ann. Phys. 369 (1921) 253–287. R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles, Taylor & Francis, 1989. S. Nose, J. Chem. Phys. 81 (1984) 511–519. S. Nose, Mol. Phys. 52 (1984) 255–268. W. Shinoda, M. Shiga, M. Mikami, Phys. Rev. B 69 (2004) 134103. H. Liang, M. Upmanyu, H. Huang, Phys. Rev. B 71 (2005) 241403. T.-Y. Zhang, M. Luo, W.K. Chan, J. Appl. Phys. 103 (2008). H. Park, J. Zimmerman, Phys. Rev. B 72 (2005). A. Roder, W. Kob, K. Binder, J. Chem. Phys. 114 (2001) 7602–7614. Y. Shi, Appl. Phys. Lett. 96 (2010) 121909. A. Pedone, G. Malavasi, M.C. Menziani, A.N. Cormack, U. Segre, J. Phys. Chem. B 110 (2006) 11780–11795. C. Moynihan, A. Easteal, M. Debolt, J. Tucker, J. Am. Ceram. Soc. 59 (1976) 12–16. M. Tomozawa, C.-Y. Li, T.M. Gross, J. Non-Cryst. Solids 356 (2010) 1194–1197. C.R. Kurkjian, P.K. Gupta, R.K. Brow, Int. J. Appl. Glass Sci. 1 (2010) 27–37. E. Le Bourhis, T. Rouxel, J. Non-Cryst. Solids 316 (2003) 153–159. T. Rouxel, J. Sangleboeuf, J. Non-Cryst. Solids 271 (2000) 224–235. N. Mott, Philos. Mag. B-Phys. Condens. Matter Stat. Mech. Electron. Opt. Magn. Prop. 56 (1987) 257–262. H. Scholze, Glass: Nature, Structure, and Properties, Springer-Verlag, 1991. P.K. Gupta, C.R. Kurkjian, J. Non-Cryst. Solids 351 (2005) 2324–2328. E. Webb, S. Garofalini, J. Chem. Phys. 101 (1994) 10101–10106. A. Carre, J. Horbach, S. Ispas, W. Kob, Europhys. Lett. 82 (2008) 17001.