Molecular dynamics simulation of morphological instabilities in solid-fluid interfaces

Molecular dynamics simulation of morphological instabilities in solid-fluid interfaces

B.V. Costa^a, P.Z. Coura^b and O.N. Mesquita^c

Departamento de Física ICEX, Universidade Federal de Minas Gerais

Caixa Postal 702, 30123-970, Belo Horizonte, MG, Brazil

[e-mail]; ^b [e-mail]; ^[e-mail]

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.

I Introduction

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 [1]. Morphological instabilities in the solid-fluid interface have been studied for a long time [2]. 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 : F_a,a, F_b,b and F_a,b=F_b,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.

III Simulation

Our simulation is carried out using molecular dynamics with all particles interacting through the potential.

where f

_i,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, r_c = 2.5 s_aa 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(r_c). Another term (df/ dr)r=r_c(r-r_c) 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.02s_aa (m_a / _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 r_v = 3 s_aa. 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 c_x × c_z = (3.5s_aa)², such that in recalculating the Verlet table we have to search only in neighbor cells. Inittialy we distribute N = n_x ×n_z = 162 × 110 particles over the two dimensional surface L_x × L_z = 162 × 2^1/6 ×110 × 2^1/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 r_f = 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 170s_aa£ z £ 151s_aa, temperature is fixed to T₀ = 0. In another region 0 £ z £ 30 s_aa, temperature is fixed to T_h = 0.7_aa/k_B, higher than the melting temperature T_f = 0.4 _aa/k_B of the pure material. The melting temperature was measured independently as described below. We let the system evolute for 1.2 × 10⁵ steps in time of size Dt = 0.02 × s_aa × (m_a / _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 v_p = 4 × 10^-3 (_aa/m_a)^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 T_f we did an independent simulation with 1.5 × 10³ particles in a box of fixed dimensions and density 7 × 10^-3s_aa, below the solid density. We obtained that T_f» 0.4 × _aa/k_B. 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 c₀ 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 s_aa. For simplicity, from now on we measure energy in units of _aa, distance in units of s_aa, mass in units of m_a so that temperature and time are measured in units of k_B/_aa and (m_as/

_aa)^-1/2 respectively. The simulation is carried out for _bb = 0.1 and _ab = 0.5. Value of s_ab are 1.00, 1.05, 1.15, 1.20. We always use s_aa = s_bb. Typically, in our simulations of directional growth, the solid phase density r_s was 5.7 times larger than the fluid phase density r_f. 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/v_f² = D/v_p² (r_f/r_s)², where D is the diffusion coefficient and v_f is the velocity of the fluid phase. Since D increases approximately linearly with decreasing density ratio, the computation time goes approximately linear with (r_f/r_s). We checked the velocity profile of the fluid phase and determined that solute transport is mainly diffusive for our simulations, with D = 1.3 s_aa(_aa/m_a)^1/2

IV Results

We studied the evolution of cellular instabilities during directional solidication for fixed temperature gradient G = 5 × 10^-4 (_aa/k_Bs_aa) and fixed pulling speed v_p = 4 × 10^-3 (_aa/m_a)^1/2, but for different values of the segregation coefficient k. From Mullins-Sekerka predictions [2], the interface is unstable if l_T < l_D, where l_T is the thermal length and l_D is the diffusion length. l_T» and l_D» » 60s_aa. Since L » 0.59 _aa (latent heat) [7], one obtains that l_T» 27s_aa. 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 s_ab as compared to s_aa, 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 s_ab = 1.05 s_aa; k » 0.04 for s_ab = 1.15 s_aa and k » 0.01 for s_ab = 1.20 s_aa, 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 ×10⁴ time steps for s_ab=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.

V Conclusions

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 s_ab becomes larger than s_aa; (3) For k smaller than 0.25 the cellular instability is triggered for the values of G, v_p and c₀ 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.

Acknowledgements

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.

• [1] D. P. Woodruff, The solid-Liquid Interface (Cambridge University Press, London, 1973).
• [2] W. W. Mullins and R. F. Sekerka, J. Appl. Phys. 34, 323 (1963); 35, 444 (1964).
• [3] J.S. Langer, Rev. Mod. Phys. 52, 1 (1980).
• [4] M. C. Cross and P. Hohenberg, Rev. Mod. Phys. 65, 851 (1993).
• [5] J. D. Weeks and G. H. Gilmer, Adv. Chem. Phys. 40, 157 (1979).
• [6] M.J.P. Nijmeijer and D. P. Landau, Comput Mater Sci. 1, 389 (1993).
• [7] P. Z. Coura, O. N. Mesquita and B. V. Costa, Phys. Rev B 59, 3408 (1999).
• [8] P. Z. Coura, O. N. Mesquita and B. V. Costa, Int. J. Mod. Phys. C 9, 857 (1998).
• [9] 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.
• [10] J.D. Hunt, K.A. Jackson and H. Brown, Rev. Sci. Instrum. 37, 805 (1966).
• [11] J. M. A. Figueired, M. B. L. Santos, L. O. Ladeira and O. N. Mesquita, Phys. Rev. Lett. 71, 4397 (1993).

Publication Dates

• Publication in this collection
  26 Mar 2002
• Date of issue
  Sept 2001

History

• Received
  01 Feb 2001
• Reviewed
  23 May 2001