A Material Point Method Formulation for Plasticity

Biswajit Banerjee

November 14, 2006

1 Governing Equations
The equations that govern the motion of an elastic-plalstic solid include the balance laws for mass, momentum, and
energy. Kinematic equations and constitutive relations are needed to complete the system of equations. Physical
restrictions on the form of the constitutive relations are imposed by an entropy inequality that expresses the second
law of thermodynamics in mathematical form.
The balance laws express the idea that the rate of change of a quantity (mass, momentum, energy) in a volume must
arise from three causes:

1. the physical quantity itself flows through the surface that bounds the volume,
2. there is a source of the physical quantity on the surface of the volume, or/and,
3. there is a source of the physical quantity inside the volume.

Let Ω be the body (an open subset of Euclidean space) and let ∂Ω be its surface (the boundary of Ω).
Let the motion of material points in the body be described by the map

x = ϕ(X) = x(X) (1)

where X is the position of a point in the initial configuration and x is the location of the same point in the deformed
configuration. The deformation gradient (F ) is given by
∂x
F = = ∇0 x . (2)
∂X

1.1 Balance Laws

Let f (x, t) be a physical quantity that is flowing through the body. Let g(x, t) be sources on the surface of the body
and let h(x, t) be sources inside the body. Let n(x, t) be the outward unit normal to the surface ∂Ω. Let v(x, t) be
the velocity of the physical particles that carry the physical quantity that is flowing. Also, let the speed at which the
bounding surface ∂Ω is moving be un (in the direction n).
Then, balance laws can be expressed in the general form ([1])
Z  Z Z Z
d
f (x, t) dV = f (x, t)[un (x, t) − v(x, t) · n(x, t)] dA + g(x, t) dA + h(x, t) dV . (3)
dt Ω ∂Ω ∂Ω Ω
1
Note that the functions f (x, t), g(x, t), and h(x, t) can be scalar valued, vector valued, or tensor valued -
depending on the physical quantity that the balance equation deals with.
It can be shown that the balance laws of mass, momentum, and energy can be written as (see Appendix):

ρ̇ + ρ ∇ · v = 0 Balance of Mass
ρ v̇ − ∇ · σ − ρ b = 0 Balance of Linear Momentum
(4)
σ = σT Balance of Angular Momentum
ρ ė − σ : (∇v) + ∇ · q − ρ s = 0 Balance of Energy.

In the above equations ρ(x, t) is the mass density (current), ρ̇ is the material time derivative of ρ, v(x, t) is the
particle velocity, v̇ is the material time derivative of v, σ(x, t) is the Cauchy stress tensor, b(x, t) is the body force
density, e(x, t) is the internal energy per unit mass, ė is the material time derivative of e, q(x, t) is the heat flux
vector, and s(x, t) is an energy source per unit mass.
With respect to the reference configuration, the balance laws can be written as

ρ det(F ) − ρ0 = 0 Balance of Mass

T
ρ0 ẍ − ∇0 · P − ρ0 b = 0 Balance of Linear Momentum
T T
(5)
F ·P =P ·F Balance of Angular Momentum
ρ0 ė − P T : Ḟ + ∇0 · q − ρ0 s = 0 Balance of Energy.

In the above, P is the first Piola-Kirchhoff stress tensor, and ρ0 is the mass density in the reference configuration.
The first Piola-Kirchhoff stress tensor is related to the Cauchy stress tensor by

P = det(F ) F −1 · σ . (6)
We can alternatively define the nominal stress tensor N which is the transpose of the first Piola-Kirchhoff stress tensor such that

N := P T = det(F ) (F −1 · σ)T = det(F ) σ · F −T . (7)

Then the balance laws become

ρ det(F ) − ρ0 = 0 Balance of Mass

ρ0 ẍ − ∇0 · N − ρ0 b = 0 Balance of Linear Momentum
(8)
F · NT = N · F T Balance of Angular Momentum
ρ0 ė − N : Ḟ + ∇0 · q − ρ0 s = 0 Balance of Energy.

The gradient and divergence operators are defined such that

3 3 3
X ∂vi X ∂vi X ∂Sij
∇v = ei ⊗ ej = vi,j ei ⊗ ej ; ∇ · v = = vi,i ; ∇ · S = ei = σij,j ei . (9)
∂xj ∂xi ∂xj
i,j=1 i=1 i,j=1

where v is a vector field, BS is a second-order tensor field, and ei are the components of an orthonormal basis in
the current configuration. Also,
3 3 3
X ∂vi X ∂vi X ∂Sij
∇0 v = Ei ⊗ Ej = vi,j Ei ⊗ Ej ; ∇0 · v = = vi,i ; ∇0 · S = Ei = Sij,j Ei (10)
∂Xj ∂Xi ∂Xj
i,j=1 i=1 i,j=1

2
where v is a vector field, BS is a second-order tensor field, and Ei are the components of an orthonormal basis in
the reference configuration.
The contraction operation is defined as
3
X
A:B= Aij Bij = Aij Bij . (11)
i,j=1

1.2 The Clausius-Duhem Inequality

The Clausius-Duhem inequality can be used to express the second law of thermodynamics for elastic-plastic
materials. This inequality is a statement concerning the irreversibility of natural processes, especially when energy
dissipation is involved.
Just like in the balance laws in the previous section, we assume that there is a flux of a quantity, a source of the
quantity, and an internal density of the quantity per unit mass. The quantity of interest in this case is the entropy.
Thus, we assume that there is an entropy flux, an entropy source, and an internal entropy density per unit mass (η)
in the region of interest.
Let Ω be such a region and let ∂Ω be its boundary. Then the second law of thermodynamics states that the rate of
increase of η in this region is greater than or equal to the sum of that supplied to Ω (as a flux or from internal
sources) and the change of the internal entropy density due to material flowing in and out of the region.
Let ∂Ω move with a velocity un and let particles inside Ω have velocities v. Let n be the unit outward normal to
the surface ∂Ω. Let ρ be the density of matter in the region, q̄ be the entropy flux at the surface, and r be the
entropy source per unit mass. Then (see [1, 2]), the entropy inequality may be written as
Z  Z Z Z
d
ρ η dV ≥ ρ η (un − v · n) dA + q̄ dA + ρ r dV . (12)
dt Ω ∂Ω ∂Ω Ω
The scalar entropy flux can be related to the vector flux at the surface by the relation q̄ = −ψ(x) · n. Under the
assumption of incrementally isothermal conditions (see [3] for a detailed discussion of the assumptions involved),
we have
q(x) s
ψ(x) = ; r=
T T
where q is the heat flux vector, s is a energy source per unit mass, and T is the absolute temperature of a material
point at x at time t.
We then have the Clausius-Duhem inequality in integral form:
q·n ρs
Z  Z Z Z
d
ρ η dV ≥ ρ η (un − v · n) dA − dA + dV . (13)
dt Ω ∂Ω ∂Ω T Ω T

We can show that (see Appendix) the entropy inequality may be written in differential form as
!
q ρs
ρ η̇ ≥ −∇ · + . (14)
T T

In terms of the Cauchy stress and the internal energy, the Clausius-Duhem inequality may be written as (see
Appendix)
q · ∇T
ρ (ė − T η̇) − σ : ∇v ≤ − . (15)
T
3
1.3 Constitutive Relations

A set of constitutive equations is required to close to system of balance laws. For large deformation plasticity, we
have to define appropriate kinematic quantities and stress measures so that constitutive relations between them may
have a physical meaning.

1.3.1 Thermoelastic Relations

Let the fundamental kinematic quantity be the deformation gradient (F ) which is given by
∂x
F = = ∇0 x ; det F > 0 .
∂X
A thermoelastic material is one in which the internal energy (e) is a function only of F and the specific entropy (η),
that is
e = ē(F , η) .
For a thermoelastic material, we can show that the entropy inequality can be written as (see Appendix)

q · ∇T
   
∂ē ∂ē
ρ −T η̇ + ρ − σ · F −T : Ḟ + ≤0. (16)
∂η ∂F T

At this stage, we make the following constitutive assumptions:

1. Like the internal energy, we assume that σ and T are also functions only of F and η, i.e.,

σ = σ(F , η) ; T = T (F , η) .

2. The heat flux q satisfies the thermal conductivity inequality and if q is independent of η̇ and Ḟ , we have

q · ∇T ≤ 0 =⇒ −(κ · ∇T ) · ∇T ≤ 0 =⇒ κ≥0

i.e., the thermal conductivity κ is positive semidefinite.

Therefore, the entropy inequality may be written as

   
∂ē ∂ē −T
ρ − T η̇ + ρ −σ·F : Ḟ ≤ 0 .
∂η ∂F

Since η̇ and Ḟ are arbitrary, the entropy inequality will be satisfied if and only if
∂ē ∂ē ∂ē ∂ē
− T = 0 =⇒ T = and ρ − σ · F −T = 0 =⇒ σ = ρ · FT .
∂η ∂η ∂F ∂F
Therefore,
∂ē ∂ē
T = and σ=ρ · FT . (17)
∂η ∂F
Given the above relations, the energy equation may expressed in terms of the specific entropy as (see Appendix)

ρ T η̇ = −∇ · q + ρ s . (18)
4
 • Effect of a rigid body rotation of the internal energy:
If a thermoelastic body is subjected to a rigid body rotation Q, then its internal energy should not change. After a rotation, the new deformation
gradient (F̂ ) is given by
F̂ = Q · F .
Since the internal energy does not change, we must have
e = ē(F̂ , η) = ē(F , η) .
Now, from the polar decomposition theorem, F = R · U where R is the orthogonal rotation tensor (i.e., R · RT = RT · R = 1 ) and U is the
symmetric right stretch tensor. Therefore,
ē(Q · R · U , η) = ē(F , η) .
Now, we can choose any rotation Q. In particular, if we choose Q = RT , we have

ē(RT · R · U , η) = ē(1 · U , η) = ẽ(U , η) .

Therefore,
ē(U , η) = ē(F , η) .
This means that the internal energy depends only on the stretch U and not on the orientation of the body.

1.3.2 Alternative strain and stress measures

The internal energy depends on F only through the stretch U . A strain measure that reflects this fact and also
vanishes in the reference configuration is the Green strain

1 1
E = (F T · F − 1 ) = (U 2 − 1 ) . (19)
2 2

Recall that the Cauchy stress is given by

∂ē
· FT . σ=ρ
∂F
We can show that the Cauchy stress can be expressed in terms of the Green strain as (see Appendix)

∂ē
σ=ρF · · FT . (20)
∂E

Recall that the nominal stress tensor is defined as

N = det F (σ · F −T ) .

From the conservation of mass, we have ρ0 = ρ det F . Hence,

ρ0
N= σ · F −T . (21)
ρ

The nominal stress is unsymmetric. We can define a symmetric stress measure with respect to the reference
configuration call the second Piola-Kirchhoff stress tensor (S):

ρ0 −1
S := F −1 · N = P · F −T = F · σ · F −T . (22)
ρ

In terms of the derivatives of the internal energy, we have

ρ0 −1
 
∂ē ∂ē
S= F · ρF · ·F T
· F −T = ρ0
ρ ∂E ∂E
5
and
ρ0
 
∂ē ∂ē
N= ρF · · FT · F −T = ρ0 F · .
ρ ∂E ∂E
That is,
∂ē ∂ē
S = ρ0 and N = ρ0 F · . (23)
∂E ∂E

1.3.3 Stress Power

The stress power per unit volume is given by σ : ∇v. In terms of the stress measures in the reference
configuration, we have  
∂ē
σ : ∇v = ρ F · ·F T
: (Ḟ · F −1 ) .
∂E
Using the identity A : (B · C) = (A · C T ) : B, we have

ρ
    
∂ē T −T ∂ē
σ : ∇v = ρF · ·F ·F : Ḟ = ρ F · : Ḟ = N : Ḟ .
∂E ∂E ρ0

We can alternatively express the stress power in terms of S and Ė. Taking the material time derivative of E we
have
1
Ė = (F˙T · F + F T · Ḟ ) .
2
Therefore,
1
S : Ė = [S : (F˙T · F ) + S : (F T · Ḟ )] .
2
Using the identities A : (B · C) = (A · C T ) : B = (B T · A) : C and A : B = AT : B T and using the symmetry
of S, we have
1 1
S : Ė = [(S · F T ) : Ḟ T + (F · S) : Ḟ ] = [(F · S T ) : Ḟ + (F · S) : Ḟ ] = (F · S) : Ḟ .
2 2

Now, S = F −1 · N . Therefore, S : Ė = N : Ḟ . Hence, the stress power can be expressed as

σ : ∇v = N : Ḟ = S : Ė . (24)

If we split the velocity gradient into symmetric and skew parts using

∇v = l = d + w

where d is the rate of deformation tensor and w is the spin tensor, we have

σ : ∇v = σ : d + σ : w = tr(σ T · d) + tr(σ T · w) = tr(σ · d) + tr(σ · w) .

Since σ is symmetric and w is skew, we have tr(σ · w) = 0. Therefore, σ : ∇v = tr(σ · d). Hence, we may also
express the stress power as
tr(σ · d) = tr(N T · Ḟ ) = tr(S · Ė) . (25)

6
1.3.4 Helmholtz and Gibbs free energy

Recall that
∂ē
S = ρ0 .
∂E
Therefore,
∂ē 1
= S.
∂E ρ0
Also recall that
∂ē
=T .
∂η
Now, the internal energy e = ē(E, η) is a function only of the Green strain and the specific entropy. Let us assume,
that the above relations can be uniquely inverted locally at a material point so that we have

E = Ẽ(S, T ) and η = η̃(S, T ) .

Then the specific internal energy, the specific entropy, and the stress can also be expressed as functions of S and T ,
or E and T , i.e.,

e = ē(E, η) = ẽ(S, T ) = ê(E, T ) ; η = η̃(S, T ) = η̂(E, T ) ; and S = Ŝ(E, T )

We can show that (see Appendix)

d 1 dψ 1
(e − T η) = −Ṫ η + S : Ė or = −Ṫ η + S : Ė . (26)
dt ρ0 dt ρ0
and
d 1 1 dg 1
(e − T η − S : E) = −Ṫ η − Ṡ : E or = Ṫ η + Ṡ : E . (27)
dt ρ0 ρ0 dt ρ0
We define the Helmholtz free energy as

ψ = ψ̂(E, T ) := e − T η . (28)

We define the Gibbs free energy as

1
g = g̃(S, T ) := −e + T η + S:E. (29)
ρ0

The functions ψ̂(E, T ) and g̃(S, T ) are unique. Using these definitions it can be showed that (see Appendix)

∂ ψ̂ 1 ∂ ψ̂ ∂g̃ 1 ∂g̃
= Ŝ(E, T ) ; = −η̂(E, T ) ; = Ẽ(S, T ) ; = η̃(S, T ) (30)
∂E ρ0 ∂T ∂S ρ0 ∂T
and
∂ Ŝ ∂ η̂ ∂ Ẽ ∂ η̃
= −ρ0 and = ρ0 . (31)
∂T ∂E ∂T ∂S

7
1.3.5 Specific Heats

The specific heat at constant strain (or constant volume) is defined as

∂ê(E, T )
Cv := . (32)
∂T
The specific heat at constant stress (or constant pressure) is defined as

∂ẽ(S, T )
Cp := . (33)
∂T
We can show that (see Appendix)
∂ η̂ ∂ 2 ψ̂
Cv = T = −T (34)
∂T ∂T 2
and
∂ η̃ 1 ∂ Ẽ ∂ 2 g̃ ∂ 2 g̃
Cp = T + S: =T + S : . (35)
∂T ρ0 ∂T ∂T 2 ∂S∂T
Also the equation for the balance of energy can be expressed in terms of the specific heats as (see Appendix)

ρ
ρ Cv Ṫ = ∇ · (κ · ∇T ) + ρ s + T βS : Ė
ρ0
! (36)
1 ρ
ρ Cp − S : αE Ṫ = ∇ · (κ · ∇T ) + ρ s − T αE : Ṡ
ρ0 ρ0

where
∂ Ŝ ∂ Ẽ
βS := and αE := . (37)
∂T ∂T
The quantity βS is called the coefficient of thermal stress and the quantity αE is called the coefficient of thermal
expansion.
The difference between Cp and Cv can be expressed as
!
1 ∂ Ŝ ∂ Ẽ
Cp − Cv = S−T : . (38)
ρ0 ∂T ∂T

However, it is more common to express the above relation in terms of the elastic modulus tensor as (see Appendix
for proof)
1 T
Cp − Cv = S : αE + αE : C : αE (39)
ρ0 ρ0
where the fourth-order tensor of elastic moduli is defined as
∂ Ŝ ∂ 2 ψ̂
C := = ρ0 . (40)
∂ Ẽ ∂ Ẽ∂ Ẽ
For isotropic materials with a constant coefficient of thermal expansion that follow the St. Venant-Kirchhoff
material model, we can show that (see Appendix)
1
α tr(S) + 9 α2 K T .

Cp − Cv =
ρ0
8
1.3.6 Kinematic Relations

Additive split of the rate of deformation tensor: The rate of deformation tensor (d) is given by

1
∇v + (∇v)T

d= (41)
2

where v is the velocity. We assume that the rate of deformation can be additively decomposed into an elastic part, a
plastic part, and a thermal part:
d = de + dp + dth . (42)
The thermal part of the rate of deformation is computed separately and subtracted from d to give

de := d − dth = ee + dp .

For simplicity, we use the symbol d instead of de in the following development.

We split the rate of deformation tensors into volumetric and deviatoric parts:
1 1 1
d= θ̇ 1 + η ; de = tr(de ) 1 + η e ; dp = tr(dp ) 1 + η p (43)
3 3 3
1
where θ̇ = tr(d), the deviatoric part of the rate of deformation tensor is η = d − 3 tr(d) 1 , and 1 is the
second-order identity tensor.

1.3.7 Stress

To deal with nearly incompressible behavior, we introduce an isomorphic split of the Cauchy stress tensor into a
volumetric and a deviatoric part of the form ([4]):

σ = σm 1
b + σs sb ⇐⇒ σ =p1 +s (44)

where
1 √ 1 1 1 s √
σm = √ tr(σ) = 3 p ; σs = ksk ; p = tr(σ) ; s = σ− tr(σ) 1 ; 1
b= ; sb = ; kAk = A : A .
3 3 3 k1 k ksk
(45)
Note that the following identities hold:

σm = σ : 1
b ; σs = σ : sb ; 1
b :1
b = 1 ; sb : sb = 1 ; 1 b : sb˙ = 0 ; sb : sb˙ = 0 .
b : sb = 0 ; 1 (46)

The rate form of equation (44) is

b + σ̇s sb + σs sb˙
σ̇ = σ̇m 1 ⇐⇒ σ̇ = ṗ 1 + ṡ . (47)

In the following development, we enforce frame indifference while evaluating the constitutive relation by assuming
that the stresses and rate of deformation have been rotated to the reference configuration using

s = RT · s · R ; η = RT · η · R

where the rotation tensor R is obtained from a polar decomposition of the deformation gradient.
9
1.3.8 Elastic Relations:

The constitutive model is assumed to consist of two parts. The first is a hypoelastic model for the deviatoric part of
the stress and the second is a equation of state for the mean stress.
The deviatoric stress is given by the rate equation

ṡ = C(σm , T ) : η e ⇐⇒ b T ) : ηe
ṡ = C(p, (48)

where ṡ is the deviatoric stress rate, C is a pressure and temperature dependent fourth order elastic modulus tensor,
and η e is the deviatoric part of de (the elastic part of the rate of deformation tensor).
Let us assume that the deviatoric elastic response of the material is isotropic, i.e.,

ṡ = 2 µ(σm , T ) η e ⇐⇒ b(p, T ) η e .
ṡ = 2 µ (49)

Remark: Note that since the shear modulus depends on σm , the rate of s should have a contribution that contains
σ̇m . We ignore this contribution in the subsequent development.
The mean stress (σm ) is calculated using a Mie-Grüneisen equation of state of the form ([5, 6])
√
√ 3 ρ0 C02 (1 − J e )[1 − Γ0 (1 − J e )] √
σm = 3p=− − 2 3 Γ0 e ; J e := det F e (50)
[1 − Sα (1 − J e )]2

where C0 is the bulk speed of sound, ρ0 is the initial mass density, 2 Γ0 is the Grüneisen’s gamma at the reference
state, Sα = dUs /dUp is a linear Hugoniot slope coefficient, Us is the shock wave velocity, Up is the particle
velocity, and e is the internal energy density (per unit reference volume), F e is the elastic part of the deformation
gradient. For isochoric plasticity,
ρ0
J e = J = det(F ) = .
ρ

The change in internal energy is computed using

Z
e = ρ0 Cv dT ≈ ρ0 [Cv (T ) T − Cv (T0 ) T0 ] (51)

where T0 is the reference temperature and Cv is the specific heat at constant volume. We assume that Cp and Cv are
equal.
The rate equation (47) has three rate terms. We need constitutive relations for each of these three terms.
The rate form of the pressure equation is obtained by taking the material time derivative of (50) to get
√ !
√ 3 ρ0 C02 J e [1 + (Sα − 2 Γ0 )(1 − J e )] J˙e √
σ̇m = 3 ṗ = e 3 e
− 2 3 Γ0 ė (52)
[1 − Sα (1 − J )] J

Now,
J˙e √
e
= tr(de ) = θ˙e = de : 1 = 3 de : 1
b ; ė = ρ0 Cv (T ) Ṫ .
J
For isochoric plasticity, the above relations become

J˙ √
= tr(de ) = tr(d) = θ̇ = d : 1 = 3 d : 1
b.
J
10
Hence,
√ √
σ̇m = 3 ṗ = 3 K(J e ) de : 1
b − 2 3 ρ0 Cv Γ0 Ṫ (53)
where
ρ0 C02 J e [1 + (Sα − 2 Γ0 )(1 − J e )]
K(J e ) = . (54)
[1 − Sα (1 − J e )]3
The stress rate σ̇s is given by

∂ √
σ̇s = ( s : s) = ṡ : sb = 2 µ(σm , T ) η e : sb .
∂t
Therefore,
σ̇s = 2 µ(σm , T ) η e : sb . (55)
The rate of sb is obtained by noting that !
ṡ ṡ : sb
sb˙ = − sb .
ksk ksk
Therefore,
2 µ(σm , T ) e
sb˙ = [η − (η e : sb)b
s] . (56)
σs
A slightly different, but equally valid, decomposition of the stress rate can be obtained by setting

σ̇s = 2 µ(σm , T ) de : sb . (57)

and
2 µ(σm , T ) h e i
sb˙ = d − (de : 1 b − (de : sb)b
b )1 s . (58)
σs
We use the above decomposition in the following development.

1.3.9 Plastic Relations:

We consider yield functions of the form

f (σm , σs , ˙, p , T ) ≤ 0 (59)
where ˙ is the total strain rate, p is the plastic strain, and T is the temperature. Here the strain rates and the plastic
strain are defined as s s
Z t
2 2 p p 0
˙ := d : d ; ˙p := d : d ; p := ˙p dt . (60)
3 3 0
As a particular case, we consider a class of isotropic yield functions of the von Mises form

3 2 µ2 (σm , T )
f (σm , σs , ˙, p , T ) = σs − σy2 (˙, p , T ) (61)
2 µ20

where σy is the yield stress, µ is the shear modulus, and µ0 is a reference shear modulus.

11
Remark: This form of the yield function can be used for deformation induced anisotropic plasticity by converting
it to the form used by Maudlin and Schiferl ([7]). In addition, the effect of porosity can be incorporated using
Brannon’s approach ([4]).
At yield, we have
3 2 µ2 (σm , T )
σs − σy2 (˙, p , T ) =0.
2 µ20
Taking a time derivative of the above equation gives us

µ2 ∂σy µ ∂µ
3 σs σ̇s − 2 σy 2 − 2 σy2 2 σ̇m = 0 . (62)
µ0 ∂t µ0 ∂σm

We assume that the plastic part of the rate of deformation can be obtained from a plastic potential in the form
(associated flow rule)
∂f
dp = λ̇ = λ̇ fσ (63)
∂σ
The unit outward normal to the yield surface is given by

1 ∂f 1 ∂f
n
b= = fσ ; ξ = k k = kfσ k ; n
b:n
b=1 =⇒ dp = λ̇ ξ n
b. (64)
ξ ∂σ ξ ∂σ
Therefore, s s
2 2
Z t
0
˙p := λ̇ ξ ; p := λ̇ ξ dt .
3 3 0
Note that
∂σm b ; ∂σs = sb .
=1 (65)
∂σ ∂σ
Assuming that ,
˙ p , and T remain unchanged during the elastic part of the deformation, we have
∂f ∂f b ∂f p
b + fs sb ; ξ = (fm )2 + (fs )2 .
fσ = = 1+ sb = fm 1 (66)
∂σ ∂σm ∂σs
For the yield function (61), we have

µ(σm , T ) ∂µ
fm = −2 σy2 (˙, p , T ) ; fs = 3 σs . (67)
µ20 ∂σm

Therefore, equation (62) can be written as

µ2 ∂σy
fs σ̇s − 2 σy + fm σ̇m = 0 . (68)
µ20 ∂t
Also,
" #
µ(σ m , T ) ∂µ
dp = λ̇ −2 σy2 (˙, p , T ) 1
b + 3 λ̇ σs sb . (69)
µ20 ∂σm

This implies that

µ(σm , T ) ∂µ
dp = (dp : 1 b + (dp : sb)b
b )1 s ; dp : 1
b = −2 λ̇ σy2 (˙, p , T ) ; dp : sb = 3 λ̇ σs ;
µ20 ∂σm
12
Recall,

b + σ̇s sb + σs sb˙
σ̇ = σ̇m 1
√
σ̇m = 3 K(J e ) de : 1b − 2 3 ρ0 Cv Γ0 Ṫ
σ̇s = 2 µ(σm , T ) de : sb
2 µ(σm , T ) h e i
sb˙ = d − (de : 1 b − (de : sb)b
b )1 s
σs
Using the decomposition de = d − dp , we get
√
σ̇m = 3 K(J e ) (d − dp ) : 1
b − 2 3 ρ0 Cv Γ0 Ṫ
σ̇s = 2 µ(σm , T ) (d − dp ) : sb
2 µ(σm , T ) h i
sb˙ = d − dp − (d : 1 b + (dp : 1
b )1 s + (dp : sb)b
b − (d : sb)b
b )1 s
σs
or,
b + σ̇s sb + σs sb˙
σ̇ = σ̇m 1
" #
µ(σ m , T ) ∂µ √
σ̇m = 3 K(J e ) d : 1 b + 2 λ̇ σy2 (˙, p , T )
2 − 2 3 ρ0 Cv Γ0 Ṫ
µ0 ∂σm
(70)
σ̇s = 2 µ(σm , T ) (d : sb − 3 λ̇ σs ) = 2 µ(σm , T ) (η : sb − 3 λ̇ σs )
2 µ(σm , T ) 2 µ(σm , T )
sb˙ = [η − (d : sb)b
s] = [η − (η : sb)b
s]
σs σs
Now, the consistency condition requires that

∂
[f (σm , σs , ˙, p , T )] = 0
∂t
or,
∂f ∂f ∂f ∂ ˙ ∂f ∂f
σ̇m + σ̇s + + ˙ + Ṫ = 0
∂σm ∂σs ∂ ˙ ∂t ∂p ∂T
or,
fm σ̇m + fs σ̇s + fd ¨ + fp ˙ + ft Ṫ = 0 (71)
where
s
∂f µ2 ∂σy 2 d˙ : d 1
fd := = 2 σy 2 ; ¨ = ; d˙ = RT · [∇a + (∇a)T ] · R
∂ ˙ µ0 ∂ ˙ 3 kdk 2
∂f µ2 ∂σy ∂f µ2 ∂σy
fp := = 2 σy 2 ; ft := = 2 σy 2 .
∂p µ0 ∂p ∂T µ0 ∂T

Note that the above equation is the same as equation (68). Also note that,

b − λ̇ fm ] − G(T ) Ṫ ; σ̇s = 2 µ [η : sb − λ̇ fs ]
σ̇m = 3 K [d : 1

where √
G(T ) := 2 3 ρ0 Cv (T ) Γ0 .

13
Plugging into the consistency condition, we get
b − λ̇ fm ] − fm G(T ) Ṫ + 2 µ fs [η : sb − λ̇ fs ] + fd ¨ + fp ˙ + ft Ṫ = 0
3 K fm [d : 1

or,
λ̇ (3 K fm2 + 2 µ fs2 ) = 3 K fm d : 1
b + 2 µ fs η : sb + fd ¨ + fp ˙ + [ft − fm G(T )] Ṫ = 0
or,
b + 2 µ fs η : sb + fd ¨ + fp ˙ + [ft − fm G(T )] Ṫ
3 K fm d : 1
λ̇ = . (72)
3 K fm2 + 2 µ fs2

1.4 Mass Balance:

The mass balance equation is

ρ̇ + ρ ∇ · v = 0 (73)
where ρ is the current mass density, and ρ̇ is the material time derivative of ρ.

1.5 Momentum Balance:

The balance of linear momentum can be written as

∇ · σ + ρ b = ρ v̇ (74)

where σ is the Cauchy stress, b is the body force per unit volume, and v̇ is the acceleration.
The momentum equation can be written as

∇ · (p 1 + s) + ρ b = ρ v̇

or,
∇ · (p 1 ) + ∇ · s + ρ b = ρ v̇ . (75)

1.6 Energy Balance:

The balance of energy can be written as

σ : d − ∇ · q + ρ ḣ = ρ ė (76)

where q is the heat flux per unit area, h is the heat source per unit mass, and e is the specific internal energy. A dot
over a quantity represents the material time derivative of that quantity.
We convert the energy equation into an approximate equation for heat conduction that includes heat source term
due to plastic dissipation ([8]):
ρ Cv Ṫ − ∇ · (κ · ∇T ) − ζ σ : dp = ρ ḣ (77)
where Cv is the specific heat at constant volume, κ is a second order tensor of thermal conductivity, and ζ is the
Taylor-Quinney coefficient. Expressing the stress and the plastic part of rate of deformation into volumetric and
deviatoric parts, we get
 
1 p p
ρ Cv Ṫ − ∇ · (κ · ∇T ) − ζ (p 1 + s) : tr(d ) 1 + η = ρ ḣ
3
14
or,  
1 1
ρ Cv Ṫ − ∇ · (κ · ∇T ) − ζ p tr(dp )1 : 1 + tr(dp ) s : 1 + p 1 : η p + s : η p = ρ ḣ
3 3
or,
ρ Cv Ṫ − ∇ · (κ · ∇T ) − ζ (p tr(dp ) + s : η p ) = ρ ḣ (78)

2 Weak Forms
2.1 Rate of deformation:

The weak form of the kinematic relation is given by

Z  
∗ 1 T

W : d− ∇v + (∇v) dΩ = 0 (79)
Ω 2

where W ∗ is an arbitrary second-order tensor valued weighting function. We write the rate of deformation and the
weighting function in terms of their volumetric and deviatoric components to get
Z    
∗ 1 ∗ 1 1 T

dev(W ) + tr(W ) 1 : η + θ̇ 1 − ∇v + (∇v) dΩ = 0
Ω 3 3 2
or,
Z 
1 1
dev(W ∗ ) : η + θ̇ dev(W ∗ ) : 1 − dev(W ∗ ) : ∇v + (∇v)T +
 
Ω 3 2

1 ∗ 1 ∗ 1 ∗
 T

tr(W ) 1 : η + θ̇ tr(W )1 : 1 − tr(W )1 : ∇v + (∇v) dΩ = 0
3 9 6
or,
Z  
1  1 1
dev(W ∗ ) : η − dev(W ∗ ) : ∇v + (∇v)T + θ̇ tr(W ∗ ) − tr(W ∗ ) tr(∇v) + tr((∇v)T )
  
dΩ = 0
Ω 2 3 6
or, Z  
∗ 1 ∗ 1 ∗
dev(W ) : (η − d) + θ̇ tr(W ) − tr(W ) ∇ · v dΩ = 0
Ω 3 3
or, Z  
1 ∗ 1 ∗ 1 ∗
− θ̇ dev(W ) : 1 + θ̇ tr(W ) − tr(W ) ∇ · v dΩ = 0
Ω 3 3 3
or, Z  
1 1
θ̇ tr(W ∗ ) − tr(W ∗ ) ∇ · v dΩ = 0
Ω 3 3
Define w∗ := 1/3 tr(W ∗ ). Then the weak form becomes
Z  
w∗ θ̇ − ∇ · v dΩ = 0 (80)
Ω

15
2.2 Constitutive Relations:

The weak forms of the constitutive relations are given by

Z ◦ 
W ∗∗ : s − C(p, T ) : η e dΩ = 0 (81)
Ω

and
ρ0 C02 (1 − J)[1 − Γ0 (1 − J)]
Z  
∗∗
w p− − 2 Γ0 e dΩ = 0 (82)
Ω [1 − Sα (1 − J)]2
where W ∗∗ is a tensor valued weighting function and w∗∗ is a scalar weighting function.

2.3 Momentum Balance:

To derive the weak form of the momentum equation, we multiply the momentum equation with a vector-valued
weighting function (w) and integrate over the domain (Ω). The weighting function (w) satisfies velocity boundary
conditions on the parts of the boundary where velocities are prescribed. Then we get,
Z Z
w · [∇ · (p 1 ) + ∇ · s + ρ b] dΩ = ρ w · v̇ dΩ . (83)
Ω Ω
Using the identity v · (∇ · S) = ∇ · (S T
· v) − S : ∇v, where S is a second-order tensor valued field and v is a
vector valued field, we get (using the symmetry of the stress tensor)
Z Z
{∇ · [(p 1 ) · w] − p 1 : ∇w + ∇ · (s · w) − s : ∇w + ρ w · b} dΩ = ρ w · v̇ dΩ
Ω Ω
or, Z Z
{∇ · (p w) − p ∇ · w + ∇ · (s · w) − s : ∇w + ρ w · b} dΩ = ρ w · v̇ dΩ . (84)
Ω Ω
Applying the divergence theorem, we have
Z Z Z Z
{−p ∇ · w − s : ∇w + ρ w · b} dΩ + n · (p w) dΓ + n · (s · w) dΓ = ρ w · v̇ dΩ
Ω Γ Γ Ω
or, Z Z Z Z
{ρ w · v̇ + p ∇ · w + s : ∇w} dΩ = ρ w · b dΩ + p n · w dΓ + n · (s · w) dΓ
Ω Ω Γ Γ
or, Z Z Z Z
{ρ w · v̇ + p ∇ · w + s : ∇w} dΩ = ρ w · b dΩ + p n · w dΓ + (sT · n) · w dΓ
Ω Ω Γ Γ
or,
Z Z Z Z
{ρ w · v̇ + p ∇ · w + s : ∇w} dΩ = ρ w · b dΩ + p n · w dΓ + (σ T · n − p n) · w dΓ (85)
Ω Ω Γ Γ
or, Z Z Z
{ρ w · v̇ + p ∇ · w + s : ∇w} dΩ = ρ w · b dΩ + (σ T · n) · w dΓ . (86)
Ω Ω Γ
In the above, n is the outward normal to the surface Γ.
If the applied surface traction is t̄ = σ T · n, we get (since w is zero on the part of the boundary where velocities
are specified)
Z Z Z
{ρ w · v̇ + p ∇ · w + s : ∇w} dΩ = ρ w · b dΩ + t̄ · w dΓ . (87)
Ω Ω Γt

16
2.4 Energy Balance:

The weak form of the heat conduction equation is given by

Z h i Z
p p
b ρ Cv Ṫ − ∇ · (κ · ∇T ) − ζ (p tr(d ) + s : η ) dΩ =
w w
b ρ ḣ dΩ .
Ω Ω

Using the identity φ(∇ · v) = ∇ · (φ v) − (∇φ) · v we get

Z Z Z Z
p p
b ρ Cv Ṫ dΩ −
w {∇ · [w(κ
b · ∇T )] − (∇w)b · κ · ∇T } dΩ − w
b ζ [p tr(d ) + s : η ] dΩ = w
b ρ ḣ dΩ .
Ω Ω Ω Ω

Using the divergence theorem, we have

Z n o Z n o Z
p p
w b · κ · ∇T dΩ =
b ρ Cv Ṫ + (∇w) b ρ ḣ + ζ [p tr(d ) + s : η ] dΩ + n · [w(κ
w b · ∇T )] dΓ .
Ω Ω Γ

If a heat flux (q = n · (κ · ∇T )) is specified on part of the boundary (Γq ), the weak form becomes
Z n o Z n o Z
p p
w b · κ · ∇T dΩ =
b ρ Cv Ṫ + (∇w) w
b ρ ḣ + ζ [p tr(d ) + s : η ] dΩ + b q dΓ .
w (88)
Ω Ω Γq

3 MPM Discretization
3.1 Rate of Deformation:

The weak form of the rate of deformation equation is

Z  
w∗ θ̇ − ∇ · v dΩ = 0 . (89)
Ω

We assume a discontinuous approximation for θ̇ over each cell and a weighting function w∗ of the form
nc
p nc
p
X X
∗
θ̇(x) ≈ θ̇q χq (x) ; w (x) = wp∗ χp (x) (90)
q=1 p=1

where ncp is the number of particles in a grid cell, wp∗ are the weighting functions at the particle locations, and
χp (x) are the basis functions. The velocities are approximated in the standard MPM way using
ng
X
v(x) ≈ vg Ng (x) (91)
g=1

where ng is the number of grid points, vg are the velocities at the grid vertices, and Ng (x) are the grid basis
functions.

17
Remark: The approximations for w∗ and θ̇ need not be a sum over the particles. We could alternatively assume a
constant value over each cell or a sum over the grid nodes with the same basis functions as the grid basis functions.
We have chosen a sum over particles in a cell to keep our formulation consistent with the standard MPM procedure.
Plugging the approximations into the weak form, we get the following equations that are valid over each cell
 c  c  c 
Z np Xnp ng 
X X
 wp∗ χp (x) θ̇q χq (x) − ∇ ·  vg Ng (x) dΩ = 0
Ωc p=1

q=1

g=1

or,  c 
nc nc
 
p Z np
X g 
X X
∗
wp χp (x) θ̇q χq (x) − vg · ∇Ng dΩ = 0 .
p=1 Ωc 
q=1 g=1


The abitrariness of wp∗ gives us a system of ncp equations

 c 
Z X np ncg 
X
χp (x) θ̇q χq (x) − vg · ∇Ng dΩ = 0 ; p = 1 . . . ncp
Ωc 
q=1

g=1

or,
nc
p Z  nc
g Z 
X X
χp (x)χq (x) dΩ θ̇q − χp (x)∇Ng dΩ · vg = 0 ; p = 1 . . . ncp
q=1 Ωc g=1 Ωc

or,
nc
p Z  nc
g Z 
X X
χp (x)χq (x) dΩ θ̇q = χp (x)∇Ng dΩ · vg ; p = 1 . . . ncp .
q=1 Ωc g=1 Ωc

If we identify the basis functions χp with the particle characteristic functions and use the property ([9], p. 485)
(
1 if x ∈ Ωp
χp (x) = (92)
0 otherwise

we get
"Z # ncg "Z #
X
χp (x)χp (x) dΩ θ̇p = χp (x)∇Ng dΩ · vg ; p = 1 . . . ncp .
Ωp ∩Ωc g=1 Ωp ∩Ωc

Define the volume of the particle inside the cell as

Z
Vpc := χp (x)χp (x) dΩ . (93)
Ωp ∩Ωc

We then have
ncg "Z #
X
Vpc θ̇p = χp (x)∇Ng dΩ · vg ; p = 1 . . . ncp .
g=1 Ωp ∩Ωc

Define the gradient weighting function as

Z
∇Npg := χp (x)∇Ng dΩ . (94)
Ωp ∩Ωc

18
Then,
nc
g
X
Vpc θ̇p = ∇Npg · vg ; p = 1 . . . ncp .
g=1
or,
ncg
1 X
θ̇p = ∇Npg · vg ; p = 1 . . . ncp (95)
Vpc
g=1

3.2 Constitutive Relations:

The weak form of the constitutive relation for pressure is

ρ0 C02 (1 − J(x))[1 − Γ0 (1 − J(x))]

Z  
∗∗
w (x) p(x) − − 2 Γ0 e(x) dΩ = 0 . (96)
Ω [1 − Sα (1 − J(x))]2

We use a discontinuous approximation for the pressure field (p) and a weighting function w∗∗ of the form
nc
p nc
p
X X
∗∗
p(x) ≈ pq χq (x) ; w (x) = wp∗∗ χp (x) . (97)
q=1 p=1

After plugging these approximations into the weak form, we get the following relations over each cell
 c  c 
np np
Z X X 2
ρ0 C0 (1 − J(x))[1 − Γ0 (1 − J(x))] 
∗∗
 wp χp (x)  pq χq (x) − − 2 Γ0 e(x) dΩ = 0
Ωc  [1 − Sα (1 − J(x))]2 
p=1 q=1

or,
 c 
nc
 
p np
ρ0 C02 (1 − J(x))[1 − Γ0 (1 − J(x))]
X Z X 
∗∗ 
wp χp (x) pq χq (x) − − 2 Γ0 e(x) dΩ = 0 .
Ωc  [1 − Sα (1 − J(x))]2 
p=1 q=1

Invoking the arbitrariness of wp∗∗ we get

 c 
np
Z X 2
ρ0 C0 (1 − J(x))[1 − Γ0 (1 − J(x))] 
χp (x) pq χq (x) − − 2 Γ0 e(x) dΩ = 0 ; p = 1 . . . ncp
Ωc  [1 − Sα (1 − J(x))]2 
q=1

or,
nc
p Z
ρ0 C02 (1 − J(x))[1 − Γ0 (1 − J(x))]
X  Z  
χp (x) χq (x) dΩ pq = χp (x) 2
− 2 Γ0 e(x) dΩ ; p = 1 . . . ncp
Ωc Ωc [1 − S α (1 − J(x))]
q=1

Once again, if we identify the basis functions χp with the particle characteristic functions and use the property
(
1 if x ∈ Ωp
χp (x) = (98)
0 otherwise

19
and we get
"Z #
ρ0 C02 (1 − J(x))[1 − Γ0 (1 − J(x))]
Z  
χp (x) χp (x) dΩ pp = χp (x) − 2 Γ0 e(x) dΩ ; p = 1 . . . ncp
Ωp ∩Ωc Ωp ∩Ωc [1 − Sα (1 − J(x))]2

or,
ρ0 C02 (1 − J(x))[1 − Γ0 (1 − J(x))]
Z  
Vpc pp = χp (x) 2
− 2 Γ0 e(x) dΩ ; p = 1 . . . ncp .
Ωp ∩Ωc [1 − S α (1 − J(x))]
Identifying the volume averaged J(x) and e(x) as θp and ep , respectively, we get

ρ0 C02 (1 − θp )[1 − Γ0 (1 − θp )]
Vpc pp = Vpc − 2 Γ0 ep ; p = 1 . . . ncp
[1 − Sα (1 − θp )]2
or,
ρ0 C02 (1 − θp )[1 − Γ0 (1 − θp )]
pp = − 2 Γ0 ep ; p = 1 . . . ncp . (99)
[1 − Sα (1 − θp )]2
The weak form of the constitutive relation for the deviatoric stress rate is given by
Z h◦ i
W ∗∗ (x) : s(x) − C(p, T, x) : η e (x) dΩ = 0 . (100)
Ω

We assume that the weighting function is

np
X
∗∗
W (x) = Wq∗∗ χq (x) (101)
q=1

and plug these into the weak form to get the following relations
 
Z X np h◦ i
 Wq∗∗ χq (x) : s(x) − C(p, T, x) : η e (x) dΩ = 0
Ω q=1

or,
np Z h◦ i 
X
Wq∗∗ : e
χq (x) s(x) − C(p, T, x) : η (x) dΩ = 0 .
q=1 Ω

The arbitrariness of Wq∗∗ gives

Z h◦ i
χq (x) s(x) − C(p, T, x) : η e (x) dΩ = 0 ; q = 1 . . . np
Ωq ∩Ω

Using arguments similar to those used for the pressure constitutive relation, we get

◦
sq = Cq (p, T ) : ηqe ; q = 1 . . . np (102)

20
3.3 Momentum Balance:

The weak form of the momentum balance relation is

Z Z Z
[ρ(x) w(x) · v̇(x) + p(x) ∇ · w + s(x) : ∇w] dΩ = ρ(x) w(x) · b(x) dΩ + t̄(x) · w(x) dΓ . (103)
Ω Ω Γt

The first step in the MPM discretization is to convert the integral over Ω into a sum of integrals over particles. To
achieve this we assume that the particle characteristic functions have the form
(
1 if x ∈ Ωp
χp (x) = (104)
0 otherwise

Then the weak form of the momentum balance equation can be written as
np Z
X
χp (x) [ρ(x) w(x) · v̇(x) + p(x) ∇ · w + s(x) : ∇w] dΩ =
p=1 Ωp ∩Ω
np Z Z
X
χp (x) ρ(x) w(x) · b(x) dΩ + t̄(x) · w(x) dΓ .
p=1 Ωp ∩Ω Γt

The weighting function is

ng
X
w(x) = wg Ng (x) . (105)
g=1

The velocity is approximated as

ng
X
v(x) ≈ vh Nh (x) . (106)
h=1
The material time derivative of v is also approximated in a similar manner (see [10]). Thus,
ng
X
v̇(x) ≈ v̇h Nh (x) . (107)
h=1

Plugging these into the left hand side of the momentum equation we get
   (  
np Z X ng  ng )
X ng 
X X
LHS = χp (x) ρ(x) wg Ng (x) · v̇h Nh (x) + p(x) ∇ · wg Ng (x)
p=1 Ωp ∩Ω
   
g=1 h=1 g=1
 
X ng 
+s(x) : ∇ wg Ng (x)  dΩ
 
g=1

or,
 )
ng np ( ng !
X X X Z
LHS = wg ·  χp (x) ρ(x) Ng (x) Nh (x) dΩ v̇h  +
g=1 p=1 h=1 Ωp ∩Ω
   
ng np Z ng np Z
X X X X
wg ·  χp (x) p(x) ∇Ng dΩ + wg ·  χp (x) s(x) · ∇Ng dΩ
g=1 p=1 Ωp ∩Ω g=1 p=1 Ωp ∩Ω

21
or,

ng np "( ng ! )
X X X Z
LHS = wg ·  χp (x) ρ(x) Ng (x) Nh (x) dΩ v̇h +
g=1 p=1 h=1 Ωp ∩Ω
Z Z ##
χp (x) p(x) ∇Ng dΩ + χp (x) s(x) · ∇Ng dΩ
Ωp ∩Ω Ωp ∩Ω

Similarly, the right hand side of the weak form of momentum balance can be written as
   
np Z ng Z ng
X X X
RHS = χp (x) ρ(x)  wg Ng (x) · b(x) dΩ + t̄(x) ·  wg Ng (x) dΓ
p=1 Ωp ∩Ω g=1 Γt g=1

or,  
ng np Z ng Z 
X X X
RHS = wg ·  χp (x) ρ(x) Ng (x) b(x) dΩ + wg · t̄(x) Ng (x) dΓ
g=1 p=1 Ωp ∩Ω g=1 Γt

or,   
ng Xnp Z  Z
X
RHS = wg ·  χp (x) ρ(x) Ng (x) b(x) dΩ + t̄(x) Ng (x) dΓ
g=1
 Ωp ∩Ω
p=1
 Γt

Combining the left and right hand sides and invoking the arbitrariness of wg , we get
np "( ng ! )
X X Z
χp (x) ρ(x) Ng (x) Nh (x) dΩ v̇h +
p=1 h=1 Ωp ∩Ω
Z Z #
χp (x) p(x) ∇Ng dΩ + χp (x) s(x) · ∇Ng dΩ =
Ωp ∩Ω Ωp ∩Ω
 
Xnp Z  Z
χp (x) ρ(x) Ng (x) b(x) dΩ + t̄(x) Ng (x) dΓ ; g = 1 . . . ng

p=1 Ωp ∩Ω  Γt

If ρ(x), p(x), s(x), and b(x) are assumed to be constant over each particle, we can replace these quantities with
the values at the particles to get
np "( ng Z ! )
X X
ρp χp (x) Ng (x) Nh (x) dΩ v̇h +
p=1 h=1 Ωp ∩Ω
Z Z #
pp χp (x) ∇Ng dΩ + sp · χp (x)∇Ng dΩ =
Ωp ∩Ω Ωp ∩Ω
 
Xnp Z  Z
ρp bp χp (x) Ng (x) dΩ + t̄(x) Ng (x) dΓ ; g = 1 . . . ng .

p=1 Ωp ∩Ω  Γt

Then, using the definition of the gradient weighting function (94), and defining
Z
Npg := χp (x) Ng (x) dΩ (108)
Ωp ∩Ω
22
we get
 !    
ng  np Z  Xnp  Xnp 
X X
ρp χp (x) Ng (x) Nh (x) dΩ v̇h + pp ∇Npg + sp · ∇Npg =
h=1

p=1 Ωp ∩Ω  
p=1
 
p=1

 
X np  Z
ρp bp Npg + t̄(x) Ng (x) dΓ ; g = 1 . . . ng

p=1
 Γt

Define the mass matrix as

np Z
X
mgh := ρp χp (x) Ng (x) Nh (x) dΩ . (109)
p=1 Ωp ∩Ω

Define the internal force vector as

np
X
fgint := [pp ∇Npg + sp · ∇Npg ] . (110)
p=1

Define the body force vector as

np
X
fgbod := ρp bp Npg . (111)
p=1

Also define the external force vector as Z

fgext := t̄(x) Ng (x) dΓ (112)
Γt
Then we get the semidiscrete system of equations
ng
X
mgh v̇h + fgint = fgbod + fgext ; g = 1 . . . ng (113)
h=1

3.4 Energy Balance:

The weak form of the energy balance equation is

Z n o
w(x)
b ρ(x) Cv (x) Ṫ (x) + ∇wb · κ(x) · ∇T dΩ = (114)
Ω
Z n o Z
p p
w(x)
b ρ(x) ḣ(x) + ζ [p(x) tr(d (x)) + s(x) : η (x)] dΩ + w(x)
b q(x) dΓ .
Ω Γq

As before, we express this equation as a sum over particles to get

np Z n o
X
χp (x) w(x)
b ρ(x) Cv (x) Ṫ (x) + ∇w
b · κ(x) · ∇T dΩ =
p=1 Ωp ∩Ω
np Z n o Z
X
p p
χp (x) w(x)
b ρ(x) ḣ(x) + ζ [p(x) tr(d (x)) + s(x) : η (x)] dΩ + w(x)
b q(x) dΓ
p=1 Ωp ∩Ω Γq

We approximate T using
ng
X
T (x) ≈ Tg Ng (x) . (115)
g=1
23
The material time derivative of T is then given by
ng
X
Ṫ (x) ≈ Ṫg Ng (x) . (116)
g=1

The weighting functions are also chosen to be of the same form:

ng
X
w(x)
b = w
bg Ng (x) . (117)
g=1

Plugging these into the weak form, we get

 
np Z ng ng !
X X X
χp (x)  bg Ng (x) ρ(x) Cv (x)
w Ṫh Nh (x) dΩ
p=1 Ωp ∩Ω g=1 h=1
 
np Z ng ng !
X X X
+ χp (x) ∇  bg Ng (x) · κ(x) · ∇
w Th Nh (x) dΩ
p=1 Ωp ∩Ω g=1 h=1
 
np Z ng n o
X X
= χp (x)  w
bg Ng (x) ρ(x) ḣ(x) + ζ [p(x) tr(dp (x)) + s(x) : η p (x)] dΩ
p=1 Ωp ∩Ω g=1
 
Z ng
X
+  w
bg Ng (x) q(x) dΓ
Γq g=1

For simplicity, let us assume that κ(x) is isotropic, i.e., κ(x) = κ(x). Then,
 
ng X np Z ng !

X X
wbg χp (x) Ng (x) ρ(x) Cv (x) Ṫh Nh (x) dΩ
p=1 Ωp ∩Ω
 
g=1 h=1
 
ng X np Z ng !

X X
+ wbg χp (x) κ(x) ∇Ng · Ṫh ∇Nh dΩ
p=1 Ωp ∩Ω
 
g=1 h=1
 
ng np Z n o
X X
= wbg  χp (x) Ng (x) ρ(x) ḣ(x) + ζ [p(x) tr(dp (x)) + s(x) : η p (x)] dΩ
g=1 p=1 Ωp ∩Ω
ng (Z )
X
+ w
bg Ng (x) q(x) dΓ
g=1 Γq

Invoking the arbitrariness of w

bg , we get
np Z ng ! np Z ng !
X X X X
χp (x) Ng (x) ρ(x) Cv (x) Ṫh Nh (x) dΩ + χp (x) κ(x) ∇Ng · Th ∇Nh dΩ
p=1 Ωp ∩Ω h=1 p=1 Ωp ∩Ω h=1
np Z n o
X
= χp (x) Ng (x) ρ(x) ḣ(x) + ζ [p(x) tr(dp (x)) + s(x) : η p (x)] dΩ
p=1 Ωp ∩Ω
Z
+ Ng (x) q(x) dΓ ; g = 1 . . . ng
Γq

24
or,
   
ng np Z ng np Z
X X X X
 χp (x) Ng (x) Nh (x) ρ(x) Cv (x) dΩ Ṫh +  χp (x) κ(x) ∇Ng · ∇Nh dΩ Th
h=1 p=1 Ωp ∩Ω h=1 p=1 Ωp ∩Ω
np
XZ n o
p p
= χp (x) Ng (x) ρ(x) ḣ(x) + ζ [p(x) tr(d (x)) + s(x) : η (x)] dΩ
p=1 Ωp ∩Ω
Z
+ Ng (x) q(x) dΓ ; g = 1 . . . ng
Γq

If we assume that ρ(x), Cv (x), κ(x), ḣ(x), p(x), s(x), and dp (x) are constant over each particle, we get
   
ng np Z ng np Z
X X X X
 ρp Cvp χp (x) Ng (x) Nh (x) dΩ Ṫh +  κp χp (x) ∇Ng · ∇Nh dΩ Th
h=1 p=1 Ωp ∩Ω h=1 p=1 Ωp ∩Ω
np Z np Z
X X
= ρp ḣp χp (x) Ng (x) dΩ + ζ pp tr(dpp ) χp (x) Ng (x) dΩ
p=1 Ωp ∩Ω p=1 Ωp ∩Ω
np Z Z
X
+ ζ sp : ηpp χp (x) Ng (x) dΩ + Ng (x) q(x) dΓ ; g = 1 . . . ng
p=1 Ωp ∩Ω Γq

Using the definition of the particle weighting functions (108), we have

   
ng np Z ng np Z
X X X X
 ρp Cvp χp (x) Ng (x) Nh (x) dΩ Ṫh +
  κp χp (x) ∇Ng · ∇Nh dΩ Th
h=1 p=1 Ωp ∩Ω h=1 p=1 Ωp ∩Ω
np np np Z
X X X
= ρp ḣp Npg + ζ pp tr(dpp ) Npg + ζ sp : ηpp Npg + Ng (x) q(x) dΓ ; g = 1 . . . ng
p=1 p=1 p=1 Γq

or,
   
ng np Z ng np Z
X X X X
 ρp Cvp χp (x) Ng (x) Nh (x) dΩ Ṫh +
  κp χp (x) ∇Ng · ∇Nh dΩ Th
h=1 p=1 Ωp ∩Ω h=1 p=1 Ωp ∩Ω
np Z
Xh i
= ρp ḣp + ζ pp tr(dpp ) + ζ sp : ηpp Npg + Ng (x) q(x) dΓ ; g = 1 . . . ng .
p=1 Γq

Define the thermal inertia matrix as

np Z
X
Mgh := ρp Cvp χp (x) Ng (x) Nh (x) dΩ , (118)
p=1 Ωp ∩Ω

the thermal stiffness matrix as

np Z
X
Kgh := κp χp (x) ∇Ng · ∇Nh dΩ , (119)
p=1 Ωp ∩Ω

25
the thermal source vector as
np h i
X
Sg := ρp ḣp + ζ pp tr(dpp ) + ζ sp : ηpp Npg , (120)
p=1

and the external heat flux vector as Z

qgext := Ng (x) q(x) dΓ . (121)
Γq

Then we get the following semidiscrete system of equations:

ng ng
X X
Mgh Ṫh + Kgh Th = Sg + qgext ; g = 1 . . . ng . (122)
h=1 h=1

4 Appendix
1. The integral
Z b(t)
F (t) = f (x, t) dx
a(t)

is a function of the parameter t. Show that the derivative of F is given by

!
dF d
Z b(t) Z b(t) ∂f (x, t) ∂b(t) ∂a(t)
= f (x, t) dx = dx + f [b(t), t] − f [a(t), t] .
dt dt a(t) a(t) ∂t ∂t ∂t

This relation is also known as the Leibnitz rule.

The following proof is taken from [11].
We have,
dF F (t + ∆t) − F (t)
= lim .
dt ∆t→0 ∆t
Now,
"Z #
F (t + ∆t) − F (t) 1 b(t+∆t) Z b(t)
= f (x, t + ∆t) dx − f (x, t) dx
∆t ∆t a(t+∆t) a(t)

1
»Z b+∆b Z b –
≡ f (x, t + ∆t) dx − f (x, t) dx
∆t a+∆a a
1
» Z a+∆a Z b+∆b Z b –
= − f (x, t + ∆t) dx + f (x, t + ∆t) dx − f (x, t) dx
∆t a a a
1
» Z a+∆a Z b Z b+∆b Z b –
= − f (x, t + ∆t) dx + f (x, t + ∆t) dx + f (x, t + ∆t) dx − f (x, t) dx
∆t a a b a
Z b b+∆b a+∆a
f (x, t + ∆t) − f (x, t) 1 1
Z Z
= dx + f (x, t + ∆t) dx − f (x, t + ∆t) dx .
a ∆t ∆t b ∆t a

Since f (x, t) is essentially constant over the infinitesimal intervals a < x < a + ∆a and b < x < b + ∆b, we may write

F (t + ∆t) − F (t)
Z b f (x, t + ∆t) − f (x, t) ∆b ∆a
≈ dx + f (b, t + ∆t) − f (a, t + ∆t) .
∆t a ∆t ∆t ∆t

Taking the limit as ∆t → 0, we get

" # "Z # " # " #
F (t + ∆t) − F (t) b f (x, t + ∆t) − f (x, t) ∆b ∆a
lim = lim dx + lim f (b, t + ∆t) − lim f (a, t + ∆t)
∆t→0 ∆t ∆t→0 a ∆t ∆t→0 ∆t ∆t→0 ∆t

or,
dF (t)
Z b(t) ∂f (x, t) ∂b(t) ∂a(t)
= dx + f [b(t), t] − f [a(t), t] . 
dt a(t) ∂t ∂t ∂t

26
2. Let Ω(t) be a region in Euclidean space with boundary ∂Ω(t). Let x(t) be the positions of points in the region and let v(x, t) be the velocity field in
the region. Let n(x, t) be the outward unit normal to the boundary. Let f (x, t) be a vector field in the region (it may also be a scalar field). Show that

Z ! Z Z
d ∂f
f dV = dV + (v · n)f dA .
dt Ω(t) Ω(t) ∂t ∂Ω(t)

This relation is also known as the Reynold’s Transport Theorem and is a generalization of the Leibnitz rule.
This proof is taken from [12] (also see [13]).
Let Ω0 be reference configuration of the region Ω(t). Let the motion and the deformation gradient be given by

x = ϕ(X, t) ; F (X, t) = ∇0 ϕ .

Let J(X, t) = det[F (x, t)]. Then, integrals in the current and the reference configurations are related by
Z Z Z
f (x, t) dV = f [ϕ(X, t), t] J(X, t) dV = f̂ (X, t) J(X, t) dV .
Ω(t) Ω0 Ω0

The time derivative of an integral over a volume is defined as

! !
1
Z Z Z
d
f (x, t) dV = lim f (x, t + ∆t) dV − f (x, t) dV .
dt Ω(t) ∆t→0 ∆t Ω(t+∆t) Ω(t)

Converting into integrals over the reference configuration, we get

!
1
Z „Z Z «
d
f (x, t) dV = lim f̂ (X, t + ∆t) J(X, t + ∆t) dV − f̂ (X, t) J(X, t) dV .
dt Ω(t) ∆t→0 ∆t Ω0 Ω0

Since Ω0 is independent of time, we have

! " #
f̂ (X, t + ∆t) J(X, t + ∆t) − f̂ (X, t) J(X, t)
Z Z
d
f (x, t) dV = lim dV
dt Ω(t) Ω0 ∆t→0 ∆t
Z
∂
= [f̂ (X, t) J(X, t)] dV
Ω0 ∂t
Z „ «
∂ ∂
= [f̂ (X, t)] J(X, t) + f̂ (X, t) [J(X, t)] dV
Ω0 ∂t ∂t

Now, the time derivative of det F is given by (see [13], p. 77)

∂J(X, t) ∂
= (det F ) = (det F )(∇ · v) = J(X, t) ∇ · v(ϕ(X, t), t) = J(X, t) ∇ · v(x, t) .
∂t ∂t

Therefore, !
Z Z „ «
d ∂
f (x, t) dV = [f̂ (X, t)] J(X, t) + f̂ (X, t) J(X, t) ∇ · v(x, t) dV
dt Ω(t) Ω0 ∂t
Z „ «
∂
= [f̂ (X, t)] + f̂ (X, t) ∇ · v(x, t) J(X, t) dV
Ω0 ∂t
Z “ ”
= ḟ (x, t) + f (x, t) ∇ · v(x, t) dV
Ω(t)

where ḟ is the material time derivative of f . Now, the material derivative is given by

∂f (x, t)
ḟ (x, t) = + [∇f (x, t)] · v(x, t) .
∂t

Therefore, !
Z Z „ «
d ∂f (x, t)
f (x, t) dV = + [∇f (x, t)] · v(x, t) + f (x, t) ∇ · v(x, t) dV
dt Ω(t) Ω(t) ∂t
or, !
Z Z „ «
d ∂f
f dV = + ∇f · v + f ∇ · v dV .
dt Ω(t) Ω(t) ∂t
27
 Using the identity
∇ · (v ⊗ w) = v(∇ · w) + ∇v · w

we then have ! „ «
Z Z
d ∂f
f dV = + ∇ · (f ⊗ v) dV .
dt Ω(t) Ω(t) ∂t

Using the divergence theorem and the identity (a ⊗ b) · n = (b · n)a we have

Z ! Z Z Z Z
d ∂f ∂f
f dV = dV + (f ⊗ v) · n dV = dV + (v · n)f dV .
dt Ω(t) Ω(t) ∂t ∂Ω(t) Ω(t) ∂t ∂Ω(t)

3. Show that the balance of mass can be expressed as:

ρ̇ + ρ ∇ · v = 0

where ρ(x, t) is the current mass density, ρ̇ is the material time derivative of ρ, and v(x, t) is the velocity of physical particles in the body Ω
bounded by the surface ∂Ω.
Recall that the general equation for the balance of a physical quantity f (x, t) is given by
»Z – Z Z Z
d
f (x, t) dV = f (x, t)[un (x, t) − v(x, t) · n(x, t)] dA + g(x, t) dA + h(x, t) dV .
dt Ω ∂Ω ∂Ω Ω

To derive the equation for the balance of mass, we assume that the physical quantity of interest is the mass density ρ(x, t). Since mass is neither
created or destroyed, the surface and interior sources are zero, i.e., g(x, t) = h(x, t) = 0. Therefore, we have
»Z – Z
d
ρ(x, t) dV = ρ(x, t)[un (x, t) − v(x, t) · n(x, t)] dA .
dt Ω ∂Ω

Let us assume that the volume Ω is a control volume (i.e., it does not change with time). Then the surface ∂Ω has a zero velocity (un = 0) and we get
Z Z
∂ρ
dV = − ρ (v · n) dA .
Ω ∂t ∂Ω

Using the divergence theorem Z Z

∇ · v dV = v · n dA
Ω ∂Ω

we get Z Z
∂ρ
dV = − ∇ · (ρ v) dV.
Ω ∂t Ω
or, Z » –
∂ρ
+ ∇ · (ρ v) dV = 0 .
Ω ∂t
Since Ω is arbitrary, we must have
∂ρ
+ ∇ · (ρ v) = 0 .
∂t
Using the identity
∇ · (ϕ v) = ϕ ∇ · v + ∇ϕ · v

we have
∂ρ
+ ρ ∇ · v + ∇ρ · v = 0 .
∂t
Now, the material time derivative of ρ is defined as
∂ρ
ρ̇ = + ∇ρ · v .
∂t
Therefore,
ρ̇ + ρ ∇ · v = 0 .

28
4. Show that the balance of linear momentum can be expressed as:

ρ v̇ − ∇ · σ − ρ b = 0

where ρ(x, t) is the mass density, v(x, t) is the velocity, σ(x, t) is the Cauchy stress, and ρ b is the body force density.
Recall the general equation for the balance of a physical quantity
»Z – Z Z Z
d
f (x, t) dV = f (x, t)[un (x, t) − v(x, t) · n(x, t)] dA + g(x, t) dA + h(x, t) dV .
dt Ω ∂Ω ∂Ω Ω

In this case the physical quantity of interest is the momentum density, i.e., f (x, t) = ρ(x, t) v(x, t). The source of momentum flux at the surface is
the surface traction, i.e., g(x, t) = t. The source of momentum inside the body is the body force, i.e., h(x, t) = ρ(x, t) b(x, t). Therefore, we have
»Z – Z Z Z
d
ρ v dV = ρ v[un − v · n] dA + t dA + ρ b dV .
dt Ω ∂Ω ∂Ω Ω

The surface tractions are related to the Cauchy stress by

t=σ·n.
Therefore, »Z – Z Z Z
d
ρ v dV = ρ v[un − v · n] dA + σ · n dA + ρ b dV .
dt Ω ∂Ω ∂Ω Ω

Let us assume that Ω is an arbitrary fixed control volume. Then,

Z Z Z Z
∂
(ρ v) dV = − ρ v (v · n) dA + σ · n dA + ρ b dV .
Ω ∂t ∂Ω ∂Ω Ω

Now, from the definition of the tensor product we have (for all vectors a)

(u ⊗ v) · a = (a · v) u .

Therefore, Z Z Z Z
∂
(ρ v) dV = − ρ (v ⊗ v) · n dA + σ · n dA + ρ b dV .
Ω ∂t ∂Ω ∂Ω Ω

Using the divergence theorem Z Z

∇ · v dV = v · n dA
Ω ∂Ω

we have Z Z Z Z
∂
(ρ v) dV = − ∇ · [ρ (v ⊗ v)] dV + ∇ · σ dV + ρ b dV
Ω ∂t Ω Ω Ω
or, Z » –
∂
(ρ v) + ∇ · [(ρ v) ⊗ v] − ∇ · σ − ρ b dV = 0 .
Ω ∂t
Since Ω is arbitrary, we have
∂
(ρ v) + ∇ · [(ρ v) ⊗ v] − ∇ · σ − ρ b = 0 .
∂t
Using the identity
∇ · (u ⊗ v) = (∇ · v)u + (∇u) · v
we get
∂ρ ∂v
v+ρ + (∇ · v)(ρv) + ∇(ρ v) · v − ∇ · σ − ρ b = 0
∂t ∂t
or, » –
∂ρ ∂v
+ρ∇·v v+ρ + ∇(ρ v) · v − ∇ · σ − ρ b = 0
∂t ∂t
Using the identity
∇(ϕ v) = ϕ ∇v + v ⊗ (∇ϕ)
we get » –
∂ρ ∂v
+ρ∇·v v+ρ + [ρ ∇v + v ⊗ (∇ρ)] · v − ∇ · σ − ρ b = 0
∂t ∂t
From the definition
(u ⊗ v) · a = (a · v) u
29
 we have
[v ⊗ (∇ρ)] · v = [v · (∇ρ)] v .
Hence, » –
∂ρ ∂v
+ρ∇·v v+ρ + ρ ∇v · v + [v · (∇ρ)] v − ∇ · σ − ρ b = 0
∂t ∂t
or, » –
∂ρ ∂v
+ ∇ρ · v + ρ ∇ · v v + ρ + ρ ∇v · v − ∇ · σ − ρ b = 0 .
∂t ∂t
The material time derivative of ρ is defined as
∂ρ
ρ̇ = + ∇ρ · v .
∂t
Therefore,
∂v
[ρ̇ + ρ ∇ · v] v + ρ + ρ ∇v · v − ∇ · σ − ρ b = 0 .
∂t
From the balance of mass, we have
ρ̇ + ρ ∇ · v = 0 .
Therefore,
∂v
ρ + ρ ∇v · v − ∇ · σ − ρ b = 0 .
∂t
The material time derivative of v is defined as
∂v
v̇ = + ∇v · v .
∂t
Hence,
ρ v̇ − ∇ · σ − ρ b = 0 .

5. Show that the balance of angular momentum can be expressed as:

σ = σT

We assume that there are no surface couples on ∂Ω or body couples in Ω. Recall the general balance equation
»Z – Z Z Z
d
f (x, t) dV = f (x, t)[un (x, t) − v(x, t) · n(x, t)] dA + g(x, t) dA + h(x, t) dV .
dt Ω ∂Ω ∂Ω Ω

In this case, the physical quantity to be conserved the angular momentum density, i.e., f = x × (ρ v). The angular momentum source at the surface
is then g = x × t and the angular momentum source inside the body is h = x × (ρ b). The angular momentum and moments are calculated with
respect to a fixed origin. Hence we have
»Z – Z Z Z
d
x × (ρ v) dV = [x × (ρ v)][un − v · n] dA + x × t dA + x × (ρ b) dV .
dt Ω ∂Ω ∂Ω Ω

Assuming that Ω is a control volume, we have

Z » – Z Z Z
∂
x× (ρ v) dV = − [x × (ρ v)][v · n] dA + x × t dA + x × (ρ b) dV .
Ω ∂t ∂Ω ∂Ω Ω

Using the definition of a tensor product we can write

[x × (ρ v)][v · n] = [[x × (ρ v)] ⊗ v] · n .

Also, t = σ · n. Therefore we have

Z » – Z Z Z
∂
x× (ρ v) dV = − [[x × (ρ v)] ⊗ v] · n dA + x × (σ · n) dA + x × (ρ b) dV .
Ω ∂t ∂Ω ∂Ω Ω

Using the divergence theorem, we get

Z » – Z Z Z
∂
x× (ρ v) dV = − ∇ · [[x × (ρ v)] ⊗ v] dV + x × (σ · n) dA + x × (ρ b) dV .
Ω ∂t Ω ∂Ω Ω

To convert the surface integral in the above equation into a volume integral, it is convenient to use index notation. Thus,
»Z – Z Z Z
x × (σ · n) dA = eijk xj σkl nl dA = Ail nl dA = A · n dA
∂Ω i ∂Ω ∂Ω ∂Ω
30
where [ ]i represents the i-th component of the vector. Using the divergence theorem
Z Z Z Z
∂Ail ∂
A · n dA = ∇ · A dV = dV = (eijk xj σkl ) dV .
∂Ω Ω Ω ∂xl Ω ∂xl

Differentiating,
Z Z » – Z » – Z
∂σkl ∂σkl ˆ ˜
A · n dA = eijk δjl σkl + eijk xj dV = eijk σkj + eijk xj dV = eijk σkj + eijk xj [∇ · σ]l dV .
∂Ω Ω ∂xl Ω ∂xl Ω

Expressed in direct tensor notation, Z Z h i

A · n dA = [E : σ T ]i + [x × (∇ · σ)]i dV
∂Ω Ω

where E is the third-order permutation tensor. Therefore,

»Z – Z h i
x × (σ · n) dA == [E : σ T ]i + [x × (∇ · σ)]i dV
∂Ω i Ω

or, Z Z h i
x × (σ · n) dA == E : σ T + x × (∇ · σ) dV .
∂Ω Ω

The balance of angular momentum can then be written as

Z » – Z Z h Z
∂ i
x× (ρ v) dV = − ∇ · [[x × (ρ v)] ⊗ v] dV + E : σ T + x × (∇ · σ) dV + x × (ρ b) dV .
Ω ∂t Ω Ω Ω

Since Ω is an arbitrary volume, we have

» –
∂
x× (ρ v) = −∇ · [[x × (ρ v)] ⊗ v] + E : σ T + x × (∇ · σ) + x × (ρ b)
∂t

or, » –
∂
x× (ρ v) − ∇ · σ − ρ b = −∇ · [[x × (ρ v)] ⊗ v] + E : σ T .
∂t
Using the identity,
∇ · (u ⊗ v) = (∇ · v)u + (∇u) · v
we get
∇ · [[x × (ρ v)] ⊗ v] = (∇ · v)[x × (ρ v)] + (∇[x × (ρ v)]) · v .
The second term on the right can be further simplified using index notation as follows.

∂
[(∇[x × (ρ v)]) · v]i = [(∇[ρ (x × v)]) · v]i = (ρ eijk xj vk ) vl
∂xl
» –
∂ρ ∂xj ∂vk
= eijk xj v k v l + ρ v k v l + ρ xj vl
∂xl ∂xl ∂xl
„ « „ «
∂ρ ∂vk
= (eijk xj vk ) vl + ρ (eijk δjl vk vl ) + eijk xj ρ vl
∂xl ∂xl
= [(x × v)(∇ρ · v) + ρ v × v + x × (ρ ∇v · v)]i
= [(x × v)(∇ρ · v) + x × (ρ ∇v · v)]i .

Therefore we can write

∇ · [[x × (ρ v)] ⊗ v] = (ρ ∇ · v)(x × v) + (∇ρ · v)(x × v) + x × (ρ ∇v · v) .
The balance of angular momentum then takes the form
» –
∂
x× (ρ v) − ∇ · σ − ρ b = −(ρ ∇ · v)(x × v) − (∇ρ · v)(x × v) − x × (ρ ∇v · v) + E : σ T
∂t

or, » –
∂
x× (ρ v) + ρ ∇v · v − ∇ · σ − ρ b = −(ρ ∇ · v)(x × v) − (∇ρ · v)(x × v) + E : σ T
∂t
or, » –
∂v ∂ρ
x× ρ + v + ρ ∇v · v − ∇ · σ − ρ b = −(ρ ∇ · v)(x × v) − (∇ρ · v)(x × v) + E : σ T
∂t ∂t
31
 The material time derivative of v is defined as
∂v
v̇ = + ∇v · v .
∂t
Therefore,
∂ρ
x × [ρ v̇ − ∇ · σ − ρ b] = −x × v + −(ρ ∇ · v)(x × v) − (∇ρ · v)(x × v) + E : σ T .
∂t
Also, from the conservation of linear momentum
ρ v̇ − ∇ · σ − ρ b = 0 .

Hence,
∂ρ
0=x× v + (ρ ∇ · v)(x × v) + (∇ρ · v)(x × v) − E : σ T
∂t
„ «
∂ρ
= + ρ∇ · v + ∇ρ · v (x × v) − E : σ T .
∂t
The material time derivative of ρ is defined as
∂ρ
ρ̇ = + ∇ρ · v .
∂t
Hence,
(ρ̇ + ρ ∇ · v)(x × v) − E : σ T = 0 .

From the balance of mass

ρ̇ + ρ ∇ · v = 0 .

Therefore,
E : σT = 0 .

In index notation,
eijk σkj = 0 .

Expanding out, we get

σ12 − σ21 = 0 ; σ23 − σ32 = 0 ; σ31 − σ13 = 0 .

Hence,
σ = σT

6. Show that the balance of energy can be expressed as:

ρ ė − σ : (∇v) + ∇ · q − ρ s = 0

where ρ(x, t) is the mass density, e(x, t) is the internal energy per unit mass, σ(x, t) is the Cauchy stress, v(x, t) is the particle velocity, q is the
heat flux vector, and s is the rate at which energy is generated by sources inside the volume (per unit mass).
Recall the general balance equation
»Z – Z Z Z
d
f (x, t) dV = f (x, t)[un (x, t) − v(x, t) · n(x, t)] dA + g(x, t) dA + h(x, t) dV .
dt Ω ∂Ω ∂Ω Ω

In this case, the physical quantity to be conserved the total energy density which is the sum of the internal energy density and the kinetic energy
density, i.e., f = ρ e + 1/2 ρ |v · v|. The energy source at the surface is a sum of the rate of work done by the applied tractions and the rate of heat
leaving the volume (per unit area), i.e, g = v · t − q · n where n is the outward unit normal to the surface. The energy source inside the body is the
sum of the rate of work done by the body forces and the rate of energy generated by internal sources, i.e., h = v · (ρb) + ρ s.
Hence we have
»Z „ « – Z „ « Z Z
d 1 1
ρ e+ v·v dV = ρ e+ v · v (un − v · n) dA + (v · t − q · n) dA + ρ (v · b + s) dV .
dt Ω 2 ∂Ω 2 ∂Ω Ω

Let Ω be a control volume that does not change with time. Then we get
Z » „ «– Z „ « Z Z
∂ 1 1
ρ e+ v·v dV = − ρ e + v · v (v · n) dA + (v · t − q · n) dA + ρ (v · b + s) dV .
Ω ∂t 2 ∂Ω 2 ∂Ω Ω

Using the relation t = σ · n, the identity v · (σ · n) = (σ T · v) · n, and invoking the symmetry of the stress tensor, we get
Z » „ «– Z „ « Z Z
∂ 1 1
ρ e+ v·v dV = − ρ e + v · v (v · n) dA + (σ · v − q) · n dA + ρ (v · b + s) dV .
Ω ∂t 2 ∂Ω 2 ∂Ω Ω

32
 We now apply the divergence theorem to the surface integrals to get
Z » „ «– Z » „ « – Z Z Z
∂ 1 1
ρ e+ v·v dV = − ∇ · ρ e + v · v v dA + ∇ · (σ · v) dA − ∇ · q dA + ρ (v · b + s) dV .
Ω ∂t 2 Ω 2 Ω Ω Ω

Since Ω is arbitrary, we have

» „ «– » „ « –
∂ 1 1
ρ e+ v·v = −∇ · ρ e + v · v v + ∇ · (σ · v) − ∇ · q + ρ (v · b + s) .
∂t 2 2

Expanding out the left hand side, we have

» „ «– „ « „ «
∂ 1 ∂ρ 1 ∂e 1 ∂
ρ e+ v·v = e+ v·v +ρ + (v · v)
∂t 2 ∂t 2 ∂t 2 ∂t
„ «
∂ρ 1 ∂e ∂v
= e+ v·v +ρ +ρ ·v.
∂t 2 ∂t ∂t

For the first term on the right hand side, we use the identity ∇ · (ϕ v) = ϕ ∇ · v + ∇ϕ · v to get
» „ « – „ « » „ «–
1 1 1
∇· ρ e+ v·v v =ρ e+ v·v ∇·v+∇ ρ e+ v·v ·v
2 2 2
„ « „ « „ «
1 1 1
=ρ e+ v·v ∇·v+ e+ v · v ∇ρ · v + ρ ∇ e + v · v · v
2 2 2
„ « „ «
1 1 1
=ρ e+ v·v ∇·v+ e+ v · v ∇ρ · v + ρ ∇e · v + ρ ∇(v · v) · v
2 2 2
„ « „ «
1 1
=ρ e+ v·v ∇·v+ e+ v · v ∇ρ · v + ρ ∇e · v + ρ (∇vT · v) · v
2 2
„ « „ «
1 1
=ρ e+ v·v ∇·v+ e+ v · v ∇ρ · v + ρ ∇e · v + ρ (∇v · v) · v .
2 2

For the second term on the right we use the identity ∇ · (S T · v) = S : ∇v + (∇ · S) · v and the symmetry of the Cauchy stress tensor to get

∇ · (σ · v) = σ : ∇v + (∇ · σ) · v .

After collecting terms and rearranging, we get

„ «„ « „ « „ «
∂ρ 1 ∂v ∂e
+ ρ ∇ · v + ∇ρ · v e+ v·v + ρ + ρ ∇v · v − ∇ · σ − ρ b · v + ρ + ∇e · v + −σ : ∇v + ∇ · q − ρ s = 0 .
∂t 2 ∂t ∂t

Applying the balance of mass to the first term and the balance of linear momentum to the second term, and using the material time derivative of the
internal energy
∂e
ė = + ∇e · v
∂t
we get the final form of the balance of energy:
ρ ė − σ : ∇v + ∇ · q − ρ s = 0 .

7. Show that the Clausius-Duhem inequality in integral form:

q·n ρs
„Z « Z Z Z
d
ρ η dV ≥ ρ η (un − v · n) dA − dA + dV .
dt Ω ∂Ω ∂Ω T Ω T

can be written in differential form as !

q ρs
ρ η̇ ≥ −∇ · + .
T T

Assume that Ω is an arbitrary fixed control volume. Then un = 0 and the derivative can be take inside the integral to give

q·n ρs
Z Z Z Z
∂
(ρ η) dV ≥ − ρ η (v · n) dA − dA + dV .
Ω ∂t ∂Ω ∂Ω T Ω T

Using the divergence theorem, we get

!
q ρs
Z Z Z Z
∂
(ρ η) dV ≥ − ∇ · (ρ η v) dV − ∇· dV + dV .
Ω ∂t Ω Ω T Ω T
33
 Since Ω is arbitrary, we must have !
∂ q ρs
(ρ η) ≥ −∇ · (ρ η v) − ∇ · + .
∂t T T

Expanding out !
∂ρ ∂η q ρs
η+ρ ≥ −∇(ρη ) · v − ρ η (∇ · v) − ∇ · +
∂t ∂t T T
or, !
∂ρ ∂η q ρs
η+ρ ≥ −η ∇ρ · v − ρ ∇η · v − ρ η (∇ · v) − ∇ · +
∂t ∂t T T
or, !
q ρs
„ « „ «
∂ρ ∂η
+ ∇ρ · v + ρ ∇ · v η+ρ + ∇η · v ≥ −∇ · + .
∂t ∂t T T

Now, the material time derivatives of ρ and η are given by

∂ρ ∂η
ρ̇ = + ∇ρ · v ; η̇ = + ∇η · v .
∂t ∂t

Therefore, !
q ρs
(ρ̇ + ρ ∇ · v) η + ρ η̇ ≥ −∇ · + .
T T

From the conservation of mass ρ̇ + ρ ∇ · v = 0. Hence,

!
q ρs
ρ η̇ ≥ −∇ · + .
T T

8. Show that the Clausius-Duhem inequality !

q ρs
ρ η̇ ≥ −∇ · +
T T

can be expressed in terms of the internal energy as

q · ∇T
ρ (ė − T η̇) − σ : ∇v ≤ − .
T

Using the identity ∇ · (ϕ v) = ϕ ∇ · v + v · ∇ϕ in the Clausius-Duhem inequality, we get

! !
q ρs 1 1 ρs
ρ η̇ ≥ −∇ · + or ρ η̇ ≥ − ∇·q−q·∇ + .
T T T T T

Now, using index notation with respect to a Cartesian basis ej ,

!
1 ∂ ` −1 ´ ´ ∂T 1
ej = − T −2
`
∇ = T ej = − 2 ∇T .
T ∂xj ∂xj T

Hence,
1 1 ρs 1 1
ρ η̇ ≥ − ∇·q+ q · ∇T + or ρ η̇ ≥ − (∇ · q − ρ s) + q · ∇T .
T T2 T T T2
Recall the balance of energy

ρ ė − σ : ∇v + ∇ · q − ρ s = 0 =⇒ ρ ė − σ : ∇v = −(∇ · q − ρ s) .

Therefore,
1 1 q · ∇T
ρ η̇ ≥ (ρ ė − σ : ∇v) + q · ∇T =⇒ ρ η̇ T ≥ ρ ė − σ : ∇v + .
T T2 T
Rearranging,
q · ∇T
ρ (ė − T η̇) − σ : ∇v ≤ − .
T

34
 9. For thermoelastic materials, the internal energy is a function only of the deformation gradient and the temperature, i.e., e = e(F , T ). Show that, for
thermoelastic materials, the Clausius-Duhem inequality

q · ∇T
ρ (ė − T η̇) − σ : ∇v ≤ −
T

can be expressed as
q · ∇T
„ « „ «
∂e ∂e
ρ −T η̇ + ρ − σ · F −T : Ḟ ≤ − .
∂η ∂F T
Since e = e(F , T ), we have
∂e ∂e
ė = : Ḟ + η̇ .
∂F ∂η
Therefore,

q · ∇T q · ∇T
„ « „ «
∂e ∂e ∂e ∂e
ρ : Ḟ + η̇ − T η̇ − σ : ∇v ≤ − or ρ −T η̇ + ρ : Ḟ − σ : ∇v ≤ − .
∂F ∂η T ∂η ∂F T

Now, ∇v = l = Ḟ · F −1 . Therefore, using the identity A : (B · C) = (A · C T ) : B, we have

σ : ∇v = σ : (Ḟ · F −1 ) = (σ · F −T ) : Ḟ .

Hence,
q · ∇T
„ «
∂e ∂e
ρ −T η̇ + ρ : Ḟ − (σ · F −T ) : Ḟ ≤ −
∂η ∂F T
or,
q · ∇T
„ « „ «
∂e ∂e
ρ −T η̇ + ρ − σ · F −T : Ḟ ≤ − .
∂η ∂F T

10. Show that, for thermoelastic materials, the balance of energy

ρ ė − σ : ∇v + ∇ · q − ρ s = 0 .

can be expressed as
ρ T η̇ = −∇ · q + ρ s .

Since e = e(F , T ), we have

∂e ∂e
ė = : Ḟ + η̇ .
∂F ∂η
Plug into energy equation to get
∂e ∂e
ρ : Ḟ + ρ η̇ − σ : ∇v + ∇ · q − ρ s = 0 .
∂F ∂η
Recall,
∂e ∂e
=T and ρ = σ · F −T .
∂η ∂F
Hence,
(σ · F −T ) : Ḟ + ρ T η̇ − σ : ∇v + ∇ · q − ρ s = 0 .

Now, ∇v = l = Ḟ · F −1 . Therefore, using the identity A : (B · C) = (A · C T ) : B, we have

σ : ∇v = σ : (Ḟ · F −1 ) = (σ · F −T ) : Ḟ .

Plugging into the energy equation, we have

σ : ∇v + ρ T η̇ − σ : ∇v + ∇ · q − ρ s = 0

or,
ρ T η̇ = −∇ · q + ρ s .

35
11. Show that, for thermoelastic materials, the Cauchy stress can be expressed in terms of the Green strain as

∂e
σ =ρF · · FT .
∂E

Recall that the Cauchy stress is given by

∂e ∂e ∂e
σ=ρ · FT =⇒ σij = ρ FT = ρ Fjk .
∂F ∂Fik kj ∂Fik

The Green strain E = E(F ) = E(U ) and e = e(F , η) = e(U , η). Hence, using the chain rule,

∂e ∂e ∂E ∂e ∂e ∂Elm
= : =⇒ = .
∂F ∂E ∂F ∂Fik ∂Elm ∂Fik

Now,
1 T 1 T 1
E= (F · F − 1 ) =⇒ Elm = (F Fpm − δlm ) = (Fpl Fpm − δlm ) .
2 2 lp 2
Taking the derivative with respect to F , we get

∂F T
„ « „ «
∂E 1 ∂F ∂Elm 1 ∂Fpl ∂Fpm
= · F + FT · =⇒ = Fpm + Fpl .
∂F 2 ∂F ∂F ∂Fik 2 ∂Fik ∂Fik

Therefore,

∂F T
» „ «– » „ «–
1 ∂e ∂F 1 ∂e ∂Fpl ∂Fpm
σ= ρ : · F + FT · · FT =⇒ σij = ρ Fpm + Fpl Fjk .
2 ∂E ∂F ∂F 2 ∂Elm ∂Fik ∂Fik

Recall,
∂A ∂Aij ∂AT ∂Aji
≡ = δik δjl and ≡ = δjk δil .
∂A ∂Akl ∂A ∂Akl
Therefore, » – » –
1 ∂e ` ´ 1 ∂e
σij = ρ δpi δlk Fpm + Fpl δpi δmk Fjk = ρ (δlk Fim + Fil δmk ) Fjk
2 ∂Elm 2 ∂Elm
or, " #
» – „ «T
1 ∂e ∂e 1 ∂e ∂e
σij = ρ Fim + Fil Fjk =⇒ σ= ρ F · +F · · FT
2 ∂Ekm ∂Elk 2 ∂E ∂E
or, "„ #
«T
1 ∂e ∂e
σ= ρF · + · FT .
2 ∂E ∂E

From the symmetry of the Cauchy stress, we have

σ = (F · A) · F T and σ T = F · (F · A)T = F · AT · F T and σ = σ T =⇒ A = AT .

Therefore,
„ «T
∂e ∂e
=
∂E ∂E
and we get
∂e
σ = ρF · · FT .
∂E

12. For thermoelastic materials, the specific internal energy is given by

e = ē(E, η)

where E is the Green strain and η is the specific entropy. Show that

d 1 d 1 1
(e − T η) = −Ṫ η + S : Ė and (e − T η − S : E) = −Ṫ η − Ṡ : E
dt ρ0 dt ρ0 ρ0

where ρ0 is the initial density, T is the absolute temperature, S is the 2nd Piola-Kirchhoff stress, and a dot over a quantity indicates the material time
derivative.
36
 Taking the material time derivative of the specific internal energy, we get

∂ē ∂ē
ė = : Ė + η̇ .
∂E ∂η

Now, for thermoelastic materials,

∂ē ∂ē
T = and S = ρ0 .
∂η ∂E
Therefore,
1 1
ė = S : Ė + T η̇ . =⇒ ė − T η̇ = S : Ė .
ρ0 ρ0
Now,
d
(T η) = Ṫ η + T η̇ .
dt
Therefore,
d 1 d 1
ė − (T η) + Ṫ η = S : Ė =⇒ (e − T η) = −Ṫ η + S : Ė .
dt ρ0 dt ρ0

Also, !
d 1 1 1
S:E = S : Ė + Ṡ : E .
dt ρ0 ρ0 ρ0

Hence,

! !
d d 1 1 d 1 1
ė − (T η) + Ṫ η = S:E − Ṡ : E =⇒ e−T η− S:E = −Ṫ η − Ṡ : E .
dt dt ρ0 ρ0 dt ρ0 ρ0

13. For thermoelastic materials, show that the following relations hold:

∂ψ 1 ∂ψ ∂g 1 ∂g
= Ŝ(E, T ) ; = −η̂(E, T ) ; = Ẽ(S, T ) ; = η̃(S, T )
∂E ρ0 ∂T ∂S ρ0 ∂T

where ψ(E, T ) is the Helmholtz free energy and g(S, T ) is the Gibbs free energy.
Also show that
∂ Ŝ ∂ η̂ ∂ Ẽ ∂ η̃
= −ρ0 and = ρ0 .
∂T ∂E ∂T ∂S

Recall that
ψ(E, T ) = e − T η = ē(E, η) − T η .

and
1
g(S, T ) = −e + T η + S:E.
ρ0
(Note that we can choose any functional dependence that we like, because the quantities e, η, E are the actual quantities and not any particular
functional relations).
The derivatives are
∂ψ ∂ē 1 ∂ψ
= = S; = −η .
∂E ∂E ρ0 ∂T
and
∂g 1 ∂S 1 ∂g
= :E= E; =η.
∂S ρ0 ∂S ρ0 ∂T
Hence,
∂ψ 1 ∂ψ ∂g 1 ∂g
= Ŝ(E, T ) ; = −η̂(E, T ) ; = Ẽ(S, T ) ; = η̃(S, T )
∂E ρ0 ∂T ∂S ρ0 ∂T

From the above, we have

∂2ψ ∂2ψ ∂ η̂ 1 ∂ Ŝ
= =⇒ − = .
∂T ∂E ∂E∂T ∂E ρ0 ∂T

37
 and
∂2g ∂2g ∂ η̃ 1 ∂ Ẽ
= =⇒ = .
∂T ∂S ∂S∂T ∂S ρ0 ∂T
Hence,
∂ Ŝ ∂ η̂ ∂ Ẽ ∂ η̃
= −ρ0 and = ρ0 .
∂T ∂E ∂T ∂S

14. For thermoelastic materials, show that the following relations hold:

∂ê(E, T ) ∂ η̂ ∂ 2 ψ̂
=T = −T
∂T ∂T ∂T 2

and
∂ẽ(S, T ) ∂ η̃ 1 ∂ Ẽ ∂ 2 g̃ ∂ 2 g̃
=T + S: =T 2
+S : .
∂T ∂T ρ0 ∂T ∂T ∂S∂T
Recall,
ψ̂(E, T ) = ψ = e − T η = ê(E, T ) − T η̂(E, T )
and
1 1
g̃(S, T ) = g = −e + T η + S : E = −ẽ(S, T ) + T η̃(S, T ) + S : Ẽ(S, T ) .
ρ0 ρ0
Therefore,
∂ê(E, T ) ∂ ψ̂ ∂ η̂
= + η̂(E, T ) + T
∂T ∂T ∂T
and
∂ẽ(S, T ) ∂g̃ ∂ η̃ 1 ∂ Ẽ
=− + η̃(S, T ) + T + S: .
∂T ∂T ∂T ρ0 ∂T
Also, recall that
∂ ψ̂ ∂ η̂ ∂ 2 ψ̂
η̂(E, T ) = − =⇒ =− ,
∂T ∂T ∂T 2
∂g̃ ∂ η̃ ∂ 2 g̃
η̃(S, T ) = =⇒ = ,
∂T ∂T ∂T 2
and
∂g̃ ∂ Ẽ ∂ 2 g̃
Ẽ(S, T ) = ρ0 =⇒ = ρ0 .
∂S ∂T ∂S∂T
Hence,
∂ê(E, T ) ∂ η̂ ∂ 2 ψ̂
=T = −T
∂T ∂T ∂T 2

and
∂ẽ(S, T ) ∂ η̃ 1 ∂ Ẽ ∂ 2 g̃ ∂ 2 g̃
=T + S: =T 2
+S : .
∂T ∂T ρ0 ∂T ∂T ∂S∂T

15. For thermoelastic materials, show that the balance of energy equation

ρ T η̇ = −∇ · q + ρ s

can be expressed as either

ρ ∂ Ŝ
ρ Cv Ṫ = ∇ · (κ · ∇T ) + ρ s + T : Ė
ρ0 ∂T
or !
1 ∂ Ẽ ρ ∂ Ẽ
ρ Cp − S: Ṫ = ∇ · (κ · ∇T ) + ρ s − T : Ṡ
ρ0 ∂T ρ0 ∂T

where
∂ê(E, T ) ∂ẽ(S, T )
Cv = and Cp = .
∂T ∂T
If the independent variables are E and T , then

∂ η̂ ∂ η̂
η = η̂(E, T ) =⇒ η̇ = : Ė + Ṫ .
∂E ∂T
38
 On the other hand, if we consider S and T to be the independent variables

∂ η̃ ∂ η̃
η = η̃(S, T ) =⇒ η̇ = : Ṡ + Ṫ .
∂S ∂T

Since !
∂ η̂ 1 ∂ Ŝ ∂ η̂ Cv ∂ η̃ 1 ∂ Ẽ ∂ η̃ 1 1 ∂ Ẽ
=− ; = ; = ; and = Cp − S:
∂E ρ0 ∂T ∂T T ∂S ρ0 ∂T ∂T T ρ0 ∂T

we have, either
1 ∂ Ŝ Cv
η̇ = − : Ė + Ṫ
ρ0 ∂T T
or !
1 ∂ Ẽ 1 1 ∂ Ẽ
η̇ = : Ṡ + Cp − S: Ṫ .
ρ0 ∂T T ρ0 ∂T

The equation for balance of energy in terms of the specific entropy is

ρ T η̇ = −∇ · q + ρ s .

Using the two forms of η̇, we get two forms of the energy equation:

ρ ∂ Ŝ
− T : Ė + ρ Cv Ṫ = −∇ · q + ρ s
ρ0 ∂T

and
ρ ∂ Ẽ ρ ∂ Ẽ
T : Ṡ + ρ Cp Ṫ − S: Ṫ = −∇ · q + ρ s .
ρ0 ∂T ρ0 ∂T
From Fourier’s law of heat conduction
q = −κ · ∇T .
Therefore,
ρ ∂ Ŝ
− T : Ė + ρ Cv Ṫ = ∇ · (κ · ∇T ) + ρ s
ρ0 ∂T
and
ρ ∂ Ẽ ρ ∂ Ẽ
T : Ṡ + ρ Cp Ṫ − S: Ṫ = ∇ · (κ · ∇T ) + ρ s .
ρ0 ∂T ρ0 ∂T
Rearranging,
ρ ∂ Ŝ
ρ Cv Ṫ = ∇ · (κ · ∇T ) + ρ s + T : Ė
ρ0 ∂T

or,
!
1 ∂ Ẽ ρ ∂ Ẽ
ρ Cp − S: Ṫ = ∇ · (κ · ∇T ) + ρ s − T : Ṡ .
ρ0 ∂T ρ0 ∂T

16. For thermoelastic materials, show that the specific heats are related by the relation
!
1 ∂ Ŝ ∂ Ẽ
Cp − Cv = S−T : .
ρ0 ∂T ∂T

Recall that
∂ê(E, T ) ∂ η̂
Cv := =T
∂T ∂T
and
∂ẽ(S, T ) ∂ η̃ 1 ∂ Ẽ
Cp := =T + S: .
∂T ∂T ρ0 ∂T
Therefore,
∂ η̃ 1 ∂ Ẽ ∂ η̂
Cp − Cv = T + S: −T .
∂T ρ0 ∂T ∂T
Also recall that
η = η̂(E, T ) = η̃(S, T ) .
39
 Therefore, keeping S constant while differentiating, we have

∂ η̃ ∂ η̂ ∂E ∂ η̂
= : + .
∂T ∂E ∂T ∂T

Noting that E = Ẽ(S, T ), and plugging back into the equation for the difference between the two specific heats, we have

∂ η̂ ∂ Ẽ 1 ∂ Ẽ
Cp − Cv = T : + S: .
∂E ∂T ρ0 ∂T

Recalling that
∂ η̂ 1 ∂ Ŝ
=−
∂E ρ0 ∂T
we get
!
1 ∂ Ŝ ∂ Ẽ
Cp − Cv = S−T : .
ρ0 ∂T ∂T

17. For thermoelastic materials, show that the specific heats can also be related by the equations

1 ∂2ψ 1 T ∂E
„ « „ «
∂E ∂E ∂E ∂E ∂S ∂E
Cp − Cv = S: + : : = S: + : : .
ρ0 ∂T ∂T ∂E∂E ∂T ρ0 ∂T ρ0 ∂T ∂E ∂T

Recall that
∂ψ
S = ρ0 = ρ0 f (E(S, T ), T ) .
∂E
Recall the chain rule which states that if
g(u, t) = f (x(u, t), y(u, t))

then, if we keep u fixed, the partial derivative of g with respect to t is given by

∂g ∂f ∂x ∂f ∂y
= + .
∂t ∂x ∂t ∂y ∂t

In our case,
u = S, t = T, g(S, T ) = S, x(S, T ) = E(S, T ), y(S, T ) = T, and f = ρ0 f .

Hence, we have
S = g(S, T ) = f (E(S, T ), T ) = ρ0 f (E(S, T ), T ) .

Taking the derivative with respect to T keeping S constant, we have

0 1
2 3
∂g ∂S7
 6 ∂f ∂E ∂f ∂T7
7
=  = ρ0 6
4 ∂E : ∂T + ∂T ∂T 5
7
∂T ∂T
 

or,
∂f ∂E ∂f
0 = : + .
∂E ∂T ∂T
Now,
∂ψ ∂f ∂2ψ ∂f ∂2ψ
f = =⇒ = and = .
∂E ∂E ∂E∂E ∂T ∂T ∂E
Therefore,
∂2ψ ∂2ψ
„ « „ «
∂E ∂ ∂ψ ∂E ∂ ∂ψ
0 = : + = : + .
∂E∂E ∂T ∂T ∂E ∂E ∂E ∂T ∂T ∂E
Again recall that,
∂ψ 1
= S.
∂E ρ0
Plugging into the above, we get
∂2ψ ∂E 1 ∂S 1 ∂S ∂E 1 ∂S
0 = : + = : + .
∂E∂E ∂T ρ0 ∂T ρ0 ∂E ∂T ρ0 ∂T

40
 Therefore, we get the following relation for ∂S/∂T :

∂S ∂2ψ ∂E ∂S ∂E
= −ρ0 : =− : .
∂T ∂E∂E ∂T ∂E ∂T

Recall that
1
„ «
∂S ∂E
Cp − Cv = S−T : .
ρ0 ∂T ∂T
Plugging in the expressions for ∂S/∂T we get:

1 ∂2ψ 1
„ « „ «
∂E ∂E ∂S ∂E ∂E
Cp − Cv = S + T ρ0 : : = S+T : : .
ρ0 ∂E∂E ∂T ∂T ρ0 ∂E ∂T ∂T

Therefore,
1 ∂2ψ 1 T
„ « „ «
∂E ∂E ∂E ∂E ∂S ∂E ∂E
Cp − Cv = S: +T : : = S: + : : .
ρ0 ∂T ∂E∂E ∂T ∂T ρ0 ∂T ρ0 ∂E ∂T ∂T
Using the identity (A : B) : C = C : (A : B), we have

1 ∂2ψ 1 T ∂E
„ « „ «
∂E ∂E ∂E ∂E ∂S ∂E
Cp − Cv = S: +T : : = S: + : : .
ρ0 ∂T ∂T ∂E∂E ∂T ρ0 ∂T ρ0 ∂T ∂E ∂T

18. Consider an isotropic thermoelastic material that has a constant coefficient of thermal expansion and which follows the St-Venant Kirchhoff model,
i.e,
αE = α 1 and C = λ 1 ⊗ 1 + 2µ I

where α is the coefficient of thermal expansion and 3 λ = 3 K − 2 µ where K, µ are the bulk and shear moduli, respectively.
Show that the specific heats related by the equation

1 ˆ
α tr(S) + 9 α2 K T .
˜
Cp − Cv =
ρ0

Recall that,
1 T
Cp − Cv = S : αE + αE : C : αE .
ρ0 ρ0
Plugging the expressions of αE and C into the above equation, we have

1 T
Cp − Cv = S : (α 1 ) + (α 1 ) : (λ 1 ⊗ 1 + 2µ I) : (α 1 )
ρ0 ρ0
α α2 T
= tr(S) + 1 : (λ 1 ⊗ 1 + 2µ I) : 1
ρ0 ρ0
α α2 T
= tr(S) + 1 : (λ tr(1 ) 1 + 2µ 1 )
ρ0 ρ0
α α2 T
= tr(S) + (3 λ tr(1 ) + 2µ tr(1 ))
ρ0 ρ0
α 3 α2 T
= tr(S) + (3 λ + 2µ)
ρ0 ρ0
α tr(S) 9 α2 K T
= + .
ρ0 ρ0

Therefore,
1 ˆ
α tr(S) + 9 α2 K T .
˜
Cp − Cv =
ρ0

41
References
[1] T. W. Wright. The Physics and Mathematics of Adiabatic Shear Bands. Cambridge University Press,
Cambridge, UK, 2002.

[2] R. C. Batra. Elements of Continuum Mechanics. AIAA, Reston, VA., 2006.

[3] G. A. Maugin. The Thermomechanics of Nonlinear Irreversible Behaviors: An Introduction. World Scientific,
Singapore, 1999.

[4] R. M. Brannon. The Consistent Kinetics Porosity (CKP) model: A theory for the mechanical behavior of
moderately porous solids. Sandia Report SAND2000-2696, Sandia National Laboratories, Albuquerque, New
Mexico, 2000.

[5] M. L. Wilkins. Computer Simulation of Dynamic Phenomena. Springer-Verlag, Berlin, 1999.

[6] M. A. Zocher, P. J. Maudlin, S. R. Chen, and E. C. Flower-Maudlin. An evaluation of several hardening

models using Taylor cylinder impact data. In Proc. , European Congress on Computational Methods in
Applied Sciences and Engineering, Barcelona, Spain, 2000. ECCOMAS.

[7] P. J. Maudlin and S. K. Schiferl. Computational anisotropic plasticity for high-rate forming applications.
Comput. Methods Appl. Mech. Engrg., 131:1–30, 1996.

[8] P. Perzyna. Constitutive equations for thermoinelasticity and instability phenomena in thermodynamic flow
processes. In Stein E., editor, Progress in Computational Analysis of Inelastic Structures: CISM Courses and
Lectures No. 321, pages 1–78. Springer-Verlag-Wien, New York, 1993.

[9] S. G. Bardenhagen and E. M. Kober. The Generalized Interpolation Material Point Method. Comp. Model.
Eng. Sci., 5(6):477–495, 2004.

[10] D. Sulsky, S. Zhou, and H. L. Schreyer. Application of a particle-in-cell method to solid mechanics.
Computer Physics Communications, 87:236–252, 1995.

[11] M. D. Greenberg. Foundations of Applied Mathematics. Prentice-Hall Inc., Englewood Cliffs, New Jersey,
1978.

[12] T. Belytschko, W. K. Liu, and B. Moran. Nonlinear Finite Elements for Continua and Structures. John Wiley
and Sons, Ltd., New York, 2000.

[13] M. E. Gurtin. An Introduction to Continuum Mechanics. Academic Press, New York, 1981.

42