Atomistic simulation of mass transport phenomena in inhomogeneous intermetallic systems

(1)

Atomistic simulation of mass

transport phenomena in

inhomogeneous intermetallic

systems

Piotr Sowa

Thesis presented for the degree of

Doctor of Physical Sciences

Faculty of Physics, Astronomy

and Applied Computer Science

Jagiellonian University

in Kraków, Poland

(2)

(3)

Wydział Fizyki, Astronomii i Informatyki Stosowanej Uniwersytet Jagielloński

Oświadczenie

Ja niżej podpisany Piotr Sowa (nr indeksu: 1015966) doktorant Wydziału Fizyki, Astronomii i Informatyki Stosowanej Uniwersytetu Jagiellońskiego oświadczam, że przedłożona przeze mnie rozprawa doktorska pt. „Atomistic simulations of diffusion in inhomogeneous intermetallic systems” jest oryginalna i przedstawia wyniki badań wykonanych przeze mnie osobiście, pod kierunkiem prof. Rafała Abdank-Kozubskiego. Pracę napisałem samodzielnie.

Oświadczam, że moja rozprawa doktorska została opracowana zgodnie z Ustawą o prawie autorskim i prawach pokrewnych z dnia 4 lutego 1994 r. (Dziennik Ustaw 1994 nr 24 poz. 83 wraz z późniejszymi zmianami).

Jestem świadom, że niezgodność niniejszego oświadczenia z prawdą ujawniona w dowolnym czasie, niezależnie od skutków prawnych wynikających z ww. ustawy, może spowodować unieważnienie stopnia nabytego na podstawie tej rozprawy.

Kraków, dnia . . . .

(4)

(5)

Atomistyczne symulacje transportu masy w

niejednorodnych układach międzymetalicznych

Piotr Sowa

Streszczenie

Opracowanie spójnej metodyki modelowania zjawisk transportu masy w skali atomowej w układach krystalicznych stanowi ogólny cel niniejszej pracy. W szczególności uwaga została poświęcona badaniu dwóch procesów kontrolowanych przez migrację atomów: samo- i interdyfuzji oraz relaksacji porządek-porządek uporządkowania atomowego w związkach międzymetalicznych.

Praca składa się z trzech części. W pierwszej części zostały szczegółowo omówione atomistyczne podstawy procesów dyfuzji i relaksacji. W części drugiej przedstawiono metodologię przeprowadzonych badań, gdzie opisane zostały przede wszystkim modele badanych układów oraz zastosowane techniki symulacji Monte Carlo. W ostatniej części zaprezentowano uzyskane wyniki dotyczące (i) samo-dyfuzji i relaksacji porządek-porządek w układach jednorodnych oraz (ii) procesu interdyfuzyji w układach heterogenicznych (pary dyfuzyjne).

Wyprowadzona zależność między termodynamicznymi energiami aktywacji dla samo-dyfuzji i relaksacji uporządkowania atomowego pozwoliła przedstawić mierzalne doświadczalnie współczynniki samo-dyfuzji i czasu relaksacji za pomocą elementarnych wielkości, takich jak stężenie defektów i częstotliwość skoków atomowych. Poprawność znalezionej formuły sprawdzono poprzez symulowanie obu procesów w dwóch różnych systemach ze strukturą krystaliczną typu B2. Wykorzystano model energetyczny w ramach Hamiltonianu Isinga. Podczas modelowania kinetyki procesów uwzględniono zależność barier migracji dla przeskoków atomów do wakancji od lokalnej konfiguracji atomowej.

Zastosowano dwie techniki symulacji Monte Carlo. Równowagowe konfiguracje, które uwzgledniały temperaturową zależność koncentracji wakancji, wygenerowano za pomocą symulacji Semi-Grand Canonical Monte Carlo (SGCMC). Natomiast migrację atomową za pośrednictwem mechanizmu wakancyjnego symulowano za pomocą algorytmu Kinetic Monte Carlo (KMC). Odtworzono i wyjaśniono doświadczalnie zaobserwowaną w związku międzymetalicznym NiAl zależność między energiami aktywacji dla samodyfuzji i relaksacji porządek-porządek.

Z kolei badanie dyfuzji w stopach heterogenicznych zostało zrealizowane przez symulację transportu masy w parze dyfuzjnej. Główny nacisk położono na opracowanie fizycznie ugruntowanego modelu źródeł wakancji, które mogłyby być zastosowane w metodologii Monte Carlo do symulacji. Nowy algorytm został zaimplementowany w standardowym algorytmie KMC. Bazuje on na równoważeniu stężenia wakancji powstałych w wyniku procesu interdyfuzji, zgodnie z lokalną konfiguracją atomową w próbce. Testy zaproponowanego modelu źródeł wakancji przeprowadzono w oparciu o współczynniki interdyfuzji, efekty korelacyjne oraz zjawisko Kirkendalla wynikające z symulacji układów niewykazujących uporządkowania, jak i układach z uporządkowaniem atomowym. Wyniki uzyskano pośrednio w oparciu o teorię Darkena-Manninga na podstawie współczynników korelacji zmierzonych przez symulacje KMC oraz bezpośrednio za pomocą metody Boltzmanna-Matano zastosowanej do wygenerowanych metodą KMC profili koncentracji składników. Otrzymane rezultaty pozwoliły ocenić prawidłowość opracowanego modelu bezpośrednich symulacji pary dyfuzyjnej poprzez porównanie rezultatów wynikających z obu podejść.

(6)

(7)

Atomistic simulation of mass transport phenomena in

inhomogeneous intermetallic systems

Piotr Sowa

Abstract

The general objective of the present Thesis is the elaboration of a consistent methodology for simulation-based modelling of atomic-migration-controlled phenomena in crystalline systems. Of interest are two atomic-migration controlled processes: self- and interdiffusion and order-order relaxations in chemically ordering systems.

The Thesis consists of three parts. The theoretical background of the atomistics of diffusion and diffusion-controlled processes is discussed in detail in the Part I. Part II contains the description of the applied methodology of the research. The description covers the applied models of the investigated systems and the applied techniques of Monte Carlo simulations. The obtained results concerning (i) self-diffusion and order-order relaxations in homogeneous systems and (ii) interdiffusion process in heterogeneous systems (diffusion couples) are presented in the Part III. The relation between thermodynamic activation energies for self- diffusion and the order-order relaxation was derived. The experimentally measurable self-diffusion coefficient and the relaxation time of order-order relaxations – the two parameters quantifying the rates of the processes – have been expressed in terms of more elementary quantities, such as defect concentration and atomic jump frequencies. The validity of the derived dependency was checked by simulating both processes in two different binary systems showing B2-ordering. The systems were modelled with Ising Hamiltonians completed with migration barriers (saddle-point energies) partially dependent of the local atomic configurations. Two Monte Carlo techniques were used. While the equilibrium atomic configurations of the systems including the temperature-dependent equilibrium vacancy concentration were generated by means of the Semi-Grand Canonical Monte Carlo (SGCMC) simulation, the vacancy-mediated atomic migration was simulated by the Kinetic Monte Carlo (KMC) algorithm. The relationship between the activation energies for self-diffusion and order-order relaxations experimentally observed in the NiAl intermetallic was reproduced and explained. The study of diffusion in heterogeneous alloys was realised by the simulation of mass transport in a diffusion couple. The main emphasis was put on the development of the physically grounded model of vacancy sources and sinks, which could be applied in the Monte Carlo methodology of simulations. The elaborated algorithm was implemented in the standard KMC algorithm. Tests were performed on various alloy systems. The vacancy concentration generated during the process of interdiffusion was on-line locally equilibrated, according to the virtual atomic configuration in the sample. The evaluated interdiffusion coefficients, as well as the correlation and Kirkendall effects resulting from the simulation of the binary disordered and ordered systems, were analysed either indirectly by using KMC-determined correlation factors, or directly by means of the Boltzmann-Matano formalism applied to KMC-generated concentration profiles. The correctness of the proposed model of direct interdiffusion simulations was verified by comparing the results following from both tests.

(8)

(9)

Acknowledgements

I would like to express my sincerest gratitude to my supervisor Prof. Rafał Abdank-Kozubski for the guidance and assistance throughout my thesis work. He was very patient in his mentorship and spent considerable time working with me and giving me precious advises. I also especially thank Prof. Graeme Murch and Prof. Irina Belova for scientific discussions and friendliness. I also wish to thank my family and friends for their support and motivation.

I gratefully acknowledge financial assistance from Polish National Science Center [Grant No. 2015/16/T/ST3/00501], Polish Ministry of Science and High Education [Grant No. 3135/7. PR/2014/2], European Community’s Seventh Framework Programme (FP7 PEOPLE-2013-IRSES) [Grant No. EC-GA 612552] and Endeavour Fellowship and the Australian Research Council Discovery Grants schemes [Grant No. ERF-RDDH-5049-2016].

(10)

(11)

List of Figures

2.1 Scheme of a tracer diffusion experiment. Deposited on the surface isotope atoms diffuse into the sample. . . 26 2.2 Cross section of the experimental set–up used by Kirkendall and co-workers [15].

Black circles denote Mo wires. The dotted line shows initial contact plane. Incorporated Mo markers moved towards each other. . . 29 2.3 Scheme of diffusion between neighbouring crystallographic planes. . . 34 2.4 The elementary unit cell of the B2 superstructure. Empty circles denote the

α-sublattice, and filled circle indicates β-sublattice. . . . 39

3.1 Variation of the system configuration energy due to a jump performed by a p-type atom initially occupying the i-th lattice site and moving to a nearest neighbour vacancy residing on the j-th lattice site. . . . 48 3.2 An example of vacancy concentration dependence on the relative chemical potential

∆µAV for AB-V system. The dotted line indicates the coexisting phases and equilibrium relative chemical potential ∆µ(eq)_AV for A atoms. . . 53 3.3 An example of a developed compositional profile (a), atomic fluxes (b) and dJV/dx

(c) for the case of D∗_A> D∗_B. . . 57 3.4 Generation (a) and annihilation (b) of the vacancy due to the dislocation climbing

down and climbing up respectively. . . 58 3.5 An example of shrinking of the sample due to the annihilation of the vacancies.

Dashed line marks the initial size of the sample. dx is a total change of the size when the dislocation is removed from the sample. . . 58 3.6 Scheme of vacancy annihilation. (a) Change in the distortion in the real material

due to the dislocation climbing up; (b) A real sample is mapped onto the perfect lattice. The path of a hypothetical dislocation sliding is shown; (c) The modelled sample after annihilation of the vacancy. . . 59 3.7 Cross section of the diffusion couple with embedded compositional profile developed

in the sample (bold solid line). Dashed line indicate the regions where local vacancy concentration CV and composition CA is monitored during simulation with the equilibration algorithm (Alg. 6). Initial diffusion couple was wrapped with empty sites on the left-hand and right-hand side of the sample which is marked by the dotted lines. . . 61

5.1 Equilibrium relative chemical potentials ∆µAV, ∆µBV simulated for ’DIS’ (a) and ’ORD’ (b) models. . . 104 5.2 Concentration of C_A(β) antysites (_{) and C}_B(α) antysites (_{#) in the ’DIS’ system (a)}

and in the ’ORD’ system (b). . . 105 5.3 Equilibrium vacancy concentration C_V(eq)for the ’DIS’ system (a) and for the ’ORD’

system (b). . . 105 5.4 Thermodynamic factors simulated for ’DIS’ (a) and ’ORD’ (b) models: circles (_#)

denote φA; squares () denote φB. . . 106 5.5 The ratio wA/wB of the jump frequencies for the ’DIS-ASY’ system at T = 1000K. 107 5.6 (a) Jump frequencies of A atoms wA indicated by filled circles ( ) and jump

frequencies of B atoms wB denoted by empty circles (#) for ’ORD-SYM’ system at T = 1400K; (b) The ratio wA/wB of the jumps frequencies for ’ORD-SYM’ system at T = 1400K. . . . 108

(14)

List of Figures

5.7 (a) A atom jump frequencies, (b) B-atom jump frequency and (c) the ratio wA/wB of the jumps frequencies for ’ORD-ASY’ system at T = 1400K. . . . 108 5.8 Tracer correlation factors fifor the ’SYM-ASY’ system at T = 1000K: (a) A-atom

tracer correlation factor and (b) B-atom tracer correlation factor. . . 109 5.9 Tracer correlation factors fifor the ’ORD-SYM’ system at T = 1400K: filled circles

(_{) denote A-atom tracer correlation factor and empty circles (#) denote B-atom} tracer correlation factor. . . 109 5.10 Tracer correlation factors fifor the ’ORD-ASY’ system at T = 1400K: (a) A-atom

tracer correlation factor and (b) B-atom tracer correlation factor. . . 110 5.11 Self diffusion coefficients for the ’DIS-ASY’ system at T = 1000K: (a) A-atoms

self-diffusion coefficient and (b) B-atoms self-diffusion coefficient. . . 110 5.12 (a) Self diffusion coefficients for the ’ORD-SYM’ system at T = 1400K. Filled

circles (_{) denote A atom jump frequencies and empty circles (#) denote B-atom} jump frequency. (b) The ratio of the self diffusion coefficients for the ’ORD-SYM’ system at T = 1400K. . . . 111 5.13 Self diffusion coefficients for the ’ORD-ASY’ system at T = 1400K. (a) Self diffusion

coefficient for A atoms. (b) Self diffusion coefficient for B-atoms. (c) The ratio of the self-diffusion coefficients. . . 111 5.14 Collective correlation factors for ’DIS-SYM’ (a) and ’DIS-ASY’ (b) systems at

T = 1000K. Filled circles (_{) denote f}AA; Empty circles (#) denote fBB; Filled squares (_{) denote f}AB; Empty squares () denote fBA. . . 112 5.15 Collective correlation factors for ’ORD-SYM’ (a) and ’ORD-ASY’ (b) systems at

T = 1400K. Filled circles (_{) denote f}AA; Empty circles (#) denote fBB; Filled squares (_{) denote f}AB; Empty squares () denote fBA. . . 112 5.16 Interdiffusion coefficient (a) and vacancy wind factor S (b) for the ’DIS-ASY’ system

at T = 1000K. The dotted line denotes the theoretical prediction from the MAA theory (Sec. 2.4.3). . . 113 5.17 Inter diffusion coefficient (a) and vacancy wind factor S (b) for the ’ORD-SYM’

system at T = 1400K. . . . 114 5.18 Inter diffusion coefficient (a) and vacancy wind factor S (b) for the ’ORD-ASY’

system at T = 1400K. . . . 114 5.19 A diffusion couple for the DIS system built out of two samples with 45 × 25 × 25 bcc

unit cells each. Diffusion couple was wrapped by the vacancies (light grey). Sample were generated with the equilibrium number of vacancies. A atoms are marked as a dark grey colour. B atoms are marked by a grey colour. . . 115 5.20 Scheme of the moving average. The y value for xn is calculated on the base of the

Eq. 5.3 by taking points X2, X3, X4, X5. . . 115 5.21 Integration procedure over time and space of concentration. . . 116 5.22 Time evolution of the interdiffusion coefficient ˜D at the stoichiometry composition

CA= 0.5 for the ’DIS-SYM’ system at T = 1000K. . . . 116 5.23 Averaged ˜D(x) for the ’DIS-SYM’ system at T = 1000K. Open circles (_{#) show}

the raw data of BM analysis of profiles from the saturated region. Filled circles (₎ are the averaged values. . . 117 5.24 Interdiffusion coefficient for ’DIS-ASY’ system simulated with standard KMC

methods without equilibration mechanism at T = 1000K. Open circles (_{#) shows} results obtained from BM analysis applied to the profiles from the saturated region (Fig. 5.22). Averaged inter-diffusion coefficient is denoted with filled circles (_{). . . 117} 5.25 Vacancy distribution developed in the diffusion couple of ’DIS-ASY’ system at

T = 1000K without equilibration mechanism. . . . 118 5.26 Interdiffusion coefficient for ’DIS-ASY’ system at T = 1000K without vacancy

equilibration mechanism. . . 118 5.27 Interdiffusion coefficient for disordered ’DIS-SYM’ system at T = 1000K. . . . 119 5.28 Intediffusion coefficient for the case of the active source/sink on the A-rich side (a)

and B-rich side (b) for ’DIS-SYM’ system at T = 1000K. Empty circles denote (_#) results obtained form B-M analysis for compositional profile from saturated region. Filled circles (_{) denote averaged interdiffusion coefficient. . . 120} 5.29 Results obtained for the ’DIS-ASY’ system at T = 1000K. . . . 120

(15)

List of Figures

5.30 Interdiffusion coefficient for ’DIS-ASY’ system simulated with equilibration mechanism (Alg. 6) at T = 1000K. . . . 121 5.31 Vacancy distributions developed in the diffusion couple of ’DIS-ASY’ system at

T = 1000K simulated with the equilibration mechanism (see Alg. 6). . . . 122 5.32 Interdiffusion coefficient for ’ORD-ASY’ system simulated with equilibration

mechanism (Alg. 6) at T = 1400K. . . . 122 5.33 Snapshots of the evolving (moving) sample in the [0, 1, 0] direction. The initial

position of the sample (reference frame) is marked by horizontal solid lines. The x-axis shows the distance from the initial contact plane with constant lattice a0 units. Grey colour denotes A atoms; Black colour denotes B atoms; White colour denotes V atoms (vacancies). . . 123 5.34 Time evolution of the Kirkendall plane position xK measured for ’DIS-ASY’ system

at T = 1000K. The initial sample was built from pure A sample on the left-hand side and pure B sample on the right-hand side (see Fig. 5.33a). . . 124

(16)

(17)

List of Tables

5.1 The interaction energy parameters Vpq (see Sec. 3.1.1). . . 103 5.2 Migration barrier parameters E_bar+ (Sec. 3.1.2) for A and B atoms. . . 104 5.3 Monitored parameters. . . 107

(18)

(19)

List of Algorithms

1 Scheme of the direct exchange algorithm. . . 50

2 Scheme of the Semi-Grand Canonical Monte Carlo algorithm. . . 51

3 Scheme of the Residence Time algorithm. . . 51

4 The algorithm of path localisation. . . 59

5 The algorithm for defining the path direction K. . . . 60

(20)

(21)

List of Symbols

a0 : Lattice constant.

C : Volume based concentration.

c : Mole fraction (site concentration), Eq. (2.6). N : Number of particles.

M : Total number of sites.

J : Flux of atoms across the plane, Eq. (2.1). D : Diffusion coefficient, Eq. (2.1).

D∗ : Tracer (Self) diffusion coefficient, Sec. 2.2.1.

D(I) _{: Intrinsic diffusion coefficient, Eq. (2.3).}

µ : Chemical potential.

φ : Thermodynamic factor, Eq. (2.9). ˜

D : Interdiffusion coefficient, Sec. 2.2.3. η : Long range order parameter, Eq. (2.99). τ : Relaxation time, Eq. (2.25).

Π : Frequency of an atom jumps (probability per unit time), Eq. (2.87).

Π0 : Timescale factor, Eq. (3.3).

Γ : Frequency of an atom jumps to neighbouring vacancy in the crystalline structure (probability per unit time), Eq. (2.94).

ν : Jump frequency of the atoms. Total number of atoms jumps per unit time, Eq. (2.82). w : Average jump frequency of an atom equal to the total number of atoms jumps per unit

time per number of atoms, Eq. (2.81).

λ : Length of the jump.

r : Displacement of an atom in a single jump.

R : Total displacement of an atom after n jumps, Eq. (2.62). ∆R : Total displacement of all atoms in the system, Eq. (2.61).

fi : Tracer correlation factor, Eq. (2.65).

(22)

(23)

Chapter 1

Introduction

1.1 Mass transport phenomena

The phenomenon of mass transport in materials is of crucial importance in the theoretical and applied research because many different processes and properties of the material are controlled and affected by the behaviour of single atoms. The other term naming the transport of matter is diffusion. Diffusion is an essential mechanism of any structural transformation including chemical ordering, phase transitions or crystal growth. It is relevant for the kinetics of many microstructural changes that occur during the processing of materials. Typical examples are nucleation, phase transformations, homogenisation and recrystallisation of alloys. The knowledge on the atomistic mechanisms and time scales of diffusion and diffusion controlled processes governing both the fabrication and the stability of the material is essential for the material manufacturers.

Diffusion processes can be categorised into two main groups due to the way how diffusion can be observed. The flow of matter is described as a chemical diffusion in the case of inhomogeneous systems, for example, two materials of the different composition being in contact (welded together). Such a system, where atoms diffuse across the contact plane is often called a diffusion couple. In contrary, the flow of matter in a chemically homogeneous environment is referred to as a self-diffusion process. There are various experimental methods for studying diffusion in solids. The standard technique used in studies of self-diffusion is radiotracer method, often combined with techniques of depth profiling such as Secondary Ion Mass Spectrometry (SIMS) or Electron MicroProbe Analysis (EMPA) [1]. Other types of methods that can be used to study diffusion over a broad spectrum of time scales are Nuclear Magnetic Relaxation (NMR), Resonance Mössbauer spectroscopy and X-ray Photon Correlation Spectroscopy (XPCS) [2].

Each technique has its own advantages. Among all of the methods, the NMR covers the widest range of diffusivities. With nuclear resonant reflectivity, it is possible to study the motion of selected atoms but it requires samples with a one dimensional superstructure, preferably of a period just about one order of magnitude above the lattice constant. Very low diffusion rates can be measured by XPCS, which are in the order of 10−23m2/s−1. It based on a monitoring of changes in a distribution of speckles obtained by scattering of coherent X-rays on the material. A distinct correlation between the arrangement of atoms and the speckle pattern gives information about the atomic dynamics. This knowledge enables one to establish a model of the atomistic mechanisms of diffusion in a material.

Similarly, as in other domains of natural sciences, the primary method of research is an experiment. However, despite being the most direct approach, it is sometimes difficult to carry out or too expensive and often requires additional knowledge to interpret the results. Hence, various theories of diffusion were developed which allow modelling of fluxes of matter or evolution of a diffusion couple on the macroscopic level. The most common theoretical models for the case of a disordered alloy are for example those elaborated by Darken [3], Manning [4] or Moleko [5].

Atomistic aspects of diffusion and diffusion-controlled processes can be successfully modelled using simulation techniques [6]. By means of appropriate Monte Carlo (MC) or Molecular Dynamics (MD) algorithms, recently accompanied by the atomistic variants of the Phase Field method (Phase Field Crystal, Atomic Density Function), diffusion may be directly modelled as a collective process composed of elementary atomic jumps, whose particular mechanism is often implemented (especially in the case of standard MC algorithms). In the case of intermetallics and

(24)

1.2. Motivation for the study

high ordered alloys, this is predominantly the vacancy mechanism.

Application of the simulation techniques requires that a proper model of the examined system is available – i.e. its Hamiltonian must be determined on the ground of a particular theory. Nowadays still more and more common models based on Density Functional Theory (DFT) are used in ab-initio MD or MC simulations of systems. Such techniques are, however, still applicable only to small systems (clusters built of maximum 103 atoms). Modelling of nanostructured systems still most often involves effective (quasi-empirical) potentials determined either by fitting to experimental results, or recently by fitting to either experimental, or recently to DFT calculation results.

Atomistic simulation is very often used as a validation tool for the macroscopic theories of diffusion. Moreover, nowadays, they are an integrated part of advanced experimental techniques. The significant advantage of the atomistic simulations of diffusion is a direct availability of a full analysis of almost all features of the process including correlation effects, jump frequencies, energy barriers or diffusion mechanism.

The general objective of the present Thesis is the elaboration of a consistent methodology for simulation-based modelling of atomic migration-controlled phenomena in crystalline systems. The method is based on:

• An atomistic model of a system parameterized by atomic interactions and the geometry of its crystalline structure.

• Monte Carlo algorithms for the simulation of both the equilibrium observables (Semi Grand Canonical MC) and the atomic migration controlled kinetics (Kinetic Monte Carlo, implemented with a particular mechanism of atomic movement).

Of interest are two atomic-migration controlled processes: (i) self- and interdiffusion and (ii) order-order relaxation in chemically order-ordering systems.

The first part of the Thesis is focused on physical aspects of self-diffusion and order-order relaxations being examples of the steady-state and non-steady-state processes. The experimentally observed differences between the features of the two phenomena are elucidated regarding the properties of the collectivity of the migration of individual atoms. The novelty of this approach consists of the comparative atomistic analysis of the two phenomena showing the effect of the steady- and non-steady-state on the elementary processes controlling both of them.

The second part of the Thesis is devoted to the design of an atomistic model of the simulation of diffusion couple experiments. A diffusion couple is a chemically inhomogeneous system. The study of diffusion directly in atomic scale in such conditions is crucial in understanding the nature of the interdiffusion mechanism on the most basic level. In particular, it allows to explain experimental results, as well as to test new theories concerning multicomponent systems and systems where analytical methods cannot be used due to complex boundary conditions.

In the following sections, details of the problem raised in this Thesis will be discussed along with the motivation and scheme of the study.

1.2 Motivation for the study

Our goal is to simulate the dynamic evolution of atoms. The principal tool in this class of atomistic simulation methods is molecular dynamics, which is grounded on a solving of classical equations of motion. However, its serious limitation is that short time steps (10−15s) are used to do an accurate calculation. As a result, the time achieved in a simulation is limited to microseconds, while diffusion processes usually occur on much longer scales.

Monte Carlo techniques overcome this limitation by treating dynamics of alloys as a series of diffusive jumps from state to state instead of following the trajectory of all atoms in every vibrational step, as it is explained in the Sec. 3.2. Thanks to this, Monte Carlo simulations can reach much longer times. Within this approach it possible to simulate entire processes and to interpret the results in terms of statistical thermodynamics. The evaluation of thermodynamic activation energy of various diffusion controlled processes is one of the examples.

Relation between the thermodynamic activation energy for diffusion and structure relaxation

Self-diffusion and order-order relaxations are fundamental processes in many technologies aiming at the generation of the desired microstructure of advanced materials based on high ordered alloys [7].

(25)

1.2. Motivation for the study

Since both phenomena need thermally activated atomic migration, it is, in general, expected that when occurring in the same system they would give rates that are similarly temperature dependent. However, the recently reported resistometric measurements of Partyka, Nikiel [8], [9] revealed that it is not always the case and that in B2-NiAl intermetallic the efficiency of thermal activation of ordering/disordering is lower than the efficiency of thermal activation of diffusion kinetics. The result is non-trivial due to the same atomic jumps being involved in both processes of diffusion and ordering/disordering.

It can be understood in the way that in case of order–order transformation system is under driving force coming from non-equilibrium (caused by the change of the temperature) condition which can activate additional mechanism for atomic jumps, which is not the case of diffusion process occurring in condition where degree of order is preserved. The most likely key to understanding the observed effect is the fact that the vacancy-mediated atomic migration may involve elementary atomic jumps showing specific correlations depending on whether the process runs in stationary (tracer diffusion) or non-stationary (ordering/disordering) conditions.

The purpose of the present studies is to provide a systematic and consistent analysis of the dynamics of vacancy-mediated atomic jumps in the above two sorts of conditions. The relationships between the dynamics of diffusion and chemical ordering in systems modelled with different Hamiltonians will be analysed by performing appropriate atomistic simulations. Such investigations addressing systems that show, or not, a tendency for the triple-defect formation will be the subject of the first part of the Thesis.

Vacancy thermodynamics and direct simulations of a diffusion couple experiment As already mentioned in the introduction, diffusion couple experiment is a common way to study mass transport phenomena in the solid state. Depending on the special experimental conditions, different macroscopic diffusion parameters can be measured by analyzing the measured concentration profiles. The development of the contemporary advanced technologies created, however, a definitive demand for detailed atomistic description of the diffusion mechanisms which, despite the availability of particular experimental techniques (e.g. XPCS), is the most straightforwardly provided by atomistic simulations.

The underlying assumption of the atomistic simulations is that the heterogeneous sample can be approximated as a series of homogeneous subsystems. Accordingly, the process of atomic migration resulting from the concentration gradient of the components is expressed in terms of the parameters of diffusion occurring in an environment free from any gradients. Such an approach limitates of course, the scope of diffusional processes which can be studied – e.g. it is commonly assumed that the vacancies stay in thermodynamic equilibrium. The latter is possible due to the activity of various vacancy sources and sinks, everywhere in the crystal during diffusion.

The reasoning is, however, not generally applicable as one can imagine the situations with an insufficient density of vacancy sources and sinks causing too slow vacancy relaxation. This effect is increased especially in the cases where the strong compositional dependence of vacancy concentration is observed like it is in the case of the intermetallics [10], [11]. Dynamic evolution of composition and accompanied vacancy concentration during interdiffusion can quite significantly reduce the reliability of the used assumptions. Experimentally this effect may be observed as a difference in the results for chemical and tracer diffusion because the two types of diffusion experiment might refer to the different distribution of vacancy concentration. In this context it appears very useful to consider vacancies as an additional component of the system [12] whose equilibrium concentration is determined in terms of specific phase equilibria and include effects of vacancies directly in the simulations.

The aim of this part of the Thesis is to develop a physically grounded algorithm for creating and annihilating vacancies in a diffusion couple during diffusion running in the concentration gradient. The algorithm will enable a direct simulation of a diffusion couple experiment, which in turn makes it possible to tests phenomenological theories and study the effect in more complex systems where non-equilibrium effects are important. The new model of vacancy generation/annihilation is developed starting with the assumption of fast vacancy relaxation. By using the atomistic simulation techniques, the correctness of the proposed model is tested on disordered and ordered alloys systems. The work addresses the fundamental problem of chemical interdiffusion in binary alloys in the case where vacancies with non-equilibrium concentrations are generated during the interdiffusion process. The elaborated model can by applied to study effects of non-equilibrium distribution and localised sources/sinks of vacancies.

(26)

1.3. General scheme of the study

1.3 General scheme of the study

The Thesis consists of three parts. The theoretical background of the atomistic of diffusion and diffusion-controlled processes is discussed in detail in the Part I. Part II contains the description of the applied methodology of the performed research. The description covers the used models of the investigated systems and the applied techniques of Monte Carlo simulations. The obtained results corresponding to the two goals of the study are presented in the Part III. The description of the obtained results is followed by a chapter containing their discussion and conclusions. The full bibliography of the works studied when preparing the Thesis is presented at the end.

(27)

Part I

(28)

(29)

Chapter 2

Elements of the theory of diffusion

The subject of diffusion consists of a wide range of topics covering different aspects of mass transport properties. Only a basic description of the bulk diffusion in binary systems is given in this chapter. Other aspects of diffusion in more complex systems are not necessary for the understanding of problems raised in this Thesis.

In this chapter, the ground terminology is introduced in chronological order. The Section 2.2 covers a general overview of diffusion in the solid state. Starting from a macroscopic level, we introduce the reader to the basic definitions of experimentally measurable observables which include various types of diffusion coefficients. Next, in Section 2.3 a generalised thermodynamic approach is presented. In the second part of this chapter (Sec. 2.4) the atomistic underlay of defined macroscopic parameters is explained on the basis of the random walk theory. Finally, atomic-scale details of diffusional processes with a non-equilibrium environment such as structural relaxation are discussed in the Sec. 2.4.5.

2.1 Diffusion mechanisms

The diffusion is considered as a phenomenon of matter transport. This means at the atomistic scale a movement of atoms/particles whose nature definitely differs between fluids and solids. While in fluids all degrees of freedom are available for the particle movement, in the solid state the movements of atoms are dominated by atomic vibrations. Diffusive jumps are usually single-atom jumps showing fixed distances, which in crystalline systems are of the order of the lattice parameter.

In solids, translational atomic displacements are strongly inhibited and occur only in specific conditions. Thus, only a specific diffusion mechanisms are allowed. The most common are direct exchange mechanism, interstitial mechanism, vacancy mechanism or ring mechanism. The simplest one is the exchange mechanism, which consists of a simultaneous direct interchange of the positions of two neighbouring atoms. This type of mechanism dominates in very loosely packed crystals. The interstitial mechanism operates in systems built of atoms with strongly different sizes – i.e., containing small impurity atoms. The smaller atoms can then move by jumping between the interstitial positions. The vacancy mechanism dominates in most of metallic systems such as substitutional solid solutions and intermetallic phases. As in thermal equilibrium metallic crystals always contain a certain number of vacant sites-vacancies, which provide an easy path for diffusion, the process of diffusion can proceeds by the jumps of atoms to these vacancies.

2.2 Phenomenological description of diffusion

The Fick’s laws (Eq. (2.1)) constitute the essential foundation of diffusion processes. These equations for the first time appeared in the original work of Adolf Fick in 1855 [13] where they were used to represent a continuum description of diffusion of a salt-water system. The concept of the diffusion coefficient which was introduced by Fick suggested a linear response between the concentration gradient and the rate of mixing. Fick’s laws describe the diffusion as an empirical fact, without relying on any mechanism.

(30)

2.2. Phenomenological description of diffusion

It is known from observation that when heterogeneous monophasic systems are annealed the material flows in such a way that the gradients are smoothed. After a long enough time, the sample becomes homogeneous and the net flow of matter is no longer observed. In the case of this type of problem, Fick’s first law relates the components fluxes, defined as the number of particles passing through an area unit per unit of time with its concentration gradient on that plane. In the case of a one-dimensional system, if the x axis is taken parallel to the concentration gradient of the component i the flux Ji of the component along the gradient is given by

Ji= −Di ∂Ci

∂x . (2.1)

Here Ci denote the flux and concentration of the of the i-type respectively. The proportionality coefficient, Di, is named a diffusion coefficient or diffusivity of the species i, it has the dimension of the units [m2_s−1_].

The Fick’s law is valid for different diffusional processes, but the diffusion coefficient may have different meanings depending on the particular circumstances.

2.2.1 Tracer diffusion

A particularly simple situation where Fick’s law has an application is an experiment with a diffusion of tracer isotope in an otherwise homogeneous crystal.

A concentration gradient of the radiotracer is created by deposition of a thin layer of the tracer element onto the surface of a sample. The only factor that can lead here to the net flux of tracer is the concentration gradient of the tracer itself provided the total tracer concentration is small enough to keep the composition of the sample unchanged and to make sure that the tracer atoms do not interact one with another. Such a conditions lead to the diffusion independent on the local concentration of isotope. After deposition, an isothermal annealing is performed for some time, and the tracer concentration-depth profile is analysed. A typical configuration of the experiment for a tracer/self-diffusion is illustrated in the Fig. 2.1. The most important task is to derive the diffusion coefficient by comparing the depth profile with the appropriate solution of Fick’s law.

Described experiment set-up allows studying a diffusion in a chemically homogeneous environment. In this case, the diffusion of an atom is called tracer diffusion and is expressed by the tracer diffusion coefficient D_i∗, where i stands for the different atomic species. The term self-diffusion coefficient refers to the particular case where the tracer atoms are the same species as the non-tracer atoms in the crystal.

Figure 2.1: Scheme of a tracer diffusion experiment. Deposited on the surface isotope atoms diffuse into the sample.

2.2.2 Intrinsic diffusion coefficient

In the preceding section, the diffusion was limited to the situation where there were no additional effects which could influence the movement of tracer atoms. Nevertheless, in the case when for example external physical forces are present, it can lead to the situation where the diffusing particles experience extra drift motion. The same effect can arise not only from the applied external forces but also from other internal factors such as interactions between diffusing atoms. All those additional factors are called driving forces (see Sec. 2.3). Those effects are included in the diffusion equation for atoms (Eq. (2.1)) now written as

Ji= −D∗i ∂Ci

∂x + CihviF, (2.2)

The term CihviF has a dimension of the flux and originates from the fact that the individual atoms feel some non-zero velocity – the drive velocity hvi_F – in x direction due to the presence

(31)

of driving forces. In cases when the drift velocity is proportional to the concentration gradient ∂Ci/∂x Eq. (2.2) can be rewritten as

Ji= −D (I) i

∂Ci

∂x . (2.3)

Here D(I)_i is the intrinsic diffusion coefficient of the species i.

Driving forces which are proportional to the concentration gradient are those following from the diffusion induced gradient of an electric field in ionic crystal or from the excess of chemical potential due to the presence of non-zero interaction terms between components.

Regular solution

As already mentioned in the previous section, in the case of a reaction occurring between the constituents and leading to the change of the local-surrounding-dependent properties of the atoms, it is a gradient of the chemical potential ∂µi/∂x of the component i which acts as a driving force. The average drift velocity (Eq. (2.2)), in this case, can be expressed in terms of the Nernst-Einstein relation [1], which enters to the equation for the flux of component i in the form:

hvi_F = D ∗ iFi kBT (2.4) and Ji= −Di∗ ∂Ci ∂x + CiD∗iFi kBT , (2.5)

where Fi = −∂µ0i/∂x. Here µ0i has a meaning of the contribution of all factors other than the entropy of mixing.

The general relation between the chemical potential and concentration is given by

µi= µ0i+ kBT lnci, (2.6)

where ci = Ci/PiCi is a mole fraction. It can be further expressed in terms of the activities ai and the activity coefficient γi of the component

µi= µ (0)

i (p, T ) − kBT lnai= µ (0)

i (p, T ) − kBT lnγici, (2.7)

where µ(0)_i (p, T ) is a standard chemical potential for given pressure p and temperature T . Equations (2.6) and (2.7) yield

Fi= ∂µ0_i

∂x = kBT ∂lnγi

∂x . (2.8)

Substituting Eq. (2.8) into Eq. (2.5), replacing the volume base concentration by the mole fraction and remembering that the total concentration of all components is constant, one obtains the relation between the intrinsic and the tracer diffusion coefficients

D_i(I)= D_i∗(1 +∂lnγi ∂lnci

) = D_i∗φi. (2.9)

The quantity φi is called a thermodynamic factor. In ideal solid solutions φi= 0. Additional effect appears in the case of φi6= 0 (non-ideal solid solutions). Then, according to the Eq. (2.9) the rate of mixing of the components can be increased, what is the case of a solid solution with a positive heat of mixing or decreased in case of a negative heat of mixing.

2.2.3 Interdiffusion

So far, described was diffusion in homogeneous systems and defined were self- and intrinsic diffusion coefficients to describe different effects connected with mass transport of atoms. Let us now focus on binary AB alloys.

In general, the diffusivity of atoms may differ for different components. Hence, if two semi-infinite bars with the various concentrations of A and B atoms are joined (forming a diffusion couple) and annealed, atoms start to penetrate the bar volumes with different ratios and according to Eq. (2.9) two fluxes JAand JBare generated in the sample. As a single diffusional process – the

(32)

intermixing (homogenization) of the binary AB system, is effectively observed, the total rate of this homogenisation can be described by a single chemical (interdiffusion) coefficient ˜D(C). Because of the concentration gradient, the composition permanently changes in particular positions of the couple and the interdiffusion coefficient also changes along the sample.

The important theoretical problem of the relation of this single interdiffusion coefficient ˜D(C) with the intrinsic and self-diffusion coefficients of the A and B components is discussed in detail in the next sections.

2.2.4 Nernst-Planck equation

An example of the operation of a driving force in an ideal solid solution case is an intermixing of two species of cations A and B in an ionic crystal where the anions remain fixed in their positions. If the diffusion of the cations is vacancy mediated then at the absence of driving forces, the fluxes JA and JB are related to the concentration gradients via Eq. (2.1). With only two species of cations we have that ∂CA/∂x = −∂CB/∂x.

When D_A∗ and D_B∗ are not equal, a net flux of electric charge JA6= JBwould develop. However, the condition of electrical neutrality requires that this net flux must be zero. In practice, an electric field ED arises as an effect of diffusion of cations and decrease the further diffusion of cation. The driving force F resulting from this field equals qED, where q is the cation charge.

According to the Nernst-Einstein relation (Eq. (2.4)) and Eq. (2.2) the drift term is given by

Ji= −Di∗ ∂Ci ∂x + qCiD∗i kBT ED. (2.10)

Setting JA= JB and solving for ED it follows that

ED= kBT q D∗ A− D∗B CAD∗A+ CBDB∗ ∂CA ∂x . (2.11)

Inserting Eq. (2.11) to Eq. (2.10) and comparing with Eq. (2.3) so called Nernst-Planck equation for ideal solution is obtained

˜ D = D(I)_A = D(I)_B = D ∗ AD ∗ B NADA∗ + NBDB∗ . (2.12)

This equation was elaborated on the example of electric driving force, but it applies to any situation where diffusion is constrained, i.e. where the bottleneck for diffusion is present. It can be seen in the case when D∗_A>> D_B∗. From Eq. (2.12) we have that the overall rate of mixing is still limited by the B component diffusivity.

As reported by Svoboda [14], an extra term CV/C eq

V also appears when off-equilibrium vacancy concentration is present. Finally

˜ DN ernst−P lanck= CV C_Veq D∗_AD∗_B NADA∗ + NBDB∗ . (2.13)

For the case of a non-ideal solid solution D(I)_A and D(I)_B should be additionally multiplied by the appropriate thermodynamic factor φA(B) (see Eq. (2.9)).

2.2.5 The Kirkendall effect

Another relation between interdiffusion coefficient and self-diffusion coefficients is obtained in the case of a binary diffusion couple consisting of the species showing different values of intrinsic diffusion coefficients D(I)_A 6= D(I)_B .

In this case, a net atom flux occurs across any plane in the diffusion zone. It causes the crystal to swells on one side of the diffusion plane and shrinks on the other side. Kirkendall and coworkers discovered this effect (Kirkendall shift) in a copper-brass diffusion couple in the 1947 [15].

The original Simglskas and Kirkendall experiment was performed using the Cu/Zn diffusion couple and showed a shrinking of the inner brass core (see Fig. 2.2) due the difference between the rates of the Zn- and Cu-atoms diffusions (D_Zn(I) > D_Cu(I)). The Kirkendall shift was measured by incorporating inert inclusions, called markers (e.g., Mo or W wires, ThO2 particles), at the interface where the diffusion couple was initially joined.

(33)

2.2. Phenomenological description of diffusion Cu + 30% Zn 100% Cu 100% Cu Cu Brass Cu

Figure 2.2: Cross section of the experimental set–up used by Kirkendall and co-workers [15]. Black circles denote Mo wires. The dotted line shows initial contact plane. Incorporated Mo markers moved towards each other.

This experiment was a break-through in the studies of diffusion. It proved that the diffusion in crystals occurs mostly by the vacancy mechanism. Before it was commonly believed that diffusion in solids takes place via direct exchange or ring mechanism (see Sec. 2.1), which implied that the diffusivities of both components of a binary alloy must be equal. The fact that in a diffusion process running in a solid-state the species may diffuse with different rates gave rise to fundamental reformulations of the existing atomistic models on the solid-state diffusion.

From an atomistic point of view the Kirkendall shift can be explained as a result of a generation and annihilation of vacancies. When the flux of one component is higher than the fluxes corresponding to the remaining elements a vacancy flux in a direction opposite the flux of more mobile particles is generated to compensate the difference between the A and B fluxes. Consequently, regions of vacancy under- and super-saturation are created. As the system itself tends to maintain equilibrium vacancy concentration the vacancies must be either created, or annihilated on the respective sides of the contact surface. Each generation (annihilation) of a vacancy generates localised displacement of sample which in overall cause lattice drift. In the case of insufficient plastic relaxation during the process, the vacancies can cluster and form pores or voids called Kirkendall porosity and Kirkendall voids. The initial contact plane (marker plane) is referred to as the Kirkendall plane.

2.2.6 Darken equation

In 1948 Darken proposed a phenomenological analysis of the Kirkendall’s and co-worker’s results and gave the first mathematical description of interdiffusion [16].

The correct description of the diffusion requires that a good reference system is chosen. Darken assumed that the marker wire identifies a given lattice plane. Thus, the Fick’s first law for the diffusion fluxes with respect to the local crystal lattice are

JA= D (I) A ∂CA ∂x , JB = D (I) B ∂CB ∂x . (2.14)

If the plane is moving with the velocity vK (Kirkendall velocity) with the respect to the ends of the diffusion couple (laboratory reference frame), the flux J_i0 of species i with respect to the ends of the sample is Ji0 = −D (I) i ∂Ci ∂x + CivK. (2.15)

vK is a velocity appearing in addition to the driving velocity vF. The total flux of atoms Jtot0 with respect to the laboratory frame is

J_tot0 = JA0 + JB0 = vK(CA+ CB) − (D (I) A − D (I) B ) ∂Ci ∂x . (2.16)

If the ends of the diffusion couple are far away from the interface so that there are no concentration gradients and no diffusion occurs at these points, then Jtot0 must be equal to zero, as the ends of the sample stay fixed. Then the velocity vK can be expressed as

vK = (D (I) A − D (I) B ) ∂cA ∂x , (2.17)

(34)

where cA= CA/(CA+CB) and ∂cA/∂x denotes the molar concentration gradient at the Kirkendall plane.

Substituting Eq. (2.17) to Eq. (2.15) one arrives at a general expression for the interdiffusion coefficient ˜ DDarken= cAD (I) B + cBD (I) A , (2.18)

which gives the rate at which the concentration gradient tends to smooth itself out.

2.2.7 Vacancy flow therm – Manning factor

The Darken model presented above relates the coefficients for inter-, intrinsic- and self-diffusion of the diffusion couple components and provides a valid approximation of the physical reality. However, it may be regarded as only a first approximation which assumes that there is no contribution of correlation effects to the experiment in which ˜D is determined (see Sec. 2.3 and Eq. (2.40)).

In the case of the vacancy mechanism Manning showed that the atomic fluxes may be indirectly affected by the net flux of vacancies. It is because an extra flow of vacancies perturbs the equilibrium distribution of vacancies in the neighbourhood of the atom. Effectively, more vacancies approach a given atom from one side. It causes an additional drift of particles in the direction opposite to the vacancy flux. As a result, the flux of either component has an additional part defined by Manning as a vacancy-wind factor.

Assuming that vacancies do not group up and are not bounded by the material, for the case of binary alloy Manning elaborated the following modification for intrinsic diffusion coefficients

D(I)_i = D∗iφiri, (2.19)

where riis a vacancy-wind corrections. An approximate expressions for ri in the framework of an alloy showing no order, are

rA= 1 + 1 − f0 f0 cA(D∗A− D ∗ B) cADA∗ − cBD∗B , rB= 1 − 1 − f0 f0 cB(D∗A− D ∗ B) cADA∗ − cBD∗B . (2.20)

f0 is the geometric tracer correlation factor, a characteristic constant depending on the structure and diffusion mechanism [17].

From Eqs. (2.19) and (2.20) we have that when D∗_A > D∗_B, then the vacancy flow increases (rA > 1) the flux of A atoms and reduces the net flux of B atoms down its gradient (rB < 1). It is understandable as the effective vacancy flux, created due to the difference of diffusivity (Eqs. (2.15), (2.17)) has a direction opposite to the direction of the flux of A atoms and same direction as the flux of B atoms.

Rewriting Darken relation (Eq. (2.18)) in the form that includes the effect of vacancy flow, Manning provided an expression for the total vacancy-wind factor S:

S = 1 +1 − f0 f0

cAcB(D∗A− DB∗)2

(cAD∗B+ cBD∗A)(cADA∗ + cBD∗B)

. (2.21)

The final Darken-Manning relation for inter-diffusion coefficient which gives the rate for smoothing the concentration gradient is:

˜

DDarken−M anning= (cADB∗ + cBD∗A)φS, (2.22)

where the constrain (Eq. (2.24)) yielded by Gibbs-Duhem relations were applied. For the case of binary alloy we have that

cAdµA+ cBdµB+ cVdµV = 0. (2.23)

With equilibrium vacancy concentration (dµV ≈ 0) and dilute number of vacancies (dcA≈ −dcB), the Eq. (2.23) and Eq. (2.9) give

φA= φB= φ. (2.24)

A transparent derivation of Eq. (2.20) can be found in [17], [18]. Vacancy-wind corrections for chemical diffusion in intermetallic compounds depend on the structure, the type of disorder and the diffusion mechanism. The idea of the vacancy-wind has been extended upon multicomponent and ordered systems by Alnatt, Murch and Belova [19], [20].

(35)

2.2.8 Diffusion-controlled processes – order-order relaxation

The preceding sections were devoted to diverse phenomena controlled by diffusion acting as a transport of matter. Introduced and defined were related diffusion coefficients being parameters reflecting the rate of this processes. Diffusion (migration) of atoms controls, however, various structural transformations among which the phenomenon of chemical ordering is crucial in materials science.

Chemical ordering may be investigated directly using various diffraction techniques and indirectly by measuring different chemical-order-dependent material properties (e.g. electrical resistivity). In any case, of interest are so-called ‘order-order’ or ‘order-disorder’ relaxations during which the system evolves from the non-equilibrium towards the equilibrium state (degree) of the chemical order. Such processes naturally run due to atomic migration (displacements). Thus, are diffusion controlled and show particular rates whose quantification needs a parameter (different from the diffusion coefficient). Definition of the new parameter – the ‘relaxation time’ – stems from the experimental observations.

The typical experimental procedure has the following steps:

1. A sample is equilibrated at a given temperature Ti < TC – such equilibration leads to equilibrium degree of the chemical order corresponding to T = Ti with a particular arrangement of atoms (and resistivity of a sample). TC means the critical temperature at which order-disorder transition occurs.

2. The temperature T is increased or decreased to arbitrary chosen temperature Tj < TC. Time-dependence of the measured chemical-order-dependent property X is monitored. The measurement is continued until the time-dependence X(t) saturates at a newly attained equilibrium degree of the chemical order.

Detailed analysis of X(t) curves gives an opportunity to evaluate the time scale coefficient τ (relaxation time) by fitting the phenomenological equation to the experimental results [21]

X(t) − X(t → ∞) X(t = t0) − X(t → ∞) =X i A(τi)exp −t τi ; X i A(τi) = 1, (2.25)

where A(τi) represents a distribution of the particular exponential processes, with the relaxation time τi.

2.2.9 Temperature dependence of diffusion coefficient and relaxation

time

In general, the rate of diffusion and so on following diffusional processes depend on thermodynamic variables such as temperature, pressure, and composition [21]. In this section, we concentrate on the dependence of the diffusivity and relaxation time on temperature.

The self-diffusion coefficient D and the reciprocal relaxation time τ−1depend rather strongly on temperature, increasing with growing temperatures correspondingly. They are frequently, but by no means always, found to obey the Arrhenius formula, at least within finite temperature ranges.

Experimentally measured diffusion coefficient often fit a relation [22], [23]

D = D0exp(−E_a(D)/kBT ), (2.26)

where both Ea(D) and the pre-exponential factor D0 are independent of temperature. Equation (2.26) is called Arrhenius equation for diffusion. By definition, the experimental quantity Ea(D) is given by

E_a(D)= −kB(∂lnD/∂(1/T )) (2.27)

and is called the thermodynamic activation energy for diffusion.

Similarly Arrhenius relation is obeyed by the relaxation times τi [24]–[26]

τ−1 = τ₀−1exp(−E_a(o−o)/kBT ), (2.28)

where Ea(o−o)is the thermodynamic activation energy for disordering(ordering) in case when τ was measured as a result of the Eq. (2.25) applied for increased(decreasing) temperature.

(36)

2.3. Non-equilibrium thermodynamic approach

The exponential dependence of D and τ on temperature, which follows naturally from the kinetic theory, was explained in details in the Sec. 2.4.5. The general increase of D and τ with temperature is a consequence of the particular characters of self-diffusion and order-order relaxation. The physical interpretation of the activation parameters E(D)a ,E

(o−o)

a and of D0, τ0 depends on the diffusion mechanism, on the type of diffusion process, diffusing element and on the lattice geometry. It should be noted that deviations from linearity are often observed. It is an effect connected with the presence of multiple diffusion mechanisms, impurities or grain boundaries.

2.3 Non-equilibrium thermodynamic approach

In the preceding chapters, basic parameters describing diffusion processes have been considered on the base of Fick’s first law. In this section an alternative approach will be given as together with previous it will give a better understanding of different aspects of diffusion phenomena.

The continuum theory of Fick can successfully explain many diffusion phenomena. However, the conditions ∂C/∂x = 0 is sometimes insufficient for describing the diffusion (see Sec. 2.2.2). Equation (2.1) does not recognise others driving forces which can influence diffusion. For example, it does not consider that a concentration gradient of one species gives rise to a flux of another.

The more general approach is based on irreversible non-equilibrium thermodynamic. In the theory of transport developed by Onsager, the flux of matter is related to the driving force Xk by equations [27], [28]

Ji = X

k

LikXk, (2.29)

where the Likare the phenomenological coefficients.

The definition of the driving force includes not only actual forces generated, for example by an electric field, but also forces having any other sources arising from the tendency of the system to return to the equilibrium condition after having fluctuated. It includes, for example, gradients of chemical potential, temperature gradient or stress gradient.

The phenomenological coefficients Likdepends on the thermodynamic parameters (concentra-tion, temperature, etc.) but are not a function of the Xk. The general meaning of the summation, in Eq. (2.29), over all forces, is that each transport process affects all the others processes. While the diagonal terms Lii are connected with a direct impact on the process, the off-diagonal ones determine the indirect influence of each force. At the absence of a magnetic field the following relation, known as the Onsager reciprocity theorem [29] holds for all off-diagonal coefficients

Lik= Lki. (2.30)

Directly access to the phenomenological coefficients through the experiment is rather complicated. Fortunately, however, they can be related to the tracer diffusion coefficient D∗ defined in Sec. 2.2.1.

Consider a one-dimensional diffusion via vacancies in a binary alloy and suppose that the system is isothermal and that external forces are absent. In such conditions, the system consists of two atomic components A and B and vacancies (index V) on a single lattice. In general, we do not erase vacancy influence on the flux of atoms. Then the Onsager equation for this problem are

JA= LAAXA+ LABXB+ LAVXV, JB = LBAXA+ LBBXB+ LBVXV, JV = LV AXA+ LV BXB+ LV VXV.

(2.31)

Assuming vacancy mechanism and the mass conservation it must hold that

JA+ JB+ JV = 0. (2.32)

If this is true for arbitrary values of the forces Xi, substitution of Eqs. (2.31) into Eq. (2.32) gives

LAA+ LBA+ LV A= 0, LBA+ LBB+ LBV = 0, LAV + LBV + LV V = 0.

(37)

2.3. Non-equilibrium thermodynamic approach

Above equations show that the kinetic coefficients of the vacancy flux are related to those of the atomic species (see Sec. 2.2.7). By using the reciprocity theorem (Eq. (2.30)) and by combining Eqs. (2.33) with the flux equations (Eq. (2.31)), the following expressions for the fluxes of the atomic species are obtained

JA= LAA(XA− XV) + LAB(XB− XV), JB= LAB(XA− XV) + LBB(XB− XV), JV = LAV(XA− XV) + LBV(XB− XV).

(2.34)

Equations (2.34) apply to the situation where a vacancy mechanism is operating. If a vacancy mechanism is not present or the vacancies are everywhere in thermal equilibrium, there will be no driving forces coming from vacancies, and XV as well as LiV will equal zero

JA= LAAXA+ LABXB, JB = LABXA+ LBBXB.

(2.35)

The definition of an equilibrium state, saying that the chemical potential of the constituents µi must be the same all over the system, gives that in our case Xi is proportional to the gradient of the chemical potential. System slightly out of equilibrium will tend to return with the rate proportional to the departure from equilibrium, thus

Xi= −  ∂µi ∂x T . (2.36)

By using the definition of the chemical potential (Eq. (2.7)) of a real solid solution, after simple algebra, we arrive at the formula

Xi= −  ∂µi ∂x T = −kBT ∂lnγici ∂x = −kBT φic −1 i ∂ci ∂x. (2.37)

Here φi is a thermodynamic factor introduced in Eq. (2.9).

To relate phenomenological coefficients Lik with the diffusion coefficient D (I)

A within Darken solution (Eq. (2.14)) one has to assume additionally that all off-diagonal coefficients LABand LBA are equal to zero. Then Equations (2.35) are

JA= LAAXA, JB = LBBXB.

(2.38)

Then we need to combine Eq. (2.37) and Eq. (2.38) and compare two expressions for JA, e.g. Eqs. (2.38) and (2.14) remembering that mole fraction is related with volume based concentration via c = CV /M , where M is the total number of sites

JA= −LAAkBT φic−1i V M ∂Ci ∂x = −D (I) A ∂Ci ∂x = −D ∗ Aφi ∂Ci ∂x . (2.39)

Finally, ones arrives at the equation for phenomenological coefficient [16]

Lii=

D_i∗ciM V kBT

; Lij= 0. (2.40)

The Manning relation for Lij coefficients based on the Eq. (2.19) is more accurate as it takes into account all diagonal and off-diagonal parts in Eq. (2.35). The detailed elaboration can be followed in [17]. By separating appropriate contributions to fluxes Ji (Eq. (2.14)) arising from XA and XB, Manning expression for the phenomenological coefficients are found in the form

Lik= D∗_iciM V kBT δik+ 1 − f0 f0 ckD∗k cAD∗A+ cBD∗B . (2.41)

The general relation between the diffusivities and the phenomenological coefficients follows directly from Eqs. (2.35), (2.3) and the expressions for Xi (Eq. (2.37)). Remembering that c = CV /M once obtain D(I)_A = kBT φA V M  LAA cA −LAB cB , D_B(I)= kBT φB V M  LAA cB −LBA cA , (2.42)

Atomistic simulation of mass transport phenomena in inhomogeneous intermetallic systems

Atomistic simulation of mass

transport phenomena in

inhomogeneous intermetallic

systems

Piotr Sowa

Thesis presented for the degree of

Doctor of Physical Sciences

Faculty of Physics, Astronomy

and Applied Computer Science

Jagiellonian University

in Kraków, Poland

Oświadczenie

Atomistyczne symulacje transportu masy w

niejednorodnych układach międzymetalicznych

Piotr Sowa

Streszczenie

Atomistic simulation of mass transport phenomena in

inhomogeneous intermetallic systems

Piotr Sowa

Abstract

Acknowledgements

Contents

I

Theoretical background

23

II

Methodology

45

III

Results and Conclusion

63

List of Figures

List of Tables

List of Algorithms

List of Symbols

Chapter 1

Introduction

1.1

Mass transport phenomena

1.2

Motivation for the study

1.3

General scheme of the study

Part I

Chapter 2

Elements of the theory of diffusion

2.1

Diffusion mechanisms

2.2

Phenomenological description of diffusion

2.2.1

Tracer diffusion

2.2.2

Intrinsic diffusion coefficient

2.2.3

Interdiffusion

2.2.4

Nernst-Planck equation

2.2.5

The Kirkendall effect

2.2.6

Darken equation

2.2.7

Vacancy flow therm – Manning factor

2.2.8

Diffusion-controlled processes – order-order relaxation

2.2.9

Temperature dependence of diffusion coefficient and relaxation

time

2.3

Non-equilibrium thermodynamic approach