Molecular dynamics simulations of reactive wetting in metal–ceramic systems

Acta mater. 48 (2000) 4419–4424 www.elsevier.com/locate/actamat

MOLECULAR DYNAMICS SIMULATIONS OF REACTIVE WETTING IN METAL–CERAMIC SYSTEMS T. P. SWILER1 and R. E. LOEHMAN2* 1

Advanced Materials Laboratory, University of New Mexico, Albuquerque, NM, USA and 2Sandia National Laboratories, Advanced Materials Laboratory, 1001 University Boulevard SE, Suite 100, Albuquerque, NM 87106, USA

Abstract—We have simulated the wetting of sapphire by aluminum using the ES+ potential. This potential is unique among empirical potentials because it applies a combination of the embedded atom method (EAM) and Coulomb potentials on a system of atoms with defined electronegativities so that the atom charges as well as the interatomic forces and system energies found by typical empirical potentials can be determined. Because we do not have to specify atom bond character a priori, we can directly model reactive wetting where atoms from the metal droplet may replace cations in the ceramic substrate and undergo the appropriate change from metallic to ionic bonding. Also, by using this empirical potential, we can simulate large dynamic systems that cannot be simulated by first-principles methods, albeit with some loss in detail and accuracy. In our simulations of high-temperature wetting we observe the formation of an oxygen-deficient reaction layer between the liquid and the substrate. The driving force for the creation of this layer is the partial oxidation of the metallic aluminum, which results in partial reduction of the aluminum ions in the substrate and diffusion of oxygen from the substrate to the reaction layer. We believe that this change in stoichiometry at the metal–ceramic interface may influence diffusivities in the surface of the ceramic and may therefore have an important effect on the kinetics of the evolution of surface features observed experimentally during reactive wetting.  2000 Acta Metallurgica Inc. Published by Elsevier Science Ltd. All rights reserved. Keywords: Wetting; Ceramics; Metals; Surfaces & interfaces; Computer simulation

1. INTRODUCTION

Classical wetting involves the behavior of a bulk liquid on a flat, non-deformable solid in equilibrium with a vapor phase. Here, the equilibrium state of the system is simply defined by the balance of interacting surface tensions that results in a well-defined equilibrium liquid contact angle. In reality, the situation is more complex, and wetting may typically include other processes such as adsorption and reaction. For example, a precursor film can form in the absence of a bulk liquid, making this part of wetting fundamentally an adsorption process. Similarly, the reaction of a liquid with a substrate can form microscopic surface features on the substrate that influence wetting on the macroscopic scale. These deviations from perfect classical wetting can confuse attempts to describe the behavior of a system with a classical wetting model. Thus, we are studying wetting behavior in a complex system using a fundamental simulation model to provide insight into these different wetting phenomena.

* To whom all correspondence should be addressed. Tel.: 001 505 272 7601; fax: 001 505 272 7304.

Molecular dynamics (MD) is a simulation method by which we can investigate wetting at the atomic scale. Given appropriate interaction potentials, one can follow all of the atomic motions that lead to continuum phenomena. The use of MD in simulations of wetting is not new. Many existing MD simulation studies of wetting have used Lennard-Jones potentials [1–9]. However, monoatomic Lennard-Jones liquids are known to have unrealistically high vapor pressures, so it is sometimes difficult to distinguish the liquid–vapor interface due to the high density of the vapor phase in these systems. Additionally, spreading kinetics may be unduly influenced by the dense vapor phase in these systems. More realistic embedded atom method (EAM) potentials [10] have been formulated for several metallic systems [11] and their applicability for modeling liquid metals has been demonstrated [12]. These potentials couple an atom-densitydependent embedding potential with a simple pairwise potential to model the energy of an atom in a metallic system where each atom interacts with the overall electron gas as well as the atoms in its vicinity. One advantage of EAM potentials over LennardJones potentials is that the vapor pressures of

1359-6454/00/$20.00  2000 Acta Metallurgica Inc. Published by Elsevier Science Ltd. All rights reserved. PII: S 1 3 5 9 - 6 4 5 4 ( 0 0 ) 0 0 2 2 8 - 7

4420

SWILER and LOEHMAN: MOLECULAR DYNAMICS SIMULATIONS

monoatomic EAM liquids are lower and more realistic. We have previously performed simulations of metal–metal wetting using EAM potentials and have obtained liquid droplets with low vapor pressures, clearly defined equilibrium wetting angles, and realistic spreading rate dependencies [13]. Interactions of metals with ceramic or glass substrates have also been modeled using MD. Typically, these models have used an ionic potential to model interactions within the ceramic, an EAM or LennardJones potential to model interactions within the metallic phase, and a Lennard-Jones potential to model interactions between the neutral metal atoms and the ions in the ceramic [14]. The aim of those simulations was primarily to study the structural effects resulting from the weak interactions between the metallic and ceramic or glass phases, and issues regarding the chemistry at the interfaces were mainly ignored. Because bonding in metals and ceramics is significantly different, the use of empirical potentials in MD simulations to accurately model ceramic–metal interactions has not been particularly successful. Firstprinciples calculations have typically been the method of choice for obtaining reliable predictions of structures and transition energies along specific paths for processes at ceramic–metal interfaces [15, 16]. Firstprinciples molecular dynamics simulations are now possible, but the temporal and spatial scales of these simulations are still quite limited. One notable approach used to explain the interactions between a ceramic and a metal phase is the concept of image potentials [17]. Because metals have large dielectric constants, it is proposed that a surface charge is induced in the metal as a response to the partially ionic ceramic surface it contacts. Since this surface charge is opposite in sign but otherwise mirrors the charges at the surface of the ceramic, the term “image potentials” is derived, and these potentials are always attractive. Although this concept gives insight into the nature of attraction between a metal and a ceramic, it assumes that the metal is a continuous medium and therefore cannot be expressed as a potential suitable for MD simulations. In this work, we sought a method by which we could model the reactive wetting of alumina (Al2O3) by metallic aluminum. We have used empirical potentials that allow us to build on the techniques used in simulations of wetting in EAM metallic systems. We first considered the modified EAM (MEAM) potential of Baskes [18, 19] as applied to the Al–O system [20]. Although this potential may be successful for modeling the structure of low-temperature interfaces, we did not find it useful for modeling reactive wetting. We wanted to capture the change in behavior of metallic aluminum when it oxidizes in the presence of oxygen, but the MEAM potential does not allow for the electron transfer that is central to the concept of oxidation. We then implemented the ES+ potential of Streitz and Mintmire [21–23], which couples an EAM potential with an ionic potential and allows for

charge transfer between atoms. The advantage of this potential is that it captures the essence of how we view bonding in metals and ceramics while not forcing formal ionic charges that may be not be appropriate over the whole system. It also allows for reaction and exchange of metal atoms between the oxide and metallic phases. 2. PROCEDURE

We used molecular dynamics to simulate reactive wetting in the Al–O system. An appropriate set of potential functions is required to give a simulated system an identity. In these simulations, we used the ES+ potential of Streitz and Mintmire for the Al–O system [21]. This potential is the sum of an electrostatic potential and an EAM potential, E = Ees + EEAM. To understand the charge-transfer mechanism that this potential permits, it is instructive to summarize the functional dependencies for the electrostatic portion of the potential, Ees = E0 +

冘

qici +

i

冘

1 qqV , 2 i,j i j ij

(1)

where qi is the sum of a core charge zi and a variable charge distributed symmetrically about atom i with the distribution function →

→

fi(兩 r ⫺ r i兩) =

E0 =

冘

→ → z3i exp(⫺2zi兩 r ⫺ r i兩), π

Ei(0) +

i

ci = c0i +

冘

(2)

1 Ecc(z , z , r , f , f ), 2 i⫽j ij i j ij i j

(3)

冘

(4)

cvc ij (zj, rij, fi, fj)

j

and Vij = J0i dij + Vvv ij (rij, fi, fj). →

(5)

→

The charge distribution function fi(兩 r ⫺ r i兩) of equation (2) could be constructed from Slater 1s orbitals and was chosen for mathematical convenience. Note that ci of equation (4) is akin to the electronegativity of atom i and is dependent on the environment about atom i, while J0i of equation (5) is the self-Coulomb repulsion or atomic hardness of atom i. The charges, qi, are calculated to minimize the total energy of the electrostatic system as given by equation (1). The reader is referred to the original reference for the full vc vv definition of the functions Ecc ij , cij and Vij , as well as

SWILER and LOEHMAN: MOLECULAR DYNAMICS SIMULATIONS

the parameters c0i , J0i , zi and zi used for the Al–O system. This potential reproduces the formal ionic charges of about +3 for Al and ⫺2 for O in the bulk Al2O3 structure, while Al in the bulk metallic phase remains essentially neutral. Our implementation of this potential involves determining atom charges every time step by solving a set of N linear equations, where N is the number of atoms. Although this approach is the most robust, the amount of computation required scales as N3, limiting the size of the systems that we can simulate to about 3000 atoms. The EAM portion of the total potential includes a Finnis–Sinclair [24] type EAM potential and a pairwise potential that handles core repulsion as well as a pairwise attractive term. As in earlier works, the EAM portion of the potential, which primarily comes into play when modeling the metallic liquid, results in a low vapor pressure so that the wetting line can be clearly delineated. This potential includes Coulombic interactions with 1/r dependence. Streitz and Mintmire performed a full Ewald sum [25, 26] to calculate Coulombic interaction energies. Briefly, an Ewald sum combines a real-space sum and a reciprocal-space sum to obtain rapid convergence of 1/r interactions in a periodic system. The relative extents of the real-space and reciprocal-space portions of the sum can be altered by a parameter, b, which only affects the rate of convergence and not the final answer. The real-space portion of the sum takes precedence over the reciprocalspace portion when b is small and vice versa. In our implementation of this potential, we wanted to avoid a full Ewald sum because of the assumptions in the dimensionality of the periodicity inherent in the reciprocal-space portion of the sum. Because the reciprocal-space portion of the sum can be made arbitrarily small by making b sufficiently small, we were able to accurately determine Coulombic interactions by using just the real-space portion of the sum with a small value for b (0.20) and a large real-space cutoff ˚ ). distance (12 A The desire to avoid the problems associated with such a long-range potential prompted us to investigate other methods to treat 1/r Coulombic interactions. A simple screening function has recently been shown to give a good convergence in summing Coulombic ˚ ) cutoff distances [27]. interactions at short (苲6 A Although this approach is successful in determining sums of the Coulomb interactions, we found that it was not successful in accurately recreating the matrix of interactions between individual atom pairs that the ES+ potential requires and therefore could not be simply substituted for the Ewald sum approach. In our simulations, we used 1440 atoms in the Al2O3 substrate and 252 atoms in the Al liquid for a total of 1692 atoms. The dimensions of the substrate ˚ long (x) by 14 A ˚ wide (y) were approximately 66 A ˚ thick (z). The simulation cell was periodic by 13 A in the x- and y-dimensions to prevent edge effects that would have revealed the limited size of the substrate

4421

to the spreading liquid. To simulate wetting on a bulk material, the atoms in the bottom half of the Al2O3 substrate were constrained while the atoms in the top half were free to move and react with the Al liquid placed on top. The thickness of the substrate was sufficient to ensure that the minimum distance between the dynamic atoms and the bottom of the substrate was beyond the range of the short-range portion of the interaction potential. We then put the liquid on top of the substrate. To minimize the amount of unused substrate area in the simulation and thereby reduce the number of atoms required, we simulated the spreading of a line in the x-dimension rather than of a circular droplet in both x- and y-dimensions. We formed this line of liquid by creating a crystalline block of the liquid component that spanned the y-periodic dimension and then relaxing the block at a temperature in excess of the melting temperature. In all subsequent discussion, we shall refer to this line of liquid as the drop. The resulting geometry (after assumed spreading) is shown in Fig. 1. With this geometry, we found our system size to be sufficient to provide both the area of substrate and the volume of liquid for unconstrained spreading on realizable time scales. We performed these wetting simulations at the temperatures 1000 K and 1600 K. The experimental melting points of Al metal and Al2O3 are 933 K and 2323 K, respectively, and the simulated range represents the temperatures at which wetting experiments in this system are often performed. The entire system was thermostatted using a Nose–Hoover thermostat, as necessitated by the size of the simulation. Although this can perturb relevant atomic motions while the system is evolving, the perturbation becomes small as equilibrium is approached, and we believe that its effects on the system are small. 3. RESULTS AND DISCUSSION

Our simulations make several predictions that are in agreement with current experimental and theoretical descriptions of the Al2O3 surface and its behavior during reactive wetting. It is generally accepted that the basal plane (0001)

Fig. 1. The simulation geometry used in these spreading simulations. The system is periodic in the horizontal plane, effectively simulating an array of infinitely long lines of liquid on an infinite substrate plane.

4422

SWILER and LOEHMAN: MOLECULAR DYNAMICS SIMULATIONS

Fig. 2. Location of atomic layers in the z-dimension of the substrate relaxed at 300 K. Negative coordinate values correspond to the constrained, unrelaxed portion of the substrate, while positive coordinates correspond to the dynamic portion of the substrate.

of α-Al2O3 is terminated by a half layer of aluminum atoms that has relaxed inward so as to provide bulk neutrality while also minimizing the formation of a dipolar layer at the surface. As demonstrated by firstprinciples calculations by Manassidis and Gillan [28], the outer aluminum half-layer relaxes inward, the oxygen layer immediately beneath stays put, and the next two aluminum half-layers relax towards each other. We see a similar behavior for our simulated substrate at 300 K as shown in Fig. 2. Here, the outer Al layer splits so that its average location is near the location predicted by Manassidis and Gillan, the O layer immediately beneath stays near its unrelaxed position, and the next two Al half-layers appear to relax into a single broadened layer. Because of the charge-transfer mechanism in the ES+ potential, the actual charge developed on the atoms in the relaxed state is of interest. In Fig. 3, we show the system with atoms shaded to identify atom charge. More quantitatively, we find that the atom charges approach their formal values of +3 for Al and ⫺2 for O, as is shown in Fig. 4, but the magnitudes of the atom charges are reduced near the surface as compared with bulk values. Note that the some scatter in the atom charges is observed near the center of the ˚ ). This indicates that bulk behavior substrate (兩z兩⬍2 A has not yet been fully reached in a system of this size. We investigated the density of the aluminum drop-

Fig. 3. Dynamic portion of the relaxed Al2O3 substrate relaxed at 300 K with atom charges indicated by shade. The larger, darker atoms are oxygen while the smaller, lighter atoms are aluminum.

Fig. 4. Atom charge as a function of z-coordinate in the dynamic portion of the substrate.

let by counting the number of atoms in rectangular cells wholly contained within the liquid. We found that the Al liquid density was in the range of 2.60– 2.73 g/cm3 at 1000 K and 2.61–2.68 g/cm3 at 1600 K. These values are both significantly higher than the experimentally determined density for liquid aluminum of 2.35 g/cm3 at 1000 K [29], but are slightly less than the density for face-centered cubic (fcc) Al of 2.7 g/cm3. At this scale of simulation, we cannot determine whether this deviation is the result of limitations of the potential function or the result of nonbulk liquid properties in the very small drop. When we simulated the wetting of Al on Al2O3 at 1600 K, we found that some of the oxygen atoms from the substrate diffused into the aluminum droplet and subsequently formed a reaction layer that developed and extended with time, as shown in Fig. 5. This resulted in the formation of a suboxide of aluminum in the surface region of the substrate in contact with the droplet, with an associated decrease in atom charge, as shown in Fig. 6 for the system after 5 ps of wetting. In simulations of wetting in this system at 1000 K, we found that the droplet was more compact, i.e., its surface was better defined, than at 1600 K, as shown in Fig. 7. The compactness of the droplet resulted in fewer atomic excursions on to the substrate surface, limiting the extension of the reaction layer. In contrast to 1600 K results, where we found that the reaction layer continued to develop during the entire length of our simulation, at 1000 K the reaction layer changed little after wetting for 3 ps. Ordering in the reaction layer can be seen in Fig. 8, which shows that oxygen from the Al2O3 substrate extends into the aluminum droplet. We see that the region at the base of the droplet that was initially free of oxygen clearly contains correlated layers of oxygen after 5 ps of simulated wetting at both 1600 K and 1000 K. When oxygen leaves the Al2O3 substrate to partially oxidize aluminum in the droplet, the substrate becomes partially reduced, forming a suboxide phase. The stability of the Al2O3 at the interface is therefore reduced and the formation of oxygen vacancies pro-

SWILER and LOEHMAN: MOLECULAR DYNAMICS SIMULATIONS

4423

Fig. 7. Image of wetting system at 1000 K after 5 ps.

Fig. 8. Density of oxygen in the system after 5 ps of wetting by z-coordinate. Fig. 5. Picture of wetting system at 1600 K after (a) 1 ps, (b) 3 ps and (c) 5 ps of wetting. Atom charges are indicated by shade as given by the legend, with the larger darker atoms being oxygen and the smaller lighter atoms being aluminum.

Fig. 6. Charges as a function of z-coordinate in the wetting system after 5 ps at 1600 K.

vides a mechanism for enhanced near-surface diffusion. We believe that this mechanism increases the response kinetics of the substrate surface to the tensions of the droplet and therefore allows the formation of ridges in a reactive wetting system as observed by Saiz et al. [30]. An important question is whether or not a suboxide reaction layer really exists. It is known experimentally that Al2O3 does not form a suboxide in bulk [31]. However, the occurrence of such a suboxide at a high-temperature interface between Al and Al2O3

cannot be dismissed. To settle this question with more confidence, we are planning to investigate the energetics of a suboxide at an Al2O3 surface using firstprinciples methods. Clearly the present system is quite reducing, so it is worth addressing the issue of oxygen content on the non-wetted portion of the Al2O3 surface. One might expect that the terminating species of the Al2O3 substrate is dependent on the chemical potential of oxygen in the system and that, with excess aluminum and no excess oxygen, the surface should be aluminum-rich. In reality, the vapor pressure for oxygen is sufficiently low that desorption of oxygen will not occur on the length and time scales that can be observed in our simulations. A vapor pressure of 10⫺4 atm would average to about one liberated or impinging atom per 500 ns over the simulated surface area, but our simulations were 100,000 times shorter. Thus, although equilibration of surfaces by vapor transport may occur rapidly in experiments, this process cannot be simulated using molecular dynamics. In contrast, the initial stages of reactive wetting can be modeled using MD because this process involves atomistic motions that can be immediately observed in the simulation. 4. CONCLUSIONS

By making use of the ES+ potential that gives realistically low liquid vapor pressures and enables

4424

SWILER and LOEHMAN: MOLECULAR DYNAMICS SIMULATIONS

charge transfer between atoms, we have performed the first molecular dynamics simulation study of reactive wetting in a metal/oxide system. Despite the current limitations on system size using this potential, our simulations have given us insight into the reactive wetting of Al2O3 by Al, and have yielded the following specific results.

앫 Relaxation of the alumina (0001) surface is in qualitative agreement with predictions made by first-principles methods. 앫 Reactive wetting in the simulated Al–Al2O3 system commences by diffusion of oxygen out of the surface of the Al2O3 substrate into the metallic aluminum liquid, forming an aluminum suboxide reaction layer that extends into the droplet. 앫 Depending on the temperature of the system, this reaction layer may extend beyond the edge of the droplet as aluminum atoms from the liquid make excursions on the Al2O3 substrate that pull oxygen atoms out of the substrate. The reaction layer extends at 1600 K but does not extend at 1000 K, where the droplet shape remains more compact. Because the reaction layer can extend outward at 1600 K but not at 1000 K, the system continues to evolve over the 5 ps duration of the simulations at the higher temperature, but the evolution of the system stops by 3 ps at the lower temperature. 앫 No desorption of oxygen from the surface of the Al2O3 substrate occurs over the limited duration of the simulations in accordance with expectations given the low vapor pressure of oxygen in this system. 앫 An aluminum suboxide reaction layer forms at the liquid–substrate interface. Future work will include better validation of the ES+ potential for modeling reactive wetting in the Al/Al2O3 system and simulation of phenomena that require larger systems and longer times. Validation will involve comparing the potential energies of specific configurations obtained from the ES+ potential and first-principles methods. Simulation of larger systems and longer times will be accomplished using algorithms that scale better with increased system size, such as the those used by Campbell et al. [32] in their implementation of the ES+ potential. Such simulations will better define the equilibrium state of the current system and will enable observation of the processes that produce the surface features on polycrystalline Al2O3 substrates that are observed experimentally.

Acknowledgements—This work was supported by the DOE Office of Basic Energy Sciences and Contract DE-AC0494AL85000 at Sandia National Laboratories. REFERENCES 1. Sikkenk, J. H., Indekeu, J. O. and van Leeuwen, J. M. J. Phys. Rev. Lett., 1987, 59, 98. 2. Nijmeijer, M. J. P., Bruin, C., Bakker, A. F. and van Leeuwen, J. M. J. Phys. Rev. A, 1990, 42, 6052. 3. Nijmeijer, M. J. P., Bruin, C., Bakker, A. F. and van Leeuwen, J. M. J. Mol. Phys., 1991, 72, 927. 4. Yang, J., Koplik, J. and Banavar, J. R. Phys. Rev. Lett., 1991, 67, 3539. 5. Nieminen, J. A., Abrahan, D. B., Karttunen, M. and Kaski, K. Phys. Rev. Lett., 1992, 69, 124. 6. Bruin, C., Nijmeijer, M. J. P. and M Crevecoeur, R. J. Chem. Phys., 1995, 102, 7622. 7. Philips, J. M. Phys. Rev. B, 1995, 51, 7186. 8. Burlatsky, S. F., Oshanin, G., Cazabat, A. M., Moreau, M. and Reinhardt, W. P. Phys. Rev. E, 1996, 54, 3832. 9. Jin, W., Koplik, J. and Banavar, J. R. Phys. Rev. Lett., 1997, 78, 1520. 10. Daw, M. S., Foiles, S. M. and Baskes, M. I. Mater. Sci. Rep., 1992, 9, 251. 11. Foiles, S. M., Baskes, M. I. and Daw, M. S. Phys. Rev. B, 1986, 33, 7983. 12. Foiles, S. M. Phys. Rev. B, 1985, 32, 3409. 13. Swiler, T. P., Acta. mater., 2000, in press. 14. Webb, E. B. and Garofalini, S. H. Surf. Sci., 1994, 319, 381. 15. Duffy, D. M., Harding, J. H. and Stoneham, A. M. Acta metall. mater., 1992, 40, S11. 16. Harding, J. H., Stoneham, A. M. and Venables, J. A. Phys. Rev. B, 1998, 57, 6715. 17. Tasker, P. W. J. Chem. Soc., Faraday Trans., 1990, 86, 1311. 18. Baskes, M. I. Phys. Rev. Lett., 1987, 59, 2666. 19. Baskes, M. I. Phys. Rev. B, 1992, 46, 2727. 20. Baskes, M. I. in Proceedings of the 1996 International Workshop on Computer Modelling and Simulations for Materials Design, Tsukuba, Japan, ed. S. Nishijima and H. Onodera, National Research Institute for Metals, Tsukuba, Japan, 1996, p. 1. 21. Streitz, F. H. and Mintmire, J. W. Phys. Rev. B, 1994, 50, 11996. 22. Streitz, F. H. and Mintmire, J. W. J. Adhesion Sci., 1994, 8, 853. 23. Streitz, F. H. and Mintmire, J. W. Langmuir, 1996, 12, 4605. 24. Finnis, M. W. and Sinclair, J. E. Phil. Mag. A, 1984, 50, 45. 25. Ewald, P. P. Annln. Phys., 1921, 64, 253. 26. Smith, E. R. Proc. Roy. Soc. Lond. A, 1981, 375, 475. 27. Wolf, D., Keblinski, P., Philpot, S. R. and Eggebrecht, J. J. Chem. Phys., 1999, 110, 8254. 28. Manassidis, I. and Gillan, M. J. J. Am. Ceram. Soc., 1994, 77, 335. 29. Smith, P. M., Elmer, J. W. and Gallegos, G. F. Scripta mater., 1999, 40, 937. 30. Saiz, E., Tomsia, A. P., and Cannon, R. M. Acta mater., 1998, 46, 2349. 31. Yanagida, H. and A Kroeger, F. J. Am. Ceram. Soc., 1968, 51, 704. 32. Campbell, T., Kalia, R. K., Nakano, A., Vashishta, P., Ogata, S. and Rogers, S. Phys. Rev. Lett., 1999, 82, 4866.