versão impressa ISSN 0103-9733
Braz. J. Phys. v.31 n.3 São Paulo set. 2001
Molecular dynamics simulation of morphological instabilities in solid-fluid interfaces
B.V. Costaa, P.Z. Courab and O.N. Mesquitac
Departamento de Física ICEX, Universidade Federal de Minas Gerais
Caixa Postal 702, 30123-970, Belo Horizonte, MG, Brazil
a email@example.com; b firstname.lastname@example.org; c email@example.com
Received on 1 February, 2001. Revised version received on 23 May, 2001
We perform molecular dynamics simulations of directional growth of a binary alloy. We fix the temperature gradient, pulling speed, impurity concentration and only vary the impurity segregation coefficient. By changing the range of the Lennard Jones potential of impurity atoms as compared to the range of the potential of solvent atoms, the elastic energy cost causes a decreases of solubility of the impurity atoms into the matrix (solvent) solid phase and consequently a decrease in the impurity segregation coefficient. Within certain range of segregation coefficients, the growing interface is planar; below it, the interface becomes unstable and a cellular structure emerges.
Rapid growth of a solid-fluid interface under directional solidification conditions is a common method of investigating solidification morphologies as well as segregation during crystalization . Morphological instabilities in the solid-fluid interface have been studied for a long time . It is known that during growth, the crystal-fluid (melt, vapor or solution) interface, is not at thermodynamic equilibrium. The moving interface is a dynamical system, that can display a variety of dynamical instabilities and pattern formation [3, 4].
The rapid expansion of the use of high quality crystalline materials in optical and eletronic devices during the last decades has stimulated research, both theoretical and experimental, on dynamics of crystallization. A better understanding on solidification of metals and eutectic fibers are of unquestionable technological interest. Computer simulations have played an important role on the development and understanding of recent models of crystal growth [5-9]. A crystal can grow from the adjacent fluid (melt, vapor or solution) by different mechanisms, depending on the structure of the interface (rough or smooth), material purity, growth rates, temperature gradients and related factors. These requirements can be met in a controlled way in experiments of directional growth, where a sample in an appropriate furnace is submitted to a temperature gradient and pulled with a fixed speed towards the colder region of the furnace. For practical crystal growth, the sample can be cast into a quartz tube with chosen diameter and length. This is a 3-dimension Bridgman growth arrangement. However, detailed studies about dynamics of crystal growth have been conducted in very thin transparent samples (sandwiched between glass slides), such that the crystal-fluid interface can be visualized and recorded with the use of videomicroscopy techniques [10, 11]. The results of such experiments have been comparared with results of two-dimensional models of crystal growth. Our computer simulations are also carried out in two dimensions. With the use of molecular dynamics techniques we simulate the solidification of a two-component system consisting of solvent (atoms a) and solute (atoms b) interacting via Lennard-Jones(LJ) potentials : Fa,a, Fb,b and Fa,b=Fb,a which has been shown to describe very well the directional growth of a binary mixture [7, 8, 9]. the method is describe in details in section III. By tuning the parameters of the LJ potential, we can change the structure of the interface and the segregation coefficient such that a planar interface can go unstable.
In this paper, we simulate a binary system with a rough crystal-fluid interface and with several values of the LJ interaction potential as discussed below. This work is organized as follow. In section II we briefly discuss the directional growth of binary mixtures. In section III we discuss the simulation method we have used. In section IV we show our results and in section V we present our conclusions.
II Directional Growth of Binary Mixtures
A great deal of experimental observations in directional growth have been done in thin samples of transparent materials, where presumably the third dimension is not important, and the crystal-fluid interface can be followed in time, by using videomicroscopy with digital image analysis. The growth furnace consists of two metal blocks at controlled temperatures:one block with temperature above the sample melting temperature and the other below. A thin sample of the material to be studied is sandwiched between glass slides with spacers, whose thickness is in general of a few micra, to keep the system as close as possible of a 2d geometry, and guarantee that the transport of mass is mainly diffusive with no convection in the fluid phase. The sample is then put on top of the metal blocks, with good thermal contact. A temperature gradient appears along of the sample cell. The sample is then pulled towards the cold block in a very precise and controlled way by a pulling system and starts to solidify. After an initial transient the system achieves a steady-state where the interface position becomes fixed in the laboratory frame, the growth speed is the same, but opposite, as the pulling speed and the solute concentration profile becomes steady in the laboratory frame. Morphological instabilities of the planar interface are inhibited during directional growth of pure materials. However, for binary systems, depending on concentration of solute, temperature gradient and growth speed, morphological instabilities can occur: The planar interface becomes cellular and eventually dendritic. This is the so called Mullins-Sekerka instability (MS). In our simulations we are able to observe both regimes: planar and cellular interfaces. In this work we will focus our attention on the evolution of an interface where we clearly see segregation, transport of solute at the interface and morphological instabilities.
Our simulation is carried out using molecular dynamics with all particles interacting through the potential.
where fi,j is the LJ potential
The indexes i and j stands for particles in the positions i and j respectively, and 0 £ i,j £ N, where N is the total number of particles and i,j = | i-i |. A cutoff, rc = 2.5 saa is introduced in the potential in order to accelerate the simulation. If the force on a particle is found by summing contributions from all particles acting upon it, then this truncation limits the computational effort to an amount proportional to the total number of particles N. Of course this truncation introduces discontinuities both in the potential and force. To smooth this discontinuity we introduce a constant term f(rc). Another term (df/ dr)r=rc(r-rc) is introduced to remove the force discontinuity. Particles in our simulation move according the Newton's law, that generate a set of 2N coupled equations of motion which are solved by increasing forward in time the physical state of the system in small time-steps of size Dt = 0.02saa (ma / aa). The resulting equations are solved by using the Beeman's method of integration. In order to improve the method we use a Verlet and a celular table. The Verlet table consists of an address vector which contains the number and position of each particle inside a circle of radius rv = 3 saa. After some steps in time, the neighborhood of each particle changes, so that, we have to refresh the Verlet table. This refreshment process can take a long time. In order to make it shorter we divide the system in cells of size cx × cz = (3.5saa)2, such that in recalculating the Verlet table we have to search only in neighbor cells. Inittialy we distribute N = nx ×nz = 162 × 110 particles over the two dimensional surface Lx × Lz = 162 × 21/6 ×110 × 21/6. We assume periodic boundary conditions in the x direction. In the z direction we divide the system in two distinct regions, a solid and a fluid one. In the solid region particles stand initially in their equilibrium position in a total of 162 × 30 particles. On the fluid region the density is initially rf = 0.5s, giving a total of 162 × 80 particles, randomly distributed in a triangular lattice and slightly dislocated from their equilibrium position. We impose a temperature gradient, G, along the z direction using a velocity renormalization approach [7, 8, 9]. In one side of the system defined by 170saa £ z £ 151saa, temperature is fixed to T0 = 0. In another region 0 £ z £ 30 saa, temperature is fixed to Th = 0.7aa/kB, higher than the melting temperature Tf = 0.4 aa/kB of the pure material. The melting temperature was measured independently as described below. We let the system evolute for 1.2 × 105 steps in time of size Dt = 0.02 × saa × (ma / aa)1/2, which seems to be enough to equilibrate the system. Once the equilibrium is reached we start pulling the system in the +z direction, at a pulling velocity vp = 4 × 10-3 (aa/ma)1/2. The z = 0 layer works as a source of particles, mantaining a constant particle flow to the material. In order to obtain an estimate for Tf we did an independent simulation with 1.5 × 103 particles in a box of fixed dimensions and density 7 × 10-3saa, below the solid density. We obtained that Tf » 0.4 × aa/kB. The considered system consists of two different types of particles. One is the solvent ( a particles) and the solute (b particles). We define three types of interactions, solute-solute (b - b), solvent-solvent (a - a), solute-solvent (b - a). The initial solute concentration c0 is 5%. This concentration is defined as the ratio between the number of particles in the solute and the total number of particles in the solvent. To calculate the density along the growing crystal we use strips of size Dz = 20 saa. For simplicity, from now on we measure energy in units of aa, distance in units of saa, mass in units of ma so that temperature and time are measured in units of kB/aa and (ma s/aa)-1/2 respectively. The simulation is carried out for bb = 0.1 and ab = 0.5. Value of sab are 1.00, 1.05, 1.15, 1.20. We always use saa = sbb. Typically, in our simulations of directional growth, the solid phase density rs was 5.7 times larger than the fluid phase density rf. For a given pulling speed the time for the system to achieve steady-state, and consequently the simulation time, decreases as one increases the ratio between the solid and fluid densities. The reason is that the typical time for reaching steady-state is of the order of ~ D/vf2 = D/vp2 (rf/rs)2, where D is the diffusion coefficient and vf is the velocity of the fluid phase. Since D increases approximately linearly with decreasing density ratio, the computation time goes approximately linear with (rf/rs). We checked the velocity profile of the fluid phase and determined that solute transport is mainly diffusive for our simulations, with D = 1.3 saa(aa/ma)1/2
We studied the evolution of cellular instabilities during directional solidication for fixed temperature gradient G = 5 × 10-4 (aa/kBsaa) and fixed pulling speed vp = 4 × 10-3 (aa/ma)1/2, but for different values of the segregation coefficient k. From Mullins-Sekerka predictions , the interface is unstable if lT < lD, where lT is the thermal length and lD is the diffusion length. lT » and lD » » 60saa. Since L » 0.59 aa (latent heat) , one obtains that lT » 27saa. Therefore the interface is stable for k > 0.25. By decreasing k one can trigger the cellular instability. The value of k estimated from Fig. 1 is about 0.3. The interface is stable as predicted above. By increasing the value of sab as compared to saa, the impurity to be able to incorporate into the crystal has to deform its lattice. The high elastic energy cost causes a decrease of the impurity solubility into the crystal phase, resulting in a smaller k. From Figs. 2 and 3 we estimate k » 0.22 for sab = 1.05 saa; k » 0.04 for sab = 1.15 saa and k » 0.01 for sab = 1.20 saa, all values are in the unstable region. As k gets smaller, the crystal becomes cleaner and impurities concentrate into the grooves of the cellular structure.
|Figure 1. Final configuration of the simulated crystal after 8 ×104 time steps for sab=1.0. Grey and black circles represent solvent and solute atoms respectively. The solvent (grey) crystaline structure is hexagonal as expected for a Lennard-Jones solid. Pulling direction is along Z. Morphological instability is not observed. There is solute segregation at the solid-fluid interface, with k » 0.3.|
|Figure 2. From top to botton the plots show the surface's time evolution for sab = 1.05. The sequence starts at t = 1.6 × 104 following in time steps of size Dt = 4.0 × 104. Symbols are the same as in figure 1. The estimated k » 0.22 indicates that the interface is in the unstable region as observed.|
|Figure 3. Final configuration of the simulated crystal after 8 × 104 time steps. From top to botton sab = 1.05, 1.15 and 1.20. Symbols are the same as in figure 1. As sab increases, segragation increases (k decreases) facilitating the appearence of cellular instabilities.|
Our computer simulations of directional growth of a binary mixture clearly show: (1) The effect of impurity segregation during growth; (2) Segregation increases (k decreases) as sab becomes larger than saa; (3) For k smaller than 0.25 the cellular instability is triggered for the values of G, vp and c0 used. The impurities segregated at the interface concentrate into the grooves of the cellular structure resulting in impurity inclusions that cause large defects in the crystal matrix.
This work was partially supported by CNPq, FAPEMIG and FINEP. Numerical work was done in the LINUX parallel cluster at the Laboratório de Simulação Departamento de Física - UFMG.
 D. P. Woodruff, The solid-Liquid Interface (Cambridge University Press, London, 1973). [ Links ]
 W. W. Mullins and R. F. Sekerka, J. Appl. Phys. 34, 323 (1963); 35, 444 (1964). [ Links ]
 J.S. Langer, Rev. Mod. Phys. 52, 1 (1980). [ Links ]
 M. C. Cross and P. Hohenberg, Rev. Mod. Phys. 65, 851 (1993). [ Links ]
 J. D. Weeks and G. H. Gilmer, Adv. Chem. Phys. 40, 157 (1979). [ Links ]
 M.J.P. Nijmeijer and D. P. Landau, Comput Mater Sci. 1, 389 (1993). [ Links ]
 P. Z. Coura, O. N. Mesquita and B. V. Costa, Phys. Rev B 59, 3408 (1999). [ Links ]
 P. Z. Coura, O. N. Mesquita and B. V. Costa, Int. J. Mod. Phys. C 9, 857 (1998). [ Links ]
 B.V.Costa, P.Z. Coura and O.N. Mesquita in "Computer Simulation in Condensed Matter Physics XII", Springer Proceedings in Physics 85, Eds. D.P. Landau, S.P. Lewis and H.-B Schuetler (Springer Verlag, 2000)p. 191. [ Links ]
 J.D. Hunt, K.A. Jackson and H. Brown, Rev. Sci. Instrum. 37, 805 (1966). [ Links ]
 J. M. A. Figueired, M. B. L. Santos, L. O. Ladeira and O. N. Mesquita, Phys. Rev. Lett. 71, 4397 (1993). [ Links ]