Heat transfer rates in sheared beds of inertial particles at high Biot numbers

We investigate heat conduction rates through the contact network of sheared granular materials that consist of particles with a density much larger than that of the ambient fluid. The tool ParScale, a finite difference solver for intra-particle transport processes, is coupled to the discrete element method solver LIGGGHTS to carry out the simulations. Heat transfer to the surrounding fluid is considered via a fixed heat transfer coefficient. We identify a combination of Biot and Peclet number as the key non-dimensional influence parameter to describe the non-dimensional heat flux through the particle bed and provide an analytical solution to calculate mean particle temperature profiles. Simulations over a wide range of particle volume fractions, Biot and Peclet numbers are then performed in order to develop a continuum model for the heat flux. We show that for Biot numbers below 10-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-3}$$\end{document} our results agree with a simple model based on a uniform particle temperature. However, for Biot numbers above a certain threshold, the conductive heat flux decreases substantially depending on the flow situation, and temperature gradients within the particle should be considered in order to correctly predict the heat transfer rate to the surrounding fluid.


Introduction
The reliable prediction of local temperatures in reactors for, e.g., CO 2 absorption, or biomass combustion is a research field of high interest. In these reactors particles undergo thermal or chemical transformation processes such as heating, drying, gasification, combustion, or reduction. The parallel sequence of these phenomena is still not fully understood [1,2], and often highly dependand on the local temperature gradient inside the particle [3]. When using materials with a low heat conductivity, these gradients become extremely important in order to predict, e.g., reaction fronts. These fronts may lead to a significant change of the particle properties. Also, the surrounding flow field is typically also affected.
In addition, an accurate understanding of heat flow through static (or slowly moving) granular materials is needed in various heat treatment processes, e.g., for clayed soils, or thermal granular compaction. These systems have been investigated in detail by Smart et al. [4], Shi et al. [5], or Vargas et al. [6]. Also, a hotly debated question is the heating of brittle rocks in fault gouges during earthquakes in the field of geophysics, for instance as discussed in Rice [7]. The rising temperature in these systems has a significant effect on the overall behaviour since it leads to a pressurization of the fluid, and might cause fluidization [7,8]. In the field of siliceous materials it is well known that rising temperatures even lead to grain melting as shown by Otsuki et al. [9]. A temperature rise in such particle beds can be seen even without any kind of external heat source, but just through the dissipation of mechanical energy, and the inability of a layer of particles to transfer heat to the environment [10].
Many researchers (e.g., [11,12]) aimed towards the development of a rheological model in the context of a continuum description of granular materials. The study of Chialvo et al. [11] can be seen as a breakthrough in these studies, since different stress regimes (and the associated transitions between these regimes) for granular flow were quantitatively described for the first time. Mohan et al. [13] showed the existence of different regimes in the context of conductive and convective liquid transport in a sheared particle bed. In this latter work similarities between the transport of thermal energy and liquid adhering to the surface of the particles have been already discussed. Thereby, a dimensionless shear rate and a Peclet number were identified as the key dimensionless parameters. Unfortunately, the work of Mohan et al. [13] relied on the assumption of zero heat flux to the ambient fluid, critically limiting the applicability of their conclusions. Recently, Morris et al. [14] made a step towards more efficient simulations by introducing correction terms that account for the artificial softening of the particles in order to correctly calculate the conductive thermal flux.
The above mentioned previous work relied on the discrete element method (DEM), and mostly focussed on the rheological behaviour of granular matter. In case heat transfer was studied, typically a rather rudimentary approach for the local particle temperature distribution was adopted. Thus, often a spatially homogeneous particle temperature (i.e., no intraparticle temperature gradient) is considered. A good example for such an approach is the work of Tsory et al. [15]. Also, heat transfer to the ambient fluid is often not accounted for. For example, Rognon et al. [16,17] did not even attempt to model the transferred heat from the particle to the ambient fluid (e.g., air), in line with the assumptions of Mohan et al. [13]. A more advanced study was perfromed by Shi et al. [5], who performed a thermally coupled CFD-DEM simulation in a rotary kiln under conduction-dominated flow conditions. Unfortunately, they used strongly simplified models for the fluid and solid phase. Most important, they also neglected intra-particle property profiles, which clearly limits the predictive capabilities of their model.
In related fields, Batchelor and O'Brien [18] obtained results for the effective conductivity of the granular material in the case of non-flowing uniform spherical particles embedded in a non-flowing matrix material. Thus, particles and the surrounding matrix were not allowed to move, and hence their analytical solution is rather limited. Shimizu [19] carried out thermal simulation using an Euler-Lagrange approach. He relied on a "thermal pipe network" to represent heat conduction due to particle-particle interactions, and failed to resolve intra-particle temperature gradients. Cheng et al. [20] were able to model the effective conductivity of a non-moving packing by performing simulations that take (i) conduction through the solid particles, (ii) conduction through the contact area between contacting particles, and (iii) the conduction via the ambient (but stagnant) fluid into account.
Only some research groups include more refined models for the temperature distribution inside individual particles. Relevant studies include the work of Feng et al. [21,22] or Oschmann et al. [23]. The former propose a so-called discrete thermal element model (DTEM), which can be readily integrated with the DEM. Unfortunately, also this previous work made strong simplifications, e.g., by neglecting direct heat transfer in the contact zones, or by relying on a pipenetwork structure to compute thermal fluxes. The model of Oschmann et al. can be seen as the ultimate model for heat conduction within a particle bed: they discretize each particle with ca. 20,000 computational nodes, and are hence able to predict the full three-dimensional temperature distribution within each single particle. Clearly, such an approach is computationally very demanding. This critically limits the size of the particle bed to be studied. Consequently, the number of parameter variations considered in the work of Oschmann et al. is too small to draw more general conclusions.
In summary, there is still very little quantitative knowledge about the thermal behaviour of particle beds that consist of particles with a non-uniform temperature distribution. This is especially true for beds that are sheared, and for which convective heat transport (i.e., transport due to individual particle motion) might be the dominating transport mechanism. Clearly, a thorough understanding of the relative rates of (i) the transferred heat to the ambient fluid, (ii) the convective heat transport, and (iii) the heat that is transported via conduction is still missing in literature.
When introducing thermal fluxes in the context of the granular flow simulations considered in our present contribution, it is important to consider the exact meaning of their definition. Critically, the term "convective flux" is fundamentally different from its meaning in classical continuum models: convective fluxes in the context of the present work are caused by random motion of individual particles. This is in contrast to the convective transport of thermal energy due to macroscopic (i.e., particle-average) motion: this "macroscopic" convective transport can be modeled with relative ease in continuum models [24]. Hence, this mode of thermal transport is not discussed in the present work. Another term is the transferred heat to the ambient fluid: in our present contribution this refers to the heat exchanged between the fluid and individual particles. At this point we reference the reader to our presentation in Sect. 2.4 which provides a more holistic view of all relevant thermal fluxes. Most important, these fluxes are the main output of our simulations, and can be subsequently used in continuum models to describe thermal transport in granular materials and fluid-particle systems. Unfortunately, rigorous closures for these fluxes are still missing in literature, currently limiting the application of continuum models to engineering problems. In this study, we focus on suspensions characterized by a large Stokes number (St = 2/9 ρ p ud p /μ f ) which are relevant for a number of practical applications involving gas-particle suspensions [25].
It is our goal to improve this situation by coupling DEMbased simulations to a tool that is capable of resolving intra-particle temperature profiles: our tool ParScale (see [26,27]) allows us to efficiently solve a one-dimensional heat conduction equation within each particle. This enables us to provide answers to questions related to the occurence of situations in which intra-particle temperature profiles are significant. Clearly, this is of utmost importance for the correct selection of an appropriate simulation model for a certain fluid-particle heat transport problem. Specifically, we analyse the relevance of intra-particle temperature gradients by probing an adequate space of dimensionless parameters, most important the Biot and Peclet number. This allows us to identify critical dimensionless system parameters for which heat conduction inside the particle becomes limiting. Most important, the use of dimensionless parameters allows us to draw conclusions that are applicable to a wider range of real-world applications compared to previous work. Specifically, our contribution is meant to extend the work of Rognon et al. [16] to (i) a wider range of volume fractions, (ii) to fluidparticle systems with heat transfer to the ambient fluid, and (iii) includes the temperature distribution inside the particles. Same as Mohan et al. [13] we attempt to collapse the predicted conductive flux into one curve. This is done to lay the foundation for the development of a unified continuum model for predicting thermal transport in high-Biot number fluid-particle suspensions.
We note in passing that including heat transfer to the ambient fluid requires the specification of a heat transfer coefficient (details will be discussed in Sect. 2.4), or its dimensionless counterpart the Nusselt Number. Ranz and Marshall [28] developed one of the first Nusselt number correlations for thermal transfer between spheres and ambient fluid. Clearly their correlation is not applicable in dense particle flow which we focus on in the current manuscript. Gunn [29] proposed a still widely-used correlation for the Nusselt number in the porosity range of 0.35-1.0 for various Reynolds numbers. Thus, using Gunn's results it is straight forward to estimate the heat transfer rate, and consequently also the Biot number, in any given fluid-particle flow systems.
Our paper is structured as follows: in Sect. 2 we describe the simulation method, including the governing equations for the DEM and the intra-particle heat transport model. Furthermore, we provide details about the coupling sequence between LIGGGTHS and ParScale. Section 3 summarizes our theoretical analysis of transport of thermal energy in a sheared particle bed. Section 3.2 considers effects due to heat conduction with the particles, and Sect. 3.3 aims on presenting results for limiting cases in order to better interpret our numerical data in Sect. 4. These subsections hence constitute the first part of our results. The second part of our results is placed in Sect. 4 where the sheared particle simulation setup is described in detail, including the effect of all relevant dimensionless numbers. Section 5 views our results in light of previous findings, and comments on possible future applications of ParScale. A consistent list summarizing the nomenclature used, including all quantities of secondary importance which are not mentioned after their first appearance in the paper, is available in "List of symbols".

Simulation method
An overview of all relevant transport equations used in our simulations is provided in Table 1. Also, we have summarized the corresponding unknown simulation variables in Table 2. Note that some of the equations are discussed in more detail in the next Section. We want to emphasize that we solve the transport equations for the fluid only in Sect. 3.2, while in all other parts of our contribution we consider (i) the fluid to be held at a fixed temperature, and consider (ii) a fixed particle-fluid heat transfer coefficient. Also, as implied by our title, we consider inertial particles, i.e., particle motion is not affected by fluid motion. This is done to avoid the large computational effort needed to resolve local fluid flow and temperature, which would make our study unaffordable considering our current computational resources.

Particle flow model
The open-source software package LIGGGHTS [13,30] is a well established discrete element method-based tool and Table 1 Covering equations for particle flow (Newtons equation of motion), particle flow (Continuity equation, Navier-Stokes equation, scalar transport for e.g., temperature) and equations regarding internal thermal transport in spherical coordinates in one-dimensional form Internal thermal transport Particle temperature T was used in combination with a spring-dashpot model [30] in the present work. Note that we do not solve any transport equations for the fluid, but consider (i) the fluid to be held at a fixed temperature, and consider (ii) a fixed particle-fluid heat transfer coefficient. Specifically, the following contact force models in the tangential (i.e., f t i,j ) and normal direction (i.e., f n i,j ) were applied:

Reference system Equation
with: Thereby δ ij indicates the normal overlap, which is positive in case of overlapping particles. Below we summarize the most important quantities that are relevant for the interpretation of our work. For the meaning of other symbols we refer to the Nomenclature. The characteristic contact time t co of the above linear spring-dashpot model is defined as: with This contact time t co limits the time step t for the integration of Newtons equation of motion. The coefficient of restitution e n is defined as and was adjusted to 0.9. In our simulations, the spring stiffness and damping coefficient in the normal and tangential direction are set equal. In order to obtain a certain restitution coefficient, and a certain dimensionless shear rate γ * , we adjust η and k. For more information on the used parameter ranges and routines to calculate, e.g., the effective Young's modulus, we refer the interested reader to Sect. 1 and "Appendix".

Transport within a particle
In order to calculate intra-particle properties, a tool called ParScale (short for Particle Scale) [31] was developed. It is capable of resolving individual property profiles (e.g., that of the temperature, or a dissolved chemical species) inside a particle. Specifically, ParScale relies on the so-called Method of Lines (MoL) as outlined in the next paragraph.
The resulting system of equations is solved by an already established multi-step solver called CVODE. This enables ParScale to perform flexible and stable simulations, even in case of fast changing environmental conditions, or (numerically) stiff systems that are typical for non-isothermal reactive porous particles. The next sections outline the governing equations, provide more detailed information on the solution procedure, and give an overview of the parallel coupling strategy.

Governing equations
In this paper only heat conduction inside non-porous particles is taken into account. Note, the current development state of ParScale already allows the treatment of mass transfer, multiple chemical reactions, as well as phase change phenomena (e.g., evaporation) taking place inside porous particles. However, these phenomena are not accounted for in our present study.
The basic approach of the MoL is to (i) discretize the transport equation related to all intra-particle properties in space, and (ii) then integrate the discretized equations in time. At this stage, a one-dimensional spatial discretization in a radial coordinate system centered at the particle position has been employed. Therefore, every particle is spatially split in a fixed number of grid points. At each of these points, the differential heat balance in spherical coordinates is solved (see Eq. 8). We want to amplify that the thermal boundary condition of the heat transfer equation solved in ParScale is a combination of the transferred heat flux to the surrounding q trans and the per-particle heat conduction rate (which is calculated by summing over all particle-particle contact points, similar to what is done in Eq. 10). Since we are only tracking solid, non-porous particles, we can assume for the particles' mean heat conductivity λ p = const. After introducing the thermal diffusivity a, the differential heat balances reads: For discretisation in physical space, a second order central differencing scheme is used. A verification case, i.e., the convective cooling of a sphere, was used to ensure a correct implementation of the discretization schemata, and the integration of ParScale with the solver CVODE. For details of the differencing scheme, as well as results of the verification studies, the reader is referred to "Appendix" and Sect. 1.

Code architecture and parallel coupling strategy
ParScale is developed in a C++ environment, relies on a hierarchical class structure, and is packaged in a library that can be easily integrated with other particle-based simulation tools. Most important, ParScale is equipped with wellstructured data containers for the easy transfer of particle information during a parallel simulation run. Note, that data containers are separated from an external particle simulator, in order to ensure a maximum level of compatibility with different simulators and parallelization strategies. In summary, ParScale is scalable, flexible, and can be adapted to the users needs by the simple addition of user-defined classes. ParScale is currently hosted on a public Github repository [31], together with detailed installation instructions and a comprehensive user documentation. ParScale has two different run modes: (i) a stand-alone mode, and (ii) a coupling mode in which ParScale acts as a slave to an external DEM or CFD-DEM simulator. At the current development state ParScale features coupling to LIGGGHTS, as well as CFDEMcoupling. Example cases placed in the public repository demonstrate the usage of ParScale, contain several verification cases to prove the correct implementation, and provide some post-processing and plotting routines. Since ParScale is combined with the open-source solver LIGGGHTS in the current work, an overview of a single time step of a coupled LIGGGHTS-ParScale simulation is shown in Fig. 1. At the beginning of every time step LIGGGHTS is running as the master program. It solves for the particle movement (i.e., all DEM-based calculations) and hands over surface temperatures, actual heat transfer coefficients, and fluxes due to an eventual conductive heat transfer between contacting particles. Some other coupling options are available, e.g., setting a specific environment temperature for each particle, which is used in the present contribution. Furthermore, it is possible to reset the whole particle temperature to a constant value. After all calculations of LIGGGHTS are completed, ParScale initializes from the last know time step and reacts to changes due to the coupling of the particle motion solver. Specifically, it updates the surface temperature, calculates source terms (e.g., due to reactions or phase change phenomena), and pushes these particle properties back to LIGGGHTS for the next time step. For visualisation and post-processing purposes, additional information (e.g., core, or the volume-averaged temperature, as well as conductive, or convective heat fluxes) can be exchanged with LIGGGHTS.

Thermal fluxes
The definition of thermal fluxes is taken from the work of Mohan et al. [13] and a simplified example assemble including a description of coordinates can be found in Fig. 7. Physically seen, the convective flux (Eq. 9) is defined as the redistribution of thermal energy of the particles due to its macroscopic movement inside the simulation domain. In the contrary, the conductive thermal flux (Eq. 10) accounts for the heat exchanged between two particles in case of a collision. Drawing back to classical thermodynamics, the transferred heat flux (Eq. 11) describes the exchange of thermal energy between fluid and particle. This transferred amount of heat is mainly a function of the heat transfer coefficient α and the temperature difference between particle and ambient fluid. Note that q cond is defined for particles with equal physical properties regarding thermal conductivity and is calculated for the sum of all contacts in the domain with the volume V .
The conditional claim of a fully periodic, sheared bed with Lees-Edwards boundary conditions is that the domainaveraged particle velocity equals zero. When considering the sheared particle ensemble used in our work (see Sect. 4), it is important to note that we consider three different "thermal" regions as shown in Fig. 7 in Sect. 4. Most important, we only consider the center region (i.e., region B in Fig. 7) to calculate the fluxes. In such a case one has to correct the calculation of the convective flux as discussed in the next paragraph. Even though the above conditional claim related to the domain-averged particle velocity is fulfilled for the full box, this does not apply for our center region in which the thermal calculations are performed. Since we fix the temperature of the particles located inside the top and bottom region to enforce a predefined temperature gradient, the number of particles in these regions fluctuates over time. Also, dividing the domain into three regions causes a temporal fluctuation of the number (as well as the average speed) of the particles in each region. Thus, we see a drift convective flux due to the average macroscopic movement of the particles in the center region, which is characterized by the velocity v y,drift .
or explicitly calculated As mentioned above, the number of particles in our center region is not fixed and can change due to the deformation of the box. Consequently, we need to weight our already corrected convective flux with the number of particles in the region when computing a time-average flux: For the reference flux we choose the conductive flux in case our computational domain is completely filled with solid material. Thus, the reference heat flux is calculated as: where ∂ T /∂ y is the temperature gradient in the gradient direction of the shear flow.

Dimensional analysis
The key dimensionless flow parameters, i.e., the coefficient of restitution, the friction coefficient, the dimensionless stresses, the particle volume fraction, and the dimensionless shear rate have been well documented by Chialvo et al. [11], and Mohan et al. [13]. In agreement with Chialvo et al. [11] we define a scaled shear rate aṡ According to the experiments by Nordstrom et al. [32], and recent simulations [10,33,34], simple expressions that correlate the above defined dimensionless quantities can be found. Specifically, the scaled stress-shear rate data can be collapsed into curves by appropriately adjusting the critical volume fraction (φ c ), as well as a set of exponents as shown below. This scaling will be an inspiration for the re-scaling of the conductive heat flux in the present contribution (see Sect. 4). Thus, of key relevance are the three resulting stress regimes (here illustrated by the pressure), which are characterized by the following expressions: Previous work only considered the conductive and convective heat transport rate in a sheared, but uncooled, particle bed (see, e.g., [17]). As pointed out by Mohan et al. [13], the Peclet number is the most influential parameter for heat transport in granular materials and reads: Clearly, Pe characterizes the relative rates of convective and conductive heat transport within the particles in the absence of heat transfer to the ambient fluid. In the present work, however, we consider exactly this heat transfer rate to the ambient fluid, which is characterized by the heat transfer coefficient α. Therefore, we must introduce a new dimensionless quantity, which we choose to be the Biot number defined as: The Biot number can be interpreted as the ratio of (i) the intra-particle and (ii) the external resistance to heat transport. Thus, for high Biot numbers we expect intra-particle temperature profiles to be significant. We will see that also the Peclet number affects intra-particle temperature profiles, since conductive transport between particles (caused by particle-particle collisions) leads to a heat up of the outer shell of the particles.
As an alternative to the Peclet number, one could define a dimensionless convective heat transport rate by using the transferred heat flux (to the ambient fluid) as the reference, i.e., Π can be interpreted as a scaled length over which convective transport is able to balance the cooling (or heating) caused by the ambient fluid. This dimensionless number is useful in situations in which one wants to probe the effect of the the particles' heat conductivity on the results: for fixed values of Π one can illustrate the effect of heat conduction within the particles by varying the Biot number. Note, that for fixed Pe a variation of Bi does not lead to such a transparent illustration, since also Pe is a function of the particles' heat conductivity, and both dimensionless numbers are zero in the limit of infinitely fast heat conduction within the particles. Since Π is proportional to Pe and Bi, it is straight forward to interpret our results also in light of this new dimensionless quantity. We have already attempted this in most parts of the present contribution.

Effect of heat conduction within the particle
The goal of this Section is to theoretically investigate the cooling of a sphere under different conditions, and justify the assumption of spherical symmetry made in the current work. As pointed out by Schmidt et al. [35], as well as Oschmann et al. [23], the computational cost for a fully resolved particle is enormous: calculation time and data storage increase by a factor of about 1000 when using three-dimensionally resolved simulations. Fortunately, a one-dimensional discretization (as followed in the current work) is often an excellent approximation (in case of a sphere) when considering heat exchange in a certain Reynolds number range (see Nikrityuk et al. [36]). To support this finding, we study the transient cooling of a sphere, and consider different Reynolds numbers. We then compare the temperature change over the radius with the data provided by Nikrityuk et al. [36]. The latter performed a fully three-dimensional simulation, including a resolved flow and temperature profile simulation in the fluid surrounding the particle. For the simulation details we refer to Nikrityuk et al. [36]. We have summarized the key dimensionless parameters in the caption of Fig. 2.
In general, an excellent agreement between the DNS data and simulations performed with ParScale over a wide range of dimensionless times, i.e., Fo numbers, can be observed. This is especially true for large Re numbers, where only marginal differences between DNS data and predictions by ParScale are observed. As already pointed out by Nikrityuk et al. [36], the maximum difference between a fully three-dimensional simulation and a one-dimensional simulation is between 2% and ≈ 5%. To support this finding, we investigated conjugate heat transfer associated with the flow around a single sphere under significant internal heat transport limitation (see Fig. 3 for details). We listed all boundary conditions in Table 3 and governing equations in Table 1. As can be seen, the temperature profiles inside the particles approach spherical symmetry for higher Reynolds numbers. The reason for this behaviour is that the the temperature distribution in the ambient fluid becomes more uniform, since   Fig. 3) a certain temperature gradient in the fluid in the main flow direction is observed. This temperature gradient leads to a shift of the maximum intra-particle temperature towards the rear stagnation point, causing a significant temperature gradient in the azimuthal direction. Fortunately, most industrial applications are concerned with improving external transport limitations, and hence operate at Reynolds numbers larger than unity. Also, a closer inspection of the global temperature difference in the left panel of Fig. 3 indicates a temperature fluctuation of about 0.17 dimensionless units. This is because Bi < 1, i.e., the external heat exchange rate is typically much smaller than the heat conduction rate in the particle in case of low Reynolds number flows. In summary, the results in Figs. 2 and 3 support our assumption of a spherically-symmetric temperature profile within particles for most gas-particle systems. Of course, this statement only holds for flows involving single spherical particles, and cannot be easily generalized to non-spherical particles, or dense particle ensembles. Thus, we currently cannot provide a rigorous justification for spherical symmetry in dense granular flows where particle-particle heat conduction occurs. Fortunately, this mode of thermal transport (i.e., heat conduction through the contact areas) is typically small: an analysis of the amount of heat exchanged due to conduction (denoted as Q cond ) and transfer to the surrounding fluid (i.e., Q trans ) on a per particle basis helps to support this argument. Specifically, one finds that the ratio of these heat exchange rates is (based on the expressions for the heat fluxes introduced in Eqs. (10) and (11): Here we have introduced as a parameter that characterizes a typical overlap area, which is typically much smaller than unity. Also, we have denoted with T co a mean temperature difference of the contacting particles. If we now accept that the term T co /T i − T f,i is of order unity, we find Thus, the relative importance of the conductive heat flux decreases with increasing Biot number. It is now obvious that the heat exchange rate due to particle-particle contacts will not impact the shape of intra-particle temperature profiles, since the latter are significant only for Bi > 1. Thus, one can readily transfer the findings related to intra-particle temperature profiles in a single particle to a cooled bed of particles in case one considers Biot numbers of order unity or larger.

Analytical solution for heat conduction in a cooled particle bed
In the following we present an analytical approach to calculate the temperature distribution in a sheared bed with a constant ambient fluid temperature T f , and study the limiting cases with respect to the Biot number. To the best of our knowledge no analytical solution for the mean temperature profile in a cooled bed of moving particles exists. To approach this issue, we first calculate the mean particle temperature distribution inside a static bed of cooled particles. Therefore, we use the particle heat conductivity λ p , as well as simple geometrical factors that describe the morphology of the particle bed. Specifically, we assume that the bed has an average effective cross-sectional area A, and an average effective cross-sectional perimeter U (see Fig. 4 that also illustrates their local counterparts). We note that the calculation of the effective area and perimeter cannot be done by taking a simple arithmetic mean of their local counterparts. This is due the fact that in the governing equation for heat conduction in the fin (discussed in the next paragraph) the cross-sectional area A comes into play in a non-linear fashion. Consequently, A will be also impacted by the average particle-particle contact area, which depends on φ p , and hence the applied stress.
Fortunately, in what follows it is not necessary to provide an expression for the effective cross-sectional area and perimeter. We will see that considering these effective quantities is sufficient to explain the qualitative features of the mean particle temperature distribution. In case we accept the above simplification, we can reduce the problem of heat conduction inside a static bed to the problem of heat transport in and around a convectively-cooled fin (see Fig. 4, left panel). The steady-state solution of this classical heat transport problem is straightforward, and relies on the following differential heat balance: Note that the above expression only considers thermal transport through the particles, and hence avoids complexities arising due to, e.g., heat conduction via the ambient fluid as mentioned in Kuipers et al. [24]. We feel that such complexities should be incorporated by a separate model for thermal transport in the surrounding fluid, and hence can be added in the future. Normalization of the above heat balance equation by introducing leads to The first term on the right hand side of the above equation is equal to the Biot number. The second term, i.e., (g * ), represents a dimensionless metric for the bed morphology. The latter fluctuates in case the bed of particles is sheared (see our discussion below). However, it is reasonable to assume that this is irrelevant for the steady-state solution to a first approximation. We next consider thermal transport in a sheared bed. Clearly, the rate of heat transport through the bed will be increased due to the motion of particles in such a case. This will not change the structure of the above differential equation as long as we assume that particle motion in the gradient direction leads to a diffusive transport of thermal energy. In other words, we assume that the additional heat flux caused by particle motion linearly scales with the temperature gradient. Consequently, only the magnitude of the heat conductivity will change, and one could simply introduce an effective heat conductivity of the sheared particle bed, i.e., λ eff in Eq. 30. This effective heat conductivity will be some function of φ p ,γ and the particle parameters. Naturally, one would express λ eff /λ p as some function of the Peclet number, the particle concentration and the dimensionless shear rate. Hence, one can easily split off the particles' heat conductivity from λ eff and remain the general structure of Eq. 30. Thus, we simply lump the effect of particle flow, i.e., the convective heat transport due to random particle motion, and the consequent change of A into the (unknown) function g * (Pe, φ p ,γ ). Clearly, g * will decrease with increasing shear rate and increasing Peclet number (i.e., softer parti-cles moving faster), since the effective heat conductivity and the area available for conduction will increase with increasing speed of shearing (or increasing particle softness). While g * lumps a large number of phyical effects, we have not attempted to decouple these effects in the present contribution. This is done since such a decoupling is not essential to understand the qualitative features of the particle temperature profiles discussed in the next paragraphs. Also, we note in passing that in the following, we have set T f to be constant as indicated in Fig. 4.
We now use the following Ansatz for the solution of the temperature profile: with the boundary conditions This leads us to the following constants: Combining these constants with Eq. (31) one obtains the following mean particle temperature distribution: This equation does not allow us to explicitly calculate the effect of the shear rate on the resulting temperature profile. However, it provides us with an analytical solution to study certain limiting cases, namely sheared beds under the influence of high and low thermal heat exchange rates to the ambient fluid. Also, this analytical solution allows us to draw qualitative conclusions on the effect of the particle shear rate as mentioned above: faster shearing (i.e., increasing Π or Pe) will lead to a more uniform temperature profile, since convective heat transport (due to random particle motion) becomes the dominating mode of thermal transport compared to heat transfer to the ambient fluid. We note in passing that after calculating the temperature gradient from the above temperature profile, as well as integrating over the bed height, the mean heat flux based on the above analytical solution bed can be calculated as: This heat flux has been normalized with the reference flux (i.e., q s ), and requires the knowledge of λ eff /λ p or, alternatively, g * . Unfortunately, the latter quantities are unknown, and hence must be modelled (e.g., via models for the conductive and convective fluxes developed in the present contribution). In the following we study certain limiting cases, and only compare the qualitative features of the temperature profiles inside the sheared bed between the simulation and the above analytical solution. Note, that g * has been adjusted in the following such that the difference between analytical solution and simulation data is a minimum. Thus, the following comparison is only helpful to judge whether the shape of the temperature profiles can be predicted. Models for the effective heat flux (due to conduction and convection) are proposed in Sect. 4, which can be used to compute the effective heat conductivity in a straight forward manner.

Small Biot number limit
The first limiting case we study is a sheared bed without any heat transfer to the ambient fluid. This is in line with the work of Rognon et al. [16] and Mohan et al. [13]. Therefore, it is clear that for a steady-state situation the temperature profile is linear, as seen in Fig. 5. A very good agreement can be observed from the above figure. Volume-averaged and surface temperature are collapsing in one point since the Biot number is small (and in the limit of thermally insulated particles, i.e., no heat exchange with the ambient fluid, Bi = 0). Due to the low Biot number the solution is independent of the fluid temperature, the dimensionless shear rate and the volume fraction. Unfortunately, this solution does not lead to a deeper understanding of the total fluxes inside the bed and the relative contributions of convective and conductive transport rates. Thus, in the following we study a set of dimensionless parameters for which the temperature distribution is mainly governed by the heat transferred to the ambient fluid, and hence might help us to better understand the associated heat fluxes.

High Biot number limit
The second limiting case is a sheared bed which is strongly cooled by the ambient fluid, i.e., heat transfer to the ambient fluid is fast. It is obvious that in this limiting case the temperature of the ambient fluid is of key importance. Also, we expect that for Π approaching zero (at finite Bi) the convective and conductive flux will vanish, since the particle temperature approaches the fluid temperature. To illustrate this, we study high Biot number flows to highlight the mean particle temperature profiles in a sheared particle bed. For the shown set of dimensionless parameters a clear difference to Fig. 5 can be observed. The heat transfer to the ambient fluid leads to a fast relaxation of the bed temperature to the fluid temperature. Also, for even higher Biot numbers we observe a steeper gradient near the top and bottom of the bed (data not shown). This is explained by the corresponding decrease of Π = Pe/Bi, i.e., the dimensionless distance over which convective heat transport outbalances heat transfer to the ambient fluid decreases. In contrast, this gradient decreases with increasing Peclet number: fast shearing in combination with a low thermal conductivity of the particles smoothens out temperature gradients. Another key observation is that in case of low Π -number flows, the temperature gradient (and hence the thermal fluxes through the particle bed) is no longer constant over the simulation domain. Consequently, the size of the considered domain may affect the predicted mean thermal fluxes in that domain. Hence, in what follows it is essential to consider the size of the simulation domain for a correct interpretation of the results. Specifically, we now explore a wider range of dimensionless parameters to establish a more quantitative understanding of thermal fluxes through sheared particle beds (Fig. 6).

Simple shear flow
We now investigate the thermal behaviour of a sheared particle bed under different flow and thermal conditions. By applying Lees-Edwards boundary condition to a fully periodic box driven by a homogeneous shear flow, a combination of LIGGGHTS and ParScale can be used to study various influence parameters. The dimensions of the box are chosen to be H/d p = 15 with a wide range of particle volume fractions (φ p = 0.30-0.64) studied. The z-direction serves as the spanwise direction, whereas in x-direction the homogeneous shear flow is driven. In y-direction we apply a temperature gradient by fixing the temperatures in the top and bottom layer. Figure 7 summarizes this setup in the y − x plane, in which region A and C are the hot and cold region, respectively. In region B the temperature of the particles is allowed to evolve freely. Consequently we do not take the top and bottom layer into account during our post processing procedures, and just analyse thermal transport in the center region, i.e., region B. In order to study the transferred heat to the ambient fluid (note that [13,16] considered particles in vacuum, i.e., perfectly insulated), the fluid temperature is set equal to the average of the top and bottom temperature. Note, that the number of particles, as well as their mean velocity is not constant in the center region, a fact that needs to be accounted for in the analysis. Wide ranges of various parameters and their influence on the thermal fluxes (corresponding to [13,16]) are studied in what follows. An overview of these ranges is given in Table 4. These parameter ranges in Table 4 are chosen according to realistic scenarios involving solids particles in a gas. Relevant industrial applications for the high Biot number are, for example, thermally thick biomass particles (Mehrabian et al. [37]). The volume fraction covers more dilute flows (φ p = 0.3), but also includes the densest packing possible under the assumption of mono-disperse spheres (φ p = 0.64). Dimensionless shear rate and Peclet number relate to Rognon et al. [16] and account for soft and hard spheres and particle movement that can typically be found in earthquakes, landslides and general granular lubrication, respectively. Biot numbers up to Bi = 25 are typical for granular flows involving highly porous particles with low heat conductivity (e.g., wood, or polymeric materials), a typical size of d p ∼5-10 mm, and a Reynolds number of Re = 100 or higher.
In our simulations all particles are initialized with a uniform temperature, are forced to the temperatures in the top and bottom layer, and are sheared according to a predefined set of dimensionless parameters. After a certain shear deformation, the flow has converged into a statistical quasi-steady state, in which convective and conductive fluxes through the bed (i.e., region B) can be averaged and compared. Figure 8 shows typical results at such a state. In order to visualize the variety of flow situations, Fig. 8 summarizes results for certain limiting cases. Specifically, Fig. 8a shows the lowest Peclet and Biot number studied (compare Table 4). Figure 8b shows results for a flow that is dominated by the transferred flux to the ambient fluid, whereas Fig. 8c, d show results for the highest Peclet number studied and the extrema in the Biot number. As can be seen, a non-trivial distribution of particle temperatures develops in the simulation domain, depending on the combination of Peclet and Biot number. Interestingly, intra-particle temperature gradients appear for the situations characterized by a high Peclet number (see Fig. 8c, d). As expected, these gradients are most significant for the combination of a high Pe and Bi.
We next study the influence of the particle stiffness (i.e., the dimensionless shear rate) on conductive fluxes in the sheared bed over a wide range of Biot and Pelcet numbers. This analysis is motivated by the regimes observed for the evolution of the stresses in granular materials.

Effect of particle stiffness on conductive fluxes
Due to the large difference of the conductive flux depending on the particle concentration, Fig. 9 only depicts the influence of the particle stiffness for high particle loadings. Namely we study (a) φ p = 0.58 and (b) φ p = 0.64 for a Peclet number of 0.01. Thereby, higher dimensionless shear ratesγ indicate softer particles (compare Eq. 16). Note, for lower particle concentrations the conductive flux dramatically decreases, and hence is of no significance for the practical application.
As seen on the left side of the figure, the dimensionless conductive heat flux drops to 10 −5 (low Biot limit) and 10 −9 (high Biot limit) for situations involving very stiff particles, i.e., atγ = 10 −4 . With increasing dimensionless shear rate the dimensionless conductive heat flux increases with a constant slope of two in double logarithmic space. This is observed over a wide range of Biot numbers. Thus, the transferred heat to the environment has no influence on this scaling law. Only at the highest studied shear rate this is not true any more, since the flow approaches an intermediate regime of heat conduction, similar to the one observed for stresses in sheared granular matter. This trend can also be observed for smaller particle volume fractions (data not shown), with the only difference that conductive fluxes are generally much smaller in such a flow. In the intermediate regime the contribution of the conductive heat flux reaches a noticeable value. Most important, this is also the case for intermediate Biot numbers which are relevant for some industrial applications (i.e., Bi ≈ 1). In contrast, the contribution for higher Biot numbers remains in the negligible range (e.g., for Bi = 10, the conductive flux is only 1/1000 of the reference flux). Investigating a higher particle volume fraction (i.e., φ p = 0.64, see Fig. 9 (panel b)), the dimensionless conductive heat flux is not affected by the dimensionless shear rate anymore, and can be assumed to be constant for every shear rate. The results for the higher volume fraction and the dimensionless value are in agreement with Rognon et al. [16], even though their study was limited to one volume fraction. In our present work the flow reached the quasi-static regime where the ratio of contact forces and stiffness, which can be interpreted as a mean particle-particle overlap, does not change when changing the particles stiffness. Thus, the widely accepted fact [11] that for dense quasi-static granular (i) the dimensional contact stress scales linearly with the particles stiffness, and (ii) that the normalized stress is independent of the (dimensionless) shear rate in fixed volume simulations is of crucial importance: since contact stress and contact force are directly proportional, this also means that a typical inter-particle (contact) force scales linearly with the particles stiffness. Hence, the ratio of contact forces and stiffness, which can be interpreted as a mean particle-particle overlap, does not change when changing the particles stiffness in a fixed volume quasistatic flow. Since particle-particle conductive heat flux is affected by the overlapping area between particles, one would also expect that for dense granular flow in the quasi-static regime the conductive heat flux is not affected by the particles stiffness. Exactly this is what we observe in panel b of Fig. 9.
Again, as already seen in Fig. 9 (panel a), the conductive flux near a shear rate ofγ = 10 −1 has an exceptional character also for the very dense regime: a slight increase of the conductive heat flux can be observed, indicating the shift to the intermediate regime. After evaluating the influence of a variety of Biot numbers and particle stiffness on the conductive flux, we will now investigate their effect on the convective fluxes as defined in Eq. (12).

Effect of the Peclet number and shear rate on convective fluxes
A number of research groups already made an attempt to model the convective heat flux in sheared beds. Mohan et al. [13] investigated the convective transport in combination with wet particles including liquid bridge formation, but never studied the influence on the resulting thermal flux. Rognon et al. [16] investigated the scaling of the convective heat flux in the gradient direction, and found that q conv,y ≈ 10 −2 q s τ . The thermal number τ introduced by Rognon et al. scales linearly with the Peclet number defined in our contribution. Specifically, by assuming that the stress scales with k n /d p , we find from a comparison with Rognon et al.'s paper that τ =∼γ Pe. Most important, since Rognon et al. [16] found that the particle stiffness has very little influence on the thermal transport rate for high Peclet numbers, τ and Pe are expected to be directly proportional and not affected byγ . The previous work of Rognon et al. showed a regime map where the convective flux exceeds the conductive flux for τ ≥ ≈ 20. We extend this previous analysis of Rognon et al. [16] to a variety of particle volume fractions, and have summarized our results in Fig. 10. For the moment we focus on a low Biot number in order to be consistent with Rognon et al. [16]. Both particle volume fractions considered  Fig. 11 Influence of the Peclet number and dimensionless shear rate on the dimensionless convective flux for a variety of particle volume fractions aγ = 10 −1 , bγ = 10 −2 in the present work show similar values for the dimensionless convective heat flux, indicating that random particle motion is only mildly influenced by the particle concentration in our simulations. For high particle concentrations (i.e., flow in the quasi-static regime, see Fig. 10b) our results are in good agreement with Rognon et al. [16]. This is true for both the slope with respect to the Peclet number, as well as the range of values. However, Rognon et al. [16] presented results in which an increase of the inertial number (i.e., I ∼γ ) does not necessarily cause an increase in the convective flux. This was reported to be especially the case at high Peclet numbers. We find that this is not always the case for somewhat lower particle concentrations: the dimensionless convective heat flux derived from our simulations indicated a significantly different, i.e., 1.5 times larger heat flux when increasingγ from 10 −4 to 10 −1 . In other words: our results indicate that softer particles clearly lead to higher convective fluxes in dilute flows characterized by a low Peclet number.
We now have a closer look at the dependency of the convective fluxes near the critical particle concentration, i.e., the jamming point. Specifically, we consider the convective heat flux for the same particle stiffness (equalγ ) over a wide range of particle volume fractions. Figure 11a summarizes our findings for comparably stiff particles, specifically, we observe a minimum for the convective flux at φ p ≈ 0.58. We can only speculate on the origin of this minimum, which is clearly a result of the dependency of the particles' random motion on the particle concentration. A tentative explanation is that below the critical volume fraction the velocity in Eq. (9) is the leading term. For φ p>0.58 the internal energy of the particles (m eff c p T vol,avg ) is increasing, simply because more particles fit into a certain volume. However, the convective flux increases by approximately one order of magnitude in case we increase the particle volume fraction by just 5 percent. Hence, our simple explanation cannot fully describe such a strong increase in the thermal transport rate.
We finally note that according to Chialvo et al. [11] the critical particle volume fraction for transition to the quasistatic regime (in the context of granular rheology) is system dependant, and mainly influenced by the inter-particle friction coefficient. Specifically, this previous study identified a minimum critical volume fraction for μ = 1.0 of (φ c=0.581 ), and a maximum for μ = 0.0 at (φ c = 0.636). Thus, when applying the results of our study, one needs to consider the inter-particle friction coefficient used in our present contribution.
Next, we consider the effect of the Biot number on the conductive fluxes. Figure 12 depicts the influence of the Biot number on the conductive flux over a wide range of particle volume fractions. Thereby we focus on the highest and smallest studied Peclet number in order to gain some understanding how the speed of shearing affects the heat fluxes. As can be seen, the dimensionless conductive heat flux is constant for very low Biot numbers, and monotonically decreases with increasing Biot number (there are some fluctuations for Pe = 100, situations involving stiff particles, and comparably dilute particle beds. Conductive fluxes are very small for these systems, and hence not of primary interest). This trend is also seen for other Peclet numbers and dimensionless shear rates (data not shown). Most important, the point at which the the Biot number influences the conductive heat flux shows a minor dependence on the particle volume fraction φ p and shear rateγ . For Pe = 0.01 this point is Bi ≈ 10 −2 , while for Pe = 100 this point is located at Bi ≈ 1. As seen before the scattering of conductive fluxes is inversely proportional to the dimensionless shear rate: higher dimensionless shear rates lead to less scattering (compare Fig. 12a, b). A major dependency of the dimensionless conductive heat flux on the Peclet number is observed, especially for high Bi. According to Eq. (20), at the same dimensionless shear rate and cooling rate (i.e., fixed α), an increase in the particles' thermal conductivity results in a lower Peclet number. However, in such a scenario the Biot number would also decrease inversely proportional to the heat conductivity. Thus, it is instructive to plot the conductive heat flux as a function of Π −1 , since then the particles' heat conductivity does not influence Π . Figure 13 illustrates these data, where we can now easily identify the effect of intra-particulate heat transfer: for a fixed value of Π and a fixed dimensionless shear rate, a decrease in Pe reflects an increase of the particles' heat conductivity, and hence smaller intra-particle temperature gradients.

Impact of the Biot number on conductive fluxes
If we now consider the data forγ = 10 −2 in Fig. 13 (panel a), we observe that the conductive heat flux decreases upon an increase of the Peclet number (i.e., increasing intra-particle temperature gradients). Thus, the surface temperature for the case of Pe = 100 is closer to that of the ambient fluid, leading to a suppressed heat conduction. Surprisingly, this effect is less pronounced for stiff particles, i.e.,γ = 10 −4 . The conductive heat transfer rate for such stiff and rapidly cooled particles is, however, small, and hence of no technological relevance. Considering a more dense particle bed (i.e., the data shown in Fig. 13, panel b), we see similar trends, illustrating the importance to account for intra-particle temperature gradients in situations in which dense particle beds are rapidly cooled or slowly sheared (i.e., Π is small). Interestingly, in particle beds undergoing fast shearing (characterized by large values for Π ), the effect of intra-particle temperature gradients is rather small for φ = 0.58. However, as illustrated in Fig. 13, panel b, for φ = 0.64 and stiff particles, an increase of the Biot number (for Π = const) results in a large drop of the conductive flux. We have not systematically explored the origin of this effect. We speculate that it is due to the drop of the surface temperature in high Binumber flows that greatly affects the conductive flux through the bed. After evaluating the influence of the Biot number, we now look more closely on the influence of the Peclet number in dense packed beds for a fixed Bi number. This is to illustrate the important effect of the speed of shearing on the conductive heat flux, which has been identified in the previous paragraph.

Effect of the Peclet number on conductive flux at several Biot numbers
In Fig. 14 we plot the dimensionless conductive heat flux over various Peclet numbers with different dimensionless shear rates and Biot numbers ranging from 10 −5 to 25. As mentioned before, a higher Peclet number reflects a fast shearing of the particle bed. From Fig. 14 (panel a to d) we see that a change in the Peclet number has very little influence on the conductive heat flux for situations characterized by low Biot numbers (e.g., Bi ≤ 0.1). This is in good agreement with Rognon et al. [16], and simply indicates that cooling caused by the ambient fluid is irrelevant since particles are quickly dispersed. In contrast, the dimensionless heat flux for higher Biot numbers strongly depends on the Peclet number: in some cases (e.g., Fig. 14, panel a, Bi = 25) we observe changes up to three orders of magnitude in the conductive heat flux. In other flow situations (Fig. 14b-d) the enhancement of the conductive heat flux is only seen for relatively high Peclet numbers. In summary, in dense particulate flows that are characterized by a larger Peclet number, a smaller difference of the conductive heat transfer rate is observed upon a change of the Biot number. The physical argument is, that in these quickly sheared beds heat is re-distributed by random particle motion, and hence intra-particle temperature gradients are of secondary importance. As mentioned before, Rognon et al. [16] showed the existence of a general scaling law for the heat flux in sheared beds, and Mohan et al. [13] adapted and extended this to the flux of liquid adhering to the particles. We next make an attempt to establish a regime map that also accounts for the cooling by the ambient fluid.

Regimes of conductive heat transport
As demonstrated by Mohan et al. [13], a regime map for the conductive fluxes can be constructed by appropriate rescaling of the dimensionless shear rate and fluxes. The idea is to perform simulations for a variety of particle concentrations (in our case φ was varied between 0.50 and 0.64), and then hope for a collapse of the data after the re-scaling operation. We also follow such and idea in the present contribution, and use for the re-scaling operation the exponents a = 1/3, b = 7/5 and φ c = 0.613 as shown in the axis labels of Fig. 15. Most important, in addition to the particle concentration, we have also varied the Peclet and Biot number in our study. A collapse of the regimes for this variation in the thermal parameters would indicate that the (normalized) conductive flux is only a function of the flow parameters of the system (i.e.,γ and φ). We now probe if this is truly the case.
First, we investigate the influence of the Peclet number on the inertial and quasi-static regime at fixed Biot numbers. Figure 15 (panel a) reveals the dependency of the re-scaled fluxed on the Peclet number for a slow rate of cooling, i.e., Bi = 10 −4 . Our data indicates a reasonable agreement of the re-scaled data for all studied Peclet numbers (note that we studied Pe = 0.01 in most detail). Since we focus on particle volume fractions with φ p ≥ 0.50, all three flow regimes known from granular rheology can be identified: (i) an interial regime characterized by a quadratic increase of the conduc-tive flux with the shear rate, (ii) a quasi-static regime where fluxes are independent of the shear rate, and (iii) an intermediate regime that is only relevant for very soft particles. Most important, our data is also in good agreement with the previous work of Mohan et al. [13]. Clearly, this is due to the small amount of heat transferred to the ambient fluid: in Fig. 15 (panel a) we considered a very small Biot number, and large values for Π ranging from 10 2 to 10 6 . The latter hints to the fact that the rate of cooling is small compared to convective heat transport, i.e., the redistribution of thermal energy due to particle motion in our simulation domain is fast.
As one might anticipate, a different behaviour will develop at higher Biot numbers, i.e., in case particles are rapidly cooled by the ambient fluid. Specifically, at a Biot number of Bi = 2, and as shown in Fig. 15 (panel b), we observe that the Peclet number has a critical effect on the conductive heat flux. First, considering the data for Pe = 100 (i.e., Π = 50) we observe that conductive heat fluxes are still of comparable magnitude to that observed for Bi = 10 −4 . Also, the three regimes can be identified, all with the same scaling with respect to the shear rate. The physical interpretation is that the redistribution of thermal energy is still fast compared to the cooling rate as indicated by the large value for Π . Thus, particles cannot cool down significantly while they move through the shear flow in our simulations. Consequently, the Biot number (characterizing the intra-particle temperature distribution) has little effect for large Peclet number flows as long Π is large. This is due to the fact that the overall cooling of the particles is small, and hence intra-particle temperature gradients are (relative to the temperature gradient applied to the simulation domain) also small.
Most important, the situation changes as we decrease the Peclet number: for Pe = 1 (i.e., Π = 0.5) we already see a substantial decrease of the conductive heat flux, and with decreasing Pe (and hence Π 1) thermal transport due to conduction is increasingly suppressed. Also, the slope of the curve characterizing the normalized conductive heat flux in the inertial regime appears to gradually increase with decreasing Π . However, this increase is rather small. We do not attempt to speculate on the physical origin of this finding due to the lack of more data that can support our arguments.
The physical interpretation of the observed decrease of the conductive heat flux for intermediate and small Pe numbers is that a layer of "equilibrated" particles forms: these particles have a temperature close to that of the ambient fluid. Consequently, this layer of equilibrated particles chokes off the conductive heat flux through the simulation box, since heat is primarily transferred to the ambient fluid. Clearly, this layer becomes thicker (and the layer in which particles still having an appreciable temperature gradient becomes thinner) when decreasing Π . As already discussed in our theoretical analysis presented in Sect. 3.3, the size of the simulation domain is important in these situations: clearly, the size of the simulation domain determines the relative thickness of the layer containing equilibrated particles, and hence affects the computed (domain-averaged) heat fluxes.

Conclusions
We investigated the heat transfer through a sheared bed of inertial particles that is cooled with a fixed rate by an ambient fluid. Also, we presented the newly developed library ParScale which we coupled to the well established DEM solver LIGGHHTS. This library, which models intra-particle transport processes, was used to resolve radial temperature profiles within every single particle. This enabled us to model the effect of the transferred heat flux (to the ambient fluid) on intra-particle temperature profiles. Hence, we accounted for the competition of heat transferred to the surrounding fluid, and the rate of heat conduction within the particles.
Next, we presented a simple analytical solution for the temperature profile in a sheared bed, and proved the consistency of this solution with our simulation data. Our solution is applicable over a wide range of flow situations, and requires the specification of only one additional system-dependent parameter if applied to simple shear flow.
Furthermore, we showed that our approach to discretize the particle in only one spatial direction (i.e., radially) is in good agreement with data obtained from direct numerical simulations (DNS) that utilize a full three-dimensional discretiztion. Specifically, our simulations indicate that the difference between DNS data, and the data obtained from a one-dimensional discretization is decreasing with rising Reynolds numbers. Also, the ratio of the thermal conductiv-ities of the particle and the surrounding fluid is of a certain significance. Often, intra-particle temperature gradients are of industrial relevance only in high Reynolds number flows (i.e., large particles), and less important in case fluid-particle relative motion is slow (i.e., small Reynolds numbers). This suggests that our approach that relies on ParScale is an affordable, yet accurate way to account for intra-particle temperature gradients.
We then studied thermal transport rates for granular flows characterized by a low Biot number (i.e., slow heat transfer to the surrounding fluid), and found good agreement with literature data that considered particles in vacuum. Also, we extended the analysis by studying a wide range of particle volume fractions. Most important, in situations characterized by a small dimensionless shear rate we found that the convective flux shows a global minimum at the critical particle volume fraction. This finding coincides with the well known maximum in the fluctuations of the (normalized) pressure at the critical state. For softer particles (i.e., flows characterized by a large dimensionless shear rate) this minimum vanishes. The convective heat flux is then nearly unaffected by the particle concentration.
Subsequently, we focused on higher Biot numbers, and studied the influence of the Biot and Peclet number on the conductive heat flux through a sheared particle bed. We found that the conductive heat flux remains almost constant for relevant Biot numbers (i.e., speeds of cooling) at high Peclet numbers. For dense flows only a weak dependency is seen with respect to the particle stiffness, whereas for more dilute flows the conductive heat flux drops by more than one order of magnitude in case the particle stiffness is increased. This clearly hints to different regimes of heat conduction in dense and dilute beds, in analogy to the three stress regimes that are well accepted in granular rheology.
In order to make one step towards a unified model to predict thermal transport through sheared particle beds, we made an attempt to collapse the scaled conductive heat fluxes. We find that for some of the investigated flow situations the conductive heat fluxes can be collapsed into a single curve. This is a strong indication that a unified regime map that characterizes the mode of heat conduction in a cooled bed of particles can be constructed. Specifically, we showed that for high Peclet and small Biot numbers (i.e., large values for the parameter Π ) the map characterizing the conductive heat flux is nearly independent of the Biot number. Thus, our data can be used to make reliable predictions of conductive heat transfer rates for these situations. In contrast, for small Peclet and high Biot numbers still challenges remain to close our model. Clearly, in order to create a rigorous regime map that allows a precise calculation of conductive, convective and transferred heat transfer rates, a deeper analysis of the data presented in this study is necessary. Most important, such an analysis has to account for the (normalized) size of the simulation domain. This is due to the nature of particle cooling in a granular shear flow, which leads to a non-uniform temperature gradient. We expect that the newly introduced parameter Π , as well as the newly derived expression for the local temperature gradient (see Sect. 3.3), will be inspiration for future work that may account for these effects.
The work of Morris et al. [14] discussed the effect of stiffness on the results, and provides a stiffness correction. Similarly, our results could be used to propose such a stiffness correction. Indeed, and as outlined in Morris et al. [14], developing such a correction for dense granular flows remains a delicate task.
Our study certainly has some limitations: especially in dense flows of soft particles, our approach to model the transferred heat flux should be extended by accounting for the surface area of the particle that is covered due to particle deformation during a collision. For example, this could be done using well established expressions for the particleparticle contact area, or (in more extreme situations with respect to particle concentration and softness) once could adopt a Monte Carlo area integration method. Such a method is already implemented in LIGGGHTS. Another important assumption of the present analysis was that the ambient fluid is kept at a fixed temperature. This assumption is valid as long as the heat transferred to the fluid is quickly removed, e.g., convection in the ambient fluid is sufficiently strong. Specifically, this means that the spatial fluid temperature variation is small over a characteristic distance (i.e., a few particle diameters), which is the case for high fluid Peclet numbers, and depends on the particle concentration (Munichi et al. [38]). This assumption of fast convective transport holds for most applications involving comparably large (i.e., mmsized) particles suspended in a gas. However, future thoughts may extend our present work to applications that are characterized by extreme temperature gradients, i.e., systems involving smaller particles and low fluid Peclet numbers (i.e., slow fluid-particle relative motion). Also, since we fixed the ambient fluid temperature and the heat transfer coefficient, both quantities must be provided from a supplementary fluid model (e.g., that used in Kuipers et al. [24]). This exceeded the current study, however, can be implemented in a straight forward manner. Future considerations could then establish an even more rigorous understanding of thermal transport in flowing granular matter and the ambient fluid encountered in a number of engineering applications.