1 Introduction
The design and operation of various industrial processes is highly dependent upon the transport of momentum and thermal energy within a gas–solids flow. In many systems, domain walls or immersed surfaces are utilized as the primary energy source to heat a particle-laden mixture (Syamlal & Gidaspow Reference Syamlal and Gidaspow1985; Kuipers, Prins & Van Swaaij Reference Kuipers, Prins and Van Swaaij1992; Nijemeisland & Dixon Reference Nijemeisland and Dixon2004; Chaudhuri, Muzzio & Tomassone Reference Chaudhuri, Muzzio and Tomassone2006; Patil et al. Reference Patil, Smit, van Sint Annaland and Kuipers2006; Cong et al. Reference Cong, He, Chen, Ding and Wen2007; Shi, Vargas & McCarthy Reference Shi, Vargas and McCarthy2008; Martinek & Ma Reference Martinek and Ma2014; Morris et al. Reference Morris, Pannala, Ma and Hrenya2015, Reference Morris, Ma, Pannala and Hrenya2016; Yohannes et al. Reference Yohannes, Emady, Anderson, Paredes, Javed, Borghard, Muzzio, Glasser and Cuitiño2016; Ansart et al. Reference Ansart, García-Triñanes, Boissière, Benoit, Seville and Simonin2017; Bongo Njeng et al. Reference Bongo Njeng, Vitu, Clausse, Dirion and Debacq2018). Under such conditions, the drag and heat transfer occurring near a wall will be of primary significance. Despite prevalent use of such flows in industry, fundamental explorations on wall-to-particle heat transfer or near-wall particle drag have not been emphasized in the literature. While a variety of drag and convection correlations have been reported for unbounded gas–solids flows (no walls) (Ranz & Marshall Reference Ranz and Marshall1952; Whitaker Reference Whitaker1972; Haider & Levenspiel Reference Haider and Levenspiel1989; Feng & Michaelides Reference Feng and Michaelides2000; Richter & Nikrityuk Reference Richter and Nikrityuk2012), they inherently do not account for boundary effects. By and large, these unbounded correlations are applied as is to the near-wall region, where their validity is expected to deteriorate. On many occasions, direct numerical simulation (DNS) has been employed to probe the drag and heat transfer occurring within an unbounded gas–solids system (Feng & Michaelides Reference Feng and Michaelides2000, Reference Feng and Michaelides2008, Reference Feng and Michaelides2009; Dan & Wachs Reference Dan and Wachs2010; Deen et al. Reference Deen, Kriebitzsch, van der Hoef and Kuipers2012, Reference Deen, Peters, Padding and Kuipers2014; Richter & Nikrityuk Reference Richter and Nikrityuk2012; Tavassoli et al. Reference Tavassoli, Kriebitzsch, van der Hoef, Peters and Kuipers2013; Tavassoli, Peters & Kuipers Reference Tavassoli, Peters and Kuipers2015; Kruggel-Emden et al. Reference Kruggel-Emden, Kravets, Suryanarayana and Jasevicius2016; Kravets & Kruggel-Emden Reference Kravets and Kruggel-Emden2017). However, works to date which account for boundary effects are far less inclusive (Nijemeisland & Dixon Reference Nijemeisland and Dixon2004; Radl, Municchi & Goniva Reference Radl, Municchi and Goniva2016; Municchi & Radl Reference Municchi and Radl2017).
The heat transfer occurring between a particle and a wall is comprised of convective, conductive and radiative mechanisms. For the case of moderate system temperatures 
$(T<700~\text{K})$, radiation is often neglected since it is not a significant contribution to the overall heat transfer (Wen & Chang Reference Wen and Chang1967; Flamant & Menigault Reference Flamant and Menigault1987). Under these circumstances, the relevant heat transfer mechanisms may be simplified to convection and conduction only. Typically, correlations for unbounded systems (Ranz & Marshall Reference Ranz and Marshall1952; Whitaker Reference Whitaker1972; Feng & Michaelides Reference Feng and Michaelides2000; Richter & Nikrityuk Reference Richter and Nikrityuk2012) are utilized without modification to approximate fluid–particle heat transfer in the near-wall region. To account for wall–particle heat transfer, which has been believed to be mostly conductive, particle-scale theories have been developed and employed. Specifically, the conduction occurring between a particle and wall is made up of two contributions: (i) direct conduction through the particle–wall contact area (Cooper, Mikic & Yovanovich Reference Cooper, Mikic and Yovanovich1969; Batchelor & O’Brien Reference Batchelor and O’Brien1977) and (ii) indirect conduction between a particle and wall separated by a thin layer of fluid (Rong & Horio Reference Rong and Horio1999). In many practical cases, indirect conduction tends to dominate over direct conduction – i.e. when the ratio of thermal resistances associated with direct and indirect conduction is much greater than unity, or 
$\unicode[STIX]{x1D6FD}=R_{p}k_{g}{\hat{h}}_{PFW}(0)/2k_{p}R_{c}\gg 1$ (Lattanzi & Hrenya Reference Lattanzi and Hrenya2017), where 
$R_{p}$ is the particle radius, 
$k_{g}$ is the gas thermal conductivity, 
${\hat{h}}_{PFW}(0)$ is the solution to the indirect conduction integral at a particle–wall separation distance of zero, 
$R_{c}$ is the radius of contact between the particle and wall and 
$k_{p}$ is the particle thermal conductivity. While particle-scale theories for indirect conduction (Delvosalle & Vanderschuren Reference Delvosalle and Vanderschuren1985; Cheng, Yu & Zulli Reference Cheng, Yu and Zulli1999; Rong & Horio Reference Rong and Horio1999; Vargas & McCarthy Reference Vargas and McCarthy2002) have been applied to a wide variety of systems (Zhou, Flamant & Gauthier Reference Zhou, Flamant and Gauthier2004; Di Maio, Di Renzo & Trevisan Reference Di Maio, Di Renzo and Trevisan2009; Zhou, Yu & Zulli Reference Zhou, Yu and Zulli2009; Morris et al. Reference Morris, Pannala, Ma and Hrenya2015, Reference Morris, Ma, Pannala and Hrenya2016; Lattanzi & Hrenya Reference Lattanzi and Hrenya2016), the theories themselves have only been validated for static systems (Mishra et al. Reference Mishra, Lattanzi, LaMarche, Morris and Hrenya2019). Most commonly, indirect conduction theory assumes that each particle is surrounded by a static-fluid lens 
$(R_{Lens})$, as denoted by the dashed line in figure 1. When the fluid lens overlaps with the wall, one-dimensional conduction is assumed to occur through the fluid lens. Therefore, the fluid lens thickness is the key length scale that establishes distances over which particle–wall conduction will occur. For dynamic systems, indirect conduction theory has been shown to be most sensitive to the fluid lens thickness parameter, which is traditionally set according to the particle size 
$(R_{Lens}\propto R_{p})$ (Lattanzi & Hrenya Reference Lattanzi and Hrenya2017). The state of the art for modelling heat transfer to a particle in the near-wall region still involves the use of unbounded heat transfer correlations in conjunction with particle-scale theories for indirect conduction (Morris et al. Reference Morris, Pannala, Ma and Hrenya2015). For gas–solids flows at moderate temperatures (dominated by convection and indirect conduction), further work is required to assess the accuracy of these approximations in the near-wall region.
Figure 1. An illustration of the static-fluid lens (dashed line) theory employed by indirect conduction theory (left particle). The heat transfer occurring between a particle and wall when the lens overlaps with a wall – i.e. the particle–wall separation distance 
$(\unicode[STIX]{x1D6FF})$ is less than the fluid lens thickness 
$(R_{Lens}-R_{p})$ (right particle).
In the present work, uniform flow of a fluid past a hot plate and a static, cold particle was simulated by a hybrid lattice Boltzmann–random walk particle tracking (LBM-RWPT) DNS code (Wang et al. Reference Wang, Koch, Yin and Cohen2009; Metzger, Rahli & Yin Reference Metzger, Rahli and Yin2013; Lattanzi, Yin & Hrenya Reference Lattanzi, Yin and Hrenya2019a,Reference Lattanzi, Yin and Hrenyab) to examine the heat and momentum transfer to a spherical particle in the near-wall region. The presence of a hot wall in this work allowed boundary effects on particle drag force and wall-to-particle heat transfer to be quantified. Particle drag forces and heat rates obtained from LBM-RWPT are compared to predictions made from unbounded correlations (Ranz & Marshall Reference Ranz and Marshall1952; Haider & Levenspiel Reference Haider and Levenspiel1989) coupled with indirect conduction (Rong & Horio Reference Rong and Horio1999) closures commonly employed within the discrete element method (DEM) framework. Use of an unbounded drag correlation with a free-stream fluid velocity was found to over-predict the drag force in the near-wall region since the effect of the slow moving fluid adjacent to the wall is not accounted for. To capture effects resulting from the particle’s local environment, a local Reynolds number was utilized to develop a local, unbounded drag correlation. The drag correlation with local Reynolds number was observed to capture the drag force in the near-wall region as well as in the limit of an unbounded system (i.e. large particle–wall separation distances; 
$\unicode[STIX]{x1D6FF}$). Similarly, unbounded convection correlations were found to under-predict the heat transfer to a particle in the near-wall region since they do not account for the presence of the wall. By contrast, the combination of unbounded convection and indirect conduction considerably improves agreement with LBM-RWPT. Specifically, indirect conduction theory is observed to capture the main physics associated with heat transfer enhancement in the near-wall region. However, further heat transfer enhancement is observed in LBM-RWPT at particle–wall separation distances 
$(\unicode[STIX]{x1D6FF})$ that are not captured by indirect conduction theory. Namely, indirect conduction theory sets the fluid lens thickness according to geometric arguments based upon the particle size 
$(R_{Lens}=1.4R_{p})$ (Vargas & McCarthy Reference Vargas and McCarthy2002; Morris et al. Reference Morris, Pannala, Ma and Hrenya2015, Reference Morris, Ma, Pannala and Hrenya2016; Lattanzi & Hrenya Reference Lattanzi and Hrenya2017) and predicts near-wall heat transfer will only occur when the fluid lens intersects the wall 
$(\unicode[STIX]{x1D6FF}\leqslant 0.4R_{p})$. However, setting the fluid lens thickness in this manner neglects the thermal length scale associated with the fluid near the wall (boundary layer thickness). An excess wall Nusselt number was developed to account for such near-wall heat transfer enhancement. Superposition of the local, unbounded convection correlation and excess wall correlation is observed to accurately capture the DNS data while the excess wall correlation asymptotically decays to zero in the limit of large particle–wall separation.
2 Background: indirect conduction theory
To account for the indirect conduction occurring between a particle and wall, we employ the theory proposed by Rong and Horio (Rong & Horio Reference Rong and Horio1999; Morris et al. Reference Morris, Pannala, Ma and Hrenya2015). In this theory, particles are assumed to be surrounded by a static-fluid lens (dashed line in figure 1). When the lens overlaps with the wall, one-dimensional conduction through the fluid lens is assumed to occur between the particle and wall. Motivation for describing the fluid lens as ‘static’ is guided by the effect of no-slip boundary conditions on the particle and wall. As the separation distance 
$(\unicode[STIX]{x1D6FF})$ between the particle and wall becomes small, the fluid velocities between the particle and wall are dramatically reduced from the free-stream velocity. The rate of heat transfer due to indirect conduction is found by integrating Fourier’s law over the area of overlap between the fluid lens and wall (Morris et al. Reference Morris, Pannala, Ma and Hrenya2015)
$$\begin{eqnarray}\left.\begin{array}{@{}rcl@{}}\dot{Q}_{PFW}\, & \equiv \, & h_{PFW}[T_{w}-T_{p}]=\displaystyle \int _{r_{in}}^{r_{out}}{\displaystyle \frac{2\unicode[STIX]{x03C0}k_{g}r}{\text{Max}(l,s)}}(T_{w}-T_{p})\,\text{d}r,\\[15.0pt] r_{in}\, & =\, & \left\{\begin{array}{@{}ll@{}}r_{s}=\sqrt{R_{p}^{2}-(s-R_{p}-\unicode[STIX]{x1D6FF})^{2}} & \unicode[STIX]{x1D6FF}\leqslant s,\\ 0 & \unicode[STIX]{x1D6FF}>s,\end{array}\right.\\[17.0pt] r_{out}\, & =\, & \left\{\begin{array}{@{}ll@{}}\sqrt{R_{p}^{2}-(R_{p}+\unicode[STIX]{x1D6FF})^{2}} & \hspace{33.5001pt}\unicode[STIX]{x1D6FF}>\sqrt{R_{Lens}^{2}-R_{p}^{2}}-R_{p},\\ R_{p} & \hspace{33.5001pt}\unicode[STIX]{x1D6FF}\leqslant \sqrt{R_{Lens}^{2}-R_{p}^{2}}-R_{p},\end{array}\right.\end{array}\right\}\end{eqnarray}$$ where 
$\dot{Q}_{PFW}$ is the rate of heat transfer due to indirect conduction between the wall and particle, 
$h_{PFW}$ is the particle–fluid–wall heat transfer coefficient, 
$T_{w}$ is the wall temperature, 
$T_{p}$ is the particle temperature, 
$r$ is the radial position of the fluid lens overlap, 
$l(r)$ is the conduction distance at a radial position of 
$r$, 
$s$ is the minimum conduction distance, 
$\unicode[STIX]{x1D6FF}$ is the particle–wall separation distance and 
$R_{Lens}$ is the fluid lens radius. To evaluate the integral in (2.1), a fluid lens radius 
$(R_{Lens})$ and minimum conduction distance 
$(s)$ must be specified. An upper bound for 
$R_{Lens}$ is generally determined from geometric arguments and is given by 
$R_{Lens}=\sqrt{2}R_{p}\approx 1.41R_{p}$. Namely, the maximum fluid lens radius is set such that the upper bound of integration in (2.1) 
$(r_{out})$ does not exceed the particle radius at the point of solid body contact 
$(\unicode[STIX]{x1D6FF}=0)$ – i.e. the conduction distance 
$(l)$ remains well defined. The fluid lens radius utilized in this work matches that commonly employed in other works 
$(R_{Lens}=1.4R_{p})$ (Vargas & McCarthy Reference Vargas and McCarthy2002; Morris et al. Reference Morris, Pannala, Ma and Hrenya2015, Reference Morris, Ma, Pannala and Hrenya2016; Lattanzi & Hrenya Reference Lattanzi and Hrenya2016, Reference Lattanzi and Hrenya2017). The minimum conduction distance 
$(s)$ in (2.1) acts as a lower bound for the conduction distance 
$(l)$. The minimum conduction distance can be physically interpreted as corresponding to either the size of surface asperities (roughness) or the mean free path of the gas (perfectly smooth). For the former case, large-scale asperities on the surface of a particle or wall will result in finite separation distances even at contact. For the latter case, as the particle and wall tend to solid body contact 
$(\unicode[STIX]{x1D6FF}=0)$, the conduction distance 
$(l(r))$ becomes small with respect to the mean free path of the gas and rarefaction effects become non-negligible. By setting the minimum conduction distance to the mean free path of the gas (air at standard temperature and pressure; 
$2.75\times 10^{-8}~\text{m}$), the integration in (2.1) avoids conduction lengths where rarefaction effects are significant. Here, the particle and wall will be assumed to be perfectly smooth and the range of separation distances 
$(\unicode[STIX]{x1D6FF})$ considered will be significantly larger than the mean free path of the gas. Therefore, the lower bound of integration in (2.1) 
$(r_{in})$ will always be 0 in the present work (i.e. particle–wall contact will not be considered). The integral in (2.1) may be non-dimensionalized and directly evaluated (Lattanzi & Hrenya Reference Lattanzi and Hrenya2017)
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle {\hat{h}}_{PFW}=2\unicode[STIX]{x03C0}\left[(1+\hat{\unicode[STIX]{x1D6FF}})\ln \left[\left|\frac{\sqrt{1-\hat{r}_{out}^{2}}-(1+\hat{\unicode[STIX]{x1D6FF}})}{\hat{\unicode[STIX]{x1D6FF}}}\right|\right]+\sqrt{1-\hat{r}_{out}^{2}}-1\right],\\ \displaystyle \hat{r}_{out}=\left\{\begin{array}{@{}ll@{}}\sqrt{{\hat{C}}^{2}-(1+\hat{\unicode[STIX]{x1D6FF}})^{2}} & \hat{\unicode[STIX]{x1D6FF}}>\sqrt{{\hat{C}}^{2}-1}-1,\\ 1 & \hat{\unicode[STIX]{x1D6FF}}\leqslant \sqrt{{\hat{C}}^{2}-1}-1,\end{array}\right.\end{array}\right\} & & \displaystyle\end{eqnarray}$$ where 
$\hat{~}$ denotes normalization by the particle radius, and 
${\hat{C}}=R_{Lens}/R_{p}=1.4$ is the fluid lens proportionality constant. The rate of heat transfer at a given dimensionless separation distance 
$(\hat{\unicode[STIX]{x1D6FF}}=\unicode[STIX]{x1D6FF}/R_{p})$ then becomes 
$\dot{Q}_{PFW}=k_{g}R_{p}{\hat{h}}_{PFW}(\hat{\unicode[STIX]{x1D6FF}})[T_{w}-T_{p}]$.
3 Numerical techniques
3.1 Lattice Boltzmann method
The DNS framework is a hybrid scheme based on two coupled methods. The first is the lattice Boltzmann method (LBM), which is utilized to resolve the fluid phase – i.e. solve the Navier–Stokes (NS) equations. The LBM scheme employed here matches that developed by Ladd and co-workers (Ladd Reference Ladd1994a,Reference Laddb; Ladd & Verberg Reference Ladd and Verberg2001). Due to a foundation in statistical mechanics, LBM discretizes the continuous Boltzmann equation rather than the NS equations. Since the Boltzmann equation governs the evolution of the molecular velocity distribution, LBM utilizes discrete velocity distributions (population densities) as opposed to the hydrodynamic variables. The discrete velocity distributions are updated in this work according to the classic streaming and collision process
$$\begin{eqnarray}n_{i}(\boldsymbol{r}+\boldsymbol{c}_{i}\unicode[STIX]{x0394}t,\unicode[STIX]{x0394}t)\equiv n_{i}^{\ast }(\boldsymbol{r},t)=n_{i}(\boldsymbol{r},t)+\unicode[STIX]{x1D6FA}_{i}(\boldsymbol{n}(\boldsymbol{r},t)),\end{eqnarray}$$ where 
$n_{i}$ is the discrete velocity distribution associated with molecular velocity 
$\boldsymbol{c}_{i}$, 
$\unicode[STIX]{x0394}t$ is the LBM time step, 
$\unicode[STIX]{x1D6FA}_{i}$ is the collision operator (function of all the velocity distributions at a node 
$\boldsymbol{n}(\boldsymbol{r},t)$) and 
$n_{i}^{\ast }$ is the post-collision distribution function (expanded about the local equilibrium; 
$\boldsymbol{n}^{eq}$). The hydrodynamic quantities are given by the moments of the velocity distribution functions
$$\begin{eqnarray}\unicode[STIX]{x1D70C}=\mathop{\sum }_{i}n_{i},\quad \boldsymbol{j}\equiv \unicode[STIX]{x1D70C}\boldsymbol{u}=\mathop{\sum }_{i}n_{i}\boldsymbol{c}_{i},\quad \unicode[STIX]{x1D72B}=\mathop{\sum }_{i}n_{i}\boldsymbol{c}_{i}\boldsymbol{c}_{i},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D70C}$ is the density, 
$\boldsymbol{j}$ is the momentum, 
$\boldsymbol{u}$ is the macroscopic velocity and 
$\unicode[STIX]{x1D72B}$ is the fluid stress tensor. The update scheme given in (3.1) may ultimately be shown to recover the incompressible Navier–Stokes equations in the low Mach limit with the following closures for the shear 
$(\unicode[STIX]{x1D702})$ and bulk 
$(\unicode[STIX]{x1D702}_{b})$ viscosities (Ladd & Verberg Reference Ladd and Verberg2001)
$$\begin{eqnarray}\unicode[STIX]{x1D702}=-\unicode[STIX]{x1D70C}c_{s}^{2}\left[\frac{1}{\unicode[STIX]{x1D706}}+\frac{1}{2}\right],\quad \unicode[STIX]{x1D702}_{b}=-\frac{2\unicode[STIX]{x1D70C}c_{s}^{2}}{3}\left[\frac{1}{\unicode[STIX]{x1D706}_{b}}+\frac{1}{2}\right],\end{eqnarray}$$ where 
$c_{s}^{2}=1/3$ is the square of the speed of sound, and 
$\unicode[STIX]{x1D706}$ and 
$\unicode[STIX]{x1D706}_{b}$ are eigenvalues of the collision matrix. Here, 
$\unicode[STIX]{x1D706}$ corresponds to the relaxation of the off-diagonal portion of the non-equilibrium stress tensor while 
$\unicode[STIX]{x1D706}_{b}$ corresponds to the relaxation of the diagonal portion of the non-equilibrium stress tensor. Coupling between the fluid phase and solid particles is completed by imposing a no-slip boundary condition at the particle surface. The net force and torque applied to a particle by the fluid is given by surface integration of the interphase momentum transfer (resulting from the no-slip boundary condition). The force and torque due to interphase momentum transfer and particle collisions may be utilized to find the new particle velocity and position (solid body mechanics). However, the particle in this work is held static and at finite particle–wall separation distances (no particle collisions occur).
3.2 Random walk particle tracking
The second method within the DNS framework is random walk particle tracking (RWPT). RWPT is employed here to solve the advection–diffusion equation for thermal energy (Gardiner Reference Gardiner1986; Salamon, Fernandez-Garcia & Gomez-Hernandez Reference Salamon, Fernandez-Garcia and Gomez-Hernandez2006; Wang et al. Reference Wang, Koch, Yin and Cohen2009; Metzger et al. Reference Metzger, Rahli and Yin2013; Lattanzi et al. Reference Lattanzi, Yin and Hrenya2019a,Reference Lattanzi, Yin and Hrenyab),
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}T)=\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}(T),\end{eqnarray}$$ where 
$T$ is the thermal temperature and 
$\unicode[STIX]{x1D6FC}$ is the thermal diffusivity. Similar to LBM, RWPT does not directly involve the continuum equation ((3.4) for RWPT), but instead RWPT monitors the positions of many tracers as they undergo displacement. The movement of each tracer depends upon the local velocity field obtained via LBM as well as random fluctuations. An explicit time integration scheme is utilized within the present work to update the position of each tracer
$$\begin{eqnarray}\boldsymbol{r}(t+\unicode[STIX]{x0394}t)=\boldsymbol{r}(t)+\boldsymbol{u}(\boldsymbol{r},t)\unicode[STIX]{x0394}t+\unicode[STIX]{x1D743}(t)\sqrt{2\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}t},\end{eqnarray}$$ where 
$\boldsymbol{r}$ is the position of a tracer, 
$\boldsymbol{u}$ is the velocity at the tracer position before the step (found via trilinear interpolation of the LBM velocity field), 
$\unicode[STIX]{x1D743}$ is a random vector whose entries are sampled from a Gaussian distribution with zero mean and unit variance, 
$\unicode[STIX]{x1D6FC}$ is the thermal diffusivity of the gas, 
$\unicode[STIX]{x0394}t$ is the random walk time step. The thermal temperature in RWPT is proportional to the local tracer concentration. In the present work, we impose a temperature gradient 
$(\unicode[STIX]{x0394}T=T_{1}-T_{0})$ by utilizing two tracer types. Tracers labelled as type ‘1’ correspond to the higher temperature 
$(T_{1})$ while tracers labelled as type ‘0’ correspond to the lower temperature 
$(T_{0})$. Alternatively, the temperature field may be resolved by a single tracer type; however, this would require dynamically re-allocating tracer arrays as the tracers are generated and eliminated at respective constant-temperature boundaries. By utilizing two tracer types, the tracer count remains constant and we only require a conversion of tracer type to recover the Dirichlet boundary condition. The local temperature and dimensionless temperature are given as
$$\begin{eqnarray}\displaystyle & \displaystyle T(\boldsymbol{r},t)=T_{1}\frac{C_{1}(\boldsymbol{r},t)}{C_{t}}+T_{0}\frac{C_{0}(\boldsymbol{r},t)}{C_{t}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}(\boldsymbol{r},t)\equiv \frac{T(\boldsymbol{r},t)-T_{0}}{T_{1}-T_{0}}=\frac{C_{1}(\boldsymbol{r},t)}{C_{t}}, & \displaystyle\end{eqnarray}$$ where 
$C_{1}$ is the concentration of type 1 tracers, 
$C_{0}$ is the concentration of type 
$0$ tracers, 
$C_{t}=C_{1}+C_{0}$ is the total concentration of tracers, a constant throughout the domain, and 
$\unicode[STIX]{x1D703}$ is the dimensionless temperature. For 
$\unicode[STIX]{x1D703}=0$, the system is at a temperature of 
$T_{0}$ and for 
$\unicode[STIX]{x1D703}=1$ the system is at a temperature of 
$T_{1}$. The minimum temperature (
$T_{0}$) and maximum temperature (
$T_{1}$) here will correspond to the particle temperature and wall temperature, respectively (discussed in § 4). For the case of mass transfer, additional tracer types may be introduced to track the concentrations of multiple species, thereby solving a system of advection–diffusion equations.
The RWPT method presented above is often applied to advection dominated transport, since it does not suffer from numerical dispersion. It can be used to simulate both heat transfer and mass transfer that do not actively affect fluid flow. Our method has been applied to both freely evolving multiphase flows (Metzger et al. Reference Metzger, Rahli and Yin2013) and static multiphase systems with conjugate heat transfer (Lattanzi et al. Reference Lattanzi, Yin and Hrenya2019b). While the flow and heat transfer presented in this study may be readily simulated with classic computational methods that employ immersed boundaries or body-fit grids, future work that examines moderately dense particle packings or freely evolving particle suspensions with inter- and intraparticle temperature gradients are expected to benefit from the simple manner in which RWPT addresses these challenges.
4 Systems simulated
Uniform flow past a hot wall and a static, cold particle was considered; see figure 2. Due to the presence of the hot wall, the steady-state fluid flow will be characterized by the development of a hydrodynamic and thermal boundary layer near the bottom plate. The centre of the particle was located 5 particle diameters 
$(D_{p})$ away from the leading edge of the plate 
$(L=5D_{p})$ in all simulations, while the particle–wall separation distance 
$(\unicode[STIX]{x1D6FF})$ and the particle Reynolds number 
$(Re_{Part}\equiv |U_{f}-U_{s}|D_{p}/\unicode[STIX]{x1D708}=U_{\infty }D_{p}/\unicode[STIX]{x1D708})$ were varied. A range of 
$Re_{Part}\in [1~10]$ was selected since it is representative of values present in applications concerned with wall-to-particle heat transfer (Morris et al. Reference Morris, Ma, Pannala and Hrenya2016; Yohannes et al. Reference Yohannes, Emady, Anderson, Paredes, Javed, Borghard, Muzzio, Glasser and Cuitiño2016; Ansart et al. Reference Ansart, García-Triñanes, Boissière, Benoit, Seville and Simonin2017) but also allowed for a straightforward assessment of the drag and Nusselt number corrections that are necessary in the near-wall region. The range for 
$\unicode[STIX]{x1D6FF}$ was chosen such that the particle resides completely within the boundary layers as well as completely outside the boundary layers and is given by 
$\unicode[STIX]{x1D6FF}/R_{p}\in [0.07~12]$. Since the distance from the leading edge 
$(L)$ was fixed, the resulting plate Reynolds number 
$(Re_{Plate}\equiv U_{\infty }L/\unicode[STIX]{x1D708}=5Re_{Part};Re_{Plate}\in [5~50])$ will lie in the intermediate regime and the flow will be laminar 
$(Re_{Plate}<O(10^{6}))$ (White Reference White2005). The particle diameter and fluid Prandtl number 
$(Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D6FC})$ were fixed and set to 
$600~\unicode[STIX]{x03BC}\text{m}$ and 0.7, respectively. The particle diameter was resolved by 15 LBM nodes 
$(D_{p}/\unicode[STIX]{x0394}x_{LBM}=15)$ and a tracer concentration of 1.0 
$(C_{t}=1.0)$ was used in all simulations. The selection of these resolutions will be discussed below in the Grid Convergence section (§ 5). A complete overview of the simulation conditions is given in tables 1–2 while the fluid and particle properties are contained within table 3.
Figure 2. The geometry and boundary conditions utilized to simulate uniform flow past a hot plate (bottom wall in red) and a static, cold particle (blue sphere). The particle–wall separation distance 
$(\unicode[STIX]{x1D6FF})$ is the distance between the bottom of the particle and the wall (varied) while 
$L$ is the distance from the leading edge to the centre of the particle (fixed).
Table 1. LBM-RWPT geometry and mesh.
Table 2. Simulation conditions.
To impose the required boundary conditions given in figure 2, a variety of methods were employed. The hydrodynamic boundary conditions were imposed in the LBM framework. Namely, the no-slip and uniform inflow boundary conditions were achieved via the bounce-back method (Ladd & Verberg Reference Ladd and Verberg2001). The free-slip and outflow boundary conditions were completed by way of the anti-bounce-back method (Jansen & Krafczyk Reference Jansen and Krafczyk2011) and extrapolation (Yang Reference Yang2013; Kruggel-Emden et al. Reference Kruggel-Emden, Kravets, Suryanarayana and Jasevicius2016), respectively. The thermal boundary conditions were imposed in the RWPT framework. Specifically, the constant-temperature boundary at the inflow 
$(\unicode[STIX]{x1D703}=0.2)$ and bottom wall 
$(\unicode[STIX]{x1D703}=1)$ was achieved by a two-step process. First, all tracers that cross the boundary are specularly reflected back into the domain. Second, a number is sampled from 
${\mathcal{U}}(0,1)$. If the sampled value is less than or equal to 
$\unicode[STIX]{x1D703}$, the reflected tracer type is set to 1; else, the tracer type is set to 0. The inlet fluid temperature boundary condition 
$(\unicode[STIX]{x1D703}=0.2)$ is chosen such that it is less than the wall temperature 
$(\unicode[STIX]{x1D703}=1)$ but greater than the particle temperature 
$(\unicode[STIX]{x1D703}=0)$. By setting the inlet temperature boundary condition in this manner a thermal gradient between the particle and fluid will be sustained at large particle–wall separation distance 
$(\unicode[STIX]{x1D6FF})$ and the particle heat transfer will approach the Nusselt correlation for unbounded spheres. By contrast, as the particle approaches the wall 
$(\unicode[STIX]{x1D6FF}\rightarrow 0)$, the inlet temperature is of less significance since the fluid near the wall will be dictated by the wall temperature. Therefore, the effect of the flow on wall–particle heat transfer can be directly evaluated. The constant particle temperature 
$(\unicode[STIX]{x1D703}=0)$ is achieved by setting all tracers that enter the particle to type 0. The adiabatic boundary is imposed by specularly reflecting tracers back into the domain (no alteration of type). The thermal outflow boundary is achieved by a semi-reflecting barrier (Lattanzi et al. Reference Lattanzi, Yin and Hrenya2019a). If a tracer reaches the outflow plane, the probability of being specularly reflected back into the domain 
$(P^{\ast })$ is calculated as in Lattanzi et al. (Reference Lattanzi, Yin and Hrenya2019a). A number is then sampled from 
${\mathcal{U}}(0,1)$. If the value is less than 
$P^{\ast }$, the tracer is specularly reflected back into the domain; otherwise, the tracer is re-seeded at the inflow plane and its type is set according to the temperature boundary condition at the inflow plane 
$(\unicode[STIX]{x1D703}=0.2)$.
Table 3. Hydrodynamic and thermal properties.
Since the particle–wall separation distance 
$(\unicode[STIX]{x1D6FF})$ will become small in the present work, some further comments on the bottom wall boundary condition (constant temperature) and its interaction with the particle are in order. As discussed in Lattanzi et al. (Reference Lattanzi, Yin and Hrenya2019a,Reference Lattanzi, Yin and Hrenyab), the impenetrable boundary is valid if the velocity in the direction normal to the wall tends to zero, which will occur for a no-slip boundary (i.e. the bottom wall here). The basis for the specular reflection treatment is that it recovers the diffusive heat flux emanating from the wall while confining operations to only tracers that cross the boundary during a time step. A subtle, but key, distinction must be made about the difference in which stochastic and discretization methods quantify fluxes. Discretization methods commonly employ low-order polynomial approximation of the continuum equation (3.4) on a grid. Thus, for the fluxes to be well approximated, gradients in the solution variable must be sufficiently resolved by a fine grid. By contrast, the flux in a stochastic method (random walk) is obtained from many tracer trajectories (
$\unicode[STIX]{x0394}\boldsymbol{r}$ in (3.5)) that are not confined to any grid. For the case of a particle and wall separated by a small distance, properly capturing the heat flux with a discretization method would require that the temperature field be resolved between the particle and wall to a high degree of accuracy. For the same system, a stochastic method would require that the trajectories and number of tracers emanating from the wall be statistically meaningful. The wall heat flux and its interaction with the particle was strictly enforced in a sequential manner here. Specifically, a tracer identified as crossing the bottom wall was treated as noted in the previous paragraph and then subsequently checked for entry into the particle. This treatment ensures that, even for small separation distances, the number of tracers crossing the particle surface (flux) is properly counted. Additionally, as noted in Lattanzi et al. (Reference Lattanzi, Yin and Hrenya2019b), the random walk time step may be set independent of the LBM time step but should not exceed it. The random walk time step allows for control of the tracer displacement (
$\unicode[STIX]{x0394}\boldsymbol{r}$). For small separation distances, the tracer displacement between the particle and wall will be strongly governed by diffusion. The random walk time step was set to ensure that the diffusive displacement remained small with respect to the particle–wall separation distance 
$\unicode[STIX]{x1D6FF}/\sqrt{2\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}t}=4.0$ – i.e. 99.98 % of tracers will sample a diffusive displacement that is less than 
$\unicode[STIX]{x1D6FF}$ and 68 % will sample a diffusive displacement that is less that 
$\unicode[STIX]{x1D6FF}/4$. Thus, the particle heat flux is resolved with many tracer trajectories but without the requirement that many tracers reside between the particle and wall; as will be shown in § 5, the addition of tracers acts to suppress the random fluctuations but has little effect on the mean flux.
5 Grid convergence
To gauge the effect of resolution on the present simulations, the LBM grid size 
$(D_{p}/\unicode[STIX]{x0394}x_{LBM})$ and tracer concentration 
$(C_{t})$ were varied for a system with fixed particle Reynolds number 
$(Re_{Part}=10)$ and two separation distances 
$(\hat{\unicode[STIX]{x1D6FF}}=0.07,0.5)$. The conditions were selected since they are expected to be a good representation of the grid insensitivity requirements. Namely, the largest Reynolds number was considered as well as separation distances that place the particle near the edge of the hydrodynamic boundary layer and very close to the wall. Results from LBM-RWPT were also compared to the commercial COMSOL software which employs a body-fit unstructured grid.
When utilizing the link bounce-back method in LBM, the particle boundary nodes are placed half-way along links that connect internal and external nodes; thereby achieving a discrete representation for the spherical surface. The location of the boundary nodes (link midpoints) can be motivated by work completed on Poiseuille flow where the bounce-back method was shown to recover the no-slip boundary condition at the half-way position between two nodes if the channel was aligned with the LBM lattice (Ginzbourg & Adler Reference Ginzbourg and Adler1994). However, for channels at arbitrary angles, the no-slip boundary condition does not lie exactly at the midpoint between two nodes and is a function of the inclination angle and fluid viscosity (Ginzbourg & Adler Reference Ginzbourg and Adler1994). For a curved surface superimposed on a rectangular grid (particle on LBM lattice), a variety of plane angles will be present and the boundary nodes (link midpoints) will not always align with the curved solid surface. To account for the misalignment of boundary nodes and a curved solid surface, interpolation (Chun & Ladd Reference Chun and Ladd2007) or multireflective (Ginzburg & dHumières Reference Ginzburg and dHumières2003) treatments may be employed. Alternatively, since the hydrodynamic diameter of a sphere represented by the link bounce-back method is a function of particle size and viscosity, the effect of misalignment may be accounted for by introducing a hydrodynamic or effective radius that is displaced by 
$\unicode[STIX]{x0394}H$ outward in the radial direction from the link midpoints, and 
$\unicode[STIX]{x0394}H$ is obtained from calibration simulations (Ladd & Verberg Reference Ladd and Verberg2001). We did not apply a hydrodynamic calibration during the grid convergence study but provide a discussion on its effects below. In our simulations, we hold the LBM viscosity fixed 
$(\unicode[STIX]{x1D708}_{LBM}\cong 1/50)$.
Figure 3. The drag coefficient (a) and lift coefficient (b) obtained with the free-stream velocity.
The drag 
$(C_{D})$ and lift 
$(C_{L})$ coefficients
$$\begin{eqnarray}\displaystyle & \displaystyle C_{D,\infty }=\frac{8F_{D}}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}D_{p}^{2}U_{\infty }^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle C_{L,\infty }=\frac{8F_{L}}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}D_{p}^{2}U_{\infty }^{2}} & \displaystyle\end{eqnarray}$$ were computed with the free-stream velocity 
$(U_{\infty })$ and plotted as a function of grid size; see figure 3. For 
$D_{p}/\unicode[STIX]{x0394}x_{LBM}\geqslant 15$, the drag and lift coefficients become nearly constant and exhibit a nearly constant over-prediction of the drag force by ∼7 % and lift force by ∼5 % when compared to COMSOL (black markers in figure 3), which is consistent with boundary displacement. Non-monotonic trends in the computed drag and lift force are a result of changes to the boundary displacement with particle size, which become nearly constant for higher resolutions (Ladd & Verberg Reference Ladd and Verberg2001), as well as the location of the particle with respect to the underlying grid at a specified separation distance. At a fixed resolution of 
$D_{p}/\unicode[STIX]{x0394}x_{LBM}=15$, the tracer concentration was varied to assess its impact on the heat flux 
$(q^{\prime \prime })$; see figure 4. The heat fluxes are plotted with 95 % confidence intervals that correspond to the stochastic (temporal) fluctuations at steady state. The mean heat flux is observed to be nearly constant while the fluctuations (confidence intervals) decrease as 
$C_{t}$ is increased. RWPT is in very good agreement with COMSOL for both separation distances and captures the heat flux to within ∼1 %. We note that the advective portion of the tracer flux in RWPT employs a fluid velocity that is interpolated. Therefore, velocity fluctuations in the near-surface region, arising from the numerically rough particle, will be attenuated. Additionally, the Prandtl number in this study is less than unity, making the diffusive portion of the tracer flux more significant. Thus, the fluid velocity fluctuations arising from a numerically rough particle do not have a primary effect on the tracer heat flux due to interpolation of the tracer velocity and the considered Prandtl number, which is why we observe better agreement between RWPT and COMSOL. From these results, a resolution of 
$D_{p}/\unicode[STIX]{x0394}x_{LBM}=15$ and 
$C_{t}=1.0$ was chosen.
As the particle resolution increases, the links will become a better approximation of the physical surface of the particle; however, the effect of the lattice viscosity will not be removed. The error between the present LBM simulations and COMSOL are ∼5–7 % for drag and lift, with LBM yielding an over-prediction. To gauge the effect of boundary displacement, a calibrated hydrodynamic radius 
$(\unicode[STIX]{x0394}H=0.25)$ was utilized for a simulation at the chosen resolution 
$(D_{p}/\unicode[STIX]{x0394}x_{LBM}=15;C_{t}=1.0)$ and a separation distance of 
$\hat{\unicode[STIX]{x1D6FF}}=0.5$. The resulting drag/lift forces are observed to decrease by ∼5 % and the heat flux increases by ∼3 %; where it is noted that the heat flux in RWPT is always computed with the physical particle size 
$(R_{p}=7.5)$. The reduction in drag/lift force, due to the use of a calibrated radius, would largely account for the observed discrepancy between LBM and COMSOL. We have conducted further uncalibrated LBM simulations and obtain similar convergence properties for pressure driven flow with and without walls. Additionally, other studies in the literature show very similar behaviour when a calibration is not employed with the link bounce-back method (Kruggel-Emden et al. Reference Kruggel-Emden, Kravets, Suryanarayana and Jasevicius2016). The application of 
$\unicode[STIX]{x0394}H$ will shift the link midpoints (no-slip boundary; 
$R_{p}=7.25$) inward from the physical particle surface 
$(R_{p}=7.5)$. However, the shift of the effective no-slip boundary may lead to regions near the physical particle surface (where RWPT counts the flux) that have finite fluid velocities in the direction normal to the surface. The fluid velocities will contribute to non-physical enhancement of the heat flux due to advective transport of the tracers. In an effort to keep the no-slip boundary location consistent with RWPT, it was chosen to not apply a hydrodynamic calibration in the present work.
Figure 4. The heat flux obtained from LBM-RWPT.
6 Results
6.1 Boundary layer considerations
In the present work, a hydrodynamic and thermal boundary layer develops near the bottom wall. From boundary layer theory, the ratio of the wall thermal boundary layer thickness 
$(\unicode[STIX]{x1D6FF}_{T})$ to the wall hydrodynamic boundary layer thickness 
$(\unicode[STIX]{x1D6FF}_{H})$ is found to scale as 
$\unicode[STIX]{x1D6FF}_{T}/\unicode[STIX]{x1D6FF}_{H}=Pr^{-1/3}$ (White Reference White2005; Schlichting & Gersten Reference Schlichting and Gersten2017). For the Prandtl number considered in this work (0.7), the thermal boundary layer thickness will be larger than the hydrodynamic boundary layer thickness by approximately 12 %. The extent to which the particle interacts with the thermal boundary layer depends upon the separation distance 
$(\unicode[STIX]{x1D6FF})$; see figures 5–6 for velocity and temperature profiles, respectively. For the case of a small separation distance (figures 5a–6a), the particle is within the boundary layers and will interact with the wall to a great degree. By contrast, for large separation distances (figures 5b–6b), the particle is outside the boundary layers and will have a small interaction with the wall.
Figure 5. The dimensionless streamwise velocity 
$(\hat{U} _{x})$ profile for 
$Re_{Part}=10$ and a separation distance of 
$\unicode[STIX]{x1D6FF}/R_{p}=0.5$ (a) and 
$\unicode[STIX]{x1D6FF}/R_{p}=10.0$ (b). The black sphere indicates the location of the particle.
Figure 6. The dimensionless temperature 
$(\unicode[STIX]{x1D703})$ profile for 
$Re_{Part}=10$ and a separation distance of 
$\unicode[STIX]{x1D6FF}/R_{p}=0.5$ (a) and 
$\unicode[STIX]{x1D6FF}/R_{p}=10.0$ (b). The black sphere indicates the location of the particle.
When considering the task of developing correlations from highly resolved methods for use in unresolved computational fluid dynamics discrete element methods (CFD-DEM), the definition adopted for the driving force becomes significant. The driving forces utilized to develop drag and convection correlations in unbounded systems, such as Ranz & Marshall (Reference Ranz and Marshall1952), Whitaker (Reference Whitaker1972), Haider & Levenspiel (Reference Haider and Levenspiel1989), Feng & Michaelides (Reference Feng and Michaelides2000) and Richter & Nikrityuk (Reference Richter and Nikrityuk2012), are the difference between the free-stream fluid velocity and particle velocity 
$(\unicode[STIX]{x0394}U=U_{\infty }-U_{p})$ and the difference between the free-stream fluid temperature and particle temperature 
$(\unicode[STIX]{x0394}T=T_{\infty }-T_{p})$, respectively. In the present work, we show that utilizing the classic definition for the thermal driving force inherently neglects the effect of the wall temperature 
$(T_{w})$ and makes the resulting 
$h$ values specific to the thermal boundary condition. By contrast, the local field surrounding the particle (fluid temperature 
$T_{Loc}$; fluid velocity 
$U_{Loc}$) is a function of all system boundary conditions. Accordingly, our analysis indicates that the use of local mean variables in the present work allows results to be extended to other systems and is a natural choice for spatially varying flows.
Local mean variables may be approximated by the integral of a point variable 
$(T_{f};U_{f})$, with respect to a weighting function 
$(g(r))$, over the fluid volume surrounding the particle 
$(\unicode[STIX]{x1D734}_{y})$ (Deen et al. Reference Deen, Kriebitzsch, van der Hoef and Kuipers2012; Capecelatro & Desjardins Reference Capecelatro and Desjardins2013; Tavassoli et al. Reference Tavassoli, Kriebitzsch, van der Hoef, Peters and Kuipers2013),
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}T_{Loc}={\displaystyle \frac{\displaystyle \iiint g(|\boldsymbol{r}_{y}-\boldsymbol{r}_{p}|)T_{f}(\boldsymbol{r}_{y})\,\text{d}\unicode[STIX]{x1D734}_{y}}{\displaystyle \iiint g(|\boldsymbol{r}_{y}-\boldsymbol{r}_{p}|)\,\text{d}\unicode[STIX]{x1D734}_{y}}},\\[30.0pt] U_{Loc}={\displaystyle \frac{\displaystyle \iiint g(|\boldsymbol{r}_{y}-\boldsymbol{r}_{p}|)U_{f}(\boldsymbol{r}_{y})\,\text{d}\unicode[STIX]{x1D734}_{y}}{\displaystyle \iiint g(|\boldsymbol{r}_{y}-\boldsymbol{r}_{p}|)\,\text{d}\unicode[STIX]{x1D734}_{y}}},\\[30.0pt] g(r)=\exp \left(-{\displaystyle \frac{r}{R_{p}}}\right),\quad {\displaystyle \frac{r}{R_{p}}}\in [0~5],\end{array}\right\}\end{eqnarray}$$ where 
$T_{f}(\boldsymbol{r}_{y})$ is the fluid temperature, 
$U_{f}(\boldsymbol{r}_{y})$ is the fluid velocity, 
$g(r)$ is the integration kernel, 
$\boldsymbol{r}_{p}$ is the location of the particle centre and 
$\unicode[STIX]{x1D734}_{y}$ is the fluid volume within a sphere of radius 
$5R_{p}$ whose centre coincides with 
$\boldsymbol{r}_{p}$. The fluid volume (
$\unicode[STIX]{x1D734}_{y}$) is contained within 
$r/R_{p}\in [1~5]$ while the compact support (
$\unicode[STIX]{x1D734}$) of 
$g(r)$ is contained within 
$r/R_{p}\in [0~5]$. When deriving the volume-averaged equations of motion, Anderson & Jackson (Reference Anderson and Jackson1967) note that local mean variables (6.1) are not uniquely determined unless there is a sufficient separation of scales between the particle–particle spacing and variations in the macroscopic system; which makes results insensitive to the chosen integration kernel 
$(g(r))$ and its characteristic width (
$\unicode[STIX]{x1D70E}_{1/2}$). Deen et al. (Reference Deen, Kriebitzsch, van der Hoef and Kuipers2012) utilized the weighting function given here to approximate the local fluid temperature in DNS simulations of flow through a stationary array of spheres and suggests that a kernel support of radius 
$5R_{p}$ or greater leads to constant heat transfer coefficients. Additionally, the form of 
$g(r)$ given here displays properties that are consistent with the assumptions commonly employed when deriving the volume-averaged equations of motion for a gas–solids mixture (Anderson & Jackson Reference Anderson and Jackson1967; Jackson Reference Jackson1997; Capecelatro & Desjardins Reference Capecelatro and Desjardins2013). Namely, for 
$r>0$, 
$g(r)$ monotonically decreases, is differentiable for all degrees of freedom 
$(C^{\infty })$, and has a characteristic width of 
$\unicode[STIX]{x1D70E}_{1/2}$; which physically corresponds to the radial coordinate at which the normalized integral of 
$g(r)$ is equal to one half (Anderson & Jackson Reference Anderson and Jackson1967) (see (A 1)). Direct computation of local variables, by applying (6.1) to DNS point data, raises the question of whether or not results are sensitive to the integration kernel parameters. The effect of kernel width 
$(\unicode[STIX]{x1D70E}_{1/2})$ and compact support (
$\unicode[STIX]{x1D734}$) are discussed in detail in appendix A for the case of unbounded flow past a cold, static sphere. Local variables computed from DNS point data are shown to depend upon the kernel width but are insensitive to the kernel support, so long that the kernel support is large enough that the truncated tails of 
$g(r)$ do not contribute significantly to the integral (see (A 4)). Therefore, correlations derived with local mean variables will be specific to the kernel width employed 
$(\unicode[STIX]{x1D70E}_{1/2})$. Similarly, in unresolved CFD-DEM methods, the size of the fluid grid will dictate the filtering length scale for interphase transport as well as the magnitude of resolved and unresolved fluid stresses. For volume filtered approaches like those given by Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013), the filter length scale is directly set by the Gaussian kernel width. The integration kernel employed here has a 
$\unicode[STIX]{x1D70E}_{1/2}\approx 2.5R_{p}$, which is similar to the grid size commonly employed in classic CFD-DEM as well as the mollification kernel employed by Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013), a normalized integral on 
$\unicode[STIX]{x1D734}$ of ≈0.90, and yields results very similar to the Gaussian kernel employed by Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013) (see figure 22 versus 24 in appendix A for comparison). Here, we define 
$\unicode[STIX]{x0394}T_{Loc}=T_{Loc}-T_{p}$ as the thermal driving force for convection 
$(h=q^{\prime \prime }/\unicode[STIX]{x0394}T_{Loc})$, and 
$\unicode[STIX]{x0394}U_{Loc}=U_{Loc}-U_{p}$ as the hydrodynamic driving force for drag 
$(C_{D}=8F_{D}/\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}D_{p}^{2}U_{Loc}^{2})$.
Figure 7. (a) Comparison between unbounded Nusselt numbers when the free-stream (
$Nu_{\infty }$) or local fluid temperature (
$Nu_{Loc}$) was utilized. (b) Comparison between the unbounded drag coefficients when the free-stream (
$C_{D,\infty }$) or local fluid velocity (
$C_{D,Loc}$) was utilized.
Physically speaking, as 
$\unicode[STIX]{x1D6FF}$ becomes large with respect to the wall boundary layer thicknesses (figures 5b–6b), the boundary effects on particle heat and momentum transfer will become negligible. In the limit of 
$\unicode[STIX]{x1D6FF}\rightarrow \infty$, the resulting drag coefficient 
$(C_{D})$ and Nusselt number 
$(Nu=hD_{p}/k_{g})$ should converge to those obtained for an unbounded system (Ranz & Marshall Reference Ranz and Marshall1952; Whitaker Reference Whitaker1972; Haider & Levenspiel Reference Haider and Levenspiel1989; Feng & Michaelides Reference Feng and Michaelides2000; Richter & Nikrityuk Reference Richter and Nikrityuk2012). However, use of a local variable 
$(T_{Loc};U_{Loc})$ will prevent 
$C_{D}$ and 
$Nu$ from converging to the classic unbounded correlations even in the limit of 
$\unicode[STIX]{x1D6FF}\rightarrow \infty$. To illustrate this concept, simulations of unbounded, uniform flow past a sphere were completed; see figure 7. The disagreement between quantities computed with local variables 
$(C_{D,Loc};Nu_{Loc})$ and existing correlations for unbounded systems that utilize 
$U_{\infty }$ and 
$T_{\infty }$ can be attributed to a reduction in driving force based on local variables 
$(\unicode[STIX]{x0394}T_{Loc}<\unicode[STIX]{x0394}T;U_{Loc}<U_{\infty })$. Similarly, unresolved CFD-DEM frameworks experience the same phenomena, where fluid disturbances due to the presence of the particle (interphase drag and heat transfer) yield interpolated fluid quantities that differ from the undisturbed (free-stream) values. When examining two-way coupled point-particle methods, it was observed that the undisturbed fluid velocity can be of greater significance than using the most appropriate drag force model (Mehrabadi et al. Reference Mehrabadi, Horwitz, Subramaniam and Mani2018). Efforts to improve unresolved CFD-DEM drag have been focused on velocity correction schemes that are derived in the Stokes limit (Horwitz & Mani Reference Horwitz and Mani2016; Ireland & Desjardins Reference Ireland and Desjardins2017; Horwitz & Mani Reference Horwitz and Mani2018) or the Oseen approximation (Balachandar, Liu & Lakhote Reference Balachandar, Liu and Lakhote2019). By contrast, the correlations obtained here employ local variables for consistency with the driving forces available to an unresolved CFD-DEM simulation that utilizes the same kernel width employed here. Specifically, the Nusselt numbers obtained here with 
$\unicode[STIX]{x0394}T=T_{\infty }-T_{p}~(Nu_{\infty })$ agree with classic, unbounded convection correlations while the Nusselt numbers obtained with 
$\unicode[STIX]{x0394}T_{Loc}=T_{Loc}-T_{p}~(Nu_{Loc})$ are larger. Similarly, the drag coefficients obtained with 
$U_{Loc}$
$(C_{D,Loc})$ are larger than the correlation of Haider & Levenspiel (Reference Haider and Levenspiel1989) while the drag coefficients obtained with 
$U_{\infty }~(C_{D,\infty })$ agree with Haider & Levenspiel (Reference Haider and Levenspiel1989). Shifts in the abscissa of figure 7 when computing local quantities 
$(C_{D,Loc};Nu_{Loc})$ are a result of also converting 
$Re_{Part}$ to a local quantity 
$(Re_{Part,Loc}=U_{Loc}D_{p}/\unicode[STIX]{x1D708})$. Figure 7 shows that local variables lead to self-similar shifts in the Nusselt number and drag coefficient data, which implies that the functional forms for the unbounded correlations still hold but with different scaling. Therefore, the unbounded correlations of Richter & Nikrityuk (Reference Richter and Nikrityuk2012) and Haider & Levenspiel (Reference Haider and Levenspiel1989) are re-fit to the 
$Nu_{Loc}$ and 
$C_{D,Loc}$ data to obtain
$$\begin{eqnarray}\displaystyle & \displaystyle Nu_{Loc,UB}\equiv \frac{h_{Loc,UB}D_{p}}{k_{g}}=3.4+0.03Re_{Part,Loc}^{0.4}+0.18Re_{Part,Loc}^{0.7}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle C_{D,Loc,UB}\equiv \frac{8F_{D}}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}D_{p}^{2}U_{Loc}^{2}}=\frac{38.5}{Re_{Part,Loc}}(1~+0.11Re_{Part,Loc}^{0.70}), & \displaystyle\end{eqnarray}$$ where 
$Nu_{Loc,UB}$ is the local, unbounded Nusselt number, 
$C_{D,Loc,UB}$ is the local, unbounded drag coefficient and 
$Re_{Part,Loc}$ is the local particle Reynolds number. Since the correlations obtained in (6.2)–(6.3) employ local variables, they are more consistent with modified correlations that account for particle disturbances (Municchi, Goniva & Radl Reference Municchi, Goniva and Radl2016; Ireland & Desjardins Reference Ireland and Desjardins2017). Modified drag correlations are obtained from filtering the solution for Stokes flow past a sphere with respect to an integration kernel. For the case of a modified Nusselt correlation, the thermal temperature must also be filtered. In appendix B we give the modified drag correlations arising from the local fluid velocities of Ireland & Desjardins (Reference Ireland and Desjardins2017) and Municchi et al. (Reference Municchi, Goniva and Radl2016). In appendix C we show that the hydrodynamic disturbances in appendix B may be utilized to obtain thermal disturbances but the thermal disturbances will correspond to a sphere undergoing steady diffusion in the radial direction. Additionally, in appendix C we also consider filtering the asymptotic solution of Acrivos & Taylor (Reference Acrivos and Taylor1962) to estimate the effect of finite Reynolds numbers on the thermal disturbances. This case is referred to as ‘thermal Stokes filtered’ and contains an ‘inner’ and ‘outer’ solution that correspond to the near-surface region of the particle and the far field, respectively. The modified drag laws of Ireland & Desjardins (Reference Ireland and Desjardins2017) and Municchi et al. (Reference Municchi, Goniva and Radl2016) agree with the correlation developed here (6.3) at lower Reynolds numbers but tend to predict larger 
$C_{D}$ as the Reynolds number increases; see figure 8. The same effect was reported in figure 7 of Ireland & Desjardins (Reference Ireland and Desjardins2017) where their filtered fluid velocity (based on Stokes flow) led to over-predictions for the residual drag force at Reynolds numbers of 
$O(1)$ but was highly accurate for Reynolds numbers up to 
$O(10^{-1})$. At higher Reynolds numbers, the local velocity will be under-predicted by the filtering of Stokes flow and the filtered drag correlation will be shifted upward too far, leading to the over-prediction of the drag coefficient. Therefore, the departure of (6.3) from the modified correlations of Ireland & Desjardins (Reference Ireland and Desjardins2017) and Municchi et al. (Reference Municchi, Goniva and Radl2016) at higher Reynolds numbers has a physical basis. Additionally, the modified Nusselt correlations of Ireland & Desjardins (Reference Ireland and Desjardins2017) and Municchi et al. (Reference Municchi, Goniva and Radl2016) show reasonable agreement with the correlation developed here (6.2) as the particle Reynolds number tends to zero but yield significantly larger 
$Nu$ values at finite Reynolds numbers; see figure 9. The larger disagreement between (6.2) and the Nusselt correlations of Ireland & Desjardins (Reference Ireland and Desjardins2017) and Municchi et al. (Reference Municchi, Goniva and Radl2016) can be attributed to under-predictions in both the disturbed fluid velocity and the disturbed fluid temperature at elevated Reynolds numbers.
Figure 8. Comparison between filtered drag correlations with the same kernel attributes as that employed here over (a) the simulation conditions and (b) outside the simulation conditions.
Figure 9. Comparison between filtered Nusselt correlations with the same kernel attributes as that employed here over (a) the simulation conditions and (b) outside the simulation conditions.
In the opposite limit of separation distance 
$(\unicode[STIX]{x1D6FF}\rightarrow 0)$, a choice must be made in terms of the definition for 
$\unicode[STIX]{x1D734}_{y}$ appearing in (6.1). Since the radius of 
$\unicode[STIX]{x1D734}_{y}$ is 
$5R_{p}$ (significantly larger than the particle), a subset of 
$\unicode[STIX]{x1D734}_{y}$ will overlap with the wall (
$\unicode[STIX]{x1D734}_{w}$). For this case, the volume of 
$\unicode[STIX]{x1D734}_{y}$ overlap with the wall (
$\unicode[STIX]{x1D734}_{w}$) as well as the fluid volume (
$\unicode[STIX]{x1D734}_{f}$) was incorporated into the volume integration performed in (6.1) (
$\unicode[STIX]{x1D734}_{y}=\unicode[STIX]{x1D734}_{f+w}$) and the variables within 
$\unicode[STIX]{x1D734}_{w}$ were set to the boundary condition values 
$(\unicode[STIX]{x1D703}=1;\boldsymbol{u}=\mathbf{0})$. An alternative approach for handling the near-wall region that involves mirroring the integration kernel across the wall was utilized by Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013) to obtain a Neumann boundary condition. By contrast, we employ a Dirichlet boundary condition here. If the integration kernel is mirrored across the bottom wall in the present work, then the variable to be averaged 
$(T_{f};U_{f})$ must still be specified within the wall volume. Setting variables within the wall volume to their boundary condition quantities 
$(T_{w};U_{w})$ results in a treatment that is quite similar to the one proposed here (integration in (6.1) over 
$\unicode[STIX]{x1D734}_{f+w}$). However, the mirroring of the kernel across a Dirichlet boundary may lead to drastic changes in the local mean variable as the integration kernel begins to intersect the wall – i.e. the original kernel will average the temperature around the particle but the mirrored kernel will solely average the wall temperature. The choice to employ (6.1) with an integration over 
$\unicode[STIX]{x1D734}_{f+w}$ was motivated by the interpolation techniques employed within CFD-DEM (Garg et al. Reference Garg, Galvin, Li and Pannala2012), to which our correlations are to be applied. In CFD-DEM, the fluid velocity and temperature are found by interpolation of the CFD numerical grid to the centre of the particle. For a numerical cell adjacent to a wall, the fluid variables will be interpolated between the known wall nodes 
$(T_{w};U_{w})$ and the solution values at the adjacent nodes within the domain. By including 
$\unicode[STIX]{x1D734}_{w}$ into the calculation of local variables, the resulting values smoothly approach the boundary condition values as the separation distance decreases and are consistent with quadratic interpolation techniques in CFD-DEM; see figure 10 for the effect on 
$T_{Loc}$.
Figure 10. The local, dimensionless fluid temperature 
$(\unicode[STIX]{x1D703}_{Loc})$ found via integration over 
$\unicode[STIX]{x1D734}_{y}$ that includes the wall volume 
$(\unicode[STIX]{x1D734}_{f+w})$ (solid lines) as well as just the fluid volume 
$(\unicode[STIX]{x1D734}_{f})$ (dashed lines) versus the dimensionless distance between the particle centre and the wall 
$(\hat{\unicode[STIX]{x1D6FF}}+1)$. Note, for 
$\hat{\unicode[STIX]{x1D6FF}}+1\geqslant 5$, 
$\unicode[STIX]{x1D734}_{y}$ does not intersect the wall 
$(\unicode[STIX]{x1D734}_{w}=0\rightarrow \unicode[STIX]{x1D734}_{y}=\unicode[STIX]{x1D734}_{f})$. The wall temperature 
$(\unicode[STIX]{x1D703}_{w})$ is included for reference.
Figure 11. The simulated drag force for each particle–wall separation distance and Reynolds number.
Figure 12. The drag coefficient computed with 
$U_{\infty }$ (a) and 
$U_{Loc}$ (b).
Figure 13. The normalized drag coefficient computed with 
$U_{\infty }$ (a) and 
$U_{Loc}$ (b).
6.2 Near-wall drag
The primary drag force on the particle acts in the streamwise direction (positive x-direction) and was extracted from each simulation at steady state; see figure 11. As the dimensionless separation distance 
$(\hat{\unicode[STIX]{x1D6FF}}=\unicode[STIX]{x1D6FF}/R_{p})$ decreases, the particle begins to interact with the slower moving fluid contained within the wall hydrodynamic boundary layer and the drag force decreases. Additionally, the drag reduction in the near-wall region occurs at smaller separation distances for increasing Reynolds number due to compression of the wall hydrodynamic boundary layer. The drag coefficient is computed from the drag force with the free-stream fluid velocity 
$(U_{\infty };C_{D,\infty })$ as well as the local fluid velocity 
$(U_{Loc};C_{D,Loc})$; see figure 12. Use of the free-stream velocity results in drag coefficients that decrease in the near-wall region, due to the reduced drag forces in the near-wall region. By contrast, use of the local fluid velocity results in the exact opposite behaviour since the local fluid velocity also reduces as the particle becomes closer to the wall. To gauge the drag reduction in the near-wall region, with respect to an unbounded system, the drag coefficients (
$C_{D,\infty }$ and 
$C_{D,Loc}$) are normalized by the unbounded correlations (see Haider & Levenspiel (Reference Haider and Levenspiel1989) and (6.3), respectively); see figure 13. Normalized drag coefficients 
$({\hat{C}}_{D})$ that employ the free-stream velocity are a function of 
$\unicode[STIX]{x1D6FF}$ and decrease in the near-wall region (figure 13a). By contrast, use of the local fluid velocity yields normalized drag coefficients that are nearly unity for all separation distances – i.e. the drag force for a particle in the near-wall region is well approximated by an unbounded particle with the same average environment (figure 13b). Local maxima in the data of figure 13(b) correspond to particles near the edge of the wall hydrodynamic boundary layer and persist at separation distances where the averaging volume does not intersect the wall. The bumps are a product of volume averaging (6.1) into the boundary layer and can be attenuated by reducing the kernel width in (6.1), thereby confining the averaging to more localized regions about the particle. Despite small mismatches, equation (6.3) (unbounded correlation modified for local velocity field) yields very reasonable predictions in both the near-wall region as well as far from the boundary. By contrast, use of the Haider & Levenspiel (Reference Haider and Levenspiel1989) correlation with the free-stream fluid velocity will lead to large over-predictions in the near-wall region and the errors will become more significant as the separation distance is reduced (figure 13a). Therefore, the local variable treatment may permit simple modifications to classic correlations in order to account for wall proximity, rather than redevelopment.
Figure 14. Normalized drag coefficients (computed with 
$U_{\infty }$) compared to (6.5)
Since the shift to local variables is a new concept and previous works have by and large been completed with free-stream variables, we also consider the task of constructing a compression of the data that were normalized with free-stream variables (figure 13a). From the physical arguments given in the Boundary Layer Considerations section (§ 6.1), there is an expectation that the drag reduction is due to the interaction of the particle with the slower moving fluid in the wall boundary layer and that the drag coefficient will asymptotically approach the unbounded correlation for large particle–wall separation distances. Therefore, the ratio of the separation distance to the hydrodynamic boundary layer 
$(\hat{\unicode[STIX]{x1D6FF}}_{H}=\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}_{H})$ is expected to be a crucial length scale, where 
$\unicode[STIX]{x1D6FF}_{H}$ may be approximated from boundary layer theory as (White Reference White2005)
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{H}\approx 5.0\frac{x}{Re_{Plate}^{1/2}}.\end{eqnarray}$$ Plotting the normalized drag coefficient with free-stream fluid velocity against 
$\hat{\unicode[STIX]{x1D6FF}}_{H}$ shows that the data roughly compress on a single curve; see figure 14. The compressed data show a shape that is consistent with the analytical solution for flow past the leading edge of a plate (Jessee Reference Jessee2015). More specifically, an adjusted form of the solution in Jessee (Reference Jessee2015) is observed to replicate the data well and asymptotically approaches unity as 
$\hat{\unicode[STIX]{x1D6FF}}_{H}\rightarrow \infty$ (the physically correct limit)
$$\begin{eqnarray}{\hat{C}}_{D,\infty }(\hat{\unicode[STIX]{x1D6FF}}_{H})=0.4+\frac{0.6}{\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\right)^{2}-1}\left[\frac{\unicode[STIX]{x03C0}}{2}\left(\tan ^{-1}(\hat{\unicode[STIX]{x1D6FF}}_{H})+\frac{0.7\hat{\unicode[STIX]{x1D6FF}}_{H}}{(0.7\hat{\unicode[STIX]{x1D6FF}}_{H})^{2}+1}\right)+\frac{0.7}{(0.7\hat{\unicode[STIX]{x1D6FF}}_{H})^{2}+1}-1\right].\end{eqnarray}$$ It must be mentioned that the 0.4 and 0.6 constants in (6.5) govern the drag reduction at particle–wall contact 
$({\hat{C}}_{D,\infty }(0)=0.4)$ and in the limit of large separation distance 
$({\hat{C}}_{D,\infty }(\infty )=0.4+0.6)$ and are expected to be a function of 
$R_{p}/\unicode[STIX]{x1D6FF}_{H}$. The robustness of (6.5) (wall function for the normalized drag coefficient with free-stream fluid velocity) and (6.3) (new, unbounded drag correlation with local fluid velocity) will be tested in the Model Assessment section (§ 7) by conducting simulations at larger leading-edge lengths and thus larger hydrodynamic boundary layer thicknesses (i.e. different 
$R_{p}/\unicode[STIX]{x1D6FF}_{H}$).
Figure 15. The lift coefficients computed with 
$U_{\infty }$ (a) and 
$U_{Loc}$ (b).
6.3 Near-wall lift
For uniform, unbounded flows at the Reynolds numbers considered here 
$(Re_{Part}<210)$, the lift force on a sphere is essentially zero due to the axisymmetric velocity field (Bagchi, Ha & Balachandar Reference Bagchi, Ha and Balachandar2000). By contrast, shear and vortex flows have been shown to generate considerable lift forces (Bagchi & Balachandar Reference Bagchi and Balachandar2002). Segre and Silberberg showed that particles in Poiseuille flows tend to move away from walls and reach an equilibrium separation distance (Segre & Silberberg Reference Segre and Silberberg1962a,Reference Segre and Silberbergb). For asymptotically small Reynolds numbers, Saffman and Mcloughlin derived closed forms for the lift force (Saffman Reference Saffman1965; McLaughlin Reference McLaughlin1993). For a stationary sphere in a linear shear flow, the lift force will always act towards the side with higher fluid velocity (Kurose & Komori Reference Kurose and Komori1999). Due to the shearing within the hydrodynamic boundary layer, it is expected that significant lift forces will be present here and that they will act in the positive y-direction (towards the faster moving free-stream fluid). The lift coefficients in figure 15 show that the lift forces in the present simulations are of substantial magnitude and that they will act to push particles away from the wall (force is in the vertical y-direction, positive 
$C_{L}$), which is consistent with Kurose & Komori (Reference Kurose and Komori1999). Use of the free-stream velocity leads to lift coefficients that increase with increasing separation distance but reach a maximum at some critical separation distance 
$(\unicode[STIX]{x1D6FF}_{crit})$, after which the lift coefficients decrease. Use of the local fluid velocity leads to trends that are loosely similar to those obtained with the free-stream fluid velocity but the maxima are shifted closer to the wall and 
$\unicode[STIX]{x1D6FF}_{crit}$ does not increase to the same degree with decreasing Reynolds number. The maxima in the lift force are expected to correspond to particles near the edge of the boundary layer where the product of shear stress and fluid velocity is greatest. Plotting the lift coefficients against 
$\hat{\unicode[STIX]{x1D6FF}}_{H}$ shows that the maxima in 
$C_{L,\infty }$ do indeed correspond to 
$\hat{\unicode[STIX]{x1D6FF}}_{H}\cong 1.0$; see figure 16. While interesting, the task of developing a near-wall correlation for the lift force is beyond the scope of the present work and is not attempted here.
Figure 16. Lift coefficients computed with 
$U_{\infty }$ when the separation distance is normalized by the hydrodynamic boundary layer thickness.
6.4 Near-wall heat transfer
For each LBM-RWPT simulation, the heat rate to the particle 
$(\dot{Q})$ and local fluid temperature 
$(T_{Loc})$ was extracted at steady state. The heat rates obtained from LBM-RWPT 
$(\dot{Q})$ are directly compared to unbounded convection correlations 
$(\dot{Q}_{conv})$ and the indirect conduction theory 
$(\dot{Q}_{PFW})$ commonly employed in CFD-DEM methods. First, the unbounded convective correlation of Ranz & Marshall (Reference Ranz and Marshall1952) with a local fluid temperature 
$(\dot{q}_{conv}=h_{conv}A_{p}\unicode[STIX]{x0394}T_{Loc})$ is compared to LBM-RWPT; see figure 17(a). As the particle–wall separation distance becomes small, the heat transfer coefficient grows rapidly (note logarithmic 
$x$-axis) and the unbounded convection correlation fails to characterize the heat transfer enhancement that occurs in the near-wall region. This behaviour is expected since the correlation given in Ranz & Marshall (Reference Ranz and Marshall1952) (unbounded system) does not account for the thermal source associated with the boundary. Note that the dimensionless heat rate 
$(\hat{q})$ does not decay to unity as the separation distance becomes large. This behaviour is solely a result of utilizing 
$\unicode[STIX]{x0394}T_{Loc}$ as the thermal driving force (see 
$Nu_{Loc}$ in figure 7a) and 
$\hat{q}$ would tend to unity if 
$\unicode[STIX]{x0394}T=T_{\infty }-T_{p}$ were utilized for the thermal driving force. Furthermore, 
$\hat{q}$ decreases with increasing 
$Re_{Part}$, and indicates that temperature gradients around the particle are intensified by increased fluid velocity, which makes 
$\unicode[STIX]{x0394}T_{Loc}\rightarrow \unicode[STIX]{x0394}T$.
Inclusion of the indirect particle-fluid-wall (PFW) conduction mechanism (Rong & Horio Reference Rong and Horio1999; Lattanzi & Hrenya Reference Lattanzi and Hrenya2017) into the total heat rate 
$(\dot{Q}_{conv}+\dot{Q}_{PFW}=h_{conv}A_{p}\unicode[STIX]{x0394}T_{Loc}+k_{g}R_{p}{\hat{h}}_{PFW}(\hat{\unicode[STIX]{x1D6FF}})[T_{w}-T_{p}])$ is observed to agree markedly better with LBM-RWPT than the convection correlation alone; see figure 17(b). In contrast to the unbounded convection correlation, indirect conduction theory accounts for the effect of a boundary by assuming that one-dimensional conduction occurs through a stagnant layer of fluid between the particle and wall 
$(R_{Lens})$. However, heat transfer enhancement due to the hot wall is still observed at length scales not predicted by indirect conduction theory (peaks in figure 17b). The length scale for indirect conduction theory is the fluid lens thickness 
$(R_{Lens}-R_{p})$ and is set according to the particle size 
$(R_{Lens}=1.4R_{p})$ (Lattanzi & Hrenya Reference Lattanzi and Hrenya2017) – i.e. 
$\dot{Q}_{PFW}$ only contributes to the total heat rate when 
$\hat{\unicode[STIX]{x1D6FF}}<(R_{Lens}-R_{p})/R_{p}=0.4$ to the left of the peaks in figure 17(b).
Figure 17. The heat rate obtained via LBM-RWPT 
$(\dot{Q})$ normalized by the correlation of Ranz & Marshall (Reference Ranz and Marshall1952) (a) and the sum of Ranz & Marshall (Reference Ranz and Marshall1952) and indirect conduction theory (Rong & Horio Reference Rong and Horio1999) 
$(\dot{Q}_{conv}+\dot{Q}_{PFW})$ (b) versus dimensionless separation distance 
$(\hat{\unicode[STIX]{x1D6FF}})$.
Physically speaking, heat transfer enhancement due to the boundary should occur at a length scale associated with the thermal boundary layer thickness 
$(\unicode[STIX]{x1D6FF}_{T})$ of the plate, rather than the particle radius; see figure 18. For example, if a particle that is large with respect to 
$\unicode[STIX]{x1D6FF}_{T}$ is considered (right particle in figure 18), the onset of indirect conduction (fluid lens just intersects the wall; 
$\unicode[STIX]{x1D6FF}=0.4R_{p}$) would correspond to a particle outside the thermal boundary layer. For this case, the inherent assumptions of indirect conduction theory (static, one-dimensional conduction) are violated since the hot fluid contained within the thermal boundary layer is advected between the particle and wall. This scenario may correspond to large particle sizes and/or compression of the thermal boundary layer with the Reynolds and Prandtl numbers. The advection of fluid between the particle and the wall acts to reduce the thermal gradients near the particle surface from those predicted by indirect conduction theory, and thus, the heat transfer to the particle in this case is over-predicted by indirect conduction theory. By contrast, if a particle small with respect to 
$\unicode[STIX]{x1D6FF}_{T}$ is considered (left particle in figure 18), the onset of indirect conduction (fluid lens just intersects the wall; 
$\unicode[STIX]{x1D6FF}=0.4R_{p}$) corresponds to a particle that is fully immersed in the boundary layer. Therefore, the heat transfer enhancement occurring when the particle is within the thermal boundary layer 
$(\unicode[STIX]{x1D6FF}<\unicode[STIX]{x1D6FF}_{T})$ but not within the fluid lens thickness 
$(\unicode[STIX]{x1D6FF}>0.4R_{p})$ cannot be captured by indirect conduction theory – i.e. the particle may reside in the thermal boundary layer where heat transfer enhancement occurs but the fluid lens does not intersect with the wall. In this case, the heat transfer to the particle is under-predicted by indirect conduction theory. Note that the ratio of the particle size to thermal boundary layer thickness considered in the LBM-RWPT simulations here is most analogous to the ‘small’ case in figure 18, which is why the combination of convection and indirect conduction tends to under-predict the overall heat transfer (figure 17b).
Figure 18. An illustration of particles (solid black lines) and their fluid lenses (dashed grey line) overlaid on the scaling of the thermal boundary layer thickness 
$(\unicode[STIX]{x1D6FF}_{T})$. For particles small with respect to 
$\unicode[STIX]{x1D6FF}_{T}$ (left particle), the onset of indirect conduction occurs when the particle is inside the boundary layer. For particles large with respect to 
$\unicode[STIX]{x1D6FF}_{T}$ (right particle), the onset of indirect conduction occurs when the particle is outside the boundary layer.
Clearly the presence of the thermal source at the wall leads to thermal gradients that act to enhance the heat transfer over what would be obtained with a particle in an unbounded system. To quantify the heat transfer enhancement due to interaction between the particle and wall thermal boundary conditions, the heat flux in excess of the local, unbounded convection correlation (6.2) is considered 
$(q_{w}^{\prime \prime }=q^{\prime \prime }-h_{Loc,UB}\unicode[STIX]{x0394}T_{Loc})$. The thermal driving force for the excess heat flux 
$(q_{w}^{\prime \prime })$ is the thermal gradient between the two solid bodies 
$(\unicode[STIX]{x0394}T_{w}=T_{w}-T_{p})$ and is utilized to obtain an excess heat transfer coefficient that is associated with the wall 
$(h_{w}=q_{w}^{\prime \prime }/\unicode[STIX]{x0394}T_{w})$. The excess wall Nusselt number 
$(Nu_{w}=h_{w}D_{p}/k_{g})$ characterizes the additional heat transfer that a particle in the near-wall region will experience when compared to a particle in an unbounded system with the same average conditions. As the separation distance is reduced, large thermal gradients that persist between the particle and wall, where the fluid velocity is smallest, yield heat fluxes that are much larger than those experienced by a particle undergoing steady, unbounded convection and this is why (6.2) with a local fluid temperature is not capable of capturing the overall heat transfer. While the presence of the hot wall shows significant departure from an unbounded system, the excess wall Nusselt numbers compress well onto a single curve (see figure 19) and may be approximated by
$$\begin{eqnarray}Nu_{w}(\hat{\unicode[STIX]{x1D6FF}})\equiv \frac{h_{w}(\hat{\unicode[STIX]{x1D6FF}})D_{p}}{k_{f}}=1.42\exp (-20.12\hat{\unicode[STIX]{x1D6FF}})+1.54\exp (-2.67\hat{\unicode[STIX]{x1D6FF}}).\end{eqnarray}$$By making use of (6.2) and (6.6), the total heat flux to the particle then becomes
$$\begin{eqnarray}q^{\prime \prime }=h_{Loc,UB}\unicode[STIX]{x0394}T_{Loc}+h_{w}(\hat{\unicode[STIX]{x1D6FF}})\unicode[STIX]{x0394}T_{w}.\end{eqnarray}$$ Note that 
$Nu_{w}$ asymptotically decays to zero as 
$\hat{\unicode[STIX]{x1D6FF}}\rightarrow \infty$, which is the physically correct behaviour 
$(Nu\rightarrow Nu_{Loc,UB})$. The under-shoot of the data at 
$\hat{\unicode[STIX]{x1D6FF}}\approx 4.0$ again corresponds to averaging effects near the edge of the boundary layer and are more pronounced with increasing 
$Re_{Part}$, due to the reduction in hydrodynamic and thermal length scales as the system becomes more advection dominated. While not considered here, it is expected that the exponential relaxation rates in (6.6) will display a dependence upon the Prandtl number. As the thermal boundary layers (wall and particle) compress with increasing Prandtl number, the thermal interaction between particle and wall (heat transfer in excess of convection) will be confined to increasingly smaller length scales, leading to larger relaxation rates. The robustness of (6.7) will also be tested in § 7 by modifying the leading-edge length and thermal boundary conditions.
Figure 19. The excess wall Nusselt number obtained via LBM-RWPT versus the non-dimensional separation distance on a logarithmic axis.
Due to the restrictions on parameter space, the formal accuracy of indirect conduction theory for a generic system is outside of the scope of the present work. However, by identifying the thermal boundary layer thickness as the key length scale, some general trends may be noted. For particles that are large with respect to 
$\unicode[STIX]{x1D6FF}_{T}$ (right particle in figure 18), the current indirect conduction theories within DEM are expected to over-predict the heat transfer to the particle. This can be traced back to the violation of the static-fluid lens assumption over a length scale of 
$0.4R_{p}$. Note that the boundary layer thickness can vary spatially and will compress with increasing Reynolds and Prandtl numbers. For the case of a particle that is small with respect to 
$\unicode[STIX]{x1D6FF}_{T}$ (left particle in figure 18), indirect conduction theories are expected to under-predict the heat transfer to the particle (observed here in figure 17b). In this case, the particle is well within the boundary layer (where heat transfer enhancement occurs) at the onset of indirect conduction 
$(\unicode[STIX]{x1D6FF}\leqslant 0.4R_{p})$.
7 Model assessment
Additional simulations were conducted with larger leading-edge lengths and different thermal boundary conditions to test the robustness of the drag ((6.3) (unbounded correlation with local fluid velocity) and (6.5) (wall function with free-stream fluid velocity)) and convection ((6.2) (unbounded Nusselt correlation with local fluid temperature) and (6.6) (excess wall Nusselt number correlation)) correlations developed here from DNS data; see tables 4–5. Due to the additional computational demand imposed by larger simulation domains and longer run times, only 
$Re_{Part}=10$ was considered and the tracer concentration was reduced to 
$C_{t}=0.7$. It was shown in § 5 that a reduction in tracer concentration does not appreciably impact the mean heat flux but does increase stochastic fluctuations to a small degree.
Table 4. LBM-RWPT geometry and mesh.
Table 5. Simulation conditions.
The steady-state drag forces at 
$L/D_{p}=10$ again decreases with decreasing particle–wall separation distance and the observed trend is quite similar to that obtained with 
$L/D_{p}=5$ (figure 11); see figure 20(a). However, the drag reduction at 
$L/D_{p}=10$ is observed at larger separation distances than those at 
$L/D_{p}=5$ due to the growth in the thickness of the wall hydrodynamic boundary layer with increasing leading-edge length 
$(\unicode[STIX]{x1D6FF}_{H}(x))$. The drag force obtained here from DNS simulations was normalized by the drag force predicted by the correlations (
$F_{D,Cor}$; equation (6.3) (new, unbounded drag correlation with local fluid velocity) and (6.5) (wall function for the normalized drag coefficient with free-stream fluid velocity)) in figure 20(b) to assess the accuracy and robustness of both methods. Both correlations perform reasonably well at very small and very large separation distances but (6.3) with the local fluid velocity is observed to better capture the drag occurring at intermediate separation distances that are near the edge of the wall boundary layer.
The steady-state heat fluxes at 
$L/D_{p}=10$ show significant growth as the separation distance is reduced and the trends are again similar to the 
$L/D_{p}=5$ results (figure 17a); see figure 21(a). The heat fluxes were also normalized by the correlation predictions (
$q_{Cor}^{\prime \prime }$; equation (6.7)) and show strong agreement for all thermal boundary conditions and separation distances; see figure 21(b). These observations suggest that the proposed superposition of local, unbounded convection and a wall interaction may sufficiently capture the heat transfer occurring in systems beyond that utilized to construct the correlation.
Figure 21. The heat flux versus non-dimensional separation distance (a) and the heat flux normalized by (6.7) (b).
8 Conclusions
Direct numerical simulation was utilized to examine the effect of a hot boundary on heat and momentum transfer in a gas–solids mixture. Drag forces obtained via LBM-RWPT show significant drag reduction in the near-wall region that is not captured by unbounded drag correlations based on a free-stream fluid velocity. By contrast, the use of an unbounded drag correlation with local fluid velocity is observed to capture the data reasonably well for large and small particle–wall separation distances.
The heat rates obtained via LBM-RWPT show that unbounded convection correlations are not sufficient in the near-wall region but the combination of such correlations with indirect conduction theory agrees markedly better with DNS results. Nonetheless, such particle-scale theories still exhibit discrepancies with DNS, which can be traced to the effect of the thermal length scales; specifically, the length scale associated with near-wall heat transfer enhancement is proportional to the thermal boundary layer thickness and not the particle radius, which is the length scale utilized by indirect conduction theory. The use of an unbounded convection correlation with local fluid temperature is also observed to under-predict the heat transfer in the near-wall region due to sustained temperature gradients arising due to the proximity of the particle and wall thermal boundary conditions. The heat transfer in excess of the unbounded convection correlation with local fluid temperature is utilized to define an excess wall Nusselt number that only depends upon separation distance. Superposition of the local, unbounded Nusselt number and excess wall Nusselt number defines a new correlation which is valid in the near-wall region but asymptotically decays to the unbounded convection correlation in the limit of large particle–wall separation distance.
The volume averaging process utilized here to define local variables is shown to depend upon the width of the integration kernel. As the kernel width is decreased, the point variables near the surface of the particle are more heavily weighted by the filtering process and the resulting local variables depart from the free-stream variables to a greater extent. Therefore, drag and heat transfer correlations based upon volume-averaged quantities will undergo shifts away from correlations based upon undisturbed quantities and the shift will become greater as the kernel width is reduced. The same phenomenon is also observed in unresolved CFD-DEM simulations when the width of the mollification kernel is reduced (Ireland & Desjardins Reference Ireland and Desjardins2017).
While not considered here, the particle(s) may translate in space as well as rotate (angular velocity). Furthermore, the diameter of the particle, Prandtl number and thermal wall boundary condition may be altered. The impact of each parameter on particle heat and momentum transfer is not known a priori but will be the subject of future work.
Acknowledgements
The authors gratefully acknowledge funding support provided by the National Science Foundation under grant CBET 1512630.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Kernel attributes
The computation of local variables (as given in (6.1)) requires that a functional form for the integration kernel 
$(g(r))$ be specified. To maintain consistency with the methods utilized to derive the volume-averaged equations of motion, 
$g(r)$ should monotonically decrease for 
$r>0$, be differentiable for all degrees of freedom 
$(C^{\infty })$ and have a characteristic width of 
$\unicode[STIX]{x1D70E}_{1/2}$. The characteristic width 
$(\unicode[STIX]{x1D70E}_{1/2})$ is defined as the radial coordinate at which the integral of the normalized kernel 
$({\hat{g}}(r))$ is equal to one half (Anderson & Jackson Reference Anderson and Jackson1967)
$$\begin{eqnarray}\int _{\unicode[STIX]{x1D70E}_{1/2}}^{\infty }{\hat{g}}(r)r^{2}\,\text{d}r\equiv \frac{\displaystyle \int _{\unicode[STIX]{x1D70E}_{1/2}}^{\infty }g(r)r^{2}\,\text{d}r}{\displaystyle \int _{0}^{\infty }g(r)r^{2}\,\text{d}r}=\frac{1}{2},\end{eqnarray}$$ where 
${\hat{g}}(r)$ is normalized over the entire space and not just the fluid volume (
$\unicode[STIX]{x1D734}_{y}$), which is why the lower bound of the integral in the denominator of (A 1) starts at 0. Two kernels that meet these criteria are the exponential and Gaussian function (utilized by Deen et al. (Reference Deen, Kriebitzsch, van der Hoef and Kuipers2012) and Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013), respectively)
$$\begin{eqnarray}\displaystyle & \displaystyle g(r)=\exp \left(-a\left(\frac{r}{R_{p}}\right)\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle g(r)=\exp \left(-\frac{1}{2a^{2}}\left(\frac{r}{R_{p}}\right)^{2}\right), & \displaystyle\end{eqnarray}$$ where 
$a$ is a constant that will scale the decay rate of 
$g(r)$ and allow 
$\unicode[STIX]{x1D70E}_{1/2}$ to be set. Since the volume integration within a DNS simulation domain must occur at finite radial distances 
$(r)$, a compact support must be imposed upon the integration kernel (
$g(r)$ is truncated after some threshold 
$r$). The compact support defines the union of the fluid volume (
$\unicode[STIX]{x1D734}_{y}$) and the particle volume. For the integral in (6.1) to converge, the compact support must be sufficiently large with respect to the kernel width (
$\unicode[STIX]{x1D70E}_{1/2}$) - i.e. the normalized integral on 
$\unicode[STIX]{x1D734}$ is close to unity
$$\begin{eqnarray}\int {\hat{g}}(r)\,\text{d}\unicode[STIX]{x1D734}\approx 1.\end{eqnarray}$$ For the integration kernels presented in (A 2)–(A 3), the kernel width 
$(\unicode[STIX]{x1D70E}_{1/2})$ and normalized integral (A 4) are given for varying values of 
$a$ and compact support 
$(\unicode[STIX]{x1D734})$; see tables 6–7. For 
$a=1$ in (A 2) (
$\unicode[STIX]{x1D70E}_{1/2}\approx 2.5R_{p}$), the integration kernel of Deen et al. (Reference Deen, Kriebitzsch, van der Hoef and Kuipers2012) is obtained and the normalized integral (A 4) is greater than 0.90 for 
$\unicode[STIX]{x1D734}>5R_{p}$ (see table 6). This result is in agreement with the findings of Deen et al. (Reference Deen, Kriebitzsch, van der Hoef and Kuipers2012) where it was reported that an averaging box of radius 
${\geqslant}5R_{p}$ led to convergence of the local heat transfer coefficients. However, convergence of the integral in (6.1) does not rectify the ambiguity associated with the inherent length scale of the integration kernel 
$(\unicode[STIX]{x1D70E}_{1/2})$ – i.e. an integration kernel with 
$\unicode[STIX]{x1D70E}_{1/2}=1.5R_{p}$ and 
$\unicode[STIX]{x1D734}=4R_{p}$ could alternatively be selected and while the volume averaging integral (6.1) has converged, the resulting local variables 
$(T_{Loc};U_{Loc})$ will differ from the case 
$\unicode[STIX]{x1D70E}_{1/2}=2.5R_{p}$ and 
$\unicode[STIX]{x1D734}=6R_{p}$ (also converged integral). When considering a fluidized bed, Anderson & Jackson (Reference Anderson and Jackson1967) argue that the dependence upon 
$g(r)$ and 
$\unicode[STIX]{x1D70E}_{1/2}$ becomes insignificant when 
$\unicode[STIX]{x1D70E}_{1/2}$ is large with respect to the particle–particle spacing and small with respect to the variations of the complete system. For dilute particle flows, this criterion would imply that there is not a unique averaging method since the particle spacing does not display a separation of scales from the macroscopic system. While dilute flows do not display a separation of scales in the particle spacing, they will experience a separation of scales between their local hydrodynamic environment and the macroscopic system – i.e. two particles in shear flow separated by a distance large enough that they do not impact each other.
Table 6. Exponential kernel attributes.
Table 7. Gaussian kernel attributes.
Figure 22. Nusselt numbers for unbounded flow past a sphere obtained from local variables computed with the kernel attributes given in table 6 (exponential kernel, equation (A 2)). The kernel width 
$(\unicode[STIX]{x1D70E}_{1/2})$ decreases from left to right. The open circle markers denote a support of 
$\unicode[STIX]{x1D734}=4R_{p}$ while the open diamond markers denote a support of 
$\unicode[STIX]{x1D734}=5R_{p}$ and the open square markers denote a support of 
$\unicode[STIX]{x1D734}=6R_{p}$. The solid red line denotes the free-stream correlation of Richter & Nikrityuk (Reference Richter and Nikrityuk2012) and the dashed blue line denotes the filtered correlation given in (C 1).
Figure 23. Drag coefficients for unbounded flow past a sphere obtained from local variables computed with the kernel attributes given in table 6. The kernel width 
$(\unicode[STIX]{x1D70E}_{1/2})$ decreases from left to right. The open circle markers denote a support of 
$\unicode[STIX]{x1D734}=4R_{p}$ while the open diamond markers denote a support of 
$\unicode[STIX]{x1D734}=5R_{p}$ and the open square markers denote a support of 
$\unicode[STIX]{x1D734}=6R_{p}$. The solid red line denotes the free-stream correlation of Haider & Levenspiel (Reference Haider and Levenspiel1989) and the dashed blue line denotes the filtered correlation given in (B 13).
Figure 24. Nusselt numbers for unbounded flow past a sphere obtained from local variables computed with the kernel attributes given in table 7 (Gaussian kernel, equation (A 3)). The kernel width 
$(\unicode[STIX]{x1D70E}_{1/2})$ decreases from left to right. The open diamond markers denote a support of 
$\unicode[STIX]{x1D734}=5R_{p}$. The solid red line denotes the free-stream correlation of Richter & Nikrityuk (Reference Richter and Nikrityuk2012) and the dashed blue line denotes the filtered correlation given in (C 1).
We apply the averaging kernels in (A 2)–(A 3) (tables 6–7) to data obtained from simulation of uniform flow past a static, cold particle and compare the results to free-stream correlations (Haider & Levenspiel Reference Haider and Levenspiel1989; Richter & Nikrityuk Reference Richter and Nikrityuk2012) (solid red lines) as well as filtered correlations (Ireland & Desjardins Reference Ireland and Desjardins2017) (dashed blue lines); see figures 22–25. As the kernel width is reduced, the averaging is more heavily weighted to the slow, cold fluid near the particle surface and the resulting 
$Nu$ and 
$C_{D}$ values undergo shifts away from the free-stream correlations of Haider & Levenspiel (Reference Haider and Levenspiel1989) and Richter & Nikrityuk (Reference Richter and Nikrityuk2012). The shifts away from the free-stream correlations are due to interphase transport that reduce the fluid velocity and temperature in the vicinity of the particle. Furthermore, the shifted data obtained with local variables are self-similar to the free-stream correlations and suggests that the local fluid environment contains sufficient information for correlating the drag and heat transfer, provided the averaging volume is large with respect to the particle boundary layer. Note, that for a given 
$\unicode[STIX]{x1D70E}_{1/2}$, the data in each panel of figures 22–23 converge to a trend as the compact support 
$(\unicode[STIX]{x1D734})$ of the integration kernel increases ((A 4) becomes satisfied as the upper integration bound tends to 
$\infty$). For similar kernel widths 
$(\unicode[STIX]{x1D70E}_{1/2})$, equations (A 2) and (A 3) yield very comparable results for 
$Nu$ and 
$C_{D}$; see figures 22 versus 24 and 23 versus 25. For integration kernels of finite width, an offset will be present between correlations that utilize local variables and correlations that utilize free-stream variables and the offset will become larger as 
$\unicode[STIX]{x1D70E}_{1/2}$ decreases. Therefore, a correlation derived with local variables will be specific to the 
$\unicode[STIX]{x1D70E}_{1/2}$ length scale, assuming 
$\unicode[STIX]{x1D734}$ is large enough to achieve convergence. Similarly, the simulation of flows via unresolved CFD-DEM involves a characteristic length scale associated with the fluid grid size 
$(\unicode[STIX]{x0394}x)$. For methods like those derived by Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013), the length scale is associated with the mollifier width (
$\unicode[STIX]{x1D70E}_{1/2}=3R_{p}$ in that work) since the source terms due to interphase transport are distributed among adjacent nodes according to the mollification kernel. For both classic CFD-DEM and the strategy of Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013), the length scale (
$\unicode[STIX]{x0394}x$ and 
$\unicode[STIX]{x1D70E}_{1/2}$, respectively) is often 
${\approx}2-3R_{p}$. For this reason, we employ (A 2) with 
$\unicode[STIX]{x1D70E}_{1/2}=2.5R_{p}$ and 
$\unicode[STIX]{x1D734}=5R_{p}$; which is an approximately converged integral that has been previous employed with DNS data (Deen et al. Reference Deen, Kriebitzsch, van der Hoef and Kuipers2012) and whose kernel width is consistent with unresolved CFD-DEM methods (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013).
Figure 25. Drag coefficients for unbounded flow past a sphere obtained from local variables computed with the kernel attributes given in table 7. The kernel width 
$(\unicode[STIX]{x1D70E}_{1/2})$ decreases from left to right. The open diamond markers denote a support of 
$\unicode[STIX]{x1D734}=5R_{p}$. The solid red line denotes the free-stream correlation of Haider & Levenspiel (Reference Haider and Levenspiel1989) and the dashed blue line denotes the filtered correlation given in (B 13).
Appendix B. Filtered drag correlations
B.1 Ireland & Desjardins (Reference Ireland and Desjardins2017)
Discrete particles interacting with a continuum fluid phase will experience an interphase exchange of momentum and thermal energy. The interphase exchange leads to source terms within unresolved CFD-DEM frameworks that effectively couple the phases. The source term due to drag will alter the fluid velocity in the vicinity of the particle, which is subsequently interpolated to calculate the particle drag force for the next time step. The disturbance to the fluid velocity, caused by the drag source term, leads to errors in the computed drag force for two-way coupled frameworks. Ireland & Desjardins (Reference Ireland and Desjardins2017) filter the solution for Stokes flow past a sphere to account for the fluid velocity disturbance introduced by the particle. For completeness, and to establish a procedure for appendix C, we re-derive the solution of Ireland & Desjardins (Reference Ireland and Desjardins2017) here.
The Stokes velocity field in spherical coordinates is given by (Happel & Brenner Reference Happel and Brenner1983)
$$\begin{eqnarray}\displaystyle & \displaystyle u_{r}=U_{\infty }\left[1-\frac{3}{2}\left(\frac{R_{p}}{r}\right)+\frac{1}{2}\left(\frac{R_{p}}{r}\right)^{3}\right]\cos (\unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{\unicode[STIX]{x1D703}}=U_{\infty }\left[-1+\frac{3}{4}\left(\frac{R_{p}}{r}\right)+\frac{1}{4}\left(\frac{R_{p}}{r}\right)^{3}\right]\sin (\unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$ which may be mapped to Cartesian coordinates via rotation matrices. Sequential rotations about the 
$\unicode[STIX]{x1D719}$ and 
$\unicode[STIX]{x1D703}$ axes lead to a total rotation matrix of
$$\begin{eqnarray}\unicode[STIX]{x1D64D}=\left[\begin{array}{@{}ccc@{}}\cos (\unicode[STIX]{x1D703}) & -\sin (\unicode[STIX]{x1D703}) & 0\\ \sin (\unicode[STIX]{x1D703})\cos (\unicode[STIX]{x1D719}) & \cos (\unicode[STIX]{x1D703})\cos (\unicode[STIX]{x1D719}) & -\sin (\unicode[STIX]{x1D719})\\ \sin (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D719}) & \cos (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D719}) & \cos (\unicode[STIX]{x1D719})\end{array}\right].\end{eqnarray}$$ Multiplication of 
$\unicode[STIX]{x1D64D}$ by (B 1)–(B 2) in vector format yields
$$\begin{eqnarray}\displaystyle & \displaystyle u_{x}=U_{\infty }\left[1-\frac{3}{4\hat{r}}\left(1+\cos ^{2}(\unicode[STIX]{x1D703})\right)-\frac{1}{4\hat{r}^{3}}(1-3\cos ^{2}(\unicode[STIX]{x1D703}))\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{y}=U_{\infty }\left[\sin (\unicode[STIX]{x1D703})\cos (\unicode[STIX]{x1D703})\cos (\unicode[STIX]{x1D719})\left(-\frac{3}{4\hat{r}}+\frac{3}{4\hat{r}^{3}}\right)\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{z}=U_{\infty }\left[\sin (\unicode[STIX]{x1D703})\cos (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D719})\left(-\frac{3}{4\hat{r}}+\frac{3}{4\hat{r}^{3}}\right)\right], & \displaystyle\end{eqnarray}$$ where 
$\hat{r}=r/R_{p}$. The filtered fluid volume fraction 
$(\unicode[STIX]{x1D700}_{g,Loc})$ is obtained by integrating the Gaussian kernel over the fluid volume
$$\begin{eqnarray}\unicode[STIX]{x1D700}_{g,Loc}\equiv \frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a^{3}}\int _{1}^{\infty }\exp \left(-\frac{1}{2a^{2}}\hat{r}^{2}\right)\hat{r}^{2}\,\text{d}\hat{r}=1-\text{erf}\left(\frac{1}{a\sqrt{2}}\right)+\frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a}\exp \left(-\frac{1}{2a^{2}}\right).\end{eqnarray}$$The phase-averaged Stokes velocity is obtained by applying the Gaussian filter to the fluid velocity in the x-direction
$$\begin{eqnarray}\unicode[STIX]{x1D700}_{g,Loc}u_{x,Loc}=\frac{\sqrt{2/\unicode[STIX]{x03C0}}}{2a^{3}}\int _{0}^{\unicode[STIX]{x03C0}}\int _{1}^{\infty }\exp \left(-\frac{1}{2a^{2}}\hat{r}^{2}\right)u_{x}\hat{r}^{2}\sin (\unicode[STIX]{x1D703})\,\text{d}\hat{r}\,\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$ where the 2 present in the denominator of (B 8) is a result of having integrated over 
$\unicode[STIX]{x1D719}$ but not 
$\unicode[STIX]{x1D703}$. Terms in (B 4) without 
$\unicode[STIX]{x1D703}$ dependence simplify to
$$\begin{eqnarray}U_{\infty }\frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a^{3}}\int _{1}^{\infty }\exp \left(-\frac{1}{2a^{2}}\hat{r}^{2}\right)\left[\hat{r}^{2}-\frac{3}{4}\hat{r}-\frac{1}{4\hat{r}}\right]\text{d}\hat{r},\end{eqnarray}$$ while terms with 
$\unicode[STIX]{x1D703}$ dependence 
$(\cos ^{2}(\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703}))$ are integrated to obtain
$$\begin{eqnarray}U_{\infty }\frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a^{3}}\int _{1}^{\infty }\exp \left(-\frac{1}{2a^{2}}\hat{r}^{2}\right)\left[-\frac{1}{4}\hat{r}+\frac{1}{4\hat{r}}\right]\text{d}\hat{r}.\end{eqnarray}$$ The first term in brackets in (B 9) is exactly 
$\unicode[STIX]{x1D700}_{g}U_{\infty }$. Combining 
$\hat{r}$ and 
$1/\hat{r}$ terms in (B 9)–(B 10), we arrive at
$$\begin{eqnarray}U_{\infty }\frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a^{3}}\int _{1}^{\infty }\exp \left(-\frac{1}{2a^{2}}\hat{r}^{2}\right)(-\hat{r})\,\text{d}\hat{r}=-U_{\infty }\frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a}\exp \left(-\frac{1}{2a^{2}}\right).\end{eqnarray}$$Dividing by the solids’ volume fraction (B 7), we obtain the solution reported by Ireland & Desjardins (Reference Ireland and Desjardins2017)
$$\begin{eqnarray}u_{x,Loc}=\left[1-\frac{{\displaystyle \frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a}}\exp \left(-{\displaystyle \frac{1}{2a^{2}}}\right)}{1-\text{erf}\left({\displaystyle \frac{1}{a\sqrt{2}}}\right)+{\displaystyle \frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a}}\exp \left(-{\displaystyle \frac{1}{2a^{2}}}\right)}\right]U_{\infty },\end{eqnarray}$$ where 
$u_{x,Loc}$ is the local fluid velocity for Stokes flow around a sphere, and 
$a$ is the standard deviation in (A 3). We note that the local fluid velocity is the difference between the undisturbed fluid velocity 
$(U_{\infty })$ and the velocity correction 
$\unicode[STIX]{x1D701}_{u_{f}}$ (equation 36 in Ireland & Desjardins (Reference Ireland and Desjardins2017)).
The local fluid velocity given in (B 12) may be utilized to modify existing drag correlations (Haider & Levenspiel Reference Haider and Levenspiel1989) for use with local (disturbed) fluid velocities obtained in unresolved CFD-DEM methods. Specifically, equation (B 12) has the general form 
$u_{x,Loc}=C_{H}U_{\infty }$ where 
$C_{H}$ (term in brackets in (B 12)) is a scaling factor that depends upon the filter width and is less than unity. The drag force is given by 
$F_{D}=\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}D_{p}^{2}u^{2}C_{D}$. Substituting in 
$u_{x,Loc}=C_{H}U_{\infty }$ for 
$u$ illustrates that 
$C_{D}$ must be scaled by 
$1/C_{H}^{2}$ in order to obtain the same drag force with an undisturbed velocity (
$U_{\infty }$). The drag correlation of Haider & Levenspiel (Reference Haider and Levenspiel1989) may be modified as follows
$$\begin{eqnarray}C_{D}=\frac{24}{C_{H}Re_{Part,Loc}}(1+0.18C_{H}^{-0.65}Re_{Part,Loc}^{0.65}),\end{eqnarray}$$ where 
$C_{H}$ is the bracketed term in (B 12) evaluated for 
$a=1.625$ and 
$Re_{Part,Loc}=C_{H}Re_{Part}$.
B.2 Municchi et al. (Reference Municchi, Goniva and Radl2016)
Additionally, Municchi et al. (Reference Municchi, Goniva and Radl2016) considered a similar procedure with a top-hat kernel and obtained
$$\begin{eqnarray}u_{x,Loc}=\left[\frac{2f^{2}-f-1}{2(f^{2}+f+1)}\right]U_{\infty },\end{eqnarray}$$ where 
$f=\unicode[STIX]{x1D734}/R_{p}$ is the dimensionless ratio of the filter support to the particle radius (5 in the present study). The same procedure as described above is utilized to arrive at a modified drag law but with 
$C_{H}$ equal to the bracketed term in (B 14).
Appendix C. Filtered heat transfer correlation
When considering heat transfer, the Nusselt correlation must be modified to account for the disturbed fluid temperature and the disturbed Reynolds number. To obtain a local fluid temperature, we begin by filtering the analytical solution for radial diffusion from a sphere. Subsequently, we consider thermal Stokes flow past a sphere (Acrivos & Taylor Reference Acrivos and Taylor1962). In contrast to the hydrodynamics, a simple solution for this case is not readily available and the asymptotic method employed by Acrivos & Taylor (Reference Acrivos and Taylor1962) yields ‘inner’ and ‘outer’ solutions that are valid near the particle surface and far away, respectively.
The rate of interphase heat transfer is given by 
$\dot{Q}=h(T_{\infty }-T_{p})$, which simplifies to 
$\dot{Q}=hT_{\infty }$ for zero particle temperature. Employing a local fluid temperature (
$T_{Loc}=C_{T}T_{\infty }$) implies that the Nusselt number should be scaled by 
$1/C_{T}$ to conserve the heat rate. However, the Nusselt number also depends upon the particle Reynolds number, which also experiences a disturbance. Modifying the Nusselt correlation of Richter & Nikrityuk (Reference Richter and Nikrityuk2012) we arrive at
$$\begin{eqnarray}Nu=\frac{1}{C_{T}}[1.76+0.55C_{H}^{-1/2}Re_{Part,Loc}^{1/2}Pr^{1/3}+0.014C_{H}^{-2/3}Re_{Part,Loc}^{2/3}Pr^{1/3}],\end{eqnarray}$$ where 
$C_{T}$ is the thermal analogue of 
$C_{H}$ and is closed in the following subsections.
C.1 Diffusion filtering
For the case of a cold particle placed in a quiescent bath of hot fluid, the transient diffusion equation in the radial coordinate governs the evolution of the fluid temperature. The solution to this problem may be found via combination of variables (Lattanzi et al. Reference Lattanzi, Yin and Hrenya2019b)
$$\begin{eqnarray}T=\left[1-\frac{1}{\hat{r}}\text{erfc}\left(\frac{\hat{r}-1}{\sqrt{4\unicode[STIX]{x1D6FC}t}}\right)\right]T_{\infty },\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FC}$ is the thermal diffusivity of the fluid. For long time scales (steady state), the erfc term will tend to unity and the solution simplifies to 
$1-1/\hat{r}$. Filtering of a constant and 
$1/\hat{r}$ term will yield
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a^{3}}\int _{1}^{\infty }\exp \left(-\frac{1}{2a^{2}}\hat{r}^{2}\right)\hat{r}^{2}\,\text{d}\hat{r}=\unicode[STIX]{x1D700}_{g,Loc}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a^{3}}\int _{1}^{\infty }\exp \left(-\frac{1}{2a^{2}}\hat{r}^{2}\right)\hat{r}\,\text{d}\hat{r}=\frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a}\exp \left(-\frac{1}{2a^{2}}\right). & \displaystyle\end{eqnarray}$$ Collecting terms and dividing by 
$\unicode[STIX]{x1D700}_{g,Loc}$ we obtain
$$\begin{eqnarray}T_{Loc}=\left[1-\frac{{\displaystyle \frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a}}\exp \left(-{\displaystyle \frac{1}{2a^{2}}}\right)}{1-\text{erf}\left({\displaystyle \frac{1}{a\sqrt{2}}}\right)+{\displaystyle \frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a}}\exp \left(-{\displaystyle \frac{1}{2a^{2}}}\right)}\right]T_{\infty },\end{eqnarray}$$ which is exactly the solution obtained by Ireland & Desjardins (Reference Ireland and Desjardins2017) for the hydrodynamics (B 12). Therefore, we compute 
$C_{T}$ in the same manner as 
$C_{H}$ for both the Ireland & Desjardins (Reference Ireland and Desjardins2017) correlation and the Municchi et al. (Reference Municchi, Goniva and Radl2016) correlation.
C.2 Thermal Stokes filtering
To assess the first-order effects arising from finite fluid flow, we first consider the inner expansion of Acrivos & Taylor (Reference Acrivos and Taylor1962), which holds near the surface of the particle. The truncated inner solution is given by
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}T_{In}=\left[1-{\displaystyle \frac{1}{\hat{r}}}+\left({\displaystyle \frac{Pe}{4}}+{\displaystyle \frac{Pe^{3}}{16}}\ln \left({\displaystyle \frac{Pe}{2}}\right)\right)T_{1}+{\displaystyle \frac{Pe^{2}}{8}}\ln \left({\displaystyle \frac{Pe}{2}}\right)\left(1-{\displaystyle \frac{1}{\hat{r}}}\right)\right]T_{\infty },\\[12.0pt] T_{1}=\left(1-{\displaystyle \frac{1}{\hat{r}}}\right)-\left(1-{\displaystyle \frac{3}{2\hat{r}}}+{\displaystyle \frac{3}{4\hat{r}^{2}}}-{\displaystyle \frac{1}{4\hat{r}^{3}}}\right)\cos (\unicode[STIX]{x1D703}).\end{array}\right\}\end{eqnarray}$$ We may again utilize (C 3)–(C 4) to filter the constant and 
$1/\hat{r}$ terms. Due to the orthogonality of 
$\sin (\unicode[STIX]{x1D703})$ and 
$\cos (\unicode[STIX]{x1D703})$, the second set of terms in 
$T_{1}$ will be zero when the 
$\unicode[STIX]{x1D703}$ direction is integrated, and thus will not contribute to 
$T_{Loc}$. Combining terms we obtain a solution of
$$\begin{eqnarray}\displaystyle T_{Loc,In} & = & \displaystyle \left[1-\frac{{\displaystyle \frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a}}\exp \left(-{\displaystyle \frac{1}{2a^{2}}}\right)}{1-\text{erf}\left({\displaystyle \frac{1}{a\sqrt{2}}}\right)+{\displaystyle \frac{\sqrt{2/\unicode[STIX]{x03C0}}}{a}}\exp \left(-{\displaystyle \frac{1}{2a^{2}}}\right)}\right]\nonumber\\ \displaystyle & & \displaystyle \times \,\left[1+\frac{Pe}{4}+\frac{Pe^{2}}{8}\ln \left(\frac{Pe}{2}\right)+\frac{Pe^{3}}{16}\ln \left(\frac{Pe}{2}\right)\right]T_{\infty }.\end{eqnarray}$$ Filtering of the inner solution is valid for smaller Péclet numbers 
$(O(10^{-1}))$ where the boundary layer is not confined to a thin region near the particle surface. By contrast, filtering of the outer solution is valid for larger Péclet numbers 
$(O(10^{1}))$ since the near-surface region will not contribute significantly to the integral. The outer solution to first order is given by
$$\begin{eqnarray}T_{Out}=\left[1-\frac{2}{Pe~\hat{r}}\exp \left(\frac{Pe~\hat{r}}{4}(\cos (\unicode[STIX]{x1D703})-1)\right)\right]T_{\infty }.\end{eqnarray}$$ The outer solution is numerically filtered and thus no closed form solution is given for 
$T_{Loc,Out}$. For the case of thermal Stokes filtering, we consider a local velocity (
$C_{H}$) given by (B 12) with a local temperature (
$C_{T}$) given by (C 7) for 
$Pe<1$ and numerical integration of (C 8) for 
$Pe>10$.
