Available online at www.sciencedirect.com
NIM B Beam Interactions with Materials & Atoms
Nuclear Instruments and Methods in Physics Research B 266 (2008) 2707–2710 www.elsevier.com/locate/nimb
Molecular dynamics study of structural changes versus deposited energy dose in a sodium borosilicate glass G. Bureau a,*, J.-M. Delaye a, S. Peuget a, G. Calas b a
DEN/DTCD/SECM, CEA Marcoule, BP 17171, 30207 Bagnols-sur-Ce`ze Cedex, France b IMPMC, 140 rue de Lourmel, 75015 Paris, France Available online 1 April 2008
Abstract The accumulation of cascades modeled by molecular dynamics in a sodium borosilicate glass allowed us to simulate the evolution of various macroscopic and structural properties up to the level of a stabilization plateau for the highest deposited nuclear energy doses. Marples’ model was used to fit the glass volume expansion to the deposited energy dose, giving the damaged volume per projectile. The volume parameter from this model approximates the cascade core volume, suggesting that the underlying mechanisms of volume expansion are contained in the cascade core and are thus related to the highest-energy events: atom ejection and thermal quenching. Ó 2008 Elsevier B.V. All rights reserved. PACS: 78.55.Qr; 61.82.Ms; 71.15.Pd Keywords: Sodium borosilicate glass; Molecular dynamics; Macroscopic properties; Structural properties; Saturation plateau; Marples’ model
1. Introduction High-level radioactive waste produced by the French nuclear industry is immobilized in a borosilicate glass matrix. Assessing the long-term behavior of nuclear glass implies evaluating the impact of cumulative alpha decay in the minor actinides it contains. When subjected to the cumulative effects of alpha decay (244Cm-doped glass specimens) or to external ion irradiation, some macroscopic properties vary appreciably with the dose [1]. Above a given dose level, the properties are observed to stabilize. To improve our understanding of these modifications, studies are carried out on glass compositions modeled by molecular dynamics (MD) [2,3] in which irradiation effects are simulated by accelerating uranium projectiles. These models are capable of qualitatively reproducing the density variations observed in actinide-doped materials [4].
*
Corresponding author. Tel.: +33 04 66 79 63 78; fax: +33 04 66 79 77
08. E-mail address:
[email protected] (G. Bureau). 0168-583X/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.nimb.2008.03.217
A molecular dynamics analysis of the property stabilization period is described here in order to understand its mechanisms. 2. Method The structure of a nuclear containment glass is chemically and physically too complex to be modeled. A simplified matrix with only three oxides (SiO2, B2O3, Na2O) in the same proportions as in the actual glass [5] was used instead. 2.1. Simulating the pristine glass by molecular dynamics A glass molar composition of 67.7% SiO2, 18.0% B2O3, and 14.2% Na2O (+0.05% UO2) was prepared by classical molecular dynamics with empirical potentials. The applicable simulation code is ‘‘homemade” code developed for this purpose. This glass is designated ‘‘CJ1” in the remainder of this text. The atomic interaction potentials have already been described in previous articles [6,7]; Born–Mayer–Huggins
2708
G. Bureau et al. / Nucl. Instr. and Meth. in Phys. Res. B 266 (2008) 2707–2710
potentials were used, supplemented by three-body terms (and by a 1/r6 term for O–U interactions). Integer charges were used to represent interactions between ions (Si4+, B3+, Na+, O2). The values used were taken from work available in the literature; some of them were fitted to the structure of glasses similar to CJ1 glass [8,9]. Three-body terms such as the following [10] were applied to the O–Si–O, O–B–O and Si–O–Si triplets to improve the local angular order of the covalent bonds. The different parameters were dependent on the nature of the atomic triplet. The values were taken from earlier work [8]. Ziegler–Biersack–Littmark (ZBL) potentials were applied to simulate short-range interactions capable of occurring during a displacement cascade [11]. The ZBL potentials were related to the pair potentials by polynomial expressions ensuring continuity of energy, forces and force derivatives. The polynomial expressions were applied ˚ and 1.0 A ˚ for atomic pairs except O–U, between 0.9 A ˚ ˚ for the O–U pair. and between 1.1 A and 1.3 A 2.2. Construction of the simulation cell The cubic simulation cell contained 64,000 atoms. The glass was prepared from a random configuration with a ˚ between atoms to avoid minimum initial distance of 1 A excessively violent forces at the outset. During the initial relaxation phase, the liquid was equilibrated at high-temperature between 4000 K and 6000 K. The atom velocities were corrected to 4000 K by a multiplication factor whenever the temperature exceeded 6000 K. This procedure allowed the temperature to reach equilibrium gradually at about 5400 K. The initial relaxation phase lasted 10,000 time steps. The liquid at equilibrium was then quenched to 293 K at a rate of 1014 K s1. A final relaxation was applied in two stages. For 10,000 time steps the glass was relaxed by regularly readjusting the atomic velocities by a multiplication factor relative to the ambient temperature, followed by relaxation for 10,000 additional time steps with unrestricted atomic velocities. Each time step was 1.0 fs throughout the glass preparation phases. The edge of the simulation cell was adjusted to obtain a final pressure near 0 Pa. The glass density at room temperature was thus 2.16 g cm3. 2.3. Molecular dynamics study of cumulative displacement cascades A series of 4 keV displacement cascades was simulated in the glass at equilibrium at room temperature to study structural changes and volume expansion resulting from ballistic collisions. Three different calculations were performed for each cascade: calculation of the cascade itself, relaxation at constant pressure to determine the new equilibrium volume of the irradiated structure, and final constant-volume relaxation to equilibrate the structure before the next cascade.
Each cascade was simulated by accelerating an initial projectile (in this case a uranium atom) to an energy of 4 keV. Translational movement was applied with periodic conditions to all the atoms in the cell to position the displacement cascade in the center of the simulation box. The cascades were simulated with variable time steps fitted along a hyster˚ per esis loop to prevent any displacement exceeding 0.034 A step. Several levels were used, and as soon as the maximum atomic velocity dropped below a level the time step was increased; conversely, when the maximum atomic velocity increased above a level (e.g. when a heavy atom collided with a light atom), the time step was diminished accordingly. ˚ thick was conDuring the cascade an outer layer 3 A strained at room temperature to adjust the atomic velocities in order to dissipate the thermal wave as it reached the edges of the cell. At the end of the cascade, the time step was 1.0 fs and the temperature had practically returned to room temperature. The typical duration of cascade is equal to about 25 ps. After being subjected to the cascade the structure was relaxed in the isothermal isobaric system at room temperature and at atmospheric pressure to determine the new equilibrium volume using the algorithm proposed by Martyna et al. [12]. The new equilibrium volume was determined by calculation over 20,000 time steps of 1.0 fs before constantvolume relaxation using the volume determined by the calculation at constant pressure. This step lasted 5000 time steps of 1.0 fs and resulted in complete relaxation of the structure before simulating the next cascade. Separate calculations were necessary for the cascade simulation and the determination of the new equilibrium volume because the pressure is not equilibrated during a cascade. Locally very high pressure zones coexist with low pressure zones, and a constant-pressure calculation of such a system would inevitably lead to significant expansion of the structure. 3. Results We analyzed the macroscopic and structural changes in CJ1 glass subjected to a series of 4 keV displacement cascades. Ninety-five cascades have been accumulated to date, with a deposited nuclear energy dose of about 4.2 1020 KeV/cm3. About 66% of the simulation cell has been touched by the core of at least one displacement cascade. 3.1. Volume expansion and depolymerization The structural changes were quantified between each cascade. Volume expansion of about 5% was observed in the CJ1 glass after 95 cumulative cascades (Fig. 1). The depolymerization index was defined as follows: X NPOLYðnÞ ¼ NC i ð0Þ NC i ðnÞ ð1Þ i¼Si;B
G. Bureau et al. / Nucl. Instr. and Meth. in Phys. Res. B 266 (2008) 2707–2710
Fig. 1. Volume expansion of the CJ1 glass versus deposited nuclear energy dose. Comparison with Marples’ model for saturation density values of 9% ± 1%.
NPOLY(n) represents the number of former-oxygen bonds that disappeared between the initial instant and the nth cascade, NCi(0) and NCi(n) are the coordination numbers of the atom i (i = Si and B) at the initial instant and after n cascades. A positive index value corresponds to depolymerization of the structure.
2709
With the cumulative cascade, the density and the volume expansion versus the deposited nuclear energy dose tends toward a saturation plateau. Depolymerization of the glass network has already been shown to be due to structural changes [3]. A mechanism by which 4-coordinate boron is converted to 3-coordinate boron can be observed as the deposited nuclear energy dose increases, resulting in a larger number of nonbridging oxygen atoms (Fig. 2(a) and (b)). The number of networkmodifying Na atoms logically increases with the number of nonbridging oxygen atoms (Fig. 2(c)), while the number of charge-compensating Na atoms diminishes, reflecting the decreasing number of Na atoms necessary to stabilize BO 4 groups. As for volume expansion, these structural changes tend to level off at a saturation plateau. The glass structural changes thus appear to cease above a certain deposited nuclear energy dose corresponding to the level at which the macroscopic properties stabilize and confirming the link between structure and properties. 3.2. Increasing structural disorder In addition to depolymerization, the irradiated CJ1 glass structure is also subjected to progressively increasing disorder. This generally results in broader curves and lower peak
Fig. 2. Number of 3-coordinate boron atoms (a), nonbridging oxygen atoms (b) and network-modifying sodium atoms (c) versus deposited energy in CJ1 glass subjected to 95 4 keV cascades.
2710
G. Bureau et al. / Nucl. Instr. and Meth. in Phys. Res. B 266 (2008) 2707–2710
values for the spatial, angular, and ring size distributions, reflecting the increasing structural disorder as the 4 keV cascades accumulate. These various distributions gradually leveled off toward a saturation plateau. 4. Discussion The volume expansion obtained using the model (Fig. 1) was greater by a factor of about 2.5 than for the experimental data. This discrepancy could be attributable to atoms ejected from the cascade core toward the periphery, resulting in a lower atomic density at the core of the cascade [13]. Annealing by atomic diffusion of low-density zones would require computing time beyond the limits of molecular dynamics. This type of transient phenomenon could account for the overestimated simulated volume expansion. Volume expansion follows an exponential curve that was fitted using Marples’ model [13]. Volume expansion in this model is proportional to the damaged volume fraction, based on the hypothesis that the passage of a projectile is sufficient to fully damage the irradiated zone. Developing the model gives: Dq ¼ Að1 eq0 VD Þ q0
ð2Þ
where q0 is the glass density, Dq the glass density variation, A the maximum density variation (%), D the cumulative alpha decay (projectiles/g) and V the glass volume damaged by one projectile. The density variation at saturation, A, is 9% ± 1% and volume damaged per projectile, V, is about 8.4 nm3 ± 1.5 nm3. We attempted to interpret the volume parameter of Marples’ model by calculating the simulated cascade volumes. The morphology of a displacement cascade is characterized by a core, in which very high-energy events occur (head-on collisions between atoms, atoms stripped from their equilibrium sites, individual displacements, strong thermal agitation), and a periphery forming an interface between the irradiated zone and the remaining volume [6]. The peripheral zone is characterized by a cloud of displaced Na atoms accompanying lower-energy events such as collective displacements of polymerized chains. We estimated the mean volumes of the core of the 4 keV cascades and of the core + periphery by a method based on drawing a large number of random points [14]. The proximity of displaced atoms with respect to a random point indicates whether it lies within the cascade or outside it. The cascade volume was then calculated from the ratio between the number of points inside the cascade and number of external points. The total number of Si, B and O displacements was used to estimate the volume of the cascade core, and the total number of all atom displacements (Si, O, B, Na) to estimate the volume of the core + periphery.
The mean cascade core volume was estimated in this way at 11.8 nm3 ± 4 nm3 and that of the core + periphery at 18.7 nm3 ± 7 nm3. The volume parameter from Marples’ model, nearer the cascade core volume than that of the core + periphery, suggests that the mechanisms responsible for volume expansion are contained in the cascade core. The cascade core is the seat of the highest-energy events, and notably the thermal quenching already cited to account for irradiation effects [3], which requires a temporary ‘‘molten” state and which can occur only within the cascade core. Similarly, the ejection of atoms from the cascade core to the periphery also creates a low-density zone in the core. The volume parameter from Marples’ model thus suggests that the causes of volume expansion must be sought partly in the thermal quenching effect – which is itself responsible for the depolymerization and increasing disorder of the glass – and partly in the mechanism by which a low-density zone is created. 5. Conclusion The accumulation of cascades modeled by molecular dynamics in a sodium borosilicate glass allowed us to simulate the evolution of various macroscopic and structural properties tending towards a stabilization plateau for the highest deposited nuclear energy doses. Marples’ model was used to fit the glass volume expansion to the deposited energy dose, giving the damaged volume per projectile. The volume parameter from this model approximates the cascade core volume, suggesting that the underlying mechanisms of volume expansion are contained in the cascade core and are thus related to the highest-energy events: atom ejection and thermal quenching. References [1] S. Peuget, J.-N. Cachia, C. Jegou, X. Deschanels, D. Roudil, V. Broudic, J.-M. Delaye, J.-M. Bart, J. Nucl. Mater. 354 (2006) 1. [2] J.M. Delaye, D. Ghaleb, Nucl. Instr. and Meth. B 191 (2002) 10. [3] G. Bureau, J.-M. Delaye, S. Peuget, G. Calas, N. Deladerrie`re, in: Proceedings of Mate´riaux 2006, Dijon, France, 13–17 October, 2006. [4] S. Peuget, P.-Y. Noel, J.-L. Loubet, S. Pavan, P. Nivet, A. Chenet, Nucl. Instr. and Meth. B 246 (2006) 379. [5] D. Petit-Maire, PhD thesis, University of Paris VI, 1998. [6] J.M. Delaye, D. Ghaleb, Phys. Rev. B 61 (2000) 14481. [7] J.M. Delaye, D. Ghaleb, J. Non-Cryst. Solids 195 (1996) 239. [8] B.P. Feuston, S.H. Garofalini, J. Chem. Phys. 89 (1988) 5818. [9] J.-M. Delaye, V. Louis-Achille, D. Ghaleb, J. Non-Cryst. Solids 210 (1997) 232. [10] F.H. Stillinger, T.A. Weber, Phys. Rev. B 31 (1985) 5262. [11] J.F. Ziegler, J.P. Biersack, U. Littmark, The Stopping Range of Ions in Matter, Pergamon, New York, 1985. [12] G.J. Martyna, M.E. Tuckerman, J.T. Douglas, M.L. Klein, Mol. Phys. 87 (1996) 1117. [13] J.-M. Delaye, D. Ghaleb, J. Nucl. Mater. 348 (2006) 243. [14] J.-M. Delaye, Private communication, 2006.