Gottstein G., Shvindlerman L.S. Grain Boundary Migration in Metals: Thermodynamics, Kinetics, Applications

Подождите немного. Документ загружается.

456 5 Computer Simulation of Grain Boundary Motion

the equations of motion. The Verlet time integration scheme [491] deﬁnes for

an atom i

(t + δt)=r

(t)+

(t)

· δt +

(t)

δt

+ ... (5.15)

(t − δt)=r

(t) −

(t)

· δt +

(t)

δt

± ... (5.16)

A second equation for the position of atom i, r

(t + δt), is found when adding

Eqs. (5.15) and (5.16) and taking into account terms up to third order. Such

equation does not explicitly depend on the velocity of atom i,

(t)

,attime

t and reads

⇒ r

(t + δt)=2 · r

(t) − r

(t − δt)+

(t)

· δt

. (5.17)

Eq. (5.17) renders the temporal evolution of the real space positions of each

atom. The velocity of atom, i,attimet,

(t)

, is found when substracting

Eq. (5.15) from (5.16)

(t)

(t + δt) − r

(t − δt)

2 · δt

(5.18)

For an adequate choice of the time step δt one uses the energy scale, length

scale, or the atomic mass. For instance, the x component of force, F

,has

the SI unit [F

]=N =

. The units of J and m are not the natural units

when dealing with atomic-level systems. These are rather the units eV and a

Hence the unit of the x component of force is [F





. The contribution

of a force, F

(t), exerted on atom i at position x at time t + δt within a

Taylor series expansion, allows one to formulate an equation of all entering

units which eventually deﬁnes the unit of the simulation time step, δt:

[x]=a

]

· [dt]

· u

· [dt]

⇔ [dt]=t

unit

= a



=3.68 · 10

−14

s (5.19)

where u is the atomic mass unit which is equal to 1.66055 ·10

−27

kg, and a

is the zero Kelvin lattice parameter of Cu. It is noted that δt depends on the

studied material since it is proportional to its zero Kelvin lattice parameter.

A diﬀerent class of integration schemes only uses the information of the atom position,

velocity and higher time derivatives at time t to predict the t + δt situation. The predicted

atom position, velocity and higher time derivatives at t + δt are then reﬁned by a corrector

procedure. This class of integration schemes is called predictor-corrector scheme. A standard

scheme is the Gear-predictor-corrector algorithm [491] which can be used in diﬀerent forms

with respect to the order of time derivatives used.

5.1 Introduction 457

A typical time step is δt =0.05 · t

unit

=1.84fs for Cu potentials. A proper

choice of the used time step can be veriﬁed by checking conservation laws

within the utilized MD method like energy conservation and/or conservation

of the linear momentum of the overall system. From a physical point of view

thetimestepshouldbeatleastonlyafraction of the minimal period related

to the highest vibrational frequencies of the system. The Debye frequency, ν

can be considered as a good ﬁrst reference to assess this aspect. Since ν

the order of 10

/s, δT 

5.1.4 Boundary Concepts

Diﬀerent boundary concepts are utilized in atomistic systems depending on

the physical properties to be studied. Commonly periodic boundary concepts

are applied to simulation boxes to suppress ﬁnite size eﬀects of the atomistic

systems due to the limited number of atoms that can be considered in atom-

istic simulations. 3D periodic boundary conditions are common in atomistic

simulations of crystalline materials. They are most adequate for ideal crystal

simulations and for bicrystal simulation boxes where a strong interaction of

the adjacent GBs, which exist in such systems due to the periodicity, can be

ruled out.

For pure twist GBs 3D periodic boundary conditions can be used when

the GBs spacing is suﬃciently large and no massive GB sliding occurs. For

low-angle twist GBs GB sliding becomes increasingly important so that a 3D

periodic border concept may not be the best choice. All tilt GBs have a strong

tendency for GB sliding.

If sliding at a GB becomes an issue, 3D periodic boundary conditions may

give rise to artiﬁcal, unphysical interactions of the two GBs. In such case a free

surface normal to the GB would oﬀer an alternative if the used interatomic

potential is able to represent the free surfaces realistically and an interaction

of the free surfaces and the GB can be neglected. Another concept is the

frozen block approach [507] (Fig. 5.2) or its variants [524].

5.1.5 Finite Temperature MD Simulations

To perform ﬁnite temperature MD simulations one has to ﬁrst deﬁne whether

the temperature T is a thermodynamic variable or not. In the case of mi-

crocanonical simulations, T is not a thermodynamical variable, but an MD

system will equilibrate to a certain temperature, which depends on the initial

choice of kinetic energy introduced to the system [491].

In a canonical ensemble the temperature is a thermodynamic variable be-

sides the number of atoms N , and the volume of the simulation box V .They

are kept thermodynamically ﬁxed during the complete course of a simulation

run. Two major MD methods exist, namely an isokinetic method [491] which

does not really represent a canonical ensemble, and ﬁnite temperature ther-

458 5 Computer Simulation of Grain Boundary Motion

Region I

Region II

3D frozen upper block

3D frozen lower block

initial position of the

3D frozen upper block

initial position of the

3D frozen lower block

FIGURE 5.2

Region I/II 3D frozen block border concept according to Lutsko and Wolf.

Here the blocks of region II are rigid frozen atomistic blocks which are shifted

about the translations vectors, t

and t

respectively, according to the net

force exerted by region I on either of the region II blocks.

5.1 Introduction 459

Region I

Region II

upper block with

frozen s

eqs. of motion

lower block with

frozen s

eqs. of motion

FIGURE 5.3

Restricted free surface border concept. Here the free surfaces are stabilized in

the normal direction of the surfaces since during performed MD simulations

the s

component of region II atoms is held ﬁxed. According to this procedure

even for simple pair potentials stable surfaces can be obtained even in the

presence of a DF.

460 5 Computer Simulation of Grain Boundary Motion

mostats [491, 493].

In the case of isokinetic simulations, the desired temperature T can be es-

tablished by rescaling the atoms’ velocities after each time step according to

the ratio of the desired kinetic energy in the overall system to the actual ki-

netic energy in the overall system [491], which, however, does not represent

the canonical ensemble [491]. Hoover [498] and Nose [494, 497] introduced a

diﬀerent approach which utilizes frictional forces, a kinetic term for the vari-

able associated with the heat bath and a logarithmic potential function to

introduce the heat bath.

5.2 Driving Force Concepts

5.2.1 Introduction

GB migration is caused by a driving force. The introduction and quantiﬁcation

of driving forces for GB migration play a key role for the atomistic modeling

of the migration of planar GBs. The study of planar GBs oﬀers the advan-

tage that the grain boundary structure is unambiguously deﬁned whereas in

curved grain boundaries the grain boundary character changes with location.

However, planar GBs do not experience a curvature driving force.

To exert a driving force (DF) on a planar GB an energy gradient across

the boundary has to be introduced into the system. In the following three

diﬀerent forms of DFs will be introduced to drive planar GBs. Bicrystal sim-

ulation boxes will be used exclusively in the simulations. The ﬁrst approach

utilizes the anisotropy of elastic constants which will cause an elastic energy

diﬀerential across the boundary upon elastic loading [508]. A second approach

makes use of a Peach-Koehler force owing to the coupling of an external stress

state with the GB structure itself. A third DF scheme introduces an excess

energy term to the orientation dependence of the atoms.

5.2.2 Elastic Driving Force

To introduce an elastic driving force one has to apply an appropriate strain or

stress to a bicrystal simulation box that will cause an elastic energy diﬀerence

in the two grains. Such an approach will only work in materials which behave

elastically anisotropically since then the orientation will determine the elastic

energy of a grain. Furthermore, only certain strain or stress states will give

rise to an elastic driving force on a GB. For individual planar GBs such states

need to be identiﬁed. This can be achieved by linear-elastic energy calculations

for distinct GB geometries. Due to the gradient of stored elastic energy across

a GB, the GB will migrate toward the grain with the higher elastic energy.

The elastic driving force itself is then simply given by the diﬀerence in the

5.2 Driving Force Concepts 461

stored elastic energy density of the two grains. From the orientation of each

grain with respect to the coordinate system (x, y, z)ofthesimulationbox

the rotation matrix can readily be derived that relates each grain orientation

to the principal coordinate system of the crystal structure studied. These

rotation matrices can then be used to transform the applied strain or stress

tensor into the principal coordinate system of grain A and grain B in order to

calculate the elastic energy density of each grain. The reference state of a cubic

material in its principal coordinate system is deﬁned by the situation where

the 100 directions of the crystal are aligned with the coordinate system that

deﬁnes the atom coordinates within the grain and where the primitive cell of

the system is given by a cube as illustrated in Fig. 5.4(a). The misorientation

relationship between the MD simulation box coordinate system (x, y, z)and

those of either grain, A and B, are illustrated in Fig. 5.4(b) for the case of a

[001] planar twist GB.

We will conﬁne our consideration to cubic materials. In the following the

elastic driving force for a cubic crystal structure will be derived. The elastic

energy density for a cubic crystal in its principal cubic coordinate system

reads

elast

ijkl

1111

(ε

+ ε

)+2C

4444

(ε

+ ε

1122

(ε

· ε

+ ε

· ε

+ ε

· ε

)

ijkl

1111

(σ

+ σ

)+2S

4444

(σ

+ σ

1122

(σ

· σ

+ σ

· σ

+ σ

· σ

), (5.20)

where the Einstein summation convention is implied. Transformation of the

strain or stress tensor into the principal coordinate system of grains A and

B, and a calculation of the diﬀerence between the elastic energy densities in

both grains, A and B, ﬁnally yields the elastic driving force p

elast

= E

elast

− E

elast

(5.21)

ijkl

(ε

− ε

ijkl

(σ

− σ

)

For [001] twist or tilt GBs this reads

elast

= E

elast



[001] , +



− E

elast



[001], −



(5.22)

ijkl

(ε

− ε

)

=sin(2θ)

: ;8 9

(2C

4444

− C

1111

+ C

1122

) ε

(ε

− ε

)

ijkl

(σ

− σ

)

=sin(2θ)(2S

4444

− S

1111

+ S

1122

)

9: ;

(σ

− σ

)

462 5 Computer Simulation of Grain Boundary Motion

if the [001] axis is aligned with the z-direction of the simulation box coordinate

system. The angle θ represents the rotation angle about the [001] axis. This

result implies a choice of the strain tensor as follows:

⎛

⎝

000

⎞

⎠

, p =(2C

4444

−C

1111

1122

)·ε

·(ε

−ε

)·sin(2θ)

(5.23)

where C

ijkl

are the elastic constants. To apply a coupled xy shear state seems

reasonable since no further components of the strain tensor do appear in

the elastic driving force equation (5.23). Hence, any additional strain tensor

component would not result in a change of magnitude of the elastic driving

force. To further simplify the situation and to avoid a volume change of the

simulation box, the following strain tensor can be chosen for [001] twist GB

migration simulations:

⎛

⎝

ε 0

ε −

000

⎞

⎠

, leading to p = −(2C

4444

−C

1111

+ C

1122

)·ε

·sin(2θ)

(5.24)

It is important to realize that the simulations can be performed either

strain or stress controlled. In some cases it is crucial to chose appropriate

boundary conditions for observing continuous GB migration. From a practical

point of view, the strain-controlled MD simulations are robust to set up, since

the MD simulation box remains ﬁxed and does not need to be adapted to

an imposed external stress state. From a physical point of view, the stress

controlled simulations are more favorable since the simulation box can adjust

its size and shape to the applied external stress tensor. This is a very important

feature for ﬁnite temperature GB simulations if one pays attention to the fact

that GBs can change their structure and character with temperature.

5.2.3 Peach-Koehler Driving Force

The well-established Peach-Koehler force exerted on single lattice dislocations

due to internal or external stress states can also be used to study the eﬀect

of external stresses on GBs. This implies that the structure of such GBs still

contains discrete GB dislocations. It is well known that low-angle GBs are

composed of discrete GB dislocations. Their spacing D is large but decreases

with rising misorientation angle Θ according to sin(Θ/2) = b/2D where b is

the Burgers vector. Frequently, high-angle GBs are assumed to have a disor-

dered structure, but in principle, any grain boundary can be represented by

a superposition of primary and secondary grain boundary dislocations.

It has been shown experimentally [511] that low-angle GBs can be driven

by Peach-Koehler forces. The force F

Disl.

, exerted on a single dislocation of

5.2 Driving Force Concepts 463

simulation box :

h = (a,b,c)

(a)

y || [010]

x || [100]

z || [001]

twist

z, z

A,B

|| [001]

(b)

FIGURE 5.4

Geometry of MD simulation box: (a) shows a schematic of a MD simulation

box of a perfect crystal while (b) presents a bicrystalline MD simulation box

containing two grains labelled A and B for the case of an arbitrary [001]

twist GB. The direction of grain rotation and the primitive cells for a simple

cubic material are illustrated in the ﬁgure as well. Note that both grains are

symmetrically misoriented with respect to the MD simulation box coordinate

system (x, y, z). The box vectors a, b and c that span out the simulation

box and which are essential to apply any form of periodic border conditions

are given in the ﬁgures as well. More details on the general aspects of the

generation scheme and the geometry of the box are given in the text.

464 5 Computer Simulation of Grain Boundary Motion

length L, Burgers vector b

Disl.

and line vector s

Disl.

by a stress state σ at the

dislocation is given by

Disl.



· b

Disl.



× s

Disl.

(5.25)

Recently it was shown experimentally that even high-angle tilt GBs can be

driven by an external stress [512]. The simulations show that appropriately

chosen pure shear stress states will result in Peach-Koehler forces on pure

edge GB dislocations at least in the case of low-angle tilt GBs.

5.2.4 Orientation-Correlated Driving Force

To avoid artifacts due to the coupling of the stress state with the grain

boundary structure an artiﬁcial driving force for GB migration can be used

[396, 524]. The driving force itself is only artiﬁcial in the sense that it does

not specify its source but can be considered like a magnetic or electrical driv-

ing force in materials with anisotropic magnetic susceptibility or dielectric

constant. An atom of a speciﬁc crystal will experience an anisotropic excess

energy due to an externally applied ﬁeld. Therefore, one can construct an ex-

cess energy term for each atom as a function of the structure factor |S(k

associated with the orientations present in the system and which are calcu-

lated in the vicinity of each atom. The vicinity of an atom is considered a

sphere of a certain radius centered around the atom. The radius will be the

same as the one used to evaluate the interatomic forces on an atom, or its

potential energy. Thus, for each atom j an excess energy term, e

corr

,canbe

obtained by a simple analytical expression. The sum over all atoms can then

simply be added to the potential energy of the overall system. By doing so this

orientation-correlated driving force (OCDF) term will enter the Lagrangian

or Hamiltonian and, hence, new equations of motion for each atom will be

obtained. The OCDF will be identiﬁed by p

corr

In the following the respective equations of motion will be derived. The

structure factor of an orientation α represented by a respective reciprocal

lattice vector k

is given by

|S(k



m=1



l=1

exp (i · k

· (r

− r

))

⎛

⎝





m=1

cos(k

· r

)







m=1

sin(k

· r

)



⎞

⎠

(5.26)

The overall excess energy contribution to the system E

corr.DF

,aswellasthe

excess energy e

corr

for each atom j in a bicrystal consisting of an orientation,

5.2 Driving Force Concepts 465

1 and 2 is, therefore, given by

corr.DF



j=1

corr



j=1



α=1

·|S(k



j=1

·|S(k

+ a

·|S(k

(5.27)

where a

and a

are coeﬃcients which determine the strength of the excess

energy due to an externally applied ﬁeld, and are chosen such that the magni-

tude of the driving force is represented, for instance, 8.0 ·10

−3

eV/atom, and

2.0 ·10

−3

eV/atom, respectively. In order to evaluate Eq. (5.27), the value N

in Eq. (5.26) will be replaced by N

which is the number of neighbors within

the cut oﬀ sphere around atom j.

The concept of an OCDF can be introduced into the utilized Parrinello-

Rahman method as follows, e.g. using pair potentials. The adapted Lagrangian

now reads

L({s

}, {˙s

}, {h}, {

h})=K({

}, {h}) − V ({s

}, {h}) − E

corr.DF

({s

}, {h})

+ K

box

({

h}) − V

elastic

({h}) (5.28)

Utilizing the Lagrangian equation of motion for the coordinates s

i,x

i,y

and

i,z

one obtains new equations of motion with the incorporated OCDF

= −



j=i

dV (r

)

·s

−G

−1

corr.DF



}, {h}



−m

·G

−1

G·˙s

(5.29)

There are some advantages of the OCDF over the elastic DF. Both DFs share

a temperature dependence. The temperature dependence of the OCDF for

ﬁxed coeﬃcients a

and a

, is equivalent to the one for an elastic DF for a

ﬁxed strain, namely it decreases with increasing temperature. This behav-

ior is associated with increasing thermal vibrations with rising temperature,

which cause the structure factor to decrease. Accordingly, the temperature

dependence of the OCDF is realistic. The main advantage of the OCDF over

the elastic DF, however, is that now elasticity and/or processes that lead to

plasticity can be separated from the DF.

Of course, the calculation of the additional force contributions are more de-

manding and complex compared to a standard pair potential procedure and,

therefore, computation time goes up.

Even though the computational routines are more complicated and require

more computation time, this approach has many advantages over the elastic

DF approaChapter In order to obtain noticeable grain boundary migration

in the short time frame of nanoseconds a large DF has to be applied to the

simulation box, i.e. large strains when introducing an elastic DF. This may

cause artifacts in grain boundary behavior and is limited by the theoretical

shear stress. In case of an OCDF the driving force can easily be set to any

desirable value by choosing appropriate coeﬃcients a

and a

, and this is a

convenient way to overcome unphysical elastic strain magnitudes.