Weinan E. Principles of Multiscale Modeling

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

5.2. EXTENSIONS OF THE CAUCHY-BORN RULE 227

neighboring atoms via Taylor expansion:







j,1

j,2

j,3













j,1

j,2













j,1

, x

j,2

)

j,1

, x

j,2

)

j,1

, x

j,2

)







≈







˜y

j,1

˜y

j,2

˜y

j,3







= x

+ u(x

i,1

, x

i,2

) +

∂u(x

i,1

, x

i,2

)

∂x

j,1

− x

i,1

)

∂u(x

i,1

, x

i,2

)

∂x

j,2

− x

i,2

) +

∂

u(x

i,1

, x

i,2

)

∂x

j,1

− x

i,1

)

∂

u(x

i,1

, x

i,2

)

∂x

j,1

− x

i,1

)(x

j,2

− x

i,2

)

∂

u(x

i,1

, x

i,2

)

∂x

j,2

− x

i,2

)

(5.2.19)

where u = (u

, u

)

is the displacement ﬁeld, x

= (x

i,1

, x

i,2

) is the reference position

of the i-th atom. Substituting this into the atomistic model, we obtain the stored energy

density at x

LCB

(∇u, ∇

u) =

j,k

V (

) (5.2.20)

where D

is the area of a unit cell,

= (x

i,1

, x

i,2

, 0)

Note that this procedure can be used for both simple and complex lattices.

As an example, let us consider the grapheen sheet which is a two dimensional complex

lattice obtained from the union of two triangular lattices. In equilibrium, one set of basis

vectors for the triangular lattice is given by:

= a(1, 0), A

= a

√

(5.2.21)

where a is the lattice constant, the shift vector is:

p = a

√

. (5.2.22)

This example is treated in [1] and [46] using the Tersoﬀ potential as the atomistic model.

Figure 5.2.3 shows the results of a twisted nanotube, m odeled using the atomistic model

as well as the continuum model obtained using the exponential Cauchy-Born rule. Figure

5.2.3 shows the resu lts of the nanotube under compression, modeled using the atomistic

228 CHAPTER 5. EXAMPLES OF MULTI-PHYSICS MODELS

model and th e continuum model obtained using the local Cauchy-Born rule. Very good

agreement is found in both cases. This demonstrates that mechanical deformation of

nanotubes can be very well described by continuum theories in this regime. The value of

such a viewpoint has also been demonstrated using much simpliﬁed continuum models

[45].

220

300

Figure 5.4: Comparison of the results of the MD simulation and the continuum model for

twisted nano-tube using the exponential Cauchy-Born rule. Dots represent the results of

the MD simulation. The grey background is computed from the exponential Cauchy-Born

rule ([1], courtesy of Ted Belytschko).

5.3 The moving contact line problem

When water or other liquid drops are placed on a solid surface, a contact line forms

at the intersection between the liquid, the air and the solid. More generally, the same

situation arises when considering the dynamics of immiscible ﬂuids in a channel [12, 14],

as in the case of secondary oil recovery when oil is drained out of the pores y some other

liquid. Like the problem of crack propagation, contact line dynamics has often been used

as an example to illustrate the need for formulating coupled continuum-atomistic models,

since it has been realized for some time that classical continuum models are inadequate

near the contact line and atomistic models should give a more accurate representation [20,

33]. While developing coupled atomistic-continuum models is certainly a very important

avenue of research, we will illustrate in this section that multiscale ideas can also be used

in a diﬀerent way, to help formulating better continuum models.

5.3. THE MOVING CONTACT LINE PROBLEM 229

−10 −5 0 5 10

x (

z (

−10 −5 0 5 10

x (

z (

(a) (b)

Figure 5.5: Comparison of the results of the MD simulation (left) and the continuum

model (r ight) for compressed nano-tube using the local Cauchy-Born rule (from Yang

and E [46]).

5.3.1 Classical continuum theory

The setup is shown in Figure 5.3. The three phases are separated by three interfaces

, Γ

and Γ

respectively. We denote the su rface tension coeﬃcient of Γ

by γ

, i = 1, 2, 3.

The equilibrium problem and the static conﬁguration of the contact line was under-

stood a long time ago in the work of Laplace, Young and Gauss [11]. At equilibrium, the

total energy of the system is the sum of the interfacial energies

E =

i=1,2,3

|Γ

| (5.3.1)

We will only consider the case of partial wetting, i.e.

|γ

− γ

| < γ

(5.3.2)

In this case, the Euler-Lagrange equation of the functional (5.3.1) is given by:

− γ

+ γ

cos θ

= 0. (5.3.3)

230 CHAPTER 5. EXAMPLES OF MULTI-PHYSICS MODELS

Γx

Ω

Figure 5.6: A two-phase (in Ω

and Ω

respectively) ﬂuid on the substrate. x

is the

contact line (here a p oint) and Γ

(i = 1, 2, 3) are the interfaces. θ

is the contact angle

(courtesy of Weiqing Ren).

This is Young’s relation for the case of partial wetting, and it deﬁnes the static contact

angle θ

Now let us turn to the dynamic problem. To be speciﬁc, let us consider the situation

of two immiscible ﬂuids in a channel separated by an interface. To model the dynamics

of such a system, the most obvious approach is to use the Navier-Stokes equation in the

ﬂuid phases:

(∂

u + u · ∇u) = −∇p + µ

∆u, in Ω

∇ · u = 0, (5.3.4)

together with the no-slip boundary condition on the solid surface:

u = u

. (5.3.5)

Here ρ

and µ

, i = 1, 2 are the density and viscosity of the two ﬂuid phases respectively,

is the velocity of the wall. At the ﬂuid-ﬂuid interface, we use the standard Laplace-

Young interfacial condition:

κ = n · [τ ] · n, t · [τ ] · n = 0 (5.3.6)

where t and n denote the un it tangent and normal vectors to the interface respectively, τ

denotes the total stress (viscous str ess and pressure) in the ﬂuids, κ is the mean curvature

of the ﬂuid-ﬂuid interface, [·] denotes the jump of the quantity inside the bracket across

the interface.

5.3. THE MOVING CONTACT LINE PROBLEM 231

It was discovered long ago that in the case of partial wetting, this model predicts a

singularity at the contact line [21], known as the Huh-Scriven singularity. This is most

clearly seen in two dimension. Let (r, φ) be the polar coordinates centered at the contact

line. Then a simple Taylor expansion at the contact line gives, for r ≪ 1,

ψ ∼ r ((Cφ + D) cos φ + (Eφ + F ) sin φ) (5.3.7)

for the streamfunction ψ, where C, D, E, F are some constants. The velocity is given in

terms of the stream function via u = (−∂

ψ, ∂

ψ). Therefore we obtain

∇u ∼

|∇u|

dx = +∞ (5.3.8)

For this reason, questions have been raised on the validity of the conventional hydrody-

namics formalism, the Navier-Stokes equation with no-slip boundary condition, for this

problem. Many proposals have been made to remove this singularity (for a review, see

[33]). While they have all succeeded in doing so, it is not clear whether any of them has

a solid physical foundation.

Detailed molecular dynamics studies have revealed the fact that near the contact line,

there is a slip region in which the no-slip boundary condition is violated [24, 38]. The

existence of such a slip region is not surprising. In fact in a Couette ﬂow setting (i.e.

channel ﬂow driven by walls moving in the opposite direction), if the system reaches a

steady state, then the contact line has to exp erience complete slip, i.e. the relative slip

velocity between the ﬂuids and the solid surface should be equal to the velocity of the

solid surface. What is more important is that these molecular dynamics studies have

opened up the door to a detailed analysis of what is going on near th e contact line at the

atomistic level.

Even though molecular dynamics is able to oﬀer much needed details about the con-

tact line problem, it is n ot very attractive as a theoretical tool, compared with the

continuum models. It is much harder for us to think about the behavior of the individual

molecules instead of the continuum ﬁelds. Therefore, from a theoretical viewpoint, it is

still of interest to develop better continuum models. We will demonstrate that multiscale

ideas are indeed quite useful in this respect.

232 CHAPTER 5. EXAMPLES OF MULTI-PHYSICS MODELS

5.3.2 Improved continuum models

The strategy that we will discuss consists of two steps. The ﬁrst step is to ﬁnd

the simplest form of the constitutive relations including the boundary conditions, from

thermodynamic considerations. The second step is to measure the detailed functional

dependence of the constitutive relation using molecular dynamics. Such a strategy was

pursued in [35].

There are three constitutive relations that we need:

1. The constitutive relation in the bulk, which describes the response of the ﬂuid in

the bulk. It has been argued that even for simple ﬂuids, the linear constitutive

relation that is usually quite accurate for bulk ﬂuids, may not be accurate enough

near the contact line, due to the presen ce of large velocity gradients there [38, 40].

The careful molecular dynamics studies carried out by Qian et al. do not support

this argument [31]. Instead, it is found that, at least for simple ﬂuids, the linear

constitutive relation holds quite well right around the contact line. In other words,

even though the velocity gradient is quite large, it is not large enough to induce

signiﬁcant nonlinear eﬀects in the bulk response. We will assume for simplicity that

linear constitutive relations are suﬃcient in the bulk.

2. The constitutive relation at the ﬂuid-solid interface. It has been observed that the

standard no-slip boundary condition no longer holds in a region near the contact

line [24, 25, 31, 32, 33, 38, 40]. This region is the slip region. Linear friction law at

the surface gives rise to the Navier boundary condition. However, nonlinear eﬀects

might be important there.

3. The constitutive relation at the contact line, which determines the relationship

between the contact angle and the contact line velocity. One obvious possibility

is to set the dynamic contact angle to be the same as the equilibrium contact

angle. However, this eﬀectively suppresses energy dissipation at the contact line.

Dissipative processes at the contact line can be quite important [5]. Therefore in

general, one has t o allow the deviation of the dynamic contact from its equilibrium

value.

To derive the needed constitutive relation, we follow the principles discussed in Chap-

ter 4 and ask the following question: What are the simplest forms of boundary conditions

5.3. THE MOVING CONTACT LINE PROBLEM 233

that are consistent with the second law of thermodynamics?

For simplicity, we will focus on the case of two-dimensional ﬂows. The total (free)-

energy of the system is given by

E =

i=1,2

Ω

|u|

dx +

i=1,2,3

|Γ

|. (5.3.9)

A straightforward calculation gives

= −

i=1,2

Ω

|∇u|

dx +

i=1,2

u · τ · n dS + γ

(cos θ

− cos θ

) u

, (5.3.10)

where θ

is the dynamic contact angle, θ

is the static contact angle deﬁned in (5.3.3),

is the velocity of the contact line, n is the unit normal vector of the solid surface.

Let u

and τ

be the slip velocity and the shear stress respectively,

= (u − u

) · t, τ

= t · τ · n. (5.3.11)

Let

= γ

(cos θ

− cos θ

) . (5.3.12)

is the so-called unbalanced or uncompensated Young’s stress [5, 12, 32, 33]. It represents

the force on the contact line that arises from the deviation of the ﬂuid interface from its

equilibrium conﬁguration. With these notations, we can rewrite (5.3.10) as

= −

i=1,2

Ω

|∇u|

dx +

i=1,2

dS + u

(5.3.13)

According to t he second law of thermodynamics, the rate of energy dissipation dE/dt

is necessarily non-positive for arbitrary realizations of u. It follows easily that each of

the three contributions has to be non-positive. We now examine the implication of this

basic requirement and the constraints it places on the form of the constitutive relations.

First, we look at the ﬂuid-solid interfaces Γ

and Γ

. There we must have

u · n = 0 (5.3.14)

This implies that the shear stress should depend only on u

. Assumin g that this depen-

dence is local, then τ

is a function of the local values of u

and its derivatives. To the

lowest order, we simply have:

= f(u

) (5.3.15)

234 CHAPTER 5. EXAMPLES OF MULTI-PHYSICS MODELS

Since

dS 6 0 (i = 1, 2) for an arbitrary realization of u, the fun ction f must

satisfy wf(w) ≤ 0

Next, we turn to the contact line. From thermodynamic consideration, τ

is neces-

sarily a functional of the contact line velocity u

. The simplest assumption is that it is

a local function of u

= f

) (5.3.16)

The function f

has to satisfy uf

(u) ≤ 0.

The Navier-Stokes equation together with the boundary conditions (5.3.6), (5.3.14),

(5.3.15) and (5.3.16) constitute our continuum model for the moving contact line problem.

The constitutive relations in (5.3.15) and (5.3.16) are not yet fully speciﬁed. The exact

functional dependence of f and f

has to be obtained through other means. For this

continuum model, the energy dissipation relation has a very clean form:

= −

i=1,2

Ω

|∇u|

dx +

i=1,2

f(u

) dS + u

) (5.3.17)

The three terms at the right hand side represent the energy dissipated in the bulk ﬂuids,

the solid surface and the contact line, respectively.

So far, we have followed the standard procedure in generalized thermodynamics or

non-equilibrium thermodynamics [13, 30]. At this point, we have two options. One is to

continue following the standard procedure and write down linear relations between the

generalized ﬂuxes (here the u

and u

) and the generalized forces (here the τ

and τ

The other is to use microscopic models to measure f and f

directly.

If we continue to follow the standard procedure in non-equilibrium thermodynamics,

then the next step is to propose linear constitutive relations of the form (see Section 4.1):

= −β

, i = 1, 2 (5.3.18)

and

(cos θ

− cos θ

) = −β

(5.3.19)

, i = 1, 2 and β

are the friction coeﬃcient at the solid surface for the two ﬂuid phases

and the contact line, respectively. Equation (5.3.18) is the well-known Navier boundary

condition. Putting things together, we obtain the following continu um model for the

moving contact line problem:

5.3. THE MOVING CONTACT LINE PROBLEM 235

1. In the bulk, we have







(∂

u + (u · ∇)u) = −∇p + µ

∇

u + f

∇ · u = 0

(5.3.20)

2. The ﬂuid interface Γ is advected by the velocity ﬁeld, and

n · [τ] n = γ

t · [τ ] n = 0

(5.3.21)

where τ = −pI + µ(∇u + ∇u

) is the stress tensor; κ is the mean curvature of the

interface.

3. On the solid surface, we have

= t · τn, u · n = 0 (5.3.22)

4. At the contact line, we have

= −γ

(cos θ

− cos θ

) . (5.3.23)

This is the model put forward by Ren and E [33].

Some remarks about the physical signiﬁcance of the parameter β

are in order:

1. It is a three-phase friction coeﬃcient.

2. It has th e dimension of viscosity. This should be contrasted to β

which has

the dimension of viscosity/length and this length has to be a microscopic length.

Conequently, the Navier slip boundary condition is only signiﬁcant in a slip region

of microscopic size. In comparison, β

should be of macroscopic signiﬁcance.

3. It also arises in the sharp interface limit of the diﬀusive interface model proposed

in [23, 31, 32], as was demonstrated in [44].

One can write (5.3.23) as

cos θ

− cos θ

= −(β

/γ

(5.3.24)

236 CHAPTER 5. EXAMPLES OF MULTI-PHYSICS MODELS

This can be viewed as a kinetic correction to the condition

= θ

(5.3.25)

(5.3.25) states that the dynamic contact angle is equal to the equilibrium contact angle.

This is the analog of the Gibbs-Thomson condition in the theory of crystal growth [19]. In

this regard, (5.3.24) is the analog of the corrected Gib bs-Thomson relation with kinetic

eﬀects taken into account [19].

It is interesting to note that if the condition (5.3.25) is used together with the Navier

boundary condition, then the Huh-Scriven singularity is removed, b ut there is still a weak

singularity [37]. Naturally one expects that if we use (5.3.24) instead of (5.3.25), then

this weak singularity is also removed. This is still an open question.

Since any function is approximately linear in a small enough interval, these linear

constitutive relations will be adequate in appropriate regimes. What is important is how

big these regimes are and what happens outside of them. To answer these questions, we

resort to microscopic models.

5.3.3 Measuring the boundary conditions using molecular dy-

namics

Thermodynamic consideration tells us what the boundary conditions should depend

on, but not the detailed functional dependence. This latter information should be ob-

tained from experiments, or, as we explain now, numerical experiments using microscopic

simulations, such as molecular dynamics.

The results of the previous section tell us that f and f

depend on the slip velocity.

Therefore, we need to set up an MD system such that, at the statistical steady state, it

has a constant slip velocity. This can be done using a MD system in the Couette ﬂow

geometry, shown in Figure 5.7, where the ﬂuid particles are conﬁned between two solid

walls, and the walls move in opposite directions at a constant speed U [24, 31, 33, 38, 39].

A standard way of modeling immisible ﬂuids is to use a modiﬁed form of the Len nard-

Jones (LJ) potential:

(r) = 4ε







− ζ







(5.3.26)

and use ζ = 1 for particles of the same species and ζ = −1 for particles of diﬀerent species.

Diﬀerent values of ε can be used to model the ﬂuid-ﬂuid and ﬂuid-solid interaction