Weinan E. Principles of Multiscale Modeling

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

4.1. CONTINUUM MECHANICS 147

One main pu rpose of multiscale, multi-physics modeling is to establish a microscopic

foundation for the constitutive relations. This is not a new idea. There is a long history

of work along this line (see for example [4, 10]). Back in the 1950’s, Tsien initiated the

subject of Physical Mechanics with the objective of providing the microscopic foundation

for continuum mechanics [39]. The practical interest was to calculate material properties

and transport co eﬃcients using principles of statistical and quantum mechanics, particu-

larly at extreme conditions. For various reasons, the subject did not become very popular

after its initial success. The interest on multiscale, multi-physics modeling has given it a

new impetus, and the right time seems to have come to try reviving the subject.

4.1.1 Stress and strain in solids

Given a solid material, denote by Ω the region it occupies in the absence of any

external loading. Our objective is to ﬁnd its new conﬁguration after some external load

is applied to the material. We will refer to the conﬁguration of the material before and

after the loading is applied as the reference and deformed conﬁgurations respectively.

Let x be the position of a material point in the reference conﬁguration, and y be

the position of the same material point in the deformed conﬁguration. Obviously y is a

function of x : y = y(x). u(x) = y(x) − x is called the displacement ﬁeld,

F = ∇u (4.1.1)

is called the deformation gradient tensor.

Let δx be a small change of x in the reference con ﬁguration and δy be the corre-

sponding change in the deformed conﬁguration. To leading order, we have

δy = (I + F)δx, (4.1.2)

where I is the identity tensor,

kδyk

− kδxk

= (δx)



(I + F)

(I + F) − I



δx (4.1.3)

= 2(δx)

Eδx.

The strain tensor E is deﬁned as

E =



(I + F)

(I + F) − I



(4.1.4)

+ F) +

148 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

In component form, we have



∂u

∂x

∂u

∂x



∂u

∂x

∂u

∂x

. (4.1.5)

The tensor E = (e

) characterizes the geometric change as a result of the deformation.

It is most relevant for solids since for ﬂuids, the deformation is usually quite large and

therefore it is more convenient to discuss the rate of strain, which is deﬁned as the

symmetric part of the gradient of the velocity ﬁeld:

D =

(∇v + (∇v)

) (4.1.6)

where v is the velocity ﬁeld.

Imagine an inﬁnitesimal surface with area ∆S at a point inside the material, with

normal vector n. Let g∆S be the force acting on one side of the surface (the side opposite

to n) by the other side (i.e. the side in the direction of n). It is easy to see that there is

a linear relation between n and g, which can be expressed as:

g = σn

σ is called the stre ss tensor.

There are two ways of deﬁning the normal vector and the area of the inﬁnitesimal

surface, depending on whether we use the Lagrangian or t he Eulerian coordinate system.

When the Eulerian coordinate system is used, the resulting stress tensor is called the

Cauchy stress tensor. When the Lagrangian coordinate system is used, the resulting

stress tensor is called the Piola-Kirkhoﬀ stress tensor. This distinction is very important

for solids, since the Lagrangian coordinate system is typically more convenient for solids,

but the Cauchy stress is easier to use and to think about. Fortunately, there is a simple

relation between the two notions.

Lemma 5. Let y = y(x) be a nondegenerate map on R

, i.e. J(x) = det(∇

y(x)) 6= 0.

Denote by dS

an in ﬁni tesimal surface element at x, and dS its image under the mapping

y(·). Let n

and n be the normal to dS

and dS respectively. Then

ndS = J(x)(∇

y(x))

−T

(4.1.7)

For a proof of this lemma, see [29]. We will use σ to d enote the Cauchy stress tensor,

and σ

the Piola-Kirkhoﬀ stress tensor.

4.1. CONTINUUM MECHANICS 149

Let S

be an arbitrary closed surface in the reference conﬁguration and let S be its

image under the mapping y(·). From this lemma, we have:

σ(y)n(y)dS(y) =

J(x)σ(y(x))(∇

y(x))

−T

(x)dS

(x) (4.1.8)

Hence

(x) = J(x)σ(y(x))(∇

y(x))

−T

(4.1.9)

Note that since σ

= σndS, the force vectors σ(y)n(y) and σ

(x) are colinear.

Their magnitude can be diﬀerent as a result of the diﬀerence between dS

and dS. They

are the same physical stress, i.e. the stress (force per area) in the deformed body at

y = y(x). Their diﬀerence is the result of the diﬀerent ways that area is measured.

In ﬂuids, due to the large deformation, one usually uses the Cauchy stress tensor.

4.1.2 Variational principles in elasticity theory

The static conﬁguration of a deformed elastic body is often characterized by a varia-

tional principle of the type:

min

Ω

(W ({u}) − f · u) dx, (4.1.10)

where f is the external force density, and W ({u}) is called the stored energy density,

or elastic energy density. It is the potential energy or free energy of the material. The

dependence on {u} means W can depend on u, ∇u, ···. This function plays the role of

the constitutive relation for this type of models.

One of the simplest forms of the constitutive relation is obtained when W only depends

on E, i.e.

W ({u}) = W

(E). (4.1.11)

This is often enough for the analysis of the elastic deformation of bulk materials.

In the absence of any additional information, we can always study ﬁrst the case when

the deformation is small, small enough so that we can further assume that the material

response is linear, namely, W

(E) is a quadratic function of the strain tensor

(E) =

i,j,k,l

ijkl

, (4.1.12)

150 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

where {e

} is deﬁned by (4.1.5). The constant term is neglected and the linear term

should vanish since the undeformed state is at equilibrium (therefore has no residual

stress). The coeﬃcients {C

αβγδ

} are called the elastic constants. For a three dimensional

problem, C has 3

= 81 components. We can always require C to be symmetric, i.e.

ijkl

= C

klij

. (4.1.13)

In addition, since E is a symmetric tensor, we have

ijkl

= C

ijlk

. (4.1.14)

These relations reduce the number of independent elastic constants to 21 [27].

In most cases, the material under consideration has additional symmetry which allows

us to reduce the number of independent parameters even further. Here we will discuss

two cases, isotropic materials and materials with cubic symmetry.

If a m aterial has cubic symmetry, its elastic properties are invariant under the ex-

change of any two coordinates. Let T be such a transformation:

• T (1) = 2, T (2) = 1, T (3) = 3

• T (1) = 3, T (2) = 2, T (3) = 1

• T (1) = 1, T (2) = 3, T (3) = 2

then

T (i)T (j)T (k)T (l)

= C

ijkl

(4.1.15)

In addition, W also has inversion symmetry: Let I be any of the following transforma-

tions:

• I(x) = −x, I(y) = y, I(z) = z

• I(x) = x, I(y) = −y, I(z) = z

• I(x) = x, I(y) = y, I(z) = −z

then

I(i)I(j)







, if i = j;

−e

, if i 6= j.

(4.1.16)

4.1. CONTINUUM MECHANICS 151

This implies that all elastic constants with any subscript appearing an odd number of

times must vanish. For instance, C

1112

has to be 0.

These symmetry considerations allow us to reduce the number of independent elastic

constants for a material with cubic symmetry to only 3. They are:

1111

= C

2222

= C

3333

1212

= C

2323

= C

3131

1122

= C

1133

= C

2233

Using this information, we can write W

(E) as

(E) =

1111



+ e



+ C

1122

+ e

) (4.1.17)

+ 2C

1212



+ e



In solid state physics, we often adopt the Voigt notation:

= e

, e

= e

, e

= e

, (4.1.18)

= 2e

, e

= 2e

, e

= 2e

= C

1111

, C

= C

1122

, C

= C

1212

For isotropic materials, we have the additional property that W is invariant under a

degree rotation around any coordinate axis. It is easy to verify that this implies the

following relation

= C

+ 2C

(4.1.19)

This reduces the number of independent elastic constants to 2:

λ = C

, µ = C

. (4.1.20)

They are called the Lam´e constants. Now W

(E) can be written as

(E) =

+ µ

. (4.1.21)

Let us now consider brieﬂy nonlinear isotropic materials. We will limit ourselves to

the case of incompressible materials, i.e. materials for which volume is invariant during

152 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

the deformation. Since W is invariant under rotation, it is only a function of the principal

values of the matrix C = (I + F)

(I + F) (recall th at F = ∇u):

W = W (λ

, λ

) (4.1.22)

where λ

, λ

are the eigenvalues of C. Examples of the empirical mo dels of W include

1. The Neo-Hookean models

W = a(λ

+ λ

− 3) (4.1.23)

where a is a constant.

2. The Mooney-Rivlin models

W = a(λ

+ λ

− 3) + b(λ

+ λ

− 3) (4.1.24)

where a and b are constants.

More examples can be found in [29].

So far we have considered the case when W depends only on ∇u. This is no longer

adequate, for example, for:

• small systems, for which the dependence on the strain gradient should be included;

• plates, sheets and rods, for which dependence on the curvature should be included.

Including the dependence on the strain gradient in W increases the complexity of the

models substantially. But in principle, the same kind of consideration still applies.

4.1.3 Conservation laws

We now turn to dynamics. The ﬁrst step toward studying dynamics is to write down

the three universal conservation laws of mass, momentum and energy.

Conservation of mass. We will ﬁrst discuss mass conservation in Lagrangian coor-

dinates. Let D

be an arbitrary subdomain in the reference conﬁguration, and D be its

image in th e deformed conﬁguration. The mass inside D must be the same as the mass

4.1. CONTINUUM MECHANICS 153

inside D

. If we denote by ρ

the density function in the reference conﬁguration and ρ

the density function in the deformed conﬁguration, then we have

(x)dx =

ρ(y)dy (4.1.25)

Since

dy = J(x)dx (4.1.26)

we obtain

(x) = ρ(y(x))J(x) (4.1.27)

In a time-dependent setting, this can be written as:

(x) = ρ(y(x, t))J(x, t) (4.1.28)

We now turn t o the Eulerian form of mass conservation which is more convenient for

ﬂuids. Let D be an arbitrary subdomain in the deformed conﬁguration, ﬁxed in time.

The instantaneous change of the mass inside D is equal to the mass ﬂux through the

boundary of D:

ρ(y, t)dy = −

∂D

· ndS(y) (4.1.29)

Here n is the unit outward normal of ∂D, J

is the mass ﬂux. It is easy to see that

= ρv (4.1.30)

where v is the velocity ﬁeld. Using the divergence theorem, we get

(∂

ρ(y, t) + ∇ · (ρv))dy = 0 (4.1.31)

Since this is true for any domain D, we have

∂

ρ + ∇

· (ρv) = 0 (4.1.32)

For the next two conservation laws, we will ﬁrst discuss the Eulerian formulation.

Their Lagrangian formulation will be discussed in the next sub-section.

Conservation of momentum. The rate of change of the total momentum in D is

the result of: (1) materials moving in and out of D due to convection and (2) the forces

acting on D. In general there are two kinds of forces acting on D: body force such as

154 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

gravity and surface force acting on the boundary of D, i.e. the stress. We will neglect

body force since it can always be added back later. Putting these together, we have

ρ(y, t)v(y, t)dy = −

∂D

ρv(v · n)dS(y) +

∂D

σ · ndS(y) (4.1.33)

In a diﬀerential form, this gives

∂

(ρv) + ∇

· (ρv ⊗ v − σ) = 0 (4.1.34)

We can also write this as

∂

(ρv) + ∇

· J

= 0 (4.1.35)

where J

is the total momentum ﬂux given by

= ρv ⊗ v − σ (4.1.36)

Conservation of angular momentum. Similarly, the rate of change of the total

angular momentum in D is given by:

ρ(y, t)y × v(y, t)dy = −

∂D

ρy × v(v · n)dS(y) +

∂D

y × σ · ndS(y) (4.1.37)

This states that the rate of change of the angular momentum inside D is due to angular

momentum ﬂux through ∂D and the torque generated by the stress on ∂D. The main

consequence of this is that

σ = σ

(4.1.38)

i.e. th e Cauchy stress tensor is symmetric. For a proof, see [29].

Conservation of energy. In the same way, we can write down the energy conser-

vation law. The total energy density E is the sum of the kinetic energy density and the

potential energy density:

E =

ρ|v|

+ ρe (4.1.39)

where e is the internal energy per unit mass. Given a volume D, the rate of change of

energy inside D is the result of energy been convected in and out of D, the work done

by the external force acting on D and the heat ﬂux:

Edy = −

∂D

E(v · n)dS(y) +

∂D

v · (σ · n)dS(y) +

∂D

q · ndS(y) (4.1.40)

4.1. CONTINUUM MECHANICS 155

Here q is the heat ﬂux. In a diﬀerential form, this gives:

∂

E + ∇

· (vE − σv + q) = 0 (4.1.41)

∂

E + ∇

· J

= 0 (4.1.42)

where

= vE − σv + q (4.1.43)

4.1.4 Dynamic theory of solids and thermoelasticity

Dynamic models of continuous media are derived from the following principles:

1. Conservation laws of mass, momentum and energy.

2. Second law of thermodynamics, the Clausius-Duhem inequality.

3. Linearization and symmetry. In the absence of additional information, a linear

constitutive relation is assumed and symmetry is used to reduce the number of

independent parameters.

In this and the next sub-section, we will demonstrate how such a program can be carried

out for simple ﬂuids and solids.

As we remarked earlier, for solids, it is more convenient to use the Lagrangian coor-

dinate system. Therefore our ﬁrst task is to express the conservation laws in Lagrangian

coordinates.

Let us begin with mass conservation. As was derived earlier, we have

ρ(y(x, t), t) =

det(I + F)

(x). (4.1.44)

Similar to mass conservation, momentum conservation requires for any subdomain in

the reference conﬁguration D

(x)∂

y(x, t)dx =

∂D

σ(y, t)

n(y, t)dS(y) =

∂D

(x, t)

(x, t)dS(x).

(4.1.45)

Here D

is the image of D

under the mapping x → y(x, t) at time t,

n(y, t) is the unit

outward normal vector of ∂D

at y,

(x) is the normal vector of ∂D

at x, σ and σ

are

156 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

respectively the Cauchy and the Piola-Kirkhoﬀ stress tensor. We have made use of the

relation (4.1.7). Applying Green’s formula to the right hand side of (4.1.45), we get

(x)∂

y(x, t) = ∇

· σ

(x, t). (4.1.46)

Similarly, energy conservation implies



(x)|∂

y(x, t)|

+ ρ

(x)e(x, t)



dx (4.1.47)

∂D

∂

y(x, t) · σ

(x, t)dS(x) −

∂D

q(x, t) ·

(x, t)dS(x).

Applying Green’s formula again, we get

∂

∂t



(x)|∂

y(x, t)|

+ ρ

(x)e(x, t)



= ∇

·(∂

y(x, t) · σ

(x, t))−∇

·q(x, t), (4.1.48)

Taking (4.1.46) into consideration, this equation can be rewritten as

∂

e = (∇

∂

y) : σ

− ∇

· q (4.1.49)

Here we have used the notation A : B to denote the scalar product of two matrices:

A : B =

i,j

, where A = (a

), B = (b

Additional information comes from the second law of thermodynamics. Denote by

s(x, t) the entropy per unit mass, T (x, t) the temperature of the system, the second law

states that the rate at which entropy increases can not be slower than the contribution

due to the exchange of heat with the environment:

(x)s(x, t)dx ≥ −

∂D

q(x, t)

T (x, t)

(x)dS(x). (4.1.50)

In a diﬀerential form, this gives

∂

s ≥ −∇





(4.1.51)

This is referred to as the Clausius-Duhem inequality.

Deﬁne the Helmholtz free energy per unit mass as

h(x, t) = e(x, t) − T (x, t)s(x, t). (4.1.52)