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

Mohamed-Amine Chadil, Stéphane Vincent, Jean-Luc Estivalèzes

► To cite this version:

Mohamed-Amine Chadil, Stéphane Vincent, Jean-Luc Estivalèzes. Gas-Solid Heat Transfer Computation from Particle-Resolved Direct Numerical Simulations. Fluids, 2021, 7 (1), pp.15. 10.3390/fluids7010015. [https://arxiv.org/abs/2007.01507](https://arxiv.org/abs/2007.01507)

HAL Id: hal-04390037

https://hal.science/hal-04390037

Submitted on 5 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

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

Mohamed-Amine Chadil <fnref n="1" />,* [ID], Stéphane Vincent <fnref n="1" /> and Jean-Luc Estivalèzes <fnref n="2" />

**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.

**Keywords:** gas-solid flows; Particle-Resolved DNS; fictitious domain; Viscous Penalty Method; heat transfer; assemblies of spheres

<figure><img src="image_1.png" /></figure>

check for
updates

**Citation:** Chadil, M.-A.; Vincent, S.; Estivalezes, J.-L. Gas-Solid Heat Transfer Computation from Particle-Resolved Direct Numerical Simulations. *Fluids* **2022**, *7*, 15. https://doi.org/10.3390/fluids7010015

Academic Editors: Mehrdad Massoudi and Sourabh V. Apte

Received: 10 September 2021

Accepted: 9 December 2021

Published: 30 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_2.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 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–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]

<footnote n="1">¹ MSME, CNRS UMR 8208, University Gustave Eiffel, University Paris Est Creteil, F-77474 Marne-la-Vallée, France; stephane.vincent@univ-eiffel.fr</footnote>

<footnote n="2">² Institut de Mécanique des Fluides de Toulouse, IMFT, Université de Toulouse, CNRS, F-31400 Toulouse, France; Jean-Luc.Estivalezes@onera.fr</footnote>

<footnote>* Correspondence: amine.chadil@cnrs.fr</footnote>
Article

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

Mohamed-Amine Chadil <fnref n="1" />,* [ID], Stéphane Vincent <fnref n="1" /> and Jean-Luc Estivalèzes <fnref n="2" />

**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.

**Keywords:** gas-solid flows; Particle-Resolved DNS; fictitious domain; Viscous Penalty Method; heat transfer; assemblies of spheres

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Chadil, M.-A.; Vincent, S.; Estivalezes, J.-L. Gas-Solid Heat Transfer Computation from Particle-Resolved Direct Numerical Simulations. *Fluids* **2022**, *7*, 15. https://doi.org/10.3390/fluids7010015

Academic Editors: Mehrdad Massoudi and Sourabh V. Apte

Received: 10 September 2021

Accepted: 9 December 2021

Published: 30 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 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–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]

<footnote n="1">¹ MSME, CNRS UMR 8208, University Gustave Eiffel, University Paris Est Creteil, F-77474 Marne-la-Vallée, France; stephane.vincent@univ-eiffel.fr</footnote>

<footnote n="2">² Institut de Mécanique des Fluides de Toulouse, IMFT, Université de Toulouse, CNRS, F-31400 Toulouse, France; Jean-Luc.Estivalezes@onera.fr</footnote>

<footnote>* Correspondence: amine.chadil@cnrs.fr</footnote>
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:

$$ \nabla \cdot \mathbf{u} = 0 \qquad (1) $$

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) = -\nabla p + \rho \mathbf{g} + \nabla \cdot \left[ \mu (\nabla \mathbf{u} + \nabla^t \mathbf{u}) \right] + \mathbf{F}_m \qquad (2) $$

$$ \rho C_v \left( \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T \right) = \nabla \cdot \left( k_f \nabla T \right) \qquad (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 sense, specific arithmetic $(\rho = C\rho_s + (1-C)\rho_f, C_v = C C_{v_s} + (1-C) C_{v_f})$, harmonic $(\mu = \frac{\mu_s\mu_f}{C\mu_f + (1-C)\mu_s}, k_f = \frac{k_{fs}k_{ff}}{Ck_{ff} + (1-C)k_{fs}})$ 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 $\frac{dx_i}{dt} = \mathbf{V}_i$, where $\mathbf{V}_i$ is the velocity of the $i$th particle deduced from the Eulerian velocities inside the particle ($x_i$ being its position). This procedure is called VOF-Lag approach in contrast with the classical VOF method [25] where the equation $\frac{\partial C}{\partial t} + \mathbf{u} \cdot \nabla C = 0$ is solved to update the color function C.
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:

$$ \nabla \cdot \mathbf{u} = 0 \qquad (1) $$

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) = -\nabla p + \rho \mathbf{g} + \nabla \cdot \left[ \mu (\nabla \mathbf{u} + \nabla^t \mathbf{u}) \right] + \mathbf{F}_m \qquad (2) $$

$$ \rho C_v \left( \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T \right) = \nabla \cdot \left( k_f \nabla T \right) \qquad (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 sense, specific arithmetic $(\rho = C\rho_s + (1-C)\rho_f$, $C_v = C C_{v_s} + (1-C) C_{v_f}$), harmonic $(\mu = \frac{\mu_s \mu_f}{C\mu_f + (1-C)\mu_s}$, $k_f = \frac{k_{fs} k_{ff}}{Ck_{ff} + (1-C)k_{fs}})$ 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 $\frac{dx_i}{dt} = V_i$, where $V_i$ is the velocity of the $i$th particle deduced from the Eulerian velocities inside the particle ($x_i$ being its position). This procedure is called VOF-Lag approach in contrast with the classical VOF method [25] where the equation $\frac{\partial C}{\partial t} + \mathbf{u} \cdot \nabla C = 0$ is solved to update the color function C.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Projected phase function C on staggered grids. Pressure in black circles, the velocities as arrows and viscosity in white circles [9].</figcaption></figure>

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 $\mu \to +\infty$ and based on the equivalence [9] $(\forall \mathbf{P} \in \Omega_i, \nabla \mathbf{u} + \nabla^T \mathbf{u} = 0 \Leftrightarrow \mathbf{V}(\mathbf{P}) = \mathbf{V}_i + \mathbf{O}_i \mathbf{P} \wedge \omega_i)$ 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]:

<caption>**Algorithm 1 Augmented Lagrangian (AL) method [27]**</caption>

for each iteration n do
    $\mathbf{u}^{*,0} = \mathbf{u}^n$ and $p^{*,0} = p^n$
    $m = 0$
    **while** $||\nabla \cdot \mathbf{u}^{*,m}|| > \epsilon_{AL}$ **do**
        $m = m + 1$
        solve: $\rho \left( \frac{\mathbf{u}^{*,m} - \mathbf{u}^{*,0}}{\Delta t} + \mathbf{u}^{*,m-1} \cdot \nabla \mathbf{u}^{*,m} \right) - \nabla(r\nabla \cdot \mathbf{u}^{*,m}) = -\nabla p^{*,m-1} + \rho \mathbf{g} + \nabla \cdot [\mu(\nabla \mathbf{u}^{*,m} + \nabla^T \mathbf{u}^{*,m})]$
        $p^{*,m} = p^{*,m-1} - r\nabla \cdot \mathbf{u}^{*,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} \cdot \nabla \mathbf{u} \approx (2\mathbf{u}^n - \mathbf{u}^{n-1}) \cdot \nabla \mathbf{u}^{n+1} \quad (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].
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Projected phase function C on staggered grids. Pressure in black circles, the velocities as arrows and viscosity in white circles [9].</figcaption></figure>

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 $\mu \to +\infty$ and based on the equivalence [9] $(\forall \mathbf{P} \in \Omega_i, \nabla \mathbf{u} + \nabla^T \mathbf{u} = 0 \Leftrightarrow \mathbf{V}(\mathbf{P}) = \mathbf{V}_i + \mathbf{O}_i \mathbf{P} \wedge \omega_i)$ 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]:

<caption>**Algorithm 1 Augmented Lagrangian (AL) method [27]**</caption>

for each iteration n do
    $\mathbf{u}^{*,0} = \mathbf{u}^n$ and $p^{*,0} = p^n$
    $m = 0$
    **while** $||\nabla \cdot \mathbf{u}^{*,m}|| > \epsilon_{AL}$ **do**
        $m = m + 1$
        solve: $\rho \left( \frac{\mathbf{u}^{*,m} - \mathbf{u}^{*,0}}{\Delta t} + \mathbf{u}^{*,m-1} \cdot \nabla \mathbf{u}^{*,m} \right) - \nabla(r\nabla \cdot \mathbf{u}^{*,m}) = -\nabla p^{*,m-1} + \rho \mathbf{g} + \nabla \cdot [\mu(\nabla \mathbf{u}^{*,m} + \nabla^T \mathbf{u}^{*,m})]$
        $p^{*,m} = p^{*,m-1} - r\nabla \cdot \mathbf{u}^{*,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} \cdot \nabla \mathbf{u} \approx (2\mathbf{u}^n - \mathbf{u}^{n-1}) \cdot \nabla \mathbf{u}^{n+1} \quad (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, $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} \qquad (5)$$

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

$$H_p = \frac{Q_p}{(T_s - T_f) \pi d^2} \qquad (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 \cdot \mathbf{n} \, dS \qquad (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 (\nabla T)_l \cdot \mathbf{n}_l \Delta S_l \qquad (8)$$

where the ongoing normal to the $l$th 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 $(\nabla T)_l$ is not known at the $l$th 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 $\bar{\bar{\sigma}}_l$ components on the Lagrangian mesh in the hydrodynamic force computation $\mathbf{F} \approx \sum_{l}^{N} \bar{\bar{\sigma}}_l \cdot \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:

$$(\nabla T)(\mathbf{C}_l) = \sum_{i=1}^{3} (\nabla T)(\mathbf{P}_i) L_i(\mathbf{C}_l) \quad , \text{where} \quad L_i(\mathbf{C}_l) = \prod_{j \neq i}^{3} \frac{|\mathbf{C}_l - \mathbf{P}_j|}{|\mathbf{P}_i - \mathbf{P}_j|} \qquad (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, $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} \qquad (5)$$

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

$$H_p = \frac{Q_p}{(T_s - T_f) \pi d^2} \qquad (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 \cdot \mathbf{n} \, dS \qquad (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 (\nabla T)_l \cdot \mathbf{n}_l \Delta S_l \qquad (8)$$

where the ongoing normal to the $l$th 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 $(\nabla T)_l$ is not known at the $l$th 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 $\mathbf{F} \approx \sum_{l}^{N} \overline{\overline{\sigma}}_l \cdot \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:

$$(\nabla T)(\mathbf{C}_l) = \sum_{i=1}^{3} (\nabla T)(\mathbf{P}_i) L_i(\mathbf{C}_l) \quad , \text{where} \quad L_i(\mathbf{C}_l) = \prod_{j \neq i}^{3} \frac{|\mathbf{C}_l - \mathbf{P}_j|}{|\mathbf{P}_i - \mathbf{P}_j|} \qquad (9)$$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Notations for the third order Lagrange extrapolation (**left**) and the third order Taylor interpolation (**right**), used in heat flux computations.</figcaption></figure>

Given that the temperature gradient is computed on the Eulerian mesh and that the extrapolation points $P_1, P_2$ and $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 $P_k$ ($k = 1, 2, 3$) using the third-order Taylor interpolation:

$$ 
\begin{aligned}
f(\mathbf{P}_k) = f(\mathbf{E}) &+ \frac{\partial f}{\partial x_i}(\mathbf{E})((\mathbf{C}_l)_i - \mathbf{E}_i) + \frac{1}{2} \frac{\partial^2 f}{\partial x_i^2}(\mathbf{E})((\mathbf{C}_l)_i - \mathbf{E}_i)^2 \\
&+ \frac{\partial^2 f}{\partial x_i \partial x_j}(\mathbf{E})((\mathbf{C}_l)_i - \mathbf{E}_i)((\mathbf{C}_l)_j - \mathbf{E}_j) + O(||\mathbf{C}_l - \mathbf{E}||^3)
\end{aligned} 
\quad (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:

$$ Nu_{loc}(\theta) = \frac{-\nabla T \cdot \mathbf{n} \, d}{(T_s - T_f)} \quad (11) $$

where for each Lagrangian mesh element center ($\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}, \quad \beta = \arctan2(z, y) \frac{180}{\pi} \quad (12) $$

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Spherical coordinate system around a particle. The flow direction is represented by the undisturbed velocity $\mathbf{U}_\infty$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Notations for the third order Lagrange extrapolation (**left**) and the third order Taylor interpolation (**right**), used in heat flux computations.</figcaption></figure>

Given that the temperature gradient is computed on the Eulerian mesh and that the extrapolation points $P_1, P_2$ and $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 $P_k$ ($k = 1, 2, 3$) using the third-order Taylor interpolation:

$$f(\mathbf{P}_k) = f(\mathbf{E}) + \frac{\partial f}{\partial x_i}(\mathbf{E})((\mathbf{C}_l)_i - \mathbf{E}_i) + \frac{1}{2} \frac{\partial^2 f}{\partial x_i^2}(\mathbf{E})((\mathbf{C}_l)_i - \mathbf{E}_i)^2 \\ + \frac{\partial^2 f}{\partial x_i \partial x_j}(\mathbf{E})((\mathbf{C}_l)_i - \mathbf{E}_i)((\mathbf{C}_l)_j - \mathbf{E}_j) + O(||\mathbf{C}_l - \mathbf{E}||^3) \quad (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:

$$Nu_{loc}(\theta) = \frac{-\nabla T \cdot \mathbf{n} \, d}{(T_s - T_f)} \quad (11)$$

where for each Lagrangian mesh element center ($\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}, \quad \beta = \arctan2(z, y) \frac{180}{\pi} \quad (12)$$

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Spherical coordinate system around a particle. The flow direction is represented by the undisturbed velocity $U_\infty$.</figcaption></figure>
### 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 [(2d, 3d, 3d); (6d, 5d, 5d)] 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 $\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 $\Delta t = 0.5 \min\left(\frac{\Delta x}{U_\infty}, \frac{\Delta x^2}{k_f/C_v}\right)$.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Streamlines and temperature field in a slice for a uniform flow past a sphere at $Re = 100$. The bottom figure is a zoom on the particle.</figcaption></figure>

#### 3.1. *Effect of the Extrapolation Distance* ($\delta$)

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 $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 $\delta = \Delta x$.
### 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 [$(2d, 3d, 3d)$; $(6d, 5d, 5d)$] 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 $\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 $\Delta t = 0.5 \min\left(\frac{\Delta x}{U_\infty}, \frac{\Delta x^2}{k_f/C_v}\right)$.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Streamlines and temperature field in a slice for a uniform flow past a sphere at $Re = 100$. The bottom figure is a zoom on the particle.</figcaption></figure>

#### 3.1. *Effect of the Extrapolation Distance* ($\delta$)

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 $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 $\delta = \Delta x$.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Drag force relative error to Schiller and Naumann [30] correlation (**left**) and Nusselt number relative error to Ranz&Marshall [31] correlation (**right**), as a function of the distance $δ$ at $Re = 100$.</figcaption></figure>

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 $δ$ (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 $δ = 0$ does not reflect the flow symmetry as illustrated in (Figure 6a), and the error between our result ($Nu_{loc}$) and Massol's ($Nu_{locMassol}$) given by $\frac{\sum_{i=1}^{N} |(Nu_{loc})_i - (Nu_{locMassol})_i|}{\sum_{i=1}^{N} |(Nu_{locMassol})_i|}$ is about 25%. This error decreases as the distance $δ$ increases until being lower than 10% for $δ \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 ($β$) plans (Figure 6b). So, as for the hydrodynamic force computation [19] and in the rest of this work, the extrapolation distance will always be $δ = 1$.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Local Nusselt coefficient for a uniform flow past an isolated sphere at $Re = 100$, as function of the azimuthal angle $θ$, for all points on the particle surface (for all $β$ plans) at (**a**) $δ = 0$ and (**b**) $δ = 1$: (.) present work, (□) Massol [10]. (**c**) is local Nusselt number relative error to Massol's result [10] as function of the extrapolation distance $δ$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Drag force relative error to Schiller and Naumann [30] correlation (**left**) and Nusselt number relative error to Ranz&Marshall [31] correlation (**right**), as a function of the distance $δ$ at $Re = 100$.</figcaption></figure>

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 $δ$ (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 $δ = 0$ does not reflect the flow symmetry as illustrated in (Figure 6a), and the error between our result ($Nu_{loc}$) and Massol's ($Nu_{locMassol}$) given by $\frac{\sum_{i=1}^{N} |(Nu_{loc})_i - (Nu_{locMassol})_i|}{\sum_{i=1}^{N} |(Nu_{locMassol})_i|}$ is about 25%. This error decreases as the distance $δ$ increases until being lower than 10% for $δ \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 ($β$) plans (Figure 6b). So, as for the hydrodynamic force computation [19] and in the rest of this work, the extrapolation distance will always be $δ = 1$.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Local Nusselt coefficient for a uniform flow past an isolated sphere at $Re = 100$, as function of the azimuthal angle $θ$, for all points on the particle surface (for all $β$ plans) at (**a**) $δ = 0$ and (**b**) $δ = 1$: (.) present work, (□) Massol [10]. (**c**) is local Nusselt number relative error to Massol's result [10] as function of the extrapolation distance $δ$.</figcaption></figure>
## 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 ($Re = 10, 20, 40, 60, 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.6 Re^{0.5} Pr^{\frac{1}{3}} \qquad (13)$$

* Feng and Michaelides [32]:

$$Nu = 0.992 + (RePr)^{\frac{1}{3}} + 0.1 Re^{\frac{2}{3}} Pr^{\frac{1}{3}} \qquad (14)$$

* Whitaker [33]:

$$Nu = 2. + (0.4 Re^{0.5} + 0.06 Re^{\frac{2}{3}}) Pr^{0.4} \qquad (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.

<figure><img src="image_9.png" /><figcaption>Figure 7. Nusselt coefficient for the uniform flow past a hot sphere at different Reynolds numbers: (—) Ranz and Marshall [31], (- - -) Feng and Michaelides [32], (· · ·) Whitaker [33], (●) present work.</figcaption></figure>

## 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

$α_d$: $L = (\frac{2}{3} \frac{\pi}{\alpha_d})^{1/3}$ d. Given a Reynolds number $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

$$F_m = -\rho \frac{\langle \mathbf{u}_f \rangle - \mathbf{u}_d}{\Delta t} + \frac{\sum_{i}^{N_p} F_i^D}{V_f} \qquad (16)$$

inserted in the momentum conservation Equation (2). $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. $F_m$ is thus adjusted until the stationary state is reached when $\langle \mathbf{u}_f \rangle = \mathbf{u}_d$. As for the uniform flow past an isolated sphere, the grid resolution is fixed by imposing
## 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 ($Re = 10, 20, 40, 60, 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.6 Re^{0.5} Pr^{\frac{1}{3}} \qquad (13)$$

* Feng and Michaelides [32]:

$$Nu = 0.992 + (RePr)^{\frac{1}{3}} + 0.1 Re^{\frac{2}{3}} Pr^{\frac{1}{3}} \qquad (14)$$

* Whitaker [33]:

$$Nu = 2. + (0.4 Re^{0.5} + 0.06 Re^{\frac{2}{3}}) Pr^{0.4} \qquad (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.

<figure><img src="image_10.png" /><figcaption>Figure 7. Nusselt coefficient for the uniform flow past a hot sphere at different Reynolds numbers: (—) Ranz and Marshall [31], (- - -) Feng and Michaelides [32], (· · ·) Whitaker [33], (●) present work.</figcaption></figure>

## 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

$α_d$: $L = (\frac{2}{3} \frac{\pi}{\alpha_d})^{1/3}$ d. Given a Reynolds number $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

$$F_m = -\rho \frac{\langle \mathbf{u}_f \rangle - \mathbf{u}_d}{\Delta t} + \frac{\sum_{i}^{N_p} F_i^D}{V_f} \qquad (16)$$

inserted in the momentum conservation Equation (2). $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. $F_m$ is thus adjusted until the stationary state is reached when $\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: $\Delta x = \frac{d}{5\sqrt{Re}}$, and the time step is given by

$$ \Delta t = 0.5 \min \left( \frac{\Delta x}{U_{\infty}}, \frac{\Delta x^2}{k_f/C_v} \right). $$

<figure><img src="image_10.png" /><figcaption>Figure 8. Streamlines and temperature field for a steady flow along the x-axis of a Face-Centered Cubic array of spheres at $Re = 50$, $\alpha_d = 0.15$ (left), and at $Re = 300$, $\alpha_d = 0.15$ (right).</figcaption></figure>

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 $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} Nu_i \qquad \langle Nu_{loc} \rangle = \frac{1}{N_p} \sum_{i}^{N_p} (Nu_{loc})_i \qquad (17) $$

where $N_p = 4$ for a FCC array.
five Cartesian cells in the boundary layer: $\Delta x = \frac{d}{5\sqrt{Re}}$, and the time step is given by

$$ \Delta t = 0.5 \min \left( \frac{\Delta x}{U_{\infty}}, \frac{\Delta x^2}{k_f / C_v} \right). $$

<figure><img src="image_11.png" /><figcaption>Figure 8. Streamlines and temperature field for a steady flow along the x-axis of a Face-Centered Cubic array of spheres at $Re = 50$, $\alpha_d = 0.15$ (left), and at $Re = 300$, $\alpha_d = 0.15$ (right).</figcaption></figure>

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 $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} Nu_i \qquad \langle Nu_{loc} \rangle = \frac{1}{N_p} \sum_{i}^{N_p} (Nu_{loc})_i \qquad (17) $$

where $N_p = 4$ for a FCC array.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** **Left:** Time evolution of fluid temperature for a uniform flow past a Face-Centered Cubic array of spheres at $Re = 50$ and (◊) $\alpha_d = 0.05$, (△) $\alpha_d = 0.1$, (□) $\alpha_d = 0.15$, (◆) $\alpha_d = 0.3$, (●) $\alpha_d = 0.4$, (▲) $\alpha_d = 0.5$, (○) $\alpha_d = 0.6$. **Right:** Time evolution of the averaged Nusselt coefficient for a uniform flow past a Face-Centered Cubic array of spheres at $Re = 50$ and $\alpha_d = 0.4$.</figcaption></figure>

## 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 ($\beta = 0$) and ($\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 $\theta = 0$, due to the fountain effect that induces a low heat flux. With increasing $\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.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** **Left:** Local Nusselt number for a uniform flow past an isolated sphere at $Re = 50$ (▲) and for a uniform flow past a FCC array of spheres at $Re = 50$, $\alpha_d = 0.15$: $\beta = 0$ ((△) Massol [10], (◆) present work) and $\beta = 45$ ((□) Massol [10], (●) present work). **Right:** Local Nusselt number for a uniform flow past a FCC array of spheres at $Re = 300$, $\alpha_d = 0.15$: $\beta = 0$ (○), $\beta = 22.5$ (△), $\beta = 45$ (□).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** **Left:** Time evolution of fluid temperature for a uniform flow past a Face-Centered Cubic array of spheres at $Re = 50$ and (◊) $\alpha_d = 0.05$, (△) $\alpha_d = 0.1$, (□) $\alpha_d = 0.15$, (◆) $\alpha_d = 0.3$, (●) $\alpha_d = 0.4$, (▲) $\alpha_d = 0.5$, (○) $\alpha_d = 0.6$. **Right:** Time evolution of the averaged Nusselt coefficient for a uniform flow past a Face-Centered Cubic array of spheres at $Re = 50$ and $\alpha_d = 0.4$.</figcaption></figure>

## 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 ($\beta = 0$) and ($\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 $\theta = 0$, due to the fountain effect that induces a low heat flux. With increasing $\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.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** **Left:** Local Nusselt number for a uniform flow past an isolated sphere at $Re = 50$ (▲) and for a uniform flow past a FCC array of spheres at $Re = 50$, $\alpha_d = 0.15$: $\beta = 0$ ((△) Massol [10], (◆) present work) and $\beta = 45$ ((□) Massol [10], (●) present work). **Right:** Local Nusselt number for a uniform flow past a FCC array of spheres at $Re = 300$, $\alpha_d = 0.15$: $\beta = 0$ (○), $\beta = 22.5$ (△), $\beta = 45$ (□).</figcaption></figure>
## 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]) $Nu_s = 2. + 0.6 Re^{0.5} Pr^{1/3}$, noted $\frac{\langle Nu \rangle}{Nu_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 ($0.05 \le \alpha_d \le 0.6$). A good agreement is also observed with Gun [14] correlation for low solid volume fractions ($\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.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Global Nusselt coefficient for a uniform flow past a FCC, normalized by Ranz and Marshall [31] Nusselt coefficient for a uniform flow past an isolated sphere. Evolution as a function of the solid volume fraction $\alpha_d$ at Reynolds numbers (a) $Re = 10$, (b) $Re = 50$, (c) $Re = 100$: (—) Gunn [14], (□) Massol [10], (●) present work.</figcaption></figure>

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 [(2d, 0, 0); (7d, 5d, 5d)]. 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,
## 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]) $Nu_s = 2. + 0.6 Re^{0.5} Pr^{1/3}$, noted $\frac{\langle Nu \rangle}{Nu_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 ($0.05 \le \alpha_d \le 0.6$). A good agreement is also observed with Gun [14] correlation for low solid volume fractions ($\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.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Global Nusselt coefficient for a uniform flow past a FCC, normalized by Ranz and Marshall [31] Nusselt coefficient for a uniform flow past an isolated sphere. Evolution as a function of the solid volume fraction $\alpha_d$ at Reynolds numbers (a) $Re = 10$, (b) $Re = 50$, (c) $Re = 100$: (—) Gunn [14], (□) Massol [10], (●) present work.</figcaption></figure>

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 $[(2d, 0, 0); (7d, 5d, 5d)]$. 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_∞$ and a velocity $U_∞x$ deduced from Reynolds number $Re = \frac{U_∞d}{v}$. 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 $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: $\Delta x = \frac{d}{5\sqrt{Re}}$, $\Delta t = 0.5 \min \left(\frac{\Delta x}{U_∞}, \frac{\Delta x^2}{k_f/C_v}\right)$.

<figure><img src="image_14.png" /><figcaption>Figure 12. Streamlines and temperature field of a steady flow along the x-axis of a random array of spheres at $Re = 100$ and $\alpha_d = 0.1$.</figcaption></figure>

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{(Q_p)_i}{(T_s - (\langle T \rangle_f)_i) \pi d^2} \quad (18) $$

where $(Q_p)_i$ is the ith sphere heat flux given by 'Equation (8)'. Estimating the local average temperature of the fluid $(\langle T \rangle_f)_i$ is not straightforward. As suggested by Deen et al. [12], a good estimate is given by:

$$ (\langle T \rangle_f)_i = \frac{\iiint (1-C) T g_i(r)}{\iiint (1-C) g_i(r)} \quad (19) $$

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

$$ g_i(r) = \exp\left(\frac{-2r}{d}\right), \quad 0 < r < 2d \quad (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]:

$$ Nu = (7 - 10(1 - \alpha_d) + 5(1 - \alpha_d)^2) (1 + 0.7Re^{0.2}Pr^{\frac{1}{3}}) + (1.33 - 2.4(1 - \alpha_d) + 1.2(1 - \alpha_d)^2) Re^{0.7}Pr^{\frac{1}{3}} $$
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_∞$ and a velocity $U_∞x$ deduced from Reynolds number $Re = \frac{U_∞d}{v}$. 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 $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: $\Delta x = \frac{d}{5\sqrt{Re}}$, $\Delta t = 0.5 \min \left(\frac{\Delta x}{U_∞}, \frac{\Delta x^2}{k_f/C_v}\right)$.

<figure><img src="image_15.png" /><figcaption>Figure 12. Streamlines and temperature field of a steady flow along the x-axis of a random array of spheres at $Re = 100$ and $\alpha_d = 0.1$.</figcaption></figure>

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{(Q_p)_i}{(T_s - (\langle T \rangle_f)_i) \pi d^2} \quad (18) $$

where $(Q_p)_i$ is the ith sphere heat flux given by 'Equation (8)'. Estimating the local average temperature of the fluid $(\langle T \rangle_f)_i$ is not straightforward. As suggested by Deen et al. [12], a good estimate is given by:

$$ (\langle T \rangle_f)_i = \frac{\iiint (1-C) T g_i(r)}{\iiint (1-C) g_i(r)} \quad (19) $$

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

$$ g_i(r) = \exp\left(-\frac{2r}{d}\right), \quad 0 < r < 2d \quad (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]:

$$ Nu = (7 - 10(1 - \alpha_d) + 5(1 - \alpha_d)^2)(1 + 0.7Re^{0.2}Pr^{\frac{1}{3}}) + (1.33 - 2.4(1 - \alpha_d) + 1.2(1 - \alpha_d)^2)Re^{0.7}Pr^{\frac{1}{3}} $$
* Deen et al. [16]:

$$Nu = (7 - 10(1 - \alpha_d) + 5(1 - \alpha_d)^2)(1 + 0.17Re^{0.2}Pr^{\frac{1}{3}}) + (1.33 - 2.31(1 - \alpha_d) + 1.16(1 - \alpha_d)^2)Re^{0.7}Pr^{\frac{1}{3}}$$

* Sun et al. [17]:

$$Nu = \frac{(-0.46 + 1.77(1 - \alpha_d) + 0.69(1 - \alpha_d)^2)}{(1 - \alpha_d)^3} + (1.37 - 2.4(1 - \alpha_d) + 1.2(1 - \alpha_d)^2)Re^{0.7}Pr^{\frac{1}{3}}$$

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.

<figure><img src="image_15.png" /><figcaption>Figure 13. Nusselt coefficient for a uniform flow past a random arrangement of hot spheres as a function of the solid volume fraction $\alpha_d$ at various Reynolds numbers (a) $Re = 10$, (b) $Re = 50$, (c) $Re = 100$: (---) Gunn [14], (—) Deen et al. [16], (---) Sun et al. [17] and (●) present work.</figcaption></figure>

## 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.

**Author Contributions:** Conceptualization, M.-A.C., S.V. and J.-L.E.; methodology, M.-A.C. and S.V.; software, M.-A.C. and S.V.; validation, M.-A.C.; formal analysis, M.-A.C.; investigation, M.-A.C.; resources, M.-A.C.; data curation, M.-A.C.; writing—original draft preparation, M.-A.C. and S.V.;
* Deen et al. [16]:

$$Nu = (7 - 10(1 - \alpha_d) + 5(1 - \alpha_d)^2)(1 + 0.17Re^{0.2}Pr^{\frac{1}{3}}) + (1.33 - 2.31(1 - \alpha_d) + 1.16(1 - \alpha_d)^2)Re^{0.7}Pr^{\frac{1}{3}}$$

* Sun et al. [17]:

$$Nu = \frac{(-0.46 + 1.77(1 - \alpha_d) + 0.69(1 - \alpha_d)^2)}{(1 - \alpha_d)^3} + (1.37 - 2.4(1 - \alpha_d) + 1.2(1 - \alpha_d)^2)Re^{0.7}Pr^{\frac{1}{3}}$$

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.

<figure><img src="image_16.png" /><figcaption>Figure 13. Nusselt coefficient for a uniform flow past a random arrangement of hot spheres as a function of the solid volume fraction $\alpha_d$ at various Reynolds numbers (a) $Re = 10$, (b) $Re = 50$, (c) $Re = 100$: (---) Gunn [14], (—) Deen et al. [16], (---) Sun et al. [17] and (●) present work.</figcaption></figure>

## 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.

**Author Contributions:** Conceptualization, M.-A.C., S.V. and J.-L.E.; methodology, M.-A.C. and S.V.; software, M.-A.C. and S.V.; validation, M.-A.C.; formal analysis, M.-A.C.; investigation, M.-A.C.; resources, M.-A.C.; data curation, M.-A.C.; writing—original draft preparation, M.-A.C. and S.V.;
writing—review and editing, M.-A.C. and S.V.; visualization, M.-A.C.; supervision, S.V. and J.-L.E.; project administration, S.V.; funding acquisition, S.V. and J.-L.E. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by ANR grant number ANR-11-BS09-0004.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** The authors gratefully acknowledge support from CINES (National Computing Center for Higher Education), IDRIS (Institut du développement et des ressources en informatique scientifique) and CCRT (National Computing Center of CEA) under project number A0092B06115.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Bagheri, G.; Bonadonna, C.; Manzella, I.; Vonlanthen, P. On the characterization of size and shape of irregular particles. *Powder Technol.* **2015**, *270*, 141–153. [CrossRef]

2. Dartevelle, S. From Model Conception to Verification and Validation, a Global Approach to Multiphase Navier-Stoke Models with an Emphasis on Volcanic Explosive Phenomenology. Technical Report LA-14346, 948564. 2007. Available online: https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-14346 (accessed on 9 September 2021).

3. Adanez, J.; Abad, A.; Garcia-Labiano, F.; Gayan, P.; de Diego, L.F. Progress in Chemical-Looping Combustion and Reforming technologies. *Prog. Energy Combust. Sci.* **2012**, *38*, 215–282. [CrossRef]

4. Lyngfelt, A.; Leckner, B.; Mattisson, T. A fluidized-bed combustion process with inherent CO₂ separation; application of chemical-looping combustion. *Chem. Eng. Sci.* **2001**, *56*, 3101–3113. [CrossRef]

5. Mattisson, T.; Keller, M.; Linderholm, C.; Moldenhauer, P.; Rydan, M.; Leion, H.; Lyngfelt, A. Chemical-looping technologies using circulating fluidized bed systems: Status of development. *Fuel Process. Technol.* **2018**, *172*, 1–12. [CrossRef]

6. Acrivos, A.; Hinch, E.J.; Jeffrey, D.J. Heat transfer to a slowly moving fluid from a dilute fixed bed of heated spheres. *J. Fluid Mech.* **1980**, *101*, 403–421. [CrossRef]

7. Pfeffer, R.; Happel, J. An analytical study of heat and mass transfer in multiparticle systems at low Reynolds numbers. *AIChE J.* **1964**, *10*, 605–611. [CrossRef]

8. Wakao, N.; Kaguei, S. *Heat and Mass Transfer in Packed Beds Volume 1*; Gordon and Breach Science Publishers Inc.: New York, NY, USA, 1983.

9. Vincent, S.; de Motta, J.C.B.; Sarthou, A.; Estivalezes, J.L.; Simonin, O.; Climent, E. A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows. *J. Compuat. Phys.* **2014**, *256*, 582–614. [CrossRef]

10. Massol, A. Simulations Numériques d'écoulements à Travers des réseaux Fixes de Spheres Monodisperses et Bidisperses, Pour des Nombres de Reynolds Modérés. Ph.D. Thesis, Institut National Polytechnique de Toulouse, Toulouse, France, 2004. Available online: https://www.researchgate.net/publication/242075914_Adaptive_Augmented_Lagrangian_techniques_for_simulating_unsteadyFlows (accessed on 9 September 2021).

11. Schönfeld, T.; Rudyard, M. Steady and Unsteady Flow Simulations Using the Hybrid Flow Solver AVBP. *AIAA J.* **1999**, *37*, 1378–1385. [CrossRef]

12. Deen, N.G.; Kriebitzsch, S.H.; van der Hoef, M.A.; Kuipers, J. Direct numerical simulation of flow and heat transfer in dense fluid-particle systems. *Chem. Eng. Sci.* **2012**, *81*, 329–344. [CrossRef]

13. Tavassoli, H.; Kriebitzsch, S.; van der Hoef, M.; Peters, E.; Kuipers, J. Direct numerical simulation of particulate flow with heat transfer. *Int. J. Multiph. Flow* **2013**, *57*, 29–37. [CrossRef]

14. Gunn, D. Transfer of heat or mass to particles in fixed and fluidized beds. *Int. J. Heat Mass Transf.* **1978**, *21*, 467–476. [CrossRef]

15. Tenneti, S.; Sun, B.; Garg, R.; Subramaniam, S. Role of fluid heating in dense gas-solid flow as revealed by particle-resolved direct numerical simulation. *Int. J. Heat Mass Transf.* **2013**, *58*, 471–479. [CrossRef]

16. Deen, N.G.; Peters, E.A.J.F.; Padding, J.T.; Kuipers, J.A.M. Review of direct numerical simulation of fluid-particle mass, momentum and heat transfer in dense gas-solid flows. *Chem. Eng. Sci.* **2014**, *116*, 710–724. [CrossRef]

17. Sun, B.; Tenneti, S.; Subramaniam, S. Modeling average gas-solid heat transfer using particle-resolved direct numerical simulation. *Int. J. Heat Mass Transf.* **2015**, *86*, 898–913. [CrossRef]

18. Thiam, E.I.; Masi, E.; Climent, E.; Simonin, O.; Vincent, S. Particle-resolved numerical simulations of the gas-solid heat transfer in arrays of random motionless particles. *Acta Mech.* **2019**, *230*, 541–567. [CrossRef]

19. Chadil, M.A.; Vincent, S.; Estivalèzes, J.L. Accurate estimate of drag forces using particle-resolved direct numerical simulations. *Acta Mech.* **2018**, *230*, 569–595. [CrossRef]

20. Kataoka, I. Local instant formulation of two-phase flow. *Int. J. Multiph. Flow* **1986**, *12*, 745–758. [CrossRef]

21. Kataoka, I.; Ishii, M.; Serizawa, A. Local formulation and measurements of interfacial area concentration in two-phase flow. *Int. J. Multiph. Flow* **1986**, *12*, 505–529. [CrossRef]
writing—review and editing, M.-A.C. and S.V.; visualization, M.-A.C.; supervision, S.V. and J.-L.E.; project administration, S.V.; funding acquisition, S.V. and J.-L.E. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by ANR grant number ANR-11-BS09-0004.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** The authors gratefully acknowledge support from CINES (National Computing Center for Higher Education), IDRIS (Institut du développement et des ressources en informatique scientifique) and CCRT (National Computing Center of CEA) under project number A0092B06115.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Bagheri, G.; Bonadonna, C.; Manzella, I.; Vonlanthen, P. On the characterization of size and shape of irregular particles. *Powder Technol.* **2015**, *270*, 141–153. [CrossRef]

2. Dartevelle, S. From Model Conception to Verification and Validation, a Global Approach to Multiphase Navier-Stoke Models with an Emphasis on Volcanic Explosive Phenomenology. Technical Report LA-14346, 948564. 2007. Available online: https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-14346 (accessed on 9 September 2021).

3. Adanez, J.; Abad, A.; Garcia-Labiano, F.; Gayan, P.; de Diego, L.F. Progress in Chemical-Looping Combustion and Reforming technologies. *Prog. Energy Combust. Sci.* **2012**, *38*, 215–282. [CrossRef]

4. Lyngfelt, A.; Leckner, B.; Mattisson, T. A fluidized-bed combustion process with inherent CO₂ separation; application of chemical-looping combustion. *Chem. Eng. Sci.* **2001**, *56*, 3101–3113. [CrossRef]

5. Mattisson, T.; Keller, M.; Linderholm, C.; Moldenhauer, P.; Rydan, M.; Leion, H.; Lyngfelt, A. Chemical-looping technologies using circulating fluidized bed systems: Status of development. *Fuel Process. Technol.* **2018**, *172*, 1–12. [CrossRef]

6. Acrivos, A.; Hinch, E.J.; Jeffrey, D.J. Heat transfer to a slowly moving fluid from a dilute fixed bed of heated spheres. *J. Fluid Mech.* **1980**, *101*, 403–421. [CrossRef]

7. Pfeffer, R.; Happel, J. An analytical study of heat and mass transfer in multiparticle systems at low Reynolds numbers. *AIChE J.* **1964**, *10*, 605–611. [CrossRef]

8. Wakao, N.; Kaguei, S. *Heat and Mass Transfer in Packed Beds Volume 1*; Gordon and Breach Science Publishers Inc.: New York, NY, USA, 1983.

9. Vincent, S.; de Motta, J.C.B.; Sarthou, A.; Estivalezes, J.L.; Simonin, O.; Climent, E. A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows. *J. Computat. Phys.* **2014**, *256*, 582–614. [CrossRef]

10. Massol, A. Simulations Numériques d'écoulements à Travers des réseaux Fixes de Spheres Monodisperses et Bidisperses, Pour des Nombres de Reynolds Modérés. Ph.D. Thesis, Institut National Polytechnique de Toulouse, Toulouse, France, 2004. Available online: https://www.researchgate.net/publication/242075914_Adaptive_Augmented_Lagrangian_techniques_for_simulating_unsteady_multiphase_flows (accessed on 9 September 2021).

11. Schönfeld, T.; Rudyard, M. Steady and Unsteady Flow Simulations Using the Hybrid Flow Solver AVBP. *AIAA J.* **1999**, *37*, 1378–1385. [CrossRef]

12. Deen, N.G.; Kriebitzsch, S.H.; van der Hoef, M.A.; Kuipers, J. Direct numerical simulation of flow and heat transfer in dense fluid-particle systems. *Chem. Eng. Sci.* **2012**, *81*, 329–344. [CrossRef]

13. Tavassoli, H.; Kriebitzsch, S.; van der Hoef, M.; Peters, E.; Kuipers, J. Direct numerical simulation of particulate flow with heat transfer. *Int. J. Multiph. Flow* **2013**, *57*, 29–37. [CrossRef]

14. Gunn, D. Transfer of heat or mass to particles in fixed and fluidized beds. *Int. J. Heat Mass Transf.* **1978**, *21*, 467–476. [CrossRef]

15. Tenneti, S.; Sun, B.; Garg, R.; Subramaniam, S. Role of fluid heating in dense gas-solid flow as revealed by particle-resolved direct numerical simulation. *Int. J. Heat Mass Transf.* **2013**, *58*, 471–479. [CrossRef]

16. Deen, N.G.; Peters, E.A.J.F.; Padding, J.T.; Kuipers, J.A.M. Review of direct numerical simulation of fluid-particle mass, momentum and heat transfer in dense gas-solid flows. *Chem. Eng. Sci.* **2014**, *116*, 710–724. [CrossRef]

17. Sun, B.; Tenneti, S.; Subramaniam, S. Modeling average gas-solid heat transfer using particle-resolved direct numerical simulation. *Int. J. Heat Mass Transf.* **2015**, *86*, 898–913. [CrossRef]

18. Thiam, E.I.; Masi, E.; Climent, E.; Simonin, O.; Vincent, S. Particle-resolved numerical simulations of the gas-solid heat transfer in arrays of random motionless particles. *Acta Mech.* **2019**, *230*, 541–567. [CrossRef]

19. Chadil, M.A.; Vincent, S.; Estivalèzes, J.L. Accurate estimate of drag forces using particle-resolved direct numerical simulations. *Acta Mech.* **2018**, *230*, 569–595. [CrossRef]

20. Kataoka, I. Local instant formulation of two-phase flow. *Int. J. Multiph. Flow* **1986**, *12*, 745–758. [CrossRef]

21. Kataoka, I.; Ishii, M.; Serizawa, A. Local formulation and measurements of interfacial area concentration in two-phase flow. *Int. J. Multiph. Flow* **1986**, *12*, 505–529. [CrossRef]
22. Chadil, M.A.; Vincent, S.; Estivalezes, J.L. Improvement of the Viscous Penalty Method for Particle-Resolved Simulations. *Open J. Fluid Dyn.* **2019**, *09*, 168–192. [CrossRef]

23. Glowinski, R.; Pan, T.W.; Hesla, T.I.; Joseph, D.; Périaux, J. A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow. *J. Comput. Phys.* **2001**, *169*, 363–426. [CrossRef]

24. Khadra, K.; Angot, P.; Parneix, S.; Caltagirone, J. Fictitious domain approach for numerical modelling of Navier-Stokes equations. *Int. J. Numer. Methods Fluids* **2000**, *34*, 651–684. [CrossRef]

25. Hirt, C.W.; Nichols, B.D. Volume of fluid (VOF) method for the dynamics of free boundaries. *J. Comput. Phys.* **1981**, *39*, 201–225. [CrossRef]

26. Caltagirone, J.P.; Vincent, S. Tensorial penalisation method for solving Navier–Stokes equations. *Comptes Rendus De L'académie Des Sci. Paris Série IIb* **2001**, *329*, 607–613.

27. Vincent, S.; Larocque, J.; Lubin, P.; Caltagirone, J.P.; Pianet, G.; Climent, E. Adaptive Augmented Lagrangian Techniques for Simulating Unsteady Multiphase Flows. In Proceedings of the 6th International Conference on Multiphase Flow (ICMF 2007), Leipzig, Germany, 9–13 July 2007.

28. Gustafsson, I. *On First- and Second-Order Symmetric Factorisation Methods for the Solution of Elliptic Difference Equations*; Chalmers University of Technology: Gothenburg, Sweden, 1978.

29. der Vost, H.V. Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of non-symmetric systems. *SIAM J. Sci. Comput.* **1992**, *33*, 631–644. [CrossRef]

30. Schiller, L.; Naumann, A.Z. Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung. *Ver. Deut. Ing.* **1933**, 77, 318–320.

31. Ranz, W.E.; Marshall, W.R. Evaporation from drops. Parts I & II. *Chem. Eng. Progr.* **1952**, *48*, 141–146.

32. Feng, Z.G.; Michaelides, E.E. A numerical study on the transient heat transfer from a sphere at high Reynolds and Peclet numbers. *Int. J. Heat Mass Transf.* **2000**, *43*, 219–229. [CrossRef]

33. Whitaker, S. Forced convection heat transfer correlations for flow in pipes, past flat plates, single cylinders, single spheres, and for flow in packed beds and tube bundles. *AIChE J.* **1972**, *18*, 361–371. [CrossRef]
22. Chadil, M.A.; Vincent, S.; Estivalezes, J.L. Improvement of the Viscous Penalty Method for Particle-Resolved Simulations. *Open J. Fluid Dyn.* **2019**, *09*, 168–192. [CrossRef]

23. Glowinski, R.; Pan, T.W.; Hesla, T.I.; Joseph, D.; Périaux, J. A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow. *J. Comput. Phys.* **2001**, *169*, 363–426. [CrossRef]

24. Khadra, K.; Angot, P.; Parneix, S.; Caltagirone, J. Fictitious domain approach for numerical modelling of Navier-Stokes equations. *Int. J. Numer. Methods Fluids* **2000**, *34*, 651–684. [CrossRef]

25. Hirt, C.W.; Nichols, B.D. Volume of fluid (VOF) method for the dynamics of free boundaries. *J. Comput. Phys.* **1981**, *39*, 201–225. [CrossRef]

26. Caltagirone, J.P.; Vincent, S. Tensorial penalisation method for solving Navier–Stokes equations. *Comptes Rendus De L'académie Des Sci. Paris Série IIb* **2001**, *329*, 607–613.

27. Vincent, S.; Larocque, J.; Lubin, P.; Caltagirone, J.P.; Pianet, G.; Climent, E. Adaptive Augmented Lagrangian Techniques for Simulating Unsteady Multiphase Flows. In Proceedings of the 6th International Conference on Multiphase Flow (ICMF 2007), Leipzig, Germany, 9–13 July 2007.

28. Gustafsson, I. *On First- and Second-Order Symmetric Factorisation Methods for the Solution of Elliptic Difference Equations*; Chalmers University of Technology: Gothenburg, Sweden, 1978.

29. der Vost, H.V. Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of non-symmetric systems. *SIAM J. Sci. Comput.* **1992**, *33*, 631–644. [CrossRef]

30. Schiller, L.; Naumann, A.Z. Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung. *Ver. Deut. Ing.* **1933**, 77, 318–320.

31. Ranz, W.E.; Marshall, W.R. Evaporation from drops. Parts I & II. *Chem. Eng. Progr.* **1952**, *48*, 141–146.

32. Feng, Z.G.; Michaelides, E.E. A numerical study on the transient heat transfer from a sphere at high Reynolds and Peclet numbers. *Int. J. Heat Mass Transf.* **2000**, *43*, 219–229. [CrossRef]

33. Whitaker, S. Forced convection heat transfer correlations for flow in pipes, past flat plates, single cylinders, single spheres, and for flow in packed beds and tube bundles. *AIChE J.* **1972**, *18*, 361–371. [CrossRef]