Chemical Engineering Science 192 (2018) 448–466

Contents lists available at ScienceDirect

Chemical Engineering Science

journal homepage: www.elsevier.com/locate/ces

Development of novel heat conduction interaction model for solid body

thermal contact in CFD based particle flow simulations
L.A. Florio
US Army ARDEC, Technology Branch, Small Caliber Armaments Division, U.S. Army Armament Graduate School at Picatinny Arsenal, Picatinny Arsenal, NJ 07806, United States

h i g h l i g h t s g r a p h i c a l a b s t r a c t

Conductive contact modeled when

gap between computational elements
must be maintained.

Volumetric energy sources simulate
thermal exchange at flagged
contacting elements.

Method proves to be robust for range
of contact resistances and thermal
properties.

Excellent agreement with standard
finite volume solution with
contacting elements.

Facilitates study of thermal contact in
dynamic, resolved particle flow/fixed
mesh topology.

a r t i c l e i n f o a b s t r a c t

Article history: Traditional computational fluid dynamics (CFD) simulations with resolved objects progressing through
Received 12 March 2018 the fluid using a moving and deforming mesh require a fixed mesh topology. Though such simulations
Received in revised form 16 June 2018 may more readily capture conditions where large gradients in the flow field exist, particularly high speed,
Accepted 18 July 2018
compressible flow phenomena, the ability to directly model the conduction portion of the contact
Available online 25 July 2018
between particles or between particles and other solid bodies immersed in the fluid is thus restricted.
A new method was devised to allow for the maintenance of the mesh topology while simulating the con-
Keywords:
tact heat conduction between a discretized particle pair and between a discretized particle and a wall for
Particle technology
Heat conduction
the intended purpose of implementation of the technique within a computational fluid dynamics (CFD)
Thermal contact based coupled particle–fluid flow model. The method is tested for a range of contact resistances
(0–1  103 m2 K/W) and material thermal conductivity (1.60–162 W/m K) and specific heat
(50.28–5028 J/kg-K) combinations with one object initially at 700 K and the other at 300 K. The temper-
ature distributions across the mid-particle or mid-particle and wall line produced by the new method at a
given time are compared to the results for a standard finite volume simulation of contacting bodies. The
maximum difference in the local temperatures obtained is under 1%. Hence, the technique is a robust
means of simulating the resolved heat conduction contact thermal energy exchange without connected
computational elements, retaining the needed fluid gap space around the solid bodies. Future work will
build upon this technique to include the convective and radiative modes of heat transfer before simulat-
ing moving and interacting particles in a fluid flow.
Published by Elsevier Ltd.

E-mail address: laurie.a.florio.civ@mail.mil

https://doi.org/10.1016/j.ces.2018.07.042
0009-2509/Published by Elsevier Ltd.
 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466 449

Nomenclature

Vectors⁄ Description [Eq.] Lm contact related length Fig. 2 [Eq. (5)]

A vector in Fig. 3 representing wall [Eq. (7)] NC number of boundary cells [Eqs. (15) and (16)]
B vector in Fig. 3 from x1 to A [Eq. (7)] nc1, nc2 parameters [Eqs. (19a) and (19b)]
d vector from particle 1 to particle 2 [Eq. (1a)] ndot distance (m) [Eq. (13)]
n unit normal vector (Fig. 1) [Eq. (1b)] R1, R2 particle radius Fig. 1 (m) [Eq. (2)]
np unit normal vector (Fig. 2) [Eq. (4b)] Q_ volumetric heat/energy source (W/m3 K) [Eq. (20)]
nw unit normal to wall (Figs. 3 and 4) [Eq. (4b)] S particle clearance (m) [Eq. (2)]
p vector from particle centroid to outer surface Fig. 2 [Eq. T temperature (K) [Eq. (3)]
(4a)] t time (s)
p1 coordinate of outer surface of particle Fig. 2 [Eq. (4a)] Tc contact temperature (K) [Eqs. (19a) and (19b)]
r particle center of gravity to contact position vector Fig. 4 tdot distance (m) [Eq. (12)]
(m) [Eq. (8)] Vol computational cell volume (m3) [Eq. (20)]
tw unit tangential vector along wall Fig. 3 [Eq. (8)] W wall dimension Fig. 8 (m)
x1, x2 center of gravity positions of particle 1 and 2 (m) [Eq. w_p, w_t distances Fig. 4 (m) [Eqs. (11a) and (11b)]
(1a)]

Greek
Scalars q particle density (kg/m3) [Eq. (21)]
2
Ai computational cell area (m ) [Eq. (15)]
CL contact area for particles [Eq. (3)]
CR contact resistance Fig. 6 (m2 K/W) [Eq. (18)] Notation
cp specific heat (J/kg K) n associated with normal
CZ parameter for determining collision [Eq. (2)] 1, 2 associated with particle 1 and particle 2
d cell centroid normal distance (m) [Eq. (17)]
k thermal conductivity (W/m-K) [Eq. (3)] ⁄
In the manuscript, bold indicates a vector.
L wall dimension Fig. 8 (m)
la dimension (m) [Eq. (6)]

1. Introduction A sample of the diverse set of complex processes and phenom-

ena that stem from thermally dependent phenomena in particle
Temperature levels and temperature gradients trigger or drive laden flows includes chemically reacting flows, ceramics process-
many significant phenomena in particle laden flows, from changes ing, powder metallurgy, packed bed reactors, blast furnace opera-
in material properties to chemical reactions, phase change, adhe- tions, propellant burn, flow path ablation, turbine blade wear,
sion, or erosion. Resolved conduction, convection, and radiation coal pyrolysis, additive manufacturing, fouling of flow paths,
heat transfer interactions between fluid flow, particulates, and flow magneto-hydrodynamics applications, and novel solar collectors
path surfaces is often necessary to capture these local condition (Das et al., 2017; Haddad et al., 2016; Kosinski et al., 2013, 2014;
dependent phenomena. While alternative techniques such as Lattanzi and Hrenya, 2016; Liu and Li, 2016; Liu et al., 2013;
discrete-element/porosity models coupled with computational Morris et al., 2016a, 2016b; Oschmann et al., 2016; Vargas and
fluid dynamics (CFD), immersed boundary condition, or Lattice- McCarthy, 2002; Yang et al., 2015). With the wide ranging applica-
Boltzmann methods have been implemented to obtain the flow tions, advancement in the body of knowledge of underlying physics
and thermal conditions, along with any mechanical or, in some of these processes and phenomena is needed. Though physical
cases, thermal interactions among the system components these experiments can provide bulk data (Yang et al., 2015), the level of
methods may often not be adequate to capture the phenomena detail required to improve the understanding of these processes
of interest. Particularly in compressible flows or other regimes often cannot be reached through physical experiments. Numerical
with high velocity, temperature, or pressure gradients, computa- simulations, in contrast, can be used to gather such data. Individual
tionally resolved particles, flow paths, and fluid flow with a moving particle translational and rotational displacements, velocities, and
and deforming mesh may be needed. These resolved CFD-moving accelerations can be traced. Available local time and space depen-
particle models generally require a fixed mesh topology so that a dent flow and temperature field distributions can be used to esti-
fluid layer is maintained around each solid body, keeping the ‘‘con- mate fluid induced forces and moments acting on solid bodies.
tacting” bodies from being explicitly connected in the computa- The portions of objects in contact, the mechanical interaction forces
tional domain. Means to simulate the resulting mechanical and moments, and the convective, conductive, and radiative ther-
interaction being solid bodies are more advanced than means to mal interactions can then be modeled. Hence, the development
simulate the multimode thermal interactions between solid bodies and application of locally resolved particle-flow based numerical
in the flow with the most common so called soft collision mechan- modeling techniques is critical to furthering the understanding of
ical interaction models able to consider elastic, elastic-plastic, and the particle flow thermal interaction phenomena and exploiting
even strain hardening material models. The modeling technique these thermal phenomena at work in particle laden flows
developed and described in this work is the foundation of an alter- (Oschmann et al., 2016, Yang et al., 2015; Haddad et al., 2016).
native multimode thermal interaction modeling method. This Commonly implemented multiphase flow methods apply sim-
work presents a technique to simulate the conductive thermal plifying assumptions so that the models may no longer consider
interaction between solid bodies without contact between the some of the important local flow and thermal effects truly needed
computational representations of the particles or other objects, to capture the thermal interaction phenomena. Eulerian methods
allowing for the simulation of locally resolved thermal conductive treat the ‘‘particles” as a continuum and as such cannot follow
interactions. individual particles. Lagrangian based models track individual
450 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466

particles, but a wide range of treatments are employed. Low parti- Oschmann et al. (2016) began work on a more detailed modeling
cle volume fraction methods neglect particle interactions, discrete technique in the study of heat conduction within particles in the
element mode (DEM)-CFD coupling may capture particle interac- presence of a surrounding fluid for use in a DEM model. Good com-
tions, but require correlations to estimate flow induced forces parisons are made to the results from a commercial finite volume
and thermal effects as the particles are not resolved in the flow code. With the more detailed thermal resolution, Oschmann et al.
(Morris et al., 2016a; Lattanzi and Hrenya, 2016; Liu et al., 2013). note significant differences in the heat fluxes, effective thermal
Correlations have been developed mainly for spherical particles, conductivity of the packed particles, and the time to reach steady
requiring special studies for other shapes such as cylindrical or state, indicating the inclusion of these localized, resolved thermal
elliptical particles (Singhal et al., 2017; Ga et al., 2016). As the par- phenomena is important in advancing the understanding of the
ticle volume fractions increase, standard correlations for flow conditions occurring in particle laden flows. Hence, more general
induced forces on the solids and convective flow-solid thermal techniques are needed that do not rely heavily on correlations,
interactions may not be appropriate as the flow around one parti- are not tied to a particular flow regime, particle speed, particle size,
cle is increasingly affected by the flow over neighboring particles. or particle concentration, and allow for the determination and
Additionally, such a formulation inherently assumes the particle usage of the local temperature and flow field and the local conduc-
size is smaller than the computational cell size, which may not tion, convection, and radiation interactions for moving and collid-
be the case for many flow conditions. Though thermal conditions ing particles carried with a fluid. Such methods can be broadly
in solid bodies may be included, immersed boundary methods applied to particle flows to further the understanding of the ther-
(Das et al., 2017; He and Tafti, 2017; Ardekani et al., 2018) and Lat- mal interaction characteristics that occur during particle transport.
tice Boltzmann methods (McCullough et al., 2016) may not be able The aim of this paper is to report on the development of a novel
to capture the flow and temperature field gradients that develop in modeling method that extends the direct flow-particle modeling
higher speed flows or the compressibility effects that are critical to method with mechanical interaction modeling method developed
the localized temperature dependent processes. by the author to include the thermal effects within a fully compu-
Florio (2014, 2015a, 2015b, 2017) has developed a direct flow- tational fluid dynamics based environment. This publication
particle coupling method specifically focused on high speed com- focuses on the development and implementation of the conduction
pressible flow applications with the particles actually incorporated contact thermal simulation technique for replicating the effects of
in the computational domain with a moving and deforming mesh, thermal conduction contact between a particle and another parti-
allowing for particle motion through the fluid. A range of mechan- cle and a particle and a wall for a range of particle and wall thermal
ical interaction models and methods for simulating various physi- material properties and contact resistances. This work sets the
cal phenomena has been devised based on the basic modeling basis for a modeling method that includes the fluid-particle con-
concept (Florio 2014, 2015a, 2015b, 2017). However, detailed ther- vection, conduction, and radiation effects with moving and inter-
mal interactions have not been included. acting particles.
Apart from the immersed boundary/Lattice Boltzmann meth-
ods, most work in particle flow, and particle flow with thermal
2. Modeling method
effects, has been in the area of DEM-CFD coupled type flows. While
contact effects between solid bodies may be considered, the parti-
The significant components of the proposed method to model
cles are not resolved in the flows and the particle volume fraction
discretized, localized heat conduction contact between two solid
is passed to the CFD code for an essentially Eulerian treatment of
bodies without requiring computationally connected geometries
the particulates. Work continues to build convective/conductive
are described in detail. The major steps include:
correlations for the thermal interactions such as in the studies by
Singhal et al. (2017) and Ga et al. (2016). Vargas and McCarthy
a. The determination of the mechanical contact related param-
(2002) consider a lumped capacitance formulation to estimate
eters and the specific elements on each of the particles/
the heat flow between colliding particles using a correlation based
objects that are in contact
total contact conductance, noting the model ‘‘breaks down” when
b. The calculation of the heat exchange between the particles/
the solid conductivity approaches the fluid conductivity. Haddad
objects based on the local temperature conditions, the ther-
et al. (2016) developed a contact conductance based modeling
mal properties of the particles/objects, the contact resis-
method to maintain a bulk particle temperature and Kosinski
tance, and the contact geometry
et al. (2013, 2014) utilize a semi-infinite solid approximation to
c. The assignment of this local energy exchange via volumetric
estimate the thermal heat exchange between contacting particles.
sources to the proper elements on each of the contacting
Lattanzi and Hrenya (2016) developed a lumped particle tempera-
particles/objects
ture, one dimensional heat conduction contact model with differ-
d. The determination of the new temperature distribution
ent formulations for particle–fluid-particle or particle fluid-wall
through the particles/objects by implementation of the
contacts activated when a ‘‘static lens” is placed around the parti-
CFD solver.
cle, overlaps with that of another body. Morris et al. (2016b) intro-
duce a correction factor to account for the artificially low Young’s
modulus applied for most DEM simulations for numerical effi- 2.1. Mechanical contact
ciency. Yang et al. (2015) use a parameter based conduction inter-
action formulation, correlations for the fluid-particle convective The first step in modeling the thermal interaction between two
interaction, and a radiation exchange solely to the flow passage, objects is determining if the particles/objects are in contact and
not other particles, based on the surface average temperatures. then finding the specific contacting elements on the boundary sur-
As seen in this literature review, particle temperatures are often face of the objects in contact. Different procedures will be used for
treated under a bulk or lumped simplifying assumption set, with the interacting particles, here all considered circular, and an inter-
conduction and convection exchanges typically requiring correla- acting particle and a wall.
tion equations. Progress towards a more localized thermal interac-
tion treatment for particles in a fluid based on resolved flow and 2.1.1. Particle – particle contact
temperature fields using methods other than the immersed bound- Consider the two particles in Fig. 1 with particle 1 of radius R1
ary or similar overset type boundary methods is limited. and centroid position vector x1 and particle 2 of radius R2 and
 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466 451

Fig. 1. Identifying contacting cells or element on each contacting object pair.

centroid position vector of x2. The vector pointing from the cen- The length, Lm, is the minimum possible value of the compo-
troid of particle 1 to the centroid of particle 2, d, is given by: nent of the vector p in the direction n, for which a point p1 is
within the contact zone. The value for Lm can be calculated from:
d ¼ x2  x1 ð1aÞ
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
where the unit normal vector, n, is:  2
1
Lm ¼ R21  CL ð5Þ
d 2
n¼ ð1bÞ
jdj
The contact status of the boundary cells can then be determined
The parameter, CZ, defined in Eq. (2), is used to track the colli- by finding the component of the np vector in the n direction, which
sion status of the particle pair. When CZ is greater than zero, the will be called la. The point p1 will be within the contact zone if the
particle pair is defined to be colliding or in contact. S is a clearance statement in Eq. (6) is true.
space or safety zone between the particles that is maintained
throughout a simulation that may include moving particles as Lm
6 la 6 1 with la ¼ np  n ð6Þ
the mesh topology must be maintained within the CFD model to R1
prevent mesh overlap. A similar procedure is used for particle 2, but with the normal
CZ ¼ R1 þ R2 þ S  jdj ð2Þ vector taken in the negative n direction.

Bold in the equation notation indicates a vector. If the two par-

ticles are defined to be contacting, additional parameters need to 2.1.2. Particle – wall contact
be determined so that the contacting portion of each of the two Next, the procedure to check for a collision and determine con-
particles can be identified. The contact length, CL, as seen in tact locations for an interaction between a particle and a wall sur-
Fig. 2, is estimated through a Hertz contact type assumption. The face are described. Because of the different shape of the surfaces
expression for the contact length is: involved and the differences in the characteristics of the normal
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vector and normal distance formulations, some modifications to
CL ¼ 2 2R1 CZ  CZ 2 ð3Þ the procedure used for particle-particle interactions are needed
when modeling the particle-wall contact. Consider a portion of
For each of the two particles in contact, the centroids of each of
the wall as a segment reaching from Point A to Point E and the
the outer boundary face cells are checked to determine if the cen-
potential interaction with particle 1 as in Fig. 3. All segment direc-
troid is within the contact zone, as indicated in Fig. 2. A centroid on
tions around the solid body are taken proceeding in the counter-
the outer surface element of particle 1 is located at p1 with the vec-
clockwise direction. The length of the segment of the wall is LW.
tor from the centroid of particle 1, x1, to the point p1 given by p (Eq.
The vector B reaches from Point A to the centroid of particle 1 at x1.
(4a)). The unit vector associated with the vector p is np (Eq. (4b)).
p ¼ p 1  x1 ð4aÞ B ¼ x1  A ð7Þ
The vector tw is a unit vector that runs parallel to the surface of
p
np ¼ ð4bÞ the wall as shown in Fig. 3 with nw as the outward directed normal
jpj
unit vector from the solid surface. The outward directed unit

Fig. 2. Contact zone identification, with contacting bounding elements on each particle indicated by the lighter color. (For interpretation of the references to color in this
figure legend, the reader is referred to the web version of this article.)
452 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466

Fig. 3. Contact zone identification, particle-wall interactions, parameter definitions.

normal from the solid particle surface at the point p1 is np as If a collision is detected, the elements on the boundary of both
defined in Eq. (4). the particle and the wall that are involved in the collision must be
The components of the B vector in the wall normal and wall identified. First, for the particle, a vector p extends from the cen-
tangential directions are calculated first to check for conditions troid of the particle, x1, to a given centroid of an element on the
under which the particle and wall might be colliding. The tangen- outer surface of the particle, p1 (Fig. 4a). Then, following the proce-
tial component, dt, of the vector B in the tw direction and the dure and conditions for the particle-particle contact, if the condi-
normal component, dn, of the vector B in the nw direction are tions noted in Eqs. (4a)–(6) are met, the element on the particle
found in Eq. (8). boundary is involved in the contact.
For the wall, a vector, w_p is directed from the wall boundary
dt ¼ tw  B dn ¼ nw  B ð8Þ surface element centroid, p1, to the centroid of the particle, x1. A
If the dn value, the normal component of vector B, is greater vector w_t is drawn from the start of the wall segment at A to
than zero, the particle is on the ‘‘correct” side of the wall and con- the boundary cell centroid, p1. These vectors are defined in Eqs.
tact is possible. If the dt value, the tangential component of the (11a) and (11b) (Fig. 4b).
vector B, also falls between 0 and LW, the total length of the wall w t ¼ p1  A ð11aÞ
segment, then the particle centroid lies somewhere in the zone
normal to but within the wall segment length. Hence, under these w p ¼ x1  p1 ð11bÞ
conditions a collision is possible. Next, the value for CZ is calcu-
lated to determine if the distance between the surfaces is within The dot product of tw and w_t is used to determine if the particle
the allowable ‘‘tolerance” for collision. The value of CZ is calculated centroid falls within the contact length as in Eq. (12).
as in Eq. (9) where S is half of the safety zone used for the two par-
tdot ¼ tw  w t ð12Þ
ticle collisions. The normal to the collision or contact used else-
where in the force and other calculations is in the negative nw The dot product of nw and w_p is also calculated to ensure the local
direction since the normal is always defined to be from the particle bounding element normal is facing in the proper direction for
1 to the wall (or particle 2). contact.

CZ ¼ R1 þ S  dn ð9Þ ndot ¼ nw  w p ð13Þ

So, if the two expressions, Eqs. (14a) and (14b) are met, then the
A contact length as defined in Eq. (10) is also determined.
wall bounding element is in the contact area.
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
CL ¼ 2 2R1 CZ  CZ 2 ð10Þ ndot > 0 ð14aÞ

However, simply checking if a particle is situated normal to a

   
CL CL
wall segment and within the appropriate distance does not capture dt  6 tdot 6 dt þ ð14bÞ
2 2
the collisions that may be occurring near the corners of the seg-
ments. Thus, in addition to the collision check above, the distance Recall, the ‘‘corner” contact is tested and treated independently.
between the particle centroid and the corner nodes of each the
wall segment is calculated and this distance is used to replace 2.2. Thermal interaction model
the dn value for the collision check in Eq. (9). For the ‘‘corner
check”, the normal to the ‘‘collision” is approximated as the unit With the collision status and bounding surface elements
vector in the direction from the particle centroid to the corner node involved in the particle-particle or particle–wall interaction identi-
on the wall surface. fied, the thermal interaction effects can be considered. The contact
 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466 453

Fig. 4. Contact zone identification, particle-wall interaction, particle contact.

area and the area weighted average temperature at the contact for 2.2.2. Calculation of the area weighted average temperature
each of the two contacting objects are determined from the CFD The area weighted average temperature on the contacting por-
model based computational geometry and temperature. This tem- tions of each of the contacting objects is determined utilizing the
perature and surface area information, together with the local bounding cell centroid temperatures as in Eq. (16). The cell tem-
material properties of the solid and the local bounding element peratures for object 1 are Ti and for object 2 are Tj. Then, the net
normal and bounding element centroid locations, is used to esti- area weighted average temperatures are called T1 and T2.
mate the local instantaneous heat transfer rate between the con- PNC 1 1 PNC 2 1
tacting bodies. The heat transfer rate is reformulated in terms of Ai T i j¼0 Aj T j
T1 ¼ i¼0
PNC 1 1 ; T 2 ¼ PNC 1 ð16Þ
a volumetric heat generation rate which is assigned as a volumetric i¼0 Ai j¼0
2
Aj
energy source to the appropriate interacting bounding elements on
each interacting body. The basic procedure used is outlined in this
section. 2.2.3. Calculation of the area weighted average normal distance
The area weighted average normal distance between the con-
2.2.1. Calculation of the local element contacting area tacting face and the cell centroid of the bounding element, d1
The total surface area of the contacting computational elements and d2, for particle 1 and particle/object 2 respectively, is also cal-
on each of the contacting bodies must be found, where the contact- culated in a similar manner to the temperature as shown in Eq.
ing area from object 1 is A1 and from object 2 is A2. As in Eq. (15), (17) (see Fig. 5).
the net area is found by summing the surface areas of the individ- PNC 1 1 PNC 2 1
Ai di j¼0 Aj dj
ual contacting elements where Ai and Aj represent the surface d1 ¼ i¼0
PNC 1 1 ; d2 ¼ PNC 1 ð17Þ
areas of the individual contacting elements from object 1 and i¼0 Ai j¼0
2
Aj
object 2 respectively and NC1 and NC2 are the total number of con-
tacting bounding elements from object 1 and object 2 respectively.
2.2.4. Heat flow through the contacting elements
X
NC 1 1 X
NC 2 1 At a given time step, the current computational based data are
A1 ¼ Ai ; A2 ¼ Aj ð15Þ used to estimate the heat flow through the contacting zones. A
i¼0 j¼0
contact resistance, CR, allows for the consideration of non-perfect

Fig. 5. Contacting element information defined.

454 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466

Fig. 7. System One: Two contacting particles.

Fig. 6. Parameters for heat flow at contact defined.

thermal contact in the simulations. The contacting zone can be thermal contact as in Fig. 7. The two particles are of the same size
modeled as an equivalent thermal circuit with the equivalent resis- and are fixed and so maintain a static contact length or area and a
tances as shown in Fig. 6. The left most zone represents particle 1, static contact condition. Only the heat conduction through the par-
the right zone particle 2, and the mid-zone represents the contact ticles is considered in this work as this is the first component of the
resistance. Tc1 and Tc2 are the contact temperatures of the each of thermal interaction modeling method. The significant system
the two ‘‘contacting” bodies. dimensions are provided in Table 1. The second system investi-
Neglecting any storage in the small volume of the contacting gated is one particle impacting a wall as shown in Fig. 8. Significant
elements and equating the heat flow across each of the material system dimensions are provided in Table 2.
sections, with k1 and k2 as the local area weighted average thermal
conductivity for particle 1 and particle/object 2 respectively, the 3.2. System conditions studied
heat flow is given by Eq. (18). As the element size decreases, the
assumptions used become more appropriate. A range of potential thermal conditions is tested utilizing the
ðT 1  Tc1 Þ ðTc1  Tc2 Þ ðTc2  T 2 Þ new modeling method described. An elevated temperature of
d1
¼ CR
¼ d2
ð18Þ 700 K is initially applied to particle 1 and the second object (the
k1 A1 A1 k2 A2 particle for the system one cases and the wall for the system two
Using the equivalent heat flow at the contact, the contacting tem- cases) is initially set to ambient thermal conditions of 300 K. The
peratures at object 1 and object 2, Tc1 and Tc2, respectively, can contact resistance and the thermal conductivity and specific heat
be determined. of the two solid bodies are varied to test and indicate the robust-
ness of the modeling method. The specific heat values included
ðT 2  T 1 Þnc1 d1 d2 in the investigation are 50.248, 502.48⁄, and 5024.8 J/kg-K and
Tc1 ¼ T 1 þ ; nc1 ¼ ; nc2 ¼ ð19aÞ
nc1 þ nc2 þ CR A1
k1 A1 k 2 A2 the thermal conductivity values are 1.6277, 16.277, and 162.77⁄
W/m-K (⁄ indicates the baseline value). The contact resistance
ðT 1  T 2 Þnc2 parameters are 0.00, 1.00E03, 1.00E04, 1.00E05, and
Tc2 ¼ T 2 þ ð19bÞ
nc1 þ nc2 þ CR A1
Table 1
Then, based on these temperature estimations, the flow of System One: Two contacting particles: Significant geometric parameters.
energy from one object to the other through a volumetric heat
source can be calculated. The expressions for the local volumetric Parameter Value (mm)

heat sources are provided in Eq. (20), where Voli is the particular R1 1.00
boundary computational cell volume for a cell on particle 1 and R2 1.00
d 1.98
Q_ i is the local cell volumetric energy source for an element on par- CL 0.282134
ticle 1. The items with the j subscript correspond to particle 2. The
formulation in Eq. (20) allows for the variation in the heat
exchange along the contacting elements, following the local varia-
tion in the temperatures around each of the bounding elements of
solids. Simply assigning a uniform volumetric source at all of the
boundary cells tends to result in the development of anomalies
in the temperature field since the local cell temperature vary. Jus-
tification of the use of the volumetric source at the boundary cells
is provided in Section 4.
ki ðTc1  T i Þ Ai kj ðTc2  T j Þ Aj
Q_ i ¼ ; Q_ j ¼ ð20Þ
di Voli dj Volj

With this formulation, the conductive contact between solids can

be simulated.

3. Problems considered

3.1. Two major systems

To develop and demonstrate the method, two main geometric

systems were studied each involving a planar, two dimensional
model to demonstrate the modeling method. The first is a two-
particle system with a portion of each particle in mechanical and Fig. 8. System Two: Particle – wall interaction.
 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466 455

Table 2 particle and the wall. However, in order to implement the method
System Two: Particle – wall interaction: Significant geometric parameters. described, significant customization of such a general code is
Parameter Value (mm) required. This customization is carried out through the use of user
R1 1.00 defined functions, or C code written by the author, that is coupled
L 10.00 to the commercial code. Through the code, the particle positions
W 3.00 are stored, contact conditions are checked, contacting elements
d 0.99 are flagged (Fig. 9), and the mechanical contact parameters are
CL 0.282134
determined based on the geometry and material properties. The
thermal contact related parameters are also found including the
1.00E06 m2 K/W. All cases involve an object (particle or wall) surface area weighted average temperatures, the surface to cen-
with a density of 8030 kg/m3. troid distances, the resulting contact temperatures and heat flux,
For system one, with two colliding particles, two main studies and finally, the volumetric heat source conditions. Through stan-
were conducted: dard finite volume techniques, the local and time varying thermal
conditions in both of the contacting, yet computationally discon-
a. Part 1 – The two particles are of the same material at the nected objects are then calculated, numerically solving Eq. (21)
baseline property values and the set of contact resistances which can account for variable material properties where Q_ is
listed above are applied. Temperatures are compared at the volumetric heat or energy source:
0.017 s using discrete data along the y = 0 line through both
particles as well as contour plots of the temperature field
(see Fig. 1 for coordinate system location at centroid of par-
ticle 1).
b. Part 2 – Various combinations of the specific heat and ther-
mal conductivity are examined for two contact resistance
values of 1.0  104 and 1.0  105 m2 K/W. The parameters
investigated are given in Table 3. Again, the temperature dis-
tributions at 0.017 s are compared including the tempera-
ture values along the line y = 0 and temperature contour
plots.

For system two, with the particle and the wall, one set of studies
is conducted. For a contact resistance of 1.0  105 m2 K/W, a lim-
ited set of material properties is examined as listed in Table 4.

4. Implementation of the method

A standard, commercial finite volume numerical tool, FLUENTÓ,

Fig. 9. Flagged contacting elements.
is used to solve for the temperature field in the particles or the

Table 3
Particle-particle interaction cases applied to system 1.

Case cp1 k1 cp2 k2 CR

(J/kg-K) (W/m-K) (J/kg-K) (W/m-K) (m2 K/W)
1 502.48 162.7 502.48 162.7 1.00E04
2 502.48 162.7 502.48 162.7 1.00E05
3 502.48 16.27 502.48 16.27 1.00E04
4 502.48 1.627 502.48 1.627 1.00E05
5 502.48 1.627 502.48 1.627 1.00E04
6 502.48 162.7 502.48 1.627 1.00E05
7 502.48 162.7 502.48 1.627 1.00E04
8 502.48 162.7 502.48 1.627 1.00E05
9 5024.8 162.7 5024.8 162.7 1.00E04
10 5024.8 162.7 5024.8 162.7 1.00E05
11 502.48 162.7 502.48 162.7 1.00E04
12 502.48 162.7 502.48 162.7 1.00E05
13 50.248 162.7 50.248 162.7 1.00E04
14 50.248 162.7 50.248 162.7 1.00E05
15 5024.8 162.7 50.248 162.7 1.00E04
16 5024.8 162.7 50.248 162.7 1.00E05

Table 4
Particle-wall interaction cases applied to system 2.

Case cp1 k1 cp2 k2 CR

(J/kg-K) (W/m-K) (J/kg-K) (W/m-K) (m2 K/W)
W1 502.48 162.7 502.48 162.7 1.00E05
W2 502.48 1.627 502.48 1.627 1.00E05
W3 502.48 162.7 502.48 1.627 1.00E05
W4 50.248 162.7 50.248 162.7 1.00E05
W5 502.48 162.7 50.248 162.7 1.00E05
456 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466

 
@T
qcp ¼ r  ðkrTÞ þ Q_ ð21Þ
@t
Insulated boundaries are applied in the current model since the
only heat transfer mechanism being studied is heat conduction dri-
ven by the temperature differences within and between the solid
bodies and so the only external energy transfer across a boundary
is the contact conduction handled through the modeling method
with its volumetric heat sources. Convective and radiative effects
can be readily incorporated in a multimode heat transfer model
independent of the source based method. Also, if moving objects
are considered, the new contact and thermal exchange conditions
can be found before conditions are found at the next timestep.
The process is then repeated.
A ‘‘boundary” layer of quadrilateral elements is applied at the
bounding surface of the particles or the wall, while the bulk mate-
rial is filled with triangular elements (Fig. 10). The quadrilateral
elements produce a more uniform distribution of the volumetric Fig. 11. Comparison of surface flux and volumetric heat source assignment of
contact heat exchange.
sources that account for the conductive contact, throughout the
contact domain. However, a purely quadrilateral element mesh
fluxes across fluid-solid boundaries within a standard computa-
type is not allowable since the current mesh motion technique uti-
tional model include the conduction, convection, and radiation
lized by FLUENTÓ requires a triangular mesh. As the ultimate goal
heat modes based on the existing conditions and geometry. The
is to use this method with moving objects, triangular elements
volumetric sources ensure the contact based heat conduction is
must be used to mesh the bulk object geometry excluding the
independent of these other heat transfer phenomena. Second, FLU-
boundaries.
ENTÓ storage of user parameters in computational cells rather
The basic modeling method described can be applied to three
than surfaces is more stable and reliable for moving and deforming
dimensional systems as there are no special restrictions on the
meshes where individual cells may divide or collapse or move from
technique. A more general formulation of the surface normal
one computational partition to another. Applied in a thin layer of
may be used and instead of quadrilateral elements around the
elements around the object boundary, the volumetric sources pro-
boundaries of the solid bodies, hexahedral elements are needed,
duce no significant difference from a direct heat flux application as
with tetrahedral elements in the bulk of the particles/objects.
shown in Fig. 11.
Two dimensional models were selected to demonstrate the utility
For each model case, a corresponding simulation is conducted
of the model over a range of operating conditions, since more sim-
using the standard finite volume commercial code with the ‘‘con-
ulations can be conducted within a given set of computational
tacting” elements actually in contact via an interface boundary
resources.
condition, so no volumetric sources are used. The data from these
Mesh sensitivity studies with this basic mesh arrangement
standard models are used to evaluate the performance of the new
demonstrate that reducing the mesh size and accompanying time
modeling method developed in this paper.
step by a factor of two from that depicted in Fig. 10 results in less
than a 1% change in the local temperature distribution along the
line y = 0 at a time of 0.017 s, sufficient for temperature gradients
to develop in both particles. Hence, a mesh of 168 elements along 5. Results
the circumference of each particle is used (Fig. 10), approximately
0.035 mm size elements for the 1 mm radius particle. With the The temperature field results associated with the three major
particle-wall, system two arrangement, there are 252 elements parametric study sets demonstrate the utility of the novel method
along the length of the wall and 77 across its thickness. for replicating the thermal conductive contact between two bodies.
A least squares cell based discretization method is specified
with second order upwinding for the energy equation and second
order discretization in time. A time-step of 1  1006 s is used in
the current models.
A volumetric energy source rather than a heat flux is used to
model the heat exchange for two main reasons. First, the heat

Fig. 12. Temperature distribution indicating the effect of variation of contact

resistance for two particles in contact (time = 0.017 s, y = 0). Contact resistance
Fig. 10. Two particles in contact. shown in m2 K/W.
 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466 457

5.1. Particle-particle contact with varying thermal contact resistance of this paper. The contact resistance controls the rate of heat transfer
between the particles. As the contact resistance drops, the restric-
In the first set of studies, particle-particle thermal contact for tions on the heat exchange between the particles diminishes. The
the system in Fig. 7 is studied over a range of contact resistances decrease in the contact resistance leads to greater temperature
with constant particle size, material properties, and initial temper- change in a particle, a large zones experiencing a temperature
atures. The particles both have the high thermal conductivity value change, and a decreased temperature ‘‘jump” at the contact.
of 162.7 W/m K and the mid-value of the specific heat at The results of the simulations with the new modeling technique
502.48 J/kg K. Contact resistances range from 0 to 1.0  103 are consistent with the expected effects of the contact resistance
m2 K/W. The particle on the left (Fig. 7) is assigned an initial, uniform variation. Fig. 12 depicts the temperatures along the line y = 0
temperature of 700 K and that on the right 300 K. A representative (see Fig. 1) at a time of 0.017 s for the various cases using the
animation of the temperature field can be found in the online version new model (marked model) and those produced using the

a) CR=1.0x10-3 m2K/W

700K

300K
-5 2
b) CR=1.0x10 m K/W

c) CR=0.0 m2K/W
Fig. 13. Temperature field contour plots indicating the effect of variation of contact resistance for two particles in contact (time = 0.017 s). Temperature in K.
458 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466

standard model, marked actual with the dashed line, with particle thermal conductivity on the temperature distributions is
the particle interface at x = 0.001 m. At a contact resistance of investigated for two contact resistances. The specific parameters
1.0  103 m2 K/W, the temperature change of particle 1 is minimal. studied are Cases 1 through 8 in Table 3. The thermal conditions
At the 1.0  104 m2 K/W contact resistance level, the maximum predicted by the current modeling method and those predicted
temperature change is under 2% at 0.017 s. The discontinuity in using the standard finite volume method are compared along the
the temperature across the interface is sharp. With a relaxed contact line at y = 0 at a time of 0.017 s as shown in Fig. 14. Temperature
resistance of 1.0  106 m2 K/W, the maximum temperature change contours, also at the 0.017 s time, in Fig. 15 help to visualize the
is 23% and at zero contact resistance, the maximum change is 29% impact of the changes in the thermal conductivity on the temper-
with a continuous temperature distribution across the interface ature field. The thermal conductivity of the material controls the
and temperature variation from the initial conditions throughout propagation or diffusion of the contact heat transfer through the
both particles. The data for the temperature distributions in particle material. Higher thermal conductivity induces higher rates
Fig. 12 are consistent with the characteristics of the contour plots of thermal diffusion and a larger zone thermally ‘‘affected” by the
in Fig. 13. contact.
Across all of the contact resistance values tested, the results of More of the particle material experiences a temperature change
the current model show excellent correspondence to the results in Cases 1/2, with the highest thermal conductivity value and the
from the standard commercial model where the fixed particles two different contact resistances, than for any of the other cases
are actually in ‘‘computational” contact as seen in Fig. 12. The studied, though the temperature change levels are not the greatest.
new modeling method properly replicates the effects of the phys- The higher thermal contact resistance of 2  104 m2 K/W severely
ical phenomena associated with conductive heat transfer between restricts the level of temperature changes, but some slight changes
two contacting particles for a range of contact resistances. can be observed across the entire particle along the line y = 0
(Figs. 14a and 15a). Particle 1 experiences a maximum 1.1%
5.2. Particle-particle contact with varying thermal conductivity decrease, and particle 2 a maximum 1.7% increase. For the lower
contact resistance case of 2  105 m2 K/W, Case 2, the tempera-
In the next set of studies, the density and specific heat of the tures have been altered by the contact across nearly the entire par-
two particles are held fixed and the effect of the variation of ticle diameter, with the largest shift in temperature near the

Fig. 14. Temperature distribution indicating the effect of variation of thermal conductivity for two particles in contact (time = 0.017 s, y = 0).
 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466 459

700K

a)Case2 k1=162.7 W/mK k2=162.7W/mK b)Case 4 k1=16.27 W/mK k2=16.27W/mK

cp1=502.48 J/kgK cp2=502.48J/kgK cp1=502.48 J/kgK cp2=502.48J/kgK

300K

c)Case 6 k1=1.627 W/mK k2=1.627W/mK d)Case 8 k1=162.7 W/mK k2=1.627W/mK

cp1=502.48 J/kgK cp2=502.48J/kgK cp1=502.48 J/kgK cp2=502.48J/kgK
Fig. 15. Temperature field contour plots indicating the effect of thermal conductivity for two particles in contact (time = 0.017 s). CR = 1  105m2 K/W.

particle-particle interface (Figs. 14a and 15a). Particle 1 experi- change are found due to the low conductivity in particle 2. At par-
ences a maximum 7.5% decrease, and particle 2 a maximum ticle 2, with the lowest thermal conductivity, the heat exchange is
12.5% increase. restricted to the immediate vicinity of the particle-particle inter-
In Cases 3/4, the thermal conductivity of the two particles face and conditions are similar to those in Cases 5 and 6 at the min-
decreases by a factor of 10. The portion of the particle undergoing imum thermal conductivity, with minimal temperature change
a temperature change progresses no more than one-third through further from the interface.
the particle diameter for Cases 3/4, significantly lower than for the Clearly, in Fig. 14, the results from the proposed modeling
higher thermal conductivity cases (Figs. 14b and 15b). With the method match the results from the standard finite volume model
slowed heat transfer progression, the extent of the temperature to within 1%.
decrease of the hotter particle (1) and the increase of the cooler
particle (2) is slightly higher than for the higher conductivity cases. 5.3. Particle-particle contact with varying specific heat
For Case 3, with highest thermal contact resistance of 2  104 m2
K/W, the peak temperature changes are about 1.4% at particle 1 In the next set of studies, the specific heat of the particle mate-
and 1.9% for particle 2. With the lower thermal constant resistance rials are varied with all other properties/geometries fixed, for two
of 2  105 m2 K/W, the maximum temperature decrease in parti- contact resistances and the maximum thermal conductivity (Cases
cle 1 is 7.8% and the maximum particle 2 increase is 18.3%. A sim- 9 through 16 in Table 3). The local temperature distribution along
ilar trend occurs with a thermal conductivity decrease by another the y = 0 line at 0.017 s are provided in Fig. 16 and temperature
factor of 10 in Cases 5/6 with only about 15% of the particle diam- contours at 0.017 s are in Fig. 17. The specific heat value controls
eter undergoing a temperature change from the initial conditions the thermal storage capacity of the material, so that for the same
in Figs. 14c and 15c. However, the temperature change levels grow heat input, a lower specific heat will result in higher temperature
with particle 1 temperature decrease increasing to 3.3% and 26.6% changes in a system.
for Cases 5 and 6 and peak temperature increases at particle 2 ris- For the highest specific heat level sets (Case 9/10), the temper-
ing to 17.8% and 28.5% respectively. atures remain close to their initial values due to the large thermal
Finally, with the high thermal conductivity in particle 1, the hot capacity and the spread in the temperature change is only through
particle side of the system, and the minimum thermal conductivity 25% of the particle diameter. The decrease in the particle 1 temper-
in particle 2, the cold side of the system, in cases 7/8, a combina- ature and increase in the particle 2 temperature are only 0.3% and
tion of the trends described above is observed (Figs. 14d and 0.7% respectfully for Case 9 with the higher resistance (Fig. 16a).
15d). A slight temperature change (indicated by the shallow down- The lower contact resistance case produces larger temperature
ward slope in the negative x range) is seen throughout the diame- changes of a maximum of 4.5 and 11.1% for particle 1 and 2 for
ter of particle 1 (Fig. 14d), as the energy exchange is ‘‘spread” over Case 10 as in Figs. 16a and 17a.
the majority of particle 1 mass/volume and the conditions are sim- When the specific heat is decreased by a factor of ten in Cases
ilar to Cases 1 and 2 in particle 1. Lower levels of temperature 11/12, the thermal capacity of the particles drops, leading to more
460 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466

a) k1=162.7 W/mK k2=162.7W/mK cp1=5024.8 J/kgK c)k1=162.7 W/mK k2=162.7W/mK cp1=50.248 J/kgK
cp2=5024.8J/kgK (Case 9 and 10) cp2=50.248J/kgK(Case 13 and 14)

b)k1=162.7 W/mK k2=162.7W/mK cp1=502.48 J/kgK d)k1=162.7 W/mK k2=162.7W/mK cp1=5024.8 J/kgK
cp2=502.48J/kgK(Case 11 and 12) cp2=50.248J/kgK(Case 15 and 16)

Fig. 16. Temperature distribution indicating the effect of variation of specific heat for two particles in contact (time = 0.017 s, y = 0).

700K

a)Case 10 k1=162.7 W/mK k2=162.7W/mK b) Case 12 k1=162.7 W/mK k2=162.7W/mK

cp1=5024.8 J/kgK cp2=5024.8J/kgK cp1=502.48 J/kgK cp2=502.48J/kgK

300K

c) Case 14 k1=162.7 W/mK k2=162.7W/mK d) Case 16 k1=162.7 W/mK k2=162.7W/mK

cp1=50.248 J/kgK cp2=50.248J/kgK cp1=5024.8 J/kgK cp2=50.248J/kgK

Fig. 17. Temperature field contour plots indicating the effect of specific heat for two particles in contact (time = 0.017 s). CR = 1  105 m2 K/W.
 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466 461

resistance, to 7.5 and 18.3% for particles 1 and 2 respectively

(Figs. 16a and b and 17a and b).
In Cases 13/14, the particle specific heat is decreased further by
a factor of 10. A more substantial temperature change is noted
across each of the particles (Figs. 16c and 17c), with at most a
1.8% temperature change for the higher contact resistance and a
15.7% temperature change for the lower thermal contact resistance
for particle 1 and a 4.3% and 36% maximum temperature increase
at particle 2 for the two contact resistance cases.
Finally, with particle 1 at the highest specific heat and particle 2
at the lowest specific heat (Cases 15/16), particle 1 experiences lit-
tle temperature change (Figs. 16d and 17d), while the particle 2
experiences a more significant temperature change throughout
the particle. For Case 15, a maximum 0.6% percent temperature
decrease at particle 1 and 5.3% increase at particle 2 develops with
values of 3.5% and 48.5% at particle for Case 16 with the reduced
contact resistance. The conditions similar to those in Cases 13/14
with the low specific heat are seen for particle 2, though with a
greater temperature increase since the temperatures in particle 1
with the larger specific heat remain higher than in the Cases 13/14.

5.4. Particle-wall contact with varying thermal conductivity

Next, particle-wall contact studies are conducted, varying the

thermal conductivity in the three cases W1 through W3 in Table 4.
Fig. 18 depicts the temperature at y = 0 and the good agreement
between the current and the standard models. Fig. 19 shows con-
tour plots at 0.017 s. The higher thermal conductivity allows for
swifter propagation of the heat away from the interface. While,
the geometric layout and mass of the two contacting objects differ
in these cases, the basic trends relating the temperature field to the
thermal conductivity follow those for the particle-particle contact.
In general, the thermally affected zone remains closer to the
particle-wall interface as the thermal conductivity decreases and
the temperature gradients/changes in the temperature are larger
in the immediate vicinity of the particle-wall interface due to the
reduced propagation rate of the heat flow. The wall behaves nearly
as a semi-infinite solid with the temperature change restricted to a
small distance from the particle-wall interface and higher temper-
ature gradients near the interface as the conductivity decreases
(Fig. 18a–c).

5.5. Particle-wall contact with varying specific heat

Finally, the specific heat values of the particle and wall are var-
ied. The cases considered are W1, W4, and W5 in Table 4. Fig. 20
contains the local temperature distribution at the line y = 0 and
Fig. 21 contains contour plots of the temperature field, both at
0.017 s. As the specific heat decreases, the temperature changes
experienced by the wall and particle grow just as with the two par-
ticles. However, since the mass of the wall is larger than that of the
particle, the effect of the variation in the specific heat is generally
less pronounced at the wall.
First, the specific heat of both objects is reduced by a factor of
10 moving from case W1 to W4. The maximum temperature drop
for the particle 1 is 17.8% for W4 compared to 10.1% for W1, with a
14% decrease in the temperature at the outer diameter for W4 and
Fig. 18. Temperature distribution indicating the effect of variation of thermal
conductivity for particle-wall contact (time = 0.017 s, y = 0). negligible change for W1. At the wall, the maximum temperature
rise remains at about 16.6% for both cases, but for the W4 case,
the temperature rise reaches the surface at x = 0.003 m
(Figs. 20a/b and 21a/b).
extensive changes in the temperature field due to the contact with With the higher specific heat assigned to the particle and the
some, though slight, temperature changes occurring throughout lower specific heat assigned to the wall in Case W5, the tempera-
the particles. With the higher contact resistance in Case 11, the ture change of the particle is minimal, except near the interface
maximum temperature change rises to 0.9 and 1.3% for particles as in Fig. 20c, while the temperature in the wall is similar to that
1 and 2 respectively and for Case 12, with the lower contact of W4 with a slightly higher ‘‘contact” temperature. The maximum
462 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466

a) k1=162.7 W/mK k2=162.7W/mK cp1=502.48 J/kgK cp2=502.48J/kgK(Case W1)

700K

300K

b) k1=1.627 W/mK k2=1.627W/mK cp1=502.48 J/kgK cp2=502.48J/kgK( Case W2)

c) k1=162.7 W/mK k2=1.627W/mK cp1=502.48 J/kgK cp2=502.48J/kgK (Case W3)

Fig. 19. Temperature field contour plots indicating the effect of thermal conductivity for particle-wall contact (time = 0.017 s). CR = 1  105 m2 K/W.
 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466 463

Overall, the studies for both the particle-particle and the

particle-wall interactions have exhibited good comparisons
between the new model and standard model predicted
temperatures.

5.6. Complex geometry application with stack of particles between two

conducting plates

In order to further illustrate the implementation of the model-

ing technique described in this paper, the heat conduction through
the set of 11 particles and two plates in contact was simulated for
the system shown in Fig. 22a. Each row of particles begins at a dif-
ferent temperature with the top row at 1300 K, the middle row at
900 K, and the lower row at 500 K. The plates are initially at 300 K.
The conduction heat flow through the contacting bodies is simu-
a)k1=162.7 W/mK k2=162.7W/mK cp1=502.48 J/kgK lated using the method explained in this work. The temperature
cp2=502.48J/kgK (Case W1) contour at 14.5 ms in Fig. 22b shows the local temperature rise
near the contacting zones in upper plate and the corresponding
cooling of the zones in the upper row of particles from above
and from the cooler particles in the mid-row. The localized heating
of portions of the tops of the mid-particles and cooling of the lower
portions can also be observed in the contour lots. The local temper-
ature around the circumference of the second particle in the mid-
dle row at the times indicated is shown in Fig. 22c. The progression
of the temperature changes due to the contacting particles can be
clearly seen. Hence, the model can be applied to multi-particle and
multi-wall systems as well.

6. Conclusions

The results presented and discussed in this work demonstrate

that the alternative modeling technique developed specifically to
facilitate the simulation of particle – flow mechanical and thermal
interactions in particle laden flow can readily predict the temper-
b)k1=162.7 W/mK k2=162.7W/mK cp1=50.248 J/kgK
ature fields that develop due to thermally contacting particles
cp2=50.248J/kgK(Case W4) and particle – walls with heat conduction as the only mode of heat
transfer. The higher temperature gradients and the slower thermal
propagation associated with a reduced thermal conductivity and
the higher temperature changes associated with a lower specific
heats and lower thermal mass were successfully predicted by the
model. Excellent correspondence to results of simulations with
actually contacting bodies in a standard finite volume method
were presented. The method was also applied to a multi-particle,
multi-wall system.
Such basic phenomena as material property changes, phase
change, chemical reactions, erosion and particle adhesion and
resuspension are highly dependent upon local temperature values.
The commonly used numerical methods reported in the literature
cannot capture the level of detail required to properly simulate
these effects, particularly in systems with high particle concentra-
tions, large particles, compressible flow, or high pressure, temper-
ature, or velocity gradients. The capability to model the heat
conduction based thermal interaction without contacting nodes/
c)k1=162.7 W/mK k2=162.7W/mK cp1=502.48 J/kgK elements in the computational domain, but with the fully dis-
cp2=50.248J/kgK(Case W5) cretized particles is an advancement towards the end goal of mod-
eling the conduction, convection, and radiation based thermal
Fig. 20. Temperature distribution indicating the effect of variation of specific heat interactions between moving and mechanically interacting parti-
for particle-wall contact (time = 0.017 s, y = 0). cles, the flow path surfaces, and the driving fluid flow.
With confidence gained in the modeling method from the
results in this work, the technique developed can be extended next
temperature drop in the particle is about 3.5% and the maximum to include the convective and radiative modes of heat transfer, first
temperature rise in the wall increases slightly to 25%, though the with fixed particles and then with moving and mechanically
variation from the W4 results farther from the interface is minimal. interacting particles.
464 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466

a) k1=162.7 W/mK k2=162.7W/mK cp1=502.48 J/kgK cp2=502.48J/kgK(Case W1)

700K

300K

b) k1=162.7 W/mK k2=162.7W/mK cp1=50.248 J/kgK cp2=50.248J/kgK(Case W4)

c) k1=162.7 W/mK k2=162.7W/mK cp1=502.8 J/kgK cp2=50.248J/kgK(Case W5)

Fig. 21. Temperature field contour plots indicating the effect of specific heat for particle-wall contact (time = 0.017 s). CR = 1  105 m2 K/W.
 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466 465

1400K

a)Temperature contours at inial condions

300K

a)Temperature contours aer 14.5 ms

c)Local temperature distribuons around the second mid-row parcle at indicated mes
Fig. 22. Application to 11 particle, two wall more complicated system.

Declaration of interest References

None. Ardekani, M.N., Asman, L.A., Picano, F., Brandt, L., 2018. Numerical study of heat
transfer in laminar and turbulent pipe flow with finite-sized spherical particles.
Int. J. Heat Fluid Flow 71, 189–199.
Das, S., Deen, N.G., Kuipers, J.A.M., 2017. A DNS study of flow and heat transfer
Submission declaration through slender fixed-bed reactors randomly packed with spherical particles.
Chem. Eng. Sci. 160, 1–19.
Florio, L.A., 2014. Computation of fluid-particle interactions with high-speed
This work has not been published previously and is not under compressible flows and multiple particles with deformation, plasticity and
consideration for publication elsewhere. collision. Powder Technol. 268, 19–37.
Florio, L.A., 2015a. Dynamic coupling of fluid and structural mechanics for
simulating particle motion and interaction in high speed compressible gas
particle laden flow. J. Fluids Struct. 54, 171–201.
Funding Florio, L.A., 2015b. Direct modeling of coupled compressible flow and macroscopic
particle spread and particle ejection. Mecannica 50 (9), 2257–2291.
This work was funded as an U.S. Army Armament Graduate Florio, L.A., 2017. Estimation of particle impact based erosion using a coupled direct
particle-compressible gas computational fluid dynamics model. Powder
School at Picatinny Arsenal Faculty Research project, FY2017.
Technol. 305, 625–1621.
Ga, J., Zhou, Z., Ya, A., 2016. Particle scale study of heat transfer in packed and
fluidized beds of elliptical particles. Chem. Eng. Sci. 144, 201–215.
Appendix A. Supplementary material Haddad, H., Guessasma, M., Fortin, J., 2016. A DEM-FEM coupling based approach
simulating thermomechanical behavior of frictional bodies with interface layer.
Int. J. Solids Struct. 81, 203–218.
Supplementary data associated with this article can be found, in He, L., Tafti, D.K., 2017. Heat transfer in an assembly of ellipsoidal particles at low to
the online version, at http://dx.doi.org/10.1016/j.ces.2018.07.042. moderate Reynolds Numbers. Int. J. Heat Mass Transf. 114, 324–336.
466 L.A. Florio / Chemical Engineering Science 192 (2018) 448–466

Kosinski, P., Balakin, B.V., Middha, P., Hoffman, A.C., 2013. Heat conduction during Morris, A.B., Ma, Z., Pannala, S., Hrenya, C.M., 2016a. Simulations of heat transfer to
collisions of cohesive and viscoelastic particles. Int. J. Heat Mass Transf. 58, solid particles flowing through an array of heated tubes. Sol. Energy 130, 101–
107–116. 115.
Kosinski, P., Balakin, B.V., Middha, P., Hoffman, A.C., 2014. Collisions between Morris, A.B., Pannala, S., Ma, Z., Hrenya, C.M., 2016b. Development of soft-sphere
particles in multiphase flows: Focus on contact mechanics and heat conduction. contact models for thermal heat conduction in granular flows. Am. Inst. Chem.
Int. J. Heat Mass Transf. 70, 674–687. Engineers J. 62, 4526–4535.
Lattanzi, A.M., Hrenya, C.M., 2016. A coupled, multiphase heat flux boundary Oschmann, T., Schiemann, M., Kruggel-Emden, H., 2016. Development and
condition for the discrete element method. Chem. Eng. J. 304, 766–773. verification of a resolved 3D inner particle heat transfer model for the
Liu, B., Li, Y., 2016. Simulation of effect of internals on particle mixing and heat Discrete Element Method (DEM). Powder Technol. 291, 392–407.
transfer in downer reactor using discrete element method. Powder Technol. Singhal, A., Cloete, S., Radl, S., Ferreira, R.Q., Amini, S., 2017. Heat transfer to a gas
297, 89–105. from densely packed beds of cylindrical particles. Chem. Eng. Sci. 172, 1–12.
Liu, D., Bu, C., Chen, X., 2013. Development and test of CFD-DEM model for complex Vargas, W.L., McCarthy, J.J., 2002. Conductivity of granular media with stagnant
geometry: a coupling algorithm for Fluent and DEM. Comput. Chem. Eng. 58, interstitial fluids via thermal particle dynamics simulation. Int. J. Heat Mass
260–268. Transf. 45, 4847–4856.
McCullough, J.W.S., Leonardi, C.R., Jones, B.D., Aminossadati, S.M., Williams, J.R., Yang, W.J., Zhou, Z.Y., Yu, A.B., 2015. Particle scale studies of heat transfer in a
2016. Lattice Boltzmann methods for the simulation of heat transfer in particle moving bed. Powder Technol. 281, 99–111.
suspensions. Int. J. Heat Fluid Flow, Part B, 150–165.