- Email: [email protected]
A molecular dynamics study of sintering between nanoparticles
Computational Materials Science 45 (2009) 247–256
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
A molecular dynamics study of sintering between nanoparticles Lifeng Ding a, Ruslan L. Davidchack b, Jingzhe Pan a,* a b
Department of Engineering, University of Leicester, Leicester LE1 7RH, UK Department of Mathematics, University of Leicester, Leicester LE1 7RH, UK
a r t i c l e
i n f o
Article history: Received 10 July 2008 Received in revised form 16 September 2008 Accepted 18 September 2008 Available online 8 November 2008 PACS: 61.46.Df 81.20.Ev Keywords: Molecular dynamics Nanoparticles Sintering Modelling
a b s t r a c t The paper presents a molecular dynamics study on the interactions between nanoparticles at elevated temperatures. The emphasis is on the comparison between the molecular dynamics model and the continuum model using solid state physics. It is shown that the continuum model is unable to capture the sintering behaviour of nanoparticles. This is not because the continuum theory does not apply at the nano-scale but because the nanoparticles behave in so many different scenarios of the continuum theory that a meaningful model has to predict these scenarios, using the molecular dynamics for example. In the MD simulation, it is observed that the particles reorient their crystalline orientations at the beginning of the sintering and form different types of ‘‘necks” between different particles. This leads to different mechanisms of matter redistribution at the different necks. It is also observed that the particles can switch the mechanism of matter transportation half-way through the sintering process. It would be very difficult, if not impossible, to handle these complexities using the continuum model. However assuming the right scenario, the continuum theory does agree with the MD simulation for particles consisting of just a few thousands atoms. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction Nanoparticles are of a growing interest to modern industry and technology. Many nanostructured materials have shown exciting properties and applications [1,2]. The main challenge is to avoid coarsening during powder or material processing so that the nanostructure can be maintained. Sintering is the simplest and most cost effective way to fabricate solid material from nanopowders. A lot of efforts have been made to understand the process of nanoparticles sintering [1–3]. Although there are cases where considerable progress has been achieved, [3] generally the problem of grain coarsening is still a major concern for material processing. To understand the fundamental mechanisms of nanoparticle sintering, computer simulations, particularly atomistic simulations, offer a convenient and practical way to investigate such small scale phenomena. There have been many attempts to model sintering of nanoparticles using MD simulation. Zhu and Averback [4] performed molecular dynamics simulations on copper nanoparticles and they suggested that the high shear stresses in the neck could lead to rapid nanoparticle sintering, therefore the primary mechanism responsible for sintering should be plastic deformation. Zeng [5] also studied copper nanoparticle sintering via MD simulation. They noticed that surface and grain-boundary diffusion are two main mechanisms, although dislocations are present at the early * Corresponding author. E-mail address: [email protected] (J. Pan). 0927-0256/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2008.09.021
stage of the sintering. Raut and Bhagat [6] questioned the importance of plastic deformation and surface diffusion during nanoparticle sintering. Moreover, Tsuruta [7] carried out large scale MD simulation of both silicon nitride crystal and amorphous nanoclusters sintering. They observed relative rotation at early sintering stage and concluded that sintering was controlled by surface diffusion. In many ways the application of molecular dynamics to sintering modelling is still at its infancy. Most of the reported simulations of sintering were people’s first attempts. The aim was to explore the feasibility of the MD modelling instead of providing accurate predictions. Recently, Ch’ng and Pan [8] observed some interesting phenomenon in their computer simulation of nanoparticles sintering using a continuum model. They found that for nanoparticles a small particle placed between two larger neighbours may possess a much higher resistance to coarsening than micron-sized particles. This numerical finding would explain the so called ‘‘frustration” effect which has been observed by several groups, i.e. both oxide and metal nanomaterials show a strong resistance to grain coarsening [9–13]. However it may be argued that the continuum model used by Ch’ng and Pan [8] is invalid for nanoparticles because such particles should be considered as discrete assemblies of atoms. The purpose of this paper is to provide a careful comparison between the molecular dynamics model and the continuum model for nanoparticle sintering in order to gain an insight into the difference between the two models. To avoid introducing too many factors in the comparison, we limit ourselves to two-dimensional problems and to the simplest atomic
248
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
2. Molecular dynamics and continuum models for sintering
uration corresponds to the (1 1 1) oriented layer of a three dimensional face-centered-cubic (fcc) crystal. Before the three particles are placed in contact with each other, each particle is equilibrated for 10,000 t* at constant temperature. The constant temperature is realized by introducing a nonholonomic constraint into the equations of motion in order to keep the kinetic energy constant [18]. The constraint equation to ensure constant temperature is
2.1. The molecular dynamics model
NEk ¼
We use the classical molecular dynamics with a simple twodimensional Lennard-Jones model. The pairwise Lennard–Jones (LJ) potential [16] between two atoms at a distance r from each other is defined as follows:
where the average kinetic energy of an atom, Ek, is related to the system temperature: Ek = 3/2kBT. The constrained equation of motion is
12 6 r r UðrÞ ¼ 4e r r
r€i ¼
potential – the Lennard-Jones model. For the continuum model, Herring’s theory [14] for solid state diffusion is adopted and a rigorous numerical scheme developed by Pan and Cocks [15] is used. By keeping both models simple, there is strong confidence in the numerical solution schemes and therefore the comparison is built on a very solid basis.
ð1Þ
ð2Þ
where
/ðxÞ ¼ 2x3 3x2 þ 1
ð3Þ
for 0 < x < 1.
xðrÞ ¼
r rm rc rm
1 fi þ ar_i m
ð5Þ
ð6Þ
and since E_ k ¼ 0, or equivalently
where e is the depth of the potential well and r is the range of the potential. All measurements are expressed in units derived from the three basic units of energy e, length r and mass of an atom m. For example, the unit of time is t* = r(m/e)1/2, the unit of temperature is e/kB, where kB is the Boltzmann constant, and so on. In order to speed up the simulation and avoid self-interaction of particles in simulation with periodic boundary conditions, the long ‘‘tail” of LJ potential interaction is truncated by introducing a differentiable function /ðxðrÞÞ, which is equal to 1 for x 6 0 and 0 for x P 1. Then the truncated potential is given by
U tr ðrÞ ¼ UðrÞ/ðxðrÞÞ
N 1 X m r_ 2 ; 2 i¼1 i
ð4Þ
where rc = 2.5r is the cut-off radius and we choose rm = 0.95 rc. The classical Newton’s equations of motion are integrated using the leap-frog method [17]. The time step of 0.02 t* is used to guarantee stable numerical integration of the equations of motion. An example of the starting configurations is shown in Fig. 1. Each bulk atom has 6 nearest neighbour atoms. This two dimensional config-
Fig. 1. The initial configuration of the molecular dynamics model of three particles in contact with each other.
X
€r i r_ i ¼ 0
ð7Þ
i
it follows that the value of the Lagrange multiplier a is
P
r_ fi : _2 i ri
i i a¼ P
m
ð8Þ
The equilibrated particles are rotated with respect to each other in order to investigate the effect of crystallographic orientation on the process of neck formation. The particles are then placed next to each other such that the distance between the nearest atoms from two contacting particles is about 1.7r. In order to prevent the loss of any vapour atoms from the system, periodic boundary conditions are employed, so that a vapor atom leaving the simulation cell through the boundary is re-introduced into the cell with the same velocity at the opposite boundary. The size of the simulation cell is large enough to prevent atoms on the particles from interacting with their neighbor images except for the vapor atoms. Following Ch’ng and Pan, [8] three different particle arrangements are used in order to study the particle interactions. In configuration I, a small particle whose radius is 16r is placed next to a bigger particle which is twice the size of the small one. In configuration II, two same size particles are placed in contact. Their radius is 32r. In configuration III one small particle is placed along a straight line between two big particles (see Fig. 2). The larger par-
Fig. 2. The three different configurations of particles: (I) two particles of different size; (II) two particles with same size; (III) one small particle in the middle of two large particles.
249
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
ticle consists of around 3000 atoms while the smaller one consists of around 1000 atoms. The simulation results, including atom coordinates, velocities, and potential energy are recorded every 500 time steps. Thus the simulation progress could be observed directly from dynamic movie obtained by plotting all the atom coordinates. When two particles join together they create a ‘‘neck”. We monitor neck growth with time by measuring its width, defined as the distance between the outer atoms at the narrowest place. We also measure the particle center to centre distance determined by measuring the distance between mass centers of the core parts of the particles with a radius about 10r which remain crystalline throughout the simulation and thus provide a faithful measure of the relative particle displacement during sintering. 2.2. The continuum model The solid state diffusion theory originally due to Herring [19] is adopted for the continuum model. According to this theory, the diffusive flux j, defined as volume of matter passing across unit area perpendicular to the flux direction per unit time, is assumed to be linearly dependent on the gradient of the chemical potential l of the diffusing species
j¼
Dd ol ; kB T os
ð9Þ
where D is diffusivity, d is the thickness of the layer through which the material diffuses, kB is the Boltzmann constant, T is the absolute temperature and s is the local coordinate along the diffusion path. Along a grain-boundary, the gradient of the atomic chemical potential is induced by the gradient of stress r acting normal to the boundary [20]
l ¼ Xr;
ð11Þ
ð12Þ
ð13Þ
where cs is the specific surface energy and j is the principal curvature of the surface. It is positive for a concave free surface. So we have
js ¼
Ds ds Xcs oj kB T os
In a MD simulation, it is inevitable to have certain amount of vapour atoms. When a vapour atom flies across the periodic boundary, the angular momentum of the system changes causing the system to move and rotate rigidly. Sintering occurs by atoms migrating mainly along the particles surface. Atoms at the core of a particle only vibrate thermally and move rigidly. The translational and rotational rigid motions have to be eliminated, as shown in Fig. 3, during post-processing in order to filter out information such as diffusion coefficient, neck size, centre to centre distance, etc. This numerical detail is presented here to demonstrate the feature of the MD model that does not exist in the continuum model. Considering a rigid core of a particle which contains n atoms, after certain time t, the rigid core rotates by an angle a about its mass centre. Inside the rigid core, the atom coordinates before and after rotation are denoted by ri and r0i . The angle a of the rigid rotation of the particle core can be calculated as
where di = rircm is the atom coordinate with respect to the mass centre of the particle core. The rigid motion can then be subtracted from the atomic coordinates according to the following coordinate transformation equations:
Along a free surface, the gradient of the atomic chemical potential is induced by the gradient of the free surface curvature [20]
l ¼ Xcs j;
3.1. Subtraction of rigid motion of the entire system
ð10Þ
where Db is the grain-boundary diffusivity and db is the grainboundary thickness. As matter is deposited onto (or removed from) a grain-boundary, the grains on either side of the boundary move with a normal velocity vb. Matter conservation requires that
ojb þ mb ¼ 0: os
3. Molecular dynamics simulation of nanoparticles sintering
P 0 jdi di j sina ¼ Pi 0 0 ; i jdi jjdi j
where X is the atomic volume, therefore
Db db X or jb ¼ ; kB T os
numerical scheme to solve the coupled diffusion problem which is used in this study. When comparing with the molecular dynamics model it is important to bear in mind that a continuum model has to assume a ‘‘mechanism of matter redistribution”, grain-boundary diffusion for example, while the MD model does not make such an assumption. It is also useful to point out that Herring’s theory assumes thermodynamic equilibrium at every step of the microstructure evolution while the molecular dynamics model has no such limitation.
ð16Þ
x00i ¼ ðy0i ycm Þcosa þ ðy0i ycm Þsina þ xcm ;
ð17Þ
y00i ¼ ðx0i xcm Þsina þ ðy0i ycm Þcosa þ ycm ;
ð18Þ
r 0i
ðx0i ; y0i Þ
r 00i
ðx00i ; y00i Þ
where ¼ and ¼ and after rotation subtraction.
denote atom coordinates before
3.2. Volume fraction of atoms at the particle surface From the point of view of the continuum model, the driving force for particle sintering is the reduction in the free energy asso-
ð14Þ
in which Ds is the free surface diffusivity. As matter is deposited onto (or removed from) a free surface, it migrates with a normal velocity vs. Matter conservation requires that
ojs þ ms ¼ 0 os
ð15Þ
with ms assumed to be positive when matter is added to the surface. The grain-boundary and surface diffusion processes are coupled at the triple contact point where the grain-boundary meets the free surface. Furthermore, the grain-boundary stress has to satisfy an equilibrium condition. Pan and Cocks [15] developed a simple
Fig. 3. Mathematical subtraction of rigid rotation during post processing.
250
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
Fig. 4. Number of atoms as a function of its potential energy indicating their locations at different parts of the particle during sintering.
ciated with the particle surface. In a continuum model, the surface energy is simply proportional to the total surface area. It is therefore instructive to monitor the volume fraction of atoms on the particle surface in the MD simulation to check against the continuum concept. In order to do this, it is necessary to distinguish atoms from different parts of the system, i.e. on the surface, in the bulk, on the grain-boundary or in the vapour. This can be achieved by examining the potential energy of interaction of an atom with other atoms. We record the average potential energy of each atom over 500 time steps. The histogram in Fig. 4a shows the number of atoms as a function of their potential energy. It can be clearly observed that the number of atoms peaks around several values of the potential energy, each representing a group of atoms with a particular number of nearest neighbour atoms. The numbers of neighbouring atoms are indicated by colours as shown in the fig-
Fig. 5. Volume fraction of surface atoms as a function of sintering time.
ure. The peaks from left to right represent atoms with 6, 5, 4, 3, 0 neighbours respectively. The left prime peak represents the majority atoms in the bulk which have all six neighbours arranged in a correct crystal configuration. The sub-peak under the prime peak represents atoms on the grain-boundary (i.e. distorted crystal configuration). All the other peaks except for the vapour atoms represent atoms on the particle surface. This identification of different atoms can be used to monitor mass transport during the sintering process. Fig. 5 shows the volume fraction of the surface atoms (with potential energy larger than 5.0) as a function of time for the two particles sintering in configuration II. It can be seen clearly that the total number of surface atoms decreases as sintering proceeds, which is consistent with the continuum concept that the driving force for sintering is the reduction in the total surface energy.
Fig. 6. Particles with initially aligned crystal orientations maintain the alignment during sintering; only half of the particles are shown to reveal their relative rotation.
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
Fig. 7. Particles with initial misaligned crystal orientations. The small particle reorients itself at the very beginning of the sintering process; only half of the particles are shown to reveal their relative rotation.
3.3. Effect of crystalline orientation of the particles The information that the continuum model does not contain is the crystalline orientation of the particles. It turns out that this information in the MD model leads to a sintering behaviour which the continuum model is unable to capture. Fig. 6 shows two particles with perfectly aligned initial crystalline orientation while Fig. 7 shows three particles of different initial crystalline orientations. Only half of the particles are shown to reveal the relative rotation between the particles during sintering. The MD simulation shows that the particles with perfectly aligned initial crystalline orientation stay aligned during sintering. For particles of different initial crystalline orientations, a quick adjustment occurs in the
251
crystalline orientations at the very beginning of the sintering progress. The particles then remain in the adjusted orientations throughout the rest of the sintering process. This phenomenon has been observed by other researchers both experimentally [21] and in MD simulation [5]. The reorientation is believed to be driven by the minimization of the grain-boundary energy. However the consequence of the reorientation on sintering has not been fully realised in previous studies. In our MD simulation, it is observed that the reorientation can eliminate a grain-boundary from the system as shown in Fig. 8 and significantly retard the sintering process. In some other cases, a stable grain-boundary is formed after the initial adjustment, as shown in Fig. 9. Our simulations show that three factors can help the formation of a stable grain-boundary: (a) a large angle of initial misalignment between the particles, (b) big particle size and (c) higher sintering temperature. In the two dimensional MD model, the maximum misalignment angle is 30 degrees because of the 6-fold symmetry of the crystalline structure. The initial misalignment has to be large enough in order to resist the elimination of the grain-boundary by the particle reorientation. Our simulation suggests that the grain boundaries are more likely to form between two large nanoparticles if the initial misalignment angle is larger then 20 degrees. Furthermore the MD simulation shows that a grain-boundary is more likely to form between the two large nanoparticles in configuration II (see Fig. 2), but rarely form between two smaller nanoparticles. The smaller particles are relatively easy to reorient themselves to eliminate the grain-boundary. A typical evolution of the system in configuration III, is that the small particle quickly aligns with one of the large particles to form a single crystal while a grain-boundary is formed between the small particle and the other large one. 3.4. Sintering of particles The primary conclusion of the continuum modelling by Ch’ng and Pan [8] is that nanoparticles arranged in the configuration III shown in Fig. 2 seem to possess certain resistance to particle coarsening while the configuration I does not have such resistance. However using the continuum model for particles consisting of only a few thousands of atoms pushes its validity to the limit. One major feature
Fig. 8. Two nanoparticles sinter together to form a single crystal, the lines mark the crystalline orientations before and after the reorientation.
252
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
Fig. 9. Two nanoparticles sinter together to form a grain-boundary; the lines mark the crystalline orientations before and after the reorientation.
Fig. 10. MD simulation of one small and two large particles sintering at T = 0.26 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
Fig. 11. MD simulation of a small and a large particle sintering at T = 0.26 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
that the continuum model is unable to capture is the elimination of the grain-boundary as shown in Section 3.3. If a grain-boundary is eliminated, the grain-boundary diffusion process ceases to exist, leaving only surface diffusion and vapour evaporation and condensation as the mechanisms for matter redistribution. This phenomenon has a profound effect on the sintering behaviour of the particles.
253
The melting temperature of the two-dimensional Lennard– Jones crystal is Tm = 0.415 e/kB [22]. We performed MD simulations at two different temperatures: T = 0.19 e/kB and T = 0.26 e/kB, which are 46% and 63% of the melting temperature, respectively. It is clear that sintering occurs much faster at the higher temperature. For the two identical particles of initial diameter
Fig. 12. MD simulation of two large particles sintering at T = 0.26 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
Fig. 13. MD simulation of one small and two large particles sintering at T = 0.19 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
Fig. 14. MD simulation of a small and a large particle sintering at T = 0.19 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
254
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
Fig. 15. MD simulation of two large particles sintering at T = 0.19 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
Fig. 16. MD simulation of two large particles sintering with the grain-boundary at T = 0.26 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
70r (configuration II in Fig. 2) without a grain-boundary, the neck size reached 30 r at T = 0.26 e/kB after 25000 t*, but only reached 20 r at T = 0.19 e/kB by the same time. If a grain-boundary was formed, then the neck size reached 45 r at T = 0.26 e/kB after 25000 t*. Figs.10–12 show the MD simulations for the three particle configurations sintering at T = 0.26 e/kB. We have picked three cases for which the grain-boundaries were eliminated at the beginning of the simulation. In all figures showing MD simulations the following colouring scheme is used: the bulk atoms are always coloured as blue, any atom that has been on the surface or the grain-boundary of the particle is coloured as red,1 and any atom that has been in the vapour is coloured as cyan. The colouring scheme therefore traces the history rather than the current state of the atoms, providing a rough picture of the matter redistribution during the simulation run. In Fig. 10 one small particle is placed between two larger particles. At the very beginning, the particles quickly joined together and formed two necks. Although there was some time lag between the formations of the two necks, by 100 t* both necks were fully filled by surface atoms (shown in red) and the blue bulk part of the small particle has reduced in size. At the later stage, both surface diffusion and vapour condensation contributed to the growth of the neck. Fig. 11 shows one small particle and one large particle in contact with each other. Initially, the quick neck formation was sim-
1 For interpretation of the references to colors, the reader is referred to the web version of this article.
ilar to the case in Fig. 10. The neck was filled by atoms mostly coming from the small particle through surface diffusion. In the later stage of the sintering, evaporation and condensation dominated the matter redistribution. After the neck has been filled up the top of the small particle started to collapse and the whole system evolved toward a circular shape. Fig. 12 shows two large particles of the same size in contact with each other. Again the sintering process started with a quick neck formation followed by neck growth first by surface diffusion and then by vapour condensation. If one uses the continuum model, such behaviour of switching from one mechanism of matter transportation to another half-way through the sintering process would be very difficult, if not impossible, to handle. These three cases therefore highlight the difference between the MD model and the continuum model. The locations of the centres of mass of the particle cores were monitored during the simulation. In all the three cases, there was no perceptible change of the distance between the particles. This indicated the absence of grain-boundary diffusion and that any bulk deformation of the particle was negligible. The simulations were then repeated at T = 0.19 e/kB at which evaporation is much less important. Figs. 13–15 show these cases. It is obvious from these figures that at lower sintering temperatures surface diffusion dominates the matter redistribution. In all the above cases there was no grain-boundary between the particles in contact. Fig. 16 shows the case of two particles of the same size sintering at T = 0.26 e/kB for which a grainboundary is formed between the particles. By t = 100 t*, a small
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
Fig. 17. Centre to centre distance as a function of time for the case shown in Fig. 16.
255
Fig. 17 shows the centre to centre distance between the two particles as a function of time. It is clear that the particles are approaching each other during the sintering process confirming that atoms have been taken away from the grain-boundary by the grain-boundary diffusion. This can be compared with the case of Figs. 12 and 15 for which there was no grain-boundary between the two particles and the distance between particles remained unchanged. Fig. 18 shows the neck size as a function of time for the cases showing in Figs. 12, 15 and 16. It can be observed from the figure that grain-boundary diffusion promotes the necksize growth. Faster surface diffusion plus evaporation condensation at higher temperatures are also among the factors promoting faster necksize growth. Fig. 19 shows the MD simulation of one small particle in contact with two big neighbours sintering at T = 0.26 e/kB. This is perhaps the most interesting case which demonstrates the difference between the MD model and the continuum model. By t = 100 t* a grain-boundary was formed at the upper neck but the initial grain-boundary at the lower neck has been eliminated. Then two very different behaviours of the neck growth can be observed. The upper neck grows by grain-boundary diffusion while the lower one grows by surface diffusion. The grain-boundary also migrates into the small particle. It is also obvious that the existence of a grain-boundary accelerates significantly the neck growth at the upper neck. A continuum model on the other hand cannot possibly predict such a behaviour as the mechanism of matter transportation has to be assumed rather than predicted by such a model. 4. Direct comparison between MD and continuum models
Fig. 18. Comparison of neck size growth as a function of time for the cases shown in Figs. 12, 15 and 16.
neck made up of surface atoms has been formed. As sintering proceeds, the neck is filled up dominantly by the red atoms, meaning that these atoms have all been at the particle surface or the grain boundary at some stage of the sintering. The grain-boundary was migrating backward and forward between the two particles.
The most interesting outcome of the MD simulations presented so far is that asking if a MD model agrees with a continuum model is the wrong question. For nanoparticles, the sintering behaviour observed in a MD simulation is beyond the predictability by the continuum model. If the assumed mechanism of matter redistribution in the continuum model did occur in a MD simulation, then it would be interesting to see if the two models can predict the same speed of neck growth. Fig. 20 compares the neck sizes as functions of time predicted by the continuum and MD models respectively. The case is for two particles of the same size of R = 70 r sintering at T = 0.19 e/kB which is controlled by surface diffusion. Logarithmic scales are used for both axes. The straight line was obtained from the continuum model while the discrete dots were obtained from several independent MD simulations of the system initialised within the same macroscopic states (i.e. same atom positions and temperature) but different microscopic states (i.e. randomly
Fig. 19. MD simulation of one small and two large particles sintering at T = 0.26 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
256
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
uum theory of solid state physics does not apply at this scale, but because nanoparticles behave in so many different scenarios of the continuum theory that one has to predict these scenarios using the molecular dynamics model. In our MD simulation it has been observed that the particles may reorient themselves to match their crystalline orientations at the beginning of the sintering, and thus form different types of necks between different particles (either with or without the grain-boundary). This leads to different mechanisms of matter redistribution at the different necks. It has also been observed that the particles switch the mechanism of matter transportation half-way through the sintering process. None of these can be handled by the continuum model. However assuming the right scenario, the continuum theory does agree with the MD simulation for particles consisting of just a few thousand atoms. Fig. 20. Neck size divided by initial particle radius as a function of time predicted by the MD model (discrete dots) and continuum model (solid line) for the sintering of two particles of the same size, each consisting of 5000 atoms at T = 0.19 e/kB. Surface diffusion is the dominant mechanism for matter redistribution. Logarithmic scales are used for both axes.
Acknowledgements This work was supported by the EPSRC research grant S97996 which is gratefully acknowledged. One of the authors (RLD) did part of the work during his study leave granted by the University of Leicester. The computations were performed on the University of Leicester Mathematical Modelling Centre’s cluster, which was purchased through the EPSRC strategic equipment initiative. References
Fig. 21. Neck size divided by initial particle radius as a function of time predicted by the MD model (discrete dots) and continuum model (solid line) for the sintering of two particles of the same size, each consisting of 5000 atoms, at T = 0.26 e/kB. Grain-boundary diffusion is the dominant mechanism for matter redistribution. Logarithmic scales are used for both axes.
chosen initial velocities of individual atoms). Fig. 21 shows similar comparison where grain-boundary diffusion is the controlling mechanism of matter redistribution. The sintering temperature was set as T = 0.26 e/kB. All the other parameters are the same as that in Fig. 20. The MD prediction is smoother in Fig. 21 than Fig. 20 because the neck size was estimated manually from a single MD simulation. It is remarkable to see that the continuum model agrees with the MD model for particles consisting of only around 10,000 atoms. 5. Conclusions For nanoparticles, it is inappropriate to use continuum models to predict their sintering behaviour. This is not because the contin-
[1] R. Birringer, Materials Science and Engineering, A: Structural Materials: Properties, Microstructure and Processing A117 (1989) 33–43. [2] H. Gleiter, Progress in Materials Science 33 (1989) 223–315. [3] I.W. Chen, X.H. Wang, Nature (London) 404 (2000) 168–171. [4] H. Zhu, R.S. Averback, Materials Science and Engineering, A: Structural Materials: Properties, Microstructure and Processing A204 (1995) 96–100. [5] P. Zeng, S. Zajac, P.C. Clapp, J.A. Rifkin, Materials Science and Engineering, A: Structural Materials: Properties, Microstructure and Processing A252 (1998) 301–306. [6] J.S. Raut, R.B. Bhagat, K.A. Fichthorn, Nanostructured Materials 10 (1998) 837– 851. [7] K. Tsuruta, A. Omeltchenko, R.K. Kalia, P. Vashishta, Europhysics Letters 33 (1996) 441–446. [8] H.N. Ch’ng, A.C.F.C.J. Pan, Acta Materialia 55 (2007) 813–824. [9] J.R. Groza, R.J. Dowding, Nanostructured Materials 7 (1996) 749–768. [10] M.J. Mayo, D.C. Hague, D.J. Chen, Materials Science and Engineering, A: Structural Materials: Properties, Microstructure and Processing A166 (1993) 145–159. [11] R.W. Siegel, S. Ramasamy, H. Hahn, Z. Li, T. Lu, R. Gronsky, Journal of Materials Research 3 (1988) 1367–1372. [12] J.A. Eastman, Y.X. Liao, A. Naravanasamy, R.W. Siegel, Materials Research Society Symposium Proceedings 155 (1989) 255–266. [13] R.W. Siegel, Materials Research Society Symposium Proceedings 196 (1990) 59–70. [14] C. Herring, Journal of Applied Physics 21 (1950) 301–303. [15] A.C.F.C.J. Pan, Acta Metallurgica Material 43 (1995) 1395–1406. [16] J.E.C. Lennard-Jones, Proceedings of the Physical Society 43 (1931) 461. [17] J.A. Snyman, Applied Mathematical Modelling 6 (1982) 449. [18] W.G. Hoover, A.J.C. Ladd, B. Moran, Physical Review Letters 48 (1982) 1818. [19] F. Herring, Surface Tension as a Motivation for Sintering, McGraw- Hill, New York, 1955. [20] C. Herring, The Physics of Powder Metallurgy, McGraw-Hill, New York, 1951. [21] J. Rankin, B.W. Sheldon, Materials Science and Engineering A: Structural Materials: Properties, Microstructure and Processing A204 (1995) 48–53. [22] J.A. Barker, D. Henderson, F.F. Abraham, Physica A: Statistical Mechanics and Its Applications (Amsterdam, Netherlands) 106A (1981) 226–238.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
A molecular dynamics study of sintering between nanoparticles Lifeng Ding a, Ruslan L. Davidchack b, Jingzhe Pan a,* a b
Department of Engineering, University of Leicester, Leicester LE1 7RH, UK Department of Mathematics, University of Leicester, Leicester LE1 7RH, UK
a r t i c l e
i n f o
Article history: Received 10 July 2008 Received in revised form 16 September 2008 Accepted 18 September 2008 Available online 8 November 2008 PACS: 61.46.Df 81.20.Ev Keywords: Molecular dynamics Nanoparticles Sintering Modelling
a b s t r a c t The paper presents a molecular dynamics study on the interactions between nanoparticles at elevated temperatures. The emphasis is on the comparison between the molecular dynamics model and the continuum model using solid state physics. It is shown that the continuum model is unable to capture the sintering behaviour of nanoparticles. This is not because the continuum theory does not apply at the nano-scale but because the nanoparticles behave in so many different scenarios of the continuum theory that a meaningful model has to predict these scenarios, using the molecular dynamics for example. In the MD simulation, it is observed that the particles reorient their crystalline orientations at the beginning of the sintering and form different types of ‘‘necks” between different particles. This leads to different mechanisms of matter redistribution at the different necks. It is also observed that the particles can switch the mechanism of matter transportation half-way through the sintering process. It would be very difficult, if not impossible, to handle these complexities using the continuum model. However assuming the right scenario, the continuum theory does agree with the MD simulation for particles consisting of just a few thousands atoms. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction Nanoparticles are of a growing interest to modern industry and technology. Many nanostructured materials have shown exciting properties and applications [1,2]. The main challenge is to avoid coarsening during powder or material processing so that the nanostructure can be maintained. Sintering is the simplest and most cost effective way to fabricate solid material from nanopowders. A lot of efforts have been made to understand the process of nanoparticles sintering [1–3]. Although there are cases where considerable progress has been achieved, [3] generally the problem of grain coarsening is still a major concern for material processing. To understand the fundamental mechanisms of nanoparticle sintering, computer simulations, particularly atomistic simulations, offer a convenient and practical way to investigate such small scale phenomena. There have been many attempts to model sintering of nanoparticles using MD simulation. Zhu and Averback [4] performed molecular dynamics simulations on copper nanoparticles and they suggested that the high shear stresses in the neck could lead to rapid nanoparticle sintering, therefore the primary mechanism responsible for sintering should be plastic deformation. Zeng [5] also studied copper nanoparticle sintering via MD simulation. They noticed that surface and grain-boundary diffusion are two main mechanisms, although dislocations are present at the early * Corresponding author. E-mail address: [email protected] (J. Pan). 0927-0256/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2008.09.021
stage of the sintering. Raut and Bhagat [6] questioned the importance of plastic deformation and surface diffusion during nanoparticle sintering. Moreover, Tsuruta [7] carried out large scale MD simulation of both silicon nitride crystal and amorphous nanoclusters sintering. They observed relative rotation at early sintering stage and concluded that sintering was controlled by surface diffusion. In many ways the application of molecular dynamics to sintering modelling is still at its infancy. Most of the reported simulations of sintering were people’s first attempts. The aim was to explore the feasibility of the MD modelling instead of providing accurate predictions. Recently, Ch’ng and Pan [8] observed some interesting phenomenon in their computer simulation of nanoparticles sintering using a continuum model. They found that for nanoparticles a small particle placed between two larger neighbours may possess a much higher resistance to coarsening than micron-sized particles. This numerical finding would explain the so called ‘‘frustration” effect which has been observed by several groups, i.e. both oxide and metal nanomaterials show a strong resistance to grain coarsening [9–13]. However it may be argued that the continuum model used by Ch’ng and Pan [8] is invalid for nanoparticles because such particles should be considered as discrete assemblies of atoms. The purpose of this paper is to provide a careful comparison between the molecular dynamics model and the continuum model for nanoparticle sintering in order to gain an insight into the difference between the two models. To avoid introducing too many factors in the comparison, we limit ourselves to two-dimensional problems and to the simplest atomic
248
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
2. Molecular dynamics and continuum models for sintering
uration corresponds to the (1 1 1) oriented layer of a three dimensional face-centered-cubic (fcc) crystal. Before the three particles are placed in contact with each other, each particle is equilibrated for 10,000 t* at constant temperature. The constant temperature is realized by introducing a nonholonomic constraint into the equations of motion in order to keep the kinetic energy constant [18]. The constraint equation to ensure constant temperature is
2.1. The molecular dynamics model
NEk ¼
We use the classical molecular dynamics with a simple twodimensional Lennard-Jones model. The pairwise Lennard–Jones (LJ) potential [16] between two atoms at a distance r from each other is defined as follows:
where the average kinetic energy of an atom, Ek, is related to the system temperature: Ek = 3/2kBT. The constrained equation of motion is
12 6 r r UðrÞ ¼ 4e r r
r€i ¼
potential – the Lennard-Jones model. For the continuum model, Herring’s theory [14] for solid state diffusion is adopted and a rigorous numerical scheme developed by Pan and Cocks [15] is used. By keeping both models simple, there is strong confidence in the numerical solution schemes and therefore the comparison is built on a very solid basis.
ð1Þ
ð2Þ
where
/ðxÞ ¼ 2x3 3x2 þ 1
ð3Þ
for 0 < x < 1.
xðrÞ ¼
r rm rc rm
1 fi þ ar_i m
ð5Þ
ð6Þ
and since E_ k ¼ 0, or equivalently
where e is the depth of the potential well and r is the range of the potential. All measurements are expressed in units derived from the three basic units of energy e, length r and mass of an atom m. For example, the unit of time is t* = r(m/e)1/2, the unit of temperature is e/kB, where kB is the Boltzmann constant, and so on. In order to speed up the simulation and avoid self-interaction of particles in simulation with periodic boundary conditions, the long ‘‘tail” of LJ potential interaction is truncated by introducing a differentiable function /ðxðrÞÞ, which is equal to 1 for x 6 0 and 0 for x P 1. Then the truncated potential is given by
U tr ðrÞ ¼ UðrÞ/ðxðrÞÞ
N 1 X m r_ 2 ; 2 i¼1 i
ð4Þ
where rc = 2.5r is the cut-off radius and we choose rm = 0.95 rc. The classical Newton’s equations of motion are integrated using the leap-frog method [17]. The time step of 0.02 t* is used to guarantee stable numerical integration of the equations of motion. An example of the starting configurations is shown in Fig. 1. Each bulk atom has 6 nearest neighbour atoms. This two dimensional config-
Fig. 1. The initial configuration of the molecular dynamics model of three particles in contact with each other.
X
€r i r_ i ¼ 0
ð7Þ
i
it follows that the value of the Lagrange multiplier a is
P
r_ fi : _2 i ri
i i a¼ P
m
ð8Þ
The equilibrated particles are rotated with respect to each other in order to investigate the effect of crystallographic orientation on the process of neck formation. The particles are then placed next to each other such that the distance between the nearest atoms from two contacting particles is about 1.7r. In order to prevent the loss of any vapour atoms from the system, periodic boundary conditions are employed, so that a vapor atom leaving the simulation cell through the boundary is re-introduced into the cell with the same velocity at the opposite boundary. The size of the simulation cell is large enough to prevent atoms on the particles from interacting with their neighbor images except for the vapor atoms. Following Ch’ng and Pan, [8] three different particle arrangements are used in order to study the particle interactions. In configuration I, a small particle whose radius is 16r is placed next to a bigger particle which is twice the size of the small one. In configuration II, two same size particles are placed in contact. Their radius is 32r. In configuration III one small particle is placed along a straight line between two big particles (see Fig. 2). The larger par-
Fig. 2. The three different configurations of particles: (I) two particles of different size; (II) two particles with same size; (III) one small particle in the middle of two large particles.
249
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
ticle consists of around 3000 atoms while the smaller one consists of around 1000 atoms. The simulation results, including atom coordinates, velocities, and potential energy are recorded every 500 time steps. Thus the simulation progress could be observed directly from dynamic movie obtained by plotting all the atom coordinates. When two particles join together they create a ‘‘neck”. We monitor neck growth with time by measuring its width, defined as the distance between the outer atoms at the narrowest place. We also measure the particle center to centre distance determined by measuring the distance between mass centers of the core parts of the particles with a radius about 10r which remain crystalline throughout the simulation and thus provide a faithful measure of the relative particle displacement during sintering. 2.2. The continuum model The solid state diffusion theory originally due to Herring [19] is adopted for the continuum model. According to this theory, the diffusive flux j, defined as volume of matter passing across unit area perpendicular to the flux direction per unit time, is assumed to be linearly dependent on the gradient of the chemical potential l of the diffusing species
j¼
Dd ol ; kB T os
ð9Þ
where D is diffusivity, d is the thickness of the layer through which the material diffuses, kB is the Boltzmann constant, T is the absolute temperature and s is the local coordinate along the diffusion path. Along a grain-boundary, the gradient of the atomic chemical potential is induced by the gradient of stress r acting normal to the boundary [20]
l ¼ Xr;
ð11Þ
ð12Þ
ð13Þ
where cs is the specific surface energy and j is the principal curvature of the surface. It is positive for a concave free surface. So we have
js ¼
Ds ds Xcs oj kB T os
In a MD simulation, it is inevitable to have certain amount of vapour atoms. When a vapour atom flies across the periodic boundary, the angular momentum of the system changes causing the system to move and rotate rigidly. Sintering occurs by atoms migrating mainly along the particles surface. Atoms at the core of a particle only vibrate thermally and move rigidly. The translational and rotational rigid motions have to be eliminated, as shown in Fig. 3, during post-processing in order to filter out information such as diffusion coefficient, neck size, centre to centre distance, etc. This numerical detail is presented here to demonstrate the feature of the MD model that does not exist in the continuum model. Considering a rigid core of a particle which contains n atoms, after certain time t, the rigid core rotates by an angle a about its mass centre. Inside the rigid core, the atom coordinates before and after rotation are denoted by ri and r0i . The angle a of the rigid rotation of the particle core can be calculated as
where di = rircm is the atom coordinate with respect to the mass centre of the particle core. The rigid motion can then be subtracted from the atomic coordinates according to the following coordinate transformation equations:
Along a free surface, the gradient of the atomic chemical potential is induced by the gradient of the free surface curvature [20]
l ¼ Xcs j;
3.1. Subtraction of rigid motion of the entire system
ð10Þ
where Db is the grain-boundary diffusivity and db is the grainboundary thickness. As matter is deposited onto (or removed from) a grain-boundary, the grains on either side of the boundary move with a normal velocity vb. Matter conservation requires that
ojb þ mb ¼ 0: os
3. Molecular dynamics simulation of nanoparticles sintering
P 0 jdi di j sina ¼ Pi 0 0 ; i jdi jjdi j
where X is the atomic volume, therefore
Db db X or jb ¼ ; kB T os
numerical scheme to solve the coupled diffusion problem which is used in this study. When comparing with the molecular dynamics model it is important to bear in mind that a continuum model has to assume a ‘‘mechanism of matter redistribution”, grain-boundary diffusion for example, while the MD model does not make such an assumption. It is also useful to point out that Herring’s theory assumes thermodynamic equilibrium at every step of the microstructure evolution while the molecular dynamics model has no such limitation.
ð16Þ
x00i ¼ ðy0i ycm Þcosa þ ðy0i ycm Þsina þ xcm ;
ð17Þ
y00i ¼ ðx0i xcm Þsina þ ðy0i ycm Þcosa þ ycm ;
ð18Þ
r 0i
ðx0i ; y0i Þ
r 00i
ðx00i ; y00i Þ
where ¼ and ¼ and after rotation subtraction.
denote atom coordinates before
3.2. Volume fraction of atoms at the particle surface From the point of view of the continuum model, the driving force for particle sintering is the reduction in the free energy asso-
ð14Þ
in which Ds is the free surface diffusivity. As matter is deposited onto (or removed from) a free surface, it migrates with a normal velocity vs. Matter conservation requires that
ojs þ ms ¼ 0 os
ð15Þ
with ms assumed to be positive when matter is added to the surface. The grain-boundary and surface diffusion processes are coupled at the triple contact point where the grain-boundary meets the free surface. Furthermore, the grain-boundary stress has to satisfy an equilibrium condition. Pan and Cocks [15] developed a simple
Fig. 3. Mathematical subtraction of rigid rotation during post processing.
250
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
Fig. 4. Number of atoms as a function of its potential energy indicating their locations at different parts of the particle during sintering.
ciated with the particle surface. In a continuum model, the surface energy is simply proportional to the total surface area. It is therefore instructive to monitor the volume fraction of atoms on the particle surface in the MD simulation to check against the continuum concept. In order to do this, it is necessary to distinguish atoms from different parts of the system, i.e. on the surface, in the bulk, on the grain-boundary or in the vapour. This can be achieved by examining the potential energy of interaction of an atom with other atoms. We record the average potential energy of each atom over 500 time steps. The histogram in Fig. 4a shows the number of atoms as a function of their potential energy. It can be clearly observed that the number of atoms peaks around several values of the potential energy, each representing a group of atoms with a particular number of nearest neighbour atoms. The numbers of neighbouring atoms are indicated by colours as shown in the fig-
Fig. 5. Volume fraction of surface atoms as a function of sintering time.
ure. The peaks from left to right represent atoms with 6, 5, 4, 3, 0 neighbours respectively. The left prime peak represents the majority atoms in the bulk which have all six neighbours arranged in a correct crystal configuration. The sub-peak under the prime peak represents atoms on the grain-boundary (i.e. distorted crystal configuration). All the other peaks except for the vapour atoms represent atoms on the particle surface. This identification of different atoms can be used to monitor mass transport during the sintering process. Fig. 5 shows the volume fraction of the surface atoms (with potential energy larger than 5.0) as a function of time for the two particles sintering in configuration II. It can be seen clearly that the total number of surface atoms decreases as sintering proceeds, which is consistent with the continuum concept that the driving force for sintering is the reduction in the total surface energy.
Fig. 6. Particles with initially aligned crystal orientations maintain the alignment during sintering; only half of the particles are shown to reveal their relative rotation.
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
Fig. 7. Particles with initial misaligned crystal orientations. The small particle reorients itself at the very beginning of the sintering process; only half of the particles are shown to reveal their relative rotation.
3.3. Effect of crystalline orientation of the particles The information that the continuum model does not contain is the crystalline orientation of the particles. It turns out that this information in the MD model leads to a sintering behaviour which the continuum model is unable to capture. Fig. 6 shows two particles with perfectly aligned initial crystalline orientation while Fig. 7 shows three particles of different initial crystalline orientations. Only half of the particles are shown to reveal the relative rotation between the particles during sintering. The MD simulation shows that the particles with perfectly aligned initial crystalline orientation stay aligned during sintering. For particles of different initial crystalline orientations, a quick adjustment occurs in the
251
crystalline orientations at the very beginning of the sintering progress. The particles then remain in the adjusted orientations throughout the rest of the sintering process. This phenomenon has been observed by other researchers both experimentally [21] and in MD simulation [5]. The reorientation is believed to be driven by the minimization of the grain-boundary energy. However the consequence of the reorientation on sintering has not been fully realised in previous studies. In our MD simulation, it is observed that the reorientation can eliminate a grain-boundary from the system as shown in Fig. 8 and significantly retard the sintering process. In some other cases, a stable grain-boundary is formed after the initial adjustment, as shown in Fig. 9. Our simulations show that three factors can help the formation of a stable grain-boundary: (a) a large angle of initial misalignment between the particles, (b) big particle size and (c) higher sintering temperature. In the two dimensional MD model, the maximum misalignment angle is 30 degrees because of the 6-fold symmetry of the crystalline structure. The initial misalignment has to be large enough in order to resist the elimination of the grain-boundary by the particle reorientation. Our simulation suggests that the grain boundaries are more likely to form between two large nanoparticles if the initial misalignment angle is larger then 20 degrees. Furthermore the MD simulation shows that a grain-boundary is more likely to form between the two large nanoparticles in configuration II (see Fig. 2), but rarely form between two smaller nanoparticles. The smaller particles are relatively easy to reorient themselves to eliminate the grain-boundary. A typical evolution of the system in configuration III, is that the small particle quickly aligns with one of the large particles to form a single crystal while a grain-boundary is formed between the small particle and the other large one. 3.4. Sintering of particles The primary conclusion of the continuum modelling by Ch’ng and Pan [8] is that nanoparticles arranged in the configuration III shown in Fig. 2 seem to possess certain resistance to particle coarsening while the configuration I does not have such resistance. However using the continuum model for particles consisting of only a few thousands of atoms pushes its validity to the limit. One major feature
Fig. 8. Two nanoparticles sinter together to form a single crystal, the lines mark the crystalline orientations before and after the reorientation.
252
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
Fig. 9. Two nanoparticles sinter together to form a grain-boundary; the lines mark the crystalline orientations before and after the reorientation.
Fig. 10. MD simulation of one small and two large particles sintering at T = 0.26 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
Fig. 11. MD simulation of a small and a large particle sintering at T = 0.26 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
that the continuum model is unable to capture is the elimination of the grain-boundary as shown in Section 3.3. If a grain-boundary is eliminated, the grain-boundary diffusion process ceases to exist, leaving only surface diffusion and vapour evaporation and condensation as the mechanisms for matter redistribution. This phenomenon has a profound effect on the sintering behaviour of the particles.
253
The melting temperature of the two-dimensional Lennard– Jones crystal is Tm = 0.415 e/kB [22]. We performed MD simulations at two different temperatures: T = 0.19 e/kB and T = 0.26 e/kB, which are 46% and 63% of the melting temperature, respectively. It is clear that sintering occurs much faster at the higher temperature. For the two identical particles of initial diameter
Fig. 12. MD simulation of two large particles sintering at T = 0.26 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
Fig. 13. MD simulation of one small and two large particles sintering at T = 0.19 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
Fig. 14. MD simulation of a small and a large particle sintering at T = 0.19 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
254
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
Fig. 15. MD simulation of two large particles sintering at T = 0.19 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
Fig. 16. MD simulation of two large particles sintering with the grain-boundary at T = 0.26 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
70r (configuration II in Fig. 2) without a grain-boundary, the neck size reached 30 r at T = 0.26 e/kB after 25000 t*, but only reached 20 r at T = 0.19 e/kB by the same time. If a grain-boundary was formed, then the neck size reached 45 r at T = 0.26 e/kB after 25000 t*. Figs.10–12 show the MD simulations for the three particle configurations sintering at T = 0.26 e/kB. We have picked three cases for which the grain-boundaries were eliminated at the beginning of the simulation. In all figures showing MD simulations the following colouring scheme is used: the bulk atoms are always coloured as blue, any atom that has been on the surface or the grain-boundary of the particle is coloured as red,1 and any atom that has been in the vapour is coloured as cyan. The colouring scheme therefore traces the history rather than the current state of the atoms, providing a rough picture of the matter redistribution during the simulation run. In Fig. 10 one small particle is placed between two larger particles. At the very beginning, the particles quickly joined together and formed two necks. Although there was some time lag between the formations of the two necks, by 100 t* both necks were fully filled by surface atoms (shown in red) and the blue bulk part of the small particle has reduced in size. At the later stage, both surface diffusion and vapour condensation contributed to the growth of the neck. Fig. 11 shows one small particle and one large particle in contact with each other. Initially, the quick neck formation was sim-
1 For interpretation of the references to colors, the reader is referred to the web version of this article.
ilar to the case in Fig. 10. The neck was filled by atoms mostly coming from the small particle through surface diffusion. In the later stage of the sintering, evaporation and condensation dominated the matter redistribution. After the neck has been filled up the top of the small particle started to collapse and the whole system evolved toward a circular shape. Fig. 12 shows two large particles of the same size in contact with each other. Again the sintering process started with a quick neck formation followed by neck growth first by surface diffusion and then by vapour condensation. If one uses the continuum model, such behaviour of switching from one mechanism of matter transportation to another half-way through the sintering process would be very difficult, if not impossible, to handle. These three cases therefore highlight the difference between the MD model and the continuum model. The locations of the centres of mass of the particle cores were monitored during the simulation. In all the three cases, there was no perceptible change of the distance between the particles. This indicated the absence of grain-boundary diffusion and that any bulk deformation of the particle was negligible. The simulations were then repeated at T = 0.19 e/kB at which evaporation is much less important. Figs. 13–15 show these cases. It is obvious from these figures that at lower sintering temperatures surface diffusion dominates the matter redistribution. In all the above cases there was no grain-boundary between the particles in contact. Fig. 16 shows the case of two particles of the same size sintering at T = 0.26 e/kB for which a grainboundary is formed between the particles. By t = 100 t*, a small
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
Fig. 17. Centre to centre distance as a function of time for the case shown in Fig. 16.
255
Fig. 17 shows the centre to centre distance between the two particles as a function of time. It is clear that the particles are approaching each other during the sintering process confirming that atoms have been taken away from the grain-boundary by the grain-boundary diffusion. This can be compared with the case of Figs. 12 and 15 for which there was no grain-boundary between the two particles and the distance between particles remained unchanged. Fig. 18 shows the neck size as a function of time for the cases showing in Figs. 12, 15 and 16. It can be observed from the figure that grain-boundary diffusion promotes the necksize growth. Faster surface diffusion plus evaporation condensation at higher temperatures are also among the factors promoting faster necksize growth. Fig. 19 shows the MD simulation of one small particle in contact with two big neighbours sintering at T = 0.26 e/kB. This is perhaps the most interesting case which demonstrates the difference between the MD model and the continuum model. By t = 100 t* a grain-boundary was formed at the upper neck but the initial grain-boundary at the lower neck has been eliminated. Then two very different behaviours of the neck growth can be observed. The upper neck grows by grain-boundary diffusion while the lower one grows by surface diffusion. The grain-boundary also migrates into the small particle. It is also obvious that the existence of a grain-boundary accelerates significantly the neck growth at the upper neck. A continuum model on the other hand cannot possibly predict such a behaviour as the mechanism of matter transportation has to be assumed rather than predicted by such a model. 4. Direct comparison between MD and continuum models
Fig. 18. Comparison of neck size growth as a function of time for the cases shown in Figs. 12, 15 and 16.
neck made up of surface atoms has been formed. As sintering proceeds, the neck is filled up dominantly by the red atoms, meaning that these atoms have all been at the particle surface or the grain boundary at some stage of the sintering. The grain-boundary was migrating backward and forward between the two particles.
The most interesting outcome of the MD simulations presented so far is that asking if a MD model agrees with a continuum model is the wrong question. For nanoparticles, the sintering behaviour observed in a MD simulation is beyond the predictability by the continuum model. If the assumed mechanism of matter redistribution in the continuum model did occur in a MD simulation, then it would be interesting to see if the two models can predict the same speed of neck growth. Fig. 20 compares the neck sizes as functions of time predicted by the continuum and MD models respectively. The case is for two particles of the same size of R = 70 r sintering at T = 0.19 e/kB which is controlled by surface diffusion. Logarithmic scales are used for both axes. The straight line was obtained from the continuum model while the discrete dots were obtained from several independent MD simulations of the system initialised within the same macroscopic states (i.e. same atom positions and temperature) but different microscopic states (i.e. randomly
Fig. 19. MD simulation of one small and two large particles sintering at T = 0.26 e/kB. (a) t = 500 t*; (b) t = 50000 t*; (c) t = 200000 t*; (d) t = 500000 t*. See text for the explanation of the colour scheme.
256
L. Ding et al. / Computational Materials Science 45 (2009) 247–256
uum theory of solid state physics does not apply at this scale, but because nanoparticles behave in so many different scenarios of the continuum theory that one has to predict these scenarios using the molecular dynamics model. In our MD simulation it has been observed that the particles may reorient themselves to match their crystalline orientations at the beginning of the sintering, and thus form different types of necks between different particles (either with or without the grain-boundary). This leads to different mechanisms of matter redistribution at the different necks. It has also been observed that the particles switch the mechanism of matter transportation half-way through the sintering process. None of these can be handled by the continuum model. However assuming the right scenario, the continuum theory does agree with the MD simulation for particles consisting of just a few thousand atoms. Fig. 20. Neck size divided by initial particle radius as a function of time predicted by the MD model (discrete dots) and continuum model (solid line) for the sintering of two particles of the same size, each consisting of 5000 atoms at T = 0.19 e/kB. Surface diffusion is the dominant mechanism for matter redistribution. Logarithmic scales are used for both axes.
Acknowledgements This work was supported by the EPSRC research grant S97996 which is gratefully acknowledged. One of the authors (RLD) did part of the work during his study leave granted by the University of Leicester. The computations were performed on the University of Leicester Mathematical Modelling Centre’s cluster, which was purchased through the EPSRC strategic equipment initiative. References
Fig. 21. Neck size divided by initial particle radius as a function of time predicted by the MD model (discrete dots) and continuum model (solid line) for the sintering of two particles of the same size, each consisting of 5000 atoms, at T = 0.26 e/kB. Grain-boundary diffusion is the dominant mechanism for matter redistribution. Logarithmic scales are used for both axes.
chosen initial velocities of individual atoms). Fig. 21 shows similar comparison where grain-boundary diffusion is the controlling mechanism of matter redistribution. The sintering temperature was set as T = 0.26 e/kB. All the other parameters are the same as that in Fig. 20. The MD prediction is smoother in Fig. 21 than Fig. 20 because the neck size was estimated manually from a single MD simulation. It is remarkable to see that the continuum model agrees with the MD model for particles consisting of only around 10,000 atoms. 5. Conclusions For nanoparticles, it is inappropriate to use continuum models to predict their sintering behaviour. This is not because the contin-
[1] R. Birringer, Materials Science and Engineering, A: Structural Materials: Properties, Microstructure and Processing A117 (1989) 33–43. [2] H. Gleiter, Progress in Materials Science 33 (1989) 223–315. [3] I.W. Chen, X.H. Wang, Nature (London) 404 (2000) 168–171. [4] H. Zhu, R.S. Averback, Materials Science and Engineering, A: Structural Materials: Properties, Microstructure and Processing A204 (1995) 96–100. [5] P. Zeng, S. Zajac, P.C. Clapp, J.A. Rifkin, Materials Science and Engineering, A: Structural Materials: Properties, Microstructure and Processing A252 (1998) 301–306. [6] J.S. Raut, R.B. Bhagat, K.A. Fichthorn, Nanostructured Materials 10 (1998) 837– 851. [7] K. Tsuruta, A. Omeltchenko, R.K. Kalia, P. Vashishta, Europhysics Letters 33 (1996) 441–446. [8] H.N. Ch’ng, A.C.F.C.J. Pan, Acta Materialia 55 (2007) 813–824. [9] J.R. Groza, R.J. Dowding, Nanostructured Materials 7 (1996) 749–768. [10] M.J. Mayo, D.C. Hague, D.J. Chen, Materials Science and Engineering, A: Structural Materials: Properties, Microstructure and Processing A166 (1993) 145–159. [11] R.W. Siegel, S. Ramasamy, H. Hahn, Z. Li, T. Lu, R. Gronsky, Journal of Materials Research 3 (1988) 1367–1372. [12] J.A. Eastman, Y.X. Liao, A. Naravanasamy, R.W. Siegel, Materials Research Society Symposium Proceedings 155 (1989) 255–266. [13] R.W. Siegel, Materials Research Society Symposium Proceedings 196 (1990) 59–70. [14] C. Herring, Journal of Applied Physics 21 (1950) 301–303. [15] A.C.F.C.J. Pan, Acta Metallurgica Material 43 (1995) 1395–1406. [16] J.E.C. Lennard-Jones, Proceedings of the Physical Society 43 (1931) 461. [17] J.A. Snyman, Applied Mathematical Modelling 6 (1982) 449. [18] W.G. Hoover, A.J.C. Ladd, B. Moran, Physical Review Letters 48 (1982) 1818. [19] F. Herring, Surface Tension as a Motivation for Sintering, McGraw- Hill, New York, 1955. [20] C. Herring, The Physics of Powder Metallurgy, McGraw-Hill, New York, 1951. [21] J. Rankin, B.W. Sheldon, Materials Science and Engineering A: Structural Materials: Properties, Microstructure and Processing A204 (1995) 48–53. [22] J.A. Barker, D. Henderson, F.F. Abraham, Physica A: Statistical Mechanics and Its Applications (Amsterdam, Netherlands) 106A (1981) 226–238.