Gas-Solid Heat Transfer Computation from Particle-Resolved Direct Numerical Simulations

Abstract

Particle-Resolved simulations (PR-DNS) have been conducted using a second order implicit Viscous Penalty Method (VPM) to study the heat transfer between a set of particles and an incompressible carrier fluid. A Lagrange extrapolation coupled to a Taylor interpolation of a high order is utilized to the accurate estimate of heat transfer coefficients on an isolated sphere, a fixed Faced-Centered Cubic array of spheres, and a random pack of spheres. The simulated heat transfer coefficients are compared with success to various existing Nusselt laws of the literature.

1. Introduction

Gas-solid flows are widely encountered in nature, for example in volcanic eruptions [1,2]. They also participate in heat transfers in many industrial processes such as petroleum refining, blast furnaces, or chemical looping combustion [3,4,5]. Due to the scale disparity between these applications and the particle size at the industrial scale, the CFD simulations of such applications are based on statistical approaches where the average interphase transfer of energy between the fluid and the particles needs to be modeled.

Theoretical and Experimental approaches have been widely employed to model these unclosed terms with significant limitations. Indeed, theoretical results are limited to Stokes flows or moderate Reynolds number regimes [6,7]. Also, experimental measurements showed huge differences in Nusselt number at high volume fractions because of limited optical access [8]. This motivates the community to consider and develop Particle-Resolved Direct Numerical Simulations (PR-DNS), where all scales are resolved, to directly compute the fluid-particle interaction and the associated heat transfer closure laws. In this regard, The present work propose a method to extract heat transfer from PR-DNS carried out using Viscous Penalty Method [9].

In the last two decades, several research teams have conducted numerical studies to characterize, understand, and model the energy interphase exchange in gas-solid flows for numerous configurations of particulate motions. Among them, we can cite Massol [10] who studied heat transfers in a fixed array of monodisperse spheres using a body-fitted three-dimensional PR-DNS (AVBP code [11]). This configuration will be studied in this work and its results will be compared to ours. Deen et al. [12] and Tavassoli et al. [13] used Immersed Boundary Methods (IBM) for three-dimensional PR-DNS with inflow and outflow boundary conditions to compute the gas-solid Nusselt number that is successfully compared with Gunn’s correlation [14]. Note that Gunn [14] proposed his correlation based on experimental results. Tenneti et al. [15] utilized the so-called PUReIBM with periodic boundary conditions and proposed a correlation from their results for various flow configurations. Later, Deen et al. [16] proposed a correlation based on Tavassoli et al.  [13] and Tenneti et al. [15] results. Tenneti et al.  [15] simulations were carried out for a limited range of Reynolds numbers and solid volume fractions. As a consequence, Sun [17] took over this work and extended it to larger values of Reynolds numbers and solid volume fractions, and proposed a correlation that fits their PR-DNS data.

Our contribution is to propose an accurate heat flux estimate for PR-DNS simulations using Viscous Penalty Method (VPM) [9]. A previous work using VPM was conducted using an indirect method to compute the heat transfer in an assembly of spheres [18]. In this work, on the other hand, a direct approach is proposed to extract the heat flux between each particle and the surrounding fluid. This is an extension of the developments conducted to compute the hydrodynamic force in a previous study [19]. Furthermore, the proposed method manages to avoid the flow contaminated by the VPM method [19] in the Eulerian cells cut by the fluid/solid particle interface and therefore extracts the accurate heat flux needed to compute the Nusselt number.

The article is structured as follows. A presentation of the models and numerical methods is first proposed in Section 2. In Section 3, validations for flows interacting with isolated spherical particles at various Reynolds numbers are discussed. Simulations of a uniform flow past a Face-Centered cubic array of spheres are presented in the Section 4. Section 5 is devoted to simulations and validations of flows through random arrangements of monodispersed spheres. Finally, conclusions and perspectives are drawn.

2. Numerical Methodology

2.1. Viscous Penalty Method

The motion equations of incompressible anisothermal two-phase flows, involving a carrier fluid and a solid particle phase, are based on one-fluid model [20,21], as explained in [9,22]. They are given by:

$$\begin{array}{c}\hfill \nabla ·\mathbf{u}=0\end{array}$$

(1)

$$\begin{array}{c}\hfill \rho \left(\frac{\partial \mathbf{u}}{\partial t}+\left(\mathbf{u}·\nabla \right)\mathbf{u}\right)=-\nabla p+\rho \mathbf{g}+\nabla ·\left[\mu (\nabla \mathbf{u}+{\nabla }^t\mathbf{u})\right]+{\mathbf{F}}_m\end{array}$$

(2)

$$\begin{array}{c}\hfill \rho C_v\left(\frac{\partial T}{\partial t}+\mathbf{u}·\nabla T\right)=\nabla ·\left(k_f\nabla T\right)\end{array}$$

(3)

where ${\mathbf{F}}_m$ is a source term, $\mathbf{u}$ is the velocity, p the pressure, T the temperature, t the time, $\mathbf{g}$ the gravity vector, $\rho$, $\mu$, $C_v$ and $k_f$ respectively the density, the viscosity, the specific heat and the thermal conductivity of the equivalent fluid. Indeed, VPM belongs to the class of fictitious domain methods [23,24], where a fixed mesh for the whole domain (fluid+particles) is used to solve the dynamic of both phases (see Figure 1). Therefore, “equivalent fluid” has the fluid or solid thermophysical properties in the Eulerian cell fully fluid or solid respectively, and is a mix of both in the Eulerian mesh containing the fluid/particle interface. In that sens, specific arithmetic ${\displaystyle \left(\rho =C{\rho }_s+(1-C){\rho }_f,C_v=C{C}_{v_s}+(1-C)C_{v_f}\right)}$, harmonic ${\displaystyle \left(\mu =\frac{{\mu }_s{\mu }_f}{C{\mu }_f+(1-C){\mu }_s}\right.}$, ${\displaystyle \left.k_f=\frac{k_{f_s}{k}_{f_f}}{C{k}_{f_f}+(1-C)k_{f_s}}\right)}$ laws are used to define them according to a phase function C [22] (the quantities with the subscript s, f refers to those in the solid, fluid respectively). As for the phase function C (also called Color function), it is directly obtained by projecting the shape of the particles on the Eulerian mesh [9,19] (As illustrated in Figure 1). Indeed, at each time iteration, this projection is conducted after the particles positions are updated by solving ${\displaystyle \frac{d{\mathbf{x}}_i}{dt}={\mathbf{V}}_i}$, where $V_i$ is the velocity of the $i$ particle deduced from the Eulerian velocities inside the particle (${\displaystyle {\mathbf{x}}_i}$ being its position). This procedure is called VOF-Lag approach in contrast with the classical VOF method [25] where the equation ${\displaystyle \frac{\partial C}{\partial t}+\mathbf{u}·\nabla C=0}$ is solved to update the color function C.

The solid behavior in the zones occupied by the particles (where $C=1$), is ensured by imposing large values of viscosity in the Eulerian cells inside the particles. Indeed, the strain tensor in (2) tends toward 0 when ${\displaystyle \mu \to +\infty }$ and based on the equivalence [9] ${\displaystyle \left(\forall \mathbf{P}\in {\Omega }_i,\nabla \mathbf{u}+{\nabla }^T\mathbf{u}=0\iff \mathbf{V}\left(\mathbf{P}\right)={\mathbf{V}}_i+{\mathbf{O}}_i\mathbf{P}\wedge {\omega }_i\right)}$ a solid behaviour is, in this way, imposed in the zones inside the particles. Moreover, a specific model [26], based on a decomposition of the viscous stress tensor is implemented to implicitly impose not only the solid behavior (via the viscosity penalization) in the particles but also the coupling between fluid and solid motion.

To satisfy the fluid and solid constraints along with the velocity and pressure coupling, an augmented Lagrangian method is applied to (1), (2) as follows (Algorithm 1) [9]:

Algorithm 1: Augmented Lagrangian (AL) method [27]

for each iteration n do     ${\mathbf{u}}^{\ast ,0}={\mathbf{u}}^n$ and $p^{\ast ,0}=p^n$     $m=0$     while $\left|\right|\nabla ·{\mathbf{u}}^{\ast ,\mathbf{m}}\left|\right|>{ϵ}_{AL}$ do         $m=m+1$         solve: ${\displaystyle \rho \left({\displaystyle \frac{{\mathbf{u}}^{\ast ,m}-{\mathbf{u}}^{\ast ,0}}{\Delta t}+{\mathbf{u}}^{\ast ,m-1}·\nabla {\mathbf{u}}^{\ast ,m}}\right)-\nabla (r\nabla ·{\mathbf{u}}^{\ast ,m})=-\nabla p^{\ast ,m-1}+\rho \mathbf{g}+\nabla ·\left[\mu (\nabla {\mathbf{u}}^{\ast ,m}+{\nabla }^T{\mathbf{u}}^{\ast ,m})\right]}$         ${\displaystyle p^{\ast ,m}=p^{\ast ,m-1}-r\nabla ·{\mathbf{u}}^{\ast ,m}}$     end while end for

The mass, momentum and energy equations are discretized using an implicit finite volumes on structured staggered meshes (see Figure 1). The time derivative is approximated with a first order Euler scheme while the inertial, viscous and augmented Lagrangian terms are discretized with second-order centered schemes. All fluxes are written at time $(n+1)\Delta t$, except the non-linear inertial term that is linearized with a second order Adams-Bashforth scheme as follows

$$\mathbf{u}·\nabla \mathbf{u}\approx \left(2{\mathbf{u}}^n-{\mathbf{u}}^{n-1}\right)·\nabla {\mathbf{u}}^{n+1}$$

(4)

The obtained linear system is solved with a BiCGSTAB II iterative solver [28], preconditioned under a Modified and Incomplete LU approach [29] to speed-up the convergence of the solver. All the code is working on massively parallel computers by using MPI [9].

Note that the particles in this work are motionless. Numerically, their velocity is set to zero by imposing the velocity of the Eulerian cells near their centroids to zero by means of a linear term, ${10}^{40}\mathbf{u}$, added to the Equation (2), and the VPM propagates the zero velocity in the whole solid medium. The particles temperature $T_s$ being also constant, they are however imposed in the whole solid medium by means of a linear term, ${\displaystyle {10}^{40}(T-T_s)}$, added to the Equation (3).

2.2. Heat Transfer Rate Computation

The goal of this work is to extract from the temperature field, solution of the energy Equation (3), the Nusselt number given by:

$$Nu=\frac{H_{p}d}{k_f}$$

(5)

where d is the particle diameter and $H_p$ is the heat transfer rate given by:

$$H_p=\frac{Q_p}{\left(T_s-T_f\right)\pi d^2}$$

(6)

Here, $T_s$ and $T_f$ are the temperatures of the particle and the surrounding fluid respectively, and $Q_p$ is the heat flux across the particle surface S. It is given by:

$$Q_p={\oint }_S-k_f\nabla T·\mathbf{n}dS$$

(7)

The computation of the heat flux consists in discretizing S on a set of N elements called Lagrangian mesh (linear segments in 2D and triangles in 3D, see Figure 2), such that:

$$Q_p\approx \sum _l^N-k_f{\left(\nabla T\right)}_l·{\mathbf{n}}_l\Delta S_l$$

(8)

where the ongoing normal to the lth element ${\mathbf{n}}_l$ and its area $\Delta S_l$ are deduced from its nodes coordinates. However, the temperature T being solved on the Eulerian mesh, the temperature gradient ${\left(\nabla T\right)}_l$ is not known at the lth element center $C_l$. To deal with this problem, we have implemented the same approach that the one we have used to obtain the stress tensor ${\overline{\overline{\sigma }}}_l$ components on the Lagrangian mesh in the hydrodynamic force computation ${\displaystyle \mathbf{F}\approx \sum _l^N{\overline{\overline{\sigma }}}_l·{\mathbf{n}}_l\Delta S_l}$. This was detailed in [19]. In this previous work, we have observed that the velocity and the pressure values were inaccurate in the cells containing the interface. Therefore, we had to extrapolate the stress tensor from the area far from the particle to the particle surface. After a numerical study of all the extrapolation parameters, we have reached the conclusion that the third-order Lagrange extrapolation (whose notations are illustrated in Figure 2) is the best compromise between accuracy and computational cost [19]. Similarly, as for force coefficients, we have the temperature gradient at the center of each Lagrangian mesh element as follows:

$$\left(\nabla T\right)\left({\mathbf{C}}_l\right)=\sum _{i=1}^3\left(\nabla T\right)\left({\mathbf{P}}_i\right)L_i\left({\mathbf{C}}_l\right),\mathrm{where}{L}_i\left({\mathbf{C}}_l\right)=\prod _{j\ne i}^3\frac{|{\mathbf{C}}_l-{\mathbf{P}}_j|}{|{\mathbf{P}}_i-{\mathbf{P}}_j|}$$

(9)

Given that the temperature gradient is computed on the Eulerian mesh and that the extrapolation points ${\mathbf{P}}_1$, ${\mathbf{P}}_2$ and ${\mathbf{P}}_3$ are constructed in the normal direction to the particle surface (see Figure 2), the temperature gradient at these points are not known. Once again, we relayed on the work [19] to interpolate the temperature gradient from the Eulerian mesh to the extrapolation points ${\mathbf{P}}_k$ ${\displaystyle (k=1,2,3)}$ using the third-order Taylor interpolation:

$$\begin{array}{cc}\hfill f\left({\mathbf{P}}_k\right)=& f\left(\mathbf{E}\right)+\frac{\partial f}{\partial x_i}\left(\mathbf{E}\right)({\left({\mathbf{C}}_l\right)}_i-{\mathbf{E}}_i)+\frac{1}{2}\frac{{\partial }^{2}f}{\partial x_i^2}\left(\mathbf{E}\right){({\left({\mathbf{C}}_l\right)}_i-{\mathbf{E}}_i)}^2\hfill \\ \hfill & +\frac{{\partial }^{2}f}{\partial x_i\partial x_j}\left(\mathbf{E}\right)({\left({\mathbf{C}}_l\right)}_i-{\mathbf{E}}_i)({\left({\mathbf{C}}_l\right)}_j-{\mathbf{E}}_j)+\mathrm{O}\left(\right||{\mathbf{C}}_l-{\mathbf{E}\left|\right|}^3)\hfill \end{array}$$

(10)

where $\mathbf{E}$ is the nearest Eulerian point to ${\mathbf{P}}_k$.

Another pertinent value to consider when studying the distribution of heat transfers on the surface of the particle, is the local Nusselt number given in the spherical system (see Figure 3) by:

$$N{u}_{loc}\left(\theta \right)=\frac{-\nabla T·\mathbf{n}d}{\left(T_s-T_f\right)}$$

(11)

where for each Lagrangian mesh element center (${\displaystyle {\mathbf{C}}_l=(x,y,z)}$), the azimuthal angle $\theta$, and the polar angle $\beta$ are given by:

$$\theta =(\pi -arctan2(\sqrt{y^2+z^2},x))\frac{180}{\pi },\beta =arctan2(z,y)\frac{180}{\pi }$$

(12)

3. Convective Heat Transfer Forced by a Uniform Flow around a Stationary Sphere

The first case considered to validate heat transfer computation is the convection forced by a uniform flow past a hot sphere, illustrated in Figure 4. The computational domain lengths are $L_x=16d$ and $L_y=L_z=8d$. The Eulerian mesh refinement is constant in a box of extension $\left[\left(2d,3d,3d\right);\left(6d,5d,5d\right)\right]$ centered around the particle position. Outside this box, the Eulerian mesh is exponentially coarsen from the box to the boundaries of the simulation domain. This case was previously used to validate the drag force computation [19]. In the simulations carried out for this case as well as for all the cases studied in this work, we have imposed that the hydrodynamic boundary layer contains five Eulerian cells. Therefore, in the box surrounding the particle, the minimum cell size is ${\displaystyle \Delta x=\frac{d}{5\sqrt{Re}}}$. This mesh refinement ensures that the hydrodynamic boundary layer contains a sufficient number of extrapolation points to be accurate on force calculations at the particle surface. We have also chosen a Prandtl number $Pr=0.72$, which means that the thermal boundary layer is larger than the hydrodynamic boundary layer. In this way, we ensure that enough extrapolation points are available in the thermal boundary layer as well (see Figure 2). The time step is given by ${\displaystyle \Delta t=0.5min\left(\frac{\Delta x}{U_{\infty }},\frac{\Delta x^2}{k_f/C_v}\right)}$.

3.1. Effect of the Extrapolation Distance $\left(\delta \right)$

One of the extrapolation parameters studied in the hydrodynamic force computation work [19] was the distance $\delta$ between the Lagrangian mesh and the first extrapolation point ${\mathbf{P}}_1$ (see Figure 2). Indeed, during that study [19], inaccurate pressure and velocity values were observed in the cells cut by the interface, therefore we had to extrapolate the stress tensor from the fluid area far from the particle to the surface. As illustrated in (Figure 5 left), the accurate drag force is reached at ${\displaystyle \delta =\Delta x}$.

The temperature being a function of the velocity (Equation (3)), inaccurate values of the temperature gradient values are also expected in the Eulerian cells containing the interface. Thus, we have conducted the same study for the Nusselt number computation to locate the extrapolation area. But unlike the drag force (Figure 5 left), the Nusselt number computation does not seem to depend on the distance $\delta$ (Figure 5 right). Facing this unexpected result, we have compared the distribution of the local Nusselt number to Massol [10] results at $Re=100$ where the flow is symmetric with respect to its direction, as illustrated in (Figure 4 right). The local Nusselt number computed for $\delta =0$ does not reflect the flow symmetry as illustrated in (Figure 6a), and the error between our result $\left(N{u}_{loc}\right)$ and Massol’s $\left(N{u}_{locMassol}\right)$ given by ${\displaystyle \frac{{\sum }_{i=1}^N|{\left(N{u}_{loc}\right)}_i-{\left(N{u}_{locMassol}\right)}_i|}{{\sum }_{i=1}^N\left|{\left(N{u}_{locMassol}\right)}_i\right|}}$ is about $25\%$. This error decreases as the distance $\delta$ increases until being lower than $10\%$ for $\delta \ge 1$ (Figure 6c). It is worth noting that for this distance, the local Nusselt number distribution more accurately reflects the expected symmetry of the flow as it satisfies almost the same distribution for all the polar angle $\left(\beta \right)$ plans (Figure 6b). So, as for the hydrodynamic force computation [19] and in the rest of this work, the extrapolation distance will always be $\delta =1$.

3.2. Result on the Nusselt Coefficient

The numerical parameters of the heat transfer on the particle surface being set up, the global Nusselt number for a uniform flow past a hot sphere at different Reynolds number ${\displaystyle (Re=10,20,40,60,}$ ${\displaystyle 80,100,150,200,250,290)}$ is computed using the presented method and compared to the following correlations:

• Ranz and Marshall [31]:

  $$Nu=2.+0.6R{e}^{0.5}P{r}^{\frac{1}{3}}$$

  (13)
• Feng and Michaelides [32]:

  $$Nu=0.992+{\left(RePr\right)}^{\frac{1}{3}}+0.1R{e}^{\frac{2}{3}}P{r}^{\frac{1}{3}}$$

  (14)
• Whitaker [33]:

  $$Nu=2.+\left(0.4R{e}^{0.5}+0.06R{e}^{\frac{2}{3}}\right)P{r}^{0.4}$$

  (15)

Figure 7 provides the comparisons of our results to literature references. It can be observed that our simulated Nusselt values belong to the same value ensemble as literature works for all considered Reynolds numbers. If we go into more detail, the Nusselt numbers simulated with VPM show a better agreement with Feng and Michaelides and Whitaker correlations than with Ranz and Marshall.

4. Face-Centered Cubic Periodic Arrangement of Spheres

The second case studied in this work is a uniform flow past a Face-Centered Cubic (FCC) array of hot spheres which consists in a cube where three spheres are placed on the faces centers, and one sphere is located on the vertices with periodic boundary conditions, as illustrated in Figure 8. The domain length L is deduced from the solid volume fraction ${\alpha }_d$: ${\displaystyle L={\left(\frac{2}{3}\frac{\pi }{{\alpha }_d}\right)}^{1/3}d}$. Given a Reynolds number ${\displaystyle Re=\frac{\rho d{\mathbf{u}}_d}{\mu }}$, a desired mean fluid velocity $\langle {\mathbf{u}}_f\rangle ={\mathbf{u}}_d$ is enforced via the source term

$${\mathbf{F}}_m=-\rho \frac{\langle {\mathbf{u}}_f\rangle -{\mathbf{u}}_d}{\Delta t}+\frac{{\sum }_i^{N_p}{\mathbf{F}}_i^D}{V_f}$$

(16)

inserted in the momentum conservation Equation (2). ${\mathbf{F}}_i^D$ is the drag force computed over the ith particle using the method described in [19] and $V_f$ is the volume occupied by the fluid. ${\mathbf{F}}_m$ is thus adjusted until the stationary state is reached when ${\displaystyle \langle {\mathbf{u}}_f\rangle ={\mathbf{u}}_d}$. As for the uniform flow past an isolated sphere, the grid resolution is fixed by imposing five Cartesian cells in the boundary layer: ${\displaystyle \Delta x=\frac{d}{5\sqrt{Re}}}$, and the time step is given by ${\displaystyle \Delta t=0.5min\left(\frac{\Delta x}{U_{\infty }},\frac{\Delta x^2}{k_f/C_v}\right)}$.

This case was previously used to validate the hydrodynamic force computation [19]. In that previous work [19], the hydrodynamic force was computed when the steady-state of the flow was reached. This steady state obtained in an isothermal context is used in the present simulations as the initial condition for the momentum equation ‘Equation (2)’. In addition, the initial temperature conditions for the energy equation ‘Equation (3)’ are ${\displaystyle T_s=1,T_f=0.5}$, respectively the sphere and the fluid temperature. In this problem, the fluid heats up until its temperature $<T{>}_f$ (averaged on the fluid domain) reaches the sphere temperature (see Figure 9 left). In the meantime and after an initialization time of the temperature gradient around the spheres, a balance is reached between the heat flux and the difference between the fluid and the particle temperature. This results in the apparition of a plateau in the global Nusselt coefficient temporal evolution as illustrated in Figure 9 right. Therefore, the Nusselt number for a uniform flow past a FCC array of spheres is the value at this stabilized state.

Note that the global and local Nusselt coefficients for a uniform flow past an array of $N_p$ particles are an average of those computed on each sphere such as:

$$\langle Nu\rangle =\frac{1}{N_p}\sum _i^{N_p}N{u}_i\langle N{u}_{loc}\rangle =\frac{1}{N_p}\sum _i^{N_p}{\left(N{u}_{loc}\right)}_i$$

(17)

where ${\displaystyle N_p=4}$ for a FCC array.

4.1. Flow Analysis

Using the local Nusselt coefficient distribution presented in Figure 10, two regimes can be observed, i.e., attached (illustrated in Figure 8 left) and separated (Figure 8 right) flows, that govern the uniform flow past a FCC array of spheres, depending on the Reynolds number and the solid volume fraction ${\alpha }_d$:

4.1.1. Attached Flows

For low Reynolds numbers, the local Nusselt number is having the same behavior as for an isolated sphere (see Figure 10 left). The increase of its value is due to the presence of the other spheres (blocking effect) that speed up the flow and therefore flatten the boundary layer which increases the temperature gradient. We can also observe that the presence of surrounding particles breaks the symmetry of the flow. Indeed, (Figure 10 left) shows two distinct distributions of the local Nusselt number for two different plans ${\displaystyle (\beta =0)}$ and ${\displaystyle (\beta =45)}$, and these results are in good agreement with Massol body fitted simulations [10].

4.1.2. Separated Flows Downstream of the Spheres

By increasing the Reynolds number, a separated flow appears downstream of the spheres together with a fountain effect upstream of the spheres (see Figure 8 right). These phenomena can be observed also in the local Nusselt number distribution (Figure 10 right). Indeed, the higher value of the Nusselt number is no more at the upstream stagnation point ${\displaystyle \theta =0}$, due to the fountain effect that induces a low heat flux. With increasing ${\displaystyle \theta }$, i.e., leaving the fountain zone, the Nusselt number increases as the boundary layer shrinks due to the fluid acceleration in this zone. It then decreases until the separation point, where the fluid temperature is the highest. Finally, the Nusselt increases again in the recirculating zone behind the sphere.

4.2. Global Nusselt Coefficient

The Nusselt number of a uniform flow past an FCC array of spheres normalized by the isolated sphere Nusselt number (given by Ranz and Marshall [31]) ${\displaystyle N{u}_s=2.+0.6R{e}^{0.5}P{r}^{\frac{1}{3}}}$, noted ${\displaystyle \frac{\langle Nu\rangle }{N{u}_s}}$, is compared to existing literature results. This comparison, illustrated in Figure 11, shows a nice match with Massol [10] results at different Reynolds number and for a large range of solid volume fraction ${\displaystyle (0.05\le {\alpha }_d\le 0.6)}$. A good agreement is also observed with Gun [14] correlation for low solid volume fractions ${\displaystyle ({\alpha }_d\le 0.15)}$, which means that for this range of ${\alpha }_d$, the FCC configuration is a good approximation of a uniform flow past a random assembly of spheres.

This study shows the strong dependence of the global Nusselt number to the solid volume fraction and Reynolds number. Indeed, increasing the solid volume fraction increases the blocking effect which reduces the thickness of the thermal boundary layer and at the same time, increases the temperature gradient and consequently the Nusselt number. This effect accentuate significantly with the solid volume fraction. Indeed, Nusselt number for a sphere in FCC configuration can be five times higher than for an isolated sphere for high particle concentrations. Moreover, the Nusselt number of a sphere in an FFC arrangement seems to diverge from the Gunn correlation [14] for large solid volume fraction. Therefore, the FCC configuration is not a suitable approximation for a uniform flow past a random assembly of spheres at dense regime.

5. Finite Size Random Arrangement of Spheres

In order to assess the ability of the presented method to deal with more complex particulate flows, the last case studied in this work is a uniform flow past a random assembly of fixed hot spheres of diameter d (illustrated in Figure 12). The computational domain lengths are $L_x=9d$ and $L_y=L_z=5d$. Given a solid volume fraction ${\alpha }_d$, the particles are randomly distributed in the sub-domain $\left[\left(2d,0,0\right);\left(7d,5d,5d\right)\right]$. The assembly of spheres is filled by generating one sphere after an other. The position of the new sphere is randomly generated and a very straightforward overlap test is conducted with the already generated spheres. If no overlap is detected, the new particle is added to the assembly, otherwise, the sphere is rejected and a new sphere is generated in the same manner. The procedure is repeated until the desired solid volume fraction is reached. Furthermore, a cold fluid flows through the channel in the x-direction (streamwise), by imposing at the inlet a temperature $T_{\infty }$ and a velocity ${\displaystyle U_{\infty }\mathbf{x}}$ deduced from Reynolds number ${\displaystyle Re=\frac{U_{\infty }d}{\nu }}$. Periodic boundary conditions are imposed in crosswise directions. The $N_p$ spheres are maintained at a constant temperature $T_s$, the number of spheres is deduced from the solid volume fraction by ${\displaystyle N_p=\frac{6}{\pi }{\left(\frac{L_y}{d}\right)}^3{\alpha }_d}$. As for the two previous cases, the grid resolution and the time step are given by: ${\displaystyle \Delta x=\frac{d}{5\sqrt{Re}}}$, ${\displaystyle \Delta t=0.5min\left(\frac{\Delta x}{U_{\infty }},\frac{\Delta x^2}{k_f/C_v}\right)}$.

A fixed bed Nusselt number is considered here. It is defined as an average of individual particle Nusselt numbers:

$$\langle Nu\rangle =\frac{1}{N_p}\sum _i^{N_p}\frac{{\left(Q_p\right)}_i}{\left(T_s-{\left({\langle T\rangle }_f\right)}_i\right)\pi d^2}$$

(18)

where ${\left(Q_p\right)}_i$ is the $i\mathrm{th}$ sphere heat flux given by ‘Equation (8)’. Estimating the local average temperature of the fluid ${\left({\langle T\rangle }_f\right)}_i$ is not straightforward. As suggested by Deen et al. [12], a good estimate is given by:

$${\left({\langle T\rangle }_f\right)}_i=\frac{\int \int \int (1-C)T{g}_i\left(r\right)}{\int \int \int (1-C)g_i\left(r\right)}$$

(19)

where T and C are the temperature field and color function respectively, and $g_i\left(r\right)$ is the filter that covers the fluid volume in a $4d$ length box containing the $i\mathrm{th}$ sphere, both having the same center. It is given by:

$$g_i\left(r\right)=exp\left(\frac{-2r}{d}\right),0<r<2d$$

(20)

Simulations of uniform flow past a random assembly of hot spheres were performed for low solid volume fractions ${\alpha }_d=0.05,0.1$ at different Reynolds numbers $Re=10,50,100$. Their results are compared to the following correlations:

• Gun [14]:

  $$\begin{array}{cc}\hfill Nu=& \left(7-10\left(1-{\alpha }_d\right)+5{\left(1-{\alpha }_d\right)}^2\right)\left(1+0.7R{e}^{0.2}P{r}^{\frac{1}{3}}\right)\hfill \\ & +\left(1.33-2.4\left(1-{\alpha }_d\right)+1.2{\left(1-{\alpha }_d\right)}^2\right)R{e}^{0.7}P{r}^{\frac{1}{3}}\hfill \end{array}$$

  (21)
• Deen et al. [16]:

  $$\begin{array}{cc}\hfill Nu=& \left(7-10\left(1-{\alpha }_d\right)+5{\left(1-{\alpha }_d\right)}^2\right)\left(1+0.17R{e}^{0.2}P{r}^{\frac{1}{3}}\right)\hfill \\ & +\left(1.33-2.31\left(1-{\alpha }_d\right)+1.16{\left(1-{\alpha }_d\right)}^2\right)R{e}^{0.7}P{r}^{\frac{1}{3}}\hfill \end{array}$$

  (22)
• Sun et al. [17]:

  $$\begin{array}{cc}\hfill Nu=& \frac{\left(-0.46+1.77\left(1-{\alpha }_d\right)+0.69{\left(1-{\alpha }_d\right)}^2\right)}{{\left(1-{\alpha }_d\right)}^3}\hfill \\ & +\left(1.37-2.4\left(1-{\alpha }_d\right)+1.2{\left(1-{\alpha }_d\right)}^2\right)R{e}^{0.7}P{r}^{\frac{1}{3}}\hfill \end{array}$$

  (23)

Figure 13 shows that our results are in the same range of values than the results given in the literature, which is very encouraging as first results of real particle flow arrangements.

6. Conclusions

In the framework of a finite-size particle approach, an accurate procedure for the calculation of heat transfer coefficients for motions around fixed particles has been proposed based on third order Lagrange extrapolation coupled to third order Taylor interpolation methods. This method is an extension of the techniques proposed for hydrodynamic force estimates published in [19]. The procedure allowing the estimate of heat transfer coefficients at a particle surface has been validated with success on fixed isolated spheres, FCC and random arrangements of spheres. Comparisons have been presented to reference Nusselt correlations of the literature.

Future works will be devoted to applying our heat transfer calculation method to arrangements of ellipsoidal particles and extracting correlation laws for these configurations. Moreover, a more extensive study of heat transfer in a fluidized bed is now possible.