# Investigation of heat transfers with particle-resolved simulations: From Stokes flow to fluidized bed

E. Butaye, Rémy Quintana, S. Mer, F Bataille, A. Toutant

► To cite this version:

E. Butaye, Rémy Quintana, S. Mer, F Bataille, A. Toutant. Investigation of heat transfers with particle-resolved simulations: From Stokes flow to fluidized bed. International Journal of Heat and Mass Transfer, 2025, 241, pp.126687. 10.1016/j.ijheatmasstransfer.2025.126687. hal-04915288

HAL Id: hal-04915288

https://hal.science/hal-04915288v1

Submitted on 27 Jan 2025

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
<figure><img src="image_2.png" /></figure>

# Investigation of heat transfers with particle-resolved simulations: From Stokes flow to fluidized bed

E. Butaye *, R. Quintana, S. Mer, F. Bataille, A. Toutant *

PROMES Laboratory (UPR 8521), CNRS – University of Perpignan (UPVD), 66100 Perpignan, France

## ARTICLE INFO

### Keywords:

- Fluid - particles flow

- Fluidized bed

- DEM

- Fully resolved particles

- PRS

- Viscous penalization

- Wall-to-bed heat transfer

- Anisothermal flow

## ABSTRACT

Particle-Resolved Simulations (PRS) of fluid-solid particles are conducted to study fluid-particle heat transfers and wall-to-bed heat transfers in an isothermal liquid-solid fluidized bed. An overview of existing PRS methods to study anisothermal fluid-solid flows is presented. In the framework of fluidized bed simulations, the collision detection method is optimized using Verlet tables. The overall computation time is reduced by 20%. An original Lagrangian method to compute the fluid-particle heat transfers is presented for an isolate particle. A parametric study on the fluid-particle heat transfer is performed to assess well-known correlations of the literature for a settling particle in a quiescent fluid. 77 PRS are performed for Reynolds numbers between 1 and 32 and Prandtl number between 0.1 and 10 for a grid resolution of 20 meshes per particle diameter. For two of three correlations considered, the predicted heat flow is within 10% error. An anisothermal fluidized bed of 2134 particles is finally studied. Four fluidization velocities are considered for a solid fraction comprised between 0.13 and 0.35. Three grid resolutions are carried out to assess the sensitivity of the mesh for the lowest fluidization velocity (12, 24 and 36 meshes per particle diameter). Results show that the macroscopic behavior of the bed is well retrieved even with a coarser grid as the solid fraction is well predicted. However, strong effects of the grid resolution are observed on the fluid-particle and the wall-to-bed heat transfers. The study of the velocity-temperature correlation shows that the parietal heat transfer is driven by the turbulent heat flow near the wall ($x^+ \sim 30$).

## 1. Introduction

The concept of fluidization is studied since a hundred years. In a recent review, Zhang et al. [1] highlighted that the first prototype of fluidization was introduced by Winkler et al. [2] in 1922 for water-gas production. Since then, it has been used in many industrial areas such as pneumatic conveying for powder transport [3] or methane cracking for hydrogen production with fluidized beds [4]. In the field of solar energy, fluidized beds are used as a new heat transfer fluid (HTF) in concentrated solar power (CSP) plants. In solar tower power plants, a panel of heliostats reflects solar radiation to a receiving surface located at the top of a solar tower. The HTF is exposed to the concentrated solar flux at the receiving surface so it received the solar flux. The HTF is then either transported in a thermal storage or injected into a thermodynamic Rankine cycle to convert its thermal energy into electricity [5]. Current CSP plants use molten salt as HTF. However, the main drawback is the limitation of the working temperature (<560°C) due to the degradation of thermodynamic properties at high temperature as a result of chemical reactions [5-7]. To overcome this issue, the use of solid particles as a new HTF is of great interest.

In particle-in-tube solar receiver, introduced by Flamant in 1982 [8], solid particles are fluidized in a dispenser with a primary air injection. With an imposed pressure gradient, the particles flow upward the tube. Finally, a secondary injection directly into the tube controls the mass flow rate and thus the fluidization regime. [7]. In a recent study, Gueguen et al. [7] succeeded in operating at 650°C demonstrating a 50% increase of the wall-to-bed heat transfer in front of a mean operating temperature of the particles of 500°C (see Fig. 13 in [7]). With such systems, the main limitation is no longer the temperature of the HTF but the solar receiver wall temperature. When very high temperatures are reached, the inhomogeneous solid fraction at the wall leads to a highly inhomogeneous heat transfer between the wall and the bed, resulting in high thermomechanical stresses. It is therefore required to characterize the flow inside the solar receiver for two reasons. First, to predict the wall-to-bed heat transfer as a function of the flow regime and physical properties of the fluid and the particles. It implies to understand the complex dynamic and thermal behavior of the fluid-particles flow. Then, to limit the thermomechanical stresses exerted on the wall of the solar receiver.

<footnote>* Corresponding authors.
E-mail addresses: edouard.butaye@promes.cnrs.fr (E. Butaye), adrien.toutant@univ-perp.fr (A. Toutant).</footnote>
Nomenclature

<table><tr><td>α</td><td>Volumic presence rate [–]</td></tr><tr><td>Δt</td><td>Temporal discretization [s]</td></tr><tr><td>Δx</td><td>Spatial discretization [m]</td></tr><tr><td>δn</td><td>Overlapping during collision [m]</td></tr><tr><td>λ</td><td>Thermal conductivity [W m<sup>-1</sup> K<sup>-1</sup>]</td></tr><tr><td>Fc</td><td>Vector collision force [N m<sup>-3</sup>]</td></tr><tr><td>g</td><td>Vector gravity acceleration [m s<sup>-2</sup>]</td></tr><tr><td>u</td><td>Vector velocity [m s<sup>-1</sup>]</td></tr><tr><td>μ</td><td>Dynamic viscosity [Pa s]</td></tr><tr><td>Φ</td><td>Heat flux [W]</td></tr><tr><td>ρ</td><td>Density [kg m<sup>-3</sup>]</td></tr><tr><td>τ0</td><td>Collision time [s]</td></tr><tr><td>f</td><td>Fluid</td></tr><tr><td>p</td><td>Particle</td></tr><tr><td>w</td><td>Wall</td></tr><tr><td>c<sub>p</sub></td><td>Heat capacity [J K<sup>-1</sup> kg<sup>-1</sup>]</td></tr><tr><td>D</td><td>Diameter [m]</td></tr><tr><td>e<sub>dry</sub></td><td>Dry restitution coefficient [–]</td></tr><tr><td>e<sub>eff</sub></td><td>Effective restitution coefficient [–]</td></tr><tr><td>e<sub>wet</sub></td><td>Wet restitution coefficient [–]</td></tr><tr><td>h</td><td>Heat transfer coefficient [W m<sup>-2</sup> K<sup>-1</sup>]</td></tr><tr><td>I</td><td>Volume phase indicator function [–]</td></tr><tr><td>k</td><td>Spring stiffness [N m<sup>-1</sup>]</td></tr><tr><td>m</td><td>Mass [kg]</td></tr><tr><td>N<sub>f</sub></td><td>Facets number [–]</td></tr><tr><td>N<sub>p</sub></td><td>Particles number [–]</td></tr><tr><td>Nu</td><td>Nusselt number [–]</td></tr><tr><td>p</td><td>Pressure [Pa]</td></tr><tr><td>Pe</td><td>Peclet number [–]</td></tr><tr><td>Pr</td><td>Prandtl number [–]</td></tr><tr><td>Re</td><td>Reynolds number [–]</td></tr><tr><td>S</td><td>Surface area [m<sup>2</sup>]</td></tr><tr><td>St</td><td>Stokes number [–]</td></tr><tr><td>T</td><td>Temperature [K]</td></tr><tr><td>T'</td><td>Fluctuating temperature [K]</td></tr><tr><td>U'</td><td>Fluctuating velocity [m s<sup>-1</sup>]</td></tr><tr><td>CFD-DEM</td><td>Computational Fluid Dynamics - Discrete Element Method</td></tr><tr><td>CSP</td><td>Concentrated Solar Power</td></tr><tr><td>HTF</td><td>Heat Transfer Fluid</td></tr><tr><td>PR-DNS</td><td>Particle Resolved - Direct Numerical Simulation</td></tr><tr><td>PRS</td><td>Particle Resolved Simulation</td></tr></table>

Over the last decades, a number of experiments were conducted to study fluidized beds [9,10]. Flow regime transitions are now well characterized in terms of particle shape [11] and hydrodynamics [12, 13]. Regarding heat transfers, Gueguen et al. showed that the wall-to-bed heat transfer increases with temperature and strongly depend on the fluidization regime. The aforementioned experimental studies have described the dynamic and thermal behavior of fluidized beds macroscopically. However, a local description of the flow is not possible through an experimental process without invasive methods and nearly impossible at high temperature. The prediction of the inhomogeneous behavior of parietal heat transfer is required to understand the thermomechanical stresses exerted on the wall of the solar receiver. The greater the particle agitation, the greater the heat transfer. Since the regime of interest for solar applications is turbulent fluidization [7], the

inhomogeneous nature of heat transfer cannot be neglected. Therefore, numerical simulation is an excellent tool to fully characterize the flow inside the solar receiver. It enables to access local quantities of the fluid and particles such as physical properties, mean value and fluctuations of velocity, pressure and temperature, or the local solid fraction.

Numerical simulations of fluidized beds are challenging because of the multi-scale phenomenon occurring in the flow. In a particle-in-tube fluidized beds, the largest structures can be of the order of the meter and yet strongly depend on fluid-particle and particle-particle interactions at the particle scale [14]. At the particle scale - also referred to as microscopic scale-, the hydrodynamic force exerted by the fluid to the particles and particle-particle interactions operate. At the scale of several particle diameters (mesoscopic scale), the wake of the particle and the interaction with the wake of other particles must be considered. At a larger scale (macroscopic scale), particle assembly movements prevail, resulting in different flow regimes which depend on the fluidization velocity and the particle size. Different scales of resolution have emerged to study fluidized beds [1,14,15]. At the smallest scale, particle-resolved direct numerical simulation (PR-DNS) fully resolve the flow at the particle scale and thus capture all the aforementioned hydrodynamic interactions. This requires to capture the velocity and temperature gradients at the interface between the fluid and the particle. For a review of PR-DNS methods, the reader is referred to Tenneti and Subramanian [16], Maxey [17] and more recently Marchelli et al. [18]. The Eulerian grid resolution associated with PR-DNS is then more than 40 meshes per diameter for viscous flows and increase with the Reynolds number [19]. The numerical cost of PR-DNS limits its application to academic cases, thus excluding the possibility to study the collective effects of a large number of particles involving several hundreds of millions of particles. To study a large number of particles, Computational Fluid Dynamics - Discrete Element Method (CFD-DEM) simulations, also commonly referred to as Eulerian-Lagrangian methods in the literature, can be performed [20–22]. At this scale, particles are smaller than the Eulerian grid size, enabling to study the collective effects of a particle assembly. CFD-DEM resolve Navier-Stokes equations on an Eulerian grid and track particles in a Lagrangian way. Particles are transported by solving Newton equation. As the flow is not resolved at the particle scale in that case, correlations are employed to model the hydrodynamic forces exerted by the fluid to the particles (drag force, lift force, Basset's force...) [23, 24]. If the solid fraction is sufficiently high, particles are not simply advected by the fluid (one-way coupling) and the force exerted by the particles to the fluid have to be considered (two-way coupling) as well as collision forces between particles (four-way coupling) [25]. CFD-DEM presents the advantage to track individually each particle and so to accurately compute the mean velocity and fluctuations of particles and solid fraction but models are required to consider mutual interactions exerted between the fluid and the particles. Closure models are obtained from PR-DNS simulations, as the hydrodynamic forces exerted by the fluid on the interface can be computed explicitly for each particle. Experience at a laboratory scale can be reproduced numerically. To study fluidized beds in an industrial context, the Euler-Euler method is most appropriate. It considered the fluid and the particles as two distinct and continuous phases [26,27]. Phase averaged Navier-Stokes equations are solved for each medium and interphase coupling is ensured via closure terms. CFD-DEM simulations can be harnessed to inform Euler-Euler simulations [28,29].

An intermediate scale between PR-DNS and CFD-DEM is Particle Resolved Simulation (PRS). PRS uses the same methodology as PR-DNS but with a coarser mesh because of the numerical costs associated with PR-DNS for a fluidized bed. Particles are therefore represented by a dozen meshes per diameter on the Eulerian mesh, which is not sufficient to fully resolve hydrodynamics and thermal gradients at the interface between the fluid and the particles. A first attempt to correct hydrodynamic gradients at the fluid-particle interface was presented in [19] for a single-particle Stokes flow. This correction will not be
studied or applied in the present work, as extensions of the method are needed to tackle the complex case of fluidized beds.

PRS are widely employed to compute drag force to inform Euler-Lagrange simulations [24,30,31]. However, due to high computational cost, the Eulerian mesh is coarse and hydrodynamic gradients are not fully resolved in the meaning of PR-DNS. Indeed, at least 5 Eulerian mesh cells in the boundary layer are required to accurately capture the gradients [30]. In the best cases, a mesh sensitivity is performed on macroscopic quantities such as the mean velocity, the drag force and the bed height [32]. In other scenarios, hydrodynamic forces were accurately computed on fixed beds, shallow fluidized beds, two-dimensional or very small three-dimensional fluidized beds [30,33–35]. To this day, the largest particle resolved simulation of a fluidized bed was performed by Dou et al. [36]. The authors computed a 3D particle resolved simulation of 115 200 particles in a shallow fluidized bed using a Lattice Boltzmann Method - Discrete Element Method (LBM-DEM) approach and GPU-based parallel computation. It is an impressive achievement for particle resolved simulation demonstrating the ability of their method to simulate multi-scale flows with relatively high accuracy (the lattice step was 20 times smaller than the particle diameter).

In the case of fluidized beds in solar receivers, heat transfers need to be considered. If numerical methods to study the hydrodynamic of fluidized beds with PRS are now well established, the consideration of heat transfers is much more recent. The first study to consider three-dimensional heat transfer on a single particle was conducted by Bagchi et al. [37] in 2001. Later, two-dimensional studies investigated heat transfers in a particle assembly. Feng and Michaelides [38] developed a numerical method to consider heat transfers with an Immersed Boundary Method (IBM) on a four-way coupling basis. Boussinesq hypothesis was assumed to consider density variation in the buoyancy term. Gan et al. [39] used the same hypothesis to study heat transfers with the Arbitrary Lagrangian Eulerian (ALE) finite-element method. Both authors demonstrated that the sedimentation process of particles is affected by the local temperature and buoyancy force. Gan et al. identified five sedimentation regimes defined by the Grashof number — which confronts buoyant force with viscous force. The authors demonstrated that the terminal velocity of a sedimenting particle depends on the Grashof number in a complex way. Yu et al. [40] employed the fictitious domain method and retrieved Gan results. It was later extended to three-dimensional simulations by Dan and Wachs [41] who performed the first three-dimensional study of a single particle in sedimentation in a semi-infinite channel.

Two types of heat transfers can be identified in particle-in-tube solar receiver: fluid-particle heat transfer and wall-to-bed heat transfer. In a review, Deen et al. highlighted the work of Tenneti et al. [42] and Tavassoli et al. [43] for the study of heat transfer in fixed array of spheres. Both authors computed the mean Nusselt number as a function of the void fraction to assess the validity of Gunn's correlation [44] up to Reynolds 100 for a Prandtl number of 0.72 [42] and 1 [43]. Deen et al. [45] refit Gunn's correlation based on the numerical data of the work of Tenneti and Tavassoli. In a following study, the authors extended their method to consider heat and mass transfer between the fluid and a fixed array of particles [46]. Three Reynolds numbers were studied (120, 180 and 240). The results are consistent with Gunn's correlation (within 4%), even if the resolution was rather coarse, with 24 meshes per particle diameter. Other studies investigated mass transfers in fixed and dynamic beds (fluidized beds with high inertia particles) [47], heat and mass transfers in a shallow fluidized bed [48]. Feng and Musong [48] pointed out that the simulation is partially under-resolved at Reynolds numbers above 200 when the resolution is only 15 meshes per particle diameter. Furthermore, the authors showed that the average particle Nusselt number increases with the fluidization velocity and the fluid-to-particle heat transfer decreases along the bed height. More recently, Chadil et al. [49,50] computed with a high-order interpolation method the fluid-to-particle heat transfer for each

particle in a fixed assembly. Finally, one of the most recent work on heat transfer in fluidized beds was pursued by Thiam [51]. The authors realized an extensive study on the heat transfer in a small fluidized bed (1120 particles) with coarse mesh (16 meshes per particle diameter). An excellent agreement was found with Wen & Yu [52] correlation for the drag force and with Gunn's [44] correlation for heat transfers. The results obtained for the resolution of 16 meshes per particle diameter were confronted with those obtained with a resolution of 32 meshes per diameter for Reynolds numbers comprised between 15 and 50 and for a Prandtl number equals to 7. No significant deviations were observed for the hydrodynamic behavior. However, the authors showed that a mesh refinement would have been required to correctly predict heat transfers because of the high Prandtl number. In their work, periodic boundary conditions were applied so parietal heat transfers were not considered.

Wall-to-bed heat transfers were studied by Haid [53]. The author compared five correlations to predict the heat transfer coefficient on a database of more than 2500 experiments. As it gathers experiments with very different heater geometries, the mean relative error of all correlations lie in 67.2%. The best prediction is obtained with the correlation using 7 constants, for an error of 32.0%. For the simplest correlation, using 3 constants, the error is 34.4%. The author emphasizes that this correlation is a reasonable compromise between accuracy and complexity. More recent experiments on wall-to-bed heat transfers were conducted by Gueguen et al. [7] and Lee et al. [54].

To the best of our knowledge, there is no study of wall-to-bed and fluid-to-particle heat transfers, for a particle-in-tube fluidized bed with a PRS approach, in the literature. This study aims to clarify heat transfers in an anisothermal fluidized bed. The paper is organized as follows. First, the numerical and physical modeling are described in Section 2. Key aspects of the numerical analysis are presented in Section 3. Section 4 pertains a parametric study of the heat transfers on an isolate particle. Section 5 investigates fluid-particle and wall-to-bed heat transfers in an anisothermal liquid-solid fluidized bed. Finally, conclusions and outlooks are drawn in Section 6.

## 2. Physical and numerical modeling

The assumptions on which the model relies are described in Section 2.1. The simulations of the present study are performed with TrioCFD software. TrioCFD, previously named Trio\_U is developed by the French atomic energy commission (CEA) since 1993 [55]. Since then, it has been widely used to study heat transfers [56], turbulence [57] and bubbly flows [58], to name but a few. In the present work, TrioCFD is used to investigate anisothermal fluid-solid flows. The numerical modeling of the particle-resolved approach on the Eulerian grid is presented in Section 2.2. The choices of the discretization schemes and solvers are detailed in Section 2.3.

### 2.1. General hypothesis

The computation of anisothermal fluidized beds requires to resolve the flow of the fluid and the particles trajectory. The model is based on a one-fluid formulation of the Navier-Stokes equations and particles are model as a highly viscous fluid. The numerical model is based on the following assumptions:

(i) both phases are incompressible. The fluid velocity is very small in front of the speed of sound. All physical properties of fluid and particles are assumed constants.

(ii) particles are considered spherical and non-deformable. Particle-particle and wall-particle collision velocities are too small to induce a deformation. The imposed temperature is sufficiently far from the melting temperature of the particles for dilation to be neglected. Only monodisperse fluidized beds are considered in the present study so all particles have the same diameter.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Interpolation of Eulerian velocity on Lagrangian markers. $u_i^*$ represents the velocity of the i-eme Lagrangian marker. $\cdot$ : pressure and temperature nodes, ▲/▼ : y / x velocity nodes, $\cdot$ : Lagrangian marker.</figcaption></figure>

(iii) the temperature is constant inside the particles. A boundary condition imposes the constant temperature on their surfaces. The thermal conductivity of particles is several orders of magnitude larger than that of the fluid. Conduction can therefore be considered instantaneous in the particles. In aimed application, solid particles are in addition small (less than 100 µm). The Biot number is small enough to justify the hypothesis of homogeneous temperature.

(iv) the particles' temperature is constant in time and uniform for all the particles in the bed. The fluidized beds computed are not circulating. Particles are therefore considered as a thermal sink.

(v) radiative heat transfers are neglected in the present approach as well as contact thermal resistances as particles are maintained at a constant temperature (heat sink).

## 2.2. Particle-resolved numerical approach

Numerous studies have investigated the numerical modeling of PRS. One can identify, the fictitious domain methods [59], the lattice Boltzmann methods [60] and the Arbitrarian Lagrangian Eulerian - Finite Element Methods (ALE-FEM) [61]. The approach described below is a fictitious domain method. The TriocFD-1.9.1 one-fluid formulation of Navier-Stokes equation is employed to resolve the fluid flow (see Section 2.2.1). The surface of the particles is described with a Front-Tracking method [62]. Each particle is therefore represented by a set of Lagrangian markers (see Fig. 1). Hereinafter, Lagrangian facets will refer to the closed surface formed by three Lagrangian markers (in two-dimension, it refers to the segment formed by two Lagrangian markers). To enforce the solid-body motion, all Lagrangian markers of a particle are advected with the same velocity (see Section 2.2.2). Furthermore, spurious velocities in particles are damped with a viscous penalization. This method implicitly enforces fluid constraints on the interface and was early validated in [63,64]. A soft-sphere collision model was developed in [65] to tackle the collision process. The numerical approach was employed in [19,65] to study fluidized beds and to compute hydrodynamic forces exerted by the fluid to the particles.

### 2.2.1. Navier-Stokes equation and temperature energy equation

With the one-fluid formulation, the fluid and the particles movement are governed by the following mass, momentum and energy equations, respectively:

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

$$ \frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla (\rho \mathbf{u} \otimes \mathbf{u}) = -\nabla p + \nabla \cdot (\mu (\nabla \mathbf{u} + \nabla^T \mathbf{u})) $$

$$ + \rho g + \mathbf{F}_c \qquad (2) $$

$$ \frac{\partial (\rho C_p T)}{\partial t} + \nabla \cdot (\rho C_p \mathbf{u} T) = \nabla \cdot (\lambda \nabla T) \qquad (3) $$

$\mathbf{u}$ is the velocity, $\rho$ the pressure, $T$ the temperature, $t$ the time, $g$ the gravity vector, $\rho$ the density, $\mu$ the viscosity. The definition of $\rho$ and $\mu$ specific to the one-fluid formulation is defined in Section 2.2.3. A collision source term $\mathbf{F}_c$ is added to the classical formulation of the momentum equation to consider wall-particle and particle-particle interactions. $\mathbf{F}_c$ ensures the coupling between the Eulerian field on which Navier-Stokes equations are solved, and the Lagrangian approach employed to compute the collision forces. Velocity, pressure and temperature are solved on a staggered grid (see Fig. 1). Scalar quantities are computed at the center of Eulerian cells whereas vector quantities are computed component by component at the center of the faces of the element volume control (see Fig. 1).

Without phase change, there is no stress jump at the interface between the fluid and the particle [66,67]. Inside solid particles, the deformation tensor is null. This condition is achieved by imposing a high viscosity ratio between the fluid and the particles (viscous penalization).

For the energy equation, the boundary condition to impose is a Dirichlet condition ($T_p = 0$ K). The temperature is only solved for the fluid. As both phases are incompressible, physical properties are not temperature-dependent. Consequently, the value of 0 K does not refer to absolute zero, and only the temperature difference between the fluid and the solid should be considered. A ghost fluid method is employed to ensure it implicitly. The procedure is the following. First, the Lagrangian facets' normal is computed. Then, the temperature gradient is evaluated at the gravity center of all the Lagrangian facets by considering a pure diffusion model in the boundary layer of the particles. The gradient, computed in a Lagrangian manner, is distributed to purely fluid Eulerian cells close to the interface. A linear extrapolation is performed to extend the temperature gradient in solid cells. Finally, the convection and diffusion terms of the energy equation (see Eq. (3)) can be computed in the whole computational domain. A complete description of the procedure is detailed in Grosso et al. [56].

### 2.2.2. Transport of Lagrangian markers

Lagrangian markers that describe the particles' surface represent a mobile mesh overlying on the fixed Eulerian mesh. The interest of the method lies in its ability to represent accurately the position of the interface. The Front-Tracking (FT) method is often described as complex because of the remeshing procedures that are required for bubbles or drops. For solid particles – which are rigid by definition – no remeshing is required which drastically simplifies the FT algorithm. Nevertheless, it should be pointed out that, because of the data structure and the communications between Lagrangian and Eulerian meshes, the computation of the volume fraction in two-phase cells and the management of MPI operations for multiprocessor computing are rather complex.

The interface is advected by a transport equation:

$$ \frac{\partial \chi}{\partial t} + \mathbf{u} \cdot \nabla \chi = 0 \qquad (4) $$

where $\chi$ is the phase indicator function. $\chi = 0$ in the solid phase and $\chi = 1$ in the fluid phase (see Fig. 2).

To solve this equation numerically, a trilinear interpolation of the Eulerian velocity field is computed at each Lagrangian marker (see Eq. 1). The velocity is then computed at the gravity center of each Lagrangian marker. The velocity of a particle $p$ writes:

$$ \mathbf{u}_p = \frac{1}{S_p} \sum_{\text{facets}} S_i \mathbf{u}_i \qquad (5) $$

where $S_p$ is the particle surface, $S_i$ and $\mathbf{u}_i$ the surface and velocity of the facet $i$. The Lagrangian markers are transported with the velocity $\mathbf{u}_p$:

$$ \mathbf{x}_i^{t+1} = \mathbf{x}_i^t + \mathbf{u}_p \Delta t \qquad (6) $$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Phase indicator function.</figcaption></figure>

where $\mathbf{x}_i$ denotes the position of the Lagrangian marker $i$ and $\Delta t$ is the time step of the simulation. As mentioned in [19], the Lagrangian markers are not physical points of the interface. Consequently, if the particle is rotating, a given marker does not represent the same physical point at two different instants. This is possible because the particle is spherical. The projection of its shape onto the Eulerian mesh remains unchanged by pure rotation. The rotational motion of a particle is only considered when solving the Navier-Stokes equation, and would result from a local sheared flow. This phenomenon can be observed by looking at a cross-section of the velocity field inside the particle. The rotational motion of a particle cannot result directly from a collision. For non-spherical particles, the equation of angular momentum conservation must be resolved.

### 2.2.3. Defining the one-fluid properties

In the one-fluid formulation of Navier-Stokes equation, physical properties are defined by the presence rate of each phase in an Eulerian cell. In each Eulerian cell of volume $V$, the volume indicator function is defined as:

$$I = \int_V \chi dV \qquad (7)$$

The density is computed with an arithmetic average (see Eq. (8)) whereas the viscosity is computed with a harmonic average (see Eq. (9)).

$$\rho = I \rho_f + (1 - I) \rho_p \qquad (8)$$

$$\mu = \frac{\mu_f \mu_p}{I \mu_p + (1 - I) \mu_f} \qquad (9)$$

The thermal conductivity is only defined in purely fluid cells as temperature is not solved inside the particle ($T_p = 0$ K).

### 2.2.4. Collision modeling

The soft-sphere collision model of Hamidi et al. [65] is used to model solid-solid interactions. As the collision time is of the order $10^{-9}$ s, the resolution of solid-solid interaction with strain laws is impossible as it would drastically constrain the simulation time step. Thus, the collision is spread over several time steps of the fluid solver. Collision forces are computed using a Lagrangian approach. Each collision force exerted on a particle is discretized on the Eulerian field in the corresponding purely solid meshes. The coupling between the Eulerian approach and the Lagrangian approach is achieved by the number of the particle. Indeed, at each time step, a correspondence is performed between the Lagrangian and the Eulerian number of the particles. The model developed by Hamidi et al. [65] is briefly described hereinafter.

The contact force is modeled with a harmonic oscillator and writes:

$$F_c = -k \delta_n \mathbf{n} \qquad (10)$$

where $\delta_n$ is the particle-wall or particle-particle overlapping, $k$ is the spring stiffness constant and $\mathbf{n}$ is the normal vector to the plan of collision. $\delta_n$ is described by the following equation:

$$m_e \dot{\delta}_n + k \delta_n = 0 \qquad (11)$$

where $m_e$ is the effective mass. $m_e = \frac{m_p}{2}$ for a particle-particle collision and $m_e = m_p$ for a particle-wall collision.

The collision process is divided in two steps: the impact step corresponds to the decrease of the particle velocity during a collision while the rebound step starts when the particle's velocity changes sign. During the impact:

$$k = m_e \left( \frac{\pi}{\tau_0} \right)^2 \qquad (12)$$

where $\tau_0$ represents the collision time. During the rebound:

$$k = e_{eff}^2 m_e \left( \frac{\pi}{\tau_0} \right)^2 \qquad (13)$$

where $e_{eff}$ is the effective restitution coefficient and is another input parameter of the model. $e_{eff}$ is the ratio of the particle velocity after impact to the particle velocity before impact. $e_{eff} = e_d e_{wet}$. $e_d$ is the dry coefficient that models the energy loss during the collision by elastic deformation of the particle and vibration. $e_{wet}$ models the (unresolved) lubricating force that dissipates kinetic energy through viscosity. According to Legendre et al. [68]:

$$e_{wet} = \exp \left( -\frac{35}{S t} \right) \qquad (14)$$

The Stokes number is defined by $St = \frac{1}{9} \frac{\rho_n}{\rho_f} Re$ and the Reynolds number is defined by $Re = \frac{\rho U D_p}{\mu_f}$. The input parameters of the collision model are the dry coefficient and the collision time.

### 2.3. Discretization schemes and solvers

An explicit first-order Euler scheme is used to discretize the time derivatives of the mass, momentum and energy equations and to address Lagrangian markers. An implicit treatment of diffusion relaxes the time step. Spatial derivatives of convection and diffusion terms are discretized using a second-order centered scheme. A prediction correction algorithm is used to solve the momentum equation. The Poisson equation is then solved by the conjugate gradient method with a symmetrical successive overrelaxation preconditioner (ssor).

## 3. Numerical tool optimization and post-processing

PRS of fluidized beds is challenging because it requires to detect and model all collisions separately. The computation time to detect the particles in collisions increases with the square of the particle number because all particle pairs are tested. Thus, it limits the number of particles to be studied without a careful treatment of the algorithm of detection. This issue is addressed in Section 3.1. The assets of PRS-DEM approach lie in its ability to (i) track particles in their movement, (ii) compute the hydrodynamic force exerted by the fluid and the heat flux received on the particles surface. A description of the heat flux computing procedure is provided in Section 3.2. The post-processing of fluidized bed simulations is complex because of the large amount of data generated. The post-processing methodology is presented in Section 3.3.

### 3.1. Optimization of collision detection

One of the numerical challenges to fluidized beds computation is the detection of collision. The simplest algorithm to implement is to compute the distance to center for all particle pairs in the domain [65]. This procedure is realized by each CPU of the simulation and is therefore not optimized. Furthermore, the computational complexity scales
as $\mathcal{O}\left(\frac{N_p(N_p+1)}{2}\right)$. However, not all collisions are likely to occur, as two particles may be located several diameters apart. It is therefore essential to have efficient collision detection processing to handle several thousand or even tens of thousands of particles at reasonable cost.

An alternative procedure, namely the Verlet algorithm [69], is implemented in this paper. This method is widely employed for collision detection [64,70] in granular media. It is based on the Verlet table, which records the particles with which each particle may collide during a given time step. The algorithm is described in Appendix C.

This procedure is very similar to the one detailed by Fang et al. [70]. For all the simulations, $\epsilon_v$ is set to 30% of the particle diameter and $\Delta t^{max} = 50\Delta t_{simu}$ with $\Delta t_{simu}$ the computational time step. The computational efficiency of the approach is presented in for a fluidized bed of 2134 particles, three mesh resolutions (12, 24 and 36 meshes per particle diameter) and four solid fractions. One can see that for the coarser mesh, the computational time of the collision process is scaled by a factor comprised between 12.2 and 13.6. The lower the solid volume fraction, the more efficient the method as less collision occurs. For finest grid resolutions (24 and 36 meshes per particle diameter), one can note that the computational time of the collision process is comprised between 1.0% and 1.3% which shows the high scalability of the method. One should point out that the collision detection does not depend on the grid resolution. However, the finer the mesh, the longer the CPU time to complete a time step (due to increased CPU communications). As a result, the relative time devoted to collision detection over a time step decreases with grid resolution. For these resolutions, the simulations were only performed for the Verlet algorithm due to its computational cost.

For the coarser mesh, the overall simulation time is reduced by 20%. A parametric study on the distance $\epsilon_v$ and the time $\Delta t^{max}$ would be required to have the optimal pair of parameters. For the sake of simplicity, we propose to define $\Delta t^{max}$ as the time required for a particle to travel half the distance $\epsilon_v$ at the median velocity of all particles. Indeed, $\Delta t^{max}$ in only used as a safety factor. It must therefore be neither too small, as this would considerably increase computation time, nor too large, as some collisions might not be detected.

The order of the Verlet method is $\mathcal{O}\left(\frac{N_p(N_p+1)}{2}\right)$ every $\Delta t_{maj}$ time steps and $\mathcal{O}(N)$ otherwise. Some preliminary tests on a fluidized bed of 8536 particles showed that with Verlet algorithm, the computation time of the collision force represents 14% of the time it takes to complete a time step. Therefore, one suggests combining the Verlet Algorithm with the Link Cell (LC) method as proposed by Fang et al. [70]. The algorithm is currently under development and will not be presented in the present paper. The LC method was previously employed in the literature [64,71]. It consists to search for potential collision in sub-domains close to a given one. In [64] and [71], a sub-domain matches the domain of a CPU or a GPU. Therefore, potential collisions between long-distance particles are not computed.

The efficiency of collision detection does not necessarily justify the implementation of a complex algorithm, especially if only a few particles are being investigated. Furthermore, LC algorithms require communication between Central Processing Units (CPUs), which can also increase computation time. The choice of the algorithm therefore depends on the number of particles and the solid fraction of the bed. More precisely, it depends on the number of particles contained in each sub-domain and therefore on the choice of domain partitioning between CPU or GPU for a given solid fraction. The denser the bed, the larger the number of collisions to be handled. As it is difficult to define a criterion considering these two parameters, we propose the following classification:

* $N_p = \mathcal{O}(1) - \mathcal{O}(10)$: classical algorithm

* $N_p\mathcal{O}(10) - \mathcal{O}(10^3)$: Verlet algorithm [69]

* $N_p \ge \mathcal{O}(10^4)$: Fang algorithm [70]

<caption>**Table 1**</caption>

Collision force computation time as a percentage of the total time it takes to complete a time step. $\alpha$ is the solid volume fraction, $D_p$, $\Delta x$ are the particle diameter and the spatial discretization step, respectively. $\frac{D_p}{\Delta x}$ represents the spatial resolution as a function of the particle diameter.

<table><thead><tr><th></th><th>α</th><th>0.36</th><th>0.29</th><th>0.20</th><th>0.13</th></tr></thead><tbody><tr><td>Verlet</td><td>12</td><td>1.8</td><td>1.7</td><td>1.7</td><td>1.4</td></tr><tr><td>algorithm</td><td>24</td><td>1.3</td><td>1.1</td><td>1.1</td><td>-</td></tr><tr><td></td><td>36</td><td>1.0</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Classic<br>algorithm</td><td>12</td><td>22</td><td>22</td><td>21</td><td>19</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Diagram of the heat flux computing. Source: Figure extracted from [19].</figcaption></figure>

## 3.2. Fluid-particle heat flux postprocessing

The computation of the heat flux follows the method developed in [19] for the hydrodynamic force. The heat flux received by the particle writes:

$$ \Phi_p = \iint_{S_p} \lambda_f \nabla \Gamma \cdot n dS \qquad (15) $$

with $\lambda_f$ the thermal conductivity of the fluid, $n$ the normal to the interface and $S_p$ the particle surface. Eq. (15) is discretized on the Lagrangian mesh as follows:

$$ \Phi_p = \sum_{i}^{N_f} \lambda_f \nabla T_i \cdot n_i S_i \qquad (16) $$

where $N_f$ is the number of Lagrangian facets of a particle and $S_i$ is its surface area. For each Lagrangian facets, two interpolation points are defined at a distance $\delta$ and $2\delta$ along the normal (see Fig. 3). The distance $\delta$ is chosen equals to $\Delta x$ as defined in [19], where $\Delta x$ is the Eulerian mesh size. The temperature is evaluated in all points $P_1$ and $P_2$ with a tri-linear interpolation from the Eulerian mesh. The temperature gradient is computed with a second order upwind scheme as follows:

$$ \nabla T_i = \lambda \frac{-T_{P_2} + 4T_{P_1} - 3T_p}{2\delta} S_i \qquad (17) $$

where $T_p$ is the temperature of the particle. The reader is referred to [19] for a complete description of the interpolation procedure.

## 3.3. Post-processing methodology

The post-processing of the simulations of fluidized bed is rather complex because of the amount of data generated. For the finest
<figure><img src="image_6.png" /><figcaption>Fig. 4. Two-dimensional diagram of the configuration of the sedimentation of a particle in an infinite medium. The boundary conditions are as follows: -: Imposed pressure walls, -: Imposed velocity walls, -: Imposed temperature walls, -: Imposed temperature interface.</figcaption></figure>

mesh, containing 181 million meshes, more than 1 Terabyte of data was generated for 10 s of simulated flow. Post-processing such a large quantity of data requires specific treatment. A brief description of the procedure is detailed hereinafter for the Lagrangian statistics on one hand, and on the Eulerian statistics on the other hand. More details are also given in Appendices A and B.

Lagrangian statistics include the velocity and position of each particle, as well as the hydrodynamic force exerted by the fluid on the particle surface and the heat flux received by the particle. Data is recorded for all particles every 0.01 s of physical simulated time, which represents a good compromise between the accuracy of particle tracking, the simulation time spent writing the data and the computation time for post-processing the statistics. The supercomputer on which the computations were performed requires a maximum execution time of 24 h per simulation. Each complete simulation is therefore divided into 24 h simulation folders. Each folder contains a single file including the Lagrangian quantities of all the particles for each time step. For each folder, the velocity, position, hydrodynamic force and heat flux are extracted using bash commands and GNU Parallel [72]. Each variable is saved in an independent csv file. The data is then concatenated for all the folders of the simulation and saved in HDF format for faster rereading. Finally, the void fraction, bed height, position and velocity statistics are computed using Python multiprocessing tools.

Eulerian fields include the pressure, the velocity, the temperature and the solid volume fraction. The fields are written by TrioCFD software in a lata format at the same frequency as for the Lagrangian data. The files are first converted to the vtk format using GNU parallel [72] and Visit software. The Vtk files are loaded into Python using the module Pyvista [73]. The variables of interest are averaged for each 24 h simulation folder using Pyvista and Pandas modules coupled with Python multiprocessing tools and saved in HDF format. Finally, the data is concatenated for all the folders of the simulation and Eulerian statistics are computed. The computation of Eulerian averages on planes parallel and transverse to the flow is detailed in Appendix A. Finally, the wall heat transfer coefficient computation is explained in Appendix B.

## 4. Heat transfers on an isolated particle

The present method was validated in [19,65] for the hydrodynamic of the flow. A parametric study of the heat transfer computation is presented on the sedimentation of a particle at very small Reynolds

<table><caption>Table 2
Physical parameters of Stokes configuration.</caption><thead><tr><th>Re</th><th>Pr</th><th>ρ<sub>p</sub> /ρ<sub>l</sub></th><th>ρ<sub>v</sub> /ρ<sub>f</sub></th><th>D<sub>p</sub></th><th>ΔT</th></tr></thead><tbody><tr><td>7.8.10<sup>-2</sup></td><td>6</td><td>10</td><td>10<sup>3</sup></td><td>5.2.10<sup>-6</sup> m</td><td>5 K</td></tr></tbody></table>

under Stokes assumption (see Section 4.1) and on the sedimentation of a sphere in a quiescent fluid (see Section 4.2).

### 4.1. Sedimentation of a particle in an infinite medium with heat transfer

The aim of this section is to study temperature diffusion in a case where it is predominant compared to convection. For this purpose, the configuration described in [19] was reproduced and adapted to study heat transfer (see Fig. 4). A cubic domain of section 3 $D_p$ is considered. The theoretical velocity and pressure fields of Stokes sedimentation are imposed at the boundaries (see Tab. 1 in [19]). A fixed temperature $T_w$ is imposed on the lateral walls. The thermal boundary condition at the inlet is also an imposed temperature at $T_w$. The solid temperature is initialized at $T_p$ and remains constant during the simulation. The initial fluid initial is $T_w$. Due to the heat transfer with the particle, a boundary layer develops around the particle (see Fig. 5).

The heat flux received by the sphere is computed with two methods. The first one is described in Section 3.2 and the second relies on a heat balance at the boundaries of the domain. It writes:

$$ \Phi_p = \Phi_l + \Phi_o - \Phi_o \quad (18) $$

with respectively $\Phi_p$, $\Phi_l$, $\Phi_o$ and $\Phi_o$, the heat flux at the particle, the inlet, the lateral boundaries and the outlet. The Nusselt number is computed as follows:

$$ Nu_p = \frac{hD_p}{\lambda_f} \quad (19) $$

and the heat transfer coefficient writes:

$$ h = \frac{\Phi_p}{S_p(T_\infty - T_p)} \quad (20) $$

with $S_p$, the particle surface and $T_\infty$ the undisturbed temperature of the fluid ($T_\infty = T_w$).

The physical and numerical parameters of the simulation are summarized in Table 2. A mesh convergence study of the particle Nusslet number is presented in Table 3 for the method which uses the heat balance at the boundaries. This method presents the advantage of providing the flux actually received by the particle during the simulation, in contrast to the first method described in Section 3.2. The relative error to the finest mesh, $D_p/\Delta x = 50$ which is considered as the exact solution, is 3.5% for the coarser mesh while it drops down to 0.17% for a mesh resolution of $D_p/\Delta x = 20$. Snapshots of the steady-state temperature field are presented in Fig. 5 for two mesh resolutions (5 and 40 meshes per diameter). Even with the coarser mesh, the boundary layer of the particle is well captured. The particle Nusselt number computed with the finest mesh is compared to well-known correlations of the literature listed below:

* Ranz and Marshall [74]

$$ Nu = 2 + 0.6 Pr^{1/3} Re^{1/2} \quad (21) $$
for $2 \le Re \le 10^4$, $Pr \ge 0.6$

* Whitaker [75]:

$$ Nu = 2 + (0.4 Re^{1/2} + 0.06 Re^{2/3}) Pr^{2/5} \quad (22) $$
for $0.71 \le Pr \le 380$, $3.5 \le Re \le 7.6 \times 10^4$

* Feng and Michaelides [76]

$$ Nu = 0.992 + Pe^{1/3} + 0.1 Pe^{1/3} Re^{1/3} \quad (23) $$
for $0.1 \le Re \le 4000$, $0.2 \le Pe \le 2000$.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Temperature field of a Stokes flow past a sphere.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Mesh convergence of the Nusselt number computed with two methods. x: heat balance at the boundaries, x: interpolation method.</figcaption></figure>

<caption>Table 3</caption>

Particle Nusselt number computed with the heat flux obtained with a heat balance at the boundaries.

<table><caption>Table 4</caption><tr><td>D<sub>p</sub>/Δx</td><td>5</td><td>8</td><td>10</td><td>12</td><td>16</td></tr><tr><td>Nu<sub>p</sub></td><td>2.647</td><td>2.704</td><td>2.716</td><td>2.738</td><td>2.736</td></tr><tr><td>D<sub>p</sub>/Δx</td><td>20</td><td>24</td><td>30</td><td>40</td><td></td></tr><tr><td>Nu<sub>p</sub></td><td>2.740</td><td>2.740</td><td>2.742</td><td>2.742</td><td></td></tr></table>

Relative error of the prediction of the Nusselt number compared to the value computed with PR-DNS.

<table><thead><tr><td>Reference value</td><td>Nusselt number</td><td>Relative error</td></tr></thead><tbody><tr><td>PR-DNS</td><td>2.74</td><td>–</td></tr><tr><td>Ranz and Marshall</td><td>2.30</td><td>16.1%</td></tr><tr><td>Whitaker</td><td>2.25</td><td>17.8%</td></tr><tr><td>Feng and Michaelides</td><td>1.80</td><td>34.3%</td></tr></tbody></table>

The relative error of the correlations to the computed value with PR-DNS is presented in Table 4. The error is larger than 15% for all correlations. It can be explained by the fact that the Reynolds number of the simulation is too low and the correlations are not adapted for this flow regime. In addition, the error could be related to a confinement effect, as the domain is small compared to the particle diameter. The fixed temperature imposed on the boundaries may also alter the shape of the thermal boundary layer.

The heat balance method for calculating the flux received by the particle is only possible because the calculation domain contains only one particle. For multi-particle system, the use of the interpolation method (see Section 3.2) is required even though for resolution or collision reasons, the calculated flux might deviate from the exact heat

<caption>Table 5</caption>

Fluid properties of Cate's et al. experiment [77].

<table><thead><tr><td>Case number</td><td>ρ<sub>f</sub> [kg m<sup>-3</sup>]</td><td>μ<sub>f</sub> [Pa s]</td><td>Re</td></tr></thead><tbody><tr><td>E1</td><td>970</td><td>0.373</td><td>1.50</td></tr><tr><td>E2</td><td>965</td><td>0.212</td><td>4.10</td></tr><tr><td>E3</td><td>962</td><td>0.113</td><td>11.6</td></tr><tr><td>E4</td><td>960</td><td>0.058</td><td>31.9</td></tr></tbody></table>

flux received by the particle. A mesh convergence of the heat flux computed with this method is presented in Fig. 6. The interpolation method and the heat balance method both converge to the same value. The relative deviation falls below 4% for resolutions greater than 20 meshes per particle diameter. The interpolation method is therefore validated and will be employed to compute the flux received by the particle in the next section.

## 4.2. Settling sphere in a quiescent viscous fluid with heat transfer

A particle in sedimentation in a quiescent fluid is considered to assess the method at higher Reynolds numbers. The isothermal case of Cate's et al. experiment [77] is adapted to study fluid-particle heat transfer. The domain is a rectangular parallelepiped with square cross-section, sides 6.67 $D_p$ and height 10.67 $D_p$. The particle diameter is 15 mm. The fluid properties are reported in Table 5. This configuration was previously studied by [65] with the same method without considering the heat transfer. The position and velocity profiles were in good agreement with the experimental data even with a coarse mesh of 15 meshes per particle diameter.

### 4.2.1. Assessment of the sedimentation trajectory

A fixed wall condition is applied on lateral and bottom boundaries. A pressure imposed condition is applied at the top. The temperature is imposed at 10 K at the wall and at 0 K at the particle surface. Initially, the fluid is at 10 K. Three mesh resolutions are investigated: 20, 30 and 40 meshes per particle diameter. It required respectively 200, 490 and 1088 CPUs. Four fluid viscosities were considered by Cate et al. to compute different flow regimes. Reynolds numbers based on the terminal settling velocity – in an infinite medium – are 1.50, 4.10, 11.6 and 31.9. The temporal evolution of the vertical position and velocity of the particle are represented in Fig. 7. It should be noted that the simulation was terminated before the particle reached the bottom of the tank in cases where the Reynolds number is below 4.10. For all Reynolds number, the numerical results are converged with the mesh resolution of 20 meshes per diameter, as all curves completely overlap. This avoided excessive computation costs when the simulated trajectory already correctly predicted the experiment. The curves representing the different meshes are superimposed. It shows that the mesh resolution of 15 meshes per particle diameter of Hamidi et al. [65] was too coarse to fully capture the hydrodynamic force exerted by the fluid
<figure><img src="image_9.png" /><figcaption>Fig. 7. Temporal evolution of the particle vertical position (left) and velocity (right). $\cdot: \frac{D_c}{\Delta x} = 20$, $\cdot: \frac{D_c}{\Delta x} = 30$, $\cdot: \frac{D_c}{\Delta x} = 40$, $\bullet$: experimental data [77]. $\cdot: Re = 1.50$, $\cdot: Re = 4.10$, $\cdot: Re = 11.6$, $\cdot: Re = 31.9$. For a given Reynolds number, all mesh resolutions are superimposed.</figcaption></figure>

<caption>Table 6</caption>

Nusselt number for various Reynolds and Prandtl number. $Nu^{RM}$, $Nu^W$ and $Nu^{FM}$ are the Nusselt number computed with the correlation of Ranz and Marshall [74], Whitaker [75] and Feng and Michaelides [76], respectively. $Nu^{PR-DNS}$ is the simulated value of the Nusselt number. $\epsilon$ is the relative error of the correlation to the simulated value: $\epsilon = (Nu_{p}^{\text{correlation}} - Nu_{p}^{\text{PRS}})/Nu_{p}^{\text{PRS}}$.

<table><thead><tr><th>Pe</th><th>Re</th><th>Pr</th><th>Nu<sub>p</sub><sup>PRS</sup></th><th>Nu<sub>p</sub><sup>RM</sup></th><th>Nu<sub>p</sub><sup>W</sup></th><th>Nu<sub>p</sub><sup>FM</sup></th><th>ε<sup>RM</sup> [%]</th><th>ε<sup>W</sup> [%]</th><th>ε<sup>FM</sup> [%]</th></tr></thead><tbody><tr><td>0.67</td><td>1.33</td><td>0.5</td><td>2.23</td><td>2.55</td><td>2.40</td><td>1.96</td><td>14.3</td><td>7.86</td><td>-12.0</td></tr><tr><td>1.33</td><td>1.33</td><td>1</td><td>2.35</td><td>2.69</td><td>2.53</td><td>2.21</td><td>14.5</td><td>7.80</td><td>-5.86</td></tr><tr><td>1.94</td><td>3.88</td><td>0.5</td><td>2.47</td><td>2.94</td><td>2.71</td><td>2.44</td><td>18.8</td><td>9.54</td><td>-1.55</td></tr><tr><td>2.66</td><td>1.33</td><td>2</td><td>2.58</td><td>2.87</td><td>2.70</td><td>2.53</td><td>11.2</td><td>4.71</td><td>-2.04</td></tr><tr><td>3.88</td><td>3.88</td><td>1</td><td>2.77</td><td>3.18</td><td>2.94</td><td>2.81</td><td>14.7</td><td>5.87</td><td>1.33</td></tr><tr><td>5.50</td><td>11.0</td><td>0.5</td><td>3.04</td><td>3.58</td><td>3.23</td><td>3.15</td><td>17.8</td><td>6.31</td><td>3.68</td></tr><tr><td>7.76</td><td>3.88</td><td>2</td><td>3.18</td><td>3.49</td><td>3.24</td><td>3.28</td><td>9.74</td><td>1.76</td><td>3.26</td></tr><tr><td>11.0</td><td>11.0</td><td>1</td><td>3.55</td><td>3.99</td><td>3.63</td><td>3.71</td><td>12.4</td><td>2.09</td><td>4.56</td></tr><tr><td>15.0</td><td>29.9</td><td>0.5</td><td>4.01</td><td>4.61</td><td>4.10</td><td>4.22</td><td>14.7</td><td>2.04</td><td>5.15</td></tr><tr><td>22.0</td><td>11.0</td><td>2</td><td>4.21</td><td>4.51</td><td>4.15</td><td>4.42</td><td>7.12</td><td>-1.56</td><td>4.99</td></tr><tr><td>29.9</td><td>29.9</td><td>1</td><td>4.84</td><td>5.28</td><td>4.77</td><td>5.06</td><td>9.11</td><td>-1.54</td><td>4.53</td></tr><tr><td>59.8</td><td>29.9</td><td>2</td><td>5.87</td><td>6.14</td><td>5.65</td><td>6.12</td><td>4.50</td><td>-3.76</td><td>4.20</td></tr></tbody></table>

to the particle. This further explains the slight discrepancy with the experiment observed by the authors on particle velocity at the highest sedimentation velocity.

### 4.2.2. Assessment of fluid-particle heat transfer

In this section, the Nusselt number is computed with the method based on the interpolation from the nearest Eulerian fluid points. For each fluid viscosity considered in the previous section, three Prandtl numbers were studied ($Pr \in \{0.5, 1, 2\}$). In each case, one can compute the Nusselt number as a function of the time. However, the method required the boundary layer to be fully developed for the interpolation of the gradient to be accurate. Therefore, the transient phase from the start of the simulation cannot be studied. In each case, one can compute the Nusselt number of the particle when it has reached its terminal velocity. The data are presented in Table 6. The minimum relative error of the three correlations to the computed value is below 8% for all the computed cases. Ranz and Marshall correlation [74] always overestimates the particle Nusselt number and appeared less adapted than the correlations developed by Whitaker [75] and Feng and Michaelides [76]. One cannot define the most appropriate correlation between Whitaker and Feng and Michaelides, as neither appears to predict more accurately the Nusselt number than the other.

The parametric study is extended to a wider range of Peclet numbers ($Pe = RePr = D_pU_t/\alpha$, with $U_t$ the particle terminal velocity and $\alpha$ the thermal diffusivity). 77 simulations with a resolution of 20 meshes per particle diameter. 7 Reynolds numbers and 11 Prandtl numbers were considered. Reynolds numbers are defined by the viscosity of the fluid and correspond to the Reynolds number reached by the particle when it reaches its terminal sedimentation velocity ($Re \in \{1, 2, 3, 4, 5, 12, 32\}$). Prandtl numbers are defined by the thermal conductivity of the fluid

($Pr \in \{0.1, 0.2, 0.25, 0.33, 0.5, 1, 2, 3, 4, 5, 10\}$). In the context of flows in solar receivers, Reynolds numbers are of the order of unity and the Prandtl number is approximately 0.7 for air. Peclet numbers are therefore very low. The parametric study focuses on these flow regimes. In each case, one value of the particle Nusselt number is computed when the thermal boundary layer is fully developed. The value given by the three correlations is computed with the corresponding Reynolds and Prandtl numbers. One can plot the Nusselt number obtained with the correlations as a function of the Nusselt number computed with PR-DNS (see Fig. 8). The closer the points to the identity curve, the more accurate the correlation. The correlations of Whitaker [75] and Feng and Michaelides [76] are both within the 10% confidence interval of the PR-DNS values. This demonstrates that they are both appropriate for predicting heat transfer between the fluid and an isolated particle at low Peclet numbers.

An attempt was made to optimize the coefficients of Whitaker [75] and Feng and Michaelides correlations [76]. The following equations were considered:

$$Nu^{opt\_W}(Pr, Re) = a_1 + (b_1 Re^{c_1} + d_1 Re^{e_1}) Pr^{f_1} \quad (24)$$

$$Nu^{opt\_FM}(Pe, Re) = a_2 + b_2 Pe^{c_2} + d_2 Pe^{e_2} Re^{f_2} \quad (25)$$

where {$a_i, f_i$} are constant coefficients determined with the least square method (see Table 7). Both new correlations are represented on Fig. 8. It should be emphasized that the plots are superimposed, even if the equation is not in the same form. It suggests that there are too many degrees of freedom in determining the correlation coefficients. In addition, these correlations fail to predict the Nusselt value when the Reynolds number tends to zero. Indeed, for a sphere in a stationary fluid, the Nusselt number equals 2 [78]. A possible explanation is the fact that the computational domain is too small in front of the particle
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Assessment of the correlation of Whitaker and Feng and Michaelides with the data produced with PRS. Opt FM and Opt W refers to the correlations obtained with optimized coefficients of Feng and Michaelides, and Whitaker correlations, respectively.</figcaption></figure>

<caption>**Table 7**</caption>

Coefficients obtained with the least square method for Eqs. (24) and (25). W: coefficients of Whitaker's correlation [75], FM: coefficients of Feng & Michaelides's correlation [76].

<table><thead><tr><th></th><th>a<sub>i</sub></th><th>b<sub>i</sub></th><th>c<sub>i</sub></th><th>d<sub>i</sub></th><th>e<sub>i</sub></th><th>f<sub>i</sub></th></tr></thead><tbody><tr><td>i = 1</td><td>1.5</td><td>0.072</td><td>-2.5</td><td>0.73</td><td>0.44</td><td>0.37</td></tr><tr><td>i = 2</td><td>1.5</td><td>0.64</td><td>0.38</td><td>0.10</td><td>0.34</td><td>0.30</td></tr><tr><td>W</td><td>2</td><td>0.4</td><td>0.5</td><td>0.06</td><td>2/3</td><td>2/5</td></tr><tr><td>FM</td><td>0.992</td><td>1</td><td>1/3</td><td>0.1</td><td>1/3</td><td>1/3</td></tr></tbody></table>

diameter. At low Peclet numbers, the thermal diffusion is predominant in front of the thermal convection and the thermal boundary layer is thicker than the dynamic one for low Prandtl number. Thus, the influence of the domain boundaries is higher at low Peclet numbers and could explain the discrepancy.

## 5. Study of an anisothermal liquid-solid fluidized bed

In this section, an anisothermal liquid-solid fluidized bed is studied. The physical configuration is described in Section 5.1. The numerical method is validated by confrontation with experimental [79] and numerical [80] work of the literature in Section 5.2. Fluid-particle heat transfer and wall-to-bed heat transfer are addressed in Sections 5.3 and 5.4 respectively. Finally, velocity-temperature correlations are discussed in Section 5.5.

### 5.1. Description of the physical configuration

The simulation reproduces the experimental work of Corona [79, 81]. The author investigated an isothermal liquid-solid fluidized bed in a cylindrical vessel. The working fluid is a concentrated aqueous solution of Potassium Thiocyanate (KSCN, 64% w/w). The particles are 6 mm Pyrex spheres. The fluid viscosity is $3.8 \times 10^{-3}$ Pa s. A viscosity ratio of $10^4$ is imposed as in [65]. The density of the fluid and the particles is, respectively, 1400 kg m$^{-3}$ and 2230 kg m$^{-3}$. In the experiment, the cylindrical vessel is made of glass. The column is 8 cm in diameter and 60 cm height. The exponent of the Richardson-Zaki fluidization law [82] was measured experimentally by the authors as equal to 2.41. The terminal velocity of the particles in the fluid is $0.23 \text{ m s}^{-1}$ and the minimum fluidization velocity is estimated to $0.03 \text{ m s}^{-1}$ based on Fig. 2.24 of [79].

To numerically reproduce the experiment, a parallelipedic domain is defined with the same cross section area (see Fig. 9). The square cross-section is slightly larger than the circular cross-section by 3.1%. A no-slip boundary condition is imposed on the lateral walls. The fluid is injected at the bottom of the tank with a velocity $U_f$. A pressure imposed boundary condition is applied to the outlet section. The wall temperature is imposed 10 K higher than that of the particles. In the present work, all particles remain at a constant temperature. The entrance temperature of the fluid is that of the particles so only wall-to-bed heat transfer is considered. The thermal conductivity and the heat capacity of the fluid are constant in the whole domain ($\lambda_f = 3.8 \text{ W m}^{-1} \text{ K}^{-1}$, $c_{p,f} = 1000 \text{ J kg}^{-1} \text{ K}^{-1}$). The thermal properties of the fluid were selected so that the Prandtl number equals 1. Four fluidization velocities were studied, $U_f \in \{0.073, 0.090, 0.120, 1.150\} \text{ m s}^{-1}$. Three mesh resolutions were considered for the lowest fluidization regime ($\frac{D_p}{\Delta x} \in \{12, 24, 36\}$), whereas two mesh resolutions were investigated for the fluidization velocities of $0.090 \text{ m s}^{-1}$ and $0.120 \text{ m s}^{-1}$ ($\frac{D_p}{\Delta x} \in \{12, 24\}$). Finally, one mesh resolution is considered for the highest fluidization velocity ($\frac{D_p}{\Delta x} = 12$). The numerical configuration of each case is detailed in Table 8. The height of the computation domain was shortened for the lowest fluidization velocities in order to reduce the numerical cost of the simulation. Preliminary tests showed that it did not affect the bed height, the solid volume fraction or the particle velocity variances. The computational time step depends on the grid resolution and the fluidization velocity. For all the simulations, $\Delta t_{simu}$ lies between $5 \cdot 10^{-5}$ s and $1 \cdot 10^{-4}$ s.

### 5.2. Hydrodynamic behavior

The numerical method was previously validated by Hamidi et al. [65] for the same configuration of the fluidized bed with a mesh resolution of $\frac{D_p}{\Delta x} = 12$. The authors studied the bed height, the mean solid fraction as well as velocity variances and its anisotropy. The method has evolved since the TrioCFD-1.7.9 version used in [65]. A harmonic viscosity model is employed in the present work compared to the on-off model employed in [65]. The on-off viscosity model implicitly define an effective diameter smaller than the real diameter for the computation of the diffusive term in the momentum equation. For coarse meshes, the consideration of the real diameter lead to an overestimation of the drag force [64]. However, the harmonic viscosity model should be employed to compute the stress tensor as it was proved
<table><caption>**Table 8**
Numerical configuration of fluidized beds.</caption><thead><tr><th>Studied case</th><th>Dimensions (L/D<sub>p</sub>)</th><th>Resolution (D<sub>p</sub>/Δx)</th><th>Nb Eulerian cells</th><th>Nb CPU's</th><th>Physical simulated time (s/24 h computation)</th></tr></thead><tbody><tr><td rowspan="2">v &lt; 0.15 m s<sup>-1</sup></td><td rowspan="2">12 × 54 × 12</td><td>12</td><td>13.4 · 10<sup>6</sup></td><td>288</td><td>1.74</td></tr><tr><td>24</td><td>107 · 10<sup>6</sup></td><td>1152</td><td>0.8</td></tr><tr><td>v = 0.073 m s<sup>-1</sup></td><td>12 × 27 × 12</td><td>36</td><td>181 · 10<sup>6</sup></td><td>3888</td><td>0.35</td></tr><tr><td>v = 0.15 m s<sup>-1</sup></td><td>12 × 108 × 12</td><td>12</td><td>26.9 · 10<sup>6</sup></td><td>544</td><td>1.45</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Snapshot of a simulation at the highest solid fraction ($\alpha$ = 0.36, $U_f$ = 0.073 m s⁻¹).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Time evolution of the bed height for different fluidization velocities and mesh resolutions. ■: present work, ■: numerical work of Ozel [80]. —: $U_f$ = 0.15 m s⁻¹, - $\frac{D_x}{\Delta x}$ = 12, —: $U_f$ = 0.12 m s⁻¹, - $\frac{D_x}{\Delta x}$ = 12, —: $U_f$ = 0.09 m s⁻¹, - $\frac{D_x}{\Delta x}$ = 12, —: $U_f$ = 0.073 m s⁻¹, - $\frac{D_x}{\Delta x}$ = 12, —: $U_f$ = 0.073 m s⁻¹, - $\frac{D_x}{\Delta x}$ = 24, —: $U_f$ = 0.073 m s⁻¹, - $\frac{D_x}{\Delta x}$ = 36. For $U_f$ = 0.073 m s⁻¹, all mesh resolutions are superimposed.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Fluidization law. ■: experiment of Corona [79], —: Richardson-Zaki law [82], x: numerical work of Ozel [80], +: numerical work of Hamidi [65], ♦: present work for $\frac{D_x}{\Delta x}$ = 12.</figcaption></figure>

to converge with the grid resolution [67]. For the simulations performed in this paper, with Reynolds numbers above 460 for meshes of up to 36 meshes per diameter, mesh convergence is not achieved. This leads to a slight overestimation of the drag force, which consequently results in an underestimation of the presence rate compared to Hamidi et al. [65] and Corona [79] experimental data. For this reasons, the numerical method is further evaluated against the experimental work of Corona [79] and the numerical work of Ozel [80]. Ozel et al. employed the Implicit Tensorial Penalty Fictitious Domain Method (ITPM) of Vincent et al. [64] for their simulations. A parallelipedic domain was considered and a Darcy penalty method was employed to add a wall boundary condition in the momentum equation and thus consider a cylindrical column. A grid resolution of 12 meshes per diameter was defined by the authors. The temporal evolution of the bed height is represented in Fig. 10. The bed height closely matches the simulation conducted by Ozel for the lowest fluidization velocities with the same grid resolution. The higher the fluidization velocity, the larger the fluctuations in bed height. The mesh sensitivity is represented in the figure for the lowest fluidization velocity. As shown in the figure, the bed height is slightly affected by the grid resolution. The discrepancy observed for the highest fluidization velocities could be explained by the fact that the cross-section area differs of 3.1% between both computational domain. It could also be explained by the difference in the geometry of the cross-section as a square is defined in the present work whereas a disk is defined in the work of Ozel [80]. The average solid fraction is computed between 5% and 85% of the bed height to avoid any boundary effects. It is represented in Fig. 11 as a function of the fluidization velocity. The fluidization law is well predicted and closely matches the experimental work of Corona [79] and the numerical work of Ozel [80] and Hamidi [65]. The slight deviation observed is explained by the choice of the harmonic viscosity model as detailed above. A reduced diameter of $\Delta x$ was defined in [64] and employed in [80]. The viscosity model defined in [65] mimics
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Instantaneous visualization of the streamwise velocity field. $U_y = 0.073 \ \text{m s}^{-1}$. Left: $\frac{D_x}{\Delta x} = 12$, Middle: $\frac{D_x}{\Delta x} = 24$, Right: $\frac{D_x}{\Delta x} = 36$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Mean temperature profile along the bed height for different fluidization velocities $U_f$ and mesh resolutions $\frac{D_x}{\Delta x}$. —: $U_f = 0.073 \ \text{m s}^{-1}$, =: $U_f = 0.09 \ \text{m s}^{-1}$, U: $U_f = 0.12 \ \text{m s}^{-1}$, ■: $U_f = 0.15 \ \text{m s}^{-1}$. —: $\frac{D_x}{\Delta x} = 12$, ■: $\frac{D_x}{\Delta x} = 24$, ■: $\frac{D_x}{\Delta x} = 36$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Profile of the double temperature correlation along the bed height for different fluidization velocities $U_f$ and mesh resolutions $\frac{D_x}{\Delta x}$. —: $U_f = 0.073 \ \text{m s}^{-1}$, =: $U_f = 0.09 \ \text{m s}^{-1}$, ■: $U_f = 0.12 \ \text{m s}^{-1}$, =: $U_f = 0.15 \ \text{m s}^{-1}$. —: $\frac{D_x}{\Delta x} = 12$, ■: $\frac{D_x}{\Delta x} = 24$, ■: $\frac{D_x}{\Delta x} = 36$. For all mesh resolution, curves are superimposed along the bed height.</figcaption></figure>

the reduced diameter of the penalty method. In the present work, the harmonic viscosity model was employed as it is expected to converge with the grid resolution.

Instantaneous visualization of the streamwise velocity field is shown in Fig. 12. The whole domain is represented for the grid resolutions of 24 and 36 grid per particle diameter while the domain is cropped for the coarsest grid to confront the above grid resolutions. The particle wake appears to be more accurately captured for the 24-mesh resolution than for the 12-mesh resolution. On the other hand, few differences are visible between the resolutions of 24 and 36 meshes per diameter.

### 5.3. Fluid-particle heat transfer

The mean temperature profiles along the height of the column are represented in Fig. 13 for the four fluidization velocities. For details on the way this quantity is computed from the simulations results, the reader is referred to Appendix A. For hereinafter, the notation $\langle A \rangle_{i,j}$ denotes for the spatial average of the quantity $A$ in directions $i, j$ and $\bar{A}$ denotes for the time average. In the figure, one can observe that the lower the fluidization velocity, the more homogeneous the fluid

temperature in the bed. Indeed, the particles act as a heat sink in the simulations and the lower the fluidization velocity, the denser the fluidized bed and therefore the more homogeneous the temperature field. Particles thus limit the increase in fluid temperature. One can point out that the bed height can be identified by the change in the slope of the temperature profile along the height of the column. Above the later, the fluid is simply heated by the walls and the temperature variation with $y$ significantly rises. One can point out that a grid resolution of 12 meshes per particle diameter is too coarse to accurately capture the temperature evolution along the column height. A large discrepancy is observed between the evolution of the fluid temperature for grid resolutions of 12 and 24 meshes per diameter. However, the evolution of the fluid temperature for grid resolutions of 24 and 36 meshes per particle diameter, closely match. Therefore, for this configuration, at least 24 meshes per particle diameter is required to predict the mean temperature profiles.

The double correlation of temperature temporal fluctuations are represented with respect to the bed height and the wall distance in Figs. 14 and 15, respectively. The computation of $\langle (T'T')^f \rangle_x$ and
<figure><img src="image_17.png" /><figcaption>Fig. 15. Profile of the double temperature correlation as a function of $x^+$ for different fluidization velocities $U_f$ and mesh resolutions $\frac{D_p}{\Delta x}$. $\cdot$: $U_f = 0.073 \text{ m s}^{-1}$, $\cdot\cdot$: $U_f = 0.09 \text{ m s}^{-1}$, $\cdot\cdot\cdot$: $U_f = 0.12 \text{ m s}^{-1}$, $\cdot\cdot\cdot\cdot$: $U_f = 0.15 \text{ m s}^{-1}$. $\cdot\cdot\cdot\cdot\cdot$: $\frac{D_p}{\Delta x} = 12$, $\cdot\cdot\cdot\cdot\cdot\cdot$: $\frac{D_p}{\Delta x} = 24$, $\cdot\cdot\cdot\cdot\cdot\cdot\cdot$: $\frac{D_p}{\Delta x} = 36$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Relation between $x^+$ and the wall distance normalized by the particle diameter. The legend is the same as for Fig. 15.</figcaption></figure>

$\langle(T'T')_f\rangle_v$ is detailed in Appendix A. One can notice in Fig. 14 that the higher the Reynolds number, the higher the temperature variance. The influence of the lack of resolution for coarse grid has less effects than for the mean temperature profiles (see Fig. 13). In addition, when the fluid leaves the bed, temperature variance increase strongly. This can be explained by the destruction of fluid agitation by pseudo-turbulence. Once exiting the bed, turbulence develops, generating strong spatial and temporal fluctuations in the temperature field. The effect of the lack of resolution of near-wall flow is more visible in Fig. 15. In the figure, $x^+$ is a dimensionless distance to the wall. $x^+ = \frac{|x-x_{wall}|u^+}{v_f}$, $u^+ = \frac{\tau_w}{\rho_f}$ and $\tau_w = \mu_f \frac{\partial u_y}{\partial x}$. The link between $x^+$ and the wall distance scaled by the particle diameter is described in Fig. 16. In Fig. 15, the horizontal lines for low values of $x^+$ are due to the under-resolved discretization. All values of the double correlation of the temperature are 0 at the wall. Contrary to what is observed in Fig. 14, the maximum values of the correlation do not appear to depend on the Reynolds number. The correlation peak is located approximately at $x^+ = 10$, which correspond to a distance below 10% of the particle diameter. It emphasizes the necessity to refine the mesh at the wall. It is confirmed by the significant influence of the mesh resolution on the location of the correlation peak. The temperature variance falls to zero in the center of the bed, as the fluid and particles are at the same temperature.

The mean Nusselt number is defined as follows:

$$ \langle \overline{Nu_f} \rangle_{0D} = \frac{\langle \overline{h_f} \rangle_{0D} D_p}{\lambda_f} \qquad (26) $$

<figure><img src="image_19.png" /><figcaption>Fig. 17. Nusselt number as a function of the solid fraction. $\cdot$: $\frac{D_p}{\Delta x} = 12$, $\cdot\cdot$: $\frac{D_p}{\Delta x} = 24$, $\cdot\cdot\cdot$: $\frac{D_p}{\Delta x} = 36$. $\cdot\cdot\cdot\cdot$: Gunn's correlation [44].</figcaption></figure>

where $h_f$ is the heat transfer coefficient, $D_p$ is the particle diameter and $\lambda_f$ is the thermal conductivity of the fluid. $h_f$ writes:

$$ \langle \overline{h_f} \rangle_{0D} = \frac{\Phi_{f \to p}}{N_p S_p (T_p - \langle \overline{T_f} \rangle_{0D})} \qquad (27) $$

where $\Phi_{f \to p}$ is the total heat flux received by the particles. $\Phi_{f \to p}$ is computed with a heat balance at the boundaries as described in Section 4.1. $N_p$ is the number of particles in the bed, $T_p$ is the temperature of the particles and $\langle \overline{T_f} \rangle_{0D}$ is the mean fluid temperature in the bed. Therefore, one can express the Nusselt number as follows:

$$ \langle \overline{Nu_f} \rangle_{0D} = \frac{\Phi_{f \to p}}{N_p \lambda_f \pi D_p (T_p - \langle \overline{T_f} \rangle_{0D})} \qquad (28) $$

The Nusselt computed for each fluidization velocity and each mesh grid resolution is represented in Fig. 17. The computed Nusselt number increases with the solid fraction. A large discrepancy is observed with the correlation of Gunn [44] even for the mesh grid resolution of 24 and 36 meshes per particle diameter. In such configuration, the fluid is injected at the same temperature as the particles and is only heated by the lateral walls. Thus, near-wall particles are submitted to a mean temperature gradient. Gunn's correlation [44] was not developed to describe such fluid-to-particle heat transfer. To assess Gunn's correlation [44] in the present case, the mean heat flux received by the particles ($\Phi_{f \to p}$) should only be computed in the homogeneous part of the fluidized bed. However, the fluid-to-particle heat flux was not post-processed on a Lagrangian way for each particle in the current version of the code. Only the mean particle-to-fluid heat flux is available for post-processing. Nevertheless, it should be noted that Gunn's correlation always underestimates the fluid-particle heat transfer, which is expected as it was originally developed for fixed bed configurations. The same results were observed in [51]. Thus, as shown in Fig. 17, Gunn's correlation is not suitable to predict the overall convective heat transfer in the present configuration. Furthermore, even with 36 meshes per particle diameter, the mesh grid resolution is too coarse for the studied regime. Indeed, the particle Reynolds number defined as $Re = \frac{\rho_f U_f D_p}{\alpha_p \mu_f}$ is 461 for the lowest fluidization velocity and $Pr = 1$ and thus the Peclet number equals the Reynolds number. Consequently, a complete resolution of the thermal boundary layer would require more than a hundred meshes per particle diameter.

## 5.4. Wall-to-bed heat transfer

The wall-to-bed heat transfer is computed for the four walls of the parallelepipedic domain with a first order upwind scheme. The computation of the heat transfer coefficient is detailed in Appendix B. Instantaneous visualization of the heat transfer coefficient is shown in Fig. 18 for various mesh resolutions. The thermal wake is more
<figure><img src="image_20.png" /><figcaption>Fig. 18. Instantaneous wall-to-bed heat transfer coefficient. From left to right, the grid resolution is 12, 24 and 36 meshes per particle diameter.</figcaption></figure>

spread out for the coarser grid than for the others and the magnitude of the heat transfer is lower. A large discrepancy is observed between the grid resolutions of 12 and 24 meshes per particle diameter. Even though the thermal wake of the particles is thinner for the finest grid, few differences are observed between the two finest grid resolutions, suggesting mesh convergence.

The mean heat transfer coefficient is represented in Fig. 19. The wall-to-bed heat transfer decreases with the fluid fraction. This behavior is in agreement with the correlation of Haid [53]. The author defined the Nusselt number as a function of the Reynolds number, the Prandtl number and the solid fraction (see Eq. (29)).

$$Nu_p = 0.0734 Re_p^{0.75} Pr^{0.63} \alpha_p^{0.25} (1 - \alpha_p)^{-1} \quad (29)$$

To assess the correlation, the Reynolds number is computed with the velocity corresponding to a given solid fraction according to the Richardson–Zaki law [82]. One can point out an important effect of the mesh resolution on the wall-to-bed heat transfer coefficient. For the lowest fluid fraction, corresponding to the lowest fluidization velocity, the computed coefficient start to falls within the 30% confidence interval for the finest mesh. The finer the grid resolution, the less intense the heat transfer, contrary to what Fig. 18 suggests. Indeed, even if the points of contact with the wall are subject to high instantaneous values of the heat transfer coefficient, the contact zone and the collision time are very small. The higher coefficient computed for the coarser grid resolutions is therefore due to the more extensive thermal wake than for the finer resolutions. This results from the lack of resolution of pseudoturbulence dissipation. An unexpected result is that the error is smaller for the largest fluid fractions and therefore, for the highest fluidization velocities. Indeed, for the fluidization velocity of $U_f = 0.12 \text{ m s}^{-1}$ ($1 - \alpha = 0.8$), the heat transfer coefficient falls within the confidence interval for a grid resolution of 24 meshes per particle diameter. It should be noted, however, that the current configuration does not fall within the range of validity of the correlation. Indeed the correlation is valid for a Prandtl number between 1.65 and 7700, whereas in the present simulations $Pr = 1$.

## 5.5. Velocity–temperature correlation

Velocity–temperature correlations are investigated to study heat transfers. The mean and bulk temperatures are computed along the cross-section of the flow. The bulk temperature is defined in Appendix A. The mean and bulk temperature profiles along the cross-section area are represented in Fig. 21. The higher the velocity, the hotter the mean temperature of the fluid along the cross-section area. Indeed, the

<figure><img src="image_21.png" /><figcaption>Fig. 19. Space–time averaged wall-to-bed heat transfer coefficient. $x$: $\frac{D_x}{\Delta x} = 12$, $x$: $\frac{D_x}{\Delta x} = 24$, $x$: $\frac{D_x}{\Delta x} = 36$. —: Haid correlation [53] and its confidence interval of 34.4%.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 20. Fluidization velocity effect on the near wall evolution of space–time average of the velocity–temperature correlation. $\frac{D_x}{\Delta x} = 12$. —: $U_f = 0.073 \text{ m s}^{-1}$, —: $U_f = 0.09 \text{ m s}^{-1}$, —: $U_f = 0.12 \text{ m s}^{-1}$, —: $U_f = 0.15 \text{ m s}^{-1}$. —: streamwise velocity, —: normal wall velocity, —: wall tangent velocity.</figcaption></figure>

higher the velocity, the less cooled the fluid by the particles along the bed height (see Fig. 13). Strong spatial variations of the temperature are observed along the transverse direction to the flow as the fluid
<figure><img src="image_23.png" /><figcaption>**Fig. 21.** Mean and bulk fluid temperature along the cross-section area. $\frac{D_x}{\Delta x} = 12$. —: mean temperature, =: bulk temperature, $U_f = 0.073 \text{ m s}^{-1}$, =: $U_f = 0.09 \text{ m s}^{-1}$, =: $U_f = 0.12 \text{ m s}^{-1}$, =: $U_f = 0.15 \text{ m s}^{-1}$.</figcaption></figure>

is only heated at the wall. The computation of velocity–temperature correlation $\langle(U'T')_f\rangle_y$ is detailed in Appendix A.

**Fig. 20** represents this correlation, for the four fluidization velocities, along the cross-section area. The fluid Reynolds number can be defined as $Re = \frac{U_f L_x}{\nu_f \alpha_f}$. From the lowest to the highest fluidization velocity, the Reynolds number is respectively 2980, 3360, 3980 and 4970. The extremum of the double correlation is located further into the bed when the Reynolds number increases. The main component of the turbulent movement is obtained with the streamwise velocity. One can note that the turbulent agitation associated with the normal wall velocity is less than a third of the intensity of that of the streamwise velocity. The correlation peak associated with the normal wall velocity is located further into the bed than that of the streamwise component. Furthermore, the location of the correlation peak differs from that of the temperature variance (see Fig. 15) which was located at $x^+ \sim 10$. For both studied correlation, it corresponds to a distance less than half of the particle diameter. Finally, the velocity tangent to the wall does not contribute to the turbulent heat flow. Results on the average of each velocity component showed that the mean normal wall velocity is almost zero. Therefore, the mean convective heat flux, computed as $\rho_f c_{p,f} \langle u_x \rangle_{0D} \langle \overline{u} \overline{u} \rangle_{0D}$, is very small in front of the term $\rho_f c_{p,f} \langle u_x' T' \rangle_{0D}$. Wall-to-bed heat transfer is therefore driven by the turbulent convective heat flux. One should note that it is not the case for the streamwise mean convective heat flux as the mean fluid velocity is equal to the fluidization velocity divided by the mean fluid fraction. However, for heat transfer enhancement in the bed, it is the convective heat flux associated to the normal wall velocity that needs to be considered. Mesh resolution has a strong influence on the location of the agitation peak, as shown in Fig. 22. Indeed, the agitation peak associated with the streamwise component is located at $x^+ = 16$ for the mesh resolution of 12 meshes per particle diameter whereas it is located at $x^+ = 27$ for the mesh resolution of 36 meshes per particle diameter. The intensity of the turbulent heat flow is slightly affected by the grid resolution.

# 6. Concluding remarks

The objective of this study was to study the dynamic and the heat transfers in a fluidized bed. For this purpose, the numerical method was first validated on the study of an isolate particle in sedimentation in an infinite medium. A method to compute the heat flux received by the particle was developed, inspired by the method of hydrodynamic force computation of Butaye et al. [19]. Then a parametric study of the particle Nusselt number of a particle in sedimentation in quiescent fluid was performed. An excellent agreement with existing correlation of the literature was found except for very low Peclet numbers (below

<figure><img src="image_24.png" /><figcaption>**Fig. 22.** Resolution effect on the near wall evolution of space-time average of the velocity–temperature correlation. $U_f = 0.073 \text{ m s}^{-1}$, $\frac{D_x}{\Delta x} = 12$, $\frac{D_x}{\Delta x} = 24$, $\frac{D_x}{\Delta x} = 36$. The components of the velocity are described with the same representation than that of Fig. 20.</figcaption></figure>

1) which could be explained by the influence of the lateral boundaries at these regimes.

A study of a liquid-solid anisothermal fluidized bed showed that the macroscopic behavior of the bed was well predicted based on the bed height and the solid fraction. However, for coarse grid ($D_p/\Delta x = 12$), the method fails to predict the thermal behavior of the flow. A finer mesh ($D_p/\Delta x = 24$) was required to predict the heat transfer at a reasonable cost. The fluid Nusselt number increases with the solid fraction, in accordance with Gunn's correlation [44]. Finally, regarding the wall-to-bed heat transfer, the heat transfer coefficient increases with the solid fraction. The deviation from Haid's correlation [53] increases with the solid fraction, which is unexpected since, for the same grid resolution, the resolution of the flow decreases as the velocity increases. The velocity-temperature correlations are strongly affected by the grid resolution. Finally, it was shown that the wall-to-bed heat transfer is driven by the turbulent convective heat flux associated with the normal wall velocity.

The perspectives to this work would be to study anisothermal gas-solid fluidized beds. In such configurations, the particle diameter is less than a hundred micrometers and the fluid studied is air.

## CRediT authorship contribution statement

**E. Butaye:** Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.

**R. Quintana:** Formal analysis, Data curation.

**S. Mer:** Writing – review & editing, Visualization, Validation, Supervision, Methodology, Investigation, Formal analysis, Conceptualization.

**F. Baille:** Writing – review & editing, Validation, Supervision.

**A. Toutant:** Writing – review & editing, Visualization, Validation, Supervision, Methodology, Investigation, Formal analysis, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgment

This work was granted access to the HPC resources of CINES under the allocation 2023-A0142B11441 made by GENCI. The authors gratefully acknowledge the CEA for the development of the TRUST platform. The technical support of TRUST/TrioCFD team was also greatly appreciated. The authors would like to thank Dr. Ali Ozel for providing the simulation data.
<figure><img src="image_25.png" /><figcaption>Fig. A.23. Diagram of planes employed in the spatial average. Along streamwise direction (left) and along the cross-section area (right).</figcaption></figure>

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

## Appendix A. Computation of Eulerian statistics

The mean Eulerian properties of the flow are computed with a space-time average. Two configurations can be distinguished. The first is an average over planes parallel to the flow and the second is an average over the cross-section area (see Fig. A.23). For the first configuration, all the planes at the centre of the domain are considered over a thickness of two particle diameters and along the two transverse directions. Thus, for resolutions of 12, 24 and 36 meshes per diameter, we consider 48, 96 and 144 planes respectively. Three cross-sections are considered for the second configuration, at a height of 25%, 50% and 75% of the mean bed height.

For a field A, the space-time average in a two-dimensional plane $\overline{\langle A \rangle}$ is computed as follows:

$$ \overline{\langle A \rangle}(x, y) = \frac{1}{N_t} \sum_i \left( \frac{1}{N_{\text{planes}}^i} \sum_{\text{planes}} A(x, y, z, t) \right) $$

where $N_t$ is the number of timesteps considered in the simulation and $N_{\text{planes}}$ is the number of planes.

The zero-dimensional average of the quantity $\overline{\langle A \rangle}(x, y)$ is then computed on a section of the two-dimensional plane. To avoid any boundary effects, cells located at a distance of $3D_p$ near the fluid injection are not taken into account, nor are cells located at a distance of one particle diameter near the corners. In addition, as only the characteristics of the fluid inside the bed are studied, cells located above the average height of the bed are not taken into account in the average. Then, the 0D average writes:

$$ \overline{\langle A \rangle}_{0D} = \int_{3D_p}^{h_{\text{bed}}} \int_{-L_x/2+D_p}^{L_x/2-D_p} \overline{\langle A \rangle}(x, y) dx dy $$

One dimensional space-time average in the direction x or in the direction y are computed as follows:

$$ \overline{\langle A \rangle}_y = \int_{3D_p}^{h_{\text{bed}}} \overline{\langle A \rangle}(x, y) dy \\ \overline{\langle A \rangle}_x = \int_{-L_x/2+D_p}^{L_x/2-D_p} \overline{\langle A \rangle}(x, y) dx $$

The mean fluid temperature is computed over the planes parallel to the streamwise direction of the flow and writes:

$$ \overline{\langle T \rangle}_{f,0D} = \frac{\overline{\langle \alpha_f T \rangle}_{0D}}{\overline{\langle \alpha_f \rangle}_{0D}} $$

<figure><img src="image_27.png" /><figcaption>Fig. C.24. Diagram of the Verlet algorithm (not to scale).</figcaption></figure>

The bulk temperature writes:

$$ \overline{\langle T_b \rangle}_{f,0D} = \frac{\overline{\langle \alpha_f V T \rangle}_{0D}}{\overline{\langle \alpha_f V \rangle}_{0D}} $$

where V stands for the streamwise velocity.

Double correlation for field A and field B are computed as follows:

$$ \overline{\langle (A'B')_f \rangle}_y = \frac{\overline{\langle \alpha_f AB \rangle}_y}{\overline{\langle \alpha_f \rangle}_y} - \frac{\overline{\langle \alpha_f A \rangle}_y \overline{\langle \alpha_f B \rangle}_y}{\overline{\langle \alpha_f \rangle}_y^2} $$

$$ \overline{\langle (A'B')_f \rangle}_x = \frac{\overline{\langle \alpha_f AB \rangle}_x}{\overline{\langle \alpha_f \rangle}_x} - \frac{\overline{\langle \alpha_f A \rangle}_x \overline{\langle \alpha_f B \rangle}_x}{\overline{\langle \alpha_f \rangle}_x^2} $$

## Appendix B. Computation of the wall-to-bed heat transfer

The mean heat transfer coefficient is computed as follows:

$$ \bar{h} = \frac{1}{N_t} \sum_{i} \frac{1}{4} \sum_{i} \frac{1}{S_{i}^{l}} \iint_{S_{i}^{l}} h_{i}^{l}(x, y, z) dS \quad (B.1) $$

where $N_t$ is the number of timesteps considered for the time-average, $S_{i}^{l}$ is the surface of the lateral wall i and depends on the bed height,
$h_i^t(x, y, z)$ is the heat transfer coefficient at the position $(x,y,z)$ of the lateral wall $i$ at the time $t$.

## Appendix C. Verlet algorithm

1. **Construction of the Verlet tables** At the start of the calculation, for each particle $i$, the distance $d_{ij}$ with all the other particles is computed. If $d_{ij}$ is below a given threshold $\epsilon_v$, the particle $j$ is identified as a potential candidate for collision for particle $i$ and saved in its Verlet table. The algorithm is schematized in Fig. C.24

2. **Identifying the fastest particle** The speed of the fastest particle is identified as $v_{p,max}$.

3. **Computation of the update time** A characteristic time of the table update is constructed based on the time two particles at $v_{p,max}$, distant by $\epsilon_v$ and traveling in the opposite direction, need to collide $\Delta t_V = \frac{\epsilon_v}{2v_{p,max}}$. A numerical security $\Delta t^{max}$ is also considered to avoid excessively long time intervals without update. The final update time is $\Delta t_{upd} = \min(\Delta t_V, \Delta t_{max})$.

4. **Update for the tables** Every $\Delta t_{upd}$ the Verlet tables are recalculated

## Data availability

Data will be made available on request.

## References

[1] Y. Zhang, J. Xu, Q. Chang, P. Zhao, J. Wang, W. Ge, Numerical simulation of fluidization: Driven by challenges, Powder Technol. 414 (2023) http://dx.doi.org/10.1016/j.powtec.2022.118092.

[2] F. Winkler, Verfahren zum Herstellen von Wassergas, (DE437970C) 1922, https://patents.google.com/patent/DE437970C/de.

[3] T. Suppan, M. Neumayer, T. Bretterklieber, S. Puttinger, H. Wegleiter, A. model-based analysis of capacitive flow metering for pneumatic conveying systems: A comparison between calibration-based and tomographic approaches, Sensors 22 (2022) http://dx.doi.org/10.3390/s22030856.

[4] N. Muradov, Hydrogen via methane decomposition: an application for decarbonization of fossil fuels, Int. J. Hydrog. Energy 26 (2001) http://dx.doi.org/10.1016/S0360-3199(01)00073-8.

[5] G. Flamant, D. Gauthier, H. Benoit, J.-L. Sans, R. Garcia, B. Boissiere, R. Ansart, M. Hemati, Dense suspension of solid particles as a new heat transfer fluid for concentrated solar thermal plants: On-sun proof of concept, Chem. Eng. Sci. 102 (2013) http://dx.doi.org/10.1016/j.ces.2013.08.051.

[6] M. Castro-Quijada, D. Faundez, R. Rojas, A. Videla, Improving the working fluid based on a $NaNO_3-KNO_3-NACl-KCl$ molten salt mixture for concentrating solar poser energy storage, Sol. Energy 231 (2022) http://dx.doi.org/10.1016/j.solener.2021.11.058.

[7] R. Gueguen, G. Sahuquet, M. Tessoneaud, J.-L. Sans, E. Guillot, A.L. Gal, R. Garcia, S. Mer, A. Toutant, F. Bataille, G. Flamant, Heat transfer in a fluidized bed tubular solar receiver. On-sun experimental investigation, Sol. Energy 265 (2023) http://dx.doi.org/10.1016/j.solener.2023.112118.

[8] G. Flamant, Theoretical and experimental study of radiant heat transfer in a solar fluidized-bed receiver, AIChE J. 28 (1982) http://dx.doi.org/10.1002/aic.690280402.

[9] R. Gueguen, G. Sahuquet, S. Mer, A. Toutant, F. Bataille, G. Flamant, Gas-solid flow in a fluidized-particle tubular solar receiver: Off-sun experimental flow regimes characterization, Energies 14 (2021) http://dx.doi.org/10.3390/en14217392.

[10] B. Boissiere, R. Ansart, D. Gauthier, G. Flamant, M. Hemati, Experimental study of gas-particle dense suspension upward flow for application as a new heat transfer and storage fluid, Can. J. Chem. Eng. 93 (2014) http://dx.doi.org/10.1002/cjce.22087.

[11] D. Geldart, Types of gas fluidization, Powder Technol. 7 (1973) http://dx.doi.org/10.1016/0032-5910(73)80037-3.

[12] R. Gueguen, S. Mer, A. Toutant, F. Bataille, G. Flamant, Effect of temperature on the hydrodynamics of a fluidized bed circulating in a long tube for a solar energy harvesting application, Chem. Eng. Sci. 281 (2023) http://dx.doi.org/10.1016/j.ces.2023.119218.

[13] R. Gueguen, G. Sahuquet, S. Mer, A. Toutant, F. Bataille, G. Flamant, Fluidization regimes of dense suspensions of Geldart group A fluidized particles in a high aspect ratio column, Chem. Eng. Sci. 267 (2022) http://dx.doi.org/10.1016/j.ces.2022.118360.

[14] M.V. der Hoef, M. Ye, M. van Sint Annaland, A.A. IV, S. Sundaresan, J. Kuipers, Multiscale modeling of gas-fluidized beds, Adv. Chem. Eng. 31 (2006) http://dx.doi.org/10.1016/S0065-2377(06)31002-2.

[15] F. Alobaid, N. Almohammed, M.M. Farid, J. May, P. Rößger, A. Richter, B. Apple, Progress in CFD simulations of fluidized beds for chemical and energy process engineering, Prog. Energy Combust. Sci. 91 (2022) http://dx.doi.org/10.1016/j.pecs.2021.100930.

[16] S. Tenneti, S. Subramaniam, Particle-resolved direct numerical simulation for gas-solid flow model development, Annu. Rev. Fluid. Mech. 46 (2014).

[17] M. Maxey, Simulation methods for particulate flows and concentrated suspension, Annu. Rev. Fluid. Mech. 49 (2017) http://dx.doi.org/10.1146/annurev-fluid-122414-034408.

[18] F. Marchelli, L. Fiori, R. Di Felice, Cohesive particle-fluid systems: An overview of their CFD simulation, Can. J. Chem. Eng. (2024) http://dx.doi.org/10.1002/cjce.25269.

[19] E. Butaye, A. Toutant, S. Mer, F. Bataille, Development of particle resolved - subgrid corrected simulations: Hydrodynamic force calculation and flow sub-resolution corrections, Comput. & Fluids 267 (2023) http://dx.doi.org/10.1016/j.compfluid.2023.106071.

[20] H.A. Khawaja, CFD-DEM simulation of minimum fluidisation velocity in two phase medium, Int. J. Multiphys. 5 (2011) http://dx.doi.org/10.1260/1750-9548.5.2.89.

[21] Y. Dufresne, V. Moureau, E. Masi, O. Simonin, J. Horwitz, Simulation of a reactive fluidized bed reactor using CFD/DEM, in: Center for Turbulence Research Summer Program 2016, Vol. 26, 2016.

[22] L. Zhou, W. Lv, L. Bai, Y. Han, J. Wang, W. Shi, G. Huang, CFD-DEM study of gas-solid flow characteristics in a fluidized bed with different diameter of coarse particles, Energy Rep. 8 (2022) http://dx.doi.org/10.1016/j.egyr.2022.01.174.

[23] S. Tenneti, R. Garg, S. Subramaniam, Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres, Int. J. Multiph. Flow 37 (2011) http://dx.doi.org/10.1016/j.ijmultiphaseflow.2011.05.010.

[24] R. Beestra, M.V. der Hoef, J. Kuipers, Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres, AIChE J. 53 (2007) http://dx.doi.org/10.1002/aic.11065.

[25] B. Vowinckel, Incorporating grain-scale processes in macroscopic sediment transport models, Acta Mech. 232 (2022) http://dx.doi.org/10.1007/s00707-021-02951-4.

[26] A. Ozel, P. Fede, O. Simonin, Development of filtered Euler-Euler two-phase model for circulating fluidised bed: High resolution simulation, formulation and a priori analyses, Int. J. Multiph. Flow 55 (2013) http://dx.doi.org/10.1016/j.ijmultiphaseflow.2013.04.002.

[27] D. Dupuy, Y. Badran, R. Ansart, O. Simonin, Calibrating the frictional-pressure model from two-fluid simulation of fluidised beds in the defluidisation regime, Powder Technol. (2024) http://dx.doi.org/10.1016/j.powtec.2024.119776.

[28] N. Deen, M.V.S. Annaland, M.V. der Hoef, J. Kuipers, Review of discrete particle modeling of fluidized beds, Chem. Eng. Sci. 62 (2007) http://dx.doi.org/10.1016/j.ces.2006.08.014.

[29] D. Dupuy, R. Ansart, O. Simonin, Investigation of near-wall particle statistics in CFD-DEM simulations of dense fluidised beds and derivation of an Eulerian particle dynamic wall boundary condition, J. Fluid Mech. 982 (2024) http://dx.doi.org/10.1017/jfm.2024.36.

[30] M.-A. Chadil, S. Vincent, J.-L. Estivalezes, Accurate estimate of drag forces using particle-resolved direct numerical simulations, Acta Mech. 230 (2018) http://dx.doi.org/10.1007/s00707-018-2305-1.

[31] A. Esteghamatian, M. Bernard, M. Lance, A. Hammouti, A. Wachs, Micro/meso simulation of a fluidized bed in a homogeneous bubbling regime, Int. J. Multiph. Flow 92 (2017) http://dx.doi.org/10.1016/j.ijmultiphaseflow.2017.03.002.

[32] A. Esteghamatian, A. Hammouti, M. Lance, A. Wachs, Particle resolved simulations of liquid/solid and gas/solid fluidized beds, Phys. Fluids 29 (2017) http://dx.doi.org/10.1063/1.4979137.

[33] N. Deen, S. Kriebitzsch, M.V. der Hoef, J. Kuipers, Direct numerical simulation of flow and heat transfer in dense fluid-particle systems, Chem. Eng. Sci. 81 (2012) http://dx.doi.org/10.1016/j.ces.2012.06.055.

[34] S. Kriebitzsch, M.V. der Hoef, J. Kuipers, Fully resolved simulation of a gas-fluidized bed: A critical test of DEM models, Chem. Eng. Sci. 91 (2013) http://dx.doi.org/10.1016/j.ces.2012.12.038.

[35] K. Luo, J. Tuan, Z. Wang, J. Fan, Particle-resolved direct numerical simulation of gas-solid dynamics in experimental fluidized beds, AIChE J. 62 (2016) http://dx.doi.org/10.1002/aic.15186.

[36] J. Dou, L. Wang, W. Ge, J. Ouyang, Effect of mesoscale structures on solid phase stress in gas-solid flows, Chem. Eng. J. 455 (2023) http://dx.doi.org/10.1016/j.cej.2022.140825.

[37] P. Bagchi, M. Ha, S. Balachandar, Direct numerical simulation of flow and heat transfer from a sphere in a uniform cross-flow, J. Fluid Eng. 123 (2001) http://dx.doi.org/10.1115/1.1358844.
[38] Z.-G. Feng, E. Michaelides, Heat transfer in particulate flows with Direct Numerical Simulation, Int. J. Heat Mass Transfer 52 (2009) http://dx.doi.org/10.1016/j.ijheatmasstransfer.2008.07.023.

[39] H. Gan, J. Chang, J. Feng, H. Hu, Direct numerical simulation of the sedimentation of solid particles with thermal convection, J. Fluid Mech. 481 (2003) http://dx.doi.org/10.1017/S0022112003003938.

[40] Z. Yu, X. Shao, A. Wachs, A fictitious domain method for particulate flows with heat transfer, J. Comput. Phys. 217 (2006) http://dx.doi.org/10.1016/j.jcp.2006.01.016.

[41] C. Dan, A. Wachs, Direct numerical simulation of particulate flow with heat transfer, Int. J. Heat Fluid Flow 31 (2010) http://dx.doi.org/10.1016/j.ijheatfluidflow.2010.07.007.

[42] S. Tenneti, B. Sun, R. Garg, S. Subramaniam, Role of fluid heating in dense gas-solid flow as revealed by particle-resolved direct numerical simulation, Int. J. Mass. Transf. 58 (2013) http://dx.doi.org/10.1016/j.ijheatmasstransfer.2012.11.006.

[43] H. Tavassoli, S. Kriebitzsch, M.V. der Hoef, E. Peters, J. Kuipers, Direct numerical simulation of particulate flow with heat transfer, Int. J. Multiph. Flow. 57 (2013) http://dx.doi.org/10.1016/j.ijmultiphaseflow.2013.06.009.

[44] D. Gunn, Transfer of mass and heat to particles in fixed and fluidized beds, Int. J. Heat Mass Transfer 21 (1978) http://dx.doi.org/10.1016/0017-9310(78)90080-7.

[45] N. Deen, E. Peters, J. Padding, J. Kuipers, Review of direct numerical simulation of fluid-particle mass, momentum and heat transfer in dense gas-solid flows, Chem. Eng. Sci. 116 (2014) http://dx.doi.org/10.1016/j.ces.2014.05.039.

[46] N. Deen, J. Kuipers, Direct numerical simulation of fluid flow accompanied by coupled mass and heat transfer in dense fluid-particle systems, Chem. Eng. Sci. 116 (2014) http://dx.doi.org/10.1016/j.ces.2014.05.036.

[47] J. Derksen, Simulations of solid-liquid mass transfer in fixed and fluidized beds, Chem. Eng. J. 255 (2014) http://dx.doi.org/10.1016/j.cej.2014.06.067.

[48] Z.-G. Feng, S. Musong, Direct numerical simulation of heat and mass transfer of spheres in a fluidized bed, Powder Technol. 262 (2014) http://dx.doi.org/10.1016/j.powtec.2014.04.019.

[49] M.-A. Chadil, S. Vincent, J.-L. Estivalezes, Novel method to compute drag force and heat transfers for motions around spheres, Thermodyn. Interfaces Méc. Fluides 2 (2018) http://dx.doi.org/10.21494/ISTE.OP.2018.0299.

[50] M.-A. Chadil, S. Vincent, J.-L. Estivalezes, Gas-solid heat transfer computation from particle-resolved direct numerical simulations, Fluids 7 (2022) http://dx.doi.org/10.3390/fluids7010015.

[51] E. Thiam, Modélisation et Simulation Numérique Directe des Transferts de Chaleur dans les Écoulements Fortement Chargés en Particules (Ph.D. thesis), Université de Toulouse, 2018.

[52] C. Wen, Y. Yu, Mechanics of fluidization, Chem. Eng. Prog. Symp. Ser. 162 (1966).

[53] M. Haid, Correlations for the prediction of heat transfer to liquid-solid fluidized beds, Chem. Eng. Process. 36 (1997) http://dx.doi.org/10.1016/S0255-2701(96)04177-3.

[54] M. Lee, S. Kim, S. Kim, Experimental investigation of wall-to-bed heat transfer of carbon nanotubes in a fluidized bed, Int. J. Heat Mass Transfer 204 (2023) http://dx.doi.org/10.1016/j.ijheatmasstransfer.2023.123858.

[55] C. Calvin, O. Cueto, P. Emonot, An object-oriented approach to the design of fluid mechanics software, Math. Model. Numer. Anal. 36 (2002) http://dx.doi.org/10.1051/m2an:2002038.

[56] M. Grosso, G. Bois, A. Toutant, Thermal boundary layer modelling for heat flux prediction of bubbles at saturation: A priori analysis based on fully-resolved simulations, Int. J. Heat Mass Transfer 222 (2024) http://dx.doi.org/10.1016/j.ijheatmasstransfer.2023.124980.

[57] Y. Benarafa, O. Cioni, F. Ducros, P. Sagaut, RANS/LES coupling for unsteady turbulent flow simulation at high Reynolds number on coarse meshes, Comput. Methods Appl. Mech. Engrg. 195 (2006) http://dx.doi.org/10.1016/j.cma.2005.06.007.

[58] A. Toutant, E. Labourasse, O. Lebaigue, O. Simonin, DNS of the interaction between a deformable buoyant bubble and a spatially decaying turbulence: a priori tests for LES two-phase flow modelling, Comput. & Fluids 37 (2008) http://dx.doi.org/10.1016/j.compfluid.2007.03.019.

[59] R. Glowinsky, T. Pan, T. Hesla, D. Jospeh, J. Périaux, A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: Application to particulate flow, J. Comput. Phys. 169 (2001) http://dx.doi.org/10.1006/jcph.2000.6542.

[60] A. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation, J. Fluid. Mech. 271 (1994) http://dx.doi.org/10.1017/S0022112094001771.

[61] H.H. Hu, N. Patankar, M. Zhu, Direct numerical simulations of fluid-solid systems using the arbitrary Lagrangian-Eulerian technique, J. Comput. Phys. 169 (1992) http://dx.doi.org/10.1006/jcph.2000.6592.

[62] S. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, J. Comput. Phys. 100 (1992) http://dx.doi.org/10.1016/0021-9991(92)90307-K.

[63] J. Ritz, J. Caltagirone, A numerical continuous model for the hydrodynamics of fluid particle systems, Internat. J. Numer. Methods Fluids 30 (1999) http://dx.doi.org/10.1002/(SICI)1097-0363(19990830)30:8<1067::AID-FLD881>3.0.CO;2-6.

[64] S. Vincent, J.B. de Motta, A. Sarthou, J.-L. Estivalezes, O. Simonin, E. Climent, A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows, J. Comput. Phys. 256 (2014) http://dx.doi.org/10.1016/j.jcp.2013.08.023.

[65] S. Hamidi, A. Toutant, S. Mer, F. Bataille, Assessment of a coupled VOF-Front Tracking/DEM method for simulating fluid-particles flows, Int. J. Multiph. Flow 165 (2023) http://dx.doi.org/10.1016/j.ijmultiphaseflow.2023.104467.

[66] J. Delahaye, Jump conditions and entropy sources in two-phase systems. Local instant formulation, Int. J. Multiph. Flow 1 (1974) http://dx.doi.org/10.1016/0301-9322(74)90012-3.

[67] A. Toutant, Modélisation Physique des Interactions Entre Interfaces et Turbulences (Ph.D. thesis), Institut National Polytechnique de Toulouse, 2006.

[68] D. Legendre, R. Zenit, C. Daniel, P. Giraud, A note on the modelling of the bouncing of spherical droplets or solid spheres on a wall in viscous fluid, Chem. Eng. Sci. 61 (2006) http://dx.doi.org/10.1016/j.ces.2005.12.028.

[69] L. Verlet, Computer 'experiments' on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules, Phys. Rev. 159, http://dx.doi.org/10.1103/PhysRev.159.98.

[70] X. Fang, J. Tang, H. Luo, Granular damping analysis using an improved discrete element approach, J. Sound Vib. 308, http://dx.doi.org/10.1016/j.jsv.2007.07.034.

[71] A. Xu, B.-T. Li, Particle-resolved thermal lattice Boltzmann simulation using OpenACC on multi-GPUs, Int. J. Heat Mass Transfer 218 (2024) http://dx.doi.org/10.1016/j.ijheatmasstransfer.2023.124758.

[72] O. Tange, GNU parallel - The command-line power tool, USENIX Mag. 36 (1) (2011) 42-47, URL http://www.gnu.org/s/parallel.

[73] C.B. Sullivan, A. Kaszynski, PyVista: 3D plotting and mesh analysis through a streamlined interface for the visualization toolkit (VTK), J. Open Source Softw. 4 (37) (2019) 1450, http://dx.doi.org/10.21105/joss.01450.

[74] W. Ranz, W. Marshall, Evaporation from drops - Part 1, Chem. Eng. Prog. 48 (1952).

[75] S. Whitaker, 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. 18 (1972) http://dx.doi.org/10.1002/aic.690180219.

[76] Z.-G. Feng, E. Michaelides, A numerical study on the transient heat transfer from a sphere at high Reynolds and Peclet numbers, Int. J. Heat Mass Transfer 43 (2000) http://dx.doi.org/10.1016/S0017-9310(99)00133-7.

[77] A.T. Cate, C. Nieuwstad, J. Derksen, H.V. den Akker, Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity, Phys. Fluids 14 (2002) http://dx.doi.org/10.1063/1.1512918.

[78] R. Bird, W. Stewart, E. Lightfoot, Transport phenomena, 1960.

[79] A.A. Corona, Agitation des Particules dans un Lit Fluidisé (Ph.D. thesis), Institut National Polytechnique de Toulouse, 2008.

[80] A. Ozel, J.B. de Motta, M. Abbas, P. Fede, O. Masbernat, S. Vincent, J.-L. Estivalezes, O. Simonin, Particle resolved direct numerical simulation of a liquid-solid fluidized bed: Comparison with experimental data, Int. J. Multiph. Flow 89 (2017) http://dx.doi.org/10.1016/j.ijmultiphaseflow.2016.10.013.

[81] A.A. Corona, R. Zenit, O. Masbernat, Collisions in a liquid fluidized bed, Int. J. Multiph. Flow 37 (2011) http://dx.doi.org/10.1016/j.ijmultiphaseflow.2011.02.004.

[82] J. Richardson, W. Zaki, The sedimentation of a suspension of uniform sphere under conditions of viscous flow, Chem. Eng. Sci. 3 (1954) http://dx.doi.org/10.1016/0009-2509(54)85015-9.