32
33
Subsequently, in order to use this potential in the MD simulation, it is needed to find the interatomic forces. The forces are given by the gradient of the potential energy. Since the potential depends on the interatomic distances only, the force on atom i can be expressed as
ij ij i
j ij
j ij
i
i r
r r E r F E
v v
∑
≠ ⎟⎟⎠
⎞
⎜⎜
⎝
⎛
∂ +∂
∂
= ∂ , (4-2)
or
∑
≠=
i
j ij
ij ij
i r
F r F
v v
, (4-3)
where
ij j ij
i
ij r
E r F E
∂ +∂
∂
= ∂ . (4-4)
The expression for Fij is analytic and given with the notations introduced in Eq. (4-5).
[ ]
⎥⎥
⎦
⎤
⎢⎢
⎣
⎡ ⎟⎟
⎠
⎜⎜ ⎞
⎝
⎛ −
− +
⎥+
⎥⎦
⎤
⎢⎢
⎣
⎡ ⎟⎟
⎠
⎜⎜ ⎞
⎝
⎛ −
−
−
= 2 exp 1 − (2 ) − (2 ) exp 2 1
0 0
0 0
12 12
r q r q
r q q r
p r r
Fij Ap ij ξ φi φi ij
(4-5)
where
∑
⎥⎦
⎢ ⎤
⎣
⎡
⎟⎟⎠
⎜⎜ ⎞
⎝
⎛ −
−
=
l
kl
k r
r 1
exp )
(
0
λ λ
φ
The term φk is the contribution of surrounding particlesstems from the N-body character of the potential. The parameters of tight-binding potential used in this dissertation are listed in Table 4-1.
4-1-2 Simulation model and conditions
In the simulation model of ARB process, a bi-layered structure consisting of hexagonal closed-packed (HCP) Zr coupled with the face-centered cubic (FCC) Ni elemental layers,
34
measuring 20 nm in width, 20 nm in length and 12 nm in thickness for each elemental layer, was set to be the starting alloy model in order to trace the final cyclic transformation stage between nanocrystalline and amorphous phases during ARB. Meanwhile, the other model consisting of HCP Zr and HCP Ti is set up under identical condition with above case for comparison. The designed amounts of Zr and Ni atoms are 89,280 and 188,442 for the Zr-Ni system (Zr32Ni68), while 116,064 Zr atoms and 118,184 Ti atoms are designed for the Zr-Ti systems (Zr50Ti50). The processing ARB speed in the MD simulation is 0.05 nm/fs or 2×1012 s-1. The present simulation employs the verlocity verlet algorithm to calculate the trajectories of the atoms and the scaling method is adopted during the simulation to control system temperature at 300 K. The repeated ARB strain-and-stack procedure applied in the simulation is schematically illustrated in Fig. 4-1.
In the consideration to the computing efficiency caused of data communication in a large scale, the distributed computing is the straightforward to implement parallelization method into MD program[155]. There are two common parallel algorithms for MD, one is particle decomposition method, and the other is spatial decomposition. In the particle decomposition method, all particles in the simulation cell are assign numbers, partitioned, and distributed to processors, as shown in Fig. 4-2 (a) [155]. The allotted particles are governed by each processor from the beginning to the end of the calculation. The globally data such as positions of particles must be shared and stored into each processor, which will decrease scalability and calculating efficiency in the system with increasing the number of processors, resulted from the synchronous global communication of broadcasting. Nevertheless, easily adjusting load balance among processors by uniform assigning the same number of particles to each processor and simply program writing are thought as advantages.
In the spatial decomposition method, the particles are divided locally into the partitioned
35
space of simulation cell belonging to each processor, as depicted in Fig. 4-2 (b)[155]. In this method, only the communication data of particles near or across the boundary region are necessarily exchanged among processors, the others are still restricted to local areas.
Therefore, the scalability and communication cost are significantly better than the particle decomposition method. However, load imbalance will come up as extremely inequality of the distribution of particles. In this thesis, the particle decomposition method is adopted for executing MD simulation, since the change of modeling shape is required during the deformation process.
4-1-3 Analysis methods for structural properties
In the current MD simulation, one of the most common methods to measure the structure features of a system is the radial distribution function (RDF), g(r). Also, the RDF is common applied to neutron and X-ray scattering experiments on simple fluids or amorphous structures. The RDF plays a central role in theories of the liquid state. Numerical results for the RDF can be compared with theoretical predictions and then offer a criterion to test a particular theory [140]. Measuring the RDF in the MD simulation is straightforward by considering the ratio between the average number density n(r) at a distance r from any given atom and the density at a distance r from an atom in an ideal gas at the same overall density.
The equation of RDF [156] can be calculated as
∑
= Δ= n
i r r
r n N
r V g
1 2
2 4
) ) (
( π , (4-6)
where N denotes the number of atoms in the simulation cell, V is the volume of the same cell, and n(r) the number of particles which can be found in the shell from r to r＋Δr. For a binary alloy system, partial radial distribution function (PRDF) for atom α and atom β is calculated by [156]
36
∑
= Δ= α
π
β β
α αβ
N
i i
r r
r n N
N r V g
1 4 2
) ) (
( . (4-6)
Subsequently, a powerful tool proposed by Honeycutt and Anderson [75] is used to study the local symmetry of the atomic arrangement within short-range ordering during the transformation process in the simulation model. Honeycutt and Anderson had used the pair analysis technique (termed as the HA pair index) to study the local structure features in disordered systems such as liquid and glass in early years [75]. This approach provides very clear information about the local symmetry of atomic arrangement more than the radial distibution function, and has been used to simulate the grain boundary microstructure transition, local cluster structure, and glass forming of metallic alloys under the rapid cooling condition [157, 158].
In this technique, two atoms are viewed as forming a bond pair if they are within a given cutoff distance that equals the first minimum in the partial radial distribution function (for example, Zr-Zr, Ni-Ni and Ni-Zr in the Zr-Ni system). There is a sequence of four integers to characterize the local structures according to the HA index. The first integer is to identify whether or not the atoms bonded in the HA pair are the near neighbors; 1 means yes, and 2 means not. The second integer is the number of the near neighbors shared by the HA pair. The third integer denotes the number of bonds among the shared neighbors. When the first three indices are identical but the bond geometries are different, the forth integer is added.
According to this method, different HA indices will represent different local structures, as shown some examples in Fig. 4-3 [159]. For example, the 1421 and 1422 pairs would exist predominantly in the close-packed crystalline structure, such as a face-centered cubic (FCC)
37
type (including alloys, such as L10) or hexagonal closed-packed (HCP) type (including the B19 lattice). For a FCC cluster, the predominant pairs would be the 1421 ones; while in the HCP cluster, there are typically around 50% of the 1421 and 50% of the 1422 pairs. The 1441 and 1661 pairs are characteristic of the body-centered cubic (BCC or B2) structure. The 1551, 1541, and 1431 pairs are referred to the common short-range local structures of an amorphous or liquid state. Finally, the 1321 is a packing related to the rhombohedra, and this pair tends to evolve when the icosahedra 1551 packing is formed and can be viewed as the side product accompanying the icosahedra atomic packing.
Furthermore, one can identify several important polyhedra with specific symmetries through the utilization of HA technique. A cluster can be easy identified as a FCC polyhedron, if its center atom only has twelve near-neighbors consisting of twelve 1421 bonded pairs.
Similarly, one can define a HCP polyhedron, if a center atom is surrounded by twelve near-neighbors but consisting of six 1421 and six 1422 bonded pairs. If a center atom has twelve neighboring atoms consisting of twelve 1551 bonded pairs, one can define an icosahedron. A defected icosahedron is composed of only 1551, 1441, and 1661 bonded pairs.
Figure 4-4 shows four types of the above-mentioned clusters detected by the HA method in crystals and glasses [158, 160].
4-2 Atomic simulation of vitrification transformation in Mg-Cu thin film
4-2-1 Effective-medium theory potential
In this section, the effective-medium theory (EMT) potential is employed in this work to study the Mg-Cu system. The EMT model used in this simulation corresponds to a parameterization of the EMT barriers calculated by Jacobsen et al [161, 162]. The basic idea
38
of his method is straightforward as that the total energy of a given atom in the system can be determined by the effect of their surrounding atoms. The energy of an atom in an arbitrary environment should be calculated by first estimating it in a well chosen reference system, namely the effective medium, and then by evaluating the energy difference between the real and reference systems. Thus the total energy of the system is written as
∑
+ −∑
=
i i
i c i
c E E
E
E , ( ,), (4-8)
where Ec,i is the energy of atom i in the reference system. A perfect FCC crystal is taken as the reference system in the present model. How to derive the form of the correction term
∑
−
i i
Ec
E , in Eq. (4-8) from the density function theory has been presented in the reference papers [161, 162].
The total energy of the system can be written as
∑
+Δ +Δ=
i
el AS
i i
c n E E
E
E ,( ) 1 . (4-9)
The first term Ec,i(ni) is the cohesive function of system and ni is the embedding density of atom i. The second term ΔEAS is atomic-sphere correction which is the difference in electrostatic and exchange-correlation energy for the atoms in an interested system and in the reference system, and ΔE1el is one-electron correction that is the difference in the sum of one-electron energies in the two systems. The atomic-sphere correction ΔEAS can be approximated well by a difference of pair potential in the two systems, and the considerable effects of ΔE1el can be included in the adjustment of pair potential form Vij. Thus the Eq. (4-9) can be expressed as
[ ]
∑
+Δ=
i
AS i
i
c n E i
E
E ,( ) () (4-10)
∑ ∑ ∑
⎪⎭
⎪⎬
⎫
⎪⎩
⎪⎨
⎧ ⎥
⎦
⎢ ⎤
⎣
⎡ −
+
=
≠ ≠
i j i
ref
i j
ij ij ij
ij i
i
c n V r V r
E,( ) 21 ( ) ( ) . (4-11)
39
The embedding density ni describes the superimposing contributions to the electronic density at the site of the atom i from all atoms j, and is expressed the function of intermolecular distance rij and Wigner-Seitz sphere radius si (the same as neutral sphere radius).
∑
≠Δ
=
i j
ij i j
i n s r
n ( , ) (4-12)
The neutral sphere radius si is obtained by solving Eq. (4-12),
⎟⎟⎠
⎜⎜ ⎞
⎝
− ⎛
= 1 log 121,
2 0
i
i s
s σ
βη , (4-13)
and
[ ]
∑
≠−
−
=
i j
ij
i exp 2(r s0)
,
1 η β
σ . (4-14)
The cohesive function Ec,i(ni) can also be simplified as the function of Wigner-Seitz sphere radius si
[
( )]
)
(s E0f s s0
EC i = λ i− , (4-15)
) exp(
) 1 ( )
(x x x
f = + − , (4-16)
where s0 is the constant of the equilibrium neutral-sphere radius and E0 is the cohesive energy.
Finally the last term in the Eq. (4-10) can be written as Eq. (4-17) since there are 12 nearest neighbors at a distance γ = βsi determined by the radius si in the FCC reference system,
⎥⎦
⎢ ⎤
⎣
⎡ −
=
Δ
∑
≠i j
i ij
AS i V r V s
E ( ) 21 ( ) 12 (β ) . (4-17)
Note that β is a geometric factor approximating to(16π/3)1/3/ 2≈1.81. The pair potential in Eq. (4-17) is parameterized as Eq. (4-18),
[
( / )]
exp )
(r V0 r s0
V =− −κ β − . (4-18)
All the parameters mentioned above can be determined from experimental properties such as lattice constants, cohesive energies, as well as elastic constants for the interested metals and their figures are listed in Table 4-2[161, 162]. For two-component metal systems (A, B), the neutral radius si of a type A atom in Eq. (4-13) should be slightly adjusted to join the effects
40
of additional type B atoms by a weighting coefficient χAB
) 12(
log 1 1
1 1
1 1 2
0 AA AB AB
A B A A
A s
s σ χ σ
η η
βη + − +
−
= , (4-19)
where σ1ΑΑ is the sum of Eq. (4-14) running over the A neighbors to the atom A and σ1ΑΒ is the sum over the B neighbors. The coefficient χAB is given by
) )(
( 0
0 1 1 0
0 1
0 1
A A A B A A
B
B s s
s A
s B
AB e
e n
e
n − −
−
= −ηη η η
χ . (4-20)
The last exponential function in Eq. (4-20) is usually very small so that it can be taken to be 1 although the neutral radius si is contained at the same time. Similarly, the atomic sphere correction for a type A atom is written
⎭⎬
⎫
⎩⎨
⎧ + −
=
Δ
∑ ∑
≠i ≠
j j i
iA AA ij
AB AB ij
AA A
AS
A B
s V r
V r
V
E 21 ( ) χ ( ) 12 (β ) . (4-21)
This EMT potential form has previously been shown to give a reasonable overall description of the FCC metals including transition and noble metals as well as their alloys.
They have been great used in the study of mechanical properties of crystalline metals [163].
For the Mg-Cu system, Bailey et al. [164] refit the parameters of both elements not only from basic properties of the pure elements but also considering the formation energies of two intermetallic compounds, Mg2Cu and MgCu2. The reasonable formation energies are important for the amorphous alloys, otherwise the system will simply separate into regions of pure Cu and Mg, respectively.
4-2-2 Derivation of force for EMT potential
As the above-mentioned in dealing with tight-binding potential form, it is needed to find the forces in order to use this potential in the MD simulation. The forces are given by the
41
gradient of the potential energy. Since the potential also depends on the interatomic distances only, the force on atom i can be expressed as
∑
≠=
i
j ij
ij ij
i r
F r F
v v
, (4-22)
where
ij j ij
i
ij r
E r F E
∂ +∂
∂
= ∂ , (4-23)
and
{ }
∑
+Δ∂
= ∂
∂
∂
i
A AS i
i c ij ij
i E n E i
r r
E ,( ) () , (4-24)
ij ij
ij c
r x x x
r E x E f r E
∂
− ∂
−
∂ =
= ∂
∂
∂ ( ) exp( )
0
0 , (4-25)
where the x is given by
) 12(
ln 1 )
( 0 1AA AB 1AB
A A A
A
A s s C
x=λ − =−λ σ +χ σ . (4-26)
The constant CA is assumed that βη2A＋η1B－η1A = CA. Assuming the calculated atom i is a type A species and the species of its neighbors j atom is the same with A atom (e.g. i = A, j = A), thus the distance between atom i and atom j is marked rijA, else is rijB. By algebra and separation of variables from Eq. (4-26), we obtain
A
A ij
AA AB AB A
AA A
A
ij C C C r
r x
∂
∂ +
+
= −
∂
∂ 1
1 1
1 ( )
σ σ χ σ
λ , (4-27)
or
A
A ij
A A
ij r
s r
x
∂
= ∂
∂
∂ λ , (4-28)
where C1 = η1B－η1A.
The derivative of the sphere radius sA is given as
A
A ij
AA AB AB A
AA A ij
A
r C
C C
r s
∂
∂ +
+
= −
∂
∂ 1
1 1
1 ( )
1 σ
σ χ
σ , (4-29)
where
)]
(
exp[ 2 0
1 2
A ij
A A
ij
AA r s
r A
A
β η
σ =−η − −
∂
∂ . (4-30)
42
Similarly for the case of i atom = A species and j atom = B species, the differential equation of x in Eq. (4-28) becomes
B
B ij
A A
ij r
s r
x
∂
= ∂
∂
∂ λ , (4-31)
where
B
B ij
AB AB AB AB A
AA A ij
A
r C
C C
r s
∂
∂ +
+
= −
∂
∂ 1
1 1
1 ( )
1 χ σ
σ χ
σ , (4-32)
and
)]
(
exp[ 2 0
1 2
A ij
A A
ij
AB r s
r B
B
β η
σ =−η − −
∂
∂ . (4-33)
Finally, the differential equation of atomic sphere correction for the case of i = A, j = A is given by
∑
∑
≠ ≠ ∂− ∂
∂ + ∂
∂
= ∂
∂ Δ
∂
i j
iA AA ij ij
AB AB i ij
j
ij AA ij
ij A AS
B A
B A A
A A
A
s r V
r r V
r r V
r
E ( ) χ ( ) 12 (β ). (4-34)
The first term by the relation in Eq. (4-34) is
⎥⎦
⎢ ⎤
⎣
⎡− −
∂ =
∂
∑
≠
) (
exp )
( 0 A A ij 0A
i j
ij AA ij
s V r
r
r V A
A
A A
β β κ β
κ , (4-35)
the second term is
∑
∑
≠ ∂ ≠= ∂
∂
∂
i j
ij AB ij
AB i
j
ij AB AB
ij B
B B A
B A
r r V
r
r χ V ( ) χ ( ), (4-36)
and the last term is
A
A ij
A A A A A
A iA AA
ij r
s s s V
s
r V ∂
− ∂
−
∂ =
∂ (β ) 0 κ exp[ κ ( 0 )] . (4-37)
Similarly, the differential equation of atomic sphere correction for the case of i = A, j = B is rewritten as
∑
∑
≠ ≠ ∂− ∂
∂ + ∂
∂
= ∂
∂ Δ
∂
i j
iA AA ij ij
AB AB i ij
j
ij AA ij
ij A AS
B B
B A B
A B
B
s r V
r r V
r r V
r
E ( ) χ ( ) 12 (β ). (4-38)
The first term on the right side of the relation is quoted to zero; the second term is taken as
43
∑
∑
∑
≠ ≠ ∂ ≠+ ∂
∂
= ∂
∂
∂
i j
ij AB ij
AB i
j
ij AB ij
AB i
j
ij AB AB
ij B
B B B
B B B
B B
r r V
r r V
r
r χ V ( ) χ ( ) χ ( ), (4-39)
and the last term is given by
B
B ij
A A A A A
A iA AA
ij r
s s s V
s
r V ∂
− ∂
−
∂ =
∂ (β ) 0 κ exp[ κ ( 0 )] . (4-40)
4-2-3 Simulation model and conditions
In this study, an initial sandwich model is constructed by stacking an FCC Cu layer, an HCP Mg layer, and an FCC Cu layer along the z axis. The 2-D scheme of this model is illustrated in Fig. 4-5. Two interfaces between the Cu and Mg layers are set in the [111]
direction of FCC Cu and [0002] direction of HCP Mg along the z direction, respectively. The length scale in x, y, and z dimensions of an FCC Cu layer is 3.5 nm × 6 nm × 1.5 nm and of an HCP Mg layer is 3.5 nm × 6 nm × 3.5 nm. The numbers of Cu and Mg atoms in this sandwich model are 5,467 and 3,595, respectively. In addition to the z direction, the periodic boundary condition is employed on the x-y plane. Due to the lattice difference between Cu and Mg, the dimensions of the periodic boundary along the x and y directions are chosen to be large enough to minimize the mismatch effects at the Mg and Cu interface. In order to speed up the formation of a disordered state in the interface between two metals which are reasonably existing present in reality, an equal number of Cu and Mg atoms are artificially exchanged in the interface (see Fig. 4-5). After exchanging the Cu and Mg atoms, the sandwich model is annealed at 300 K for 5×106 time step by employing in the NVT ensemble and the equations of motion are solved by the velocity-Verlet algorithm.
4-3 Cyclic loading fatigue in a binary Zr-Cu metallic glass
4-3-1 Interatomic potential used in Zr-Cu
44
As so far, two of interatomic force field have been fitted to ZrCu, one is the tight-binding potential proposed by Duan et al.[156], the other is the EMT potential proposed by Paˇduraru [165] et a. Although the EMT potential of ZrCu provides a better precision on the description of the structure of the metallic glass than the former [165], the tight-binding potential provides a simple equation form to writing program and to save calculating efficiency. For that reason, the tight-binding potential is adopted again to model the mechanical properties of Zr-Cu in this part. The parameters of the tight-binding potential for Zr and Cu employed in the current simulation were proposed by Duan et al.[156], and the potentials have been shown to be in good agreement with the experiments.
4-3-2 Isothermal-isobaric ensemble and Nosé-Hoover Chain
On the other hand, a more complete ensemble of MD simulation, that is, the isothermal and isobaric (NPT) ensemble, is applied to control the temperature and stress to maintain a reasonable thermodynamic state during the MD-simulation process, as well as to bring the imposed tensile stress on the simulation model. The basic conception of NPT ensemble has already been introduced in Chapter 3. Nosé-Hoover NPT ensemble is known as one of the most popular NPT system in the MD. However, the Nosé-Hoover algorithm generates a correct canonical distribution for MD only when there is just one conserved quantity or there are no external forces and the center of mass remains fixed [140].
To alleviate the restriction for the Nosé-Hoover ensemble, Martyna et al. [166-168]
extended a new scheme which Nosé-Hoover themostate is coupled to a chain of thermostats, Nosé-Hoover chain. This modification of Nosé-Hoover dynamics has been shown to give a very good approximation to the canonical ensemble even in pathological cases. The
45
significant difference between the original Nosé-Hoover dynamics and new method is shown in Figs. 4-6 (a) and (b) [166-168]. In Fig. 4-6 (a), the trajectories of the harmonic oscillator and the projected distribution functions for the original method do not fill the space even if the initial conditions are changed, but the distribution functions for the new one have good approximations to the canonical results and the dynamics fills space, seeing in Fig. 4-6 (b).
As the chains of thermostats are increased, the results will be better.
The equations of motion for the positions and momenta that generate the isobaric-isothermal ensemble are [166-168]:
i i
i i r
W p
r• = mΡ + ε , (4-41)
i i
f i i
Q p W
p N
F d ⎟⎟ Ρ − Ρ
⎠
⎞
⎜⎜
⎝
⎛ +
−
=
Ρ• 1 ε ξ , (4-42)
where Nf is the degree of freedom for the number of particles in d dimensions. The symbols ξ, pξ, andQ are the coordinate, momentum and effective mass which are used to describe the thermostate in the Nosé-Hoover chain. The symbols, ε, pε and W, are the barostat variables similar to the version of thermostat.
W
V• = dVpε , (4-43)
( ) ∑
=
• = − + N Ρ −
i i
i f
ext p
Q p m N p d
p dV p
1 2
int ξ ε
ε , (4-44)
Q pε
ξ• = , (4-45)
( )
∑
=• = N Ρ + − +
i
f i
i N kT
W p p m
1
2 2
ε 1
ξ . (4-46)
In the equation, pext is the external pressure or imposed stress, and pint is internal pressure, which can be calculated during the simulation process from virial equation,
46
( ) ( )
⎥⎦
⎢ ⎤
⎣
⎡
∂
− ∂
⋅ Ρ +
=
∑ ∑
= =
N
i
N
i i i i i
V V dV r
F m r
p dV
1 1
2 int
,
1 φ
. (4-47)
The derivation of time-reversible integration scheme for solving the equation of motion numerically in the NPT ensemble is need to use Trotter expansion formula. The details of using this technique can be found elsewhere.
4-3-3 Simulation model and conditions
The current simulation model consisted of 47,424 atoms, in a mixture of 50% (atomic percent) Zr and 50% Cu in a box size of 9.5 nm × 9.5 nm × 9.5 nm. Periodic boundary conditions were employed in all three dimensions. The glassy structure was obtained by quenching the sample model at a cooling rate 5 K/ps (ps = 10-12 s) from 2,000 down to 300 K, followed by an equilibrium for 30,000 time steps. The glass transition was found to be about 650 K according to the volume versus temperature diagram. The average density is about 7.29 g/cm3, and the potential energy at this as-quenched state is about - 5.233 eV. Because the difference of densities between 1 K/ps and 10 K/ps is 0.096 %, the density is considered insensitive to the cooling rate of the Zr-Cu BMG [110].
Both the stress-control and strain-control modes are included in this simulation. The uniaxial stress-control mode was first applied along the Z axis to implement the tension and compression cyclic-loading fatigue tests. To induce the localization of deformation in the simulation model during the cyclic-loading process, the maximum applied stress amplitude, σmax, was chosen to be 1 and 2 GPa, and the strain level during each fatigue cycle was around 2.8% for σmax = 2GPa. Note that the 2 GPa stress amplitude might appear to be too high judging from the fact that the fracture tensile or compression stress for bulk Zr-Cu amorphous specimen may only be around 1.5-1.8 GPa. But because of the size effect [169, 170], the
47
stress level in micro- or nano-scale can increase to more than 2.5-3 GPa. The minimum stress, σmin, was set to zero, equivalent to complete unloading; the stress ratio σmin/σmax is equal to zero. The applied load was increased gradually at a rate of 0.5 MPa/fs (fs = 10-15 s) until the maximum tensile stress and then was held for 2,000 fs to reach a relatively equilibrium state, and vice versa. The tensile fatigue cycle is set to be 100 and stress in both x and y directions is maintained at zero during the whole processes. The period of one cycle is 12,000 fs. The applied cycling load is presented in Fig. 4-7 (a).
In order to further amplify the fatigue-induced deformation, the strain-control mode is also used. Several strain rates are selected, that is, at 5 × 109, 1 × 1010, and 2.5 × 1010 s-1; the maximum strain is set to 2.5%, 4%, and 10% in the z direction, respectively, and all minimum strains were set to be zero. The higher strain rate used in this MD work is due to the inherent nature of current MD simulations especially for cyclic loading work which usually needs very long time compared with the monotonic condition [171, 172]. On the other hand, the MD simulations of the transition metal alloys such as Zr-Ni, or Zr-Cu are generally insensitive to the strain rate effects [125, 173]. Again, the choice of 2.5 to 10%
strain might seem to be too high. But from literature reports [174], the metallic glasses in micro- or nano-scaled size can exhibit much higher plastic strain before failure.
Note that in order to recover the induced tensile strain back to zero, the induced load needs to be reversed to the compressive stress, so that the strain-control fatigue is actually under the tension-compression mode, different from the above stress-control fatigue under the tensile mode. When the applied strains reach to the designated values, the models are given a tolerable time steps for relaxation during the cyclic loading process. It was noted that the 2.5% maximum strain case did not show significant difference on the qualitative results from the above 2 GPa stress control mode. Thus the maximum strain was increased to 10% to
48
examine the possible more severe fatigue deformation under such a high strain amplitude.
The maximum stress achieved for this 10% strain-control mode approaches 3 GPa and the Poisson ration is about 0.43. The period of one cycle is about 80,000 fs. Since the strain accumulated during each cycle is much higher for this case, the fatigue cycle for this more severe 10% strain-control mode is set to be 25. The applied cycling-load pattern is presented in Fig. 4-7 (b). While the applied cyclic strain (10%) and the induced cyclic stress (3 GPa) appear to be exceeding the realistic levels for the Zr-Cu bulk metallic glass, the current simulation setup is simply to raise the cyclic strain and stress level to shorten the simulation cycle and time. The illustration applied in the simulation of cyclic fatigue is schematically illustrated in Fig. 4-8.
4-3-4 Observation of local shear strains
In order to quantify plastic deformation at the atomic level, there are three methods frequently used to identify where the irreversible plastic rearrangement are occurring. One is known as slip vector. The expression of slip vector is given in Eq. (4-48),
N r r v
N
i
ij ij sp
∑
=−
= 1
)2
' (v v
, (4-48)
where N is the total number of nearest neighbors of each atom i, and j is the nearest neighbor of atom i. rv and 'ij rv are the difference of two vectors between atom i and atom j, respectively, ij before and after deformation. Another is the atomic-bond rotation angle (ABRA) change [109]. The ABRA is quantified by calculating the change of angle ϕij between mutual orientations of the nearest-neighbor bond vectors. The angle ϕij is given by
49
⎟⎟
⎠
⎞
⎜⎜
⎝
= −⎛ ⋅
' cos 1 '
ij ij
ij ij
ij r r
r
ϕ r , (4-49)
where r and 'ij r are the nearest-neighbor bond vectors before and after deformation. While ij all atomic bonds rotate during shear deformation, the angle ϕij is directly related to the local atomic shear strain. The third method is the local Lagrangian strain [129, 140], ηMises. To calculate ηMises, it is required to seek a local transformation matrix Ji by minimizing,
⎟⎟⎠
⎞
⎜⎜⎝
⎛
⎟⎟⎠
⎞
⎜⎜⎝
=⎛
→
−
∑ ∑
∑
∈−
∈
∈ i i j Ni
ij T ij N
j ij T ij i
N j
ij i
ijJ r J r r r r
r ' '
1
. (4-50)
Therefore, each component of local Lagrangian strain ηMises is given by
(
JiJiT I)
i = −
2
η 1 , (4-51)
and
6
) (
) (
)
( 2 2 2
2 2
2 yy zz xx zz xx yy
xy xz yz Mises
i
η η η
η η
η η η η
η = + + + − + − + − . (4-52)
It is easy to find the fact that the local Lagrangian strain is the better than two formers in the measurement of local inelastic deformation. However, in addition to the extremes (just existing local normal strain or shear strain), there is almost no difference among the results of such three methods for the current case. Therefore, we present the ABRA results to demonstrate where the local shear strains appear under cyclic loading.
50