Continuum Mechanics

(1)

Continuum Mechanics

I-Liang Chern

Department of Mathematics

National Taiwan University

(2)

(3)

Chapter 1 Fluid Mechanics

1.1 Thermodynamics

The equations of gas dynamics are derived based on conservation of mass, momentum and energy. Before we derive these equations, let us review some thermodynamics. First, the basic thermo variables are pressure (p), specific volume (τ ), called state variables. The internal energy (e) is a function of p and τ . Such a relation is called a constitutive equation. The basic assumption are

∂e ∂p τ > 0, ∂e ∂τ p > 0

Sometimes, it is convinient to express p as a function of (τ, e).

In an adiabetic process (no heat enters or losses), the first law of thermodynamics (conservation of energy) reads

de + pdτ = 0. (1.1) This is called a Pfaffian equation mathematically. A function σ(e, τ ) is called an integral of (1.1) if there exists a function µ(e, τ ) such that

dσ = µ · (de + pdτ ).

Thus, σ = constant represents a specific adiabetic process. For Pfaffian equation with only two independent variables, one can always find its integral. First, one can derive equation for µ: from

σe = µ and στ = µp

and using σeτ = στ e, we obtain the equation for µ:

µτ = (µp)e.

This is a linear first-order equation for µ. It can be solved by the method of charac-teristics in the region τ > 0 and e > 0. The solutions of µ and σ are not unique. If σ is a solution, so does ¯σ with d¯σ = ν(σ)dσ for any function ν(σ). We can choose µ such that if two systems are in thermo-equilibrium, then they have the same value µ.

(7)

In other words, µ is only a function of emperical temperature. We shall denote it by 1/T . Such T is called the absolute temperature. The corresponding σ is called the physical entropy S. The relation dσ = µ(de + pdτ ) is re-expressed as

de = T dS − pdτ. (1.2) For ideal gas, which satisfies the laws of Boyle and Gay-Lussac:

pτ = RT, (1.3)

where R is the universal gas constant. From this and (1.2), treating S and τ as independent variables, one obtains

ReS(S, τ ) + τ eτ(S, τ ) = 0.

We can solve this linear first-order equation by the method of characteristics. We rewrite this equation as a directional differentiation:

R ∂ ∂S + τ ∂ ∂τ e = 0.

This means that e is constant along the characteristic curves

Rdτ dS = τ. These characteristics can be integrated as

τ e−S/R = φ.

Here φ is a positive constant. The energy e(τ, S) is constant when τ e−S/Ris a constant. That is, e = h(φ) for some function h. We notice that h0 < 0 because p = −(∂e_∂τ)S =

−e−S/R_h0_{(τ H) > 0. From T = (}∂e

∂S)τ = − 1 Rh

0_{(φ) · φ, we see that T is a function of φ.}

In most cases, T is a decreasing function of φ. We shall make this as an assumption. With this, we can invert the relation between T and φ and treat φ as a decreasing function of T . Thus, we can also view e as a function of T , say e(T ), and e(T ) is now an increasing function. Now, we have five thermo variables p, τ, e, S, T , and three relations:

pτ = RT e = e(T ) de = T dS − pdτ

Hence, we can choose two of as independent thermo variables and treat the rest three as dependent variables.

For instance, e is a linear function of T , i.e. e = cvT , where cv is a constant called

specfic heat at constant volume. Such a gas is called polytropic gas. We can obtain

pτ = RT and e = cvT =

pτ

(8)

or in terms of entropy, p = A(S)τ−γ T = A(S) R τ −γ+1 e = cvA(S) R τ −γ+1 where A(S) = (γ − 1) exp((S − S0)/cv) γ = 1 + R/cv

If we define dQ = T dS, it is easy to see that cv and cpare the specific heat at constant

volume and constant pressure, respectively.

cv = ∂Q ∂T τ = ∂e ∂T τ , cp := ∂Q ∂T p = ((∂e ∂τ)p+ p)/( ∂T ∂τ)p = ∂e ∂T p + p ∂τ ∂T p

In general, cp > cv. Because cp is the amount of heat added to a system per unit

mass at constant pressure. In order to maintain constant pressure, the volume has to expand (otherwise, pressure will increase), the extra amount of work due to expansion is supplied by the extra amount of heat cp− cv.

1.2 Compressible Euler Equations

Next, we derive the equations of gas dynamics. Let us consider an arbitrary domain Ω ⊂ R3. The mass flux from outside to inside per unit time per unit area dS is −ρv·, where n is the outer normal of Ω. Thus, the conservation of mass can be read as

d dt Z Ω ρ dx = Z ∂Ω [−ρv · n]dS = − Z Ω div (ρ v) dx

This holds for arbitrary Ω, hence we have

ρt+ div(ρ v) = 0. (1.5)

(9)

Now, we derive momentum equation. Let us suppose the only surface force is from pressure (no viscous force). Then the momentum change in Ω is due to (i) the momentum carried in through boundary, (ii) the pressure force exerted on the surface, (iii) the body force. The first term is −ρvv · n, the second term is −pn. Thus, we have d dt Z Ω ρv dx = Z ∂Ω −[ρvv · n + pn] dS + Z F dx = Z Ω div[−ρv ⊗ v − pI] + F dx This yields (ρv)t+ div(ρ v ⊗ v) + ∇p = F (1.6)

Here, the notation ∇ · ρv ⊗ v stands for a vector whoes ith component isP

j∂j(ρvivj).

The energy per unit volume is E = 1₂ρ v2+ ρe. The energy change in Ω per unit time is due to (i) the energy carried in through boundary (ii) the work done by the pressure from boundary, and (iii) the work done by the body force. The first term is −Ev · n. The second term is −pv · n. The third term is F · v. The conservation of energy can be read as d dt Z Ω E dx = Z ∂Ω [−Ev · n − pv · n] dS + Z Ω F · v dx By applying divergence theorem, we obtain the energy equation:

Et+ div[(E + p)v] = ρF · v. (1.7)

In one dimension, the equations are (without body force) ρt+ (ρu)x = 0 (ρu)t+ (ρu2+ p)x = 0 (1 2ρu 2_{+ e)} t+ [( 1 2ρu 2_{+ e + p)u]} x = 0.

Here, the unknowns are two thermo variable ρ and e, and one kinetic variable u. Other thermo variable p is given by the constitutive equation p(ρ, e).

1.3 Invariants

1.3.1 Compressible gas dynamics in Lagrangian coordinate

Given velocity field v(x, t), we can define the particle path x(t, X) to be the solution of the ODE

˙x(t) = v(x, t), x(0, X) = X.

The coordinate X is called the Lagrangian coordinate, or the material coordinate, where x is the orginal coordinate in the physical space. It is also called the Eulerian coordinate. The derivative

d dt :=

∂

(10)

is called the material derivative.

We shall rewite the Euler equation in terms of Lagrangian coordinate.

Rate of change of density First, we have d

dtρ + ρ∇ · v = 0. This means that

∇ · v = − dρ dt ρ = dτ dt τ (1.8)

is the relative rate of change of volume. A fluid is called incompressible if dρ

dt = 0. This is indeed equivalent to

∇ · v = 0.

Rate of change of velocity The momentum equation is

ρvt+ ρtv + ρv∇v + ρ(∇v)v + ∇p = 0.

Using the continuity equation, we get dv dt = −

∇p

ρ . (1.9)

Rate of change of internal energy and entropy We can subtract the kinetic energy from the energy equation to obtain the motion of internal energy:

ρde

dt + p∇ · v = 0.

This means that the change of internal energy is due to the volume change of the fluid. This together with

de = T dS − pdτ lead to

dS

dt = 0. (1.10)

In compressible flows, dS/dt = 0 means S is constant along particle path, S may still be different on different path. Indeed, as the particle passes through a shock front, the entropy increases.

(11)

1.3.2 Vorticity

The momentum equation can be rewritten as

ρdv

dt + ∇p = 0.

We can rewrite the convection term v · ∇v in the momentum equation by applying the following identity

v · ∇v = ∇ 1 2v

2

− v × ω.

Here, ω := ∇ × v is called the vorticity. This leads to

vt+ ω × v + ∇ 1 2v 2 + 1 ρ∇p = 0. (1.11) In the case of barotropic flows where p is only a function of ρ, the last term ∇p/ρ can be expressed as ∇w. This is an important formula. It says that the acceleration is due to (i) rotation (ω × v) and (ii) conservative force ∇(v2/2 + w).

One consequence from this formula is the Bernoulli law: For barotropic steady flows, the quantity

1 2v

2

+ w = constant, (1.12) along stream lines. Here, steady flow means that the solution is independent of time. We take inner product of (1.11) with v. We can get

v · ∇ 1 2v

2_{+ w}

= 0.

This means that 1₂v2_{+ w stay constant along the stream line.}

Another important result can be derived from (1.11) is by taking curl on both sides, we can eliminate those gradient terms and leave only kinematic variables.

ωt+ v · ∇ω = (ω · ∇) v − ω (∇ · v) +

1

ρ2∇ρ × ∇p. (1.13)

Here, we have used the identity from vector calculus

∇ × (ω × v) = ω(∇ · v) − v(∇ · ω) + (v · ∇)ω − ω · ∇v.

Now, we consider the variation of the quantity ω/ρ. d dt ω ρ = 1 ρ dω dt − 1 ρ2 dρ dtω = 1 ρ((ω · ∇)v − ω∇ · v) + 1 ρ3∇ρ × ∇p + ω ρ∇ · v = ω ρ · ∇ v + 1 ρ3∇ρ × ∇p.

(12)

This kinematic relation is due to Helmholtz. When p is only a function of ρ, ∇ρ×∇p = 0. Such flows are called barotropic flows. The isothermal flows, isentropic flows are examples of barotropic flows. In this case,

d dt ω ρ = ω ρ · ∇ v.

We claim this is equivalent to d dt ω ρ ∂X ∂x = 0. (1.14)

To prove this claim, we use Lagrange coordinate d dt ω ρ = ω ρ · ∇ v = ω ρ ∂X ∂x ∂v ∂X We differentiate the identity

∂X ∂x ∂x ∂X = I in t to get d dt ∂X ∂x ∂x ∂X + ∂X ∂x ∂v ∂X = 0. Plug this into the above equation, we get

d dt ω ρ + ω ρ d dt ∂X ∂x ∂x ∂X = 0. This leads to d dt ω ρ ∂X ∂x = 0. (1.15) We integrate it to obtain ω ρ(t) = ω ρ(0) ∂x(t, X) ∂X . This means that ω

ρ is transported along particle path.

Let us define the vortex line to be the integral curve of the vorticity field ω. Let X(α) be a vortex line at t = 0. That is,

dX dα =

ω(X, 0) ρ(X, 0).

Let us investigate this line following the flow x(t, X(α)). Its tangent is d dαx(t, X(α)) = ∂x ∂X dX dα = ∂x ∂X ω ρ(0) = ω ρ(t)

(13)

1.3.3 Irrotational flows

In the previous subsection, we see that for barotropic flows, if ω ≡ 0 initially, then ω ≡ 0 in all later time. In this case, we call such flows irrotational. That is,

∇ × v ≡ 0.

If the domain is simply connected, then we can find a scalar function φ such that v = ∇φ

Such a function is called the velocity potential. If in addition the flow is incompress-ible, i.e.

∇ · v = 0,

we call an incompressible and irrotational flow a potential flow. The potential theory can be applied.

(14)

Chapter 2 Viscous Fluids

2.1 Viscosity

2.1.1 Source of viscosity

Stress The stress is a surface force due to the impact of particles onto the surface from the other side per unit time. Suppose the small surface has area dA with (outer) normal n. The surface for from particles from outer region is F . The stress is a tensor which maps n to F :

σn = F

The stress can be decomposed into a pressure part and viscous part: σij = −pδij+ σij0

The pressure is isotropic. It exists even the veclocity is zero or a constant. The viscous part is an internal friction. It can occur only when different fluid particles move at different velocities. If v is a constant, then σ0 = 0. Thus σ0 is a function of ∇v. We assume ∇v is small, then it is reasonable to assume σ0 _{as a linear function}

of ∇v.

The rotation does not contribute any viscous force. For a rigid-body rotation:dv_dt = Ω × v, Ω is the rotation vector. The corresponding v satisfies E = 0. Here, E = ∇v + (∇v)T_{. From σ}0 _{being linear to ∇v, zero when ∇v = 0 or E = 0, we get}

σ0 = aE + b tr(E)I.

Here, a and b are two constants. We can modify it a little so that the first term is divergence free.

σ0 = η(∇v + (∇v)T − 2

3∇ · vI) + ζ∇ · vI. (2.1) The coefficient η is called the shear viscosity, while ζ the bulk viscosity. These two coefficient should satisfy

η > 0, ζ > 0. (2.2) for energy dissipation requirement in viscous systems, as we shall see later.

(15)

The viscous force is

∇ · σ0 = η(∇v + (∇v)T − 2

3∇ · vI) + ζ∇ · vI = η∇2v + (ζ + 1

3η)∇ · v. The momentum equation

ρ vt+ ω × v + ∇( 1 2v 2₎ + ∇p = ∇ · σ0

(16)

Chapter 3 Elasticity

The material of this chapter is mainly from Ciarlet’s book, Mathematical Elasticity, Vol. 1, Three dimensional elasticity.

3.1 Strain

Let us imagine an elastic material deforms from a region Ω0 at time 0 to Ω(t) at time

t. We denote such a deformation by y = y(t, x), or y = φt(x). That is, a position with coordinate x is deformed to position y at time t. The coordinate x is called the Lagrange coordinate, while y the Euler coordinate.

The gradient

F := ∇xy(t, x)

is called the deformation gradient. It is clear that F (0, x) = I. Let us denote det F by J . We have J (t) > 0 for all time from physical consideration. Thus the mapping φt_is

invertible. Let us call u := y − x the displacement. Its gradient ∇xu the displacement

gradient. We have F = I + ∇xu.

Examples.

1. Rigid body motion: A trivial deformation is the rigid body motion y = a + Q(t)x

where Q(t) is an orthogonal matrix QQT _{= 1.}

2. Let us consider

y = D(t)x where

˙

D(t) = diag(λ1, λ2, λ3)D(t)

and λi are constants with P_iλi = 0. This gives yi = eλitxi. i = 1, 2, 3. If

λ1 > 0 and λ2, λ3 < 0, then the material expands in x1 direction while squeezs

(17)

3. Shear motion: 4. Expansion

Let us characterize the rigid-body motion without proof.

Theorem 3.1 Let φt _{: Ω → R}3_{. Then φ}t_{(x) = a + Q(t)x is a rigid-body motion if}

and only if FTF ≡ I. Here, F = ∇xφt and Q(t) is an orthogonal matrix.

The geometry of deformation Suppose an infinitesimal vector dx at time 0 is deformed to dy at time t. Then

dy = ∂y

∂xdx = F dx.

The Euclidean inner product at time t induces a time-dependent Riemannian metric at t = 0:

(dy, dy) = (F dx, F dx) = (dx, FTF dx).

We call C := FTF the right Cauchy-Green strain tensor. It is also the first fundamen-tal form of the domain Ω(t). In terms of displacement gradient, C can be expressed as

C = (I + ∇xu)T(I + ∇xu) = I + ∇u + (∇u)T + (∇u)T(∇u).

If the deformation variation is small (i.e. |∇xu| << 1), then

C ≈ I + 2E, E = 1

2 ∇xu + (∇xu)

T_,

E is called the Cauchy strain tensor. Its components are

eij = 1 2 ∂ui ∂xj + ∂uj ∂xi .

Physical meaning of strain. To understand the physical meaning of e11, let us

consider an infinitesimal vector dx1 := (d`1, 0, 0)T. The vector is deformed to

dy1 = (1 + ∂u 1 ∂x1, ∂u2 ∂x1, ∂u3 ∂x1) T_d` 1.

Its length square

(d˜`1)2 = (dy1, dy1) ≈ 1 + 2∂u 1 ∂x1 (d`1)2

Here we have neglected the hight order terms. Thus,

d˜`1 ≈ 1 + ∂u 1 ∂x1 d`1 (3.1) That is e11 ≈ d˜`1− d`1 d`1 .

(18)

Thus, the physical meaning of e11 is the relative change of d` in the direction e1.

Similarly, to understand the physical meaning of e12, we consider two infinitesimal

vectors

dx1 = (d`1, 0, 0)T

dx2 = (0, d`2, 0)T.

The corresponding deformed vectors are

dy1 = (1 + ∂u 1 ∂x1, ∂u2 ∂x1, ∂u3 ∂x1) T_d` 1 dy2 = (∂u 1 ∂x2, 1 + ∂u2 ∂x2, ∂u3 ∂x2) T d`2.

The inner product

(dy1, dy2) ≈ ∂u 1 ∂x2 + ∂u2 ∂x1 d`1d`2. Hence, we get d˜`1d˜`2cos θ = 2e12d`1d`2

where θ is the angle between dy1 _{and dy}2_{. Using (3.1) and neglecting higher order}

terms, we get

e12=

1 2cos θ.

Thus, the physical meaning of e12is the half of the cosine of the relative angle between

the two deformed orthogonal vectors e1 and e2 directors.

compression and shear.

3.2 Stress

3.2.1 The stress tensor in Euler coordinate

The stress is a surface force. It is due to the molecular force of the elastic material on one side of the surface to the other side. More precisely, the stress σ = ∆P/∆S, P is the power, i.e. the force per unit time.

The time derivative of the trajectory φt_{(x), or y(t, x) is the velocity}

v := dy(t, x) dt .

We notice that the velocity defined here is in Lagrange coordinate system. We may also define v(t, y) = dy/dt(t, x) with y = y(t, x), then this is the velocity in Euler coordinate system, and y(t, x) is the solution of the ODE ˙y = v(t, y) with y(0, x) = x.

The momentum equation reads Z Ωt ρdv dt dy = Z ∂Ωt σ · ν dSy.

(19)

Here, d/dt means the material derivative (i.e. the partial derivative in t with fixed x) The stress is indeed a symmetric 2-tensor, which is based on the following three arguments:

• Cauchy stress principle: The stress is assumed to be a function of the Euler coordinate y and the normal of the surface ν.

• From conservation of linear momentum, we can derive that

σ(y, ν) = T (y)ν. (3.2)

• From conservation of angular momentum, we can get T (y) is symmetric. To show (3.2), given ν, we consider a plan with normal ν. The plan intersects the coordinate axes e1, e2, e3 at A1, A2, A3, respectively. Let us consider the perimid

OA1A2A3, call it ∆G, the surface A1A2A3 is called ∆S, and the three sides of the

perimid are ∆Si. The conservation of linear momentum on ∆G reads

Z ∆G ρdv dt dy = Z ∆S σ(y, ν) dS + 3 X i=1 Z ∆Si σ(y, −ei) dS

By dividing both sides by ∆S then take ∆S → 0, we can get

σ(y, ν) = −

3

X

i=1

νiσ(y, −ei)

On the other hand, from Newton’s third law, we can get σ(y, −ei) = −σ(y, ei). We

conclude that σ is linear in ν.

The conservation of angular momentum reads Z G ρ y × dv dt dy = Z ∂G (y × T ν)dS In differential form: ρy × dv dt = ∇y(y × T ) The last term in coordinate form is

∂l(ijkyjTkl) = ijkδjlTkl+ ijkyj∂lTkl

= ijkTkj + ijkyj∂lTkl

On the other hand, from linear momentum equation:

ρdv

(20)

we take cross product with y on both sides to get

ρy ×dv

dt = y × (∇yT ). This together with the angular equation, we get

ijkTkj = 0,

That is Tij = Tji. Thus, by assuming the conservation of linear momentum, the

conservation of angular momentum is equivalent to the symmetry of the stress tensor.

3.2.2 The stress tensor in Lagrangian coordinate–Piola stress

We want to express the linear momentum equation in Lagrangian coordinate. Let G0

is deformed to Gt, its boundary S0 to St. In Euler coordinate, the linear momentum

equation is Z Gt ρdv dt dy = Z St T ν dSy

In Lagrange coordinate, we have Z G0 ρ0 dv dt dx = Z S0 P n dS0 (3.3)

We want to find the relation between T and P . The tensor P is called the Piola stress tensor.

First, the conservation of mass reads Z Gt ρ dy = Z G0 ρ0dx.

The left-hand side isR

GtρJ dx. This leads to ρJ = ρ0. Using this, Z Gt ρdv dt dy = Z G0 ρdv dtJ dx = Z G0 ρ0 dv dt dx.

In (3.3), the velocity v at time t is the same vector, the stress σ is also the same vector in the Lagrange coordinate. We now change the surface force term to the Lagrange coordinate. We claim that

νdSt = J F−Tn dS0 (3.4)

With this, then

P = J T F−T.

Now, we prove (3.4). Let dx1 _{and dx}2 _{are two infinitesimal vectors such that}

(21)

Suppose dxi is deformed to dyi at time t, then

νdSt= dy1 × dy2.

In coordinate form, we have

νidSt = ijkdyj1dy 2 k = ijkyj,qyk,rdx1qdx 2 r

Multiplying both sides by yi,p and using

J = pqry1,py2,qy3,r, we obtain νiyi,pdSt= pqrJ dx1qdx 2 r = J npdS0 Thus, we get FTνdSt= J ndS0.

The equation of motion now reads

ρ0

dv

dt = ∇xP. From the symmetry of T , we get

Σ := F−1P = J F−1T F−T,

called the second Piola stress tensor, which is also symmetric.

3.3 Stress-strain relation

The stress in an elastic material is due to deformation. Therefore we assume the stress tensor T is a function of the deformation gradient F . Such a relation

T (y) = bT (x, F )

is called a constitutive relation and bT is called a response function. If bT is independent of x, the material is called homogeneous. Below, we will find the physical conditions for a response function.

3.3.1 Frame indifference

First, the stress should be independent of the frame of reference physically. Suppose y is the coordinate in the original frame, and y∗ is the coordinate in new frame which is a rotation of the old frame. That is, y∗ = Qy, and Q is a rotation. Then

(22)

F∗ := ∂y∗/∂x = QF ; the normal and the stress in the old and new frames are related by ν∗ = Qν, σ∗ = Qσ. In terms of the response function, it reads

σ∗ = T (QF )νb ∗ = bT (QF )Qν Qσ = Q bT (F )ν.

Hence,

b

T (QF ) = Q bT (F )QT. (3.5) This condition is called the axiom of material frame-indifference.

Next, we will show that bT is indeed a function of the Cauchy-Green strain tensor only. Let us write F in polar form

F = RU.

where R is a rotation and U is positive definite. From the axiom of material frame-indifference, we get

b

T (F ) = bT (RU ) = R bT (U )RT = F U−1T (U )U Fb T = F T (C)FT. where C = U2 = FTF is the right Cauchy-Green strain tensor.

Similarly, for the Piola tensor P , the frame-indifference reads

P (QF ) = J bT (QF )(QF )−T = J (Q bT (F )QT)Q−TF−T = QP (F ). and

P (F ) = J bT (F )F−T = J F T (C) = F P (C),

where J2 _{= detC and P (C) :=}√_{detC T (C). For the second Piola tensor it satisfies}

Σ(QF ) = (QF )−1P (QF ) = (F−1Q−1)(QP (F )) = Σ(F ), and

Σ(F ) := F−1P = P (C).

3.3.2 Isotropic material

A material is isotropic if its stress tensor is independent of a rotational of the strain. Suppose we have two identical materials. One has material coordinate x, while the other is x∗ with x∗ = Qx and Q being a rotation. If these two materials deform in the same way and result the same stresses, we call this material isotropic. Mathematically, this means

b

T (F Q) = bT (F ). (3.6) In terms of Piola stress tensor, it reads

P (F Q) = P (F )Q for all orthogonal matrix Q (3.7) This follows from

P (F Q) = J bT (F Q)(F Q)−T = J bT (F )F−TQ = P (F )Q.

In order to characterize the response function of an isotropic material, we review some spectral properties of a matrix.

(23)

3.3.3 Hyper-elastic material

A material is called hyper-elastic material if there exists a potential W (F ) such that

P (F ) = ∂W ∂F .

Notice that such materials satisfy the frame-indifference axiom if and only if W (QF ) = W (F ) for all Q.

To see this, we differentiate the above equation in F to get

P (F ) = ∂ ∂FW (F ) = ∂ ∂FW (QF ) = Q T∂W (QF ) ∂(QF ) = Q T_{P (QF ).}

This shows W (F ) = W (QF ) if and only if P (QF ) = QP (F ) for all Q. Similarly, a hyper-elastic material is isotropic (i.e. P (F Q) = P (F )Q) if and only if W (F Q) = W (F ) for all Q. This follows from

P (F ) = ∂ ∂FW (F ) = ∂ ∂FW (F Q) = ∂W (QF ) ∂(QF ) Q T = P (QF )QT.

Thus, for an isotropic hyper-elastic material, its potential energy function W satisfies

W (F Q) = W (QF ) = W (F ) for all Q (3.8)

Thus, W is only a function of the principal invariants of F , the eigenvalues of FTF . Two simple examples:

• Hook’s law: W (F ) = 1 2H|F |

2_{, where |F |}2 ₌P

ij|Fij|2.

• Linear isotropic elastic material:

W = λ 2(tr E)

2_{+ µ tr (E}2_{), where E =} 1

(24)

Chapter 4 Variation Formulation of Fluid

Mechnics and Elasticity

4.1 Flow maps

4.1.1 Lagrangian and Eulerian Coordinates

We will derive a unified theory for both fluid and elasticity. We shall use the following notations.

• Let x(·, X) be the particle path of a piece of material (or fluid) with position X initially. X is called the Lagrangian coordinate, whereas x the Eulerian coordinate. It is assumed that the mapping satisfies x(0, X) = X.

• Ωt: the region occupied by the material at time t. Thus, the mapping x(t, ·) :

Ω0 → Ωt. Or the mapping

x(·, ·) : [0, T ] × Ω0 → ∪t∈[0,T ]Ωt.

In most of case, we assume the whole region occupied by the material is un-changed. We denote this wholeregion by D. Thus, x : [0, T ] × D → D. We always assume x(t, ·) is 1-1 and onto from D to D.

• The deformation matrix F := ∂x ∂X, or F

i

α(t, X) = ∂xi(t, X)/∂Xα.

• The Jacobian J := det F . We have dx = JdX.

• The particle velocity is defined to be ˙x(t, X). We shall use dot for the time derivative with the Lagrangian coordinate fixed.

• In Eulerian coordinate, we define

u(t, x(t, X)) = ˙x(t, X). In other words, ˙x(t, X) = Z D u(t, x)δ(x − x(t, X)) dx where δ is the delta function.

(25)

• For a physical quantity, say the density ρ, we denote by ρ0 the density initially,

and the density ρ(t, x) is defined to be

ρ(t, x) = Z

ρ0(X)δ(x − x(t, X)) dX.

Thus, the mass in a region Ω is Z Ω ρ(t, x) dx = Z Ω Z ρ0(X)δ(x − x(t, X)) dX dx = Z x(t,X)∈Ω ρ0(X) dX

If we take Ω = Ωt, then we get

Z Ωt ρ(t, x) dx = Z Ω0 ρ0(X) dX.

Thus, the definition above implicitly assumes conservation of mass. By deviding the above equation by |Ω0|, and use limΩ0→0

|Ωt|

|Ω0| = J , we get

ρ(t, x(t, X))J (t, X) = ρ0(X).

Alternatively, we take Ω = [x1_{, x}1_{+ ) × [x}2_{, x}2_{+ ) × [x}3_{, x}3_{+ ). Then}

Z x1+ x1 Z x2+ x2 Z x3+ x3 ρ(t, x) dx = Z xi_≤xi_(t,X)<xi₊ ρ0(X) dX = Z xi_≤xi_(t,X)<xi₊ ρ0(X)J−1dx

By dividing both sides by 3, we get ρ(t, x) = ρ0(X)J−1(t, X), where x =

x(t, X).

• Internal energy W . The deformation F stores an internal energy W in the material. It is due to the strenching between atoms. This internal energy W is assumed to be a function of F . Examples – Hook’s law: W (F ) = k 2 X |Fi α|2 = k 2tr(F T_{F ).}

– St. Venant-Kirchhoff: we assume small variation, i.e. x(t, X) − X is small. Let E = 1₂(FTF − I).

W (F ) = λ 2(trE)

(26)

• The stress tensor. The derivative P = W0_{(F ) is called the Piola stress tensor.}

It is a surface force in Lagrange coordinate system. The stress in Eulerian coordinate system is denoted T , called the Cauchy-Green stress tensor. The surface force on the surface with normal ν is T ν. In Lagrange coordinate, it is P n, where n is the normal in Ω0 which is deformed to ν at time t. T νd2S

and P nd2_S

0 represent the same surface force. Thus, T νd2S = P nd2S0. Notice

that νd2S is a 2-form in Ωt, and nd2S0 is its pull back by the flow map x(t, X).

Thus, the pullback of νd2_{S is J}−1_FT_νd2_{S = nd}2_S

0. Thus,

T νd2S = P nd2S0 = P J−1FTνd2S,

We then obtain

T = J−1P FT.

4.1.2 Variation of Functionals w.r.t. the flow maps

Functionals Given a flow map x(·, ·), we define its kinetic energy, internal energy, action, total energy in the Lagrange coordinate as the follows.

• Kinetic energy: R Ω0 1 2ρ0(X)| ˙x(t, X)| 2_dX • Internal energy: R Ω0W (F ) dX • Action: A[x] = Z T 0 Z Ω0 1 2ρ0(X)| ˙x(t, X)| 2_{− W (F ) dX dt} • Total energy: E[x(t)] = Z Ω0 1 2ρ0(X)| ˙x(t, X)| 2_{+ W (F ) dX.}

Variation w.r.t. flow maps We shall study the variation of these functionals with respect to the flow map x(·, ·). Let us perturb the flow map by x(t, X) with

x0(t, X) = x(t, X), the original unperturbed one. We call

δx(t, X) := ∂

∂|=0x(t, X),

the variation of the flow map x(·, ·). We write the variation of x by δx. The variation of a functional I[x] in the direction of δx means that

δI[x] · δx := d

d|=0I[x]. Since, for small , x(t, ·) are flow maps, its variation

δx(t, X) = lim

→0

x− x0

is an infinitesimal variation of position, thus, can be called a pseudo-velocity. Let us thus define the corresponding pseudo-velocity in the Eulerian coordinate as

(27)

• δF = (∇v)F , where ∇ is the abbreviation of ∇x. From δx(t, X) = v(t, x(t, X)), we differentiate in X to get ∂δxi ∂Xj = ∂vi ∂xk ∂xk ∂Xj

Interchange ∂ and δ operators, we get δF = (∇v)F . • Let J = det F , then

δJ = tr (δF )F−1 J. This follows from the following lemma.

Lemma 4.1 Let A() be a smooth n × n matrix-valued function. Then d

ddet(A) = tr( ˙AA

−1

) det(A)

Proof. We write A = (aij). First, we recall the expansion formula of det A:

X

k

aikAjk = (det A)δij

where Aij is the signed cofactor of A at (i, j). That is, A(cof A)T = (det A)I.

We claim that ∂ det(A)/∂aij = Aij, To see that, we write det A as

det(A) =X

k

aikAik,

and notice that Aik does not involve aij for all k. Thus,

∂ det A ∂aij = Aij. Next, d ddet A = X ij daij d ∂ det(A) ∂aij =X ij ˙aijAij = X ij ( ˙A)ij((cof A)T)ij = X ij ( ˙A)ij(A−1)ijdet A = tr( ˙AA−1) det A

4.2 Inviscid Flows and Elasticity

4.2.1 Equation of Motion

1. Compressible elasticity We shall take variation of action w.r.t. flow map (position). This will give us the equation of motion

(28)

Lagrange formulation. We take the variation δx satisfying δx(t, X) = 0 for t = 0 and t = T . We can take integration by part for t.

δA[x] · δx = Z T 0 Z Ω0 (ρ0(X) ˙x(t, X) · δ ˙x − W0(F )δF ) dX dt = Z T 0 Z Ω0 −d dtρ0(X) ˙x(t, X) · δx − W 0 (F )δ∂x ∂X dX dt = Z T 0 Z Ω0 −d dtρ0(X) ˙x(t, X) · δx + ∂ ∂XW 0_{(F )} · δx dX dt

Here, we also choose W0(F )δx · n = 0 on the boundary of Ω0 so that no boundary

terms show up in the integration-by-part for the integration in X over Ω0. This leads

to the Euler-Lagrange equation

ρ0(X)¨x(t, X) = ∇X · P

∂x ∂X

.

Here, P = W0(F ) is the Piola stress tensor. In component form,

P_αi = ∂W ∂Xi α , (∇X · P )i = ∂ ∂XαP i α.

This is the equation of motion for an elastic material in Lagrange coordinate system. We can also use the variable u(t, X) := ˙x(t, X) and F (t, X) = _∂X∂x(t, X) as the new unknowns to write the above equations as a first-order system. Notice that the velocity u(t, X) is in Lagrange coordinate in this formulation. The equation of motion is

ρ0(X) ˙u(t, X) = ∇X · P (F ).

By differentiate ˙x(t, X) = u(t, X), we obtain the equation for F is

˙ Fij =

∂ui

∂Xj

(t, X).

Two special cases:

1. One dimension: In one dimension, let us call v = F1

1 − 1 and P (f11) = σ(v).

Then the equation of motion are

ρ0(X)ut = σ(v)X

vt = uX

Such an equation is a 2×2 hyperbolic system. The stress σ(v) satisfies σ0(v) > 0, but is a nonvex function of v in general. A shock wave attached to a rarefaction wave can be found.

(29)

2. Linear elasticity: In the Hook’s case where W (F ) = k₂|F |2_{, the stress is P (F ) =}

kF . The equation of motion becomes

ρ0(X)ut = k∇X · F

Ft = ∇Xu

Eliminating F , we get

ρ0(X)utt = k∇2Xu

This is the wave equation.

Eulerian formulation. In the Eulerian formulation, u is defined as u(t, x(t, X)) = ˙x(t, X), the perturbation (or the pseudo-velocity) v is defined as v(t, x(t, X)) = δx(t, x(t, X)). Let us define the material derivative as

D Dt :=

∂

∂t+ u · ∇.

The gradient is with respect to the Eulerian coordinate x. Thus, we have Du

Dt(t, x(t, X)) = ¨x(t, X).

We can use delta function to express δx and ¨x in terms of v and Du/Dt as the follows.

δx(t, X) = Z v(x, t)δ(x − x(t, X)) dx ¨ x(t, x(t, X)) · δx(t, X) = Z Du Dt · v(x, t)δ(x − x(t, X)) dx The force from the stress tensor is

G(t, X) = ∇X · P (t, X)

The variation formula Z (ρ0(X)¨x(t, X) − G(t, X)) · δx(t, X) dX can be re-expressed as Z Z ρ0 Du Dt(t, x) − G(t, X) · v(t, x)δ(x − x(X, t)) dx dX Let us define ρ(t, x) = Z ρ0(X)δ(x − x(t, X)) dX g(t, x) = Z G(t, X)δ(x − x(t, X)) dX

(30)

This means ρ(t, x) = lim x ∈ Ω |Ω| → 0 1 |Ω| Z x(t,X)∈Ω ρ0(X) dX

Then the above variation formula becomes Z T 0 Z ρ(t, x)Du Dt(t, x) − g(t, x) · v(t, x) dx dt = 0.

This is true for all pseudo-velocity v. We thus obtain the equation of motion in Eulerian coordinate as

ρDu Dt = g. The stress force

g(t, x) = Z G(t, X)δ(x − x(t, X)) dX = Z Ω0 (∇X · P (t, X))δ3(x − x(t, X)) dX Here, δ3_{(x) = δ(x}1_)δ(x2_)δ(x3_). gi(t, x) = Z ∂ ∂XαP i αδ3(x − x(t, X)) dX = Z P_αi 3 X k=1 Y `6=k δ(x`− x`_{(t, X))δ}0 (xk− xk_{(t, X))F}k αdX

Now, let write

δ0(xk− xk(t, X)) = lim h→0 1 h δ(x k + h − xk(t, X)) − δ(xk− xk(t, X)) and R Pi αFαk Q `6=kδ(x`− x`(t, X))δ(xk+ h − xk(t, X)) = lim→0 13 R xk_+h<xk_(t,X)<xk_+h+ R x`_<x`_(t,X)<x`₊P_αiF_αkdX = lim→013 Rx`+ x` Rxk+h+ xk_+h P i αFαkJ−1dx Similarly, we get R Pi αFαk Q `6=kδ(x `_{− x}`_{(t, X))δ(x}k_{− x}k_{(t, X))} = lim→0 13 R xk_<xk_(t,X)<xk₊ R x`_<x`_(t,X)<x`₊P_αiF_αkdX = lim→0 13 Rx`+ x` Rxk+ xk P i αFαkJ −1_dx

By taking → 0 and h → 0, we obtain

gi(t, x) = ∂ ∂xk(P i αF k αJ −1 ),

(31)

or g = ∇x· T , where T = P FTJ−1.

By differentiating the equality

˙x(t, X) = u(t, x(t, X),

in X, we get the equation for F :

˙

F = (∇u)F.

In terms of the Euler coordinate system, it is

DF

Dt = (∇u)F.

Now, ρ is a new unknown. It satisfies ρ(t, x)J (t, X) = ρ0(X). Differentiate this

in t with fixed X, we get

0 = (ρt+ ∇xρ · ˙x) J + ρ ˙J

= (ρt+ ∇xρ · u) J + ρ(∇ · u)J.

Eliminating J , we get the continuity equation

Dρ

Dt + ρ∇ · u = 0. Thus, the full equations are

   Dρ Dt = −ρ∇ · u DF Dt = (∇u)F ρDu_Dt = ∇ · T

where T = P (F )FT_J−1_{. The unknowns are ρ, F, u. These equations can also be}

written in conservation form:    ∂tρ + ∇ · (ρu) = 0 DF Dt = (∇u)F ∂t(ρu) + ∇ · (ρu ⊗ u) = ∇ · T.

2. Incompressible elasticity A flow map x(t, X) is called incompressible if

det ∂x ∂X(t, X)

= 1.

(32)

Lagrange formulation The constraint is

det F (t, X) = 1.

Thus, in the above variation of action, we should add a constraint term with a La-grange multiplier: δA[x] + δ Z T 0 Z p(t, X)(det F − 1) dX dt = 0.

Here, p is the Lagrange multiplier. The variation

δ(det F ) = tr((δF )F−1) det F = tr((δF )F−1). ((δF )F−1)i_j = ∂δx i ∂Xk(F −1 )k_j Thus, tr((δF )F−1) =X i X k ∂δxi ∂Xk(F −1 )k_i

We take integration by part in the variation form below to get

δ Z T 0 Z p(t, X)(det F − 1) dX dt = Z T 0 Z pX k X i ∂δxi ∂Xk(F −1 )k_i dX dt = − Z T 0 Z X k ∂ ∂Xk p(F −1 )k_i δxidX dt = − Z T 0 Z ∇X · p(F−1) · δx dX dt

Thus, we obtain the constraint flow equation

ρ0x = ∇¨ X · (P − pF−1).

The full set of equations are    ρ0˙u = ∇X · (P (F ) − pF−1) ˙ F = ∇Xu det F = 1

The unknowns are (F, u, p).

Eulerian formulation If u(t, x) is the velocity field which generates the flow map x(t, X), that is, ˙x(t, X) = u(t, x(t, X)), then the incompressibility of x(t, X) is equivalent to ∇ · u(t, x) = 0. This comes from the following formula

(33)

and the formula ∂ ˙x ∂X = (∇u) ∂x ∂X, we get ˙ J = tr ˙F F−1J = tr(∇u)J = (∇ · u) J. Thus, ˙ J = 0 if and only if ∇ · u = 0.

Similarly, let x(t, X) be a perturbation of x satisfying volume preserving property.

Let δx(t, X) := _∂∂|=0x(t, X). Let v(t, x(t, X)) = δx(t, X). Let F = ∂x/∂X.

δF := _∂∂|=0F. Then

δF = (∇v)F. Similarly,

δJ = (∇ · v)J.

Thus, the volume preserving of x is equivalent to δJ = 0, and is also equivalent to

∇ · v = 0.

Now, for the constraint ∇ · v = 0, we introduce a Lagrange multiplier p. Then we have the constrained variation

0 = Z T 0 Z (−ρDu Dt + ∇ · T ) · v + p(∇ · v) dx dt = Z T 0 Z (−ρDu Dt + ∇ · T − ∇p) · v dx dt This gives ρDu Dt + ∇p = ∇ · T.

Here, T (F ) = W0(F )FT. We still need an equation for F . From ˙ F = ∂x ∂X t = ∂u ∂X = ∂u ∂x ∂x ∂X = (∇u)F, we get DF Dt = (∇u)F.

Thus, the complete set of equations for incompressible inviscid elasticity are        Dρ Dt = 0 ∇ · u = 0 ρDu_Dt = −∇p + ∇ · T (F ) DF Dt = (∇u)F.

with the constitutive equation T (F ) = W0(F )FT_{. The unknowns are (ρ, p, u, F ). We}

notice that these equations can be written in conservation form        ∂tρ + ∇ · (ρu) = 0 ∇ · u = 0 ∂t(ρu) + ∇ · (ρu ⊗ u) = −∇p + ∇ · T (F ) DF Dt = (∇u)F.

(34)

3. Incompressible fluid The only difference between fluid and elastic material is the stress term. Consider the ideal fluids, where no stress appears. Then the equation of motion of an ideal fluid in Eulerian coordinate reads

ρDu

Dt + ∇p = 0. ∇ · u = 0. The continuity equation reads

ρt+ ∇ · (ρu) = 0.

Here, we have (ρ, p, u) as our unknowns. Two special cases.

• ρ ≡ 1, we have

Du

Dt + ∇p = 0, ∇ · u = 0.

• Barotropic flow: the pressure is a function of ρ. In this case, we can find a potential Φ such that Φ0(ρ) = p0(ρ)/ρ. Then the barotropic flow equation becomes

Du

Dt + ∇Φ = 0, ∇ · u = 0.

4.2.2 Boundary conditions

1. Compressible invisid elasticity The natural boundary condition can be found through the variation formulation. The natural condition is

W_kiδxi· nk _{= 0 for X ∈ ∂Ω} 0

for all perturbation δx. Here, n is the outer normal of Ω0. There are many possibilities

to achieve such condition.

• Dirichlet boundary condition:

x(t, X) = x0(X), X ∈ ∂Ω0,

or equivalently,

˙x(t, X) = 0, X ∈ ∂Ω0.

in Lagrange coordinate system. Or

u(t, x) = 0 x ∈ ∂Ωt.

in Euler coordinate system.

• In Lagrange coordinate system: Pi kn

k_{= 0.}

In Euler coordinate system: Ti

kνk= 0.

(35)

2. Incompressible invisid elasticity The natural boundary condition is to have

(P − pF−1) · nδx = 0

u · ν = 0, x ∈ ∂D. There is only one boundary condition on the boundary.

3. Incompressible invisid fluid

u · ν = 0, x ∈ ∂D.

There is only one boundary condition on the boundary.

4.2.3 Energy conservation laws

The energies E(t) defined below are conserved. • Compressible/incompressible inviscid elasticity

E(t) = Z Ω0 1 2ρ0(X)| ˙x(t, X)| 2_{+ W (} ∂x ∂X) dX or E(t) = Z Ωt 1 2ρ(t, x)|u(t, x)| 2_{+ W J}−1 dx.

• Incompressible inviscid fluid

E(t) = Z Ω0 1 2ρ0(X)| ˙x(t, X)| 2_dX or E(t) = Z Ωt 1 2ρ(t, x)|u(t, x)| 2_dx.

4.3 Visco-elasticity

4.3.1 Equation of Motion

In classical mechanics, if a system has a damping term, then we have dE(t)

dt = −γ| ˙x|

2 _{⇔ m¨}_{x + V}0_{(x) = −2γ ˙x.}

(36)

• Incompressible viscous fluid With the incompressibility ∇ · u = 0 and con-tinuity equation Dρ/Dt = 0, the equation

ρDu dt + ∇p = ∇ · (µ∇u). is equivalent to d dt Z ₁ 2ρ|u| 2_{dx = −} Z µ|∇u|2dx.

• Incompressible visco-elasticity With Dρ

Dt = 0 ∇ · u = 0

DF

Dt = (∇u)F the momentum equation

ρDu Dt = −∇p + ∇ · T + ∇ · (µ∇u) is equivalent to d dt Z ₁ 2ρ|u| 2_{+ W (F ) dx = −} Z µ|∇u|2dx.

4.3.2 Boundary conditions

Due to the viscous term in the momentum equation, it is natural to impose velocity on the boundary to be

u(x) = 0, x ∈ ∂D.

This consists of three condition because u ∈ R3. There are no boundary condition for ρ, F, p for compressible elasticity, and no boundary condition for p for incompressible elasticity.

(37)

(38)

Chapter 5 Membrane

5.1 Energy law of membrane

5.1.1 Surface geometry

We assume X = (X1, X2) be two dimensional and x = (x1, x2, x3) be three

dimen-sional. The deformation matrix F = ∂x/∂X is a 3 × 2 matrix. The metrix (dx, dx) induces a metric gαβ = ∂x ∂Xα · ∂x ∂Xβ

on the manifold Σt= {x(t, ·)}. That is,

(dx, dx) = (F dX, F dX) = (FTF dX, dX) =X

jk

gαβdXαdXβ.

The matrix FT_{F is called the first fundamental form of the surface Σ}

t. The area

spanned by the vectors vα = ∂x/∂Xα and vβ = ∂x/∂Xβ is p

(vα_{, v}α_{) (v}β_{, v}β_{) − (v}α_{, v}β₎2 ₌p_det(FT_{F ).}

Thus, we define J =pdet(FT_{F ). Then dS}

t= J dS0.

We want to define an energy W (F ) on the manifold Σt.

The frame indifference condition reads W (QF ) = W (F ),

where Q is any 3 × 3 orthogonal matrix. Taking the singular value decomposition of F , F = U ΛVT, where U ∈ O(3), V ∈ O(2) and Λ is a 3 × 2 matrix

Λ =   λ1 0 0 λ2 0 0   we then have W (F ) = W (U ΛVT) = W (ΛVT)

(39)

Isotropic condition reads The isotropy assumption is W (F Q1) = W (F ),

where Q1 is any 2 × 2 orthogonal matrix. If we represent F = U ΛVT, then W (F ) =

W (Λ). That is, The isotropy condition implies that W only depends on the eigenval-ues of g.

Two simple examples • Hook’s law: W (F ) = 1 2tr(F T_{F ) =} 1 2 2 X α=1 gαα = 1 2 2 X α=1 3 X i=1 |Fi α| 2_.

The corresponding Piola stress is

P_αi = ∂W ∂Fi α = F_αi, α = 1, 2, i = 1, 2, 3. • minimal surface W (F ) = σ√g = σpdet(FT_{F ) = σJ.}

Here, σ is the surface tension. The energy is Z Ω0 W (F ) dX = Z Ω0 σJ dX = Z Ωt σ dSt • Fabric: W (F ) = a(tr(FTF ) − 2)2+ b(pdet(FT_{F ) − 1)}2_.

Another fabric model is to decompose g = FT_{F into a divergence part and a}

divergence free part. That is,

g = 1 2tr(g)I + g − 1 2tr(g)I

The first part is related to the expansion or shrinking of the surface, while the second is related to the distorsion of the surface. We can design the energy to be W (F ) = a(tr(g) − 2)2+ b g − 1 2tr(g)I 2 F

The norm here is the Frobenius norm. • Vesicle model (Canham-Helfrich)

W (F ) = κ

2(H − H0)

2

+ ¯κK.

where H is the mean curvature and K the Gaussian curvature. In this formula-tion, W depends on the second fundamental form of the surface, which, by the intrincit property of surfaces, involves the derivative of F . 1

1_{For membrane model, a reference book is: Statistical Mechanics of Membranes and Surfaces,}

(40)

5.1.2 Isotropic membrane

Let x(X, t) be the flow map with X ∈ Σ0 being two dimensional.

Piola stress Pi α = ∂F∂Pi

α.

Relation between Cauchy stress and Piola stress The Piola stress is con-vinient in the Lagrange coordinate system. In Eulerian coordinate system, the cor-responding stress is the Cauchy stress T .

T = J−1P FT. This is to be verified.

5.1.3 Equation of Motion

Lagrange formulation The action is A[x] = Z T 0 Z Σ0 1 2ρ0(X)| ˙x(t, X)| 2_{− W} ∂x ∂X dX dt

The variation of the action with respective to the membrane motion x(t, X) gives the momentum equation. One can see the previous Euler-Lagrange equation is still the same.

ρ0(X)¨x(t, X) = ∇X · P

∂x ∂X

Let u(t, X) := ˙x(t, X). The equation of motion becomes ρ0(X) ˙ui(t, X) = _∂X∂αP i α(F ) ˙ Fi α = ∂u i ∂Xα Euler formulation

Let Σt be the membrane at time t. We assume that there is an 1-1 and onto mapping

x(t, X) between Σ0 and Σt. The Jacobian J 6= 0. We define the following quantities

on Σt.

1. Velocity: u(t, ·) is defined on Σt by ˙x(t, x(t, X)) if x = x(t, X). The velocity

u(t, x) may not lie on T Σt(x).

2. Pseudo-velocity v(t, ·) is defined on Σt by δx(t, x(t, X)) if x = x(t, X). The

pseudo-velocity may not lie on T Σt(x).

3. We define material derivative D/Dt by ∂/∂t + u · ∇Σ. Here, ∇Σ is the surface

(41)

4. The volume d3x can be expressed as dSdν, where dS is the area element of the surface Σt and dν is te arc length element in the normal direction. The

δ-function in 3-d can be expressed as

δ3(x − x(t, X)) d3x = δS(x − x(t, X))δν(x − x(t, X))dS dν. 5. We have δx(t, X) = Z Σt v(t, x)δS(x − x(t, X)) dS ¨ x(t, X) · δx(t, X) = Z Σt Du Dt · v(t, x)δS(x − x(t, X)) dS 6. The Piola stress is pushed forward as the Cauchy-Green stress T by

T = J−1P F.

7. The density ρ on the surface is defined to be

ρ(t, x) = Z

ρ0(X)δS(x − x(t, X)) dS0

With these, the equation of motion is

ρDu

Dt = ∇S · T. This is to be clearified.

5.2 Two-phase Flows

5.2.1 Inviscid flows

Consider two simple fluids occupied region Ωi, i = 1, 2 connected by an interface Γt.

In Ωi, the fluid is governed by

∇ · u_i = 0

ρi(ui,t+ ui · ∇ui) + ∇pi = 0

We assume ρi are constants. On the boundary Γt, the divergence free ∇ · ui = 0 leads

to

[u] · ν = 0,

where [u] = u2 − u1. This implies that the motion of the interface, which can be

characterized by its nomal velocity vn, satifies

(42)

For force balance on the interface, we first decompose a vector field V into Vν and

VS and the divergence is decomposed into ∇νVν + ∇SVS. Similar for a tensor. Now,

across an interface

[p]ν = T · ν where T is the stress on the surface.

Given vn, it will give the boundary condition ui· ν = vn on the bouundary ∂Ωi

and will be able to solve the Euler equation. The interface normal velocity vnin turn,

is determined by the other interface condition: [p] − σuH.

5.2.2 Viscous flows

Now we consider two viscous flows connected by an interface Γt. The governed

equa-tions in each region are the Navier Stokes equation: ∇ · u_i = 0

ρi(∂tui+ ui· ∇ui) + ∇pi = ∇ · µi∇ui.

On the boundary Γt, due to viscosity, we should use the condition

[u] = 0. This leads to

u1 = u2 = v,

where v is the surface velocity. From momentum equation, the surface force due to an elastic material on Γt is σHν (assume the elastic material is inviscid). Let

τi := µi∇ui − piI be the stress tensor of the fluid i in region Ωi. Then the force

balance law on the interface is

[τ ] · ν = σHν.

If v is given, then the condition ui = v and N-S equation can determine the flow in

Ωi. On the other hand, there are three conditions in [τ ] · ν = σHν to determine the

(43)

(44)

Chapter 6 Phase Field Models

6.1 Dynamics of multiphase materials

Order parameter A two-phase fluid model: In fluid system, we may encounter two different kinds of fluids in a fluid system. For example, oil and water, water and gas, or even have more compositions. We may use a composition function to indicate which material, which is called order parameter in physics literature. For instance, in a water-oil system, we may use φ = 1 for water and φ = −1 for oil. One can associate a free engrgy W with φ by

F [φ] = Z Ω f (φ, ∇φ) dx = Z Ω 1 2|∇φ| 2₊1 W (φ) dx

The energy W is the bulk free energy. A natural choice is a double well potential with minima at −1 and 1, for instance W (φ) = 1

4(φ

2_{− 1)}2_{. The energy |∇φ|}2 _{is the}

surface energy. It means that same phase material likes to group together. It costs some energy by putting materials with different phase together. When φ tends to an equilibrium with 1 and −1 on the two sides of an interface, ∇φ becomes the delta function on the interface and F [φ] measures the area of the interface.

Chemical potential The variation of F w.r.t. φ is called the chemical correspond-ing to φ. µ = δF δφ = −∇ · ∂f ∂∇φ+ ∂f ∂φ.

Evolution of order parameter The simplest case is ∂tφ + u · ∇φ = 0.

This is equivalent to

φ(t, x(t, X)) = φ0(X).

The other possibility is

∂tφ + ∇ · (uφ) = 0,

which yields

(45)

Variation of F w.r.t. flow motion Let xs be a one-parameter family of flow

maps. δxs= x0s. Let vs(t, x(t, X)) = δxs(t, X). The order parameter φs is defined by

φs(t, xs(t, X)) = φ0(X).

This means that the flow perturbation does not change the phase. Hence, φs satisfies

∂sφs+ vs∇φs = 0. The variation of F w.r.t. xs is δF · v = δF δφs ∂φs ∂s = µ(−vs∇φs) = −µ∇φ · v

Variation of F w.r.t. predomain In this case, we have φs(t, xs(t, X))Js(t, X) = φ0(X). This implies 0 = Js(∂sφs+ vs∇φs) + Js0φs = Js(∂sφs+ vs∇φs) + (∇ · v)Jsφs This gives ∂sφs+ ∇ · (φsvs) = 0. The variation of F w.r.t. xs _is δF · v = δF δφs ∂φs ∂s = µ(−∇(φsvs)) = φ∇µ · v

Incompressible case In this case,

δF · v = φ∇µ · v = −µ∇φ · v.

6.2 Equation of Motion

The action is defined to be

A[x] = Z T 0 Z Ω0 1 2ρ0(X)| ˙x(t, X)| 2_{− f (φ, ∇φ)(t, x(t, X))J(t, X) dX dt.}

The potential energy is indeed the functional

F (φ) = Z

Ωt

(46)

The variation of action w.r.t. flow motion gives

ρDu

Dt = µ∇φ. This together with

∂tφ + u · ∇φ = 0

∂tρ + ∇ · (ρu) = 0,

give the set of equations for (ρ, φ, u).

Conservation of energy We multiply momentum equation by u then integrate in X and in t to get d dt Z 1 2ρ|u| 2_{dx =} Z µu · ∇φ dx = − Z µ∂tφ dx = − Z δF δφ∂tφ dx = − Z dF dt dx

This gives the conservation of total energy. It also shows that the force µ∇φ is the opposite force of the convection of φ.

In the case

∂tφ + ∇ · (φu) = 0,

The force from the convection becomes

ρDu

Dt = −φ∇ · µ. We also get same result:

d dt Z 1 2ρ|u| 2_{+ f (φ, ∇φ) dx = 0.}

Dissipation One can add dissipation term in the momentum equation:

ρDu

Dt = µ∇φ + ∇ · τ where τ is the viscous stress.

The order parameter has a tendency to move from high potential energy to low potential energy. For the advection equation of φ, there are two models. The first one is Allan-Cahn, the second one is Cahn-Hilliard.

• Allan-Cahn:

∂tφ + u∇φ = −

δF δφ

(47)

• Cahn-Hilliard defines the flux to be

q = −m∇µ

where m is the mobility. The evolution of φ satisfies Dφ

Dt = ∇ · (m∇µ).

This is the Cahn-Hilliard equation. The total order parameter is conserved. The dissipation of energy can be obtained by multiplying momentum equation by u, multiplying advection equation by µ, then adding them together:

(ρDu Dt · u − µ∇φ · u) + (µ∂tφ + µu∇φ) = u · ∇ · τ − µ 2 . This gives dE dt = −τ · ∇u − µ 2_.

The stress tensor has the form ν(∇u + (∇u)T_{) + λ∇ · u. The viscous dissipation}

becomes

−ν|∇u + (∇u)T|2 − λ|∇ · u|2. For Cahn-Hilliard equation, the energy dissipation is

dE

dt = −τ · ∇u − m|∇µ|

2_.

6.3 Interface structure

In this section, we shall study the limit of various free energy. We are particularly interested in the kinetic energy part, which is related to the geometry of the interface.

F [φ] = Z f (∇φ, φ) dx 1. Let f (∇φ, φ) = 2|∇φ| 2₊ 1 W (φ) dx

where W is a double well potential. As → 0, φ → ±1. The interface energy tends to F [φ] → σ0|Γ(t)|. where σ0 is a constant. 2. Mean curvature: Z Ω(t) η 2 4φ − 1 η2W 0 (φ) 2 dx ≈ Z Γ(t) H2dS.

(48)

3. Gaussian curvature: Z Ω(t) η 4 φ − 1 ηW 0 (φ) W0(φ) dx ≈ Z Γ(t) K dS.

The bulk energy also has many choices.

1. Binary system: Consider a binary system with two components A and B with constration uA and uB, where uA+ uB = 1. The bulk energy is

W (uB) = µAuA+ µBuB+ RT (uAln uA+ uBln uB) + αuAuB.

Here, µA, µB are the chemical potential of components A and B, R the molar

gas constant, T the temperature, α the repulsion parameter between A and B. 2. Nernst-Plank-Poisson model: Consider binary charge system. Let p and n are respectively the concentration of positive and negative charges. The bulk energy is W (p, n) = Z Ω p log p + n log n + 1 2(n − p)V [n − p] dx where V [n − p] is the potential induced by n − p, i.e. V = n − p.

6.3.1 Structure of one dimensional interface

Allen-Cahn interface The order parameter satisfies

φt= − δF δφ = 4φ − W 0 (φ) where F = Z 1 2|∇φ| 2_{+ W (φ) dx} W (φ) = 1 4(φ 2_{− 1)}2 .

The interface is assumed to be steady, so φt= 0. Thus, we want to solve

4φ − W0(φ) = 0.

Multiplying φ on both sides, we get d dx 1 2φ 02_{− W (φ)} = 0.

We look for φ which connecting the two equilibria ±1 at x = ±∞ and with φ0(±∞) = 0. Thus, the interface we look for satisfies

1 2φ

(49)

We can solve this ODE: φ0 =√2W (6.1) By separation of variable dφ p2W (φ) = dx When W (φ) = 1₄(φ2_{− 1)}2_{, we get} 2dφ 1 − φ2 = ±dx

Integrate this, we get

x + C = Z ₁ 1 − φ + 1 1 + φdφ = ln 1 + φ 1 − φ Thus, let ξ = ex+C _{We have}

1 + φ 1 − φ = ±ξ. Or φ = ξ − 1 ξ + 1, or φ = − ξ + 1 −ξ + 1.

For the first solution, we have φ(±∞) = ±1, whereas the second solution satisfies φ(±∞) = ∓1. In general, φ can be expressed as A + B tanh(x + C).

Let us put the scale back. We consider the energy to be

F_[φ_{] =} Z 2|∇φ _|2₊1 W (φ _{) dx}

This gives the traveling wave equation

4 φ− 1 W

0

(φ) = 0.

Taking φ_{(x) = φ(x/), we get the rescaled equation. Notice that in F , the energy}

for the traveling wave solution φ is

F[φ] = F [φ] = σ0.

This is the interface energy in the normal direction of an interface. In multi-dimension, the interface energy is the integration of this value over the whole surface.

Remark. Let us consider a general double well potential W which has two stable equilibria φaand φb. Suppose W (φa) 6= W (φb). In this case, we don’t have a standing

interface. Instead, we have a traveling wave φ((x − ct)/), where c is the speed of the traveling wave. Plug this ansatz into equation δF /δφ = 0, we get

(50)

Cahn-Hilliard interface The Cahn-Hilliard equation is

φt= −∇ · (m∇µ), µ = −

δF φ . Again, we look for steady interface. This leads to

m(φ)∇µ = C, or

∇µ = C m(φ)

We assume for the moment that the mobility is independent of φ. Then we get µ = µ0 for all x.

Or

−φ00+ W0(φ) = µ0.

We look for solution with φ(±∞) = ±1 and φ0(±∞) = 0. Multiplying this equation by φ0, we get d dx −1 2φ 02 + W (φ) − µ0φ = 0.

Using the far field condition, we get

−1 2φ 02 + W (φ) − µ0φ = W (1) − µ0. Or φ02= 2(W (φ) − W (1) − µ0(φ − 1)) = F (φ) Thus, φ0 = ± (W (φ) − W (1) − µ0(φ − 1))1/2.

In order to have φ0(±∞) = 0, we see that F (1) = 0, we also need to require F (−1) = 0. This leads to W (−1) − W (1) − µ0(−1 − 1) = 0, which implies µ0 = (W (1) − W (−1))/2 = W (φa) − W (φb) φa− φb .

Here, φa and φb are the two equilibra of the double well potential W . Thus, the

only addimisible chemical potential connecting two equilibra is the relative energy difference between them. With this choice of chemical potential, we see the interface equation is the same as the Allen-Cahn interface equation. Thus, the structure is the same.

Cahn-Hilliard Traveling wave with nonzero speed Suppose the speed is c, then we look for traveling wave solution φ(x − ct). In this case, the

(51)

Combustion front In reaction-diffusion equation ut = 4u + f (u)

(52)

Appendix A

Stress-Strain relation

A.1 Constitutive relation: Stress-strain relation

The stress in an elastic material is due to deformation. Therefore we assume the stress tensor T is a function of the deformation gradient F . Such a relation

T (y) = bT (x, F )

is called a constitutive relation and bT is called a response function. If bT is indepen-dent of x, the material is called homogeneous. Below, we want to find the physical conditions of a response function.

A.1.1 Frame indifference

First, the stress should be independent of the frame of reference physically. Suppose y is the coordinate in the original frame, and y∗ is the coordinate in new frame which is a rotation of the old frame. That is, y∗ = Qy, and Q is a rotation. Then F∗ := ∂y∗/∂x = QF ; the normal and the stress in the old and new frames are related by ν∗ = Qν, σ∗ = Qσ. In terms of the response function, it reads

σ∗ = T (QF )νb ∗ = bT (QF )Qν Qσ = Q bT (F )ν.

Hence,

b

T (QF ) = Q bT (F )QT. (1.1) Such a condition for bT is called the axiom of material frame-indifference.

Next, we will show that bT is indeed a function of the Cauchy-Green strain tensor only. Let us write F in polar form

F = RU.

where R is a rotation and U is positive definite. From the axiom of material frame-indifference, we get

b

(53)

where C = U2 = FTF is the right Cauchy-Green strain tensor. Similarly, for the Piola tensor P , the frame-indifference reads

P (QF ) = J bT (QF )(QF )−T = J (Q bT (F )QT)Q−TF−T = QP (F ). and

P (F ) = J bT (F )F−T = J F T (C) = F P (C),

where J2 _{= detC and P (C) :=}√_{detC T (C). For the second Piola tensor it satisfies}

Σ(QF ) = (QF )−1P (QF ) = (F−1Q−1)(QP (F )) = Σ(F ), and

Σ(F ) := F−1P = P (C).

A.1.2 Isotropic material

A material is isotropic if its stress tensor is independent of a rotation of the strain. Suppose we have two identical materials. One has material coordinate x, while the other is x∗ with x∗ = Qx and Q being a rotation. If these two materials deform in the same way and result the same stresses, we call this material isotropic. Mathematically, this means

b

T (F Q) = bT (F ). (1.2) In terms of Piola stress tensor, it reads

P (F Q) = P (F )Q for all Q (1.3) This follows from

P (F Q) = J bT (F Q)(F Q)−T = J bT (F )F−TQ = P (F )Q.

In order to characterize the response function of an isotropic material, we review some spectral properties of a matrix.

• Given a matrix A, the characteric polynomial is

p(λ) := det(A − λI) = −λ3+ ι1λ2− ι2λ + ι3.

The coefficients ι1, ι2, ι3 are called the principal invariants of A.

• ι1 = trA = λ1+ λ2+ λ3 ι2 = 1 2 (trA) 2_{− tr(A}2_{) = tr(Cof A) = λ} 1λ2+ λ2λ3+ λ3λ1 ι3 = detA = λ1λ2λ3.

(54)

• The matrix Ap

, with p ∈ Z and (p ≥ 0 or p ≤ 1 if A is invertible) has the representation:

Ap = α0,p(ιA)I + α1,p(ιA)A + α2,p(ιA)A2.

where αk,p are functions of principal invariants.

• spectral mapping theorem: If B is a self-adjoint operator with eigenvalues/eigenvectors λi and pi, i = 1, 2, 3. Let β(λ) be a smooth function, then

β(B) =

3

X

i=1

β(λi)pipTi

Theorem 1.1 (Rivlin-Ericksen representation theorem) A mapping b_{T : M}3₊→ S3 satisfies

b

T (QF ) = Q bT (F )QT and bT (F Q) = b_{T (F ) for all F ∈ M}3₊_{, Q ∈ Q}3₊, (1.4)

if and only if bT (F ) = T (F FT_{) with}

T (B) = β0(ιB)I0 + β1(ιB)B + β2(ιB)B2.

Proof.

1. bT (F ) = bT (V R) = bT (V ) = T (V2_{) = T (B) = T (F F}T_).

2. T (QBQT_{) = QT (B)Q}T_{. This follows from}

T (QBQT) = T (QV (QV )T) = T (QV ) = QT (V )QT = QT (B)QT. This means that T (B) is only a function of eigenvalues of B.

3. For positive definite matrix B, the function T (B) has the representation: T (B) = b0(B)I + b1(B)B + b2(B)B2

where bi(B) are scalar functions of B. The proof of this claim is divided into

three cases: (i) λi are distinct, (ii) λ2 = λ3, (iii) all eigenvalues are equal.

(i) We have I = p1pT1 + p2pT2 + p3pT3 B = λ1p1pT1 + λ2p2pT2 + λ3p3pT3 B2 = λ2₁p1pT1 + λ 2 2p2pT2 + λ 2 3p3pT3

and the spectral mapping

T (B) = µ1p1pT1 + µ2p2pT2 + µ3p3pT3.

(55)

(ii) We have

I = p1pT1 + (p2pT2 + p3pT3)

B = λ1p1pT1 + λ2(p2pT2 + p3pT3)

and the spectral mapping

T (B) = µ1p1pT1 + µ2p2pT2 + µ3p3pT3.

From bT (QBQT_{) = QT Q}T_{, we get µ}

2 = µ3. This leads to

T (B) = b0(B)I + b1(B)B.

(iii) In this case, we have bT (B) = b0(B)I.

4. It remains to show that bi(B) is indeed only functions of its invariants. From

the independence of I, B, B2 and T (QBQT) = QT (B)QT, we get bi(QBQT) =

bi(B). This leads to bi are only functions of the three principal invariants.

Theorem 1.2 For an isotropic material, its second Piola stress-strain relation has the following representation

Σ(C) = γ0(ιC)I + γ1(ιC)C + γ2(ιC)C2.

Proof. (i) From the polar decomposition: F = RU = V R and C = FT_{F = U}2_{, B =}

F FT = V2, we get C = RTBR. Thus, ιC = ιB, F−1BF−T = I and F−1B2F−T = C.

(ii) Next, from Σ = J F−1T F−T and Rivlin-Ericksen theorem, we get Σ = J F−1 β0(ιB)I + β1(ιB)B + β2(ιB)B2 F−T

= J β0(ιC)C−1+ β1(ιC)I + β2(ιC)C .

(iii) From Caley-Hamilton theorem, it holds

C3 − ι1(C)C2+ ι2(C)C − ι3(C)I = 0.

Hence,

C−1 = ι−1₃ ι2I − ι1C + C2 .

(iv) It holds J = ι3(C)1/2.

From (ii), (iii) and (iv), we complete the proof of the theorem.

Theorem 1.3 For an isotropic elastic material, by normalizing Σ(I) = 0, the Piola stress tensor has the following representation:

Σ(C) = λ(tr E)I + 2µE + o(E). (1.5) Here, λ and µ are called the Lam´e constants, E = (C − 1)/2.

(56)

Proof. (i) At equilibrium, y = x and F = I, C = I, B = I, the stress is normalized to be bT (I) = 0 and Σ(I) = 0. In this case, we have

2 X i=0 βi(ιI) = 0 2 X i=0 γi(ιI) = 0.

(ii) We expand the principal invariants about E = 0: First, we have tr C = 3 + 2tr E

tr C2 = 3 + 4tr E + o(E) tr C3 = 3 + 6tr E + o(E). Second, the principal invariants can be expanded as

ι1(C) = tr C ι2(C) = 1 2 (tr C) 2_{− tr C}2_{= 3 + 4tr E + o(E)} ι3(C) = 1 6(tr C) 3₋ 1 2tr C tr C 2₊1 3tr C 3 = 1 + 2tr E + o(E).

(iii) We expand (1.5) in E = (C − I)/2 around E = 0 and use (i) and (ii).

A.1.3 Hyper-elastic material

A material is called hyper-elastic material if there exists a potential W (F ) such that

P (F ) = ∂W ∂F .

Notice that such materials satisfy the frame-indifference axiom if and only if W (QF ) = W (F ) for all Q.

To see this, we differentiate the above equation in F to get

P (F ) = ∂ ∂FW (F ) = ∂ ∂FW (QF ) = Q T∂W (QF ) ∂(QF ) = Q T_{P (QF ).}

This shows W (F ) = W (QF ) if and only if P (QF ) = QP (F ) for all Q. Similarly, a hyper-elastic material is isotropic (i.e. P (F Q) = P (F )Q if and only if W (F Q) = W (F ) for all Q. P (F ) = ∂ ∂FW (F ) = ∂ ∂FW (F Q) = ∂W (QF ) ∂(QF ) Q T _{= P (QF )Q}T_.

(57)

Thus, for an isotropic hyper-elastic material, its potential energy function W satisfies

W (F Q) = W (QF ) = W (F ) for all Q (1.6) Thus, W is only a function of the principal invariants of F . Two simple examples are

• Hook’s law: W (F ) = 1 2H|F | 2_{, where |F |}2 ₌P ij|Fij|2. • Neo-Hookean material: W (F ) = a|F |2+ g(det F ), a > 0. • Mooney-Rivlin materials:

W (F ) = a|F |2+ b|adjF |2+ g(det F ), a, b > 0.

•

W = λ 2(tr E)

2_{+ µ tr (E}2_{), where E =} 1

2(C − I) (1.7)

A.1.4 Linear elastic material

When there is no deformation, i.e. F = I,C = I, it is natural to set the corresponding stress bT (I) = 0. In terms of Piola stress tensor, Σ(I) = P (I) = 0. When |∇xu| =

<< 1, C = FTF = I + ∇xu + (∇xu)T + (∇xu)T(∇u) = I + 2E + o(), where E = 1 2 ∇xu + (∇xu) T_.

We can expand P about I in terms of E, or (C − I):

P (C) = P (I) + 1

2AE + o(|E|).

where A = (aijkl) is a 4-tensor. Then the Piola tensor P has the form:

P = F P (C) = (I + ∇u)AE + o() = AE + o().

Here, we have absorbed ∇uAE into o(). A material whose Piola tensor has the form P = AE is called a linear elastic material. In coordinate form:

Pij = aijklekl.

The four-tensor A should satisfy

Continuum Mechanics