- Email: [email protected]
Molecular dynamics simulation of SiO2 + B2O3 + Na2O + ZrO2 glass
IOURNAL OF
ELSEVIER
Journal of Non-Crystalline Solids 195 (1996) 239-248
Molecular dynamics simulation of SiO 2 + B 2 0 3 4- Na20 + ZrO 2 glass J.-M. Delaye
D. Ghaleb b
a, *
a CEA / DTA / DECM / SRMP: Commissariat ~ I'Energie Atomique, Direction des Technologies Avanc~es, D~partement d'Etude du Comportement des Mat~riaux, Section de Recherches de M~tallurgie Physique, Gif sur Yvette cedex, France b CEA / DCC / DRDD / SCD: Commissariat ~ l'Energie Atomique, Direction du Cycle du Combustible, D~partement du Retraitement, des D~chets et du D~mant~lement, Service de Confinement des D~chets, Gif sur Yvette cedex, France Received 30 March 1995; revised 11 July 1995
Abstract
A complex oxide glass involving four different oxide species was modeled using pair potentials completed by some three-body terms. The numerical model developed combines the major nuclear waste glass oxides (SiO 2, B203, Na20, ZrO 2) and maintains the same percentages of network formers and modifiers as in the actual glass. An overall correct structure was obtained. The model presents satisfactory agreement with experimental density, thermal expansion coefficient, diffusion of Na and a correct order of magnitude for viscosity.
1. Introduction
Oxide glasses intended for disposal of French nuclear wastes are characterized by a complex composition [1,2]. Silica represents about 45% of the total oxide weight. The other main oxides a r e B 2 0 3 ( ~ 15 wt%), N a 2 0 ( ~ 10 wt%), A1203 ( ~ 5 wt%) and CaO ( ~ 5 wt%). The remaining components each contribute less than 3% of the total. These glasses are produced for containment of radioactive nuclear wastes, such as the actinide elements, and it is, therefore, crucial to assess glass behavior under irradiation.
* Corresponding author. Tel: + 33 69 08 66 63. Telefax: + 33 69 08 68 67. E-mail: [email protected].
Molecular dynamics are a very useful technique for investigating phenomena on an atomic scale. Having fitted the potential to reproduce globally correct features of the oxide glass, we will study different phenomena related to displacement cascades in these systems. We present here values of the potential parameters fitted to the experimental structural features of a glass of the same composition as the simulated one. When speaking of the structure of an oxide glass, the network formers are generally dissociated from the network modifiers. Among the network formers, B, Si, AI... are responsible for the network skeleton. Some network modifiers, such as Na or Li, cause the network to become fragile as a result of depolymerization. The numerical model developed includes the ma-
0022-3093/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0 0 2 2 - 3 0 9 3 ( 9 5 ) 0 0 5 2 7 - 7
240
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996)239-248
jor oxides of nuclear waste glass, i.e., SiO2, B203, Na20 and ZrO 2, with a composition such that the percentages of network formers and modifiers correspond to those of the actual glass. Zirconium oxide was added in the model to allow for the influence of a metallic component on the oxide glass structure. The cohesion of this type of material is ensured by semi-covalent bonds (Si-O, B - O ) and a strong ionic bond (Na-O). The previous form of interaction potentials, based on pair potentials proposed by Born, Mayer and Huggins, generally gives correct structural features for oxide glasses [3-6]. In that approach, the Coulomb charges of the atoms are constant and cannot vary during the simulations. Adjustable parameters are chosen to reproduce the first-neighbor distances around the ions. However, simply applying pair potentials is not sufficient. By modeling SiO2, B203 and alkaline oxide mixtures, Soules and Varshneya [7] found an excessive coordination number around Si atoms. This implies that some Si atoms are not tetra-coordinated, in disagreement with experimental observations which show that all Si atoms are surrounded by a tetrahedron of O atoms. Three-body terms were, therefore, added to the simple two-body terms, to impose a precise value on the local angles. This three-body potential is founded on the covalent nature of some bonds, especially the S i - O and B - O bonds, the particular local order of which cannot be correctly reproduced by a simple pair potential. Different forms of three-body potential have been used by Feuston and Garofalini [9], Stillinger and Weber [10], Vashishta et al. [11] and Anderson et al. [12]. In these studies, the structural features of the computed models, first-neighbor distances and coordination numbers are in suitable agreement with experimental findings. Recently, new algorithms have been developed to address the problem of covalent bonds, with special attention to the Coulomb charges of the atoms entering the bonds. Alavi et al. [13] developed a molecular dynamics algorithm to deal with the problems of charge transfer in a SiO 2 + Li20 system, and they successfully distinguish between bonding and nonbonding O from the value of the equilibrium charge. Simple pair potentials can correctly reproduce the local tetrahedral order of Si atoms when the charges are given suitable attention. Balasubramanian and
Rao [14] differentiate the charges on the O atoms according to their local order (bridging or non-bridging between two Si atoms). Successful structural features emerge from their study of 2SiO 2 + K 2 0 glass. This approach is somewhat limited because the charge of a particular atom cannot evolve during the run, and no structural changes are allowed. Some experimental results are available to test the validity of the structure presented in this paper and to re-adjust some potential parameters if necessary. The most interesting for us is the work by Petit-Maire which consists of an EXAFS study of oxide glasses of the same composition as the modeled glass [15]. He investigates local order around the Zr atoms, finding ,an average first-neighbor Z r - O distance equal to 2.08 A ( + 0.04) and an average Zr atom coordination of 6.3. To complete the results of Petit-Maire, some experimental work has been done in the SCD at Marcoule to determine the density and viscosity of glass with the same composition as the modeled glass.
2. Experimental and computational procedure
2.1. Experimental procedure The experimental glass with the same composition as the model is a laboratory glass specimen. The glass components (Table I) were melted at 1200°C in a platinum crucible. The glass density was measured at room temperature on the sample in air and in pure water using a Regnault pycnometer. The average value of three different measurements for the four-oxide glass was 2.465 g cm -3. A Haake double-rotating cylinder viscometer was used in which molten glass was placed between the cylinders and sheared when the internal cylinder rotated. The measurements obtained between 1000 and 1300°C were used to define the Vogel-Fulcher-Tammann ( V - F - T ) parameters of the glass (Table 2).
Table 1 Experimental glass composition(mol%) SiO2
B203
NazO
ZrO2
67.247
17.411
13.548
1.793
J.-M. Delaye, D. Ghaleb /Journal of Non-Crystalline Solids 195 (1996)239-248 Table 2 Experimental viscosity measurements for the four-oxide glass T (°C)
7/ (Pa s)
V - F - T parameters Ln ~7= A + B / ( T - To)
1304 1255 1184 1112 1059
9.88 15.3 31.48 77.36 167.86
A = - 2.261 B = 2907.59 To= 411°C
2.2. Calculation procedure
ni
where zi and zj are the ion charges, o-i and ~rj are their atomic radii and n i and nj are the number of valence shell electrons. The first term represents a repulsive interaction between two ions and the second term is the Coulomb interaction; b and rij are adjustable parameters. The Coulomb terms were calculated using the complete Ewald sum method [17]. The summation was truncated at 13 ,~ in real space and at Nmax = 6 in reciprocal space. Only reciprocal vectors, 2 n : r / L , 2 n 2 1 r / L and 2n3Ti'//L , with n~ +n22 + n 2 < N 2 were considered (L is the length of the simulation cell). The cutoff radii were adjusted to that the neglected forces were on the same order of magnitude. We found the same values as those determined by Kieffer and Angell [18], whose adjustment was based on potential energy considerations. A shorter cutoff radius (6.5 ,~) was used for the purely shortrange repulsive term. The parameter values depended on whether or not a three-body term was applied, so different sets of adjustable parameters were employed. To constrain some local angles ( O - S i - O , S i - O Si and O - B - O ) which take particular values due to the partly covalent nature of the interactions, we applied to a triplet of atoms j - i - k a potential energy of the form --
The simulation technique involved molecular dynamics applied to oxide glasses contained in a cubic box. The total number of ions in a basic cell was equal to 1536 atoms. The details of the composition are indicated in Table 3. The four-oxide glasses roughly approximated the basic matrix of the French light water reactor waste containment glass. Proportions between the component oxides were conserved. Classical periodic conditions allowed modeling of pseudo-infinite systems. The motion equations were solved using Verlet's algorithm. In each simulation, the unit time step was 1 fs. To calculate the equilibrium volume of a system at a given temperature, we used the algorithm developed by Andersen [16], which introduces two new variables in the Hamiltonian equations, the volume and a fictive ' m a s s ' acting as an inert mass when the volume fluctuates. In each simulation, this mass was assumed equal to 1.54 × 105 g cm -3. The external pressure imposed was atmospheric pressure. To simulate the atomic interactions, we applied the B o r n - M a y e r - H u g g i n s potential. A potential energy, q~(ru), was associated with each pair of ions separated by a distance rij: tl)( rij ) = b 1 + - - +
241
exp nj ]
k
Pij
1 zi zj + - --, 47r~ 0 rij
(1)
Table 3 Computed glass composition (mol%) SiO2
B20 3
Na20
ZrO 2
67.25
17.68
13.32
1.74
q)(ru,rik,O~ik) = Aiexp - Fij -- rci
m a x ,
+ - rik -- rci
x (cos o . - cos Oo)
(2)
This type of potential has been previously developed by Feuston and Garofalini [9]. 00 represents the reference angle towards which the triplet angle is constrained. Ai, Yi and rci are adjustable parameters depending on the nature of atom i. The adjustable parameter values are indicated in Table 4. It is noteworthy that the constraint imposed on S i - O - S i triplets is weak, and the constraint imposed on the O - B - O triplets is strong. We performed three different experiments, one with no three-body terms (system 1), another with a three-body term only for the O - S i - O and S i - O - S i angles (system 2) and the last by imposing three-body terms for O - S i - O , S i - O - S i and O - B - O angles (system 3). For each experiment, the parameter sets associated with the pair potentials had to be adjusted
242
J.-M. Delaye,D. Ghaleb/ Journal of Non-Crystalline Solids 195 (1996)239-248
Table 4 Parameters used for the three-body terms Parameter
O-Si-O
Si-O-Si
O-B-O
A (10- t i erg): y (,~): r c (/~): 0o system 2 (deg): 0o system 3 (deg):
24 2.6 3 109.47 109.47
1 2.0 2.6 109.47 160.0
1920 2.27 2.1 109.47 109.47
was determined using A n d e r s e n ' s algorithm. The runs typically lasted 5000 time steps of 10 -15 s. Knowing the equilibrium volumes, the constant volume algorithm could be used again to relax the systems. Relaxations typically for 10000 time steps were performed for the structural investigations.
2.4. Method of determination of some physical quantities
because the first-neighbor distances were effected by the application of the three-body terms. Table 5 shows the sets of adjustable parameters used for the pair potential terms, according to the nature o f the three-body terms applied.
For each configuration, some physical quantities such as the density, Na diffusion coefficient and viscosity were measured. The diffusion coefficient of Na, the most highly mobile species, was measured at a high temperature ( T = 1073 K), using Einstein's relation:
2.3. Modeling of sample preparation
DNa = ( r 2 ) / 6 t ,
The initial system was a cristobalite structure in which B 2 0 3, N a 2 0 and ZrO 2 groups were introduced on random sites. The system temperature then increased quickly due to structural instabilities. The cristobalite structure contains some linear S i - O - S i bonds, unstable with regard to the B o r n - M a y e r Huggins potential. The temperature was controlled by controlling the atomic velocities, decreasing to 5500 K when it exceeded 6000 K. W h e n the temperature stabilized between 5500 and 6000 K, the liquid was relaxed for 6000 time steps of 10-15 s before quenching to room temperature at a rate o f 5 X 1014 K s-~. A series of structural configurations was stored at intermediate temperatures. A constant volume algorithm was used at each stage of the preparation with a density slightly above the equilibrium density to avoid void formation. The equilibrium volume of each configuration
where ( r 2) is the mean square displacement of the Na atoms and t is the duration of the simulation. The viscosities were determined by a method previously used by Gosling on a simple LennardJones fluid [19], in which an external perturbation was applied in the form o f a sinusoidal profile of additive forces. In response to this perturbation, atomic fluxes developed in the system, characterized by a sinusoidal profile of velocities, in phase with the force profile. The viscosity was then calculated by the formula
Table 5 List of parameters used for pair potentials Parameter Si O B q: +4 -2 +3 n: 8 8 2 o- (A) (system I): 1.148 1.42 0.74 o" (,~)(system 2): 1.1008 1.42 0.74 o- (A) (system 3): 1.1008 1.42 0.72
Na +1 8 1.17 1.17 1.17
Zr +4 8 1.45 1.4 1.425
For systems I and 2: r o = 0.29 ,~, b = 0.221 eV. For system 3: rij = 0.29 ,~ except ro_ o = 0.35 ,g,, b = 0.221 eV.
rl = ( p l k Z ) ( F / v ) ,
(3)
(4)
where F and v are the amplitudes of the force and velocity profiles, respectively, p is the atomic density of the system and k = 2 w / L is the wave vector of the perturbation. In the following, F is considered equal to ! x 10 -4 in cgs. Strictly, this last equation must be applied to the limiting case of a null perturbation and a null wave vector. For the viscosity determinations, the constant volume algorithm was always used.
3. Results
3.1. Structural properties The results are distinguished according to whether or not three-body potentials were applied around the
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 11996) 239-248
243
Table 6 Structural characteristics of system 1 (experimental data in parentheses) Property
Si-O
B-O
Na-O
Zr-O
Ncoord:
4.26 (4.0)
3.37 (3.67)
5.86
6.38 (6.3)
Rcut (,~):
2.1
1.9
3.1
3.0
_
1.40 (1.40)
2.48
2.17 (2.08)
'~
dlst neigh (,~k): 1.67 (1.60)
Si or B atoms. The first set of results concerns the structural characterization of a system when only pair potentials are used. Table 6 indicates the firstneighbor distances and coordination numbers around the cations, together with the experimental observations. Concerning first-neighbor distances, note that the S i - O and Z r - O distances are slightly too long. Moreover, the coordination number is too low around the B atoms and too high around the Si atoms. Every Si atom in the actual glass is tetra-coordinated, while, in the model, an excessive percentage of Si atoms have a coordination number equal to 5. The Zr coordination number is consistent with the experimental value. Therefore, we constructed a second set of models, adding to the pair potential a three-body term around the Si atoms. We constrained the O - S i - O angles to the value of the perfect tetrahedrai angle, i.e., 109.47 °, and a weak constraint was applied to the S i - O - S i angles, again with the most energetically favorable angle equal to 109.47 ° . For this series of models, we slightly decreased the atomic radii of Si and Zr to correct the deviations from experimental results noted above. The new structural characterizations are given in Table 7. Re-adjusting the atomic radius of Si corrected the excessive first-neighbor distance obtained with the previous set of parameters. Moreover, because of the angular constraint, every Si atom be-
Table 7 Structural characteristics of system 2 Property
Si-O
B-O
Na-O
Zr-O
Ncoord: Rcut (A):
4.04 2.1
3.48 1.9
5.47 3.1
5.88 3.0
dlst neigh ( ~ t ) :
1.60
1.39
2.49
2.12
~
z~
.... 8~1
i .... 911
m
, .... 104)
, .... I iiJ
i .... 1211
i .... 1311
i .... 140
I ,, 1511
Angle (deg)
Fig. 1. Angular distribution of atomic triplets.
came tetra-coordinated. The first-neighbor distances agreed well with experimental values. The coordination number around the B atoms was still too low, with 48% tri-coordinated B and 52% tetra-coordinated B. Fig. 1 shows the angular distributions of the O-Si-O, O-B3-O and O - B 4 - O angles, in which B 3 and B 4 are, respectively, the tri- and tetra-coordinated boron atoms. The mean O - S i - O angle was very near the perfect tetrahedral angle; the mean O-B3-O angle was near 120 °, signifying that the B 3 atoms were surrounded by a triangle of O atoms, and the m e a n O - B 4 - O angle was near the perfect tetrahedral angle. B 4 atoms were surrounded by tetrahedrons of O atoms. The mean S i - O - S i angle was equal to 158 °. The Zr coordination decreased compared with the previous value. The last stage was to correct the coordination number around the B atoms, which was just below 3.5 rather than the expected value of 3.67 [20]. This correction was obtained by applying a three-body term on the O - B - O angles. We noted that a relatively large constraint was required to increase the B coordination number, compared with what was needed to adjust the Si coordination number to 4. This is because the radius of the B coordination sphere is shorter and the mean coordination number due to geometrical packing of atoms is lower than the physical number (contrary to Si, for which a three-body term was applied to decrease the coordination number). It is more difficult to force a coordination number larger than the natural packing value. We also changed the reference angle for the weak
244
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996) 239-248
Table 8 Stluctural characteristics of system 3
3O
Property
Si-O
B-O
Na-O
Zr-O
Ncoord: Rcut (,~):
4.0 2.1 1.61
3.67 1.9 1.42
5.71 3.1 2.5
5.62 3.0 2.14
dlst neigh (m):
•
,
,
-
,
"
,
"
"
"
,
,
"
"
211 I B-O
svO
constraint imposed on the S i - O - S i angles from 109.47 ° to 160° , without any significant effect. In addition to applying a three-body term around B, we increased the ro_ o parameter associated with o x y g e n - o x y g e n repulsion, diminishing the repulsion strength between two O atoms. This made it easier to increase the coordination number around B, but modified some of the first-neighbor distances and coordination numbers. Other attempts were, therefore, made to readjust the atomic radii. The structural characteristics of the new model are given in Table 8. The first-neighbor distances are in good agreement with the experimental distances, except for a slightly overestimated Z r - O distance. The adjustable parameters assigned to the Zr species represent a trade-off between the coordination number and the first-neighbor distance. Increasing either parameter decreases the other. The final result is characterized by a low coordination number (our models contain few Zr atoms, and this value is probably subject to rather large fluctuation errors) and a large Z r - O distance. The coordination numbers around Si and B are well adjusted. The partial radial distribution functions are plotted in Fig. 2. Positions used for the calculations were averaged at 293 K over 300 time steps. No modifications of the angular distributions were observed for O - S i - O and S i - O - S i , and the mean S i - O - S i angle was 157 °, close to the previous value despite the change in the constraint angle from 109.47 ° to 160 °. In the literature, three series of experiments determined two different mean S i - O - S i angles. The first is due to Mozzi and Warren [21], who found a mean S i - O - S i angle of 144° in silica; Greaves et al. [22] reported a mean S i - O - S i angle equal to 160 °, while a recent NMR study of potassium disilicate indicated a mean S i - O - S i angle of 143 ° [23]. Our result is close to that of Greaves et al. The O - B 3 - O angles were modified by applying
.5
. . . .
i
o
. . . .
-
I
.
.
.
1
.
.
.
,
l
. . . .
,
1
R
. . . .
,
5
. . . .
6
{A}
Fig. 2. Partial radial distributionfunctions between cations and O atoms. The curve was calculated at room temperature by averaging the atomic positionsover 300 time steps.
the three-body term around the B atoms. We found a mean O - B 3 - O angle of 115 °, less than the 120 ° noted above. The triangular configurations of O atoms around the B 3 atoms are, therefore, distorted by the three-body term. 3.2. P h y s i c a l p r o p e r t i e s
Andersen's algorithm allows the equilibrium volume of a system to be measured at different temperatures. Fig. 3 shows the density versus temperature curves obtained for systems 2 and 3. An interesting although logical consequence of applying a constraint on the O - B - O angles is the densification of
O
O O
O
X45
El C]
[] []
[] t-'l
,
,
i
400
,
,
.
i
,
,
,
i
800
12oo
Temperature
(K)
,
, E I
, 16oo
2800
Fig. 3. Density versus temperature for systems 2 (O) and 3 ([3). • represents the experimentaldensity of the four-oxideglass.
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996) 239-248
system 3 as the thermal expansion coefficient increases. Some experimental results provided by the SCD laboratory provide a basis for estimating the validity of the numerical models. At room temperature, the experimental density of a glass of the same composition as the model is 2.465 g cm -3, not far from the computed value of 2.508 g cm -3. The density is still close to the experimental value when a three-body term around B is applied. No thermal expansion coefficient is available for our simplified system, so we compared the simulation results with the thermal expansion coefficient of the complete nuclear glass. For the French R7T7 and German GP 98-12 nuclear glasses, the thermal expansion coefficients are, respectively, 8.3 × 10 -6 and 8.1 × 10 -6 K - ~; the thermal expansion coefficients of the simulated systems without and with a three-body term around B are 2.41 × 10 -5 and 8.11 X 10 -6 K -1, respectively. Here again, the application of a threebody term around the B atoms provides a better approximation of the experimental thermal expansion coefficient. For system 3, we measured the vitreous transition temperature, Tg, by plotting the evolution of the atomic potential energy versus the temperature; the results are shown in Fig. 4. The intersection of the linear extrapolations at low and high temperatures corresponds to a Ts equal to roughly 1000°C, probably higher than the experimental results. The shift of Tg towards higher temperatures is a classic simulation artifact, for which very fast quench rates, com-
-5.98
10 '
.5.99
IO i
4.#10
D) I
4~AIJ
I0 t
T:
I
I
I
T
[
I
I
•
i
~
Tempcral~r¢
I
I
I
I
I
?
I
I
I
(h.~
Fig. 4. Atomic potential energy versus temperature for system 3 (Tg determined by the intersection of the two extrapolated linear
fits).
245
1114
8 •
•
0
•
0
I0 ~
o 5.5
6
6.5
7
7.$
8
8,5
9
9.5
lfr (*10') Fig. 5. Logarithm of Na atom diffusion coefficient versus reciprocal of T: 0 , points measured in system 2; A, points measured in system 3; O , experimental points measured on complete nuclear
glasses [28].
pared with the experimental rates, are imposed to hasten the creation of low-temperature structures. The Tg of the French nuclear glass is 502°C; the German nuclear glass, which contains more network modifiers, has a Tg of 543°C. We also examined the Na diffusion coefficients, DNa, at high temperatures. The Na atoms are small, and their Coulomb charge is low; the interactions binding them to the network are weak, and they can move easily. DNa is the most easily measurable diffusion coefficient. An Arrhenius diagram of DNa for systems 2 and 3 versus the reciprocal of the temperature, T, is plotted on a logarithmic scale in Fig. 5. The values measured with the German nuclear glass are shown for comparison. To obtain equilibrium values of DNa at 1073 K, longer runs are necessary to allow proper initial relaxation of the systems. Several observations may be made concerning this figure. First, the three slopes are basically similar, indicating good agreement for the activation energy between the computed and physical systems. A linear fit of the data provides an activation energy of 0.92 eV in system 2 and 1.03 eV in system 3. The corresponding activation energy in the actual nuclear glass is 1.08 eV. Second, the application of a threebody term around B significantly lowers DNa tO below experimental values. This qualitative behavior (computed results lower than experimental results) is
246
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996) 239-248
not surprising, since the complete nuclear glasses contain a larger fraction of network modifiers, and so is probably more disordered, insofar as we can speak of a level of disorder for glass. The diffusion of Na atoms is facilitated by the increased network disorder. The last physical value calculated was viscosity, using the method described above. The establishment of a dynamic flux in response to the application of the perturbation is a very long process, as the reader can see from Fig. 6. The y-axis represents the amplitude of the sinusoidal velocity profile versus the time of simulation. The velocity profile is always determined by dividing the difference between the current and initial positions by the time elapsed between the two states. The current and initial positions are averaged over 100 time steps to get rid of part of the thermal motion. After a sudden increase due to the application of the perturbation, the signal progressively decreases and continues to decrease beyond the 38 000 time steps of the simulation. The amplitude representing the equilibrium dynamic profile corresponds to the asymptotic value of the curve. In fact, two different phenomena interfere in the process, the first being the elastic deformation of the system and the second being the progressive establishment of viscous flow. Fig. 6 shows that the viscous flow was hidden by elastic deformation throughout the simulation. One means of evaluating
30011
~. m
2q)O0
I$1HI
I000
.
.
.
.
.
5000
.
.
.
.
.
.
.
IOIlO~
.
.
.
.
.
IGOI)O
.
.
.
.
.
.
.
201HIO
.
.
.
.
250(10
.
.
.
.
.
3000{I
.
i
.35000
.
.
.
.
40OOO
T i m e o f ~imulation (f;)
Fig. 6. Velocity profile amplitude versus simulation time (perturbation applied at t = 0, T = 1923 K).
tm
c
20
0
*
I
,
I loot)
,
i
*
I
,
i
,
MHH)
Time of the simulation
n 600O
,
i
,
I ~NH)
,
i
, I001~)
(fs)
Fig. 7. Velocity profile amplitude versus simulation time after release of the perturbation at t = 0 (T = 1923 K).
the amplitude of the velocity profile relevant to the viscous flow is to determine the asymptotic value of the curve by a suitable adjustment on the simulated points. Unfortunately, however, very different values are obtained when different fitting curves are used (1/t or e x p ( - a t ) in particular), while there is very little difference in the quality of the fit. Prolonging the application of the perturbation would be too time-consuming for a hypothetical stabilization of the velocity profile. Therefore, we attempted to obtain the viscosity by stopping the perturbation at a time, t, equal to 25 000 time steps and following the evolution of the velocity profile amplitude, calculated as the difference between the atomic positions of the henceforth undisturbed system and the initial positions, divided by the time during which the perturbation is applied. This new graph is shown in Fig. 7, where it may be seen that relaxation is much faster than before and the amplitude of the velocity profile reaches an equilibrium value in a rather short time. This second stage of the experiment, the release of the perturbation, furnishes an equilibrium amplitude of the velocity profile which can be taken with confidence for the viscosity calculation. The viscosity of system 3 was 3.3 Pa s at 1823 K and 2.1 Pa s at 1923 K compared with experimental viscosities estimated at 1.93 and 1.2 Pa s, respectively, on systems of the same composition as the model; the experimental values were extrapolated from V - F - T parameters. The order of magnitude of
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996)239-248
the viscosity of the simulated systems is correct, although the viscosity of the models is overestimated.
4. Discussion Computer models were constructed and demonstrated to produce predictions of the structural characteristics and some physical properties. It was necessary to introduce a three-body term applied to the O - B - O angles to increase the coordination number of the B atoms. The consequence of this modification was the deformation of the triangles of O atoms surrounding the tri-coordinated B atoms. We tried to apply a three-body potential with two local minima, at 109.47 and 120.0 °, to the O - B - O triplets to force either local triangular order or local tetrahedral order. Unfortunately, this method failed because the two minima were too close to one another, compared with the width of the O - B - O angle distribution. These two minima of the potential did not result in production of two separated O - B - O angle distributions. A secondary effect of the constraint applied on the O - B - O angles is lower Na diffusion coefficients; this lowering is a logical result, since the Na atoms travel via interstitial sites surrounded by O atoms and the added constraint decreases the mobility of the O atoms, making the Na paths more rigid and difficult to cross. System densification following the application of the three-body O - B - O term is a second possible explanation for the lower Na diffusion rate. As noted above, the computed glass transition point, Tg, was higher than the experimental temperature for an actual nuclear glass. This difference probably proceeds from the very fast quenching rate applied to terminate the computed systems. Scholze [24] proposed as an example the evolution of the for borosilicate glasses, and showed a difference of nearly 100 K between quench rate of 10 - 4 K min-l and 10 K min -~ Extrapolating the experimental curve to 1014 K s -I further increases Tg by 400 K. This difference is of the order of magnitude of that observed between the computed and experimental Nevertheless, the Na diffusion coefficients were
247
lower but not far from the experimental values measured on German nuclear glasses, which are not only less rigid than the model but also contain more than one alkaline species, suggesting that the mixed alkali effect and the influence of the other components are expected to be weak in complex nuclear glasses. The latter conclusion has previously been stated by Soules and Busbey [8] who did not detect a strong mixed alkali effect in simulations of diffusion in Na20 + K 2 0 + SiO 2 systems. The mixed alkali effect is known to decrease when the temperature increases, tending, in particular, to disappear when the diffusion mode becomes liquid-like. We observed that viscosities at 1823 and 1923 K were higher than experimental values (extrapolated from a Vogel-Fulcher-Tammann curve). Some care must be taken with these values because the exact value of the viscosity must be determined within the limit of a null perturbation. However, as the level of the perturbation we applied was approximately 2% of the mean force experienced by the atoms, this value should not be far from the exact one. It is possible to quantify the applied shear stress, s. The external perturbing force ( F ) profile is a sine wave, so the applied shear stress is not constant over the entire sample. Nevertheless, an order of magnitude can be given by considering the total force exerted on the atoms contained in the half simulation cell and experiencing a positive perturbing force, dividing by the cell cross-section: s = 6.7 × 10 9 N m - 2 . The theoretical shear strength of solid oxide glasses is typically about 10 ~° N m -2. Second, due to the low velocity profile imposed, it must be calculated by taking differences between two averaged positions separated by a rather large time step; thus, thermal agitation may have some influence. When an atom moves in a direction orthogonal to the direction of the external perturbing force, the magnitude of the perturbing force experienced by this atom changes, widening the velocity profile, and the profile amplitude is underestimated. We expect this effect to be negligible due to the low diffusion coefficient of the network formers, even at a temperature of 1923 K. The last point we would like to emphasize is the inability of the potential to represent physical characteristics of actual glasses. In particular, the cohesive energy is not well reproduced because, by taking
248
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996) 239-248
formal charges (Si 4+, 0 2-, B 3+, Na +, Z r 4 + ) , w e attribute to the atomic interactions a purely ionic feature, offset by the three-body terms representing the semi-covalent nature of some atomic bonds. The reality is different, and the covalent nature of the atomic bonds is related to the geometrical form of the electron densities of states localized in particular areas on the covalent bonds. The consequence is that the charge location is much more complex than the model used in this simulation and, consequently, the cohesive energy is lower. Several attempts to determine more exact charges have been published. Empirical potentials have been derived from ab initio calculations for water [25] and pure silica [26]. Nevertheless, it is a very challenging task to take into account these physical characteristics on mixed systems because, as has been shown experimentally [27], the charge of a core depends on its local environment. It is, therefore, incorrect to assign fixed charges to O or B atoms in particular because the surroundings of these atoms are different from one site to another.
5. Conclusion An overall correct structure of a complex oxide glass involving four different oxide species was modeled using pair potentials completed by some three-body terms. These structures show good agreement with experimental results for the density, thermal expansion coefficient and Na atom diffusion and are on the same order of magnitude for the viscosity.
Acknowledgements One of the authors (J.-M.D.) is grateful to the CEA ( D C C / D R D D / S C D ) for the financial support of this study. Dr N. van Doan kindly provided the original program for the development of the simulations and, together with Dr Y. Limoge, provided valuable insights after a critical reading of this paper.
References [1] F. Pacaud, C. Fillet and G. Baudin, in: Proe. 2nd Int. Seminar on Radioactive Waste Products, Ji~lich, Germany (1990) p. 577. [2] F. Pacand, C. Fillet and N. Jacquet-Francillon, Scientific Basis for Nuclear Waste Management, Vol. XV, 257, ed. C. Sombret (1992) p. 161. [3] T.F. Soules, J. Chem. Phys. 71 (1979) 4570. [4] T.F. Soules, J. Non-Cryst. Solids 49 (1982) 29. [5] K.V. Damodaran, V.S. Nagarajan and K.J. Rao, J. Non-Cryst. Solids 124 (1990) 233. [6] Q. Xu, K. Kawamura and T. Yokokawa, J. Non-Cryst. Solids 104 (1988) 26l. [7] T.F. Soules and A.K. Varshneya, J. Am. Ceram. Soc. 64 (1981) 145. [8] T.F. Soules and R.F. Busbey, J. Chem. Phys. 75 (1981) 969. [9] B.P. Feuston and S.H. Garofalini, J. Chem. Phys. 89 (1988) 5818. [10] F.H. Stillinger and T.A. Weber, Phys. Rev. B43 (1991) 1194. [I 1] P. Vashishta, R.K. Kalia, J.P. Rino and I. Ebbsj/5, Phys. Rev. B41 (1990) 12197. [12] D.C. Anderson, J. Kieffer and S. Klarsfeld, J. Chem. Phys. 98 (1993) 8978. [13] A. Alavi, L.J. Alvarez, S.R. Elliott and I.R. McDonald, Philos. Mag. B65 (1992) 489. [14] S. Balasubramanian and K.J. Rao, J. Phys. Chem. 98 (1994) 10871. [15] D. Petit-Maire, PhD thesis, University of Paris VI (1988). [16] H.C. Andersen, J. Chem. Phys. 72 (1980) 2384. [17] M.J.L. Sangster and M. Dixon, Adv. Phys. 25 (1976) 247. [18] J. Kieffer and C.A. Angell, J. Chem. Phys. 90 (1989) 4982. [19] E.M. Gosling, I.R. McDonald and K. Singer, Molec. Phys. 26 (1973) 1475. [20] W.J. Dell, P.J. Bray and S.Z. Xiao, J. Non-Cryst. Solids 58 (1983) 1. [21] R.L. Mozzi and B.E. Warren, J. Appl. Crystallogr. 2 (1969) 164. [22] G.N. Greaves, A. Fontaine, P. Lagarde, D. Raoux and S.J. Gurman, Nature (1981) 611. [23] I. Farnan, P.J. Grandinetti, J.H. Bultisberger, J.F. Stebbins, U. Werner, M.A. Eastman and A. Pines, Nature 358 (1992) 31. [24] H. Scholze, Glass: Nature, Structure and Properties (Springier, Berlin, 1991) (see p. 50). [25] G.G. Hall and C.M. Smith, Int. J. Quantum Chem. 42 (1992) 1237. [26] K. de Boer, A.P.J. Jansen and R.A. van Santen, Chem. Phys. Lett. 223 (1994) 46. [27] C.H. Hsieh, H. Jain, A.C. Miller and E.I. Kamitsos, J. Non-Cryst. Solids 168 (1994) 247. [28] Hj. Matzke and E. Vernaz, J. Nucl. Mater. 201 (1993) 295.
ELSEVIER
Journal of Non-Crystalline Solids 195 (1996) 239-248
Molecular dynamics simulation of SiO 2 + B 2 0 3 4- Na20 + ZrO 2 glass J.-M. Delaye
D. Ghaleb b
a, *
a CEA / DTA / DECM / SRMP: Commissariat ~ I'Energie Atomique, Direction des Technologies Avanc~es, D~partement d'Etude du Comportement des Mat~riaux, Section de Recherches de M~tallurgie Physique, Gif sur Yvette cedex, France b CEA / DCC / DRDD / SCD: Commissariat ~ l'Energie Atomique, Direction du Cycle du Combustible, D~partement du Retraitement, des D~chets et du D~mant~lement, Service de Confinement des D~chets, Gif sur Yvette cedex, France Received 30 March 1995; revised 11 July 1995
Abstract
A complex oxide glass involving four different oxide species was modeled using pair potentials completed by some three-body terms. The numerical model developed combines the major nuclear waste glass oxides (SiO 2, B203, Na20, ZrO 2) and maintains the same percentages of network formers and modifiers as in the actual glass. An overall correct structure was obtained. The model presents satisfactory agreement with experimental density, thermal expansion coefficient, diffusion of Na and a correct order of magnitude for viscosity.
1. Introduction
Oxide glasses intended for disposal of French nuclear wastes are characterized by a complex composition [1,2]. Silica represents about 45% of the total oxide weight. The other main oxides a r e B 2 0 3 ( ~ 15 wt%), N a 2 0 ( ~ 10 wt%), A1203 ( ~ 5 wt%) and CaO ( ~ 5 wt%). The remaining components each contribute less than 3% of the total. These glasses are produced for containment of radioactive nuclear wastes, such as the actinide elements, and it is, therefore, crucial to assess glass behavior under irradiation.
* Corresponding author. Tel: + 33 69 08 66 63. Telefax: + 33 69 08 68 67. E-mail: [email protected].
Molecular dynamics are a very useful technique for investigating phenomena on an atomic scale. Having fitted the potential to reproduce globally correct features of the oxide glass, we will study different phenomena related to displacement cascades in these systems. We present here values of the potential parameters fitted to the experimental structural features of a glass of the same composition as the simulated one. When speaking of the structure of an oxide glass, the network formers are generally dissociated from the network modifiers. Among the network formers, B, Si, AI... are responsible for the network skeleton. Some network modifiers, such as Na or Li, cause the network to become fragile as a result of depolymerization. The numerical model developed includes the ma-
0022-3093/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0 0 2 2 - 3 0 9 3 ( 9 5 ) 0 0 5 2 7 - 7
240
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996)239-248
jor oxides of nuclear waste glass, i.e., SiO2, B203, Na20 and ZrO 2, with a composition such that the percentages of network formers and modifiers correspond to those of the actual glass. Zirconium oxide was added in the model to allow for the influence of a metallic component on the oxide glass structure. The cohesion of this type of material is ensured by semi-covalent bonds (Si-O, B - O ) and a strong ionic bond (Na-O). The previous form of interaction potentials, based on pair potentials proposed by Born, Mayer and Huggins, generally gives correct structural features for oxide glasses [3-6]. In that approach, the Coulomb charges of the atoms are constant and cannot vary during the simulations. Adjustable parameters are chosen to reproduce the first-neighbor distances around the ions. However, simply applying pair potentials is not sufficient. By modeling SiO2, B203 and alkaline oxide mixtures, Soules and Varshneya [7] found an excessive coordination number around Si atoms. This implies that some Si atoms are not tetra-coordinated, in disagreement with experimental observations which show that all Si atoms are surrounded by a tetrahedron of O atoms. Three-body terms were, therefore, added to the simple two-body terms, to impose a precise value on the local angles. This three-body potential is founded on the covalent nature of some bonds, especially the S i - O and B - O bonds, the particular local order of which cannot be correctly reproduced by a simple pair potential. Different forms of three-body potential have been used by Feuston and Garofalini [9], Stillinger and Weber [10], Vashishta et al. [11] and Anderson et al. [12]. In these studies, the structural features of the computed models, first-neighbor distances and coordination numbers are in suitable agreement with experimental findings. Recently, new algorithms have been developed to address the problem of covalent bonds, with special attention to the Coulomb charges of the atoms entering the bonds. Alavi et al. [13] developed a molecular dynamics algorithm to deal with the problems of charge transfer in a SiO 2 + Li20 system, and they successfully distinguish between bonding and nonbonding O from the value of the equilibrium charge. Simple pair potentials can correctly reproduce the local tetrahedral order of Si atoms when the charges are given suitable attention. Balasubramanian and
Rao [14] differentiate the charges on the O atoms according to their local order (bridging or non-bridging between two Si atoms). Successful structural features emerge from their study of 2SiO 2 + K 2 0 glass. This approach is somewhat limited because the charge of a particular atom cannot evolve during the run, and no structural changes are allowed. Some experimental results are available to test the validity of the structure presented in this paper and to re-adjust some potential parameters if necessary. The most interesting for us is the work by Petit-Maire which consists of an EXAFS study of oxide glasses of the same composition as the modeled glass [15]. He investigates local order around the Zr atoms, finding ,an average first-neighbor Z r - O distance equal to 2.08 A ( + 0.04) and an average Zr atom coordination of 6.3. To complete the results of Petit-Maire, some experimental work has been done in the SCD at Marcoule to determine the density and viscosity of glass with the same composition as the modeled glass.
2. Experimental and computational procedure
2.1. Experimental procedure The experimental glass with the same composition as the model is a laboratory glass specimen. The glass components (Table I) were melted at 1200°C in a platinum crucible. The glass density was measured at room temperature on the sample in air and in pure water using a Regnault pycnometer. The average value of three different measurements for the four-oxide glass was 2.465 g cm -3. A Haake double-rotating cylinder viscometer was used in which molten glass was placed between the cylinders and sheared when the internal cylinder rotated. The measurements obtained between 1000 and 1300°C were used to define the Vogel-Fulcher-Tammann ( V - F - T ) parameters of the glass (Table 2).
Table 1 Experimental glass composition(mol%) SiO2
B203
NazO
ZrO2
67.247
17.411
13.548
1.793
J.-M. Delaye, D. Ghaleb /Journal of Non-Crystalline Solids 195 (1996)239-248 Table 2 Experimental viscosity measurements for the four-oxide glass T (°C)
7/ (Pa s)
V - F - T parameters Ln ~7= A + B / ( T - To)
1304 1255 1184 1112 1059
9.88 15.3 31.48 77.36 167.86
A = - 2.261 B = 2907.59 To= 411°C
2.2. Calculation procedure
ni
where zi and zj are the ion charges, o-i and ~rj are their atomic radii and n i and nj are the number of valence shell electrons. The first term represents a repulsive interaction between two ions and the second term is the Coulomb interaction; b and rij are adjustable parameters. The Coulomb terms were calculated using the complete Ewald sum method [17]. The summation was truncated at 13 ,~ in real space and at Nmax = 6 in reciprocal space. Only reciprocal vectors, 2 n : r / L , 2 n 2 1 r / L and 2n3Ti'//L , with n~ +n22 + n 2 < N 2 were considered (L is the length of the simulation cell). The cutoff radii were adjusted to that the neglected forces were on the same order of magnitude. We found the same values as those determined by Kieffer and Angell [18], whose adjustment was based on potential energy considerations. A shorter cutoff radius (6.5 ,~) was used for the purely shortrange repulsive term. The parameter values depended on whether or not a three-body term was applied, so different sets of adjustable parameters were employed. To constrain some local angles ( O - S i - O , S i - O Si and O - B - O ) which take particular values due to the partly covalent nature of the interactions, we applied to a triplet of atoms j - i - k a potential energy of the form --
The simulation technique involved molecular dynamics applied to oxide glasses contained in a cubic box. The total number of ions in a basic cell was equal to 1536 atoms. The details of the composition are indicated in Table 3. The four-oxide glasses roughly approximated the basic matrix of the French light water reactor waste containment glass. Proportions between the component oxides were conserved. Classical periodic conditions allowed modeling of pseudo-infinite systems. The motion equations were solved using Verlet's algorithm. In each simulation, the unit time step was 1 fs. To calculate the equilibrium volume of a system at a given temperature, we used the algorithm developed by Andersen [16], which introduces two new variables in the Hamiltonian equations, the volume and a fictive ' m a s s ' acting as an inert mass when the volume fluctuates. In each simulation, this mass was assumed equal to 1.54 × 105 g cm -3. The external pressure imposed was atmospheric pressure. To simulate the atomic interactions, we applied the B o r n - M a y e r - H u g g i n s potential. A potential energy, q~(ru), was associated with each pair of ions separated by a distance rij: tl)( rij ) = b 1 + - - +
241
exp nj ]
k
Pij
1 zi zj + - --, 47r~ 0 rij
(1)
Table 3 Computed glass composition (mol%) SiO2
B20 3
Na20
ZrO 2
67.25
17.68
13.32
1.74
q)(ru,rik,O~ik) = Aiexp - Fij -- rci
m a x ,
+ - rik -- rci
x (cos o . - cos Oo)
(2)
This type of potential has been previously developed by Feuston and Garofalini [9]. 00 represents the reference angle towards which the triplet angle is constrained. Ai, Yi and rci are adjustable parameters depending on the nature of atom i. The adjustable parameter values are indicated in Table 4. It is noteworthy that the constraint imposed on S i - O - S i triplets is weak, and the constraint imposed on the O - B - O triplets is strong. We performed three different experiments, one with no three-body terms (system 1), another with a three-body term only for the O - S i - O and S i - O - S i angles (system 2) and the last by imposing three-body terms for O - S i - O , S i - O - S i and O - B - O angles (system 3). For each experiment, the parameter sets associated with the pair potentials had to be adjusted
242
J.-M. Delaye,D. Ghaleb/ Journal of Non-Crystalline Solids 195 (1996)239-248
Table 4 Parameters used for the three-body terms Parameter
O-Si-O
Si-O-Si
O-B-O
A (10- t i erg): y (,~): r c (/~): 0o system 2 (deg): 0o system 3 (deg):
24 2.6 3 109.47 109.47
1 2.0 2.6 109.47 160.0
1920 2.27 2.1 109.47 109.47
was determined using A n d e r s e n ' s algorithm. The runs typically lasted 5000 time steps of 10 -15 s. Knowing the equilibrium volumes, the constant volume algorithm could be used again to relax the systems. Relaxations typically for 10000 time steps were performed for the structural investigations.
2.4. Method of determination of some physical quantities
because the first-neighbor distances were effected by the application of the three-body terms. Table 5 shows the sets of adjustable parameters used for the pair potential terms, according to the nature o f the three-body terms applied.
For each configuration, some physical quantities such as the density, Na diffusion coefficient and viscosity were measured. The diffusion coefficient of Na, the most highly mobile species, was measured at a high temperature ( T = 1073 K), using Einstein's relation:
2.3. Modeling of sample preparation
DNa = ( r 2 ) / 6 t ,
The initial system was a cristobalite structure in which B 2 0 3, N a 2 0 and ZrO 2 groups were introduced on random sites. The system temperature then increased quickly due to structural instabilities. The cristobalite structure contains some linear S i - O - S i bonds, unstable with regard to the B o r n - M a y e r Huggins potential. The temperature was controlled by controlling the atomic velocities, decreasing to 5500 K when it exceeded 6000 K. W h e n the temperature stabilized between 5500 and 6000 K, the liquid was relaxed for 6000 time steps of 10-15 s before quenching to room temperature at a rate o f 5 X 1014 K s-~. A series of structural configurations was stored at intermediate temperatures. A constant volume algorithm was used at each stage of the preparation with a density slightly above the equilibrium density to avoid void formation. The equilibrium volume of each configuration
where ( r 2) is the mean square displacement of the Na atoms and t is the duration of the simulation. The viscosities were determined by a method previously used by Gosling on a simple LennardJones fluid [19], in which an external perturbation was applied in the form o f a sinusoidal profile of additive forces. In response to this perturbation, atomic fluxes developed in the system, characterized by a sinusoidal profile of velocities, in phase with the force profile. The viscosity was then calculated by the formula
Table 5 List of parameters used for pair potentials Parameter Si O B q: +4 -2 +3 n: 8 8 2 o- (A) (system I): 1.148 1.42 0.74 o" (,~)(system 2): 1.1008 1.42 0.74 o- (A) (system 3): 1.1008 1.42 0.72
Na +1 8 1.17 1.17 1.17
Zr +4 8 1.45 1.4 1.425
For systems I and 2: r o = 0.29 ,~, b = 0.221 eV. For system 3: rij = 0.29 ,~ except ro_ o = 0.35 ,g,, b = 0.221 eV.
rl = ( p l k Z ) ( F / v ) ,
(3)
(4)
where F and v are the amplitudes of the force and velocity profiles, respectively, p is the atomic density of the system and k = 2 w / L is the wave vector of the perturbation. In the following, F is considered equal to ! x 10 -4 in cgs. Strictly, this last equation must be applied to the limiting case of a null perturbation and a null wave vector. For the viscosity determinations, the constant volume algorithm was always used.
3. Results
3.1. Structural properties The results are distinguished according to whether or not three-body potentials were applied around the
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 11996) 239-248
243
Table 6 Structural characteristics of system 1 (experimental data in parentheses) Property
Si-O
B-O
Na-O
Zr-O
Ncoord:
4.26 (4.0)
3.37 (3.67)
5.86
6.38 (6.3)
Rcut (,~):
2.1
1.9
3.1
3.0
_
1.40 (1.40)
2.48
2.17 (2.08)
'~
dlst neigh (,~k): 1.67 (1.60)
Si or B atoms. The first set of results concerns the structural characterization of a system when only pair potentials are used. Table 6 indicates the firstneighbor distances and coordination numbers around the cations, together with the experimental observations. Concerning first-neighbor distances, note that the S i - O and Z r - O distances are slightly too long. Moreover, the coordination number is too low around the B atoms and too high around the Si atoms. Every Si atom in the actual glass is tetra-coordinated, while, in the model, an excessive percentage of Si atoms have a coordination number equal to 5. The Zr coordination number is consistent with the experimental value. Therefore, we constructed a second set of models, adding to the pair potential a three-body term around the Si atoms. We constrained the O - S i - O angles to the value of the perfect tetrahedrai angle, i.e., 109.47 °, and a weak constraint was applied to the S i - O - S i angles, again with the most energetically favorable angle equal to 109.47 ° . For this series of models, we slightly decreased the atomic radii of Si and Zr to correct the deviations from experimental results noted above. The new structural characterizations are given in Table 7. Re-adjusting the atomic radius of Si corrected the excessive first-neighbor distance obtained with the previous set of parameters. Moreover, because of the angular constraint, every Si atom be-
Table 7 Structural characteristics of system 2 Property
Si-O
B-O
Na-O
Zr-O
Ncoord: Rcut (A):
4.04 2.1
3.48 1.9
5.47 3.1
5.88 3.0
dlst neigh ( ~ t ) :
1.60
1.39
2.49
2.12
~
z~
.... 8~1
i .... 911
m
, .... 104)
, .... I iiJ
i .... 1211
i .... 1311
i .... 140
I ,, 1511
Angle (deg)
Fig. 1. Angular distribution of atomic triplets.
came tetra-coordinated. The first-neighbor distances agreed well with experimental values. The coordination number around the B atoms was still too low, with 48% tri-coordinated B and 52% tetra-coordinated B. Fig. 1 shows the angular distributions of the O-Si-O, O-B3-O and O - B 4 - O angles, in which B 3 and B 4 are, respectively, the tri- and tetra-coordinated boron atoms. The mean O - S i - O angle was very near the perfect tetrahedral angle; the mean O-B3-O angle was near 120 °, signifying that the B 3 atoms were surrounded by a triangle of O atoms, and the m e a n O - B 4 - O angle was near the perfect tetrahedral angle. B 4 atoms were surrounded by tetrahedrons of O atoms. The mean S i - O - S i angle was equal to 158 °. The Zr coordination decreased compared with the previous value. The last stage was to correct the coordination number around the B atoms, which was just below 3.5 rather than the expected value of 3.67 [20]. This correction was obtained by applying a three-body term on the O - B - O angles. We noted that a relatively large constraint was required to increase the B coordination number, compared with what was needed to adjust the Si coordination number to 4. This is because the radius of the B coordination sphere is shorter and the mean coordination number due to geometrical packing of atoms is lower than the physical number (contrary to Si, for which a three-body term was applied to decrease the coordination number). It is more difficult to force a coordination number larger than the natural packing value. We also changed the reference angle for the weak
244
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996) 239-248
Table 8 Stluctural characteristics of system 3
3O
Property
Si-O
B-O
Na-O
Zr-O
Ncoord: Rcut (,~):
4.0 2.1 1.61
3.67 1.9 1.42
5.71 3.1 2.5
5.62 3.0 2.14
dlst neigh (m):
•
,
,
-
,
"
,
"
"
"
,
,
"
"
211 I B-O
svO
constraint imposed on the S i - O - S i angles from 109.47 ° to 160° , without any significant effect. In addition to applying a three-body term around B, we increased the ro_ o parameter associated with o x y g e n - o x y g e n repulsion, diminishing the repulsion strength between two O atoms. This made it easier to increase the coordination number around B, but modified some of the first-neighbor distances and coordination numbers. Other attempts were, therefore, made to readjust the atomic radii. The structural characteristics of the new model are given in Table 8. The first-neighbor distances are in good agreement with the experimental distances, except for a slightly overestimated Z r - O distance. The adjustable parameters assigned to the Zr species represent a trade-off between the coordination number and the first-neighbor distance. Increasing either parameter decreases the other. The final result is characterized by a low coordination number (our models contain few Zr atoms, and this value is probably subject to rather large fluctuation errors) and a large Z r - O distance. The coordination numbers around Si and B are well adjusted. The partial radial distribution functions are plotted in Fig. 2. Positions used for the calculations were averaged at 293 K over 300 time steps. No modifications of the angular distributions were observed for O - S i - O and S i - O - S i , and the mean S i - O - S i angle was 157 °, close to the previous value despite the change in the constraint angle from 109.47 ° to 160 °. In the literature, three series of experiments determined two different mean S i - O - S i angles. The first is due to Mozzi and Warren [21], who found a mean S i - O - S i angle of 144° in silica; Greaves et al. [22] reported a mean S i - O - S i angle equal to 160 °, while a recent NMR study of potassium disilicate indicated a mean S i - O - S i angle of 143 ° [23]. Our result is close to that of Greaves et al. The O - B 3 - O angles were modified by applying
.5
. . . .
i
o
. . . .
-
I
.
.
.
1
.
.
.
,
l
. . . .
,
1
R
. . . .
,
5
. . . .
6
{A}
Fig. 2. Partial radial distributionfunctions between cations and O atoms. The curve was calculated at room temperature by averaging the atomic positionsover 300 time steps.
the three-body term around the B atoms. We found a mean O - B 3 - O angle of 115 °, less than the 120 ° noted above. The triangular configurations of O atoms around the B 3 atoms are, therefore, distorted by the three-body term. 3.2. P h y s i c a l p r o p e r t i e s
Andersen's algorithm allows the equilibrium volume of a system to be measured at different temperatures. Fig. 3 shows the density versus temperature curves obtained for systems 2 and 3. An interesting although logical consequence of applying a constraint on the O - B - O angles is the densification of
O
O O
O
X45
El C]
[] []
[] t-'l
,
,
i
400
,
,
.
i
,
,
,
i
800
12oo
Temperature
(K)
,
, E I
, 16oo
2800
Fig. 3. Density versus temperature for systems 2 (O) and 3 ([3). • represents the experimentaldensity of the four-oxideglass.
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996) 239-248
system 3 as the thermal expansion coefficient increases. Some experimental results provided by the SCD laboratory provide a basis for estimating the validity of the numerical models. At room temperature, the experimental density of a glass of the same composition as the model is 2.465 g cm -3, not far from the computed value of 2.508 g cm -3. The density is still close to the experimental value when a three-body term around B is applied. No thermal expansion coefficient is available for our simplified system, so we compared the simulation results with the thermal expansion coefficient of the complete nuclear glass. For the French R7T7 and German GP 98-12 nuclear glasses, the thermal expansion coefficients are, respectively, 8.3 × 10 -6 and 8.1 × 10 -6 K - ~; the thermal expansion coefficients of the simulated systems without and with a three-body term around B are 2.41 × 10 -5 and 8.11 X 10 -6 K -1, respectively. Here again, the application of a threebody term around the B atoms provides a better approximation of the experimental thermal expansion coefficient. For system 3, we measured the vitreous transition temperature, Tg, by plotting the evolution of the atomic potential energy versus the temperature; the results are shown in Fig. 4. The intersection of the linear extrapolations at low and high temperatures corresponds to a Ts equal to roughly 1000°C, probably higher than the experimental results. The shift of Tg towards higher temperatures is a classic simulation artifact, for which very fast quench rates, com-
-5.98
10 '
.5.99
IO i
4.#10
D) I
4~AIJ
I0 t
T:
I
I
I
T
[
I
I
•
i
~
Tempcral~r¢
I
I
I
I
I
?
I
I
I
(h.~
Fig. 4. Atomic potential energy versus temperature for system 3 (Tg determined by the intersection of the two extrapolated linear
fits).
245
1114
8 •
•
0
•
0
I0 ~
o 5.5
6
6.5
7
7.$
8
8,5
9
9.5
lfr (*10') Fig. 5. Logarithm of Na atom diffusion coefficient versus reciprocal of T: 0 , points measured in system 2; A, points measured in system 3; O , experimental points measured on complete nuclear
glasses [28].
pared with the experimental rates, are imposed to hasten the creation of low-temperature structures. The Tg of the French nuclear glass is 502°C; the German nuclear glass, which contains more network modifiers, has a Tg of 543°C. We also examined the Na diffusion coefficients, DNa, at high temperatures. The Na atoms are small, and their Coulomb charge is low; the interactions binding them to the network are weak, and they can move easily. DNa is the most easily measurable diffusion coefficient. An Arrhenius diagram of DNa for systems 2 and 3 versus the reciprocal of the temperature, T, is plotted on a logarithmic scale in Fig. 5. The values measured with the German nuclear glass are shown for comparison. To obtain equilibrium values of DNa at 1073 K, longer runs are necessary to allow proper initial relaxation of the systems. Several observations may be made concerning this figure. First, the three slopes are basically similar, indicating good agreement for the activation energy between the computed and physical systems. A linear fit of the data provides an activation energy of 0.92 eV in system 2 and 1.03 eV in system 3. The corresponding activation energy in the actual nuclear glass is 1.08 eV. Second, the application of a threebody term around B significantly lowers DNa tO below experimental values. This qualitative behavior (computed results lower than experimental results) is
246
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996) 239-248
not surprising, since the complete nuclear glasses contain a larger fraction of network modifiers, and so is probably more disordered, insofar as we can speak of a level of disorder for glass. The diffusion of Na atoms is facilitated by the increased network disorder. The last physical value calculated was viscosity, using the method described above. The establishment of a dynamic flux in response to the application of the perturbation is a very long process, as the reader can see from Fig. 6. The y-axis represents the amplitude of the sinusoidal velocity profile versus the time of simulation. The velocity profile is always determined by dividing the difference between the current and initial positions by the time elapsed between the two states. The current and initial positions are averaged over 100 time steps to get rid of part of the thermal motion. After a sudden increase due to the application of the perturbation, the signal progressively decreases and continues to decrease beyond the 38 000 time steps of the simulation. The amplitude representing the equilibrium dynamic profile corresponds to the asymptotic value of the curve. In fact, two different phenomena interfere in the process, the first being the elastic deformation of the system and the second being the progressive establishment of viscous flow. Fig. 6 shows that the viscous flow was hidden by elastic deformation throughout the simulation. One means of evaluating
30011
~. m
2q)O0
I$1HI
I000
.
.
.
.
.
5000
.
.
.
.
.
.
.
IOIlO~
.
.
.
.
.
IGOI)O
.
.
.
.
.
.
.
201HIO
.
.
.
.
250(10
.
.
.
.
.
3000{I
.
i
.35000
.
.
.
.
40OOO
T i m e o f ~imulation (f;)
Fig. 6. Velocity profile amplitude versus simulation time (perturbation applied at t = 0, T = 1923 K).
tm
c
20
0
*
I
,
I loot)
,
i
*
I
,
i
,
MHH)
Time of the simulation
n 600O
,
i
,
I ~NH)
,
i
, I001~)
(fs)
Fig. 7. Velocity profile amplitude versus simulation time after release of the perturbation at t = 0 (T = 1923 K).
the amplitude of the velocity profile relevant to the viscous flow is to determine the asymptotic value of the curve by a suitable adjustment on the simulated points. Unfortunately, however, very different values are obtained when different fitting curves are used (1/t or e x p ( - a t ) in particular), while there is very little difference in the quality of the fit. Prolonging the application of the perturbation would be too time-consuming for a hypothetical stabilization of the velocity profile. Therefore, we attempted to obtain the viscosity by stopping the perturbation at a time, t, equal to 25 000 time steps and following the evolution of the velocity profile amplitude, calculated as the difference between the atomic positions of the henceforth undisturbed system and the initial positions, divided by the time during which the perturbation is applied. This new graph is shown in Fig. 7, where it may be seen that relaxation is much faster than before and the amplitude of the velocity profile reaches an equilibrium value in a rather short time. This second stage of the experiment, the release of the perturbation, furnishes an equilibrium amplitude of the velocity profile which can be taken with confidence for the viscosity calculation. The viscosity of system 3 was 3.3 Pa s at 1823 K and 2.1 Pa s at 1923 K compared with experimental viscosities estimated at 1.93 and 1.2 Pa s, respectively, on systems of the same composition as the model; the experimental values were extrapolated from V - F - T parameters. The order of magnitude of
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996)239-248
the viscosity of the simulated systems is correct, although the viscosity of the models is overestimated.
4. Discussion Computer models were constructed and demonstrated to produce predictions of the structural characteristics and some physical properties. It was necessary to introduce a three-body term applied to the O - B - O angles to increase the coordination number of the B atoms. The consequence of this modification was the deformation of the triangles of O atoms surrounding the tri-coordinated B atoms. We tried to apply a three-body potential with two local minima, at 109.47 and 120.0 °, to the O - B - O triplets to force either local triangular order or local tetrahedral order. Unfortunately, this method failed because the two minima were too close to one another, compared with the width of the O - B - O angle distribution. These two minima of the potential did not result in production of two separated O - B - O angle distributions. A secondary effect of the constraint applied on the O - B - O angles is lower Na diffusion coefficients; this lowering is a logical result, since the Na atoms travel via interstitial sites surrounded by O atoms and the added constraint decreases the mobility of the O atoms, making the Na paths more rigid and difficult to cross. System densification following the application of the three-body O - B - O term is a second possible explanation for the lower Na diffusion rate. As noted above, the computed glass transition point, Tg, was higher than the experimental temperature for an actual nuclear glass. This difference probably proceeds from the very fast quenching rate applied to terminate the computed systems. Scholze [24] proposed as an example the evolution of the for borosilicate glasses, and showed a difference of nearly 100 K between quench rate of 10 - 4 K min-l and 10 K min -~ Extrapolating the experimental curve to 1014 K s -I further increases Tg by 400 K. This difference is of the order of magnitude of that observed between the computed and experimental Nevertheless, the Na diffusion coefficients were
247
lower but not far from the experimental values measured on German nuclear glasses, which are not only less rigid than the model but also contain more than one alkaline species, suggesting that the mixed alkali effect and the influence of the other components are expected to be weak in complex nuclear glasses. The latter conclusion has previously been stated by Soules and Busbey [8] who did not detect a strong mixed alkali effect in simulations of diffusion in Na20 + K 2 0 + SiO 2 systems. The mixed alkali effect is known to decrease when the temperature increases, tending, in particular, to disappear when the diffusion mode becomes liquid-like. We observed that viscosities at 1823 and 1923 K were higher than experimental values (extrapolated from a Vogel-Fulcher-Tammann curve). Some care must be taken with these values because the exact value of the viscosity must be determined within the limit of a null perturbation. However, as the level of the perturbation we applied was approximately 2% of the mean force experienced by the atoms, this value should not be far from the exact one. It is possible to quantify the applied shear stress, s. The external perturbing force ( F ) profile is a sine wave, so the applied shear stress is not constant over the entire sample. Nevertheless, an order of magnitude can be given by considering the total force exerted on the atoms contained in the half simulation cell and experiencing a positive perturbing force, dividing by the cell cross-section: s = 6.7 × 10 9 N m - 2 . The theoretical shear strength of solid oxide glasses is typically about 10 ~° N m -2. Second, due to the low velocity profile imposed, it must be calculated by taking differences between two averaged positions separated by a rather large time step; thus, thermal agitation may have some influence. When an atom moves in a direction orthogonal to the direction of the external perturbing force, the magnitude of the perturbing force experienced by this atom changes, widening the velocity profile, and the profile amplitude is underestimated. We expect this effect to be negligible due to the low diffusion coefficient of the network formers, even at a temperature of 1923 K. The last point we would like to emphasize is the inability of the potential to represent physical characteristics of actual glasses. In particular, the cohesive energy is not well reproduced because, by taking
248
J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 195 (1996) 239-248
formal charges (Si 4+, 0 2-, B 3+, Na +, Z r 4 + ) , w e attribute to the atomic interactions a purely ionic feature, offset by the three-body terms representing the semi-covalent nature of some atomic bonds. The reality is different, and the covalent nature of the atomic bonds is related to the geometrical form of the electron densities of states localized in particular areas on the covalent bonds. The consequence is that the charge location is much more complex than the model used in this simulation and, consequently, the cohesive energy is lower. Several attempts to determine more exact charges have been published. Empirical potentials have been derived from ab initio calculations for water [25] and pure silica [26]. Nevertheless, it is a very challenging task to take into account these physical characteristics on mixed systems because, as has been shown experimentally [27], the charge of a core depends on its local environment. It is, therefore, incorrect to assign fixed charges to O or B atoms in particular because the surroundings of these atoms are different from one site to another.
5. Conclusion An overall correct structure of a complex oxide glass involving four different oxide species was modeled using pair potentials completed by some three-body terms. These structures show good agreement with experimental results for the density, thermal expansion coefficient and Na atom diffusion and are on the same order of magnitude for the viscosity.
Acknowledgements One of the authors (J.-M.D.) is grateful to the CEA ( D C C / D R D D / S C D ) for the financial support of this study. Dr N. van Doan kindly provided the original program for the development of the simulations and, together with Dr Y. Limoge, provided valuable insights after a critical reading of this paper.
References [1] F. Pacaud, C. Fillet and G. Baudin, in: Proe. 2nd Int. Seminar on Radioactive Waste Products, Ji~lich, Germany (1990) p. 577. [2] F. Pacand, C. Fillet and N. Jacquet-Francillon, Scientific Basis for Nuclear Waste Management, Vol. XV, 257, ed. C. Sombret (1992) p. 161. [3] T.F. Soules, J. Chem. Phys. 71 (1979) 4570. [4] T.F. Soules, J. Non-Cryst. Solids 49 (1982) 29. [5] K.V. Damodaran, V.S. Nagarajan and K.J. Rao, J. Non-Cryst. Solids 124 (1990) 233. [6] Q. Xu, K. Kawamura and T. Yokokawa, J. Non-Cryst. Solids 104 (1988) 26l. [7] T.F. Soules and A.K. Varshneya, J. Am. Ceram. Soc. 64 (1981) 145. [8] T.F. Soules and R.F. Busbey, J. Chem. Phys. 75 (1981) 969. [9] B.P. Feuston and S.H. Garofalini, J. Chem. Phys. 89 (1988) 5818. [10] F.H. Stillinger and T.A. Weber, Phys. Rev. B43 (1991) 1194. [I 1] P. Vashishta, R.K. Kalia, J.P. Rino and I. Ebbsj/5, Phys. Rev. B41 (1990) 12197. [12] D.C. Anderson, J. Kieffer and S. Klarsfeld, J. Chem. Phys. 98 (1993) 8978. [13] A. Alavi, L.J. Alvarez, S.R. Elliott and I.R. McDonald, Philos. Mag. B65 (1992) 489. [14] S. Balasubramanian and K.J. Rao, J. Phys. Chem. 98 (1994) 10871. [15] D. Petit-Maire, PhD thesis, University of Paris VI (1988). [16] H.C. Andersen, J. Chem. Phys. 72 (1980) 2384. [17] M.J.L. Sangster and M. Dixon, Adv. Phys. 25 (1976) 247. [18] J. Kieffer and C.A. Angell, J. Chem. Phys. 90 (1989) 4982. [19] E.M. Gosling, I.R. McDonald and K. Singer, Molec. Phys. 26 (1973) 1475. [20] W.J. Dell, P.J. Bray and S.Z. Xiao, J. Non-Cryst. Solids 58 (1983) 1. [21] R.L. Mozzi and B.E. Warren, J. Appl. Crystallogr. 2 (1969) 164. [22] G.N. Greaves, A. Fontaine, P. Lagarde, D. Raoux and S.J. Gurman, Nature (1981) 611. [23] I. Farnan, P.J. Grandinetti, J.H. Bultisberger, J.F. Stebbins, U. Werner, M.A. Eastman and A. Pines, Nature 358 (1992) 31. [24] H. Scholze, Glass: Nature, Structure and Properties (Springier, Berlin, 1991) (see p. 50). [25] G.G. Hall and C.M. Smith, Int. J. Quantum Chem. 42 (1992) 1237. [26] K. de Boer, A.P.J. Jansen and R.A. van Santen, Chem. Phys. Lett. 223 (1994) 46. [27] C.H. Hsieh, H. Jain, A.C. Miller and E.I. Kamitsos, J. Non-Cryst. Solids 168 (1994) 247. [28] Hj. Matzke and E. Vernaz, J. Nucl. Mater. 201 (1993) 295.