Weinan E. Principles of Multiscale Modeling

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

4.1. CONTINUUM MECHANICS 157

Assume that h is a function of the form:

h(x, t) = h

(F, T, x), (4.1.53)

then

∂

h = ∇

: ∂

F +

∂h

∂T

∂

T (4.1.54)

From (4.1.52), we also have

∂

h = ∂

e − s∂

T − T ∂

In terms of F, the energy conservation law (4.1.49) can be expressed as

∂

e = ∂

F : σ

− ∇

· q (4.1.55)

Using the Clausius-Duhem inequality

∂

s ≥ −

∇

· q +

q · ∇

T (4.1.56)

we get



∇

: ∂

F +

∂h

∂T

∂



= ρ

∂

h ≤ −ρ

s∂

T − q ·

∇

+ ∂

F : σ

, (4.1.57)

or equivalently,

∂

F : (ρ

∇

− σ

) + ∂



∂h

∂T

+ ρ



+ q ·

∇

≤ 0. (4.1.58)

Since F, T can take arbitrary values, the inequality (4.1.58) implies [23]

∂W

∂F

= σ

, (4.1.59)

s = −

∂h

∂T

q ·

∇

≤ 0.

where W = ρ

. These are constraints on the possible constitutive relations.

In the absence of additional information, the simplest thing to do is to assume a linear

relation between the remaining ﬂux q and the corresponding driving force ∇

T/T :

q = −k

∇

. (4.1.60)

Here k, the thermal conductivity, must be positive. The other constitutive information

needed is an expression for W , the stored energy density, which was discussed earlier.

158 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

4.1.5 Dynamics of ﬂuids

Next we turn to ﬂuids. As we said earlier, it is much more convenient to use the

Eulerian coordinates to describe ﬂuids. We will focus on dynamic models of ﬂuids. For

convenience, we now use u to denote the velocity ﬁeld (instead of the displacement ﬁeld

as was done earlier). In Eulerian coordinates, the basic principles take the form

• Conservation of mass:

∂

ρ + ∇ · J

= 0 (4.1.61)

• Conservation of momentum:

∂

(ρu) + ∇ · J

= 0 (4.1.62)

• Conservation of energy:

∂

E + ∇ · J

= 0 (4.1.63)

• The Clausius-Duhem principle, namely, the second law of thermodynamics:

ρ (∂

s + u · ∇s) + ∇ ·





≥ 0 (4.1.64)

where J

, J

are the mass, momentum and energy ﬂuxes, respectively, given by

= ρu, (4.1.65)

= ρu ⊗ u − σ,

= Eu − σ · u + q

It is convenient to introduce the material derivative. For any function ϕ, its material

derivative is deﬁned as

˙ϕ = D

ϕ =

∂

∂t

ϕ + u · ∇ϕ. (4.1.66)

Using material derivative, the conservation laws can be written as

˙ρ + ρ∇ · u = 0, (4.1.67)

u = ∇ · σ, (4.1.68)

ρ ˙e + ∇ · q = Tr



∇u



, (4.1.69)

ρ ˙s + ∇ · (

) ≥ 0. (4.1.70)

4.1. CONTINUUM MECHANICS 159

where Tr(A) denotes the trace of a matrix A. To obtain (4.1.69), multiply both sides of

(4.1.68) by u to get,

∂



|u|



+ ∇ ·



|u|

u − u · σ



= −Tr



∇u



(4.1.71)

Combining this equation and (4.1.63) gives (4.1.69).

Again, we deﬁne the Helmholtz free energy per unit mass as

h = e − T s, (4.1.72)

where T is the temperature of the system. We assume that h is only a function of T and

ρ, h = h

(T, ρ), then

h = ρ



˙e −

T s − T ˙s



(4.1.73)

Using (4.1.69), we get

h = Tr



∇u



− ∇ · q − ρs

T − ρT ˙s. (4.1.74)

On the other hand, we also have

h = ρ

∂h

∂T

T + ρ

∂h

∂ρ

˙ρ = ρ

∂h

∂T

T − ρ

∂h

∂ρ

(∇ · u) . (4.1.75)

Combining (4.1.73),(4.1.75), and (4.1.64), we get



σ +

∂h

∂ρ



∇u

− q ·

∇T

− ρ



s +

∂h

∂T



≥ 0. (4.1.76)

Since (4.1.76) should hold for any value of

T , we must have

s = −

∂h

∂T

(4.1.77)

which is a well-known result in thermodynamics. Now (4.1.76) becomes



σ +

∂h

∂ρ



∇u

− q ·

∇T

≥ 0 (4.1.78)

Let p =

∂h

∂ρ

be the thermodynamic pressure, and write

σ = −˜pI + σ

(4.1.79)

160 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

where ˜p =

Tr(σ), and Tr(σ

) = 0. (4.1.78) becomes

(p − ˜p)(∇ · u) + Tr



(σ

)

∇u



− q ·

∇T

≥ 0. (4.1.80)

This is the constraint on the constitutive relation from thermodynamic consideration.

In the absence of additional information, we again resort to a linear relation between

the ﬂuxes and the driving forces. This means that p − ˜p, σ

and q should all be linear

functions of ∇·u, ∇u and ∇T /T . It is easy to see that for isotropic ﬂ uids, the constitutive

relation must be of the form:

p − ˜p = κ∇ · u, (4.1.81)

= µ



∇u + ∇u

− (∇ · u)I



, (4.1.82)

q = −k

∇T

, (4.1.83)

Presence of other cross terms at the right hand side would violate isotropy or invariance

under reﬂection. Here κ, µ, k are constants, representing respectively the bulk viscos-

ity, shear viscosity and heat conductivity. Together with the conservation laws (4.1.67)-

(4.1.69), we obtain the familiar compressible Navier-Stokes equation in ﬂuid dynamics.

So far we have only discussed the simplest example of a ﬂu id system. For more

complicated systems, the principles are the same as described above. As an example, let

us consider a two-component ﬂuid system. Let C be the mass concentration of one of

the two components, then conservation of mass reads

C + ∇ · J

= 0, (4.1.84)

where J

is the diﬀusive mass current for the corresponding component. Now besides

T and ρ, the free energy also depends on C, i.e. h = h

(T, ρ, C). The Clausius-Duhem

principle gives

(p − ˜p)∇ · u + Tr



(σ

)

∇u



− q ·

∇T

− J

· ∇˜µ ≥ 0. (4.1.85)

4.2. MOLECULAR DYNAMICS 161

where ˜µ is the chemical potential ˜µ =

∂h

∂C

. In the linear resp onse regime, the constitutive

relation in this case is

p − ˜p = κ∇ · u, (4.1.86)

= µ



∇u + ∇u

− (∇ · u)I



, (4.1.87)

q = −A

∇˜µ − A

∇T

, (4.1.88)

= −B

∇˜µ − B

∇T

. (4.1.89)

This suggests that temperature gradient may drive mass diﬀusion, a phenomenon that

is commonly referred to as thermo-diﬀusion or the Soret eﬀect [12]. The coeﬃcients

, B

, A

, B

have to satisfy the constraint that the left hand side of (4.1.85) as a

quadratic form of (∇u, ∇˜µ, ∇T/T ) is semi-positive deﬁnite. Further constraints may

come from Onsager’s reciprocity relation [12], which in this case states that A

= B

To summarize, the basic steps for developing (dynamic) continuum models are as

follows:

1. Write down the conservation laws.

2. Find the constraints on the constitutive relations from the second law of thermo-

dynamics.

3. Find the constraints on the constitutive relations from symmetry considerations. In

most cases, as a ﬁrst approximation, one can postulate linear constitutive relations.

Such a general approach has its merits: It gives a framework that allows us to explore

all the possible leading order physical eﬀects, particularly cross eﬀects such as thermo-

diﬀusion.

4.2 Molecular dynamics

Next we turn to the atomistic models in the form of classical molecular dynamics. In

this section, we discuss the basic ideas of molecular dynamics, including the empirical

inter-atomic potentials, the diﬀerent ensembles, the continuum limits and the linear

response theory for computing transport coeﬃcients using atomistic modeling.

162 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

Consider a system of N atoms and denote by m

and y

the mass and position of the

j-th atom, respectively. Molecular dynamics models the evolution of this system using

Newton’s second law

= F

= −∇

V (y

, . . . , y

), j = 1, ··· , N (4.2.1)

where V is the potential energy function of the system. This system of equations can

also be written in Hamilton’s form (p

= m

= ∇

= −∇

H, j = 1, ··· , N (4.2.2)

where the Hamiltonian H is given by

H =

+ V (y

, . . . , y

). (4.2.3)

4.2.1 Empirical potentials

The ﬁrst issue we have to address is: What is the function V ? As we will see later when

discussing quantum mechanics, in principle V can be obtained from electronic structure

models. In practice, however, this is a rather costly procedure, and most models of

molecular dynamics still use empirical potentials that are carefully designed for speciﬁc

purposes.

It should be emphasized that at this point, there are no systematic ways of ﬁnding

empirical inter-atomic potentials. Current practice is a combination of coming up with

a good guess for the functional form of the potential and calibrating the parameters by

experimental data as well as data from ﬁrst principle-based calculations.

Empirical potentials for molecules

Potentials for molecules are usually in the following form:

V = V

bonded

+ V

non-bonded

= V

+ V

(4.2.4)

where V

is the contribution due to the covalently bonded-interactions, and V

is the

contribution due to the non-bonded interactions. Bonded interactions include contribu-

tions from the changes in bond length, bond angle and torsion ang le (also called dihedral

angle). In the example of propane (CH

− CH

, see Figure 4.2), there are 10

4.2. MOLECULAR DYNAMICS 163

bonds (2 C-C bonds, 8 C-H bonds), 18 angles (1 C-C-C angle, 10 C-C-H angle and 7

H-C-H angle), and 18 torsion angles (12 H-C-C-H angles and 6 H-C-C-C angles).

A typical form of V

bonds

− l

i,0

)

angles

(θ

− θ

j,0

)

torsion

angles

(cos ω

− cos ω

n,0

)

(4.2.5)

The ﬁrst term is the energy due to bond stretching, l is the bond length and l

is the

equilibrium bond length. The second term is due to the change in bond angles which

causes bending, θ

and θ

j,0

are respectively the j-th bond angle and its equilibrium

value. The third term is t he contribution from the change in torsion angles which causes

twisting. ω

is a torsion angle between three consecutive bonds, ω

n,0

is its equilibrium

value. Note that the functional form of the terms are rather simplistic, either quadratic

or simple cosines. The motivation is the same as before when we considered constitutive

relations for solids and ﬂuids: We always choose the simplest form that is consistent

with what we already know. Secondly, the parameters in (4.2.5) can be computed from

quantum mechanical calculations. This has become a standard practice when developing

force ﬁelds [5].

In principle, one might also consider mixed terms such as bond angle and torsion

angle interactions. But these are only included when it is absolutely necessary to do so.

Figure 4.2: An illustration of Propane.

better pics needed

Next we discuss the non-b onded interaction. There are two main contributions: elec-

164 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

Figure 4.3: An illustration of basic variables associated with covalent bonds.

trostatic and van der Waals interaction:

i6=j

4ε

(





−





)

i6=j

4πε

(4.2.6)

Here r

i,j

= |y

− y

|, q

is the partial charge of the ith atom, ε

is the permitivity

of vaccum. Their values need to be changed when solvent eﬀects, which is often an

important factor when studying biomolecules, are taken into account.

Van der Waals interaction accounts for the eﬀect of quantum ﬂuctuations on the

dipole moments. Consider two atoms, of distance r apart, where r is much larger than

the size of the atoms. Even though the atoms may not have permanent dipoles, quantum

and thermal ﬂuctuations will cause them to have temporary dipole moments. If one atom

acquires some temporary dipole moment, it induces a polarization on the other atom.

This induced polarization is of the order O(r

−3

) (see for example problem 3 in Chapter

11 of [32]). The energy for the dipole-dipole interaction is of the form

V (r) = −

(4.2.7)

where α

and α

are the polarizability constants for the two atoms.

(4.2.7) is always negative, implying t hat the van der Waals interaction is attractive. At

short distance, it has to be compensated by a repulsive force that prevents the two atoms

from collapsing into each other. The repulsive force comes either from charge-charge

interaction, or Pauli’s exclusion principle which prevents two electrons from occupying

the same quantum state. Unfortunately, there is no analogous derivation for the form of

the repulsive force. As a matter of convenience, it is often mo deled by a r

−12

term. This

gives the Lennard-Jones form of the van der Waals interaction:

V (r) = 4ε







−







(4.2.8)

4.2. MOLECULAR DYNAMICS 165

The two parameters in the Lennard-Jones potential are: ε, which is a scale of interaction

energy, and σ, which is a scale of the length. When multiple species of atoms are

encountered, a mixing rule is often used for the parameter values in the Lennard-Jones

potential, e.g.

(σ

+ σ

), ε

√

(4.2.9)

Empirical potentials for solids

Solids can be classiﬁed into several d iﬀerent types according to the nature of the

bonding between the atoms in the solid. These are: Covalently bonded, ionic, metallic,

molecular and hydrogen-bonded [4]. The last type occurs mostly in biological systems.

We will not discuss it here.

1. Molecular crystals. These occur mainly for inert gas elements in the last column of

the periodic table, for which the atoms interact mainly through the van der Waals

eﬀect, modeled by the Lennard-Jones potential.

2. Metallic systems. An important ingredient for modeling met allic systems is the

charge density, which is often deﬁned through

∗

) (4.2.10)

where ρ

∗

is a ﬁtting function. The most popular inter-atomic potential for metals

is the embedded atom model (EAM) [11]

V =

i6=j

Φ(r

) +

F (ρ

) (4.2.11)

Here the functions Φ, F are also ﬁtting functions. One special case is the Finnis-

Sinclair potential where F (ρ) = A

√

ρ.

3. Covalently-bonded systems. These are somewhat similar to the molecules discussed

earlier. Examples include the diamond structure of silicon and the graphene struc-

ture of carbon . Two most p opular potentials are the Stillinger-Weber potential and

Tersoﬀ potential. The Stillinger-Weber potential takes the form:

V =

f(r

) +

(h(r

, r

, θ

jik

) + permutations) (4.2.12)

166 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

where

h(r, r

′

, θ) = λe

r−a

′

−a



cos θ +



(4.2.13)

This potential is similar to the potentials we discussed above for molecules. The ﬁrst

term is the pair-wise interaction. The second three-body term is the contribution

of the bond angles.

Tersoﬀ potential, on the other hand, is quite diﬀerent. It takes the form:

V =

i6=j

−αr

− b

−βr

(4.2.14)

The interesting parameter here is b

, the local bond-order, which depends on the

coordination of the atoms and the bond strength [35].

4. Ionic crystals, such as rock salt. In these systems, the atoms interact mainly

through the ion-ion electrostatic interaction and the short range repulsion. Born

proposed the following inter-atomic potential [4]

V =

i6=j

4πε

(4.2.15)

where q

is the partial charge of the i-th ion, A and n are empirical parameters.

The diﬃculty with this potential and with modeling ionic systems in general is the

long range nature of the electrostatic interaction.

Empirical potentials for ﬂuids

Fluids are rarely modeled at the level of atomic detail. They are often coarse-grained.

The choice of the empirical potential depend s on how we coarse-grain the system. If

the coarse-grained particles are treated as spherically-shaped point particles, then the

Lennard-Jones potential is often used . If they are treated as ellipsoids or rods, then one

may use some anisotropic forms of the Lennard-Jones potential such as the Gay-Berne

potential [1]. If the coarse-grained particles are chains, i.e. beads connected by springs,

then bead-spring models are used. S uch a model typically contains two components:

1. The beads interact through a Lennard-Jones potential.

2. Neighboring beads on the same chain interact via a spring force.

An example of this type is treated in Chapter 6. Coarse-grained models of ﬂuids are

discussed more systematically in [40].