Computational Particle Mechanics, 6(1), 133-144 (2019) manuscript No.

(will be inserted by the editor)

Modeling heat flow across material interfaces and cracks

using the material point method
John A. Nairn

Accepted: July 2018

Abstract Heat conduction through an object with ma- volumetric expansion or compression, plastic, viscoelas-
terial interfaces or cracks is influenced by heat flow tic, or damage energy dissipation, or contact with fric-
across those discontinuities. This paper presents a nu- tion and to fully track thermodynamic state variables
merical particle method for modeling such heat flow (e.g., entropy, enthalpy, and free energy).
coupled to computational mechanics all within the ma- Standard MPM conduction methods extrapolate par-
terial point method (MPM). In brief, MPM models con- ticle temperature to a single temperature field, thereby
tact and cracks by extrapolating multiple velocity fields solving for a global temperature field. When a problem
to a grid. To model interfacial heat flow, MPM should has material interfaces or cracks, however, this single-
similarly extrapolate multiple temperature fields. Inter- temperature-field approach cannot model the influence
faces nodes that “see” more than one temperature field of those discontinuities on heat flow. This situation is
modify their heat flow to reflect interfacial physics. For analogous to MPM contact methods. When an MPM
example, interfaces in contact may transfer heat by per- code uses a single velocity field, the calculations can
fect conduction while separated interfaces may block model only stick contact. Extending MPM to model
heat flow or cause reduced heat flow by convection. Af- contact by friction [2] or an imperfect interface [10, 11]
ter some validation examples, two real-world examples requires use of multimaterial velocity fields. Nodes that
consider cooling an ingot within a crucible where cool- “see” only a single material proceed by conventional
ing causes the ingot to lose contact with the crucible MPM while “interface” nodes that see multiple materi-
walls and thermal imaging of cracks within an opaque als adjust nodal momenta to reflect some modeled con-
solid. tact mechanics.
Keywords Heat conduction, material point method, This paper adopts a similar approach for heat flow
explicit cracks, material contact calculations by extrapolating multiple temperature fields.
Single field nodes proceed by conventional conduction
methods while interface nodes adjust nodal heat flows
1 Introduction to model realistic heat flow across an interface. The
“Heat flow equations” section first describes global con-
A straightforward addition to particle-based, material duction analysis and then describes how to modify that
point method (MPM) modeling is to couple the me- analysis for multiple temperature fields. The methods
chanics analysis to thermal conduction. Coupling is done are applied both to multimaterial MPM and to MPM
by tracking particle temperature, extrapolating tem- with explicit cracks [9]. The methods were verified by
perature to a background grid, and solving the heat flow simple calculations across a contact interface. Applica-
equation on the grid along with the momentum equa- tions for the method are demonstrated with two real-
tion. Having such a feature allows MPM to model heat- world examples for cooling an ingot in contact with a
ing mechanisms caused by physical phenomena such as crucible where thermal shrinkage may cause the ingot to
Oregon State University, Wood Science & Engineering, lose contact with the crucible walls and thermal imag-
112 Richardson Hall, Corvallis, OR 97330, USA E-mail: ing of opaque objects with internal cracks. Although
[Link]@[Link] MPM provides no specific advantages over other nu-
2 J. A. Nairn

merical methods (such as finite element analysis) for heat capacity) and χp (x) is a particle basis function
solution of heat transport equations, MPM may have for particle p (which is typically 1 within the particle’s
advantages for modeling explicit cracks [9] and com- domain and zero elsewhere [3]). Next expand the weight
plex contact mechanics [11, 14]. The methods in the function and its gradient using standard, isoparametric
paper allow MPM simulations with cracks and inter- grid shape functions, Ni (x):
faces to account for effects of those discontinuities on X X
heat transport. w(x) = wi Ni (x) and ∇w(x) = wi ∇Ni (x) (8)
i i

2 Heat flow equations where wi are nodal values of w(x) on the grid. After
substituting all expansions, the weak form equation be-
2.1 Global analysis comesZ
X
wi Ni (x)q(x) · n̂ dS =
A common addition to MPM codes is a global heat i δV

conduction option that ignores material interfaces and Z (X X

cracks. Despite its use in many MPM codes, a full gen- + [wi ∇Ni (x) · qp χp (x)]
V i p
eralized interpolation material point (GIMP) [3] deriva- XX
tion is rarely in the literature [15]; it is given here. The + q̇s,p χp (x)wi Ni (x)
heat conduction equation is i p
)
∂T XX ∂Tp
ρC + ∇ · q(x) = q̇s (x) (1) − χp (x)wi Ni (x)ρp Cp dV (9)
∂t i p
dt
where ρ is density, C is heat capacity (per unit mass), T Exploiting the fact that w(x) is arbitrary, this equation
is temperature, q is heat flux (per unit area), and q̇s is transforms to a system of equations for node i:
X ∂Tp X X
volumetric heat source rate. For heat conduction, heat Mp C p Spi = Vp qp · Gpi + Vp q̇s,p Spi
flux is q = −k∇T where k is the thermal conductivity p
dt p p
tensor. Solving this equation in the MPM weak form
Z
gives − Ni (x)q(x) · n̂ dS (10)
δV
Z   where Mp and Vp are particle mass and volume and Spi
∂T
ρC + ∇ · q − q̇s (x) w(x) dV = 0 (2) and Gpi are
Z GIMP shape functions [3]:
V ∂t 1
Spi = χp (x)Ni (x) dV (11)
where V is total volume and w(x) is an arbitrary weight- Vp V
Z
ing function. Using the vector identity: 1
Gpi = χp (x)∇Ni (x) dV (12)
Vp V
w(x)∇ · q(x) = ∇ · (w(x)q(x)) − ∇w(x) · q(x) (3) Allowing particle volume and heat capacity to change
each time step (superscripted with (n)), the thermal en-
and the divergence theorem, the weak form equation
ergy on node i (with SI units J) can be defined as:
becomes:
Z
∂T (n)
X (n)
w(x)ρ(x)C(x) − ∇w(x) · q(x) τT i = Mp Cp(n) Tp(n) Sip (13)
V ∂t p
! Z
− q̇s (x)w(x) dV + w(x)q(x) · n̂ dS = 0 (4) The MPM thermal conduction equation becomes:
δV
(n)
where δV is the border of V and n̂ is a surface normal dτT i (n) (n)
vector. By GIMP methods [3], particle quantities are = qi + qi,q (14)
dt
expanded in a particle basis to get:
∂T X ∂Tp where X  
ρ(x)C(x) = ρp Cp χp (x) (5) (n)
qi =
(n) (n)
Vp(n) qp · Gip + q̇s,p Sip (15)
dt p
dt
p
X Z
q̇s (x) = q̇s,p χp (x) (6) (n)
qi,q =− Ni (x)q(x) · n̂ dS (16)
p
δV
are nodal thermal energy flows (with SI units J/sec or
X
q(x) = qp χp (x) (7) (n) (n)
p Watts). qi is internal thermal flow while qi,q is ther-
where subscript p denotes a particle property (note: mal flow at boundaries due to flux boundary conditions.
Cp is particle heat capacity and not constant-pressure Unlike the corresponding MPM momentum calculation
Contact Heat Flow 3

for force that depends on tracked particle stress, par- temperature must be updated by extrapolating temper-
ticle heat flux is not tracked on the particle; instead ature velocity instead of temperature, to the particle.
it is calculated on each time step using shape function A temperature extrapolation leads to artificial conduc-
gradients: tion, even if the material’s conductivity is zero [4].
X (n) (n) The conduction solution is easily coupled to stan-
qp = −k(n)p ∇Tp
(n)
= −k(n)
p Ti Gip (17) dard MPM mechanics analysis. Coupling occurs two
i
ways - through thermal expansion or through heat gen-
where kp is particle thermal conductivity and nodal erated by various mechanisms. Thermal expansion is
(n)
temperature, Ti , is found from: coupled by evaluating temperature increment on each
particle for input to constitutive law calculations that
(n)
(n) τT i (n)
X (n) account for thermal expansion. Because of the way par-
Ti = (n)
where ci = Mp Cp(n) Sip (18)
ci p
ticle temperatures update, some simulations exhibit vari-
(n)
ations in Tp within a cell. Despite these variations, the
is nodal heat capacity (SI units J/K) extrapolated from temperature field on the grid remains smooth and is a
particles. The temperature update on each node is better description of the current temperature field. As
(n) (n)
(n+1) (n) qi + qi,q (n) (n) a consequence, temperature increment on each particle
Ti = Ti + (n)
∆t = Ti + vT i ∆t (19) is best calculated from
ci
where nodal temperature “velocity” (SI units K/sec) is (n+1) (n) (n)
X (n) (n)
∆Tp = Tg→p − Tg→p where Tg→p = Ti Spi (25)
(n) (n) i
(n) qi + qi,q
vT i = (n)
(20)
ci (n)
rather then from Tp . To make this calculation possi-
(n) (n) (n)
Once temperature velocities are found, each parti- ble, particles should track both Tp and Tg→p . Tp is
(n)
cle’s temperature updates by: used in conduction equations while Tg→p is used to find
(n)
∆Tp or to implement any other feature that depends on
Tp(n+1) = Tp(n) + vT p ∆t (21) particle temperature.
(n)
Heat generated by constitutive laws may be adia-
where vT p is temperature velocity extrapolated to the batic heating (e.g., adiabatic deformation of an elastic
particle. Two options are possible. The first extrapo- material changes temperature by dTad = (M · dε)/ρ
lates the temperature velocity using where M is the stress-temperature tensor and dε is
(n)
vT p =
X (n) (n)
vT i Spi (22) the strain increment [5]) or any energy that is dissi-
i
pated as heat (e.g., plastic or viscoelastic energy dis-
sipation). These heat terms can be converted to tem-
while the second extrapolates heat flux and divides that perature changes (depending on particle’s current heat
result by particle heat capacity: capacity), accumulated during constitutive law calcula-
(n)
(n) (n) tions as dTp,ad , and then added to particle temperature
(n)∗ 1 X qi + qi,q (n)
vT p = (n)
Spi (23) during MPM particle updates. In effect, this approach
ρp C p i vi is solving an adiabatic particle heat equation:
where (n)
dTp
(n)
X (n) ρp C p = q̇s (26)
vi = Vp Sip (24) dt
p
with solution
is volume extrapolated to the node. These two methods
q̇s ∆t (n)
are identical when all particles have the same mass, den- Tp(n+1) = Tp(n) + = Tp(n) + dTp,ad (27)
sity, and heat capacity, but would differ in composite ρp Cp
materials and when modeling phase transitions where
Other heat generation might be frictional contact at
heat capacity varies in the melting region. In general,
material interfaces and cracks [14] or crack tip heating
the results seem similar, but unpublished results with
during crack propagation. These heats are often found
phase transitions suggest the first method (extrapola- (n)
on the grid and implemented by adding to qi in con-
tion of temperature rate) is more stable. This paper
duction calculations
uses that method (and actually either method would
work because all examples used particles with the same (n) (n)
X
qi (total) = qi + qj,i (28)
mass, density, and heat capacity). Note that particle j
4 J. A. Nairn

where the sum is over events j that each generates heat temperature rates for material α is:
qj,i on node i (SI units Watts). (n) (n) (n) (n)
The above derivation is for Cartesian coordinates. (n+1)
X qi,α + ∆qi,α + φi,α qi,q
Tp∈α = Tp(n) +∆t (n)
Spi (35)
Extension to axisymmetric MPM is trivial and details i ci,α
are provided in Ref. [15]. The only changes needed for
axisymmetry are to use axisymmetric Spi and Gpi shape and the back extrapolation to find temperature incre-
functions and to use particle mass (Mp ) and volume ment is
(Vp ) on per-radian basis. (n)
X (n) (n)
Tg→p∈α = Ti,α Spi (36)
i

2.2 Heat flow across material interfaces The implementation of MPM heat conduction anal-
ysis that accounts for material interfaces is now reduced
(n) (n)
Extending MPM to handle heat flow at material inter- to determining ∆qi,α and φi,α at each interface node.
faces requires each material to extrapolate to indepen-
dent temperature fields on the nodes and then to im-
2.3 Contact calculations
plement contact heat flow physics at all interface nodes
(defined as nodes that “see” temperature fields from
When MPM calculations are in multimaterial mode
more than one material). We define material specific
[2, 7, 11], the standard approach to contact is to adjust
quantities from sums that include only particles of ma-
momenta on all nodes in contact. In brief, whenever
terial α:X
(n) (n) a node has more than one material, mechanics calcu-
τT i,α = Mp Cp(n) Tp(n) Sip (29) lations determine if they are in contact [11]. If they
p∈α
are in contact, momenta are changed to reflect contact
(n) (n)
X
ci,α = Mp Cp(n) Sip (30) physics (e.g., frictional contact [2, 7] or an imperfect
p∈α interface [11]). When this mechanics analysis is cou-
 
(n) (n) (n) (n) pled to heat flow analysis, these contact calculations
X
qi,α = Vp(n) q0(n)
p · Gip + q̇s,p Sip (31) (n) (n)
p∈α should also calculate ∆qi,α and φi,α on the same con-
0(n)
where (importantly) qp is heat flow that would occur tact nodes. Furthermore, while contact forces can be
on particle p in the absence of interactions with other skipped when interface are separated, heat flow calcu-
materials. It is calculated from the temperature field lations are needed for interfaces both in contact and
for the material type of particle p: separated. Determination of contact or separation can
X τT(n) be used to implement different heat flow characteristics
0(n) i,α
X (n) (n) (n)
qp∈α = −k(n)
p T i,α Gip = −k (n)
p (n)
Gip (32) for each state. For example, interfaces in contact might
i i c i,α transfer heat by conduction, while separated interfaces
For any number of materials, global and material ex- might use convection.
trapolations are related by: To start, consider perfect conduction at an interface
such that multimaterial conduction analysis should re-
(n) (n) (n) (n)
X X
τT i = τT i,α and ci = ci,α (33) vert to the global analysis that ignores interfaces. In this
α α limit, the particle updates must be identical. Equating
Eq. (21) to Eq. (35) and solving for perfect conduction
In contrast,
P global and material heat flows may differ ∗(n)
(n) (n) (n) (n)
(qi 6= α qi,α ) because qi and qi,α ’s are based on heat flow change, ∆qi,α , gives
(n) 0(n)
different particle heat flows (qp vs. qp ) that are cal- (n)
ci,α
(n)
ci,α
!
∗(n) (n) (n) (n)
culated from different temperature field gradients. ∆qi,α = q
(n) i
− qi,α + (n)
− φi,α qi,q (37)
The governing equation on node i for material α is ci ci
(n) (n) (n)
(n) If we assume φi,α = ci,α /ci (i.e., assume that ex-
dτT i,α (n) (n) (n) (n)
= qi,α + ∆qi,α + φi,α qi,q (34) ternally applied flux is logically spread over available
dt materials according to their heat capacity fraction on
(n) the node), this term becomes:
where ∆qi,α is any additional heat flow that must be
added to account for material interfaces at node i and ∗(n)
∆qi,α = φi,α qi
(n) (n)
− qi,α (38)
(n)
φi,α is the fraction of any heat flux boundary condi-
tion on node i that must be applied to material α. The Next consider an interface with exactly two mate-
temperature update for particle p of material α using rials – α and β – having contact by convection with
Contact Heat Flow 5

convection coefficient h (in W/(m2 -K)). The convection [11], use the contact state to select contact prop-
heat flow applied to each material is: erties (e.g., assume heat flow by convection with
  convection coefficient h or equilibrated heat flow by
(n) (n) (n) (n) (n)
∆qi,α = hai Ti,β − Ti,α and ∆qi,β = −∆qi,α (39) conduction), and then repeat following steps 2 and
3 for each material (α) at the node (i).
where ai is interfacial contact area on the grid for node ∗(n)
2. Find ∆qi,α for heat flow under the perfect conduc-
i. In other words, heat flow induced by convection at
tion limit.
an interface is added to each material and total extra
3. If the contact state is using equilibrated conduction,
heat flow is zero. (n) ∗(n)
set ∆qi,α = ∆qi,α . If instead it is using convec-
The contact area can be calculated by the same (n)
methods used when implementing MPM contact meth- tion, find ai and ∆qi,α by Eq. (42). This heat flow,
ods that depend on contact area as described in (au- however, must not exceed conduction heat flow. If
(n) ∗(n) (n)
thor?) [11]; the needed area is: |∆qi,α | < |∆qi,α |, then use ∆qi,α ; otherwise set
(n) ∗(n)
q ∆qi,α = ∆qi,α .
(n) (n) (n) (n)
2(vi,α + vi,β ) min(vi,α , vi,β ) 4. If mechanical analysis determines interfaces in con-
ai = (40) tact by friction, the frictional sliding can be con-
t⊥
verted to heat by adding a calculated qf riction to
where qi,q (i.e., by adding φi,α qf riction to heat flow in each
(n)
vi,α =
X
Vp(n) Sip
(n)
(41) material’s temperature field).
(n)
p∈α Note that calculation of ∆qi,α requires knowledge of
(n) (n) (n)
is material volume extrapolated to the grid and t⊥ is both qi , which depends on qp , and qi,α , which de-
0(n)
an effective thickness of the contacting volume. For a pends on qp . As a consequence, MPM code to imple-
regular grid with equal element sides, ∆x = ∆y = ∆z, ment interfacial heat flow must extrapolate both global
(n)
t⊥ = ∆x or is equal to the constant cell size. For el- and material temperature fields to calculate both qp
ements with rectangular elements, t⊥ is needed to ac- 0(n)
and qp on each particle.
count for interfaces oriented in different grid directions.
More explanation and a method for finding t⊥ are given
in Ref. [11]. Importantly, ai is an effective contact area 2.4 Heat flow across cracks
that reduces to an area as a function of distance from
node i to the interface. This scaling is crucial for grid Analysis of heat flow across cracks is nearly identical to
independence of contact results as demonstrated in the heat flow at material interfaces except separate temper-
“Results and discussion” section. ature fields for each material are replaced by separate
A robust numerical implementation must allow for temperature fields for each side of the crack. In MPM
three or more materials on a single node. To handle with cracks (which is called CRAMP [9]), each parti-
this situation, replace material β with a virtual ma- cle node pair is assigned a crack velocity field, v(p, i)
terial that lumps all materials besides material α or = 1 or 0 depending whether a line from the particle to
(n) (n) (n) (n) (n) (n) (n) (n) (n)
τT i,β = τT i −τT i,α , ci,β = ci −ci,α , Ti,β = τT i,β /ci,β , the node crosses a crack (1) or does not cross a crack
(n) (n) (n) (0). For a single crack, each node will have at most two
and vi,β = vi − vi,α . Substituting into Eq. (39) and
velocity fields (0 and 1). The CRAMP method can be
eliminating all β terms, a general convection flow added
extended to handle two interacting cracks by allowing
to material α becomes
up to four velocity fields on each node (or v(p, i) = 0 to
(n) (n)
(n) Ti − Ti,α 3) corresponding to lines that cross no cracks (0), one
∆qi,α = hai (42) crack (1), a second crack (2), or both cracks (3) [13].
1 − φi,α
When modeling heat flow across cracks, the crack
This change is applied (individually) to each material velocityXfield specific thermal terms become:
(n) (n)
on the node. For two materials, this approach reduces τT i,j = Mp Cp(n) Tp(n) Sip δj,v(p,i) (43)
exactly to Eq. (39) and total added heat flow is zero. p
(n) (n)
X
For more than two materials, it gives a reasonable result ci,j = Mp Cp(n) Sip δj,v(p,i) (44)
(n)
although the sum of ∆qi,α may not be zero. p
The full contact heat flow algorithm is:
 
(n) (n) (n) (n)
X
qi,j = Vp(n) q00(n)
p · Gip + qs,p Sip δj,v(p,i) (45)
1. After updating momenta on the nodes, evaluate the p
velocity fields at each node with more then one ma- (n) (n)
X
vi,j = Vp(n) Sip δj,v(p,i) (46)
terial, determine whether or not they are in contact p
6 J. A. Nairn

where j is crack temperature field j on node i and δ over material temperature fields according to their
(n)
is the Kronecker delta. The crack heat flow term, qi,j , thermal mass fraction within the crack velocity field).
00(n)
includes yet another particle heat flow, qp , that is For all interface nodes, particle updates use the ma-
calculated each time step using terial temperature field (see Eq. (35)). For this update
X τT(n) to work on nodes having multiple crack velocity fields
i,v(p,i) (n)
q00(n)
p = −k(n)
p (n)
Gip (47) with only a single material, the total crack heat flow in
i ci,v(p,i) such fields, j, should be copied to that one material, α
(n) (n)
(i.e., qi,α (n) = qi,j + ∆qi,j ). Finally, for MPM code to
In other words, it extrapolates temperature gradient
implement combined interfacial and crack heat flow, it
from each node using the temperature appropriate for
must extrapolate global, material, and crack tempera-
that particle-node pair. (n) 0(n) 00(n)
ture fields to calculate qp , qp , and qp on each
The crack contact calculations determine if crack
particle.
surfaces are in contact and heat flow calculations are
applied to each crack temperature field instead of each
material temperature field. Crack heat flow calculations
are identical to material calculations in section 2.3 ex- 3 Results and discussion
cept they use crack terms instead of material terms.
The particle updates in Eqs. (35) and (36) replace α The algorithm presented is fully 3D and could be added
with the appropriate velocity field v(p, i). To imple- to any code that models 3D contact mechanics and/or
ment these calculations, the extrapolations must eval- 3D cracks. The verification examples given here were
(n)
uate both qp and qp
00(n)
on each particle. all 2D or axisymmetric. The first problem was a 2D
simulation, but models a 1D problem. Consider a strip
of length L = 100 mm from x0 = −50 mm to x1 =
2.5 Combining cracks and material interfaces +50 mm and width 20 mm with MPM background cell
size of 2.5 × 2.5 mm (and, like all simulations in this
To combine heat flow across both cracks and material paper, four particles per cell). The strip was comprised
interfaces, temperature fields have to be arranged in a of two separate isotropic materials, but with identi-
hierarchical structure. The approach used here was to cal (arbitrarily-selected) thermal properties: k = 2000
allow each node to have multiple crack velocity fields W/(m·K), C = 1000 J/(kg·K), and ρ = 1 g/cm3 . All
(up to two for a single crack or up to four to handle calculations were done in the MPM code OSParticu-
two interacting cracks). Each crack velocity field may las [12], which fully couples mechanical and thermal
have 1 to m material temperature fields, where m is the calculations. Because thermal conduction is typically
number of materials in the simulation. much slower than stress waves, these calculations used a
In this arrangement, an analysis may encounter three low-modulus, zero-expansion material (E = 0.01 MPa,
types of nodes requiring additional heat flow calcula- ν = 0.33, and thermal expansion coefficient = 0) to
tions — nodes with multiple materials within a single allow a larger time step for explicit integrations.
crack velocity field, nodes with multiple crack velocity At time zero, all particles were set to temperature
fields each having only a single material, and nodes with T0 = 0◦ C and boundary conditions set T = T0 = 0◦ C
multiple crack velocity fields containing multiple mate- at x0 and T = T1 = 100◦ C at x1 . The interface was at
rials. The first two are handled exactly as described in x = 0 and was modeled as conduction when in contact,
sections 2.2 and 2.4. The last one requires special treat- but convection (with convection coefficient h) when sep-
ment: arated. To induce contact or separation, the two mate-
rials were either pushed together or pulled apart by
1. Heat flow at material interfaces are handled first.
0.5 mm before the heat flow reached the interface. A
Because materials are within a crack velocity field,
(n) (n) Fourier series solution to this problem for conduction
these calculations replace global values (Ti , ci ,
(n) in the absence of an interface is [6]:
and qi ) with the corresponding crack temperature
(n) (n) (n)
field values (Ti,j , ci,j , and qi,j ) and then proceed X (−1)n
!
as described in section 2.2. −λ2n t
T (x, t) = T0 + ∆T ξ+2 e sin(nπξ)
2. Heat flow at cracks are handled second. The change n
nπ
to the crack temperature field heat flow is calculated (48)
exactly as described in section 2.4. When done, how-
(n) (n) (n)
ever, ci,α ∆qi,j /ci,j is added to each material tem- p = T1 − T0 , ξ = (x − x0 )/L, and λn =
where ∆T
perature field (i.e., the added heat flow is spread (nπ/L) k/(ρC).
Contact Heat Flow 7

100 is understood by looking at expected interface nodes.

Contact Figure 3 shows material points around an interface near
80 grid line 0 when in contact (Fig. 3A) or separated (Fig.
3B). The circled nodes show all nodes that would be in-
Temperature (C)

60 MPM
terface nodes when using GIMP shape functions (and
Analytical
non-GIMP methods used in early MPM papers [16]
40
should be avoided because they activate fewer nodes).
When in contact, the nodes at -1, 0 and +1 are all inter-
20
face nodes that adjust their heat flow and this arrange-
0
ment is able to exactly revert to perfect conduction for
-50 -40 -30 -20 -10 0 10 20 30 40 50 large h. The nodes at ±2 are not interface nodes, but
Position (mm)
they are sufficiently far from the interface that global
Fig. 1 Temperature profile in a bar with two materials in (n) 0(n)
(qp ) and material (qp ) heat flows would be identi-
contact at x = 0 mm. All material points started at T =
0◦ C and the results are after 500 ms with T set to 0◦ C at cal. Thus, when in contact all nodes can exactly recover
x = −50 mm and 100◦ C at x = +50 mm . The solid black perfect conduction for large h.
line is the MPM result. The dotted red line is the analytical In contrast, when separated (Fig. 3B), only the nodes
solution. at 0 are interface nodes. Although those nodes can re-
vert to conduction for very high h, the nodes at ±1
100 have only a single material field. Because those nodes
Separated are close to an interface, they will have a slight discrep-
80 (n) 0(n)
ancy between their global (qp ) and material (qp )
h=0 (n) (n)
Temperature (C)

h=104
heat flows (due to use of Ti,α instead of Ti , which
60
differ near an interface). Because their heat flow must
h=3X104
Analytical
40
be based on the material temperature field (to be cor-
h=105 rect for lower h), these nodes cannot recover the perfect
h=∞ conduction limit at high h.
20
Correcting the analysis to revert to perfect con-
0 duction limit at high h is relatively unimportant be-
-50 -40 -30 -20 -10 0 10 20 30 40 50
Position (mm) cause 1) the discrepancy from perfect conduction is typ-
Fig. 2 Temperature profile in a bar with two materials sep-
ically small and 2) when h is sufficiently high, the sim-
arated by 0.5 mm at x = 0 mm. All material points started ulation can simply be run using standard conduction
at T = 0◦ C and results are after 500 ms with T set to 0◦ C at that ignores interfaces. Nevertheless, a simple correc-
x = −50 mm and 100◦ C at x = +50 mm . The solid lines are tion seems to work. Realizing that heat flow at single-
MPM results for various values of interfacial convection coef-
ficient h (in W/(m2 -K)). The dashed red line is the analytical
material nodes ±1 cell from an interface may have slight-
solution for perfect conduction. ly reduced heat flow, one option is to compensate for
(n)
that by allowing the ∆qi,α on interface nodes to slightly
∗(n)
exceed the conduction limit ∆qi,α . Figure 4 magnifies
Figure 1 shows results after 500 ms with the in- the discontinuity region for very high h as a function
terface pushed into contact. Because contact was mod- of the excess heat flow allowed over conduction. As the
eled as conduction, the results should reduce to perfect allowed excess increased to 1% the simulations got very
conduction. The MPM solution (solid black line) is close to the analytical solution. Increasing beyond 1%
nearly identical to the analytical solution (dotted red did not improve the agreement and by 4%, the excess
line). Figure 2 shows the analogous results when the heat flow caused a temperature instability.
interface was separated. This case developed a temper- Although allowing 1% excess allowed high-h simu-
ature discontinuity at the interface. When h = 0 no lations to approach perfect conduction, what are the
heat is transferred across the interface and all temper- consequences on that excess heat flow when h is lower?
ature rise was to the right of the interface (x > 0). The Figure 5 shows the temperature discontinuity at the in-
left side material correctly stayed exactly at zero. As h terface in the separated 1D bar example as a function
increased, the discontinuity got smaller and the results of h for either 1% excess heat flow (solid line) or 0%
approached the analytical solution for perfect conduc- excess heat flow (dashed line). For low h, the added
tion. excess had no effect (because convection flow is always
Notice that the h = ∞ (which used very large h) less than the conduction limit), but at high h, the 1%
does not exactly match the analytical solution, which excess allowed the discontinuity to approach the cor-
8 J. A. Nairn

A 60

50

40

ΔT (C)
30

20
-3 -2 -1 0 +1 +2 +3
0% Excess
B 10
1% Excess
0
102 103 104 105 106 107 108
h (W/(m2 K))

Fig. 5 Temperature discontinuity at the material interface

for a bar using conditions identical to Fig. 2 as a function of
h and for 0% (dashed line) or 1% (solid line) excess heat flow
allowed at interface nodes.
C

All previous examples had the interface close to a

grid line in the background grid, but it is important
to verify heat flow calculations are independent of in-
terface location within that grid. Figure 3C shows an
interface through midpoints of background cells. The
Fig. 3 Material points (solid circles) on top of a background circled nodes show that twice as many interface nodes
grid for cases when the material interface between “black” are present for such an interface compared to an in-
and “gray” materials is A. in contact. B. separated. C. sep- terface near a grid line (c.f., Fig. 3B). If contact heat
arated but the interface is near the midpoint of background flow calculations do not account for total number of in-
cells. The red empty circles indicate interface nodes in the
MPM calculations that can adjust heat flow to account for terface nodes, these two interface locations would give
the interface. different results. This issue is identical to contact force
calculations in MPM [11] and solved by proper choice
34
of contact area ai in Eq. (41), which should be read as
0% an effective contact area. When the interface is near a
32 0.3% (n) (n) (n) (n)
grid line (Fig. 3B), vi,α = vi,β ≈ vi /2, where vi
1% is total volume extrapolated to node i. The effective
30
Temperature (C)

(n)
contact area reduces to ai = vi /t⊥ , which is equal to
28
Analytical the physical contact area. When the interface is through
(n) (n)
26 midpoints of background cells (Fig. 3C), vi,α ≈ 7vi /8
(n) (n)
24 and vi,β ≈ vi /8, where material α is the one that
surrounds the node [11]. The resulting effective con-
22 Contact gap (n)
tact area, ai = vi /(2t⊥ ), now correctly assigns half
20
-10 -8 -6 -4 -2 0 2 4 6 8 10 the contact heat flow to each interface node in Fig. 3C
Position (mm)
where the interface is affecting twice as many nodes.
Fig. 4 Temperature profile in a bar using conditions identical Figure 6 compares results for an interface near a grid
to Fig. 2 for very high h (effectively infinite) as a function
of the excess heat flow allowed at interface nodes. The plot
line to one through midpoints of background cells for
magnifies the region near the interface. The dashed red line is h = 3×104 W/(m2 -K). The results are nearly identical.
the analytical solution for perfect conduction at the interface. The dashed line, however, shows results for an interface
through midpoints but using the physical contact area
(n)
(≈ vi /t⊥ ) instead of the effective contact area ai . Be-
rect limit of zero. In brief, it appears reasonable to al- cause the uncorrected area is about twice ai , the added
ways allow 1% excess heat flow. Making this excess a heat flow is two times too high leading to a temperature
simulation parameter would allow it to be adjusted for distribution corresponding to doubling of h. For realis-
different problems if needed. The proper excess would tic simulations, especially when interfaces move during
be determined by comparing high-h simulation results calculations, grid independence requires use of ai . The
to a simulation that ignores interfaces. overlapping, h = ∞ curves show that the results are
Contact Heat Flow 9

100 70

60 h=0
80
Near Grid Line 50
Temperature (C)

60 40

ΔT (C)
h = 3X104
h = 3X104 30
40 Unscaled Ai
20
Analytical Near Cell Midpoint
10
20
h = ∞ (both)
0
0 -10
-50 -40 -30 -20 -10 0 10 20 30 40 50 -0.5 0.0 0.5 1.0 1.5 2.0 2.5
Position (mm) Interface Separation (cells)

Fig. 6 Temperature profile in a bar using conditions identical Fig. 7 Temperature discontinuity at a material interface
to Fig. 2 for h = 3 × 104 W/(m2 -K) and h = ∞ for simula- for a bar using conditions identical to Fig. 2 with h =
tions with the interface near a grid line or through midpoints 3 × 104 W/(m2 ·K) as a function of separation between the two
of background cells. The dashed blue line is for an inter- materials. The dotted red line is temperature discontinuity
face through midpoints of cells but using physical instead of when h = 0.
effective contact area. The dotted red line is the analytical
solution for perfect conduction at the interface.

While it is physically reasonable for effective convection

independent of interface location and are identical to coefficient h to decrease as separation increases, the sep-
the analytical solution. Although less correction was aration dependence apparent in Fig. 7 is controlled by
needed when the interface was near the midpoints of background grid cell size and not physics of heat trans-
cells (because of the higher number of interface nodes), port. This property could be changed by allowing the
each calculation agreed well with the perfect conduction interfacial heat flow (see Eq. (42)) to be any function
(n) (n)
limit by allowing the recommended 1% excess heat flow. of Ti , Ti,α , and δ (where δ is interfacial opening dis-
To validate heat flow analysis across cracks, simula- placement available in both contact [11] and explicit
tions in Fig. 2 were repeated with the material interface crack [9, 10] calculations). For example, decreased heat
being replaced by a crack within a single material or flow as a function of separation could be modeled by
by a crack part way along a material interface in mul- allowing h to decrease as separation increases. But no
timaterial mode. The first tested heat flow at cracks function can compensate for loss of thermal contact af-
instead of material interfaces while the second tested ter “numerical separation.” Such wide separations be-
problems that combine cracks and multimaterial mode. come a problem of modeling heat flow through voids
These new results were identical to the results in Fig. 2 rather than across interfaces. Void modeling will require
(therefore no new plot is shown). new methods. One approach might be to fill voids with
Figure 7 shows temperature discontinuity at an in- a gas having different thermal conductivity properties
terface for the simulation in Fig. 2 with h = 3 × 104 and then modeling contact between gas and solid at
W/(m2 ·K), but now as a function of interfacial sepa- material interface by the methods presented here.
ration. For interfaces in contact (separation < 0), the Two final examples considered real-world examples
modeling assumed conduction and ∆T was always zero. that illustrate problems where interfaces or cracks might
For separated interfaces, ∆T increased as a function affect heat flow. Vacuum arc remelting (VAR) is a tech-
of separation and approached the h = 0 limit (dotted nologically-important metals process [8]. In brief, a cy-
red line) for separation of more than two cells in the lindrical ingot is remelted by electrical currents into a
background grid. A separation of more than two cells cylindrical crucible. If the process can be sufficiently
corresponds to materials that are sufficiently far apart controlled, the re-solidified ingot can have fewer im-
that the simulation no longer has interface nodes. When purities and better properties. Because VAR is an im-
there are no interface nodes, no calculations are done portant and expensive process, any modeling methods
to adjust heat flow and thus heat transfer across the that can help control it would be beneficial. Full VAR
gap drops to zero. In other words, a separation of more modeling requires many features such as phase transi-
than two cells results in numerical separation and ther- tions, recrystallization kinetics, Lorentz forces induced
mal isolation of the materials. from electric currents, and more. A possible MPM ap-
Because of numerical separation, the modeling de- proach to this problem will be in a future publication.
scribed here works best for interfaces and cracks with This paper considers only a VAR-related example for
small separations (less than the background cell size). cooling a cylindrical ingot within a crucible under con-
10 J. A. Nairn

ditions where thermal contraction of the ingot causes 100

it to separate from the crucible walls. Because VAR is 90

conducted under vacuum, such separation is likely to 80

<Temperature> (C)
affect cooling efficiency. 70
60
Consider a cylindrical ingot (radius 40 mm; length 50
h=0
100 mm) fully surrounded by a crucible. Mechanical 40 h=10 4

and thermal properties of ingot (i) and crucible (c) 30 h=105

were set to Ei = Ec = 1 MPa, νi = νc = 0.33, αi = 20 h=106
40 × 10−6 C−1 , αc = 0 C−1 , ki = 10000 W/(m·K), 10 h=∞
kc = 50000 W/(m·K), Ci = Cc = 1000 J/(kg·K), and 0
0 500 1000 1500 2000
ρi = ρc = 1 g/cm3 . Most importantly, the ingot’s ther- Time (ms)
mal expansion coefficient was higher than the crucible’s Fig. 8 Average temperature in an ingot initially heated from
such that the ingot will pull away from the crucible as it 50 to 100◦ C followed by immersing the crucible in a 0◦ C
cools. The axisymmetric simulation used 2.5 × 2.5 mm reservoir as a function of cooling time for various values of h
cells, started with all material points at 50◦ C, and grad- (in W/(m2 -K)).
ually heated them to 100◦ C. The heating was done to
100
induce interfacial pressure such that all surfaces start in
full contact. Gradual heating (relative to the material’s 90
wave speed) was done to minimize dynamic stress ef-
<Temperature> (C)

80
fects. After reaching 100◦ C, two different cooling meth-
70 h=0
ods were used. First, the entire crucible was immersed
h=105
in a thermal bath at 0◦ C and crucible outer surfaces 60
h=107
were modeled by heat flux boundary conditions on sur- 50
h=∞
face particles of q (n) = hsurf (T − Tres ) (see Eq. (16)) h=10 6
40
where hsurf = 105 W/(m2 ·K), T is particle tempera-
ture, and Tres = 0◦ C is the reservoir temperature. The 30
0 500 1000 1500 2000
average temperatures in the ingot as a function of time Time (ms)
after reaching 100◦ C for various values of h across sep- Fig. 9 Average temperature in an ingot initially heated from
arated ingot/crucible interfaces are plotted in Fig. 8. 50 to 100◦ C followed by the bottom half of the crucible being
At first, all cooling was close to the perfect conduction immersed in a 0◦ C reservoir while the top half remained in
limit (h = ∞ curve), but lower h gave slightly slower a 100◦ C reservoir as a function of cooling time for various
values of h (in W/(m2 -K)).
cooling because thermal gradients induced some loss of
contact at the ingot/crucible interface. After the aver-
age temperature dropped below the initial temperature Transitions between slow and fast cooling indicate dy-
of 50◦ C, the entire ingot separated from the crucible namic motion of the ingot within the crucible. The role
and the cooling rate significantly slowed. For h = 0, of dynamic contact between ingot and crucible on heat
heat transfer stopped and therefore cooling stopped. flow would be difficult to model by anything other than
As h increased, cooling continued but cooling rates de- numerical methods such as those proposed here.
pended strongly on h. The second real-world example modeled cracks. One
Figure 9 shows corresponding results where the bot- method for detecting internal cracks in solid, opaque
tom half of the crucible was immersed in a 0◦ C reser- objects is by thermal imaging [1]. In brief, an object
voir while the top half remained in a 100◦ C reservoir. is heated by various methods and the surface temper-
The perfect condition limit (h = ∞ curve) gave the ature is imaged using an infrared camera. If cracks af-
expected result that average temperature approached fect heat flow, the number and location of cracks will
50◦ C with a steady state temperature gradient from alter the temperature distribution. The heat-flow meth-
cold side to hot side. The results with interfacial heat ods described here could help to interpret such exper-
transfer, however, gave non-intuitive results. The cool- iments. Figure 10 shows two examples of the effect of
ing curves show regions of reduced cooling rates and cracks on temperature distributions during heat flow.
regions of higher cooling rates. These regions are deter- The 100 × 50 mm isotropic material (E = 1 MPa, ν =
mined by which parts of the interface were in contact. 0.33, α = 0, k = 10000 W/(m·K), C = 1000 J/(kg·K),
When the hotter walls on the top half lost contact, the and ρ = 1 g/cm3 , thickness = 1 mm) started at 0◦ C.
cooling rate was higher, but when the colder walls on The bottom edge was heated using constant heat-flux
the bottom half lost contact, the cooling rate slowed. boundary condition of q (n) = 107 W/m2 , the top edge
Contact Heat Flow 11

56

54 No Cracks
One Crack

Surface Temperature (C)

52 Four Cracks

50

48

46

A. 44

42
-10 0 10 20 30 40 50 60 70 80 90 100 110
Position (mm)

Fig. 11 Surface temperature from the thermal fields in

Fig. 10 for one crack (dashed black line), four cracks (solid
blue line), or no crack (dotted red line).

4 Conclusions
B.
This paper describes a new MPM feature for modeling
Fig. 10 Temperature distribution in an object with A. one
crack or B. four cracks after 500 ms of heat flux applied to the heat flow across material interfaces, cracks planes, or
bottom surface. The image lightness indicates temperature both. The method works best for interfaces with small
from 40◦ C (dark) to 110◦ C (light) and 7◦ C per contour. separations (less than one cell size in the background
grid). The method has potential for analysis of real-
world problems that couple mechanical response with
was cooled by convection using q (n) = hsurf (T − Tres )
heat flow across interfaces or cracks where temperature
where hsurf = 105 W/(m2 ·K) and Tres = 0◦ C, and
changes may cause interfaces or crack planes to dynam-
the sides had zero heat flux. To maximize crack ef-
ically lose or regain contact.
fects, convective heat flow across cracks surfaces used
h = 0. The temperature distribution was evaluated at
the end of a 500-ms heating period. The MPM model
used 2.5 × 2.5 mm cells. Compliance with ethical standards
The first simulation had a single crack in the middle
of the object, which, as shown in Fig. 10A, influenced Funding: This work was made possible by the endow-
the temperature field. Thermal imaging can only ob- ment for the Richardson Chair in Wood Science and
serve surface temperature. Figure 11 shows a drop on Forest Products. The author also thanks Rigel Wood-
the top surface temperature above the crack compared side, Paul King, and Kevin Gartner for helpful discus-
to temperature distribution that would occur with no sions.
cracks. The length of the crack could be estimated from Conflict of interest: The corresponding author states
this drop, but the details would depend on distance of that he has no conflicts of interest.
the crack to the surface and on the convection coeffi-
cient for heat flow across the crack surface.
The second simulation used four arbitrarily sized References
and placed cracks. The temperature field was signifi-
cantly altered (see Fig. 10B) and the surface tempera- 1. Almond, D.P., Peng, W.: Thermal imaging of com-
ture showed highly convoluted effects of cracks on sur- posites. Journal of Microscopy 201(2), 163–170
face temperature distribution. For example, the tem- (2001)
perature above the last crack on the right was higher 2. Bardenhagen, S.G., Guilkey, J.E., Roessig, K.M.,
than average. That crack was too far below the sur- Brackbill, J.U., Witzel, W.M., Foster, J.C.: An im-
face to have continued effect on surface temperature. It proved contact algorithm for the material point
is unlikely the surface temperature alone can uniquely method and application to stress propagation in
locate internal cracks, but modeling coupled with ex- granular material. Computer Modeling in Engi-
periments could lead to derivation of a distribution of neering & Sciences 2, 509–522 (2001)
equivalent cracks that are consistent with observed sur- 3. Bardenhagen, S.G., Kober, E.M.: The generalized
face temperature variations. interpolation material point method. Computer
12 J. A. Nairn

Modeling in Engineering & Sciences 5, 477–496

(2004)
4. Brackbill, J., Kothe, D., Ruppel, H.: FLIP: A low-
dissipation, particle-in-cell method for fluid flow.
Computer Physics Communications 48(1), 25 – 38
(1988)
5. Carlson, D.E.: Linear thermoelasticity. In:
C. Truesdell (ed.) Mechanics of Solids, vol. II, pp.
297–345. Springer-Verlag, New York (1984)
6. Guenther, R.B., Lee, J.W.: Partial Differential
Equations of Mathematical Physics and Integral
Equations. Dover Books, New York (1996)
7. Lemiale, V., Hurmane, A., Nairn, J.A.: Material
point method simulation of equal channel angular
pressing involving large plastic strain and contact
through sharp corners. Computer Modeling in Eng.
& Sci. 70(1), 41–66 (2010)
8. Maurer, G.E.: Primary and secondary melt process-
ing — superalloys. In: J.K. Tien, T. Caulfield (eds.)
Superalloys, Supercomposites and Superceramics,
pp. 49–97. Academic Press, New York (1989)
9. Nairn, J.A.: Material point method calculations
with explicit cracks. Computer Modeling in En-
gineering & Sciences 4, 649–664 (2003)
10. Nairn, J.A.: Numerical implementation of imper-
fect interfaces. Computational Materials Science
40, 525–536 (2007)
11. Nairn, J.A.: Modeling of imperfect interfaces in the
material point method using multimaterial meth-
ods. Computer Modeling in Eng. & Sci. 92(3),
271–299 (2013)
12. Nairn, J.A.: Material point method
(NairnMPM) and finite element anal-
ysis (NairnFEA) open-source software.
[Link] (2017)
13. Nairn, J.A., Aimene, Y.: Modeling interacting and
propagating cracks in the material point method:
Application to simulation of hydraulic fractures in-
teracting with natural fractures (2017). In prepa-
ration
14. Nairn, J.A., Bardenhagen, S.G., Smith, G.S.: Gen-
eralized contact and improved frictional heating in
the material point method. Computer Particle Me-
chanics 5(2), 285–296 (2018). DOI 10.1007/s40571-
017-0168-1
15. Nairn, J.A., Guilkey, J.E.: Axisymmetric form
of the generalized interpolation material point
method. Int. J. for Numerical Methods in Engi-
neering 101, 127–147 (2015)
16. Sulsky, D., Chen, Z., Schreyer, H.L.: A particle
method for history-dependent materials. Comput.
Methods Appl. Mech. Engrg. 118, 179–186 (1994)