Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Modeling of multiscale continuum-atomistic systems using homogenization theory
(USC Thesis Other)
Modeling of multiscale continuum-atomistic systems using homogenization theory
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
MODELING OF MULTISCALE CONTINUUM-ATOMISITC SYSTEMS
USING HOMOGENIZATION THEORY
by
Karthikeyan Chockalingam
_ _
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(CIVIL ENGINEERING)
May 2010
Copyright 2010 Karthikeyan Chockalingam
ii
ACKNOWLEDGEMENTS
I would like to thank, Professor Carter Wellford for his advice, support and valuable
discussions during the entire period of this dissertation. It was because of his
continued commitment this work was completed successfully.
I would like to thank my committee member, Professor Aiichiro Nakano, for his
advice and extensive guidance in molecular dynamics simulation. I would also like to
thank my other committee members Professor Roger Ghanem, Professor Sami F.
Masri, and Professor Hung Leung Wong for their work on this dissertation.
Finally I would like to thank my family and friends for their continued support and
encouragement.
iii
TABLE OF CONTENTS
ACKNOWLEDGEMENTS……………………………………………………………….ii
LIST OF TABLES………………………………………………………………………..vi
LIST OF FIGURES…………………………………………………………………...…vii
ABSTRACT……………………………………………………………………………...xi
CHAPTER 1: INTRODUCTION………………………………………………………...1
1.1 Multiscale Approach
1.2 Modeling Approach – Macro Scale
1.3 Modeling Approach – Atomistic Scale
1.4 Methods for Linking Scales
1.5 Intended Applications
1.6 Objectives and Chapter Outline
CHAPTER 2: STATIC ANALYSIS: ATOMISITC – CONTINUUM
HOMOGENIZATION ………………………………………………………….……….13
2.1 Continuum to Spatial Discretization
2.2 Anharmonic Expansion about Equilibrium
CHAPTER 3: MOLECULAR MECHANICS: STATIC MULTISCALE
ALGORITHM………………………………………..…………………………………..30
3.1 Inter-Atomic Potential and Boundary Condition
3.2 Algorithm for Static Implementation
3.3 Computer Implementation
3.4 Results for One Dimension Lattice Cell
3.5 Results for a Two Dimension Triangular Lattice Cell Problem
3.6 Results for Static Multiscale Algorithm Implementation
CHAPTER 4: NON-EQUILBRIUM MOLECULAR MECHANICS, HEAT FLUX
AND VIRIAL STRESS …………………………………………………….…………...47
4.1 Non-Equilibrium Molecular Dynamics
4.1.1 Constant Heat Flux and Momentum Conserving Method
4.1.2 Constant Temperature Gradient and Momentum Conserved
Method
4.2 Heat Flux
4.3 Virial Stress
iv
CHAPTER 5: DYNAMIC ANALYSIS: ATOMISTIC – CONTINUUM……………...61
5.1 Dynamics Behavior in Macro and Discrete Fast Time Scales
5.2 Thermal behavior in Macro and Atomistic Scale
5.3 The Finite Temperature Problem
5.4 Results and Discussion
5.4.1 Problems with a Specified Temperature Distribution
5.4.1.1 Quasi-Static Thermal Problems
5.4.1.1.1 Quasi-Static Problems with Uniform
Temperature
5.4.1.1.2 Quasi-Static Problems with Non-Uniform
Temperature
5.4.1.2 Dynamic Thermal Problem
5.4.2 Coupled Thermo-Mechanical Model with Heat Transfer
CHAPTER 6: ATOMISTIC NON-LINEAR VIBRATIONS………………………….108
6.1 Dynamic Behavior in the Fast Time Scale
6.2 Solution for the Non-Linear Free Vibration Problem
6.3 Finite Temperature Amplitude Constraint
6.4 Incremental Temperature Algorithm
6.5 Test Problem
6.5.1 Calculation of Specific Heat at Constant Volume
6.5.2 Calculation of Thermal Expansion
CHAPTER 7: IMPLICIT TIME INTEGRATION ALGORITHM…………………....129
7.1 Weighted Residual Formulation
7.2 Approximation for Atomistic Force
7.3 Energy Preserving Algorithm
7.4 Result and Discussions
7.4.1 Stability of the Algorithm
7.4.2 Accuracy of the Algorithm
7.4.3 Additional Study
CHAPTER 8: SUMMARY AND CONCLUSIONS…………………………….…….152
8.1 Static Analysis
8.2 Temperature Dependence of Material Properties
8.3 Thermo-Mechanical Equation
8.4 Implicit Time Integration Algorithm
v
REFERENCES………………………………………………………………………....158
APPENDIX: Higher Order Expansion of L-J Potential..….…………………….…….162
vi
LIST OF TABLES
Table 3.1: Multiscale static algorithm………………………………….………………..34
Table 3.2: Assembly algorithm for
IJ
mn C
………………………………….……………...36
Table 5.1: Quasi-static thermo-mechanical multiscale algorithm………….…………..102
Table 6.1: Incremental temperature algorithm………….………………..….…………118
vii
LIST OF FIGURES
Figure 3.1: Lennard-Jones inter-atomic potential……………………..……………...….31
Figure 3.2: Molecular dynamics periodic boundary condition……….…………….……32
Figure 3.3: 1-D Point defect…………………………………………………………..….37
Figure 3.4: For fixed point defect ratio (L/r – 0.01)………………………………….….38
Figure 3.5: Equilibrium configuration of a perfect lattice – Energy/per atom = -2.9717
……………………………………………………………………………………………39
Figure 3.6: Equilibrium configuration of a point defect – Energy/per atom = -2.8523
……………………………………………………………………………………………39
Figure 3.7: Lattice under 2.5% strain in the Y – direction………………………….…...40
Figure 3.8: Lattice under 2.5% strain in the X – direction……………………….……...41
Figure 3.9: Energy/atom: Lattice under strain in Y – direction…………………….……41
Figure 3.10: Energy/atom: Lattice under strain in X – direction………………….……..42
Figure 3.11: Constitutive Constants: Lattice under strain in Y – direction……….……..43
Figure 3.12: Constitutive constants: Lattice under strain in X – direction………….…...43
Figure 3.13: 2-D Structural model consisting of Quad elements……………………..….44
Figure 3.14: Force displacement results for the macro model……………………….…..45
Figure 3.15: Gauss point evaluation of C1111 elastic constant…………………….……46
Figure 4.1: Argon 100 x 4 x 4 F.C.C structure…………………………………………..50
Figure 4.2: Non-equilibrium molecular dynamics Jund method………………….……..52
Figure 4.3: Constant temperature gradient and momentum conserving method…...……54
Figure 4.4: Lagrangian and Eulerian heat flux plots using NEMD methods……..……...56
Figure 4.5: Lagrangian, Eulerian and Zhou Virial stress plots………………….……….59
viii
Figure 5.1: One 3-D Hexahedron element model……………………………….……….93
Figure 5.2: Quasi-static strain versus dynamic time relaxation with uniform
temperature distribution….………………………………………………….94
Figure 5.3: Quasi-static stress versus dynamic time relaxation with uniform
temperature distribution….………………………………………………….95
Figure 5.4: Quasi-static displacement versus dynamic time relaxation with uniform
temperature distribution….………………………………………………….95
Figure 5.5: Four element constrained 3-D macro finite element model..………………..97
Figure 5.6: Quasi-static displacement versus dynamic time relaxation with
non-uniform temperature distribution………………………...…….……….98
Figure 5.7: Quasi-static stress versus dynamic time relaxation with non-uniform
temperature distribution….………………………………………………….98
Figure 5.8: Quasi-static strain versus dynamic time relaxation with non-uniform
temperature distribution….………………………………………………….99
Figure 5.9: Dynamic strain history plot with uniform temperature distribution…..……100
Figure 5.10: Dynamic stress history plot with uniform temperature distribution……....101
Figure 5.11: Dynamic distance history plot with uniform temperature distribution…...101
Figure 5.12: Quasi-static mid-plane nodal displace versus macro tome with
non-uniform temperature distribution……...……...……………………...105
Figure 5.13: Quasi-static nodal temperature versus macro tome with non-uniform
temperature distribution….……………………………………………….105
Figure 5.14a: Quasi-static nodal temperature versus molecular dynamic results with
non-uniform temperature distribution at 0 macro time……...…………...106
Figure 5.14b: Quasi-static nodal temperature versus molecular dynamic results with
non-uniform temperature distribution at 20 macro time…………..…….106
Figure 5.14c: Quasi-static nodal temperature versus molecular dynamic results with
non-uniform temperature distribution at 60 macro time…...…………….107
Figure 5.14d: Quasi-static nodal temperature versus molecular dynamic results with
non-uniform temperature distribution at 100 macro time..……………...107
ix
Figure 6.1: The evolution of starting frequencies (2.71, 5.69, 6.48, 7.21, 7.56) as a
function of temperature………….…………………………………………120
Figure 6.2: The specific heat calculated from Eq. (6.29) as compared to the MD
results obtained from Eq. (6.30)……………………………….…………..123
Figure 6.3: Typical frequency surface, involving data at the lattice constants 1.39,
1.49, 1.59, 1.64, and 1.71, presented in L-J units…..…………...………….125
Figure 6.4: Thermal expansion results obtained from the extended quasi-harmonic
model as compared to the experimental results obtained by
Peterson et al. [39]………………………………………………….……...127
Figure 6.5: Specific heat results obtained from the extended quasi-harmonic model
as compared to the experimental results of Peterson et al. [39]..…………...128
Figure 7.1: EPA Energy results at Density 1.1, Time Step 0.005 and Initial
Temperature 40K…...……………………………………………………...142
Figure 7.2: Velocity Verlet Energy results at Density 1.1, Time Step 0.005 and
Initial Temperature 40K………………..………………………………….143
Figure 7.3: EPA and Velocity Verlet Temperature results at Density 1.1, Time Step
0.005 and Initial Temperature 40K………..………………………………..143
Figure 7.4: EPA Energy results at Density 1.1, Time Step 0.01 and Initial
Temperature 40K……………...…………………………………………...144
Figure 7.5: Velocity Verlet Energy results at Density 1.1, Time Step 0.01 and
Initial Temperature 40K…...……………………………………………….145
Figure 7.6a: EPA and Velocity Verlet Temperature results at Density 1.1,
Time Step 0.01 and Initial Temperature 40K..……………………………145
Figure 7.6b: EPA at Time 0.005 and EPA at Time 0.01 Temperature results at
Density 1.1 and Initial Temperature 40K…………………………….…..146
Figure 7.7: EPA at Time 0.005 and EPA at Time 0.05 Temperature results at
Density 1.1 and Initial Temperature 40K………………………………….147
Figure 7.8: EPA and Velocity Verlet Temperature results at Density 1.1,
Time Step 0.005 and Initial Temperature 40K……………………………..148
x
Figure 7.9: EPA and Velocity Verlet atom displacement results at Density 1.1,
Time Step 0.005 and Initial Temperature 40K………………...…………..149
Figure 7.10: EPA and Velocity Verlet temperature results at Density 1.1,
Time Step 0.005 and Initial Temperature 80K…...…………………..…..150
Figure 7.11: EPA and Velocity Verlet temperature results at Density 0.8,
Time Step 0.005 and Initial Temperature 40K…….……………………..151
xi
ABSTRACT
The main objective of the dissertation is to develop multi-scale algorithms for
continuum-atomistic problems. The focus is on sequential multi-scale simulations. In
sequential multi-scale methods, the computations at the various scales are, in a sense,
decoupled. This means, for example, that, for a continuum/atomistic simulation, large
scale macroscopic continuum calculations rely on the results of fine scale computations
and information obtained on an atomistic cell.
While the procedures developed in this thesis could be used in conjunction with a
number of sequential multi-scale methods, the focus here is on the homogenization
technique. As has been the case in traditional finite element applications of
homogenization, one of the principal focuses in this thesis will be on the computation of
macro scale constitutive parameters; but, in this case, these constitutive representations
come from the atomistic calculations.
The thesis has four parts that develop various aspects of the theme of the work. The
dissertation focuses on the following applications:
1. Problems involving mechanical loading of solids and structures under
static load at zero temperature.
xii
The focus is on creating multi-scale continuum/atomistic simulation
methods which use the atomistic model to provide an improved material
representation including the effects of material defects. This topic could be
useful in modeling fracture and failure.
2. Computation of thermo-mechanical constitutive parameters at finite
temperature conditions.
This procedure focuses on using the atomistic scale calculation to define
constitutive parameters. It assumes that equilibrium conditions exist at the
atomistic scale. It does not attempt to track, in a time history sense, the
dynamics at the atomistic scale. It does require the solution of an atomistic
free vibration problem with natural frequencies dependent on temperature.
The procedure defines macroscopic thermo-mechanical constitutive
parameters, like the specific heat and the coefficient of thermal expansion,
as a function of temperature. These properties could be used directly in a
macroscopic continuum finite element model which would be valid at the
full range of temperatures.
3. Dynamic problems involving the simulation of the thermo-mechanical
behavior of systems at finite temperature, with and without heat transfer.
xiii
This procedure focuses on using the atomistic scale calculation to define
multi-scale, thermo-mechanical momentum and energy equations. It does
attempt to track, in a time history sense, the dynamics at the atomistic
scale. Energy equations are derived for both the scales based on first law
of thermo-dynamics. Two types of application problems are used to
demonstrate the theory. The first involves thermal stress analysis
simulation in which the temperature has no time variability and thus no
heat transfer occurs. The second involves simulations with time varying
temperatures and include heat transfer effects.
4. Implicit time integrations algorithms for atomistic momentum equations
that can be seamlessly coupled to macro models.
1
CHAPTER 1
INTRODUCTION
The principal objective of the dissertation is to create multi-scale computational
procedures for continuum-atomistic problems. In the following sections, the chosen
multiscale computational methods are characterized in terms of multiscale approach,
macro scale modeling procedure, atomistic scale modeling procedure, method for linking
the scales, uniform modeling procedure, and intended applications.
1.1 Multiscale Approach
In the case of continuum-atomistic simulations, the objective of the multiscale procedure
is to develop a bridge between scales to improve product quality by better characterizing
materials during the product design phase. The resulting refined material representations
potentially reduce product costs and improved performance.
Extensive research has been done in developing multiscale techniques in the past decade.
The concurrent multiscale methods address a class of problems wherein there is a direct
coupling between scales. Concurrent methods couple two different scales (continuum,
atomistic, quantum) to obtain multiscale solutions. One of the popular concurrent
2
multiscale techniques is the Quasicontinuum method [42], developed by the Tadmor and
Ortiz to describe the behavior of solids with defects. Other prominent concurrent
multiscale work has been by Broughton et al. [3], Wagner et al. [46] and Xiaoa et al. [50].
Typical application would be in the areas of modeling carbon nanotubes, crack
propagation and fracture mechanics.
There is another classes of models, known as decoupled multiscale methods, wherein
there is no direct coupling between scales. These methods are also known as sequential
methods. In sequential methods, the scales are decoupled, often meaning that the
underlying atomistic unit cell information is used to produce the macro material
representation. In decoupled methods, the atomistic calculation is done at a particular
spatial point, often a Gauss point in a finite element model, and the pertinent information
is transferred to the macro model. Chung et al. [8], Chen et al. [6] and Liu et al. [36] have
worked with this class of decoupled multiscale methods.
In this work, the decoupled multiscale procedure is utilized in all applications.
1.2 Modeling Approach – Macro Scale
It is necessary to have an analysis method for the continuum scale representation. The
finite element method (FE) [10] is one of the most widely used tools for the numerical
modeling of macro behavior, especially for the simulation of the performance of solids
3
and structures. The finite element method will be used in the continuum model associated
with this work.
1.3 Modeling Approach – Atomistic Scale
In certain critical applications, a refined and accurate material representation is
needed. These applications could involve, for example, structural fracture and failure. To
obtain an accurate material characterization for certain materials, atomistic models can be
creted and analyzed. Since these types of applications are the target of this work, it is
necessary that the multiscale analysis utilize a method capable of modeling the material
behavior at the atomistic (nano) scale. In this work, the atomistic model to be used will be
based on molecular mechanics and many of the ideas used in Molecular Dynamics (MD)
[1]. MD is a powerful numerical technique describing the atomistic behavior of materials
using Newton’s equation of motion. Though powerful, MD is computationally expensive
compared to macro simulation techniques such as finite elements. Thus the development
of multiscale continuum-atomistic technique is useful for the development of efficient
algorithms for certain classes of problems.
In terms of the products of the MD simulation activity, certain computed quantities are
important to this work. These include the Virial stress tensor [24] and the heat flux,
defined by the non-equilibrium molecular dynamics method (NEMD).
Most of the multiscale methods have demonstrated the working of algorithm limiting to
the usage of simple inter-atomic pair potential and not have adequately explored multi
4
body potentials. Though most researchers claim that their algorithm is applicable to
multi body potentials, it usually turns out to be a computationally unviable process. It
becomes important for researchers to address the issue of computational efficiency. Since
multiscale methods are expensive, it is important simulation happens in a realistic time
frame. This dissertation partly addresses this issue in proposing an efficient algorithm
that would be adoptable to different potentials without much overhead. There have been
some works done by Waismana et al. [45]
to improve numerical efficiency of multiscale
algorithms. Due to the availability of high performance computing [17] tools coupled
with efficient parallel algorithms [21] the simulation time can be significantly reduced.
1.4 Method for Linking the Scales
Scale linking
The critical activity of linking the scales varies with the multiscale method being utilized.
In the case of a decoupled multiscale model, especially one employing an atomistic scale
and the associated MD method, it is natural to select the homogenization method [9] for
the scale linking.
One of the principal assumptions in the homogenization method is that the underlying
microstructure is periodic. In solid mechanics, homogenization is used to obtain averaged
macroscopic constitutive equations. Homogenization seems particularly suitable for this
work that focuses on material constitutive characterization. Thus, homogenization will be
the primary method to link the continuum and atomistic scales in this work. The
5
motivation behind homogenization is to model the macro material properties for a
material that has periodic heterogeneous micro-structure. Guedes et. al. [22] coupled
continuum homogenization theory with adaptive finite element techniques to obtain the
homogenized material constants by using the underlying periodic micro-structure
information. For the case when the micro structure and macro structure are continua,
Tankano et. al. [43] used homogenization theory to solve problems under large
deformations, where the perturbed displacement field, due to the presence of the periodic
micro-structure, acts as a correction to the macro deformation gradient.
Multiscale techniques using homogenization theory have been developed by Chung et. al.
[7] and Chen et. al. [5] to obtain homogenized material constants.
Uniform Modeling Procedure
A uniform modeling procedure is important in linking the scales. This is especially true
when the macro scale involves a continuum and the micro scale involves a discrete
atomistic lattice. In this work a uniform weighted residual approximation procedure is
developed to deal with the two scales and their interaction. The continuous versus
discrete issue is handled by using both continuous and discrete weights in the weighted
residual algorithm.
6
1.5 Intended Applications
Zero Temperature Problems
Zero temperature problems can be resolved using molecular statics formulations
at the atomistic level. The emphasis will be on structural problems involving materials
with defects.
Finite Temperature Thermo-mechanical Problems
Multiscale atomistic, finite temperature problems remain one of the hardest problems to
be solved due to more variables being involved in satisfying the underlying conditions.
The Quasicontinuum finite temperature technique developed by Shenoy et al. [40] uses
Monte Carlo [32] simulation techniques with mixed atomistic and continuum scales to
study the equilibrium properties of defects. Wagner et al. [46] developed the bridge scale
decomposition techniques to bridge the atomistic and the continuum scales by
eliminating the extra degrees of freedom. There have been a number of applications of
multiscale modeling techniques to finite temperature problems. Fish et al. [15] have
introduced a finite temperature homogenization theory, resulting in a multi-scale thermo-
mechanical model.
In this work, a general homogenization framework for thermo-mechanical problems is
created. Both problems with and without heat transfer have been solved. In particular
thermal stress analysis problems, involving structures exposed to a set temperature field
7
and problems in which temperature changes occur because of heat transfer effects, are
considered for analysis.
In the thermo-mechanical simulation algorithms, there is no explicit calculation of the
material matrix, thus making algorithm more efficient for complex inter-atomic potentials.
An explicit time integration scheme was adopted in both the atomistic and continuum
momentum equations in order to avoid any force approximation. The issue of
simultaneously producing equilibrium at the disparate continuum and atomistic scales is
dealt with by using dynamic relaxation methodology which was also adopted in certain
cases to avoid any macro dynamic oscillation, thereby allowing for faster solution
convergence. Parallel computing techniques were used in the implementation although
there is no explicit discussion in the dissertation. The proposed thermo-mechanical
multiscale algorithm was tested to exhibit the heat transfer phenomenon in perfect Argon
crystals, and the results show good agreement with molecular dynamics.
Thermo-mechanical Constitutive Parameters
There has been a significant amount of work on the determination of thermo- mechanical
constitutive properties. Much of that work has utilized either the quasi harmonic method
or the local harmonic method [33]. Foiles [16] has studied the performance of these
methods, particularly focusing on the accuracy achieved in computing the free energy of
solids. Foiles has shown that the accuracy of traditional methods is significantly reduced
8
by the failures of the methods to properly model defects and anharmonic terms. Further
work has been carried out by Jiang et al. [30] who have developed thermo-mechanical
continuum relationships at finite temperature using local harmonic approximations and a
free energy formulation. Liu et al. [34] have employed a coarse-grained Helmholtz free
energy expression to create a multi-scale non-equilibrium thermo-mechanical simulation
tool.
In this work, the developed constitutive models, for example for the specific heat and the
thermal expansion, are computed utilizing an extended quasi-harmonic approximation.
The developed algorithm relies on procedures created to solve nonlinear free vibration
problems for continuous media [12, 20]. The present extension to the quasi-harmonic
method involves the explicit consideration of finite temperature and anharmonic effects.
These effects are introduced in the constitutive representations via the vibration
frequencies obtained from the inter-atomic potential. The focus of this work is on
defining an algorithm capable of determining the variation of the natural frequencies with
temperature, even at high thermally induced atomistic velocities. The developed
procedure can be characterized as an Incremental Temperature Algorithm (ITA). An
incremental loading type procedure is created by introducing a temperature amplitude
constraint, associated with the kinetic energy of the atomistic system. The resulting
algorithm allows the inclusion of the anharmonic terms in the potential and the tracking
of the nonlinear variation of the natural frequencies with temperature. Using the resulting
frequencies and statistical mechanics methods based on the Helmholtz free energy, the
9
thermo-mechanical constitutive parameters are computed, at any required temperature or
amplitude of velocity fluctuations.
Implicit MD Algorithms
The focus is on the key methodology used to integrate the Newton’s equation of motion
for the atomistic lattice, irrespective of the multiscale framework being utilized. A
majority of the time integration algorithms employed for finite element simulations are
implicit [27]. An algorithm is proposed to provide a consistent implicit time integration
algorithm for the atomistic problem as well. Velocity Verlet [44, 41] and Gear’s predictor
corrector method [19] are the widely used integration schemes for MD systems.
Performance of the developed implicit atomistic dynamics algorithm is compared to
Velocity Verlat results.
1.6 Objectives and Chapter Outline
In this thesis, as stated earlier homogenization methods are used to couple
atomistic and continuum scale to achieve the following results.
1. Develop a new homogenization formulation for the continuum-atomistic problem
and an associated weighted residual modeling approach.
2. To develop both static and dynamic versions of the multiscale weighted residual
algorithm.
10
3. To develop multiscale algorithms for static solid mechanics problems in which
the continuum constitutive equations are defined at the atomistic scale.
4. To incorporate finite temperature effects in the dynamic formulation. To develop
a continuum-atomistic thermal model.
5. To develop methods to determine constitutive parameters including the effects of
finite temperature.
Chapter 2 presents the formulation of homogenization theory to couple the atomistic and
continuum scale, by adding a correction to the Cauchy–Born [13] rules for the
deformation of the lattice. The micro solution obtained by solving the atomic cell
problem is used to determine the constitutive behavior of the material.
Chapter 3 reviews the basics of molecular mechanics and the particular inter-atomic pair
potential used to demonstrate the behavior of the lattice. The initial equilibrium
configuration of the triangular lattice with and without defect is presented. This is used as
the representative unit cell for solving problems in the static and finite temperature case.
The static case, including the cell under strain in X, Y directions and behavior of elastic
constants under strain, is presented.
Chapter 4 reviews the concepts of Non-Equilibrium Molecular Dynamics and goes on to
discuss the different types of Non-Equilibrium Molecular Dynamics currently used in
literature. The methodology used in computing heat flux from Non-Equilibrium
11
Molecular Dynamics procedures. Finally describes the concepts of Virial stress and its
interpretation in computing the atomistic stress.
Chapter 5 derives the dynamic equation of motion for the microstructure by fully
descretizing space and time scales. It incorporates the effects of finite temperature into
the formulations by using the first law of thermodynamics. Finally the thermo-
mechanical heat equations are presented.
Chapter 6 presents the formulation for the problem of predicating the lattice vibrations
resulting from finite temperature and defining the resulting changes in constitutive
parameters. The formulation uses harmonic approximation for the atom displacements. It
incorporates a higher order approximation of the potential in the hessien. A nonlinear
eigenvalue problem is set up to study the frequency dependence on temperature using an
incremental temperature formulation. Finally the results for the frequency dependence on
temperature for the lattice with and without defects are presented, and they can be used to
determine the variation of constitutive parameters with temperature.
Chapter 7 derives the atomistic momentum equation based on the Lagrangian equation of
motion. The atomistic force is approximated using a Taylor series expansion. An
unconditionally stable implicit algorithm is proposed. The results are compared and
discussed with respect to velocity Verlet algorithm.
12
Chapter 8 summarizes the objectives met in this dissertation and how some of the
proposed methodologies can be extended to other class of multiscale problems. Finally a
future study in the areas of multiscale polycrystalline structures is discussed.
13
CHAPTER 2
STATIC ANALYSIS: ATOMISTIC – CONTINUUM
HOMOGENIZATION
2.1 Continuum to Spatial Discretization
Let the upper case characters I, J, appearing as superscripts be atom number. Let
I
k
Y be the position vector of atom I in the Y
atomistic coordinate system. In the atomistic
coordinate system, we do not use a continuous integration over the
Y coordinate. Instead,
we replace integration with summation over the discrete atoms, considering the atom to
act like a discrete finite element. In this dissertation pair type potentials are considered. In
any “atom element”, a pair of atoms IJ is involved. Then the “continuous integration” is
replaced by a “finite sum” as follows:
11
11
NN
IJ
y
IJ
dy
YN
“Continuum” replaced by “atomistic” representation
In terms of the equilibrium equations, the weighted residual model should
produce 2 equations for each of the N atoms, producing 2N total equations. The
“displacement vector” for atom I in the
th
n direction (1,2) n
14
I
n u
(2.1)
The “displacement vector” for atom J in the
th
n direction is
J
n u
(2.2)
The “differential displacement” between “atom pairs” is defined by,
IJ J I
nn n uu u
(2.3)
The macro-displacement vector is continuous and is denoted
I
o
n
X
u
(2.4)
There are effectively two scales used, one is the continuous, one is discrete. The two
scales are,
Continuous Macro-Scale Discrete Atomistic-Scale
n X
I
n Y
( 1,... ) I N
To compare the two scales, we must evaluate the continuous scale X at the location of
atom I to form, the following displacement component:
I
nn
atom I
XX
15
The two scales are related by a small parameter
I
I n
n
X
Y
(2.5)
or
II
nn
X Y
Note that
1
I
n
I
n
Y
X
(2.6)
The relative positions of the atoms in the two scales is
IJ J I
nn n
X XX (2.7)
IJ J I
nn n
YY Y (2.8)
Typically the continuous macro-equation involves the evaluation at particular
macro-coordinate Gauss-point positions
. GP
n
X . At each Gauss Point position
. GP
n
X , a
discrete micro-coordinate system
I
n
Y is constructed.
16
The displacement vector for atom I is defined by asymptotic expansion in terms
of the parameter as follows:
1
....
II
Io I
nn n
XY
uu u
(2.9)
Here is the ratio of the typical co-ordinate dimension in the two scales.
Similarly for atom J, it can be seen that
1
....
JJ
Jo J
nn n
XY
uu u
(2.10)
From equation space (2.3), (2.9), (2.10), we obtain
1 IJ oIJ IJ
nn n uu u
(2.11)
where
J I
oIJ o o
nn n
X X
uu u
11 1 IJ J I
nn n uu u
Dividing (2.11) by
IJ
k
X , from (2.7) we obtain
1 IJ oIJ IJ
nn n
IJ IJ IJ
kk k
uu u
X XX
(2.12)
17
but approximately
,
IJ
n IJ
nk
IJ
k
u
u
X
(2.13)
and
,
oIJ
n oIJ
nk
oIJ
k
u
u
X
(2.14)
Introducing (2.13) and (2.14) in (2.12), the following expression is obtained:
1
,,
IJ
n IJ oIJ
nk nk
IJ
k
u
uu
X
(2.15)
Now modify (2.15) considering that
1IJ
u is a function of
I
n
Y
1
,,
IJ IJ
k n IJ oIJ
nk nk
IJ IJ
kk
u Y
uu
X Y
(2.16)
where,
1
IJ
k
IJ
k
Y
X
1
,,
IJ
n IJ oIJ
nk nk
IJ
K
u
uu
Y
(2.17)
18
Start with the hyper elastic continuous equilibrium equation as follows:
,, 0 0
oo
o
mnik n k i m
fW
Cu Wdv dv
vv
(2.18)
2
b
mnik
mi nk
E
C
FF
(2.19)
Where,
mnik C
- Material Matrix
b
E - The atomistic potential
nk F
- The deformation gradient
To create the multi-scale algorithms consider that u is discretely defined in terms
of atomistic positions and displacements. Consider that this definition is in terms of
atomistic pairs. Thus in going to a multi-scale algorithm, we do not use
1
y
dy
Y
for
averaging over a continuous micro-structure. Instead, we use a discrete summation over
the atomistic pairs IJ using
11
1
NN
IJ
IJ
N
as the averaging operation. From (2.18) the multi-scale version is derivable by first
expressing (2.17) in the discrete form, the material constitutive tensor presented as an
average over the atoms in the atomistic model is
19
11
1
NN
IJ
mnik mnik
IJ
IJ
CC
N
(2.20)
where,
IJ
IJ IJ IJ
ik mnik
mn
C YY
C
(2.21)
and
2
IJ b
mn
IJ IJ
mn
E
C
rr
(2.22)
or in another form
IJ
IJ IJ
i mnik
mnk
C Y
C
(2.23)
2
,
IJ
b
mnk
IJ
mnk
E
C
rF
(2.24)
And from (2.22) and(2.24), we note that
IJ
IJ
IJ
k
mn
mnk
Y
C
C
(2.25)
20
Using (2.20) to create a multi-scale version of (2.18)
0
,, 0 0
11 11
11
o o
NN NN
IJ
IJ IJ IJ IJ IJ
ik nk i
mn m
IJ IJ
IJ IJ
f
uW dv W dv YY
C
vv
NN
(2.26)
Special consideration must be devoted to the weight function
12
, , ,.....,
,
N
M
Wf XYY Y
Wf XY
Let,
,
II
I
II
XX Y
WW WX W
(2.27)
,
JJ
J
JJ
XX Y
WW WX W
(2.28)
then,
,,
JJII
IJ
T
XX Y X X Y
WW W
(2.29)
where,
,
,
J J
I I
J
XX Y
I
XX Y
WW
WW
Or using (2.27), (2.28) in (2.29), it be seen that
J I
IJ J I
T
WWX WX W W
(2.30)
where,
IJ J I
IJ J I
WWX WX
WW W
21
or
IJ
IJ IJ
T
WW W (2.31)
The weight functions behave like finite element global shape functions.
Divide by
IJ
i
X from (2.7)
,
IJ
IJ
IJ
i IJ IJ
ii
WW
W
X X
(2.32)
but
IJ
IJ
ii
dW W
dX X
Then from (2.32), (2.5) and (2.6)
,
1
IJ
IJ
IJ i
i IJ IJ
ii i
IJ
i
IJ
i
Y dW W
W
dX X Y
Y
X
thus,
,,
1
IJ
IJ
ii IJ
i
W
WW
Y
(2.33)
Introducing (2.17) and (2.33) in (2.26), the following weighted residual equation results:
1
,0 ,
11
0
0
11
11
.
1
o
o
IJ
IJ
NN
IJ n oIJ IJ IJ
ik nk i mn IJ IJ
IJ
k i
IJ
NN
IJ
m
IJ
IJ
W
u
W
udv YY
C
v
NY
Y
f
Wdv
v
N
(2.34)
22
Now consider the 1 O terms in the equation of equilibrium (2.34). These are
associated with the evolution of the macro-solutions.
,0 ,
11
1
0 ,
11
0
0
11
1
1
1
o
o
o
NN
IJ IJ oIJ IJ IJ
ik nk i mn
IJ
IJ
IJ
NN
IJ IJ n IJ IJ
ik i mn
IJ
IJ
k
IJ
NN
IJ
m
IJ
IJ
W
udv YY
C
v
N
u
W
dv YY
C
v
N
Y
W f
dv
v
N
(2.35)
Using the method of asymptotic expansions, various terms can be evaluated as follows,
The
2
1
term has no such terms because it was assumed that
o
ufX
only.
The
1
terms from (2.34) gives the atomistic equation defining the micro-solution
1IJ
n u
:
,
11
1
11
1
1
IJ
NN
IJ
oIJ IJ IJ
ik nk
mn IJ
IJ
i
IJ
IJ
IJ
NN
IJ n IJ IJ
ik
mn IJ IJ
IJ
ki
IJ
W
u YY
C
NY
W
u
YY
C
NYY
(2.36)
or
1
,
11 11
11
NN NN
IJ IJ
IJ IJ
IJ oIJ IJ
k nnk
mn mn
IJ I J
IJ I J
WW
uu Y
CC
NN
(2.37)
23
Using (2.25), the basic weighted residual equation for the atomistic problem is
1
,
11 1 1
11
NN NN
IJ IJ
IJ
IJ
IJ oIJ
nnk
mn
mnk
IJ I J
IJ I J
WW
uu
C
C
NN
(2.38)
Equation (2.38) gives the atomistic correction to the macro deformation gradient. If there
is no heterogeneity in the micro structure such as defects, there is no micro solution.
2.2 Anharmonic Expansion about Equilibrium
To make the solution easily computable by direct methods, the material constants
are approximated through a series expansion about the equilibrium position
() IJ E
n
r . Using
a Taylor series expansion,
2() 3 ()
() 1
,
4()
() 1 ( ) 1
,,
1
2
IJ E IJ E
bn bn
IJ oIJE IJ
lp p l mn IJ IJ IJ IJ IJ
mn m n l
IJ E
bn
oIJE IJ o IJE IJ
lp p l ic c i IJ IJ IJ IJ
mn l i
Er Er
ur u
C
rr rrr
Er
ur u u r u
rrrr
(2.39)
() 1 ( ) 1 () 1
,, ,
1
2
IJ IJ IJ
oIJE IJ oIJE IJ o IJE IJ
mn mnli mnl
lp p l l p p l ic c i
D G ur u H ur u u r u (2.40)
While using the Taylor series is an approximation, it is valid for the range of strains seen
in practice in solid mechanics problems.
24
Introducing (2.40) in (2.38), the following equation results:
() 1
,
1
() 1 ( ) 1
11
,,
,
11
1
1
2
1
IJ IJ
oIJE IJ
mn mnl
NN lp p l
IJ
IJ
n
IJ
o IJ E IJ o IJ E IJ
IJ
mnli
lp p l ic c i
IJ
NN
IJ
IJ
o
mnk
nk
IJ
IJ
DG ur u
uW
N
Hur u ur u
Cu W
N
(2.41)
or,
1
11
() 1 1
,
11
() 1 ( ) 1 1
,,
11
,
11
1
1
11
2
1
NN
IJ IJ
IJ
mn
n
IJ
IJ
NN
IJ IJ
oIJE IJ IJ
mnl
lp p l n
IJ
IJ
NN
IJ IJ
oIJE IJ o IJE IJ IJ
mnli
lp p l ic c i n
IJ
IJ
NN
IJ
IJ
o
mnk
nk
IJ
IJ
Du W
N
Gur u u W
N
Hur u ur u u W
N
Cu W
N
(2.42)
Now this equation must be linearized by setting,
1
11
IJ
IJ IJ
n
nn
uu u (2.43)
25
and introducing this expression in (2.42) to obtain:
1
1
11
11
() 1 1
,
11
1
() 1
,
11
() 1 1
,
1
1
1
1
2
NN
IJ IJ IJ
IJ
mn n
n
IJ
IJ
NN
IJ IJ IJ IJ
oIJE IJ IJ
mnl l n
lp p l n
IJ
IJ
IJ IJ
oIJE IJ
mnli l
lp p l
IJ
IJ IJ
oIJE IJ IJ J
in
ic c i n
Du u W
N
Gur u u u u W
N
Hur u u
W
N
ur uu uu
11
,
11
1
NN
I
IJ
NN
IJ
IJ
o
mnk
nk
IJ
IJ
Cu W
N
(2.44)
Expanding (2.44) and neglecting higher order,
11
()
,
11
()
,
1
11
()
,
11
() ()
,,
1
1
2
IJ IJ IJ IJ
oIJE
mn mnl l n
ns l p p ns ls
IJ IJ
oIJE
ln
si l p p
IJ
s
IJ IJ IJ
oIJE
mnli in
sl i c c
IJ IJ
o IJE o IJE
li
sn l p p i c c
DG ur u u
ur u u
u
N
Hur uu
ur u u r u
11
1
()
,
1
11
() ()
11
,,
,
11
1
1
2
1
NN
IJ
IJ
IJ
IJ IJ IJ
oIJE
mn mnl l
lp p NN
IJ IJ
n
IJ IJ IJ
o IJE o IJE
IJ
mnli li
IJ lp p ic c
NN
IJ
IJ
o
mnk
nk
IJ
IJ
W
DG ur u
uW
N
Hur u ur u
Cu W
N
(2.45)
26
This is the nonlinear atomistic weighted residual model which will be used in
implementation of the method. Using equation (2.23) in (2.45) we obtain:
11
()
,
11
()
,
1
11
()
,
11
() ()
,,
1
1
2
IJ IJ IJ IJ
oIJE
mn mnl l n
ns l p p ns ls
IJ IJ
oIJE
ln
si l p p
IJ
s
IJ IJ IJ
oIJE
mnli in
sl i c c
IJ IJ
o IJE o IJE
li
sn l p p i c c
DG ur u u
ur u u
u
N
Hur uu
ur u u r u
11
1
()
,
1
11
11
() ()
,,
()
,
11
1
1
2
1
NN
IJ
IJ
IJ
IJ IJ IJ
oIJE
mn mnl l
lp p
NN
IJ IJ
n
IJ IJ IJ
IJ
o IJE o IJE
IJ mnli li
lp p ic c
NN
IJ
IJ IJ E o
knk mn
IJ
IJ
W
DG ur u
uW
N
Hur u ur u
ru W
C
N
(2.46)
Introducing (2.40) in (2.46), the following equation results:
11
()
,
11
()
,
1
11
()
,
11
() ()
,,
1
1
2
IJ IJ IJ IJ
oIJE
mn mnl l n
ns l p p ns ls
IJ IJ
oIJE
ln
si l p p
IJ
s
IJ IJ IJ
oIJE
mnli in
sl i c c
IJ IJ
o IJE o IJE
li
sn l p p i c c
DG ur u u
ur u u
u
N
Hur uu
ur u u r u
11
1
()
,
1
11
11
() ()
,,
() 1
,
1
1
2
1
1
2
NN
IJ
IJ
IJ
IJ IJ IJ
oIJE
mn mnl l
lp p
NN
IJ IJ
n
IJ IJ IJ
IJ
o IJE o IJE
IJ mnli li
lp p ic c
IJ IJ
oIJE IJ
mn mnl
lp p l
mnli
W
DG ur u
uW
N
Hur u ur u
DG ur u
N
H
()
,
() 1 ( ) 1 11
,,
NN
IJ
IJ E o
knk
IJ
o IJ E IJ o IJ E IJ IJ
IJ lp p l ic c i
ru W
ur u u r u
(2.47)
27
and linearize the equation (2.47) results as follows,
1
() 1
,
1
() 1
,
1
11
()
,
1
() 1 ( )
,,
1
1
2
IJ IJ IJ
oIJE IJ
mn mnl n
ns l p p l ns ls
IJ
oIJE IJ
n
si l p p l
IJ
s
IJ IJ IJ
oIJE
mnli in
sl i c c
IJ
o IJE IJ o IJE
i
sn l p p l i c c
DG ur u u
ur u u
u
N
Hur uu
ur u u r u
11
1
()
,
1
11
11
() ()
,,
1
() 1
,
1
1
2
1
NN
IJ
IJ
IJ
IJ IJ IJ
oIJE
mn mnl l
lp p
NN
IJ IJ
n
IJ IJ IJ
IJ
o IJE o IJE
IJ mnli li
lp p ic c
IJ IJ IJ
oIJE IJ
mn mnl l
lp p l
W
DG ur u
uW
N
Hur u ur u
DG ur u u
N
1
() 1 ( )
,,
11
1
() 1
,
1
2
NN
IJ IJ IJ
o IJ E IJ IJ E o
mnli l
lp p l k nk
IJ
IJ
IJ
oIJE IJ
i
ic c i
Hur u u r u W
ur u u
(2.48)
28
Now neglecting the higher order terms in equation (2.48) to give,
11
()
,
11
()
,
11
()
,
11
() ()
,,
()
,
1
2
1
1
2
IJ IJ IJ IJ
oIJE
mn mnl l n
ns l p p ns ls
IJ IJ
oIJE
ln
si l p p
IJ IJ IJ
oIJE
mnli in
sl i c c
IJ IJ
o IJE o IJE
li
sn l p p i c c
IJ
IJ E o
mnl
ls k n k
IJ
mnli
DG ur u u
ur u u
Hur uu
ur u u r u
N
Gr u
H
1
11
1
()
,
()
,
1
()
,
1
()
,
1
1
2
NN
IJ
IJ
s
IJ
IJ
IJ
oIJE
i
sl i c c
IJ E o
knk
IJ
oIJE
l
si l p p
IJ IJ IJ
oIJE
mn mnl l
lp p
uW
ur u
ru
ur u
DG ur u
N
H
1
11
11
() ()
,,
1
()
,
()
,
11
() ()
,,
1
1
2
NN
IJ IJ
n
IJ IJ IJ
IJ
o IJE o IJE
IJ mnli li
lp p ic c
IJ IJ IJ
oIJE
mn mnl l
lp p
IJ E o
knk
IJ IJ IJ
o IJE o IJE
mnli li
lp p ic c
uW
ur u u r u
DG ur u
ru
N
Hur u ur u
11
NN
IJ
IJ
IJ
W
(2.49)
29
By rearranging equation (2.49),
1
()
,
1
()
,
11
() ()
,,
11
() ( )
,,
11
() ()
,,
1
1
2
IJ
oIJE
l
lp p ns
IJ IJ
mn mnl
ns
IJ
oIJE
n
ls n k k
IJ IJ
o IJE o IJE
ln
si l p p n k k
IJ IJ IJ
o IJE o IJE
mnli in
sl i c c n k k
IJ IJ
o IJE o IJE
li
sn l p p i c c
ur u
DG
ur u
ur u u r u
N
Hur uur u
ur u u r u
1
11
1
()
,
1
() ()
,,
1
()
,
11
2
NN
IJ
IJ
s
IJ
IJ
IJ IJ IJ
oIJE
mn mnl l
lp p
IJ IJ
o IJE o IJE
mnli l
lp p nk k
IJ
oIJE
i
ic c
uW
DG ur u
Hur u ur
N
ur u
1
11
NN
IJ IJ
n
IJ
IJ
uW
(2.50)
By switching indices and using symmetry in the expansion in equation (2.49) we obtain
(2.51)
1
()
,
1
11
11
() ( )
,,
1
()
,
1
()
,
()
,
2
1
3
2
11
2
IJ IJ IJ
oIJE
mn mnl l
lp p
NN
IJ
IJ
n
IJ IJ IJ
IJ
o IJE o IJE
IJ mnli li
lp p ic c
IJ IJ IJ
oIJE
mn mnl l
lp p
IJ IJ
oIJE
mnli l
lp p
oIJE
ic c
DG ur u
uW
N
Hur u ur u
DG ur u
Hur u
N
ur u
1
()
,
11
1
NN
IJ IJ
oIJE
n
nk k
IJ
IJ
IJ
i
ur u W
(2.51)
30
CHAPTER 3
MOLECULAR MECHANICS: STATIC MULTISCALE
ALGORITHM
3.1 Inter-Atomic Potential and Boundary Condition
Throughout the dissertation, the Lennard-Jones (L-J) potential is used to calculate
the inter-atomic interactions. The L-J potential is a pair potential and in reduced units can
be described as follows:
12 6
11
4
NN
b IJ IJ
IJ
IJ
E
yy
(3.1)
The parameter
IJ
y is the distance between atoms I and J . The chosen physical constants
are for Argon. The utilized constants are = 1.66 x 10
-14
erg and = 3.4 x 10
-8
cm. The
mass of an argon atom is taken to be 6.6 x 10
-23
gram. It should be noted that the problem
units are expressed in reduced L-J units (1 unit Length = 3.4 x 10
-8
cm, 1 unit time = 2.2
x 10
-12
Sec and 1 unit Temperature = 120 Kelvin) unless explicitly stated otherwise.
The potential
b
E for a pair of atoms, a distance r apart, is pictured in Figure 3.1
31
Figure 3.1: Lennard-Jones inter-atomic potential
The equilibrium distance
o
r between atoms can be calculated by setting the first
derivative of the potential to zero as follows:
01.12
o
b
o
r
dE
r
dr
(3.2)
The potential is truncated at a cut off radius
c
r of 2.5. In order to avoid sudden
discontinuity in the potential, the potential is shifted using the potential
'
b
E .
'
'
0
c
b
bb bc c c
r
bc
dE
EE r E r r r r r
dr
Er
(3.3)
32
In order to avoid free surfaces at edges of the lattice, the minimum image
convention is applied. This helps simulate the bulk properties of the lattice as effectively
the image is periodically repeated in both the directions, as shown in Figure 3.2:
Figure 3.2: Molecular dynamics periodic boundary condition
For a box dimension of
x
L and
y
L , the minimum image convection can be applied as
follows:
,,
22 2 2
,,
22 2 2
xx x x
ij ij ij ij
yy y y
ij ij ij ij
LL L L
x x SignR x SignR x
LL L L
y y SignR y SignR y
(3.4)
where,
,0
22
,0
22
xx
ii
xx
ii
LL
SignR x x
LL
SignR x x
(3.5)
33
3.2 Algorithm for Static Implementation
It is necessary to have an algorithm for linking the continuum scale and atomistic
scale solution variables in the multiscale analysis. Considering the two scales to be the
Macroscale and the Microscale, the algorithm, for the molecular statics analysis would
involve multiple steps. First model is set up and then the loads are applied in the macro
problem. The material behavior at each finite element Gauss point in the macro problem
is modeled using unit atomic cell from which the micro solutions
1
u is obtained. The
solved micro solution is used to determine the elastic material constants, which in turn
are employed to solve for the macro displacements. Then the load on the macro problem
is incremented and the process is repeated till the desired load on the macro problem is
reached. This process is shown in Table 3.1.
34
Table 3.1: Multiscale static algorithm
Macro loop:
Step (I): Apply load in the macro problem.
Micro loop:
Step (i): Compute the deformation gradient F .
Step (ii): Compute
() ij
r distance between atoms.
Step (iii): Apply the boundary condition
Step (iv): Apply Cauchy Born correction to
() 1
,
IJ o IJ E IJ
nnkk n
rur u .
Step (v): Compute the second order tensor
2
IJ b
mn
IJ IJ
mn
E
C
rr
Step (vi): Compute the third order tensor
2
b
IJ mnk
mnk
E
C
rF
Step (vii): Solve for
1
u .
Step (viii): Check for convergence residue tolerance or keep iterating from step(i).
Step (ix): After convergence, compute the elastic constants
2
b
mnik
mi nk
E
C
FF
Step (II): Assemble the macro stiffness and force matrix.
Step (III): Solve for macro displacement
o
u .
Step (IV): Check for convergence
00
old
uu or keep iterating from Step (i).
Step (IV): After convergence return to Step (I) for load increment.
35
3.3 Computer Implementation
An interesting aspect of the static analysis method, presented in equation (2.38), is
the definition of the effective atomistic stiffness matrix or Hessian. The assembly of the
stiffness matrix is based on the definition of the atomistic weight functions in the
weighted residual model. These weight functions are defined in equations 2.27-2.32 and
result in an assembly algorithm which is motivated by standard finite element assembly
procedures.
In this assembly procedure, the interactions between the each pair of atoms, is treated as a
one dimensional element with two nodes. Each node has D degrees of freedom depending
on the dimension of the problem. (1-D, D = 1; 2-D, D = 2; 3-D, D = 3). The assembly
program in the C++ language is presented in Table 3.2. In Table 3.2, the following
variables are used:
iglob = lotogo[pair][iloc] – The connectivity matrix lotogo maps the pairs local
node number to the global node number.
ieq = (iglob * ND) + idir – sets up the global equation number
signdis – is the global shape function for displacement
signwt – is the global weight function associated with the weighted residual
method.
36
Table 3.2: Assembly algorithm for
IJ
mn C
Assembly algorithm for
2
IJ b
mn
IJ IJ
mn
E
C
rr
Loop over ij pairs
for ( int iloc = 0; iloc < 2; iloc++)
for( int jloc = 0; jloc < 2; jloc++){
iglob = lotogo[pair][iloc];
jglob = lotogo[pair][jloc];
for( int idir = 0; idir < ND; idir++)
for( int jdir = 0; jdir < ND; jdir++){
ieq = (iglob * ND) + idir;
jeq = (jglob * ND) + jdir;
if(jloc == 0) signdis = -1;
if(jloc == 1) signdis = 1;
if(iloc == 0) signwt = -1;
if(iloc == 1) signwt = 1;
i = lotogo[pair][0];
j = lotogo[pair][1];
Compute second order terms;
Compute third order terms;
Compute fourth order terms;
Matrix = (a * second + 2.0 * b * third + (3.0/2.0) * c * fourth);
Stiffness[ieq][jeq] = Stiffness[ieq][jeq] + ( Matrix * signdis * signwt / nAtom);
}
}
37
3.4 Results for One Dimension Lattice Cell
The proposed homogenization theory can be demonstrated for a one
dimensional problem with a defect. The atoms within the vertical lines represent the unit
cell of the problem. The atoms on the left and right of the unit cell are imaginary atoms to
take into account the minimum image convention. The center atom of the unit cell is
displaced from its initial configuration in order to fictitiously introduce a defect. Then the
cell is subjected to 10% strain with a constant defect ratio having different numbers of
atoms in the unit cell. A 1-D lattice with a defect in a unit cell of length L is pictured in
Figure 3.3.
Figure 3.3: 1-D Point defect
The behavior of the micro solution
1
u for a constant defection ratio of
/0.01 Lr with varying number of atoms in the unit cell is presented in Figure 3.4.
There is a spike in the
1
u solution at the location of the defect showing discontinuity in
the displacement field. You will notice that as the number of atoms in the cell is
increased the spike in the solution gets sharper, meaning the defect presence is not felt by
the atoms located far away from the defect.
Defect
38
-0.00010
-0.00005
0.00000
0.00005
0.00010
0.00015
0.00020
0.00025
0 0.2 0.4 0.6 0.8 1
X / Unit Cell Lenght
Micro Solution
5 atoms
11 atoms
21 atoms
31 atoms
51 atoms
101 atoms
Figure 3.4: Micro solution for fixed point defect ratio (L/r – 0.01)
In the above one dimensional case, the heterogeneity comes from the presence of
defect in the unit cell. In a case of a perfect lattice, with no heterogeneity, there is no
micro solution and the unit cell under deformation due to macro load will undergo a rigid
body displacement.
3.5 Results for a Two Dimension Triangular Lattice Cell Problem
In the previous section we have demonstrated the behavior of the micro solution
in a one dimensional case. In this section a two dimensional version of the results is
presented. A 2-D triangular lattice of 50 atoms is pictured in Figure 3.5, and the lattice
constant that gives the minimum energy configuration is computed to be 1.116.
39
-5
-4
-3
-2
-1
0
1
2
3
4
5
-3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5
X - Coordinate
Y - Coordinate
Figure 3.5: Equilibrium configuration of a perfect lattice – Energy/per atom = -2.9717
A point defect is introduced into the perfect lattice in order to have some kind of
heterogeneity in the micro structure. The equilibrium configuration of the triangular
lattice with a point defect is pictured in Figure 3.6. The minimum energy equilibrium
position was computed by the conjugate gradient method. Conjugate gradient is a robust
technique to determine the minimum energy of the system and is frequently used in
molecular mechanics for this type of calculations.
-5
-4
-3
-2
-1
0
1
2
3
4
5
-3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5
X - Coordinate
Y - Coordinate
Figure 3.6: Equilibrium configuration of a point defect – Energy/per atom = -2.8523
40
Given the equilibrium position of the lattice, the positions under given levels of
macroscopic strain can be computed using the “micro loop” included in the algorithm for
static analysis in section 3.2. The lattice with a point defect is the atomic unit cell under
consideration. Figure 3.7 shows the undeformed lattice and the deformed lattice under
2.5% strain in the X-direction as computed using homogenization and the lattice under no
strain. The presentation gives the illustration of the movement of the atom subjected to
strains. As can be seen as the strains are applied to the system, the atoms move apart.
In figure 3.8 the atom displacement results from homogenization and those
computed by the conjugate gradient method are compared. The results are almost
identical. The displaced co-ordinate from both approaches falls on top of the other and
are almost indistinguishable. The strain level used in this test was 2.5%.
-5
-4
-3
-2
-1
0
1
2
3
4
5
-3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5
X - Coordinate
Y - Coordinate
Underformed (No Strain)
Deformed (Under Strain)
Figure 3.7: Lattice under 2.5% strain in the Y - direction
41
-5
-4
-3
-2
-1
0
1
2
3
4
5
-3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5
X - Coordinate
Y - Coordinate
Conjuage Gradient
Homogenization
Figure 3.8: Lattice under 2.5% strain in the X - direction
The overall accuracy of the microscopic homogenization calculation can be
established by computing the stored energy in the lattice. Figure 3.9 compares the energy
from homogenization and the conjugate gradient method, for a lattice under strain in the
Y-direction.
-2.86
-2.85
-2.84
-2.83
-2.82
-2.81
-2.8
-2.79
-2.78
-2.77
-2.76
01 2 3 45 6
Percent Strain
Energy
Homogenization
Conjugate Gradient
Figure 3.9: Energy/atom: Lattice under strain in Y - direction
42
Figure 3.10 compares the energy from homogenization and the conjugate
gradient method, for the lattice under stain in the X-direction.
-2.86
-2.85
-2.84
-2.83
-2.82
-2.81
-2.8
-2.79
-2.78
-2.77
01 2 3 4 56
Percent Strain
Energy
Homogenization
Conjugate Gradient
Figure 3.10: Energy/atom: Lattice under strain in X - direction
The deformations of the lattice, due to the macroscopic deformations, cause
changes in the elastic constant
ijkl
C . The behavior of the homogenized material elastic
constant
1111
C and
2222
C as a function of the applied strain is studied. After each strain step
the material constants are used to update the macro finite element problem to obtain the
macro variables such as displacement.
Figure 3.11 shows the behavior of the material elastic constants (
1111
C and
2222
C ),
when the lattice is under strain in the Y-direction, computed using homogenization.
43
0
10
20
30
40
50
60
70
80
90
02 46 8
Percent Strain
Elastic Constant
C1111
C2222
Figure 3.11: Constitutive Constants: Lattice under strain in Y - direction
Figure 3.12 shows the behavior of the material elastic constants (
1111
C and
2222
C ),
when the lattice is under strain in the Y-direction, defined using homogenization. These
results demonstrate the workability of the “micro loop” in section 3.2.
0
10
20
30
40
50
60
70
80
90
024 6 8
Precent Strain
Elastic Constant
C1111
C2222
Figure 3.12: Constitutive constants: Lattice under strain in X - direction
44
3.6 Results for Static Multiscale Algorithm Implementation
In the previous section the convergence of the micro solutions for a unit cell is shown. A
2-D structure consisting of Quad elements as shown in Figure 3.13 was chosen as test
problem to demonstrate the working of the multiscale algorithm present in Table 3.1. The
macro structural model is constrained at the left end and is free at the right end. Structural
load are applied at the free end in the positive X-direction. Each gauss point in the macro
structural model consists of a unit cell of with a point defect as described in the previous
section. The results present in section 3.5 are an outcome of the micro loop
implementation of the multiscale algorithm. In this section the micro loop is implemented
for every macro load increment till the desired strain is reached.
Figure 3.13: 2-D Structural model consisting of Quad elements
45
The forces are applied in load increments till a maximum strain of 5 percent is achieved
in the macro-model. In Figure 3.14 the relationship between applied force and the macro
free end displacement can be observed. The macro displacement is obtained after there is
a concurrent convergence between macro field variable
o
u and the micro field variable
1
u .
Since the convergence and implementation of the micro loop has been validated in
section 3.5, the primary step in the implementation of the multiscale algorithm is get a
compatible field variables belonging to both the scale converge simultaneous. To get a
qualitative comparison of the results an equivalent elastic model is solved. The results of
the elastic model are represented along with the homogenization results in Figure 3.14. It
can be observed the difference between the homogenization model and elastic model
starts to deviate as function of increasing strain.
0
20
40
60
80
100
120
01 2 3 45 6
Free end displacement
Applied Force
Non-linear Homogenization Model
Linear Elastic Model
Figure 3.14: Force displacement results for the macro model
In Figure 3.15 the calculation of the elastic constant C1111 at a particular gauss point for
every macro load increment is presented. The results follow the expected behavior as
46
demonstrated in the results presented in section 3.5. Thus the implementation of the static
multi-scale algorithm plays a vital role in realistically analyzing the homogenized
material behavior of macro model, where the inherent constitutive parameters are directly
derived from the atomistic model.
0
10
20
30
40
50
60
70
80
90
0 0.010.020.030.040.050.06
Strain Rate
C1111
Figure 3.15: Gauss point evaluation of C1111 elastic constant
47
CHAPTER 4
NON-EQUILBRIUM MOLECULAR MECHANICS, HEAT FLUX
AND VIRIAL STRESS
4.1 Non-Equilibrium Molecular Dynamics
A majority of the molecular dynamics simulations in the literature are micro-
canonical in nature. Micro-canonical ensembles are ensembles wherein the total energy,
number of atoms and volume of the system is constant. Microcanonical ensembles are
also known as NVE ensembles. In case of Canonical ensembles the total energy of the
system is not held constant. They are knows as NVT ensembles wherein the number of
atoms, the volume of the system and temperature in the system are held constant. In
addition there are Isothermal-Isobaric (NPT) ensembles wherein the number of atoms, the
pressure in the system and temperature in the system are held constant. In this section, the
discussion is restricted to Microcanonical ensembles and its extension to Non-
Equilibrium Molecular Dynamics.
In the case of Microcanonical, ensembles the temperature of the system reaches a steady
state, meaning the temperature of the system is evenly distributed while keeping the total
energy of the system constant. In Non-Equilibrium Molecular Dynamics, there is a
48
temperate profile such as a temperature gradient in the molecular dynamics system. In
case of Non-Equilibrium Molecular Dynamics simulation the average temperature of the
system and total energy of the system is held constant. There are various types of Non-
Equilibrium Molecular Dynamics algorithm available in literature. Following are some of
the popular Non-Equilibrium Molecular Dynamics techniques:
1. Constant temperature gradient and momentum not conserved (CTGMNC)
2. Constant temperature gradient and momentum conserved (CTGMC)
3. Constant heat flux and momentum conserved (CHFMC)
a. Jund method
b. Velocity exchange (VE) method
4. Constant heat flux and momentum not conserved (CHFMNC)
For a detailed description of all the above methods, the reader is referred to the work of
Hung et al. [26]. The paper discusses a good comparison among all the different methods.
In Microcanonical ensembles the total momentum of the system is conserved by simply
setting the initial velocities of the system in such a way that the total momentum is zero
The integration algorithm preserves the momentum throughout all the time steps. In case
of Non-Equilibrium Molecular Dynamics simulation preserving the total momentum of
the system is not trivial. Since in Non-Equilibrium Molecular Dynamics there is a
temperature profile, the momentum conservation can be easily disturbed. Depending
upon the velocity distribution in the Non-Equilibrium Molecular Dynamics system, it can
49
cause the entire system to swift in time. So it becomes important to preserve the total
momentum of the system. In this section, the discussion will be limited to describing the
constant heat flux and momentum conserving Jund method and the constant temperature
gradient and momentum conserved method. The Jund method is used as a sub-step in the
multi-scale thermal analysis procedure to be defined in subsequent chapter of this
dissertation. The current presentation is designed to verify the performance of the Jund
method in conjunction with the heat flux definitions resulting from the multi-scale
homogenization derivation presented in Chapter 5.
4.1.1 Constant Heat Flux and Momentum Conserving Method
This method was proposed by Jund et al. [31]. An Argon 100 x 4 x 4 F.C.C
structure as shown in Figure 4.1 was chosen to demonstrate the algorithm. The F.C.C
structure is sliced at 1/4
th
the length and 3/4
th
the length to serve as heat sink and heat
source plates respectively. Periodic Boundary Conditions (PBC) is applied on the three
directions as described in Chapter 3. By imposing a heat source and heat sink plates at
1/4
th
the length and 3/4
th
the length coupled with Periodic Boundary Condition, a
molecular dynamics loop is created. Such a set up ensures a more realistic system and not
disturb the symmetry of the system.
50
Figure 4.1: Argon 100 x 4 x 4 F.C.C structure
An appropriate width of 2 is chosen for the heat source and the heat sink plates such that
there are at least 100 atoms in the model. This ensures that a stable temperature profile is
established. After setting the initial velocities of the system to a target temperature
t
T , the
system is equilibrated for sufficient runs. The equilibration process is done by rescaling
the current velocities in the system to the target temperature. There are various
equilibration processes in the literature, which are not discussed here. Once the desired
average temperature of the system is reached, a constant energy e is input into the heat
sink and a constant energy of e us removed through the sink after each time iteration,
thus not disturbing the total energy of the system. To ensure the total momentum is
conserved the following velocity scaling of atoms in the heat source and sink is adopted.
The scaling of the velocity is carried out using the following formula:
,, Inew G I old G
vv v v (4.1)
where
, Iold
v - the current velocity of atom i before velocity scaling
51
, Inew
v - the new velocity of atom i after velocity scaling
G
v - the velocity of center of mass of the atoms in the plate
The scaling parameter which ensures the momentum conservation can stated as
follows:
1
R
C
e
E
(4.2)
and
R
C
E as
22
,
11
22
R
CIIold IG
II
E mv mv
(4.3)
where,
I
m - The mass of atom I
Figure 4.2 shows the implementation of the Non-Equilibrium Molecular Dynamics
method proposed by Jund for an average temperature of 0.375. It can be observed from
the figure that the temperature gradient is about the average temperature. Since the heat
sink and heat source plates are located at 1/4
th
the length and 3/4
th
the length of the
system, the symmetry in the set up can be easily observed.
52
0.25
0.3
0.35
0.4
0.45
0.5
0.55
0 1020 30 4050 60
Cell Number
Temperature
Figure 4.2: Non-equilibrium molecular dynamics Jund method
4.1.2 Constant Temperature Gradient and Momentum Conserved
Method
The constant temperature gradient and momentum conserved method is similar to
the Jund method in that a scaling procedure ensuring that a desired temperature gradient
is maintained in the system is used to define target temperatures at the heat sink and heat
source plates. The system considered is similar to the one described in Figure 4.1 expect
that the heat sink and heat source are placed at either end of the system and periodic
boundary conditions are imposed only along the directions perpendicular to the length of
the system. The only difference between the two methods is the scaling parameter
defined as follows:
53
tG
CG
TT
TT
(4.4)
where
t
T - The current temperature of the plate
t
T - The target temperature of the plate
and
G
T is defined as follows:
2
1
3
GiG
i
B
Tmv
NK
(4.5)
where,
N - the number of atoms in the plate
B
K - the Boltzmann Constant
Figure 4.3 shows the implementation of the Non-Equilibrium Molecular Dynamics
constant temperature gradient and momentum conserving method for an average
temperature of 0.4 with a target temperature of 0.3 at the sink and a target temperature of
0.5 at the source. As seen in the previous section, it can be observed from the figure that
the temperature gradient is about the average temperature. Since heat sink and heat
source plates are located at start and end of the system with no periodicity along the
length, the system lacks symmetry.
54
0.25
0.3
0.35
0.4
0.45
0.5
0.55
0 5 10 15 20 25
Distance
Temperature
Figure 4.3: Constant temperature gradient and momentum conserving method
4.2 Heat Flux
Heat flux is defined as the flow of energy per unit of area per unit of time. Heat fluxes
can be computed once a steady state temperature gradient is established. As described in
the pervious section temperature gradients can be established by using any of the above
mentioned Non-Equilibrium Molecular Dynamics methods. The heat fluxes can be
computed using Eulerian or Lagrangian reference framework. The heat fluxes
expressions in both the coordinate system can be found in the works done by Hoover [25]
and Li et al. [34].
55
The heat flux q in the Eulerian reference frame can be stated as follows:
2
,
11
2
IJ I IJ I I I I I
ijj i i i
IJ I I I
c
q F vv y m vv vv vv
V
(4.6)
where,
IJ
j
F - the force on atom I due to atom J in direction j
I
- the potential energy of atom I
I
i
v - the velocity of atom I in direction i
IJ
j
y - the distance between atom I and atom J in direction j
c
V - volume of the cell
v - the average velocity of a control volume
c
V
The heat flux q in the Lagrangian reference frame can be stated as follows:
,
1
IJ I IJ
ijj i
IJ I
c
QFvvY
V
(4.7)
where,
IJ
j
Y - The unreformed distance between atom I and atom J in direction j
56
To demonstrate the implementation of heat flux computation, a model as
described in Figure 4.1 is taken with heat source and heat sink plates at 1/4
th
the length
and 3/4
th
the length of the system respectively. Periodic boundary conditions were
imposed in all three directions. The system was equilibrated about an average
temperature of 0.375 as shown in Figure 4.2. The system was sub-divided into a series of
cells along the length of the system to compute the heat fluxes. The heat flux was
computed for each cell by averaging the fluxes over a sufficient period of time. Figure
4.4 shows the heat flux results computed using both Eulerian and Lagrangian reference
frame. It can also be inferred from the figure that since the gradient of the temperature
imposed is positive the heat flux computed is negative as the flow of energy is in the
opposite direction.
-0.12
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
0 1020 30 4050 60
Cell Number
Heat Flux
Lagrangian Heat Flux
Eulerian Heat Flux
Figure 4.4: Lagrangian and Eulerian heat flux plots using NEMD methods
57
The readers are referred to the works of Hafskjold et al. [23], Ikeshoji et al. [28],
Mountain et al. [27] and Bugel et al. [4] for a detailed discussion of heat fluxes and the
methodology used in the implementation of the procedure.
Since the current formulation is Lagrangian based, the heat flux expression stated
in Eq. (4.7) would be used in chapter 5, where the multi-scale thermo-mechanical model
is derived. The relevance of the current heat flux discussion involves the fact that the
multi-scale thermal analysis procedure, presented in Chapter 5, includes a Lagrangian
heat flux term resulting from the analysis. The heat flux agrees with the expression in Eq.
(4.7) described in the chapter, can be used to compute the multi-scale heat flux.
4.3 Virial Stress
Virial stress is the averaged atomistic stress of a molecular dynamics ensemble and
should not be confused with the atomistic point stress. There are many derivations for the
Virial stress expression since the initial works of Irving et al. [29]. The expression
derived by Irving et al. is difficult to implement in Molecular Dynamics simulation. The
work done by Hardy [24] in deriving the Virial stress expression overcomes this barrier.
The Virial stress
ij
S expression derived by Hardy is stated as follows:
,
1
ij IJ J I IJ I I I
ij ij
IJ I I
SFyyy mvv
Y
(4.8)
58
where,
Y - The volume of the system
Later the expression in equation (4.8) was modified to include only fluctuation part of the
atomistic velocities and was restated as follows:
,
1
ij IJ J I IJ I I I
ij i j
IJ I I
c
SFyyy mvvvv
V
(4.9)
If the volume element
c
V is taken as the entire volume of the system Y then the average
velocity v is set to zero. The virial stress expression in equation (4.9) is also referred to
as Eulerian virial stress in continuum mechanics. The work done by Li et al.[34] also
states the Lagrangian virial stress expression as follows:
,
1
ij IJ J I IJ
ij
IJ I
SFyyY
Y
(4.10)
In the work of Zhou [51], it is reasoned that the virial stress expression should not include
the kinetic terms, and the virial stress is restated as follows:
,
1
ij IJ J I IJ
ij
IJ I
SFyyy
Y
(4.11)
59
Though the expression purported by Zhou for sometime was widely accepted as the
expression for Virial stress, it has recently been refuted by researchers. The work
presented by Liu et al. [35] addresses most of the questions posed by Zhou and states that
there were inconsistencies in his derivation.
In this section, the virial stress is computed numerically from three proposed
expressions, namely – Eulerian Virial stress, Lagrangian Virial stress, Zhou Virial stress
to determine the valid expression to be adopted. An Argon 6 x 6 x 6 F.C.C structure was
chosen to calculate all three virial stress expression for different temperature ranges.
The results are shown in Figure 4.5. The Virial stress expression computed using
the Lagrangian and Eulerian formulation agree very well whereas the expression
purported by Zhou with the kinetic term missing differs from the accepted versions. Thus
using simple numerical test it is found that the Virial stress expression purported by Zhou
is incorrect.
0
1
2
3
4
5
6
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
Average Temperature
Virial Stress
Lagrangian virial stress
Eulerian virial stress
Zhou virial stress
Figure 4.5: Lagrangian, Eulerian and Zhou Virial stress plots
60
The relevance of the virial stress discussion involves the fact that the multi-scale thermal
analysis theory, developed in Chapter 5, includes a Lagrangian stress term, resulting from
the analysis. The stress expression matches the virial stress from equation (4.10). The
above analysis validates the appropriate of this stress expression in the computational
algorithms results from Chapter 5.
61
CHAPTER 5
DYNAMIC ANALYSIS: ATOMISTIC – CONTINUUM
5.1 Dynamics Behavior in Macro and Discrete Fast Time Scales
Suppose that there are two time scales, a slow timet and a fast time , used to
model the time dependent structural and atomistic deformation and that the two scales are
related through the following expressions:
t
(5.1)
The time scale t is associated with the average motion of the macroscopic structure. The
time scale is associated with the motion of the individual atoms. In the asymptotic
expansion from equation (2.9), the temporal variables must be included. The following
expansions for displacement, velocity, and acceleration result:
01
,,
II II
nn n
ut u X t u Y
(5.2)
01
,,
I
II I I ii
ii
uu
vu X t Y
tt
(5.3)
62
20 2 1
2
22 2
,,
III
I
nn
n
uX t u Y
u
tt t
(5.4)
Assuming that
1I
n
u is a function of only from equations (5.1) and (5.4), it can be seen
that
21 21 2
21
2222
,,
1
II II
I
nn
n
uY uY
u d
tdt
(5.5)
Introducing equation (5.5) in (5.4), the acceleration takes the following forms:
20 2 1
2
22 2
,,
1
III
I
nn
n
uX t u Y
u
tt
(5.6)
A multi-scale version of the Lagrangian continuum momentum equation in space,
00
,
ij i
j
VV
v
WdV S W dV
t
(5.7)
where
, j
j
W
W
X
Now consider the term on the left had side in equation (5.7)
0
0
0
1
1
1
,
i
V
i
V
Y
I N
I i
I
V
I
Y
v
WdV
t
v
WdY dV
Yt
v
m Y Y W X Y dY dV
Yt
(5.8)
63
where
,,
I
I
Y
YY W XYdY W XY
(5.9)
A useful decomposition of the two-scale weight is defined as follows:
,
,
I
I
WXY W X W Y
WXY W X W
(5.10)
where,
I
II
YY
WW Y
Introducing equation (5.10) in (5.8), it can be seen that
00
0
1
1
1
,
1
I N
ii
II
VV
I
I N
I
i
I
V
I
vv
WdV M W X Y dV
tY t
v
mWXWdV
Yt
(5.11)
or
00 0
11
11
II NN
I
ii i
II
VV V
II
vv v
WdV m W X dV m W dV
tY t Y t
(5.12)
64
Introducing equation (5.3) in equation (5.12) and averaged over fast time, we obtain
0
0
0
20 2 1
22
1
20 2 1
22
1
11 1
11 1
i
V
N
ii
I
V
I
N
I
ii
I
V
I
v
WdV
t
uu
mWXddV
Yt
uu
mWddV
Yt
(5.13)
By setting the atomistic weight
I
W equal to zero in equation (5.13), the 1 O inertial
term can be shown to take the following form:
00
20
2
1
11
N
ii
I
VV
I
vu
WdV m W X d dV
tY t
(5.14)
Similarly, by setting the macro weight WX
equal zero, the
1
O
inertial term can be
obtained
0
21
0
2
1
1
N
I
ii
I
V
I
vV u
WdV m W d
tY
(5.15)
Now consider the multi-scale stress term averaged over space and time
00
0
,,
,
11
11
11
ij ij
jj
VV
Y
NN
IJij
j
V
IJ
JI
SW dV SW dY d dV
Y
SW d dV
Y
(5.16)
65
On line
IJ
Y
, a parameter 0,1 exists such that
1
IJ
YY Y
(5.17)
Thus,
JI
dY Y Y d
(5.18)
or,
IJ
dY
d
dY
(5.19)
where,
IJ J I
kk k
YY Y (5.20)
In order to obtain a fully discrete temporal algorithm, define the macrostructure
discrete temporal variable such that, for macro time step t ,
L
tt Lt
0,, LO
Let K be the time step number associated with atomistic motion such that for atomistic
time step
k
K
0,, KM
The fully discrete displacement vector, at atom I , macro time step L , and atomistic time
step, K is
, ILK
n
u
66
and
,, ,
,, ,
IJ LK J LK I LK
nk nk nk
uu u
Similar to (2.17), the asymptotic expansion for the displacement gradient is
1,
,,
,,
IJ K
IJ LK o L n
nk nk IJ
k
u
uu
Y
(5.21)
and the asymptotic expansion for the displacement at these time points is from (5.2)
,0, 1,
II
ILK L I K
nn n
uu X u Y (5.22)
In addition the temporally discrete version of the acceleration expansion of (5.6) is
20, 2 1,
2,
22 2
1
II
LIK
ILK
nn
n
uX u Y
u
tt
(5.23)
Thus from (5.16) and (5.17)
00
1
,
,
11 1
0
11
K MN N
ij ijIJ K
j
VV
KI J
j
JI
W
SW dV S Y d dV
MY X
(5.24)
where
,, ,
1
ijIJ K IJ K J I IJ K
ij
SFyyY
Y
(5.25)
67
Which is the same at all 0,1 . The introducing (5.25) in (5.24)
00
1
,,
,
11 1
0
11
K MN N
ij IJ K J I IJ K
ji j
VV
KI J
j
JI
W
SW dV F y y Y d dV
MY X
(5.26)
The weight functions associated with atom I and fast time step K can be decomposed as
follows:
,
K
K
WX WX W
(5.27)
where
,,
1
KIK JK
WW W (5.28)
and
,,
K
JKIK
W
WW
(5.29)
The spatial derivative of the weight function can be simplified as follows:
,
K
K
j
WW W Y
W
X XY X
or
,
1
K
K
j
WW
W
X Y
(5.30)
but
KK
WW
YY
(5.31)
68
From (5.17), the following expressions result:
I
J I
YY
YY
(5.32)
and
11
J IIJ
d
dY Y Y Y
(5.33)
Thus from (5.30), (5.31), (5.29), (5.33)
,,
,
11
JK I K
K
j IJ
jj
W
WWW
X Y
(5.34)
Inserting (5.34) in (5.26), the following expressions is obtained:
0
0
,
1
,,
,,
11 1
0
11 1 1
ij
j
V
MN N
JK I K
IJ K J I IJ K
ij IJ
V
KI J
jj
JI
SW dV
W
Fy yY W W d dV
MY X Y
(5.35)
And integrating on , it can be seen that
0
0
,
,,
,,
11 1
,
1
1
1
ij
j
V
IJ K J I IJ K
ij
MN N
j
V
JK I K
KI J
IJ K J I
JI
i
SW dV
W
Fy yY
YX
dV
M
Fy y W W
(5.36)
69
The 1 O term in (5.36) takes the following form:
0
0
,
,,
11 1
11
ij
j
V
MN N
IJ K J I IJ K
ij
V
KI J
j
JI
SW dV
W
Fy yY dV
MY X
(5.37)
The
1
O
term in (5.36) is defined as follows:
0
,,
, 0
,
11 1
1
MN N
JK I K
ij IJ K J I
ji
V
KI J
JI
V
SW dV F y y W W
MY
(5.38)
The macro equation of motion is obtained by equating the 1 O terms from (5.14) to
(5.37) to give
0
0
20
2
1
,,
11 1
1
11
N
i
I
V
I
MN N
IJ K J I IJ K
ij
V
KI J
j
JI
u
MWXdV
Yt
W
Fy yY dV
MY X
(5.39)
or
00
20
,,
2
11 1
11
MN N
IJ K J I IJ K i
ij
VV
KI J
j
JI
u W
W X dV F y y Y dV
tM Y X
(5.40)
Since the weight functions
, J K
W and
, IK
W have an arbitrary dependence on K , it can be
assumed that the weights are such that
,
,
IK I
JKJ
WW
WW
forKK
70
and that,
,
,
0
0
IK
JK
W
W
for
1, , KM
KK
The micro equation of motion is obtained by equating the
1
O
terms in (5.15) to (5.38)
to give
21,
,,,
,
2
111
IK NNN
IK J K IK
IJ K J I i
Ii
IIJ
JI
u
MW F yyW W
(5.41)
Consider the typical atom L,1 LN . Then define the fine scale weights in Eq. (5.41)
using the Kroneckor Delta, as follows:
,
,
,
IK
IL
JK
JL
IJ K
JLIL
W
W
W
(5.42)
Introducing (5.42) in (5.41) , it can be seen that
21,
,
2
111
0
IK NNN
IJ K J I m
IIL i JLIL
IIJ
JI
u
MFyy
(5.43)
Or simplifying (5.43), the following atom L equation of motion is obtained:
21 ,
,,
2
11
0
LKNN
IL K J I LJ K J I m
Li i
IJ
IL J L
u
MFyy Fyy
(5.44)
1,
1, ,
LN
KM
71
5.2 Thermal behavior in Macro and Atomistic Scales
From statistical mechanics, the kinetic energy is related to the temperature T
through the following expression:
2
2
2
1
11
22
I N
I n B
I
u DK T
m
Nt
(5.45)
where D is the dimension of the system.
Introducing (5.3) into the right hand side of (5.45), it can be seen that
2
01,
1
22
1, 1 ,
1
11
2
1
2
2
IK N
I nn
n
I
II IK IK N
oo n n
I
I
uu
m
Nt
uu u u
m
Nt t
(5.46)
The velocity expression
, IK
v of atom I at particular micro time step K is defined by
,1,
,
IK o I K
IK nn n
uu u
v
tt
(5.47)
The averaged velocity
K
v over N atoms of equation (5.47) can be stated as follows:
1,
1
1,
11
1,
1
1
11
1
IK N
K
o
I
IK NN
o
II
IK N
o
I
u u
v
Nt
u u
Nt N
u u
tN
(5.48)
72
The averaged velocity v over N atoms and M micro times of equation (5.48) can be stated
as follows
1,
11 1 1
1,
11
11 1
11
IK MM M N
K
o
KK K I
IK MN
o
KI
u u
vv
Mt M N
u u
tM N
(5.49)
where
1
1,
11
11
IK MN
KI
uu
MN
Equation (5.49) can thus be restated as
1
o
u u
v
t
(5.50)
Using (5.50), equation (5.47) can be seen to take the following form:
1
1,
,
IK
IK
uu
vv
(5.51)
Squaring (5.51), the following form results:
2
11
1, 1 ,
2 2
,
2
IK IK
IK
uu uu
vv v
(5.52)
73
Averaging (5.52) over micro time and space produces the following expression:
2
,
11
1
1,
2
11 11
2
1
1,
11
1
1, 1
2
11
11
11 11
2
11
11 11
2
MN
IK
KI
IK MN MN
KI KI
IK MN
KI
IK I MN
KI
v
MN
uu
vv
MN MN
uu
MN
uu u
vv
MN MN
2
1
,
11
K MN
KI
u
(5.53)
The center term in (5.53) can be seen to be zero
1
1,
11
1
1,
11 1 1
11
11
2
11 11
22
22 0
IK MN
KI
IK MN MN
KI KI
uu
v
MN
uu
vv
MN MN
uu
vv
(5.54)
Now using (5.54), equation (5.53) becomes
2
1
1,
2 2
,
11 1 1
11 11
IK MN M N
IK
kI K I
uu
vv
MN MN
(5.55)
74
Introducing (5.55) in (5.45), the temperature can be decomposed into mean and
fluctuating parts:
2
,
11
2
1
1,
2
11 1
21 1
2
21 1 1
22
2
MN
IK
I
KI
B
IK MN N
II
KI I
B
mean fluctuation
B
Tmv
DK M N
uu
mv m
DK M N N
DK
(5.56)
Equation (5.56) can be restated using C
2
1
1,
2
11 1
11 1 1
,
22
IK MN N
II
KI I
uu
Tt m v m
CM N N
(5.57)
where
2
B
KD
C
Using (5.57), space where the temperature field , Tt can be decomposed into macro
and atomistic parts
0
Tt and
1
, Tt as follows:
01
,, Tt T t T t
2
0
11
11
2
MN
I
KI
CT t m v
MN
(5.58)
75
2
1
1,
1
11
11
,
2
IK MN
I
KI
uu
CT t m
MN
(5.59)
The total atomistic temperature
1,K
T , the average of the individual atom temperatures
1, IK
T , appearing in the atomistic term in (5.59), is defined as follows
1, 1 ,
1
1
N
K IK
I
TT
N
(5.60)
Where
2
1
1,
1,
2
IK
IK I
mu u
T
(5.61)
Introducing (5.61) in (5.59), the following expression results:
1,
1
11
11
MN
IK
KI
CT t T
MN
(5.62)
Thus temperature field T can be defined as follows:
01
1,
0
11
,
,
1
MN
IK
KI
CT t
CT t CT t
C
CT t T
MN
(5.63)
76
Taking the derivative of T in (5.63), the following result is obtained:
01 1
01 1
1
,
,,
1
, ,
1
,,
1
dT t
C
dt
T t Tt Tt
CC C
tt
Tt T t Tt
CC
t
Tt T t
CC
t
(5.64)
Equation (5.64) can be restated using (5.62) to give
,
0 1
11
11
IK MN
KI
T T dT C
CC
dt t M N
(5.65)
5.3 The Finite Temperature Problem
The variation of the temperature in the model can be assessed by considering T
to be a macroscopic quantity as defined from statistical mechanics. To obtain the
equation governing the macroscopic temperature field, the first law of thermodynamics
can be utilized:
ddE
Q
dt dt
(5.66)
77
Where the kinetic energy , the internal energy E , the power of the applied forces
and the heat Q are defined as follows:
1
2
o
o
o
o
ii
v
V
ij
ij
V
i i
V
KvvdV
EdV
SnvdA
Q q dVdA
(5.67)
In these energy expressions, is the mass density, is the internal energy density,
ij
S is
the second order stress tensor,
i
n is the normal to the boundary and q is the applied heat
flux. Transformations between surface and volume integrals can be accomplished using
Green’s Theorem:
oo
VV
Fdv F ndA
Using the Green’s Theorem, it can be seen that the heat flux through the surface can be
associated with a volume integral as follows:
oo
i
ii
i VV
q
qndA dV
x
(5.68)
78
a similar mapping can be used for the boundary stress power
oo
ij j
ij
ji
i VV
Sv
SvndA dV
X
(5.69)
Using (5.67), equations (5.69) and (5.68) can be restated as follows:
o
ij j
i V
Sv
dV
X
(5.70)
o
i
i V
q
QdV
X
(5.71)
Then using (5.70) and (5.71) in conjunction with (5.66) and (5.67) the first law of
thermodynamics can be seen to take the following form:
1
2
oo
oo
ii
VV
ij j
i
ii VV
dd
v v dV dV
dt dt
Sv
q
dV dV
XX
(5.72)
or
,,
oo o
oo o
jj
VV
i
ij j i ij i j
i VV V
d
v v dV dV
dt
q
S v dV S v dv dV
X
(5.73)
Regrouping the terms in (5.73), the following expression is obtained:
,
,
oo
oo
jiji j
VV
i
ij j i
i VV
vS vdV dV
q
Sv dV dV
X
(5.74)
79
The first term on the right hand side of (5.74) is zero due to mechanical equilibrium
conditions leaving the energy equation to be satisfied:
,
i
ij j i
i
q
Sv
X
(5.75)
Now introducing (5.75) in (5.66), it can be seen that
,
1
2
oo
oo
i
ii ij ji
i VV
ij
ji i i
VV
q d
v v dV S v dV
dt X
S v ndA q ndA
(5.76)
Inserting body force f term in (5.76), the following expression results:
,
1
2
oo o
oo o
i
ii ij ji
i VV V
ij i
ji jj i
VV V
q d
v v dV S v dV dV
dt X
Sv ndA qndA f vdV
(5.77)
Equation (5.77) can be reorganized to give
,
1
2
oo
oo o o
ii ij ji
VV
i
ij
ji i j j i
iVV V V
d
vvdV S v dV
dt
q
dV S v n dA q n dA f v dV
X
(5.78)
Using (5.69), the second term appearing on the left had side of (5.78) can be modified to
give
oo
ij
ji ij j
i VV
SvndA S v dV
X
(5.79)
80
,
1
2
oo
oo o o
ii ij ji
VV
i
ij j i j j i
ii VV V V
d
v v dV S v dV
dt
q
dV S v dV q n dA f v dV
XX
(5.80)
where the heat flux
i
q is defined as follows in terms of the total atomistic velocity j v
j
iij j
qSv v
Introducing this heat flux expression in (5.80), the resulting balance equation is
,
1
2
oo
oo oo
jj ij ji
VV
j
ij j
ij j i j j i
ii VV VV
d
v v dV S v dV
dt
Sv v
dV S v dV q n dA f v dV
XX
(5.81)
Rearranging the terms in (5.81) produces
,
1
2
oo
oo o
jj ij ji
VV
j
ij
ijj i
iVV V
d
v v dV S v dV
dt
Sv
dV q n dA f v dV
X
(5.82)
Using the symmetry in the stress tensor, indices on the second term on the right hand side
of (5.82) can be switched and the summed index can be replaced with m to give
,
1
2
oo
oo o
mm mi mi
VV
i
mi
mmm m
m VV V
d
v v dV S v dV
dt
Sv
dV q n dA f v dV
X
(5.83)
81
Using the expression forT from (5.45) in (5.83), it can be seen that
,
oo
oo o
mi m i
VV
i
mi
mmm m
m VV V
d
C TdV S v dV
dt
Sv
dV q n dA f v dV
X
(5.84)
Defining a global spatial partition of unity as follows,
1
1
1
R
L
L
PWX
(5.85)
Where L WX
are global weight functions associated with macro nodes L . If L WX
are
associated with the global shape functions
L
X
, then
L
L
WX X
(5.86)
and
1
1
1
R
L
L
PX
(5.87)
The functions in the spatial partition of unity
1
P act to enforce macro energy conservation
properties at the local level in the neighborhood of a global node L .
Similarly, define a global atomistic partition of unity as follows:
2
1
1
IL
N
L
PW
(5.88)
82
where
IL
W can be characterized as an atomistic weight function associated with micro-
atom L and an arbitrary atom I . In this formulation, the arbitrary atom I in the lattice of
the atomistic model is the analog of the arbitrary spatial position X
in the continuum
domain V of the macroscopic model. The natural choice for weight function
IL
W is
IL
IL
W (5.89)
The functions in the atomistic partition of unity
2
P enforce energy conservation of the
local atom L in the lattice of the atomistic configuration.
As is common practice in the use of homogenization techniques, the macro and
micro weight function representation are not used simultaneously. The terms in the
asymptotic analysis producing macro-governing equation, the macroscopic weights are
used. In the terms producing micro-governing equation, the microscopic weights are
utilized.
Similarly in this application of asymptotic analysis, the macroscopic partition of unity
1
P
is selected in the formulation which produces macro-equation, and the atomistic partition
of unity
2
P is selected when the development produces the governing atomistic equation.
To implement this switching, define a positive real parameter , 0 such that
83
1
1
sgn PP (5.90)
2
2
1sgn PP (5.91)
Then define the final partition of unity using (5.90) and (5.91)
12 PP P (5.92)
or
12
sgn 1 sgn PP P (5.93)
Another partition of unity is useful in this analysis. Suppose two unique atom positions
I and J are considered. The following partition of unity can be defined in a way similar
to the approach described previously.
3
1
1
1
2
N
IL JL
L
PWW
(5.94)
where
IL
IL
JL
JL
W
W
Then, in a way similar to that used in forming equation (5.91), let
3
3
1sgn PP (5.95)
and let
13 OP P (5.96)
84
or
13
sgn 1 sgn OP P (5.97)
Selectively inserting the partitions of unity from (5.93) and (5.97) in (5.84)
,
oo
oo o
mi m i
VV
i
mi
mmm m
mVV V
d
CTPdV SvOdV
dt
Sv
PdV q n PdA f v PdV
X
(5.98)
A continuum multi-scale version of (5.98) in terms of continuous micro-scale spatial co-
ordinates Y and temporal coordinate can be defined as follows:
,
111
11 11
11
oo
oo
o
mi m i
VY V Y
i
mi
m m
m VY VY
mm
VY
C
TPdY d dV S v OdY d dV
tY Y
Sv
PdY d dV q n PdY d dA
YX Y
fv PdY d dV
Y
(5.99)
Consider the first term in (5.99) and introducing (5.65)
1,
12
11
1,
12 1 2
11
1
11 1
11
o
o
oo
VY
IK MN
KI
V
IK MN
KI
VV
C
TPdY d dV
tY
TT
CPPdV
tM N
TCT
C P PdV P PdV
tMN
(5.100)
85
The proper partition of unity based weighting strategy can be defined by the specification
of parameter . Set equal to 0 in the macroscopic term, 0 in the atomistic term in
(5.100) and introduce
1
P and
2
P from (5.85) and (5.88) to give
1,
1111
1
11
o
oo
VY
IK RNMN
IL
L
LLKI
VV
C
TPdY d dV
tY
TCT
CWdV WdV
tMN
(5.101)
This modification leads to a summation and a consideration of all discrete atoms. This
then is consisted with the macroscopic term, where all possible spatial positions X
are
considered through the volume integration.
Consider now the stress power taken in (5.99) and introducing the parameter connecting
line elements 0,1
1
,
0
1
13
,
0
11
11
o
o
mi m i
V
mi m i
V
SvO YdddV
Y
Sv PPYdddV
Y
(5.102)
Substituting the atomistic expression for the stress power density from (5.25), we get
1
,
0
,,
1
1,
,
111
0
13
11
11
1
o
o
mi m i
V
IJ K J I IJ K
mi
MNN
o
IJ K
mi
m
KIJ
V
IJ IJ
i
SvO YdddV
Y
Fy yY
dV
u
u
MY
PP d
tY
(5.103)
86
It should be noted that this expression involves a summation over atom pair I and J . The
formulation represents an assembly of the contributions of the individual atom pairs to
the total stress power. Since assembly is involved, it is natural that the partition of
unity
3
P be used and the partition of unity
2
P which appears in equation (5.100).
Simplifying (5.103)
1
,
0
1
, ,,
13
111
0
,,
1,
11
0
13
11
11 1
2
11 1
1
2
o
o
mi m i
V
o
MNN
mi IJ K J I IJ K
mi
KIJ
V
IJ
IJ K J I IJ K
mi
NN
IJ K
m
IJ
IJ IJ
i
SvO YdddV
Y
u
Fy yY P Pd dV
MY t
Fy yY
u
MY
PP d
Y
1
1
o
M
K
V
dV
(5.104)
Again setting equal to zero in the macroscopic term and setting 0 in the atomistic
term in (5.104) it can be seen that
1
,
0
, ,,
11 11
,,
1,
11
11 1
2
11 1
11
2
2
o
o
mi m i
V
o
RM NN
mi IJ K J I IJ K
L
mi
LK IJ
V
IJ
IJ K J I IJ K
mi
IJ K
IL JL
m
J
IJ
i
SvO YdddV
Y
u
Fy yY W dV
MY t
Fy yY
u
MY
WW
Y
11 11
o
RM NN
LK I
V
IJ
dV
(5.105)
87
Now consider the heat flux term in (5.99) and using Green’s theorem it can be seen that
,
11
11 11
o
oo
i
mi
m VY
ii
mi m mi m
VY V Y
Sv
PdY d dV
YX
Sv PdY d dV Sv nPdY d dA
YY
(5.106)
Then consider the first term on the right had side in (5.106) and using (5.92) as the
partition of unity, it can be seen that
,1,2,
11 11
oo
ii
mi m mi m m
VY V Y
Sv PdY d dV Sv P P dY d dV
YY
(5.107)
The flux can be defined as follows:
1
0
IJ
i
i
i
mi mi
u u
Sv S
t
(5.108)
Where
1
1, 1 ,
1
2
IJ
IK J K
i
ii
uu u
or,
IJ
i
i
mi mi
u
Sv S
(5.109)
88
Where the first term in (5.107) can be rewritten using (5.109) and the symmetry in the
stress tensor as follows
,,
11 11
o o
IJ IJ
im
mi m mi i
VY VY
uu
S P dY d dV S P dY d dV
YY
(5.110)
Using the expression for the stress tensor from (5.25) and using (5.108) along with
the parameter
1
,
0
,,
1
1
0
111
0
1, 2,
11
11 1
2
o
o
IJ
m
mi i
V
IJ K J I IJ K
mi
MNN
IJ
m
m
KIJ
V
ii
IJ
u
SPYdddV
Y
Fy yY
dV
u u
MY
PP d
t
(5.111)
Then setting equal to zero in the microscopic term and setting 0 in the atomistic
term in (5.111), it can be seen that
1
,
0
1
0
,,
,
11 11
,,
1
0
11
11 1
2
11 1
2
o
o
IJ
m
mi i
V
IJ
RM NN
m
IJ K J I IJ K m
Li
mi
LK IJ
V
IJ
IJ K J I IJ K
mi
IJ
m
m
u
SPYdddV
Y
u u
Fy yY WdV
MY t
Fy yY
u u
MY
t
11 11
1
o
RM NN
IJL
LK IJ
V
IJ IJ
i
W
dV
Y
(5.112)
Where,
IJL JL IL
WW W (5.113)
89
Substituting (5.101), (5.112), and (5.105) in (5.98)
1,
1111
, ,,
11 11
1,
,,
1
11 1
2
1
11 1
2
oo
o
IK RNMN
IL
L
LLKI
VV
o
RM NN
mi IJ K J I IJ K
L
mi
LK IJ
V
IJ
IJ K
IJ K J I IJ K m
mi IJ
i
TCT
CWdV WdV
tMN
u
Fy yY W dV
MY t
u
Fy yY
Y
MY
11 11
1
0
,,
,
11 11
,,
1
2
11 1
2
11 1
2
o
o
RM NN
IL JL LK IJ
V
IJ
IJ
RM NN
m
IJ K J I IJ K m
Li
mi
LK IJ
V
IJ
IJ K J I IJ K
mi
dV
WW
u u
Fy yY WdV
MY t
Fy yY
MY
1
0
11 11
1
o
IJL
IJ
RM NN
m
m
IJ
LK IJ
i V
IJ
u uW
dV
tY
(5.114)
Using the asymptotic expansion, and considering the 1 O terms, the following macro
energy equation is obtained.
1
, ,,
11 11
1
0
,,
,
111
11 1
2
11 1
2
o
o
R
L
L
V
o
RM NN
mi IJ K J I IJ K
L
mi
LK IJ
V
IJ
IJ
MNN
m
IJ K J I IJ K m
Li
mi
KIJ
IJ
T
CWdV
t
u
Fy yY W dV
MY t
u u
Fy yY WdV
MY t
1
o
R
L
V
(5.115)
90
Using the asymptotic expansion, and considering the
1
O
terms, the following macro
energy equation is obtained.
1,
11 1
1,
,
11 11
1
0
,
1
1
11 1 1
22
11 1
2
o
o
IK NM N
IL
LK I
V
IJ K RM NN
IL JL
IJ K J I m
m
LK IJ
V
IJ
IJ
IJL
m
IJ K J I m
m
J
I
CT
WdV
MN
u
Fy y W W dV
MY
u u
Fy y WdV
MY t
11 1
o
RM NN
LK I
V
J
(5.116)
1,, LN
5.4 Results and Discussion
5.4.1 Problems with a Specified Temperature Distribution
Initially consider a class of problems in which the temperature is fixed and heat
transfer effects are precluded. In this problem class, neither the macro energy Eq. (5.115)
nor the Eq. (5.116) micro energy equations are incorporated in the model. The system
response depends only on the solution of the macro and micro equations of motion. Two
different subclasses will be considered that differ in terms of the way the macro model
behaves.
91
The first subclass involves quasi-static thermal problems in which the inertia term
in the macro equation of motion is omitted. The resulting macro problem is static, while
the atomistic problem is dynamic and is dependent on the static deformation gradients
resulting from the macro solution. Examples will be presented in which both constrained
and unconstrained thermal expansion occurs.
The second subclass involves dynamic thermal problem in which the inertia term
in the macro equation of motion is included. Both the macro and micro problems are
dynamic, while the atomistic problem is dependent on the macro deformation gradients
which vary with macro time t.
5.4.1.1 Quasi-Static Thermal Problems
The quasi-static thermal problem is governed by the following macro and micro
equations resulting from (5.40) and (5.44). By introducing the macro shape functions
and nodal force vector P in (5.40), the following equation is obtained.
0
0
20,
2
,,
11 1
11
P
PL i
V
L MN N
IJ K J I IJ K L
ij i
V
KI J
j
JI
u
dV
t
Fy yY dVP
MY X
(5.117)
92
The nonlinear coupling in the problem resulting from (5.117) and (5.44), includes
complicated interactions involving histories of atomistic displacements in fast time,
Virial stresses, and macroscopic deformations. To achieve convergence and reasonable
performance in solving the coupled equations, the following dynamic relaxation version
of the algorithm has been utilized:
0
,,
11 1
11
L L MN N
IJ K J I IJ K L i
ij i
V
KI J
j
JI
du
Fy yY dVP
dt M Y X
(5.118)
where is the dynamic relaxation parameter. The equation (5.118) can be rewritten as in
Eq. (5.119).
0
,1 ,
,,
11 1
11
Ln Ln
L MN N
IJ K J I IJ K L
ij i
V
KI J
j
JI
uu
T
Fy yY dVP
MY X
(5.119)
Where T the relaxation time step which should not be confused is should not be
confused with the macro time step t and fast time step . A number of multi-scale
thermo-mechanical test problems have been formulated and solved to demonstrate the
workability of the developed algorithms.
93
5.4.1.1.1 Quasi-Static Problems with Uniform Temperature
To demonstrate the workability of the dynamic relaxation procedure, a 3-D
problem involving unconstrained thermal expansion was initially considered. A single
element macro model, as pictured in Figure 4.6, is sufficient to solve this problem. A
single eight node hexahedron finite element was used in the macro simulation model. By
applying the appropriate macro node displacement constraints, the model was allowed to
freely expand in all three global coordinate directions due to the induced thermal stress.
Figure 5.1: One 3-D Hexahedron element model
A 2 x 2 x 2 Gauss Quadrature integration scheme was used to evaluate the macro volume
integral in (5.119). At each gauss point, the atomic model was based on a 4 x 4 x 4 FCC
atomic lattice structure. A constant temperature of 30K was maintained at all gauss points
during relaxation; therefore there was no heat transfer in the system.
1
4
2
3
5
6
7
8
94
Figure 5.2 – 5.4 shows iterative convergence behavior for strains, stresses, and
displacements at typical gauss points. From these Figures it can be seen that the three
normal strain components reach the same number as the algorithm converges. The fact
that the individual strain components behave differently during the relaxation process is
related to the particular macro displacement restraints applied at the nodes and their lack
of 3-D symmetry. Since the element was allowed to freely expand under thermal stress
and since no external mechanical forces were applied, the element expands uniformly
until it reaches the zero stress state. It can be seen that the dynamic relaxation procedure
produces a stable iterative process, a reasonable convergence rate, and an accurate macro
displacement solution. All of these aspects are important when attempting to simulate the
full coupled thermo-mechanical problem including heat transfer, as discussed in the
subsequent sections.
0
0.005
0.01
0.015
0.02
0.025
0.03
0 0.1 0.2 0.3 0.4 0.5
Dynamic Time Relax ation
Strain
Stress 11
Stress 22
Stress 33
Figure 5.2: Quasi-static strain versus dynamic time relaxation with uniform temperature
distribution
95
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
0 0.1 0.2 0.3 0.4 0.5
Dynamic Time Relaxation
Stress
Strain 11
Strain 22
Strain 33
Figure 5.3: Quasi-static stress versus dynamic time relaxation with uniform temperature
distribution
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
0 0.1 0.2 0.3 0.4 0.5
Dynamic Time Relax ation
Displacement
Displacement x
Displacement y
Displacement z
Figure 5.4: Quasi-static displacement versus dynamic time relaxation with uniform
temperature distribution
96
5.4.1.1.2 Quasi-Static Problems with Non-Uniform Temperature
The purpose of the next example is to demonstrate the satisfactory performance of
the dynamic relaxation procedure in solving quasi-static problems in which the macro
model is subjected to a non-uniform temperature distribution and the free thermal
expansion of the macro model is prevented because of displacement constraints.
As pictured in Figure 5.5, the macro model is composed of four hexahedron finite
elements. The displacements were restrained at both ends of the structure, and nodes 5
through 16 were constrained to displace only in the direction of the length of the beam.
As before, a 2 x 2 x 2 Gauss Quadrature integration scheme was used to evaluate the
macro volume integral in (5.119). As stated before the material representation at each
Gauss point utilized a 4 x 4 x 4 FCC lattice. The non-uniform nodal temperature
distribution in the model involved a linear temperature variation between the bar end
points with a temperature T1 = 0.3 at the left hand side and T2 = 0.4 at the right hand side.
97
Figure 5.5: Four element constrained 3-D macro finite element model
Figure 5.6-5.8 shows iterative convergence behavior for displacements at the center plane
and stresses and strains at the Gauss points closest to right end of the model. From these
Figures it can be seen that the nodes in the center plane displace in response to the non-
uniform temperature distribution.
T1
T2
1
16
17
18
2
3
4
5
6
7
8
9
10
11
12
13
15 19
20
14
1
2 3 4
98
-0.16
-0.14
-0.12
-0.1
-0.08
-0.06
-0.04
-0.02
0
0 100 200 300 400 500
Dynamic Time Relaxation
Displacement
Figure 5.6: Quasi-static displacement versus dynamic time relaxation with non-uniform
temperature distribution
-4.7
-4.5
-4.3
-4.1
-3.9
-3.7
-3.5
0 100 200 300 400 500
Dynamic Time Relax ation
Stress 11
Figure 5.7: Quasi-static stress versus dynamic time relaxation with non-uniform
temperature distribution
99
0
0.001
0.002
0.003
0.004
0.005
0.006
0.007
0.008
0 100 200 300 400 500
Dynamic Time Relaxation
Strain 11
Figure 5.8: Quasi-static strain versus dynamic time relaxation with non-uniform
temperature distribution
5.4.1.2 Dynamic Thermal Problems
The dynamic thermal problem is governed by the following temporally discrete
version of the macro and micro equations resulting of motion from (5.40) rewritten in
(5.120) and (5.44) respectively.
0
0
,1 , , 1
2
,,
11 1
2
11
Ln Ln Ln
PL
V
L MN N
IJ K J I IJ K L
ij i
V
KI J
j
JI
uu u
dV
t
Fy yY dVP
MY X
(5.120)
To demonstrate the workability of the dynamic relaxation procedure, the 3-D
problem involving unconstrained thermal expansion, considered above and pictured in
Figure 5.1, was solved one more, this time in conjunction with a full macro dynamic
100
model. The element is maintained at constant temperature of 30K. The initial macro
displacements and velocities were set to zero. The initial micro solution involves atom
positions consistent with the zero temperature equilibrium solution and atom velocities
equal thermal velocities corresponding to the fixed 30K macro temperature. The
computed dynamic solution is expected to oscillate about the quasi-static dynamic
relaxation solution obtained previously. Since the macro model has no damping, an
oscillatory solution without exponential decay is anticipated.
Figure 5.9 and Figure 5.10 show typical stress and strain behavior respectively, at each
gauss point. Fig. 10 it can be noticed that the initial stress state in compression does not
reach the same magnitude in tension, this is attributed to the fact the L-J potential is not
symmetrical about the equilibrium position. Figure 5.11 shows the time history
displacement of node 7.
-0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0 0.5 1 1.5 2 2.5
Macro Time
Strain
Strain 11
Strain 22
Strain 33
Figure 5.9: Dynamic strain history plot with uniform temperature distribution
101
-3
-2
-1
0
1
2
3
00.5 11.5 2 2.5
Macro Time
Stress
Stress 11
Stress 22
Stress 33
Figure 5.10: Dynamic stress history plot with uniform temperature distribution
-0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
00.5 11.5 2 2.5
Macro Time
Displacement
Displacement x
Displacement y
Displacement z
Figure 5.11: Dynamic distance history plot with uniform temperature distribution
102
5.4.2 Coupled Thermo-Mechanical Model with Heat Transfer
In this section the temperature in the model is allowed to vary and heat transfer
effects are included. In this problem class the macro energy (5.115) is incorporated in the
model and not the micro energy (5.116). The system response not only depends on the
solution of the macro and micro equations of motion but also on the macro temperature
field. In the current section, working of the coupled thermo-mechanical simulation
algorithm would be discussed and implemented. The solution obtained from the
algorithm would be benchmarked against Molecular Dynamics solution. The algorithm
for implementation is demonstrated in Table 5.1.
Table 5.1: Quasi-static thermo-mechanical multiscale algorithm
Step (I) Discretize and assemble the macro finite element domain.
Step (II) Initial nodal temperatures, displacement and velocities are assigned to the macro
domain.
Step (III) Solve the dynamic relaxation Eq.(5.119) to obtain the macro nodal
displacement and velocities at the next macro time step. The dynamical relaxation
equation is restated without nodal applied forces.
0
,1 ,
,,
11 1
11
Ln Ln L MN N
IJ K J I IJ K
ij
V
KI J
j
JI
uu
Fy yY dV
MY X T
Where n is the macro time step and K is the micro time step.
103
Step (IV) Interpolated macro parameters such as deformation gradient, velocity and
temperature at a gauss point is used to initialize the atomic lattice suited at that particular
gauss point. Follow atomistic momentum Eq. (5.44) is initialized
Where,
,1,
,
,
,1,
1, ,
,
21 , 2 1 ,
22
JKn J Kn J oIJ n
kl kl l kk
J Kn J Kn JJ JKn oIJ n
kkl kl kk
JK J K
kk
Y yy
u
YY yy
uu
uy
1J
k
y - The atomistic thermal displacement
Step (V) Thermal Eq. (5.115) is solved for the updated nodal temperature in the macro
model.
,1 ,
1
, ,,
11 11
1
0
,,
,
11
11 1
2
11 1
2
o
o
Pn Pn R
PL
L
V
o
RM NN
mi IJ K J I IJ K L
mi
LK IJ
V
IJ
IJ
NN
m
IJ K J I IJ K L m
mi i
IJ
IJ
TT
CdV
t
u
Fy yY dV
MY t
u u
Fy yY dV
MY t
11
o
RM
LK
V
The heat flux term appearing on the right had side of the equation is computed
using Non-Equilibrium Molecular Dynamics techniques by imposing the macro
temperature gradients.
Step (VI) Repeat from Step (2) till the desired temperature of the system is reached.
104
The coupled thermo-mechanical model problem chosen to solve is the one as
shown in Fig. 5. The equivalent molecular dynamic model, used in the test cases,
involves a perfect 42 x 8 x 8 FCC lattice, consisting of 10752 atoms with a density of
1.06 in L-J units. The results from both the multiscale model and the molecular dynamics
system are compared and studied.
Since there is no heat loss or heat input into the system, the temperature is
expected to equilibrate. It can be seen in Figure 5.12 that the nodal temperature from the
initial state at macro time zero gradually approaches the average temperature of the
system. The same trend is experienced by nodes on the same plane. Since the nodes in
element 3 and 4 were initially at a higher temperate than that of nodes in element 1 and 2,
the former set of elements were in a state of tension and the latter in compression. This
can be observed from Figure 5.13, which shows typical nodal displacement of the center
plane and as the system temperature equilibrates the displacement approaches zero.
Though not shown here, at the equilibrium state all the nodal displacement approached
zero. The rate of equilibration was verified using molecular dynamics. From Figure
5.14a-d it can be seen that there is a good agreement between the multiscale model and
the molecular dynamics results.
105
-0.16
-0.14
-0.12
-0.1
-0.08
-0.06
-0.04
-0.02
0
0 100 200 300 400 500 600 700
Macro Time
Mid Node Displacement
Figure 5.12: Quasi-static mid-plane nodal displace versus macro tome with non-uniform
temperature distribution
0.25
0.3
0.35
0.4
0.45
0.5
0.55
0 100 200 300 400 500 600
Macro Time
Nodal Temperature
Node 1
Node 5
Node 9
Node 13
Node 17
Figure 5.13: Quasi-static nodal temperature versus macro tome with non-uniform
temperature distribution
106
0
0.1
0.2
0.3
0.4
0.5
0.6
0 5 10 15 20 25
Distance
Temperature
Molecular Dynamics - 0 LJ Time
Nodal Temperature - 0 LJ Time
Figure 5.14a: Quasi-static nodal temperature versus molecular dynamic results with non-
uniform temperature distribution at 0 macro time
0
0.1
0.2
0.3
0.4
0.5
0.6
0 5 10 15 20 25
Distance
Temperature
Molecular Dynamics - 20 LJ Time
Nodal Temperature - 20 LJ Time
Figure 5.14b: Quasi-static nodal temperature versus molecular dynamic results with non-
uniform temperature distribution at 20 macro time
107
0
0.1
0.2
0.3
0.4
0.5
0.6
0 5 10 15 20 25
Distance
Temperature
Molecular Dynamics - 60 LJ Time
Nodal Temperature - 60 LJ Time
Figure 5.14c: Quasi-static nodal temperature versus molecular dynamic results with non-
uniform temperature distribution at 60 macro time
0
0.1
0.2
0.3
0.4
0.5
0.6
0 5 10 15 20 25
Distance
Temperature
Molecular Dynamics - 100 LJ Time
Nodal Temperature - 100 LJ Time
Figure 5.14d: Quasi-static nodal temperature versus molecular dynamic results with non-
uniform temperature distribution at 100 macro time
108
CHAPTER 6
ATOMISTIC NON-LINEAR VIBRATIONS
6.1 Dynamic Behavior in the Fast Time Scale
The purpose of this chapter is to develop a method for solving a particular
dynamic problem. This is the problem of atomistic lattice vibrations due to thermal
effects. These lattice vibrations can cause a dependence of material properties on
temperature and result in a dependence of the macro solution on temperature. A
procedure for defining the material property dependence on temperature will be
introduced here.
Essentially the problem to be solved is an atomistic free vibration problem. A multi-scale
version of the dynamic continuum equation in space and time is,
0 00
2
0000
,, 2
11 1
I
n
mnik n k i mn m
YY VV V
u
Cu WdYddV M WddV fWdV
Yt
(6.1)
Where is the period of the periodic temporal atomistic solution andY is the atomistic
region of interest, and
2
mn
mn
M
uu
109
Where,
mn
M - The mass matrix of the structure
Using (2.20), the multiscale version of (6.1) in semi-discrete form is obtained:
0 0
0
2
00
,, 2
11
00
11
11 1
1
I NN
IJ
IJ IJ IJ IJ I n
mn
ik nk i mn
IJ
Y VV
IJ
NN
IJ
m
IJ
V
IJ
u
CY Y u W d dV M W d dV
Nt
fW dV
N
(6.2)
The analysis can be started with the semi-discrete equation of motion for one period as
presented in equation
0
0
0
,,
11
2
0
2
11
1
0
NN
IJ
IJ IJ IJ IJ
mn
ik nk i
IJ
V
IJ
I
I n
mn
Y V
CY Y u W d dV
N
u
MWddV
t
(6.3)
Note that since the objective is a free vibration analysis, the body force term in (6.2) has
been omitted. Introducing (2.17), (2.33) in (6.3)
1
,0 ,
11
221
0
22
1
11 1
.
11 1
0
o
o
IJ
IJ
NN
IJ n oIJ IJ IJ
ik nk i mn IJ IJ
IJ
K i
IJ
oI
N
I
I nn
mn
I
W
u
W
uddv YY
C
v
NY
Y
uu
MWXW
ddv
v
Nt
(6.4)
110
Now consider a free vibration of the lattice about an unstrained macro-state. Separating
out the
1
terms and setting the global weight function
() WX to zero, gives the
following equation governing the atomistic free vibration problem.
21
2
1
1
11
11
11
0
I
N
I
I n
mn
I
NN
IJ
IJ
IJ
n
mn
IJ
IJ
u
MWd
N
Wd
u
C
N
(6.5)
6.2 Solution for the Non-Linear Free Vibration Problem
Then using the Taylor series expansion, for the L-J potential, from (2.39) and setting the
deformation gradient to be zero,
2() 3 () 4()
111
1
2
IJ E IJ E IJ E
bn bn bn
IJ IJ IJ IJ
lli mn IJ IJ IJ IJ IJ IJ IJ IJ IJ
mn m n l m n l i
Er Er Er
uuu
C
rr rrr rrrr
(6.6)
111
1
2
IJ IJ IJ
IJ IJ IJ
mn mnli mnl
lli
DGu H uu (6.7)
Introducing (2.40) in (6.5)
21
2
1
1111
11
11
()
11 1
()
2
0
I
N
I
I n
mn
I
NN
IJ IJ IJ IJ
IJ IJ IJ IJ
mn mnli mnl
llin
IJ
IJ
u
MW d
N
DGu H uu u W d
N
(6.8)
111
or,
21
2
1
1
11
11
11
111
11
()
11
()
11
()
11 1
(
2
I N
I
I n
mn
I
NN
IJ IJ
IJ
mn
n
IJ
IJ
NN
IJ IJ
IJ IJ
mnl
ln
IJ
IJ
IJ IJ
IJ IJ IJ
mnli
ln i
u
MW d
N
Du W d
N
Gu u W d
N
Hu u u W
N
11
)
0
NN
IJ
IJ
d
(6.9)
Now introduce a harmonic temporal variation in fast time to produce a
temporally discrete model,
11
11
sin( )
sin( )
II
nn
IJ IJ
nn
uU
uU
(6.10)
where
1I
n
u is the eigenvector and is the natural frequency. The temporally discrete
weight function is,
1
1
sin( )
sin( )
II
IJ IJ
WW
WW
(6.11)
Inserting (6.10) and (6.11) into (6.9) and considering the situation at a particular
macroscopic position and gauss point, the following equation is obtained:
112
21 2
1
12
11
11 3
11
11 1 4
11
11
sin ( )
11
sin ( )
11
sin ( )
11 1
sin ( ) 0
2
N
I
II
mn n
I
NN
IJ
IJ IJ
mn n
IJ
IJ
NN
IJ
IJ IJ IJ
mnl l n
IJ
IJ
NN
IJ
IJ IJ IJ IJ
mnli l n i
IJ
IJ
MU W d
N
DU W d
N
GU U W d
N
HU U U W d
N
(6.12)
The time integrals in (6.12) can be evaluated as follows:
2
3
4
11
sin ( )
2
1
sin ( ) 0
13
sin ( )
8
ad
bd
cd
(6.13)
Then substituting (6.13) in (6.12) and replacing
2
by
1
1
1
11
11
11
11 1
11
1
0
2
N
I
II
mn n
I
NN
IJ
IJ IJ
mn n
IJ
IJ
NN
IJ
IJ IJ IJ
mnl l n
IJ
IJ
NN
IJ
IJ IJ IJ IJ
mnli l n i
IJ
IJ
a
MU W
N
a
DU W
N
b
GU U W
N
c
H UUU W
N
(6.14)
113
Now this equation can be linearized by setting,
1
11
1
11
I
II
n
nn
IJ
IJ IJ
n
nn
UU U
UU U
(6.15)
where
1I
n U is the know eigenvector and
1I
n
U is the eigenvector increment associated with
the eigenvector,
1IJ
n U is the know eigenvector and
1IJ
n
U is the eigenvector increment
associated with the displacement difference between the atoms in the IJ pair, is the
know eigenvalue and is the eigenvalue increment.
1
1
1
1
1
11
11
11
11
111
111
11
1
0
2
N
II
II
n
mn n
I
NN
IJ IJ
IJ IJ
n
mn n
IJ
IJ
NN
IJ IJ IJ
IJ IJ IJ
ln
mnl l n
IJ
IJ
NN
IJ IJ IJ IJ
IJ IJ IJ IJ
ln i
mnli l n i
IJ
IJ
a
MU U W
N
a
DU U W
N
b
GU U U U W
N
c
H U UU UU U W
N
(6.16)
Neglecting the higher order terms in (6.16), we obtain the following matrix equation to
define the connection to the eigenvector and eigenvalue.
114
111
1
11
1
1
1
1
1111
1
13
2
2
1
2
NN
IJ IJ IJ IJ
IJ IJ IJ IJ
lli
mn mnl mnli n
IJ
IJ
N
I
II
mn n
I
N
II
I
n
mn
I
N
IJ IJ IJ IJ IJ
IJ IJ IJ
llin
mn mnl mnli
J
IJ
c
aD bG U H U U U W
N
a
MUW
N
a
MU W
N
c
aD bG U H U U U W
N
1
1
1
N
I
N
II
I
n
mn
I
a
MU W
N
(6.17)
This represents 3N equations in 31 N unknowns. This equation must be
supplemented by an amplitude constraint equation which can be associated with the
temperature T at a macroscopic point.
6.3 Finite Temperature Amplitude Constraint
The temperature is associated with the time averaged kinetic energy of the
lattice
1I
n
A u .
11
1
1
11
2
II N
II nm
nmn
I
uu
Au M d
N
(6.18)
Using the continuum in time temporal behavior association from (6.10), it can be seen
that (6.18) takes the following form:
115
111
1
,
2
N
IIII
nmnnm
I
a
AU MUU
N
(6.19)
Then the macroscopic temperature
T , considered as the time average of the atomistic
solution over one period of fast time oscillation, is
1
,
2
I
n
B
TAK U
DK
(6.20)
Where D is dimension of the problem,
B
K is the Bottzman constant, from equation (6.19)
And to replace the macroscopic temperature
T by a scaled temperature variableT .
11
1
2
N
II I
mn n m
I
Ta
TMUU
N
(6.21)
This equation acts as an amplitude constraint on the eigenvector magnitude in terms of an
average over the lattice. Essentially it normalizes the vibration mode.
Equation (6.21) must be linearized by an incremental process similar to the one
used in forming equation (6.16). The T must be incremented in accordance with the
following rule.
TT T (6.22)
116
Where T is a fixed temperature increment in the increment temperature algorithm. It
should be noted that T has a fixed value and is similar to an increment external force in
an incremental force algorithm. Introducing (6.15) and (6.22) in (6.21), we obtain,
11
11
1
2
N
II
III
nm
mn n m
I
a
TT M U U U U
N
(6.23)
Linearizing (6.23), it can be seen that,
1
1
1
11
1
11
1
2
2
N
I
II
n
mn m
I
N
II
I
mn
mn
I
N
II
I
mn
mn
I
a
MU U
N
a
MU U
N
a
M UU T T
N
(6.24)
Equation (6.24) and (6.17) represent 31 N equations in 31 N unknowns. It should be
noted that the system composed of equations (6.24) and (6.17) does not involve a
symmetric coefficient matrix. Such a system can be produced by multiplying (6.17) by
to obtain.
117
111
1
11
2
1
1
1
1
1111
1
3
2
2
2
NN
IJ IJ I IJ
IJ IJ IJ IJ
lli
mn mnl mnli n
IJ
IJ
N
I
II
mn n
I
N
II
I
n
mn
I
IJ IJ IJ IJ IJ
IJ IJ IJ
llin
mn mnl mnli
J
I
c
aD bG U H U U U W
N
a
MUW
N
a
MU W
N
c
aD bG U H U U U W
N
1
21
1
NN
I
J
N
II
I
n
mn
I
a
MU W
N
(6.25)
Equation (6.24) and (6.25) produce a 31 N x 31 N system of equations with a
symmetric coefficient matrix and the same solutions as (6.17) and (6.25).
6.4 Incremental Temperature Algorithm
The implementation of the Incremental Temperature Algorithm (Table 6.1), presented in
(6.24) and (6.25), involves a step-by-step incremental-iterative process. First the linear
eigenvalue problem for the system is solved to obtain linear eigenvalues and eigenvectors
for the atomic lattice. Then the behavior of the individual eigenvalues and associated
eigenvectors is studied, one eigenvalue at a time. In particular, the nonlinear behavior of
the individual natural frequencies of the system is determined as a function of
temperature.
118
Table 6.1: Incremental temperature algorithm
Linear eigenvalue problem:
Step (a) Compute the mass matrix and Hessian matrix
Step (b) Solve the linear eigenvalue problem
1
0 KM U
Non-Linear eigenvalue problem:
Step (I) For a chosen linear eigenvalue , scale the
1I
m
U eigenvector to a small starting
value and call it
1I
m U
Step (II) Compute
11
1
2
N
II
I
mn
mn
I
a
TMUU
N
Step (III) Assemble the coefficient matrix and residual vector from Eqs. (6.24) and (6.25)
Step (IV) Solve for
1I
n
U and
Step (V) Update the eigenvector
1
11
I
II
n
nn
UU U
Step (VI) Update the eigenvalue
Step (VII) Continue if residual tolerance , else return to Step (III)
Step (VIII) Increment T by T until maximum scaled temperature is reached
Step (IX) Scale eigenvector
1I
n U to produce scaled temperature T and return to Step (III)
Step (X) Repeat from Step (I) starting with a different linear eigenvalue
119
6.5 Test Problem
The atomistic model, used in the test cases, involves a perfect FCC lattice,
consisting of 108 atoms (N=108) with a density of 0.8 in L-J units.
In the finite temperature condition, atoms are constantly vibrating about an
equilibrium position. The associated equilibrium bond length y at finite temperature can
be determined by minimizing the Helmholtz free energy H, through the following
calculation:
0
H
r
(6.26)
3
1
() ln 2sinh
4
N
n
B
n
B
h
Hr TK
KT
(6.27)
where h is the Planck constant 6.63 x 10
-34
Js
.
The thermo-mechanical constitutive parameters can be defined through the
natural frequencies of the atomic lattice and the Helmholtz free energy. The ITA
procedure defined in Section 6 was used to define the frequency variation with
temperature. Figure 6.1 presents the evolution of a set of typical frequencies with
temperature.
120
0
0.1
0.2
0.3
0.4
0.5
0 5 10 15
Frequency
Temperature
7.56
7.21
6.48
5.69
2.71
Figure 6.1: The evolution of starting frequencies (2.71, 5.69, 6.48, 7.21, 7.56) as a
function of temperature
From the results in Figure 6.1, it can be seen that there is a significant dependence of the
vibration frequency on temperature. Using the computed variation of frequency with
temperature, it is possible to compute the specific heat
V
C and thermal expansion with
temperature.
6.5.1 Calculation of Specific Heat at Constant Volume
These are higher frequencies which harden as a function of temperature.
The specific heat
V
C at constant volume is defined as the heat energy required per unit
volume of solid per degree of temperature change. The specific heat
[11] is given by the
following expression:
121
2
3
2
1
exp
2
2
1exp
2
n
N
B n
vB
n
B
n
B
h
KT h
CK
KT
h
KT
(6.28)
In Eq.(6.28), the summation is over the 3N lattice natural frequencies. The ITA algorithm
can define the temperature variation of all 3N frequencies for use in Eq.(6.28). However,
in order to create a more computationally efficient algorithm, a modification to Eq.(6.28)
can be introduced. When the linear eigenvalues are the same for multiple modes, the
expression in Eq.(6.28) can be approximated. The developed model replaces the
summation over 3N vibration frequencies by a summation over the M unique linear
frequencies. For each of the individual frequencies
n
, from the set of M unique
frequencies, a weight
n
R , which accounts for the number of occurrence of a particular
frequency in the set of 3N lattice frequencies, is defined. This modification is particularly
useful in analyzing a perfect lattice and results in the following reduced version of
Eq.(6.28).
2
2
1
()
exp
2 ()
2
()
1exp
2
n
M
B n
vB n
n
B
n
B
hT
KT hT
CK R
KT
hT
KT
(6.29)
It was determined that the 108 atom FCC lattice has 24 unique frequencies
(M=24). Therefore instead of studying all 324 (3N) frequencies, only the frequencies in
the set of M were studied. Using the appropriate weights and frequencies obtained from
122
the ITA model, the specific heat was computed, as a function of temperature, using
equation (6.29). To validate the developed method, the specific heat was also calculated
using the Molecular Dynamics (MD) method.
The atomistic system used in the MD simulation was the same as the one used in
the extended quasi-harmonic model, with respect to number of atoms, density, and
interatomic potential. The MD specific heat [1] expression was computed, for a micro
canonical ensemble, using the following expression:
2
22
3 3
1
22
B
B
V
K
kNKT
C
(6.30)
where the following definitions apply:
- Average over operator
k - Kinetic energy of the system
2
k -
2
kk
The MD simulation was performed for 60,000 time steps at each temperature and the
results were averaged over the last 20,000 time steps.
123
Figure 6.2 compares the specific heat computed using the extended quasi-
harmonic model and the results of the MD method. As show in Figure 6.2, the results
obtained from Eq. (6.30) at higher temperatures compares well with the molecular
dynamics results, within the latter’s error bar indicated by their fluctuations. The
anharmonic terms used in the approximation of the potential in the extended quasi-
harmonic model represent a significant contribution to the computed specific heat. In the
next section, the specific heat results are also verified in comparison to experimental
results using analytical functions to model the frequency-temperature space.
0
5
10
15
20
25
020 40 60
Temperature (Kelvin)
Specific Heat (J / K / Mole)
Extended Quasi Harmonic Model
Molecular Dynamics
Figure 6.2: The specific heat calculated from Eq. (6.29) as compared to the MD results
obtained from Eq. (6.30)
124
6.5.2 Calculation of Thermal Expansion
The calculation of thermal expansion relies on the determination of the
equilibrium bond length y at finite temperature, as presented in Eqs. (6.26) and (6.27). As
before, for the perfect lattice, the 3N summation in the Helmholtz free energy H can be
replaced by a summation over the M unique frequencies of the system, with the number
of occurrences being indicated by
n
R . The reduced version of Eq. (6.27) can be written as
follows:
1
,
() ln 2sinh
4
M
n
Bn
n
B
hrT
Hr TK R
KT
(6.31)
The frequency in the above expression depends on both the lattice constant and the
temperature. To find the equilibrium bond length y for arbitrary temperatures, it would
be necessary to evaluate the frequency dependence on temperature at an infinite number
of lattice constants. Instead of calculating the temperature dependence of frequency at
every lattice constant, the frequency evolution, with respect to temperature, is calculated
at L different lattice constant values. Then a frequency surface is generated from the L
computed values using an analytical model. The needed frequency data can be
interpolated from the frequency surface for arbitrary values of the lattice constant and
temperature.
To develop a frequency surface, a regression function must be selected to fit the
generated frequency data. The behavior of the frequency , as a function of temperature
125
T , for fixed values of y, was found to have strong square root dependence on
temperature. The behavior of the frequency , as function of the lattice constant, for
fixed values of
T , was found to have quadratic dependence. Thus to model the frequency
surface, an analytical function of the following form was utilized:
2
12 3 4 5
, yTA Ay Ay A T AT (6.32)
Consider an example case in which the parameter L is set to 5, and the lattice
constants, at which the frequency dependence on temperature is to be computed, are 1.39,
1.49, 1.59, 1.64, and 1.71, presented in L-J units. A frequency surface, representing the
dependence of frequency on
T and lattice constant, can be created as shown in Figure 3.
0
0.2
0.4
0.6
0.8
1 1.3
1.4
1.5
1.6
1.7
1.8
0
20
40
60
80
Temperature
Lattice Constant
Frequency
Figure 6.3: Typical frequency surface, involving data at the lattice constants 1.39, 1.49,
1.59, 1.64, and 1.71, presented in L-J units
126
The associated constants, in equation (6.32), were determined for all frequencies in the
set M. The modeling error, associated with the difference between the frequency surface
in Eq. (6.32) and the data in the actual model, was found to be acceptable. These
analytically modeled frequency surfaces are introduced into the free energy expression
presented in Eq. (6.33). As a result, the free energy can be computed at any required
temperature. Then the thermal expansion is computed by minimizing the free energy to
find the equilibrium lattice parameter at a particular temperature
T .
The results on thermal expansion at different temperatures have been validated by
comparison with the experimental results for Argon obtained by Peterson et al. [39]. The
thermal expansion results are presented in Figure 6.4. There is a 1.5% discrepancy in the
computed lattice constant at
0 T . A complete fit could have been obtained by choosing
different lattice parameters. By shifting the curve to coincide with the experiment results
at
0 T , the lattice constants, obtained by minimizing the free energy and using the
extended quasi-harmonic method, lie within 0.15% of the experimental results.
127
1.55
1.56
1.57
1.58
1.59
1.6
1.61
1.62
1.63
0 0.2 0.4 0.6 0.8 1
Temperature
Lattice Constant
Extended Quasi Harmonic Model
Experimental
Figure 6.4: Thermal expansion results obtained from the extended quasi-harmonic model
as compared to the experimental results obtained by Peterson et al. [39]
An added advantage, associated with the creation of the frequency surface model defined
in Eq.(6.34), is that the frequencies at any required lattice constant and finite temperature
can be obtained by interpolation. To see how this frequency surface model could be used,
consider the computation of the variation of the specific heat with temperature at a fixed
lattice constant value. Using the extended quasi-harmonic method and constructing the
associated surface model, the specific heat, at a particular lattice constant, can be
determined. In Figure 6.5, the so determined computational results, for a lattice constant
of 1.56 in L-J units, are compared to the experimental results of Peterson et al. [39].
128
0
5
10
15
20
25
0 2040 6080
Temperature ( Kelvin )
Specific Heat (J/Mole/ K)
Extended Quasi Harmonic Model
Experimental
Figure 6.5: Specific heat results obtained from the extended quasi-harmonic model as
compared to the experimental results of Peterson et al. [39]
129
CHAPTER 7
IMPLICIT TIME INTEGRATION ALGORITHM
The paper presents an implicit time integration scheme for defining the dynamical
behavior of atomistic systems. The developed procedure incorporates the anharmonic
terms in the inter-atomic potential and thus can capture the nonlinear behavior of the
atomistic system. The defined time integration algorithm satisfies energy conservation
principles on a time step by time step basis.
7.1 Weighted Residual Formulation
The Lagrangian equation of motion for a three dimensional continuum
0
V is
defined in terms of the reference coordinate X
, the deformed position vector x
, the
equilibrium starting position vector
() E
X
, and the displacement vector relation to the
equilibrium position ( ) ux . The Lagrangian equation of motion is defined as follows:
22
2
qi
q
uS
X
(7.1)
130
Where
qi
S is the first Piola-kirchhoff stress tensor and is the atomistic fast time scale.
Defining a weight function WX
and assuming periodic boundary on the surface
0
V ,
the following weighted residual version of Eq.(7.1) is obtained:
00
22
2
qi
q VV
uW
WdV S dV
X
(7.2)
Consider the inertial term in Eq.(7.2) defined as follows:
0
2
1 2
i
V
u
I WdV
(7.3)
In the atomistic system, an equivalent form of Eq.(7.3) is defined as follows:
0
2
1 2
1
I N
I i
I
I
V
u
I mXXWXdV
(7.4)
Where,
N - The number of atoms in the volume
0
V
I
m - The mass of the atom I
- The Dirac Delta function
I
i
u - The displacement vector of atom I .
131
Simplifying Eq.(7.4), it can be seen that
2
1 2
1
I N
I i
I
I
u
I mW
(7.5)
Where,
II
WWX
is the atomistic weight function associated with atom I .
Consider now the stiffness term in (7.2) defined as follows:
0
2
2
qi
q V
W
I SdV
X
(7.6)
The fact that the model is atomistic can be incorporated by defining as continuum
network L of the line segment.
IJ
L connecting atom I and J for all atomistic pairs in the
lattice.
,1
N
IJ
IJ
JI
LL
The line segment where
IJ
L is parameterized by , [0,1]. On
IJ
L
2
1
1
IJ
IJ J I
IJ
rr
X XX
XX X
XX
(7.7)
132
The continuum weight functions WX
in Eq.(7.6) is arbitrary. The atomistic problem is
characterized by the following restricted form of WX
:
:;,1,
0:
IJ
IJ
IJ
WX X L IJ N
WX
XL
(7.8)
Where
1
IJ J I
WWW (7.9)
The Piola-Kirchhoff stress can be defined as follows,
1
NN
qi qiIJ
IJI
SS
(7.10)
Introducing 0.7, 0.8, and 0.10 in 0.6, it can be seen that
1
2
1
0
IJ NN
qiIJ IJ
q
IJI
W
IS X
X
(7.11)
Using a process similar to the one defined by Irving et al. [29] and Evans et al. [14], but
formulated in a weak form consistent with 0.2, it can be shown that the following viral
type stress representation defines
qiIJ
S
1
qiIJ IJ IJ
iq
IJ
SFX
X
(7.12)
Here the atom force
IJ
i
F is
133
IJ IJ IJ IJ
IJ
i IIJ
ii
x x
F
XX
(7.13)
where
IJ IJ
x
is the potential with the atomistic pairs and
IJ J I
x xx
Introducing (7.9) and (7.12) in (7.11), it can be seen that is
2
1
NN
IJ IJ J I
i
IJI
I Fx W W
(7.14)
The atomistic weighted residual equation of motion, consistent with the continuum
Lagrangian weighted residual equation of motion Eq. (7.2) , can be obtained by setting
1
I
from Eq. (7.5) equal
2
I from Eq. (7.14)to obtain:
2
2
11
0
I NNN
IIJIJJI i
Ii
IIJI
u
mW Fx W W
(7.15)
Consider the typical atom L, 1 L N . Then define the fine scale weights using the
Kroneckor delta, as follows:
,
,
,
IK
IL
JK
JL
IJ K
JLIL
W
W
W
(7.16)
134
Introducing Eq. (7.16) in Eq. (7.15), it can be seen that
21,
,
2
111
0
IK NNN
IJ K J I m
IIL i JLIL
IIJ
JI
u
mFxx
(7.17)
Or simplifying Eq.(7.17), the following atom L equation of motion is obtained. The
resulting atomistic equation of motion takes the following form:
21 ,
,,
2
11
0
LKNN
IL K J I LJ K J I m
Li i
IJ
IL J L
u
mFxx Fxx
(7.18)
1,
1, ,
LN
KM
7.2 Approximation for Atomistic Force
The atomistic force can be can be expanded in a Taylor series, including anharmonic
terms, to give
11
26
IJ IJ IJ IJ IJ IJ IJ IJ IJ IJ
iinn inlln ilmnlmn
F D u G uu H u uu (7.19)
Where,
The relative atomistic displacement vector is
IJ J I
nn n
uu u
135
Represents the change in the displacement vector relative to the displacement in the
equilibrium positions
() IJ E
X
, and
2()
3()
2()
IJ IJ E
IJ
in IJ IJ
in
IJ IJ E
IJ
inl IJ IJ IJ
il n
IJ IJ E
IJ
ilmn IJ IJ IJ IJ
il m n
x
D
XX
x
G
XXX
x
H
X XXX
Then introducing (7.19) in (7.14)
2
1
11
26
NN
IJ IJ IJ IJ IJ IJ IJ IJ IJ J I
in n inl l n ilmn l m n
IJI
I Du G u u H uu u W W
(7.20)
And setting
12
I I from (7.5) and (7.20)
2
2
11
11
26
I NNN
I IJIJ IJ IJIJ IJ IJ IJ IJ J I i
Iinninllnilmnlmn
IIJI
u
m W D u Guu H uuu W W
(7.21)
7.3 Energy Preserving Algorithm
Following the works presented by Wellford et al. [48, 49] for solving non-linear
structural dynamics problem, Eq. (7.21) is temporally discretizied as follows
2 1 , 1 1, 1, 1
22
2
IIK IK IK
i
uu u u
(7.22)
136
11,11, 1,1
11 1
42 4
IJ IJ IJ IJ K IJ K IJ K
mn n mn n n n
Du D u u u
(7.23)
1 1 1, 1, 1 1 , 1 , 1
1
3
IJ IJ IJ IJ IJ K IJ K IJ K IJ K
mnl n l mnl n l l l
Gu u G u u u u
(7.24)
1 1 1 1, 1, 1, 1 1 , 1
111
222
IJ IJ IJ IJ IJ IJ K IJ K IJ K IJ K
mnli n l i mnli n l i i
H uuu H u u u u
(7.25)
Plugging in Eq.(7.22)-(7.25) in Eq. (7.21) and obtain Eq. (7.26) which is a spatially and
temporally discretizied atomistic equation of motion.
1, 1 1 , 1 , 1
1, 1 1 , 1 , 1
2
1
1, 1, 1 1 , 1 , 1
1
1 , 1, 1, 1 1
1
2111
42 4
1
3
1
4
IK IK I K N
IIJ
IIJIJKIJKIJK
mn mn n n n
J
IJ
N
IJ
IJ IJ K IJ K IJ K IJ K
mnl n l l l
J
IJ
N
IJ IJ K IJ K IJ K IJ
mnli n l i i
J
IJ
uu u
mWDuuuW
G u uu uW
Huu u u
,1
0
IJ
K
W
(7.26)
By regrouping the terms in Eq. (7.26) such that the
1, 1 IJ K
n
u
terms only appear on the right
side of Eq. (7.27)
137
1, 1 1 , 1 , 1 , 1 , 1
2
1
1, 1 1 , 1
2
1, 1 , 1
1
1, 1, 1
1
111 1
43 4
1
2
11
24
1
3
N
I IJ
I I K IJ IJ IJK IJ IJ K IJK IJ K
mn mn mln n miln i l n
J
IJ
IJ
IIJK IJK
mn n n
N
IJ
IJ IJ K IJ K
mn l l
J
IJ
N
IJ IJ K IJ K
mnl n l
J
IJ
mu W D G u H u u u W
Mu u W
Du u W
Gu u
1, 1
1, 1 , 1, 1
1
1
0
4
IJ
IJ K
l
N
IJ
IJ IJ K IJ K IJ K
mnli n l i
J
IJ
uW
Huu u W
(7.27)
It can be shown from Eq. (7.27) that the total energy E at time step
1
2
K is equal to the
energy at time step
1
2
K by proving each term in the energy expansion (kinetic and
potential) is the same at every half time step. The total energy E at time step
1
2
K is
expressed in Eq. (7.28) as follows,
11 1 1 1
23 4
22 2 2 2
KK K K K
E
(7.28)
Where,
1
2
K
is the kinetic energy and
11 1
234
22 2
,,
KK K
is the potential energy of the second,
third and fourth order expansion of potential respectively at time step
1
2
K .
As stated earlier we can see from Eq. (7.29) that energy at time step
1
2
K is equal to the
energy at time step
1
2
K .
138
11
22
11
22
1
0
KK
KK
EE
or
EE
(7.29)
To illustrate Eq. (7.29)we prove the conservation of kinetic and third order term of the
energy. Let represent the kinetic part of the energy in Eq. (7.30)
1, 1 1 , 1 , 1
2
2
IK IK IK
I
I
mn
uu u
mW
(7.30)
Let the weight be the atom I velocity as represented in Eq. (7.31)
1, 1 1 , 1
2
IK I K
Iuu
W
(7.31)
and substituting Eq. (7.31) in Eq. (7.30) we get,
1, 1 1 , 1 , 1 , 1 1, 1 1 , 1 , 1 , 1
3
2
I
IK I K IK IK IK IK IK I K mn
nn n n n n n n
m
u u uu u u uu
(7.32)
and we rearrange the terms in Eq. (7.32) to get Eq. (7.33)
1 , 1 1 , 1 , 1 1, 1, 1, 1
3
1
2
I IKIK IKIK I IK IK
mn n n m m mn m m
Mu u u u M u u
(7.33)
139
Note:
1, 1 1 , 1 , 1 , 1 1, 1, 1 1 , 1 1 ,
0
I I K I K I K I K I I K IK IK IK
mn n n m m mn n n m m
mu u u u m u u u u
Because
II
mn nm
M M mass matrix symmetric, we can see that at time step
1
2
K is same
as at time step
1
2
K .
11
,,
22
1
IK IK
I
KK
t
1 1, 1 1 , 1 , 1 1,
,
2
1
2
IK I K IK IK
IK
I nn m m
mn
uu uu
KM
(7.34)
Likewise a similar derivation can be obtained for the third order term expressed as .
1, 1 , 1 1 , 1 , 1
1
1
3
N
IJ
IJ IJ IJ K IJ K IJ K IJ K
mnl n l l l
J
IJ
Gu u u u W
(7.35)
and substituting velocity for weight in Eq. (7.35)
1, 1 , 1
1
2
IJ
IK IK
nn
Wuu
(7.36)
we get Eq. (7.37)
1, 1, 1 1 , 1 , 1 1, 1, 1
1
1
6
N
IJ IJ IJ K IJ K IJ K IJ K I K I K
mnl n l l l n n
J
IJ
Gu u u u u u
(7.37)
140
and terms in Eq. (7.37) can be rearranged to show at at time step
1
2
K is same as at time
step
1
2
K .
11
,,
33
22
1
IJ K IJ K
IJ
(7.38)
Where,
1, 1, 1 1 , 1
1
,
3
2
1, 1, 1, 1
1
6
IJ K IJ K IJ K
nl m
IJ K
mnl
IJ K IJ K IJ K
nl m
uu u
G
uuu
(7.39)
Thus Eq. (7.29) holds true for every half time step. The proposed atomistic equation of
motion as presented in Eq. (7.27) is termed as the Energy Preserving Algorithm (EPA) in
subsequent sections.
7.4 Result and Discussions
As mention earlier EPA uses a potential approximation up to the fourth order by
using a Taylor series expansion. We hope to characterize the behavior of atoms at high
temperature by inclusion of these anharmonic terms. As demonstrated earlier the
developed EPA is an implicit algorithm that would be unconditionally stable at every half
time step. We would also like to examine the inaccuracy that may have been caused by
141
omitting the higher order terms. The stability and accuracy of the energy preservation
algorithm is studied and the results are compared to explicit Velocity Verlet algorithm.
To demonstrate the workability of the algorithm, we consider a 4 x 4 x 4 face centered
cubic lattice. The total number of atoms in the system is 256. An initial temperature of 40
K is assigned to the lattice system at rest. Lennard-Jones inter-atomic potential is used to
model the interaction among atoms.
7.4.1 Stability of the Algorithm
The stability of the implicit EPA would be the key strength over the explicit Velocity
Verlet algorithm. In order to analyze the stability, we need to compare algorithms that
would use the same approximation of the potential. The conventional explicit Velocity
Verlet force algorithm is modified to include the Taylor series force approximation up to
the fourth order. This would determine that the dynamic equations of motion solved in
both the algorithms are the same.
In order to capture the high frequency thermal oscillation of atoms, a time step of
0.005 is chosen to march forward in time. Typical properties studied for this type of
analysis would be the average kinetic energy, temperature, potential energy and the total
energy of the atomistic system. Since the model does not include any heat transfer
phenomenon such as a heat bath or sink, we would expect the energy to be conservative
at every time step.
142
Figure 7.1 and 7.2 shows the computed kinetic and strain energy of the system for
EPA and Velocity Verlet respectively using a time step of 0.005. It can be seen that the
total energy is conserved for every time step, and the potential and kinetic energy
fluctuating complimentary to one another. Figure 7.3 shows the temperature results from
EPA and Velocity Verlet compared for the same initial temperature condition and time
step. It can be observed that both algorithms produce identical results for every time step.
The temperature fluctuation can be averaged over the entire duration to determine the
temperature of the system. Considering the time step used is relatively small enough to
capture the high oscillation of atoms, the obtained results are used a datum for further
analysis purpose.
0
0.1
0.2
0.3
0.4
0.5
0.6
0 0.5 1 1.5 2 2.5
Time
Energy
Kinetic Energy
Potential Energy
Total Energy
Figure 7.1: EPA Energy results at Density 1.1, Time Step 0.005 and Initial Temperature
40K
143
0
0.1
0.2
0.3
0.4
0.5
0.6
0 0.5 1 1.5 2 2.5
Time
Energy
Kinetic Energy
Potential Energy
Total Energy
Figure 7.2: Velocity Verlet Energy results at Density 1.1, Time Step 0.005 and Initial
Temperature 40K
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0 0.5 1 1.5 2 2.5
Time
Temperature
EPA
Velocity Verlet
Figure 7.3: EPA and Velocity Verlet Temperature results at Density 1.1, Time Step 0.005
and Initial Temperature 40K
144
Next we study both the algorithms by increasing the time step to 0.01 but with the
same initial condition. Figure 7.4 and 7.5 shows the computed kinetic and strain energy
of the system for EPA and Velocity Verlet respectively using a time step of 0.01. Though
both the algorithms seem to conserve the total energy, the error in the energy computed
from Velocity Verlet increased by square of time step used. Figure 7.6a shows the
temperature results from EPA and Velocity Verlet compared for the same initial
conditions and time step. It can be observed that EPA shows a slight deviation with
respect to the datum results obtained at time step 0.005 and Verlet algorithm failing to
capture the high frequency oscillation. From Figure 7.6b it can be seen that EPA with
time step 0.01 agrees well from the results at time step 0.005.
0
0.1
0.2
0.3
0.4
0.5
0.6
0 0.5 1 1.5 2 2.5
Time
Energy
Kinetic Energy
Potential Energy
Total Energy
Figure 7.4: EPA Energy results at Density 1.1, Time Step 0.01 and Initial Temperature
40K
145
0
0.1
0.2
0.3
0.4
0.5
0.6
0 0.5 1 1.5 2 2.5
Time
Energy
Kinetic Energy
Potential Energy
Total Energy
Figure 7.5: Velocity Verlet Energy results at Density 1.1, Time Step 0.01 and Initial
Temperature 40K
0.05
0.1
0.15
0.2
0.25
0.3
00.5 1 1.5 22.5
Time
Temperature
EPA
Velocity Verlet
Figure 7.6a: EPA and Velocity Verlet Temperature results at Density 1.1, Time Step 0.01
and Initial Temperature 40K
146
0.05
0.1
0.15
0.2
0.25
0.3
00.5 1 1.5 22.5
Time
Temperature
EPA - 0.01
EPA - 0.005
Figure 7.6b: EPA at Time 0.005 and EPA at Time 0.01 Temperature results at Density
1.1 and Initial Temperature 40K
The stability of the algorithm is test by further increasing the time step to 0.05 (10
times the original time size used). From Figure 7.7 it can be found that the explicit
Velocity Verlet algorithm completely breaks down at time step 0.05 and the implicit EPA
is still stable and conserves energy at every time step. Though there is a little compromise
in the accuracy due to the use of larger time step, in this type of analysis the averaged
properties over the time integral are of importance and not the exact time history trace.
147
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0 0.5 1 1.5 2 2.5
Time
Temperature
EPA - 0.05
EPA - 0.005
Figure 7.7: EPA at Time 0.005 and EPA at Time 0.05 Temperature results at Density 1.1
and Initial Temperature 40K
7.4.2 Accuracy of the Algorithm
In order to analyze the accuracy, we need to compare EPA to an algorithm that does not
have any approximation of the potential. The conventional explicit Velocity Verlet which
uses the full potential in the force calculation is used to validate the accuracy of EPA.
This would help us determine the inaccuracy caused due to the omission of higher order
terms in the potential.
Figure 7.8 shows the temperature history results from EPA at time step 0.005
compared to the Velocity Verlet algorithm using no force approximation. It can be seen
there is a good agreement between both algorithms. The EPA’s temperature fluctuation is
in phase with Velocity Verlet, with EPA’s pecks slightly off. This determines the force
148
approximation introduced in implicit EPA is not significant enough to cause any
inaccuracy.
To further analyze the accuracy of EPA, we compare the atom displacement to the
Velocity Verlet. Figure 7.9 shows the time history trace of a particular atom obtained
from both algorithms. It can be observed the amplitude of vibration of the atom is in
correspondence. As discussed in the previous section the implicit EPA allows us to use
larger time step that would not be possible in Velocity Verlet algorithm. Given
unconditionally stable implicit nature of EPA and its accuracy, EPA can be used as a
effective method in solving atomistic momentum equation.
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0 0.5 11.522.5
Time
Temperature
EPA
Velocity Verlet
Figure 7.8: EPA and Velocity Verlet Temperature results at Density 1.1, Time Step 0.005
and Initial Temperature 40K
149
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
00.5 11.5 22.5
Time
Displacement
EPA
Velocity Verlet
Figure 7.9: EPA and Velocity Verlet atom displacement results at Density 1.1, Time Step
0.005 and Initial Temperature 40K
7.4.3 Additional Study
Additional study of EPA has revealed that the number of atoms in the system, the
lattice density and magnitude of temperature can determine the accuracy of the desired
result. In this section we demonstrate only the effect of temperature and density of the
system on EPA. Thermo-mechanical properties such as specific heat and thermal
expansion are highly sensitive to lattice constant and temperature changes. Jiang et al. [30]
have worked on studying these parameters effect on thermo-mechanical properties.
Given the same size, the face centered cubic crystal (4 x 4 x 4) with an initial
temperature increased to 80K (twice that of the temperature used before and almost close
to melting point of Argon) is studied. We compare the temperature results obtained from
150
EPA and Velocity Verlet in Figure 7.10. It can be observed that the temperature time
history trace obtained from EPA and Velocity Verlet do not agree as well as it did at 40K.
It can be inferred that higher temperature in the system is causing the neglected higher
order terms in the potential to contribute. The current state of thermal equilibrium has
larger amplitude of vibration that can be better characterized only by employing full
potential. Research work involved in the study of thermal equilibrium properties using
quasi harmonic approximation determines that the method is effective up to half the
melting point temperature of the crystal. There is definitely a temperature range, where
EPA could be effectively employed.
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
00.5 11.5 22.5
Time
Temperature
EPA
Velocity Verlet
Figure 7.10: EPA and Velocity Verlet temperature results at Density 1.1, Time Step
0.005 and Initial Temperature 80K
151
Next we try to determine the effect of lattice (density) in the working of EPA. The issue
of lattice constant in a way is related to magnitude of temperature in the system. It can be
seen in Figure 7.11 as we decreased the density of the system from 1.1 to 0.8 (but it
increases the lattice constant), for a initial temperature of 40K, the temperature results
from EPA and Velocity Verlet differ in the percentage of agreement. This is primarily
attributed to the characteristics of the potential employed; in this case it is Lennard-Jones
potential.
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0 0.5 1 1.5 2 2.5
Time
Temperature
EPA
Velocity Verlet
Figure 7.11: EPA and Velocity Verlet temperature results at Density 0.8, Time Step
0.005 and Initial Temperature 40K
152
CHAPTER 8
SUMMARY AND CONCLUSIONS
8.1 Static Analysis
The proposed algorithm demonstrates a method to link the macro finite element model
with the atomic unit cell using homogenization theory. Depending on the type of element
used in the macro model, there would be a number of Gauss integration points in each
element; and these Gauss integral points behave as the location of representative unit
atomic cells where the material elastic constants are extracted.
The dissertation already has illustrated the workability of the atomistic algorithm for
materials under strain and has obtained the material elastic constants and their behavior
as a function of strain. The results have excellent agreement with conjugate gradient
method.
Finally the multiscale algorithm for the static case links the atomistic scale calculation of
the constitutive parameter with the continuum scale finite element calculation. The key
aspect to the static multiscale algorithm is the inclusion of the anharmonic terms in the
material behavior and concurrent convergence of the solution at the macro and atomic
153
scale simultaneously. That ensures to capture the realistic behavior of materials under
strain. Another aspect of the algorithm is the weight residual formulation, which is
seamlessly implemented in the atomic scale using discrete weight functions. The
introduced weighed formulation is extended to dynamic cases as well. Though the
implementation is demonstrated for 2-D problems, its extension to 3-D is straight forward
without much modification.
8.2 Temperature Dependence of Material Properties
To obtain the material behavior under finite temperature condition, the dissertation uses a
harmonic approximation for the lattice vibration problem and incorporates anharmonic
terms in the Hessian approximation as initially introduced in the static case. A non-linear
eigenvalue problem is constructed and the lattice frequency dependence on temperature is
determined.
The obtained lattice frequency dependence on temperature is verified by studying
the variation of constitutive parameters such as co-efficient of thermal expansion and
specific heat under finite temperature condition. The presented computational procedure
uses atomistic models and fundamental statistical mechanics theories to obtain thermo–
mechanical constitutive parameters that incorporate the computed lattice frequencies.
154
The computed constitutive parameters account for anharmonic and finite
temperature effects. The developed algorithm can be used, in a stand alone mode, to
define the temperature dependence of macroscopic constitutive parameters. It can be used
as the atomistic scale solver in concurrent multi-scale thermo-mechanical simulations. It
can also be used in the preprocessor step of a serial multi-scale procedure to define
constitutive parameters which are fed to the macro-solver.
The performance of the developed algorithm has been shown to be good
in terms of comparisons to molecular dynamics and experimental results. The numerical
algorithm can be easily implemented, as it is quite similar to nonlinear incremental force
procedures employed in continuum nonlinear structural analysis. While the Lennard-
Jones potential in 3-D has been considered in the numerical examples, the incremental
temperature algorithm can be generalized to other potentials and materials.
Now consider possible extensions to this work. The computed examples have
considered a perfect atomistic lattice. This situation is a result of the availability of
experimental results for the perfect case. The algorithm and its applications could be
extended, by including the effects of defects and related uncertainty, to evaluate the effect
on the computed constitutive parameters. Could be used in probabilistic studies focusing
on effects of defects on constitutive parameters.
155
The theoretical development of the constitutive algorithm assumes a harmonic
time variation for the atom displacements. However, the developed algorithm is based on
a characterization of the nonlinearity inherent in the atomistic potential, and a more
accurate description of the time dependence of the atomistic displacement, more
consistent with nonlinear vibration theory, could be developed.
8.3 Thermo-Mechanical Equation
The proposed algorithms are particularly useful in analyzing thermo-elastic
problems, including thermal stress analysis and heat transfer cases. The algorithm has
been successfully implemented for quasi-static thermal problems in which the inertia
term in the macro equation of motion is omitted and dynamic thermal problem in which
the inertia term in the macro equation of motion is included. The dynamic relaxation
procedure introduced in the quasi-static test case produced simultaneous convergence at
both the macro and atomistic levels. It can be seen from the results presented that the
relaxation algorithm was consistent in solving constrained and unconstrained thermal
cases, as well as in the heat transfer problem.
One of the notable parts of the implementation is the Virial stress formulation
which is used in the macro momentum and energy equations and results in an algorithm
with no force approximations. At both the macro and atomistic scale, the derived energy
expressions, coupled with the relaxation algorithm, allow the modeling of heat transfer in
156
the system. The heat transfer results are in good agreement with molecular dynamics
results.
The next focus would be on developing efficient algorithms to compute
constitutive parameters coupled with parallel implementation. As a future extension of
the proposed algorithm, the study of material behavior for Graphite and Silicon, which is
has wide industrial applications, would be interesting. Though studying behavior of these
materials involves more complex inter-atomic potentials, the present formulation can be
easily adopted without much overhead. In addition, it would be interesting to formulate
the developed algorithms in a Eulerian setting, for application in solving fluid dynamics
problems.
8.4 Implicit Time Integration Algorithm
The proposed Energy Preserving Algorithm (EPA) is derived based on the continuum
Lagrangian equation of motion. As before the force term in the atomistic momentum
equation is approximated using a Taylor series expansion for the Hessian, allowing the
anharmonic effects. The algorithm is able to capture the high frequency atomistic
dynamics at elevated temperature conditions. EPA is found to be unconditionally stable
and energy preserving. The results show good comparison between EPA and the Velocity
Verlet method at small time steps. At larger time steps Velocity Verlet is shown to
breakdown while EPA is robust and stable.
157
In terms of future work, EPA is not restricted to modeling crystals characterized
by pair potentials; it can be employed for multi-body potential as well. Future study could
also incorporate the inclusion of parallel computations to increase the efficiency of the
implicit algorithms.
158
REFERENCES
[1] Allen, M., and Tildesley, D., “Computer Simulation of Liquids”, Oxford University
Press (1989).
[2] Belytschko, T., and Xiao, S.P., “Coupling Methods for Continuum Model with
Molecular Model”, International Journal for Multiscale Computational Engineering, 1(1),
115-126 (2003).
[3] Broughton, J.Q., Abraham, F.F., Bernstein, N., and Kaxiras, E., “Concurrent coupling
of length scales: Methodology and application,” Physical Review B, 60(4), 2391- 2403
(1999).
[4] Bugela, M., and Galliero, G., “Thermal conductivity of the Lennard-Jones fluid: An
empirical correlation”, Chemical Physics, 352(1-3), 49-257(2008).
[5] Chen, W., and Fish, J., “A mathematical homogenization perspective of virial stress”,
International Journal for Numerical Methods in Engineering, 67(2), 189 - 207 (2006).
[6] Chen, W., and Fish, J., “A generalized space-time mathematical homogenization
theory for bridging atomistic and continuum scales”, International Journal for Numerical
Methods in Engineering, 67(2), 253-271 (2006).
[7] Chung, P.W., and Namburu, R.R, “On a formulation for a multiscale atomistic-
continuum homogenization method”, International Journal of Solids and Structures, 40,
2563–2588 (2003).
[8] Chung, P.W., and Clayton, J.D., “Multiscale Modeling of Point and Line Defects in
Cubic Lattices”, International Journal of Multiscale Computational Engineering, 5, 203-
226 (2007).
[9] Cioranescu, D., and Donato, P., “An Introduction to Homogenization”, Oxford
University Press (2000).
[10] Cook, R., Malkus, D., Plesha, M., and Witt, R., “Concepts and Applications of Finite
Element Analysis”, Wiley (2001).
[11] De Wette, F.W., Fowler, L.H., and Nijboer, B. R. A., "Lattice dynamics, thermal
expansion and specific heat of a Lennard-Jones solid in the quasi-harmonic
approximation", Physica, 54(2), 292-304 (1971).
[12] Dib, G.M., Wellford, L.C., and Mindle, W., “Free and steady state vibration of non-
linear structures using a finite element-non-linear eigenvalue technique”, Earthquake
Engineering and Structural Dynamics, 8, 97-115 (1979).
159
[13] Ericksen, J.L., “Phase Transformations and Material Instabilities in Solids”,
Academic Press Inc. (1984).
[14] Evans, D.J. and Morriss, G., Statistical Mechanics of Nonequilibrium Liquids,
Cambridge University Press (2008).
[15] Fish, J., Chen, W., and Li, R., “Generalized mathematical homogenization of
atomistic media at finite temperatures in three dimensions”, Computer Methods in
Applied Mechanics and Engineering, 196, 908–922 (2007).
[16] Foiles, S. M., “Evaluation of harmonic methods for calculating the free energy of
defects in solids”, Physical Review B, 49, 14930-14938 (1994).
[17] Fosdick, L.D., Jessup, E. R., Schauble, C.J.C., and Domik, G., “Introduction to
High-Performance Scientific Computing”, The MIT Press (1996).
[18] Gates, T.S., Odegard, G.M., Frankland, S.J.V., and Clancy, T.C., “Computational
materials: Multi-scale modeling and simulation of nanostructured materials,” Composites
Science and Technology, 65, 2416–2434 (2005).
[19] Gear, C.W., Numerical Initial Value Problems in Ordinary Differential Equations,
Prentice-Hall (1971).
[20] Ghabriel, M. and Wellford, L.C., “An averaged Lagrangian-finite element technique
for the solution of nonlinear vibration problems”, Computers and Structures, 16, 1-4,
207-214 (1983).
[21] Golub, G., and Ortega, J.M., “Scientific computing: an introduction with parallel
computing”, Academic Press Professional, Inc. (1993).
[22] Guedes, J.M. and Kikuchi, N., “Preprocessing and postprocessing for materials
based on the homogenization method with adaptive finite element methods”, Computer
Methods in Applied Mechanics and Engineering, 83(2), 143-198 (1990).
[23] Hafskjold, B., Ikeshoji, T., and Ratkje, S.K., “On the molecular mechanism of
thermal diffusion in liquids”, Molecular Phyiscs, 80(6), 1389-1412 (1993).
[24] Hardy, R.J., “Formulas for determining local properties in molecular-dynamics
simulations: Shock waves”, Journal of Chemical Physics, 76, 622-628 (1982).
[25] Hoover, W.G., “Molecular Dynamics (Lecture Notes in Physics)”, Springer (1986).
[26] Huang, Z., and Tang, Z., “Evaluation of momentum conservation influence in non-
equilibrium molecular dynamics methods to compute thermal conductivity”, Physical B:
Condensed Matter, 373(2), 291-296(2006).
160
[27] Hughes, J. R., “The Finite Element Method: Linear Static and Dynamic Finite
Element Analysis”, Dover Publications (2000).
[28] Ikeshoji, T., and Hafskjold, B., “Non-equilibrium molecular dynamics calculation of
heat conduction in liquid and through liquid-gas interface”, Molecular Physics, 81(2),
251-261 (1994).
[29] Irving, J.H., and Kirkwood, J.G., “The statistical mechanical theory of transport
processes: IV. The equations of hydrodynamics,” The Journal of Chemical Physics, 18,
817–29 (1950).
[30] Jiang, H., Huang, Y., and Hwang, K.C., “A Finite-Temperature Continuum Theory
Based on Interatomic Potentials”, Journal of Engineering Materials and Technology,
127(4), 408-416 (2005).
[31] Jund, P., and Jullien, R., “Molecular-dynamics calculation of the thermal
conductivity of vitreous silica”, Physical Review B, 59, 13707-13711 (1999).
[32] Kalos, M.H., and Whitlock, P.A., “Monte Carlo Methods, Vol. I: Basics”, Wiley-
VCH (1986).
[33] Lesar, R., Najafabadi, R., and Srolovitz, D. J., “Finite-temperature defect properties
from free-energy minimization”, Physical Review Letters, 63, 624-627 (1989).
[34] Li, X., and E, W., “Multiscale modeling of the dynamics of solids at finite
temperature”, Journal of the Mehcanics and Physics of Solids, 53(7), 1650-1685(2005).
[35] Liu, B., and Qiu, X., “How to Compute the Atomic Stress Objectively?”,
Journal of Computational and Theoretical Nanoscience, 6, 1081–1089 (2009).
[36] Liu, X., and Li, S., “Nonequilibrium multiscale computational model”, The Journal
of chemical physics, 126, 124105 (2007).
[37] Mountain, R.D. and MacDonald, R.A., “Thermal conductivity of crystals: A
molecular-dynamics study of heat flow in a two-dimensional crystal”, Physical Review B,
28, 3022-3025 (1983).
[38] Park, H.S., Karpov, E.G., and Liu, W.K., “A temperature equation for coupled
atomistic/continuum simulations”, Computer Methods in Applied Mechanics and
Engineering, 193, 1713–1732 (2004).
[39] Peterson, O.G., Batchelder, D.N., and Simmons, R.O., “Measurements of X-Ray
Lattice Constant, Thermal Expansivity, and Isothermal Compressibility of Argon
Crystals”, Physical Review, 150, 703-711 (1966).
161
[40] Shenoy, V., Shenoy, V., and Phillips, R., Materials Research Society Symposium
Proceedings, 538, 465 (1999).
[41] Swope, W.C., Anderson, H.C., Berens, P.H., and Wilson, K.R., “A computer
simulation method for the calculation of equilibrium constants for the formation of
physical clusters of molecules: Application to small water clusters,” Journal of Chemical
Physics, 76, 637–649 (1982).
[42] Tadmor, E.B., Ortiz, M., and Philips, R., “Qusaicontinuum analysis of defects in
solids”, Philosophical Magazine A, 73(6), 1529-1563 (1996).
[43] Takanoa, N., Ohnishi, Y., Zako, M., and Nishiyabu, K., “The formulation of
homogenization method applied to large deformation problem for composite materials”,
International Journal of Solids and Structures, 37, 6517-6535 (2000).
[44] Verlet, L., “Computer "Experiments" on Classical Fluids. I. Thermodynamical
Properties of Lennard-Jones Molecules,” Physics Review, 159, 98-103 (1967).
[45] Waismana, H., and Fish, J., “A space–time multilevel method for molecular
dynamics simulations”, Computer Methods in Applied Mechanics and Engineering, 195
(44-47), 6542-6559 (2006).
[46] Wagner, G.J., Liu, W.K., “Coupling of atomistic and continuum simulations using a
bridging scale decomposition,” Journal of Computational Physics, 190, 249-274 (2003).
[47] Wang, W., Li, X., and Shu, C.W., “The Discontinuous Galerkin Method for the
Multiscale Modeling of Dynamics of Crystalline Solids”, Multiscale Modeling and
Simulation, 7(1), 94-320 (2008).
[48] Wellford, L.C, and Hamdan, S.M., “An Analysis of an Implicit Finite Element
Algorithm for Geometrically Nonlinear Problems of Structural Dynamics Part 1.
Stability,” Computer Methods in Applied Mechanics and Engineering, 14, 377-390
(1978).
[49] Wellford, L.C, and Hamdan, S.M., “An Analysis of an Implicit Finite Element
Algorithm for Geometrically Nonlinear Problems of Structural Dynamics Part 2.
Accuracy,” Computer Methods in Applied Mechanics and Engineering, 14, 391-399
(1978).
[50] Xiaoa, S.P., Belytschko, T., “A bridging domain method for coupling continua with
molecular dynamics”, Computer Methods in Applied Mechanics and Engineering, 193
(17-20), 1645-1669 (2004).
[51] Zhou, M., “A New Look at the Atomic Level Virial Stress: On Continuum-
Molecular System Equivalence Mathematical”, Physical and Engineering Sciences,
459(2037), 2347-2392 (2003).
162
APPENDIX: Higher Order Expansions of L-J Potential
2 ( )( ) ( )( )
"'
23 ()
() ()
11 1 1
11
ij ij ij ij NN NN
bsl lsls
bb IJ IJ ij
ij ij
IJ I J
ls
IJ IJ
Err rr
EE
rr N N r
rr
3
() ( ) ()
'''
3
()
11
() () ( ) () ( )
''
22 4
( ) () ()
11
() ( ) ()
''
() ( )
(
1
1
2
1
b
IJ IJ IJ
vl s
ij ij ij NN
vs l
b
ij
IJ
IJ
ij ij ij ij ij NN
ls slv
bsv lv
ij ij ij
IJ
IJ
ij ij ij
vls ls
b ij ij
ij
E
rrr
rrr
E
N
r
rr rrr
E
N
rr r
rrr
E
Nrr
r
3
)
11
() ( ) () () ( ) ()
'
33 3 5
() () () ()
11
1
3
NN
IJ
IJ
ij ij ij ij ij ij NN
vl s slv
bls sv lv
ij ij ij ij
IJ
IJ
rr r rrr
E
N
rr r r
4
b
IJ IJ IJ IJ
ls v g
E
TermITermIITermIII TermIV
rrrr
163
TermI
()
() ( ) () () ( ) ()
''' ''''
33 () ( )
() ()
11
() () ( ) () () ( )
''' '''
3
()
11
1
11
3
ij
ij ij ij ij ij ij NN
g
vs l v l s
bb ij ij
ij ij
IJ
g
IJ
ij ij ij ij ij ij
NN
vg l s lg v s gs v l
b
ij
IJ
IJ
r
rrr r rr
EE
rNr
rr
rr rr rr
EE
NN
r
() ( ) () ( )
5
()
11
ij ij ij ij
NN
vl s g
b
ij
IJ
IJ
rrrr
r
TermII
( ) () () () ( )
''
22 4 ()
( ) () ()
()
() ( ) () ( ) ()
'''
22 4 ()
( ) () ()
11
''
2
1
2
1
ij ij ij ij ij
ls slv
bsv lv ij
ij ij ij
g
ij
ij ij ij ij ij NN
g
ls slv
bsv lv ij
ij ij ij
IJ
IJ
gl
sv
b
r r rrr
E
r
rr r
r
r r rrr
E
Nr
rr r
r
E
N
() () ( ) ()
24 2 4
() () () ( )
() () () () () () ( ) () ( ) () 1
46
() ()
22
24
ij ij ij ij
lg sg s g
lv
ij ij ij ij
ij ij ij ij ij ij ij ij ij ij J
sgl v lgv s gvsl sl v g
ij ij
rr rr
rr r
rr rr rr rrrr
rr
1
NN
I
IJ
164
TermIII
()
( ) () ( ) () () ( ) ()
'' '''
324 () () ( ) ( )
() ( ) ()
11
() ( )
24
() ()
''
1
2
1
ij
ij ij ij ij ij ij ij NN
g
vls ls lsv vls
bb ij ij ij ij
ij ij ij
IJ
g
IJ
ij ij
gv v g
ls
ij ij
b
r
rrr rrrr
EE
rrr N r
rrr
rr
rr
E
N
( ) () () ( ) () () () ( ) () ( ) 11
46
() ()
4
NN
ij ij ij ij ij ij ij ij ij ij IJ
IJ vg l s lg v s sg v l s l v g
ij ij
rr rr rr rrrr
rr
TermIV
() () () () ( ) ()
'
33 3 5 ()
() () () ( )
()
() ( ) ( ) () () (
''
33 3 ()
() () ()
11
3
1
3
ij ij ij ij ij ij
vl s slv
bls sv lv ij
ij ij ij ij
g
ij
ij ij ij ij ij i NN
g
vl s slv
bls sv lv ij
ij ij ij
IJ
IJ
rr r rrr
E
r
rr r r
r
rr r rrr
E
Nr
rr r
)
5
()
() ()
35
() ()
() ( ) () ()
'
35 35
() () ( ) ( )
11
( ) () () (
3
33
1
3
j
ij
ij ij
gv v g
ls
ij ij
ij ij ij ij
NN
gv v g sg s g
bls lv
ij ij ij ij
IJ
IJ
ij ij ij ij
sg l v lg v s
r
rr
rr
rr rr
E
N
rr rr
rr rr
) ( )( ) ( ) ( )( ) ( )
57
() ( )
5
ij ij ij ij ij ij
vg s l s l v g
ij ij
rr rrrr
rr
Abstract (if available)
Abstract
The main objective of the dissertation is to develop multi-scale algorithms for continuum-atomistic problems. The focus is on sequential multi-scale simulations. In sequential multi-scale methods, the computations at the various scales are, in a sense, decoupled. This means, for example, that, for a continuum/atomistic simulation, large scale macroscopic continuum calculations rely on the results of fine scale computations and information obtained on an atomistic cell.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Multi-scale modeling of functionally graded materials (FGMs) using finite element methods
PDF
A complete time-harmonic radiation boundary for discrete elastodynamic models
PDF
An Experimental And Theoretical Approach To The Active Control Of A Shallow, Slack Cable Using Parametric Control Of Axial Tension
PDF
Homogenization procedures for the constitutive material modeling and analysis of aperiodic micro-structures
PDF
Stochastic and multiscale models for urban and natural ecology
PDF
Studies into vibration-signature-based methods for system identification, damage detection and health monitoring of civil infrastructures
PDF
Model, identification & analysis of complex stochastic systems: applications in stochastic partial differential equations and multiscale mechanics
PDF
Analytical and experimental studies of modeling and monitoring uncertain nonlinear systems
PDF
Analytical and experimental studies in modeling and monitoring of uncertain nonlinear systems using data-driven reduced‐order models
PDF
Continuum modeling techniques and their application to the physics of soil liquefaction and dissipative vibrations
PDF
Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems
PDF
Synchronization in carpets of model cilia: numerical simulations and continuum theory
PDF
Analytical and experimental studies in system identification and modeling for structural control and health monitoring
PDF
Algorithms for stochastic Galerkin projections: solvers, basis adaptation and multiscale modeling and reduction
PDF
Studies into data-driven approaches for nonlinear system identification, condition assessment, and health monitoring
PDF
Analytical and experimental studies in the development of reduced-order computational models for nonlinear systems
PDF
Characterization of environmental variability in identified dynamic properties of a soil-foundation-structure system
PDF
Towards an optimal design of feedback strategies for structures controlled with smart dampers
PDF
Studies into computational intelligence approaches for the identification of complex nonlinear systems
PDF
Hydrodynamic modeling and feasibility study of harnessing tidal power at the Bay of Fundy
Asset Metadata
Creator
Chockalingam, Karthikeyan
(author)
Core Title
Modeling of multiscale continuum-atomistic systems using homogenization theory
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Civil Engineering
Publication Date
01/25/2012
Defense Date
11/25/2009
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
atomistic-continuum,finite elements,molecular dynamics,multi-scale modeling,OAI-PMH Harvest,thermo-mechanical
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Wellford, Carter (
committee chair
), Ghanem, Roger Georges (
committee member
), Masri, Sami F. (
committee member
), Nakano, Aiichiro (
committee member
), Wong, Hung Leung (
committee member
)
Creator Email
chockali@usc.edu,karthi.chockalingam@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m2814
Unique identifier
UC1190829
Identifier
etd-Chockalingam-3442 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-288993 (legacy record id),usctheses-m2814 (legacy record id)
Legacy Identifier
etd-Chockalingam-3442.pdf
Dmrecord
288993
Document Type
Dissertation
Rights
Chockalingam, Karthikeyan
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
atomistic-continuum
finite elements
molecular dynamics
multi-scale modeling
thermo-mechanical