Modeling of deformation-enhanced diffusion mass transfer during mechanical alloying in binary metallic substitutional systems
Boris B. Khina
Physico-Technical Institute of the National Academy of Sciences of Belarus, Minsk, Belarus,
Mechanical alloying (MA) is known as a versatile method for producing in situ dispersion-strengthened composites and far-from equilibrium materials such as nanocrystalline alloys and compounds, supersaturated solid solutions and quasicrystalline and amorphous phases, which possess unique properties. Enhanced solid-state diffusion, which was observed experimentally in many systems under the action of intensive periodic plastic deformation (IPPD) during MA, has received substantial interest and was a subject of debates in literature. However, the physics of this complex phenomenon is not entirely understood yet, which impedes further development of MA and novel advanced materials by MA-based technologies. In this work, a mathematical model of enhanced diffusion in a binary substitutional system subjected to IPPD is developed. The model includes a number of physical factors such as generation of non-equilibrium point defects by gliding screw dislocations during deformation, interaction of point defects with edge components of dislocation loops during intervals between collisions, cross-term interaction of diffusion fluxes of lattice atoms and vacancies, interaction of the vacancy flux with a non-coherent phase boundary, etc. Numerical simulations performed for systems Cu-Al and Ni-Cr using realistic parameter values under the deformation conditions typical of MA in a vibratory mill have revealed the mechanism of this phenomenon.
The phenomenon of deformation-enhanced solid-state diffusion in metals and alloys is known experimentally for years . It was observed during mechanical alloying (MA), which is used as a versatile method for producing in situ dispersion-strengthened composites and far-from-equilibrium materials such as nanocrystalline alloys and compounds, supersaturated solid solutions, quasicrystalline and amorphous phases [2, 3]. Besides, enhanced diffusion takes place in repetitive cold rolling, equal-channel angular pressing (ECAP) and other metal processing methods that involve intensive periodic plastic deformation (IPPD). However, the physical mechanism of this intricate phenomenon is not well understood yet, which hinders further development and commercialization of the so-called mechanically driven alloys (MDA) produced by a variety of methods based on IPPD and, in particular, elaboration of novel advanced materials by MA-based technologies.
The role of diffusion processes in the formation of metastable phases during MA has been a subject of hot debates in literature [4–8]. Some authors completely deny the role of diffusion in MA and ascribe the observed effects to the “mechanical intermixing” of atoms under shear straining, but the physical meaning of this notion is not explained in detail nor numerical estimates for the fast intermixing rate are presented [5, 8]. A theory for decomposition of a solid solution (the demixing phenomenon) under IPPD has been developed , where only the rotational modes of deformation of nanograins in the course of MA are considered. In this case, the disclinations at triple grain boundary junctions are considered as the sources of non-equilibrium vacancies. The latter penetrate into the grains of a supersaturated solid solution (SSS) and thus cause enhanced diffusion of solute atoms towards the grain boundaries, which results in demixing. The considered situation corresponds to late stages of MA when the SSS has already formed, i.e., the mechanism of SSS formation is not considered while in most cases the goal of MA is to produce a metastable phase but not to decompose it. Besides, as outlined in Ref. , within this theory the mechanism of vacancy generation by disclinations is not described on the physical level nor a formula (or even a simple estimate) for the vacancy production rate is given. It should also be noted that, according to recent experimental observations, nanoparticles (e.g., Ni nanograins obtained by cold rolling) can contain dislocations whose density is high, = 4.71011-1.31012 cm2 . Therefore, the deformation mechanism via dislocation glide, which is preferable from the energy point of view, can take place in nanomaterials during IPPD.
In modeling attempts [10,11], only separate aspects were studied, e.g., diffusion along curved dislocation pipes  or a change in geometry of an elementary diffusion couple due to deformation  using the trivial Fickian diffusion equation without any additional terms accounting for the generation/annihilation of non-equilibrium point defects during IPPD. Besides, the whole processing time of MA was taken as the time of deformation-enhanced diffusion (from 1 h  to 50 h ) while it is known that during MA the deformation time during a ball-powder-ball or ball-powder-wall collision is much shorter than a time interval between collisions [12,13]. Also, the values of diffusion coefficients D used for calculations in Ref.  were fitted to match the experimental data on the alloying degree instead of estimating them independently.
As was outlined earlier , this situation in the MA science and technology is due to a lack of theoretical understanding of diffusion mechanisms under IPPD on the background of a large amount of experimental data accumulated in literature. In connection with the above, the goal of this work is to develop a rigorous mathematical model of enhanced diffusion in a binary substitutional system subjected to IPPD and perform computer modeling using the experimentally known (e.g., diffusion coefficients that are measured under equilibrium conditions) or independently estimated parameter values (e.g., deformation conditions in a vibratory mill during MA and the generation rate of non-equilibrium point defects). The model is based on certain ideas of the theory of diffusion in irradiated alloys [15, 16].
2. Formulation of the model
During early stages of MA, fracturing and cold welding of initial metal particles result in the formation of composite (lamellar) particles with characteristic lamella thickness L=0.5-0.1 m [2,3]. The physical situation considered in this work is an elementary diffusion couple “phase 1 (pure metal A)/phase 2 (metal B)” in a substitutional alloy system where, in equilibrium conditions, diffusion occurs via the vacancy mechanism, hence JA+JB+JV=0. As ball-powder-ball or ball-powder-wall collisions in the course of MA are chaotic [2, 3, 12], changes in the shape and size of this couple inside a composite particle is neglected. According to the aforesaid, the dislocation mechanism of plastic deformation dominates. Periodic deformation brings about the generation of point defects (non-equilibrium vacancies and interstitial atoms of two sorts: Ai and Bi) in the crystal lattice due to jog dragging by gliding screw dislocations according to the Hirsch-Mott theory . This process in non-conservative: formation of a vacancy actually denotes an increase of the lattice site number while the formation of a self- interstitial atom decreases the number of lattice sites by one. Besides, non-conservative is also the process of adsorption of non-equilibrium point defects by edge dislocation, which is accompanied by the climb of the latter . Hence we have to depart from the assumption about the constancy of the lattice site density, which is traditional in the classical solid-state diffusion theory.
Basing of our previous works [19–21], the diffusion equations for lattice species (atoms A, B and vacancies V) and interstitials of both sorts, Ai and Bi, are written as
where Ck is the concentration of k-th species, Jk is the diffusion flux, k are sink/source terms for k-th species (generation/annihilation rate), N is the lattice site density, = a03 1/N is the average volume per one lattice site, a0 is the lattice period, Dkn the are components of the matrix of diffusion coefficients. The latter can be determined using the theory of diffusion in alloys containing non-equilibrium vacancies  taking into account the fact that when the vacancy concentration CV substantially exceeds the equilibrium value CV0, the diffusion coefficients increase by the factor of CV/CV0:
DBB DB = DB*(gBB–gAB)CV/CV0, DBV = –DB*CB/CV0, DVB = DA–DB,
DAA DA = DA*(gAA–gBA)CV/CV0, DVV = DV+DA, DV = [CBDB*+CADA*]/CV0. (3)
Here DA* and DB* are the self-diffusion coefficients of atoms A and B, which are measured in quasi-equilibrium conditions (i.e. at CV=CV0), g is the thermodynamic factor: gkk = 1 + (ln k)/(ln Ck), gkn=(Ck/Cn)(ln k)/(ln Cn), k,n=A,B, kn, is the activity coefficient. For the sake of simplicity, hereinafter the ideal solid solution is assumed, i.e. gAA=gBB=1, gAB=gBA=0. The values of DA* and DB* obey the Arrhenius law D*=D0exp[–E/(RT)], where E is the activation energy and D0 is the preexponent.
The sink/source terms in the right-hand side (RHS) of Eq. (1) are due to generation of point defects by gliding screws during deformation, annihilation of excess vacancies and interstitial atoms in the bulk, and interaction of non-equilibrium vacancies and interstitials Ai and Bi with edge components of dislocation loops, which act as volume-distributed sinks. The former takes place during the deformation while the latter two processes occur mainly during a pause between sequential deformation events (collisions) in the course of MA. Then, basing on Eqs. (1)-(3), the following reaction-diffusion equations for lattice species (atoms B and vacancies) and interstitials of two sorts are formulated:
CB + CA +CV = 1, (8)
where Pi and Pv are the production rate of interstitials and vacancies, correspondingly, e is the density of edge dislocations, k are the coefficients describing the efficiency of edge dislocations as sinks for point defects of the k-th sort, k ~ 1 [15,16], iv is the annihilation rate of Frenkel pairs Ai-V and Bi-V. In Eqs. (6),(7) the cross-link terms for diffusion of interstitial atoms Ai and Bi are neglected.
The term N/t is determined by the rates of non-conservative processes, viz. generation of non-equilibrium point defects by gliding screws and interaction of the former with edge components of dislocation loops. Thus it can be defined as
The annihilation rate of Frenkel pairs is determined as [15,16]
iv = (4r0/)(Di + DV). (10)
Here Dior and r0 is the capture radius, r0 = b/2 [15,16] where b is the Burgers vector (the shortest spacing between atoms in the crystal lattice). The value of DV is defined above by Eq. (3) while the diffusion coefficients of self-interstitial atoms in substitutional alloys are not known and it is typically assumed that , ~ (101-103)DV [15,16] because the activation enthalpy for migration of self-interstitials 0.1 eV  is substantially lower than that for vacancies (e.g., =0.62 eV for vacancy migration in Al and 1.35 eV for Ni ).
The equilibrium vacancy concentration is calculated as
CV0 = exp[–HV/(kBT)], (11)
where HV is the vacancy formation enthalpy and kB is the Boltzmann constant.
The rate of point defect generation according to the Hirsch-Mott theory  is
Pk = (b/2)(fks)1/2, ki,v, (12)
where is the strain rate (during MA, =101-104 s1 [12,13]), 0.5 is the fraction of forest dislocations, s e = /2 is the density of screw dislocations, is the total dislocation density, fv and fi are the fractions of vacancy and interstitial producing jogs on a gliding screw, fv+fi = 1 and typically fv>fi .
The boundary conditions to reaction-diffusion equations (4)-(7) look as
Jk(x=0) = 0, Jk(x=L) = 0, kB,V,Bi,Ai. (13)
Besides, typically the boundary between phases 1 and 2 with coordinate L/2 is a localized sink for non-equilibrium vacancies, hence an additional boundary condition is to be posed for vacancy diffusion:
CV(x=L/2) = CV0. (14)
The initial conditions are:
, CV(x)=CV0 at 0<x<L, CB(x)=1(x–L/2), (15)
where 1(x–L/2) is the Heaviside step function.
3. results of comluter sinulations and DISCUSSION
The formulated non-linear problem (3)-(15) was solved numerically using the finite-difference method. The parameter values for computer simulations for systems Cu (metal A, phase 1)-Al (metal B, phase 2) and Ni (metal A, phase 1)-Cr (metal B, phase 2) are listed in Tables 1 and 2, correspondingly.
Table 1. Values of diffusion parameters for the Al-Cu system used in modeling .
phase 2 (Al-based)
phase 1 (Cu-based)
Hv,  eV
Table 2. Values of diffusion parameters for the Ni-Cr system at TMA=373 K [25, 26].
phase 2 (Cr)
phase 1 (Ni)
HV , eV
A cyclic process is considered: generation of point defects during collision and diffusion and point defect relaxation during the rest period between collisions. For the Cu-Al system, the parameters corresponded to MA in vibratory mill “SPEX 8000” with the oscillation frequency =20 Hz: the collision duration t~104 s , the cycle duration tc = 1/(2) = 0.025 s, and the rest period tr = tct = 0.0249 s; the deformation rate ~10 s1 was taken as a lower-level estimate. Since the local short-time heating due to a head-on collision is small [12,13], the particle temperature during MA was assumed to be TMA=373 K (100 C), and the dislocation density =1011 cm–2 was taken as a plausible estimate. For the Ni-Cr system, the MA parameters corresponded to the vibratory mill designed in the Belorussian-Russian University (Mogilev): =25 Hz, t=2.510–5 s , then tr=1,997510–2 s, and the strain rate estimated using the concept of Hertzian collision was =3.6104 s1 .
The results of simulations are presented in Figs.1-3. Modeling of enhanced diffusion in the Cu-Al system during MA was performed for two possible situations: (i) the phase boundary 2/1 cannot act as a vacancy sink (as was suggested in Ref.  on a merely qualitative level, at a certain stage of IPPD of nanomaterials, the boundaries of nanograins may be saturated with vacancies and cannot act as vacancy sinks), i.e. boundary condition (14) is not included in simulation (Fig. 1), and (ii) the interface 2/1 acts as a localized sink for non-equilibrium vacancies, i.e. Eq. (14) is taken into account (Fig. 2). It is seen that within the developed model substantial alloying can occur after a relatively short time of MA, 4000 s (line 5 in Fig. 1(a)). A steady-state concentration profile of vacancies establishes in a short time (250 s for Fig. 1 and about 1200 s for Fig. 2). The enhancement of diffusion is due to accumulation of non-equilibrium vacancies and increase in the diffusion coefficients (see Eqs. (3)), slow relaxation of vacancy concentration during pauses between collisions, and the linkage of diffusion fluxes of lattice atoms and vacancies via the off-diagonal diffusion coefficients, as seen from Eqs. (3)-(5). Simulation has revealed that interaction of vacancies with the phase boundary may have a selective effect on the diffusion of atoms A and B in the conditions of MA: strong flux of vacancies towards the interface enhances the flux of foreign atoms in the opposite direction in phase 2 (compare Fig. 1 and Fig. 2). The latter is sometimes called “the inverse Kirkendall effect”. Thus, interaction of diffusion fluxes with the phase boundary appears to be an important factor, which has not earlier received substantial attention in literature. The synergetic influence of the above factors on diffusion mass transfer brings about the formation of a wide zone of solid solution in phase 2 within a reasonably short time of MA, t=4000 s (line 5 in Fig. 2(a)). The occurrence of a small peak on the concentration profile of atoms A near the 2/1 phase boundary (Fig. 2(a), line 5) is due to the fact that the vacancy concentration in this point corresponds to the equilibrium value, hence the values of diffusion coefficients DA and DB are close to those for self diffusion D*, i.e. the phase boundary, which is considered as a sink for non-equilibrium vacancies, can act as a barrier for deformation-enhanced diffusion. The concentrations of interstitial atoms Ai and Bi is much less than that of vacancies, which is due to a lower rate of their generation by jogs on gliding screw dislocations and fast annihilation of the Frenkel pairs.
Fig. 1. Concentration profiles of atoms B (a), vacancies (b) and interstitials Ai (c) and Bi (d) obtained by modeling without boundary condition (14) for different durations of MA: 1, t=0; 2, t=250 s (10000 deformation-rest cycles); 3, t=1250 s (50000 cycles); 4, t=2000 s (80000 cycles); 5, t=4000 s (160000 cycles).
Fig. 2. Concentration profiles of atoms A (a), vacancies (b) and interstitials (c and d) for modeling with boundary condition (14). MA durations are the same as in Fig. 1.
Fig. 3. Profiles of atoms A (Ni) (a), vacancies (b) and interstitials (c) for modeling with boundary condition (14): 1, t=15 s (750 cycles); 2, t=100 s (5000 cycles); 3, t=400 s (20000 cycles); 4, t=1100 s (55000 cycles); 5, t=1800 s (90000 cycles).
Computer modeling of diffusion in the Ni-Cr system under the conditions of IPPD was performed using boundary condition (14), i.e. considering the phase boundary 2/1 as a sink for vacancies. In this situation, atoms A diffuse into phase 2, and a wide zone of supersaturated solid solution forms in a relatively short time of MA, t=1100-1800 s (55000-90000 deformation-rest cycles) (lines 4 and 5 in Fig. 3(a)). The reasons for strong enhancement of diffusion are as described above. Besides, noteworthy is a steep concentration gradient of non-equilibrium vacancies in phase 2 in the vicinity of the interface (Fig. 3(b)). This causes a strong flux of foreign atoms A from the phase boundary into phase 2 (the inverse Kirkendall effect). Far from the interface, the concentration of non-equilibrium vacancies in phase 2 is high, but their gradient is small. Hence diffusion of atoms A inside this phase occurs due to their own concentration gradient but with a high partial diffusion coefficient because CV>>CV0 (see Eqs. (3)). The steady-state profile of vacancies is established after a relatively short time, t=400 s, or 20000 cycles (lines 3-5 in Fig. 3(b)). The interstitial atoms Ai are accumulated mainly in phase 1 because in this phase the annihilation rate of Frenkel pairs, iv, appears to be higher than in phase 2.
Computer simulation within the frame of the developed model using the physically grounded parameter values have demonstrated that enhanced solid-state diffusion in binary metallic systems during IPPD in the conditions of MA is caused by a synergetic action of several physical factors: (i) generation of non-equilibrium vacancies due to jog dragging by gliding screw components of dislocation loops during deformation, (ii) slow relaxation of point defects during rest periods between collisions in the milling device, (iii) increase in the diffusion coefficients due to accumulation of non-equilibrium vacancies, (iv) interaction of diffusion fluxes with the phase boundary, which may act as a sink for vacancies thus changing the direction of their diffusion, and (v) cross-term interaction of diffusion fluxes of different species (lattice atoms and vacancies) via the off-diagonal diffusion coefficients.
The outcomes of simulation have revealed the physical mechanism of enhanced diffusion in metals under the action of IPPD and provide an explanation to numerous experimental data on the formation of supersaturated solid solutions during MA in binary substitutional systems [2-4]. The obtained results agree qualitatively with theoretical work  where the evolution of vacancy concentration during MA was studied using a different approach. Besides, the calculated concentrations of non-equilibrium vacancies generated by IPPD in the conditions of MA agree quantitatively with those measured experimentally in copper after equal channel angular pressing (ECAP), repetitive cold rolling and other deformation methods .
It should be noted that the accumulation of non-equilibrium vacancies during MA causes a distortion of the crystal lattice and increases the free energy of the alloy, and above a certain threshold concentration it can bring about a metastable phase transition, e.g., solid-state amorphization .
1. Larikov L N, Falchenko V M, Mazanko V F, Gurevich S M, Kharchenko G K, Ignatenko A I: Doklady Akademii Nauk SSSR 1975 221 (5) 1073-1075 (in Russian).
2. Suryanarayana C. Progress in Materials Science 2001 46 (1-2) 1-184.
3. Zhang D L: 'Processing of advanced materials using high-energy mechanical milling'. Progress in Materials Science 2004 49 (3-4) 537-560.
4. Skakov Yu A: Metal Science and Heat Treatment 2004 46 (3-4) 137-145.
5. Shtremel' M A:. Metal Science and Heat Treatment 2004 46 (3-4) 146-147.
6. Gapontsev V L, Koloskov V M: 'Metal Science and Heat Treatment' 2007 49 (11-12) 503-513.
7. Skakov Yu A: Metal Science and Heat Treatment 2007 49 (11-12) 514-516.
8. Shtremel' M A: Metal Science and Heat Treatment 2007 49 (11-12) 517-518.
9. Wu X-L, Ma E: 'Dislocations in nanocrystalline grains'. Applied Physics Letters 2006 88 (23) 231911 (3 pp.).
10. Rabkin E, Estrin Y: Scripta Materialia 1998 39 (12) 1731-1736.
11. Mahapatra T K, Das D, Manna I, Pabi S.K: Acta Materialia 1998 46 (10) 3501-3510.
12. Maurice D R, Courtney T H: Metallurgical Transactions A 1990 21 (2) 289-303.
13. Lovshenko G F, Khina B B: Friction and Wear 2005 26 (4) 434-445 (in Russian).
14. Khina B B, Froes F H: 'Modeling mechanical alloying: advances and challenges'. Journal of Metals (JOM) 1996 48 (7) 36-38.
15. Mansur L K: Journal of Nuclear Materials 1979 83 109-127.
16. Johnson R A, Lam N Q: Physical Review B 1976 13 (10) 4364-4375.
17. Nabarro F R N, Basinski Z S, Holt D B: Advances in Physics 1964 50 193-323.
18. Novikov I I: Defects of crystalline structure in metals. Moscow, Nauka, 1983 (in Russian).
19. Khina B B, Solpan I, Lovshenko G F: Journal of Materials Science 2004 39 (16-17) 5135-5138.
20. Khina B B, Lovshenko F G, Konstantinov V M, Formanek B: Metal Physics and Advanced Technologies 2005 27 (5) 609-623 (in Russian).
21. Khina B B, Formanek B: Defect and Diffusion Forum 2006 249 105-110.
22. Gurov K P, Kartashkin B A, Ugaste U E: Mutual diffusion in multiphase metallic systems. Moscow, Nauka, 1981 (in Russian).
23. Brandes E A, Brook G B (editors): Smithells metals reference book, 7th edition. Oxford, Butterworth-Heinemann, 1992.
24. Bokshtein B S: Diffusion in metals. Moscow, Metallurgiya, 1978 (in Russian).
25. Drits M E (editor): Properties of elements: a handbook. Moscow, Metallurgiya, 1985 (in Russian).
26. Grigor’ev I S, Meylikhov E Z (editors): Physical parameters: a handbook. Moscow, Energoatomizdat, 1991 (in Russian).
27. Nechaev Yu S: Physics 2001 65 (10) 1507-1514 (in Russian).
28. Zhang B Q, Lu L, Lai M O:. Physica B: Condensed Matter 2003 325 120-129.
29. Ungar T, Schafler E, Hanak P, Bernstorff S, Zehetbauer M: Materials Science and Engineering A 2007 462 (1-2) 398-401.
30. Fecht H J: Nature 1992 356 (6365) 133-135.