# Numerical Turbulent Flow Analysis through a Rotational Heat Recovery System

Maxime Piton, Florian Huchet, Bogdan Cazacliu, Olivier Le Corre

► To cite this version:

Maxime Piton, Florian Huchet, Bogdan Cazacliu, Olivier Le Corre. Numerical Turbulent Flow Analysis through a Rotational Heat Recovery System. Energies, 2022, 15 (18), 18p. 10.3390/en15186792. hal-04103886

HAL Id: hal-04103886

https://hal.science/hal-04103886

Submitted on 23 May 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Numerical Turbulent Flow Analysis through a Rotational Heat Recovery System

Maxime Piton<fnref n="1" />, Florian Huchet <fnref n="1" />,* , Bogdan Cazacliu <fnref n="1" /> and Olivier Le Corre <fnref n="2" />

**Abstract:** Herein, hydrodynamic analysis from a large-eddy simulation in Couette–Taylor–Poiseuille (CTP) geometry is numerically investigated. The present geometry is inspired by a previous experimental work in which heat transport phenomena were investigated in a heat recovery system devoted to a rotary kiln facility. The streamwise and spanwise components of the velocity and the Reynolds stress tensor are firstly validated using an experimental benchmark. The effect of the axial flow rates is studied at a fixed rotational velocity. It is shown that the streamwise velocity component damps the vortex flow organization known in Couette–Taylor (CT) flow. The bulk region and its wall footprint are therefore characterized by various methods (spectral and statistical analysis, Q-criterion). It is shown that the turbulent kinetic energy of the streamwise component in the near-wall region is augmented leading to a multi-scale nature of turbulence.

**Keywords:** LES; Couette–Taylor–Poiseuille; turbulence; rotary kiln; waste heat

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

### ✓ check for updates

**Citation:** Piton, M.; Huchet, F.; Cazacliu, B.; Le Corre, O. Numerical Turbulent Flow Analysis through a Rotational Heat Recovery System. *Energies* **2022**, *15*, 6792. https://doi.org/10.3390/en15186792

Academic Editor: Dmitry Eskin

Received: 24 August 2022

Accepted: 15 September 2022

Published: 16 September 2022

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

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

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

## 1. Introduction

Waste heat is considered to be a by-product of the overall thermal processes related to the usage of fossil fuel (e.g., refinery plants), domestic heating (e.g., thermal power central), building materials manufacturing (e.g., kilns in a cement plant), or any other raw resources requiring heat treatment. For over a decade, it has attracted significant interest, since the re-use of excess heat, usually transferred to the surroundings, could improve existing equipment while giving the opportunity to introduce sustainable low-cost technologies. The benefits are multiple and could tackle major issues tied to:

- The global carbon footprint mitigation of numerous human activities [1].

- The social benefit due to diversified usage of the available energy power from solar [2] or biomass [3] resources.

According to the heat excess quality, the technologies devoted to its valorisation can be divided into three sub-categories. At a weak quality ($T < 100 \,^\circ C$), the heat is transferred into various exchangers for its direct re-use, while at a medium quality ($100 \,^\circ C < T < 350 \,^\circ C$), energy storage materials are often preferred. For a higher quality ($T > 350 \,^\circ C$), the electricity conversion in various thermodynamic cycles is a relevant solution, despite some expected improvements in the design of many involved heat exchangers [4].

The present paper deals with the hydrodynamic characterization of heat recovery exchangers devoted to rotational thermal facilities. Rotary kilns are one such example of equipment whose energy consumption is very important, for example, in the cement industry. Ranging between 3 GJ/T and 5 GJ/T according to the clinkerization process [5], the rotary kiln process is able to transform the limestone into calcine materials at high temperatures (~2000 °C). Consequently, energy optimization [5] and CO₂ sequestration [6] are of paramount importance to promote energy recovery in the rotary kiln industry (capture, transport, and storage of waste heat). This traditional equipment releases the

<footnote n="1">1 MAST-GPEM, University Gustave Eiffel, 44344 Bouguenais, France</footnote>

<footnote n="2">2 IMT Atlantique Département Systèmes Énergétiques et Environnement, UMR CNRS GEPEA, 44007 Nantes, France</footnote>

<footnote>* Correspondence: florian.huchet@univ-eiffel.fr</footnote>
required heat from the reactive combustion of fossil fuel resources. Slightly inclined (~3°), the kiln can reach 80 m in length and 6 m in diameter, while its rotation is relatively slow (from 2 to 10 rpm). Despite several intermediate processes such as the pre-heating of raw materials, the exergy and energy yields can reach, respectively, 39% and 59% [5,7]. The wall heat losses represent 61% of the overall waste heat [8]. Among this wall heat loss, 19% is imputed to the hot exhaust gases, 6% into the condenser, and 15% to the wall. Energy recovery from the wall is mentioned in the literature [9] as an alternative to ensure waste heat recovery. Most of the authors consider a second shell mounted to the drum as a potential technological solution [10,11]. Indeed, the medium grade energy of the wall (~150 °C) or the exhaust gases (~400 °C) require augmentation, or at least maintenance, during heat waste recovery. Consequently, only an efficient heat exchanger using a second shell would be capable of fulfilling the objective by confining a fluid between an inner moving cylinder and an outer cylinder, as shown in Figure 1.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Position of the problem: two concentric cylindrical surfaces of internal and external radius, Ri and Re, and of length, Lz, are inclined at an angle φ with gravity. The rotation is exerted from the inner cylinder (drawn here) at an angular velocity, ω. The heat density flux can be imposed from any wall surface (the inner cylinder is chosen here) while the fluid motion is flowing into the gap: $e = Re - Ri$.</figcaption></figure>

In that context, the present paper is devoted to the numerical flow characterization of a heat recovery system preliminary studied at the pilot scale [12]. This experimental test bench was sized at a scale of 1:5 with respect to an industrial rotary dryer [13], in the perspective of a mounting on its external wall to simulate a heat exchanger. Two methods of heat recovery were proposed by so-called "natural" and "forced" convection methods. The thermal flow regimes were explored as functions of the Reynolds and Richardson numbers. The amount of heat recovered into the system is representative of wall heat fluxes measured in industrial plants [14]. A dimensionless relationship from the measurements
relying on the Nusselt number, Nu, and the effective Reynolds number, Re<sub>eff</sub>, has been proposed (Nu ~ Re<sub>eff</sub><sup>a</sup>). The exponent value found, equal to a~0.6, has been compared to previous works which established heat transport in swirling decaying flows. Starting from this finding, the present work aims to numerically study the turbulence properties in the annular duct involving the rotating inner cylinder.

**2. Materials and Methods**

**2.1. *Numerical Strategy***

The flow within the annular duct, composed of an inner moving cylinder, has been relatively well documented in the fluid mechanics field since the historical paper of Taylor in 1923 [15] until the present [16]. The role and apparition of the instabilities in the Couette–Taylor flow are known to be complex, depending on numerous parameters (gap effect, initial condition, etc.). The fine characterization of the different turbulence structures from large- to micro-scales remains an active research question [16]. The self-sustained process of the turbulence is characterized by several flow regimes from Taylor vortex up to certain wavy shapes of the flow responsible for decreasing the wall shear stress [17]. The superposition of an axial flow called the Taylor–Couette–Poiseuille flow has also been investigated numerically [18] and experimentally in the literature. Certain experimental data [19,20] are often chosen for numerical validation where the RANS (Reynolds-averaged Navier–Stokes) method [21] including different types of turbulence closure models fails to validate the simulations. Therefore, two approaches have become the consensus in the scientific community:

- A derived RANS approach based on a second-order model, taking into account the curve wall effect [22];

- Large-eddy simulation, which offers significant accuracy [23,24] in solving the vortices dynamics at the fluid/wall interface.

Consequently, the LES OneEqEddy and RSM models developed by OpenFOAM software are chosen to solve the Navier–Stokes equation in CTP flow.

Three axial flow are investigated corresponding to a rotational number ranged between $N = 0.85$ and $N = 5.12$. The CTP flow is herein representative of the air-cooled heat exchanger experimentally studied [12]. The radius of the heat exchanger is identical to those of the test bench. The drum radius $R_i$ is 0.17 m and the external radius $R_e$ is 0.21 m. The thickness of the annular area between the two cylinders is $e = 0.04$ m for a study length of $L_z = 20\delta$ with $\delta = e/2$ (Figure 2). The numerical length $L_z$ corresponds here to an aspect ratio $\Gamma = \frac{L_z}{e} = 10$ and the shape factor $\eta = \frac{R_i}{R_e} = 0.809$.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Geometry studied.</figcaption></figure>

The considered fluid is air (density: $\rho = 1.2 \text{ kg.m}^{-3}$, dynamic viscosity: $\mu = 1.85 \times 10^{-5} \text{ kg.m}^{-1}.\text{s}^{-1}$). The rotational speed of the drum is equal to 50 tr.min<sup>-1</sup>.
($\omega R_i = 0.94 \text{ m.s}^{-1}$), and the axial flow rate in the exchanger is between 30 and 200 m$^3$.h$^{-1}$ ($\overline{U}_z = 0.21 - 0.95 \text{ m.s}^{-1}$). The studied velocity ranges correspond to a turbulent regime for CTP flow as reported in Table 1.

<table><caption>**Table 1.** Parameters studied in the numerical simulations for the different cases A, B, and C.</caption><thead><tr><th></th><th>A</th><th>B</th><th>C</th></tr></thead><tbody><tr><td>ω [tr.min<sup>-1</sup>]</td><td></td><td>50<sup>-</sup></td><td></td></tr><tr><td>Re<sub>ω</sub></td><td></td><td>7488<sup>-</sup></td><td></td></tr><tr><td>Q [m<sup>3</sup>.h<sup>-1</sup>]</td><td>31.4</td><td>94.3</td><td>188.6</td></tr><tr><td>Re<sub>a</sub></td><td>1462</td><td>4388</td><td>8776</td></tr><tr><td>N = Re<sub>ω</sub>/Re<sub>a</sub></td><td>5.12</td><td>1.71</td><td>0.85</td></tr></tbody></table>

## 2.2. Meshes Parameters

The computational accuracy increases with the number of meshes, while the computational time increases with the number of equations to solve directly related to the number of meshes. Therefore, there is a trade-off between the number of meshes, the computational accuracy, and the computational time. In LES simulations, the mesh size is relatively important compared the RANS approach. Large scales larger than the mesh size are solved directly, and smaller scales are modelled based on the assumption of a sub-grid turbulence model below the smaller mesh size. Consequently, it is necessary to adjust the mesh size in order to satisfy this assumption.

Large-scale structures of the flow are dependent on geometry and boundary conditions, while the behaviour of small structures at the dissipative scales is assumed to be self-similar and dependent only on the fluid viscosity $\nu$ and the amount $\varepsilon$ of the turbulent kinetic energy dissipation rate.

The turbulent Taylor microscale $\lambda$ represents the characteristic scale of the finest spatial fluctuations in velocity responsible for energy dissipation. This Taylor scale is used to define the dimension below which turbulent vortices behave universally and independently of external causes.

$$ \lambda = \sqrt{\frac{15 \nu u_i'^2}{\bar{\varepsilon}}} \qquad (1) $$

where $u_i'$ is the standard deviation of each velocity component. $\bar{\varepsilon}$ is the averaged turbulent kinetic energy dissipation rate. In a steady-state flow, $\varepsilon$ is the power required to be delivered to the large scales of turbulence in order to continuously maintain the turbulent agitation of the flow.

In CTP flow, the mesh is performed according to the available computational power in order to refine the mesh to be closer to the Taylor microscale $\lambda$. The mesh consists of approximately 16.5 million hexahedral cells with sizes between 0.65 mm and 2 mm. In coordinates scaled by the friction velocity $V_\tau$, the cell size is $\Delta r^+ < 1.5$, $\Delta z^+ < 15$ and $r \Delta \theta^+ < 12$ for the inner and outer wall, respectively.

$$ \Delta r^+ = (\Delta r)_\omega V_\tau / \nu \qquad (2) $$

$$ \Delta z^+ = \Delta z V_\tau / \nu \qquad (3) $$

With,

$$ V_\tau = \sqrt{\frac{\tau}{\rho}} \quad \text{where } \tau = \mu \left( \frac{\partial V_z}{\partial r} \right)_\omega \qquad (4) $$

The mesh specifications are summarized in Table 2. The scaled cell sizes of the mesh were compared to those of the literature on CTP flow [23,24], as seen in Table 2.
<table><caption>**Table 2.** Parameters of the numerical simulations used in LES.</caption><thead><tr><th></th><th>[23]</th><th>[24]</th><th>A</th><th>B</th><th>C</th></tr></thead><tbody><tr><td>η</td><td>0.50</td><td>0.89</td><td>0.809</td><td>0.809</td><td>0.809</td></tr><tr><td>N</td><td>0.21–0.86</td><td>1.49–6.71</td><td>5.12</td><td>1.71</td><td>0.85</td></tr><tr><td>Δr<sup>+</sup></td><td>0.27–0.35</td><td>0.47–0.86</td><td>1.23</td><td>1.24</td><td>1.24</td></tr><tr><td>R<sub>i</sub>Δθ<sup>+</sup></td><td>8.01–10.40</td><td>94.76–44.45</td><td>12.13</td><td>12.15</td><td>12.18</td></tr><tr><td>R<sub>e</sub>Δθ<sup>+</sup></td><td>13.86–15.80</td><td>96.05–45.27</td><td>15.21</td><td>14.98</td><td>15.01</td></tr><tr><td>Δz<sup>+</sup></td><td>9.93–22.95</td><td>26.95–48.82</td><td>14.95</td><td>15.28</td><td>15.32</td></tr><tr><td>N<sub>r</sub>, N<sub>θ</sub>, N<sub>z</sub></td><td>(65, 64, 128)</td><td>(65, 144, 130)–(91, 128, 182)</td><td>(123, 672, 200)</td><td>(123, 672, 200)</td><td>(123, 672, 200)</td></tr><tr><td>δ<sub>t</sub> (s)</td><td>-</td><td>-</td><td>4 × 10<sup>-4</sup></td><td>4 × 10<sup>-4</sup></td><td>3 × 10<sup>-4</sup></td></tr><tr><td>CFL</td><td>-</td><td>-</td><td>0.20</td><td>0.21</td><td>0.27</td></tr></tbody></table>

## 2.3. Boundary Conditions

The boundaries of the geometry study domain are two solid walls and cyclic boundary conditions. Unlike other types of boundary conditions, cyclic conditions are not boundaries. Using these conditions allows all properties to be repeated around an axis of symmetry. This cyclic boundary condition is applied between the inlet and outlet of the geometry to model the flow along the infinite in the z-direction. This approach has been previously used [23,24] for LES studies in CTP flow to obtain a fully developed flow between the two cylinders. At the walls, a no-slip condition is applied, while the rotational speed, $ωR_i$, of the inner wall is imposed.

## 3. Numerical Validation

Here, LES OneEqEddy and RSM models were initially compared to experiments performed using laser Doppler anemometry (LDA) [20] on a geometry having a shape factor of $η = 0.5$ and a wide aspect ratio of $Γ = 244$. These experimental results have already been used by other authors [22] for validation of their numerical approaches to extend their work to real operating conditions ($η = 0.961$, $Γ = 77$), including non-isothermal flow for a wide range of axial and rotational Reynolds numbers ($0 ≤ Re_a ≤ 4870$, $3744 ≤ Re_t ≤ 37,443$).

Such experiments [20] have provided a large database on mean velocity components and turbulent Reynolds stresses in CTP flow. Initially, the geometric parameters are defined in a dimensionless form with $r^* = (r - R_i)/(R_e - R_i)$ and $z^* = z/L$. Under these conditions, the inner cylinder is found to be $r^* = 0$ and the outer cylinder is $r^* = 1$. The velocity components are also scaled by the mean flowing velocity $\bar{U}_z$ for the streamwise component $U_z^* = U_z/\bar{U}_z$, and by the cylinder rotation speed $ωR_i$ for the spanwise component $U_θ^* = U_θ/ωR_i$. The Reynolds stress components $R_{θθ}$ and $R_{zz}$ are also scaled by the mean flow velocity $\bar{U}_z$.

Figures 3 and 4 show the comparison of the numerical results of the RSM model and the LES modelling with the experimental data. In contrast with the RSM simulation, the LES simulation better predicts the average velocity, $U_θ^*$, in the annulus. The LES simulations appear to take better account of the rotational effects than the RANS simulation (Figure 3).

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Steamwise and spanwise velocity profiles. Comparison between numerical models (RSM and LES) and the experimental measurements of Escudier and Gouldson (1995) [20].</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Normalized Reynolds stress components $R_{\theta\theta}^{0.5*}$ and $R_{zz}^{0.5*}$. Comparison between numerical models (RSM and LES) and the experiments of Escudier and Gouldson (1995) [20]. Adapted with permission from Ref. [20]. Copyright (1995), with permission from Elsevier.</figcaption></figure>

The numerical results of the LES simulations obtained in OpenFOAM are close to the experimental data [20]. However, it is important to note the presence of a near-wall peak ($r^* = 0$) on the experimental turbulent component $R_{\theta\theta}^{1/2}$. This peak was not reproduced by the numerical simulations and has not been observed in the numerical simulations available in the literature [21,22,24]. The presence of this peak may be related to the experimental difficulty in LDA of measuring the velocity in a thin boundary layer. The use of appropriate experimental methods, such as the electrochemical method, in the near-wall region is fundamental for wall shear stress determination [25].

# 4. Results and Discussion

The predictions of the LES simulations have been extended to the study case in order to obtain a better knowledge of the hydrodynamics of the CTP flow in the industrial application [12]. In this previous work, only heat transport characterization was performed on the exchanger geometry, characterized by a shape factor of $\eta = 0.809$ and an aspect ratio of $\Gamma = 22.5$.

## 4.1. Mean Flow Fields

The rotation parameter $N$, defined in Table 1, is used in the study of CTP flow to compare inertial effects related to the rotational and axial motions of the flow [19,20,23,24,26,27]. As shown in Figure 5, the trajectory of the fluid particles behaves differently according to the $N$ values. The axial flow effect reduces their path lengths from the inlet to the outlet.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Representation of the streamlines at the centre of the annular geometry for $N = 5.12$ case (A), $N = 1.71$ case (B), and $N = 0.85$ case (C).</figcaption></figure>

The instantaneous velocity fields of the streamwise and spanwise components are shown in Figures 6 and 7, to visualize the turbulent structures of the flow. The radial velocity component has not been shown because it tends to zero throughout the system. The figures show that the flow is a combination of CT flow in the near-wall region of the rotating cylinder and Poiseuille flow in the annular space. The combination of the CT flow in the near-wall region superimposed on a Poiseuille flow explains the helical shape of the
flow streamlines around the inner cylinder. Nevertheless, the effect of rotation decreases rapidly on the mean flow with the increase in the axial flow.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Velocity fields of the streamwise (a) and spanwise (b) components in a cross-section of the annular geometry (the axis of rotation is on the right side of each figure and the flow moves from bottom to top).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Velocity fields of the spanwise component in a section of the annular geometry for *N* = 5.12 case (A), *N* = 1.71 case (B), and *N* = 0.85 case (C).</figcaption></figure>

Figure 8 shows the radial distributions of the axial mean velocity components $U_z^*$ and tangential $U_\theta^*$ for three values of *N* for $z^* = 0.5$ (see Figure 9).
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Velocity profiles of the streamwise (a) and spanwise (b) components obtained by the LES OneEqEddy model of OpenFOAM.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Position of the area of measurement in the studied geometry.</figcaption></figure>

The mean axial velocity profile is closer to a turbulent Blasius profile than a Poiseuille flow; the streamwise velocity profile loses its symmetry plan with the increase in the rotational effect at $N = 5.12$. The axial flow rate has an influence on the velocity profile of the spanwise component, $U_\theta^*$. The spanwise velocity component is composed of a rotating central flow, with two thin near-wall boundary layers. This type of profile corresponds to the velocity profiles typically observed in a turbulent regime in rotating flow configurations. In this case study, the central region rotates at approximately 40% of the inner cylinder velocity and the increase in the axial flow tends to stabilize the flow.

## 4.2. Reynolds Stress Tensor

The Reynolds stress tensor of the turbulent motion obtained by the LES simulation was analysed at $z^* = 0.5$ (Figure 9) for different values of the parameter, $N$, to study its impact on the turbulent fluctuations.

The study of the distribution of the six Reynolds stresses normalized by the average axial flow velocity $R_{zz}^* = R_{zz}^{1/2}/\bar{U}_z$ and $R_{\theta z}^* = R_{\theta z}^{1/2}/\bar{U}_z$, is displayed in Figure 10. From the Reynolds stress distribution, the turbulence rate appears to be mainly concentrated in the near-wall region of the boundary layer of the flow with relatively high turbulence peaks. The normal components $R_{\theta\theta}^*$ and $R_{zz}^*$ are lower in the bulk of the flow, while $R_{rr}^*$ is higher in the centre of the flow. The diagonal components of the Reynolds tensor of turbulent motion form the turbulent kinetic energy $= 1/2(R_{rr} + R_{\theta\theta} + R_{zz})$. From Figure 10, the $R_{\theta\theta}$ term appears to be dominant compared to the axial component of the turbulent fluctuations $R_{zz}$.
Therefore, the turbulent kinetic energy increases with N and the axial motion contributes to a large transfer of turbulent kinetic energy.

<figure><img src="image_13.png" /><figcaption>Figure 10. Distribution of the six components of the Reynolds stress tensor. The diagonal components of the dimensionless tensor $R_{rr}^{1/2*}$, $R_{\theta\theta}^{1/2*}$ and $R_{zz}^{1/2*}$ are represented by (a,c,e) and the cross components $R_{r\theta}^{*}$, $R_{rz}^{*}$ and $R_{\theta z}^{*}$ are represented by (b,d,f).</figcaption></figure>

The crossed components $R_{rz}^{*}$ and $R_{\theta z}^{*}$ show strong stress variations at the inner and outer cylinder walls. The turbulence intensity between the two cylinders increases with N and appears to be more unstable at the rotating cylinder boundary layer. This can also be attributed to the small effects of the radius of curvature.

The energy redistribution induced by the effect of the radius of curvature is defined by the parameter $K^*$ [28]. This parameter measures the relative contribution of the turbulent intensity by comparing the turbulent contributions normal to the mean flow.

$$K^* = \frac{2\overline{u'_z^2}}{(\overline{u'_r^2} + \overline{u'_\theta^2})} \qquad (5)$$

From Figure 11, the effect of the radius of curvature appears to be less pronounced than that observed in preceding works [24,25]. It appears that the value of $K^*$ is larger at the outer cylinder wall than at the inner wall. This leads to the conclusion that the energy
transferred near the wall by the streamwise component $\overline{u'_{z}^{2}}$ to the other two components is enhanced at the rotating inner cylinder.

<figure><img src="image_14.png" /><figcaption>Figure 11. Radial distribution of the energy redistribution parameter obtained by the LES OneEqEddy model of OpenFOAM.</figcaption></figure>

The parameter $K^*$ is greater than 1 for $N_C = 0.85$ and $N_B = 1.71$, and close to 1 for $N_A = 5.12$. Therefore, the contribution of the axial turbulence intensity of the turbulent kinetic energy $\overline{u'_{z}^{2}}$ is larger than the contribution of the other two normal components of the flow. With the increase in the rotational speed contribution, the values of $K^*$ decrease and are significantly lower than the values obtained by Chung et al. (2002) [26] in stationary cylinders. This decrease in the values of $K^*$ can be explained by the increase in the contribution of the turbulent intensity $\overline{u'_{\theta}^{2}}$ on the flow with rotation. As expected, the contributions of $\overline{u'_{z}^{2}}$ decrease strongly with increases in the parameter $N$. For large values of $N$, the rotational effect dominates over the axial flow.

### 4.3. Boundary Layer Properties

For confined flows, the wall treatment is important because it can develop its own turbulence properties. The flow is in a laminar regime very close to the wall in the viscous sublayer, and moving away, the turbulence increases rapidly because of the high velocity gradients. The turbulent boundary layer can be divided into three areas: a viscous sublayer near the wall, a transition zone where turbulence and viscosity are equally important, and a zone dominated by turbulence.

Figure 12 shows the radial and normalized variation of the effective velocity $U^+ = U/V_τ$ as a function of the wall coordinate $r^+ = r_p V_τ / \nu$, where $r_p$ is the distance between the wall and the measurement. The velocity $U$ and the frictional velocity $V_τ$ are defined as follows:

$$U = \sqrt{U_z^2 + (\omega R_i - U_\theta)^2} \qquad (6)$$

$$V_τ = \nu^{1/2} \left( \left( \frac{\partial U_θ}{\partial r} \right)_p^2 + \left( \frac{\partial U_z}{\partial r} \right)_p^2 \right)^{1/4} \qquad (7)$$

In confined flow, specific treatment is required in addition to the turbulence model used in the area where the flow is developed. OpenFOAM proposes two approaches to solve the boundary layer: a wall law or a wall model. The wall law is founded on an analytical formulation of the wall region, while the wall model resolves the overall boundary layer. The resolution of the boundary layer requires a refined mesh in the near wall and is therefore expensive in terms of computation time. Thus, the present simulations
were performed by a wall model in order to solve and study the boundary layer of the flow.

<figure><img src="image_15.png" /><figcaption>Figure 12. Logarithmic boundary layer profiles for the rotating wall.</figcaption></figure>

Figure 12 shows the hydrodynamic structure of the boundary layer as a function of the rotation parameter N. The wall modelling at the viscous sublayer $r^+ < 5$ appears well resolved, $r^+ = u^+$. In this viscous sublayer, the viscous stress is greater than the turbulent stress. The size of the transient region between the viscous sublayer and the outer region $5 < r^+ < 30$ is larger or smaller depending on the magnitude of the rotation parameter N. In this region, turbulent kinetic energy generation and dissipation are dominant. Finally, the logarithmic region between $30 < y^+ < 200$ is defined by:

$$U^+ = \frac{1}{\kappa} \ln r^+ + B \qquad (8)$$

In this equation, $\kappa$ represents the universal von Kármán constant (typically 0.41) and B a constant close to 5 for a turbulent flow pipe geometry. In this study case, the von Kármán constant usually found in confined turbulent flow is not found. In the literature, the effect of the radius of curvature is known to decrease the slope of the transition zone. According to recent work [29], the rotation of the flow significantly affects the turbulence by introducing a nonlinear term into the logarithmic wall law. This observation was examined by Poncet et al., (2014) [24] in CTP flow and more recently by Belkadi et al., (2022) [30] in CT flow.

These results confirm the intensification of the hydrodynamic structures due to rotational effects. These significantly decrease the size of the boundary layer, ending up below that encountered in a pipe flow. However, the increase in the axial flow has the effect of enlarging the hydrodynamic boundary layer. The hydrodynamic boundary layer approaches the size of the boundary layer encountered in the case of pipe flow (case C, Figure 12).

### 4.4. Energy Cascade in the Near-Wall Region

Spectral analyses on the different cases of the Taylor-Couette-Poiseuille flow show a transfer of kinetic energy from small to large wave numbers. This evolution is called an energy cascade.

This energy cascade is calculated from the velocity fluctuations $u'_z$. The frequency domain transition of the signal recorded during the numerical experiments is ensured by a discrete Fourier transformation. The power spectral density was obtained from a series of operations to ensure a relevant statistical representation of a random signal. Based on signal processing theory and applied to turbulent wall-bounded flows [31]:

$$E^i(U) = \int_0^N u'_z(t)_N \times F_{ha}(t) \times \exp(-j2\pi ft)dt \qquad (9)$$
$$E(U) = \frac{\sum_{i=1}^{N} E^{i}(U)}{N} \qquad (10)$$

The principle is based on the extraction of the fluctuating signal from N blocks (N = 41), N representing the number of blocks into which the signal has been divided, each block having the same number of points, $N_e = 4048$. The number of points $N_e$ depends on the temporal resolution of the study. Finally, the signal is then multiplied by a Hanning’s window $F_{ha}(t)$ in order to avoid the edge effects caused by the Fourier transform on a finite signal.

$$F_{ha}(t) = 0.5 \left( 1 + \cos \frac{\pi t}{N_e} \right) \qquad (11)$$

Figure 13 shows the energy cascade for three rotation numbers N at the rotating wall. Three distinct zones of the $E(U)$ energy spectrum are exhibited, which were interpreted by such turbulence pioneers of the last century as Kolmogorov [32] or Hinze [33].

<figure><img src="image_16.png" /><figcaption>Figure 13. Spectral analysis obtained for different values of N.</figcaption></figure>

The first zone is the zone at large scales, where turbulence production takes place. The large-scale structures are related to the mean flow fields which obey the boundary conditions of the flow system. The second zone is the so-called “inertial zone” where a large range of micro-scales of turbulence co-exist. Its decrease is still the subject of many studies and is compared to the model of an isotropic nature of turbulence at a slope equal to -5/3. In this zone, the energy spectrum depends only on k and $\epsilon$. Finally, the third zone is called the dissipation zone. It is very small in our case and consists of the smallest scales, where turbulent kinetic energy is dissipated by viscous effects.

### 4.5. Turbulent Structures

This section presents the transfer phenomena at the moving wall ($r = 0$), due to the appearance of vortex structures in the flow. These vortex structures enhance the momentum transfer in the near-wall region. In numerical fluid mechanics, there are several parameters for the detection of vortex structures. Among them, the Q-criterion has been chosen for the 3D visualization of coherent structures and the vorticity criterion $\omega_{ij}$ to track the presence of vortex structures in the boundary layer.

#### 4.5.1. Q-Criteria

The visualization of coherent structures was obtained by using the Q-criterion proposed by Hunt et al. (1988) [34]. This criterion uses the second variant of the velocity gradient tensor defined by:

$$Q = -\frac{1}{2}(\Omega_{ij}\Omega_{ij} - S_{ij}S_{ij}) \qquad (12)$$
With,

$$ \Omega_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} - \frac{\partial u_j}{\partial x_i} \right) \quad \text{and} \quad S_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \qquad (13) $$

$\Omega_{ij}$ measures the rotation of a portion of fluid and $S_{ij}$ its shear. Thus, regions for which the Q-criterion is positive characterize regions where rotation dominates shear, thus allowing visualization of coherent vortices. The Q-criterion is commonly used in the literature to identify small vortex structures.

The vortex structures in the flow were thus identified by the positive iso-surfaces of the Q-criterion. The combination of Taylor–Couette near-wall flow superimposed on Poiseuille flow explains the presence of the vertical helical structures wound along the wall (Figure 14).

<figure><img src="image_17.png" /><figcaption>Figure 14. Iso-contours of the Q-criterion in the entire computational domain. (a,d) represent the iso-contours for case A, (b,e) for case B, and (c,f) for case C on the inner cylinder (a–c), as well as on the outer cylinder (d–f). (a) $Q = 200$, (b) $Q = 500$, (c) $Q = 1000$.</figcaption></figure>

The results of the three-dimensional LES simulations highlight the presence of unstable coherent structures at the two opposite boundary layers. These structures appear as coiled spirals along the length of the cylinders. At the inner wall, the vortex structures are coiled in the direction opposite to the rotation of the cylinder, while they are coiled in the direction of the rotation of the cylinder at the near outer wall. With the increase in the number of rotations $N$, the structures tilt progressively in the tangential direction. The angle formed on the inner cylinder in the tangential direction is equal to $-18^{\circ}$ for $N = 5.12$ (case A), $-55^{\circ}$ for $N = 1.71$ (case B), and $-72^{\circ}$ for $N = 0.85$ (case C).
### 4.5.2. Vorticity

To attest to the presence of vortices in the boundary layer, the $(r, z)$ mapping of the vorticity vector is presented in Figure 15. The vorticity vector is defined by:

$$ \omega_{ij} = \left( \frac{\partial u_i}{\partial x_j} - \frac{\partial u_j}{\partial x_i} \right) \vec{u}_z \quad (14) $$

<figure><img src="image_18.png" /><figcaption>Figure 15. Vorticity fields (the axis of rotation is on the right side of each figure and the flow moves from the bottom to the top) superimposed on the skin-friction coefficient $C_f$ for the study cases (A-C), and a Couette-Taylor flow. Increasing the parameter $N$ leads to an increase in the vorticity $\omega_{ij}$ and a reduction in the coefficient $C_f$.</figcaption></figure>
Thus, the centres of the vortices can be detected by a local maximum of the vorticity modulus. For each case, at both walls, we find an alternation of counter-rotating vortices highlighted by $\omega_{ij} > 0$ and $\omega_{ij} < 0$.

The skin-friction coefficient $C_f(z)$ characterizing the wall shear stress is estimated at a time $t$, between $r^* = 0$ and the first cell of the mesh located at $r^* = 0.03$ along the $z^*$ axis.

$$C_f = \frac{\tau_p}{1/2\rho U_m^2} \qquad (15)$$

$$\tau_p = \sqrt{\tau_z^2 + \tau_\theta^2} \qquad (16)$$

where, $\tau_p$ is the wall shear stress, $\tau_z = \mu \frac{\partial U_z}{\partial r} \Big|_{r=0}$ and $\tau_\theta = \mu \frac{\partial U_\theta}{\partial r} \Big|_{r=0}$.

Temporal variations in the skin-friction coefficient $C_f$ are observed in the near-wall region. These intermittent fluctuations are explained by the presence of counter-rotating vortices in the boundary layer (Figures 15 and 16). The magnitude of the skin-friction coefficient $C_f(t)$ increases with the rotation parameter $N$ (Figure 17). This is inferred by the observation of a greater amount of energy exchanged between the counter-rotating vortices ejecting some of the fluid located in the boundary layer towards the bulk flow. However, this energy exchanged between the boundary layer and the core appears to decrease with increasing flow rate.

<figure><img src="image_19.png" /><figcaption>Figure 16. Schematic representation of the counter-rotating vortices at the wall of the Taylor-Couette flow.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 17. Skin-friction coefficient $C_f$ versus time and its average value $\overline{C_f}$ as function of the $N$ parameter.</figcaption></figure>
Figure 17 (left-hand) shows the skin-friction coefficient $C_f(t)$ for the study cases A, B, and C estimated as a function of time $t$ between $r^* = 0$ and $r^* = 0.03$ for $z^* = 0.5$. The magnitude of the skin-friction coefficient is again larger for high values of $N$ (case A) and decreases with increasing axial flow and decreasing parameter $N$ (cases B and C). However, the wavelengths are reduced which can be explained by an increase in the turbulent kinetic energy transport of counter-rotating vortices in the near-wall region.

Finally, Figure 17 (right-hand) shows the average friction coefficient as a function of the rotation parameter $N$. Thus, a slight increase in the skin-friction coefficient $\overline{C}_f$ is observed with parameter $N$, decreasing with the increase in the axial flow rate.

## 5. Conclusions

The present paper aims to characterize flow through a rotational heat recovery system used for heat transport in the frame of energy harvesting in the rotary kiln industry. The importance of better understanding the role of turbulent eddies in the regime of forced convection is a necessity at the pilot scale and the industrial scale. Solving unsteady 3D flow will help to improve process modelling in transient conditions. Consequently, the large-eddy simulation technique (LES OneEqEddy) with a no-slip condition was chosen to solve the flow in Couette–Taylor–Poiseuille flow geometry. Rigorous validation was first employed to compare this numerical method to the RANS method and the existing results of the literature in terms of mean velocity profiles and Reynolds stress tensors.

The LES was then applied to the case of interest, characterized by a shape factor and an aspect ratio respectively equal to 0.809 and 10. Three axial flow rates were investigated at a constant rotational velocity of the inner cylinder equal to 50 rpm, their dimensionless rotation parameters being equal to $R_A = 5.12$, $R_B = 1.71$, and $R_C = 0.85$.

The main results are the following:

* Streamwise components of the velocity exhibit a Blasius "type" flow profile, while the spanwise components are characterized by three different regions; the inner boundary layer, the central region, and the external boundary layer.

* The distribution of the Reynolds stress tensor is composed of two main normal components in the streamwise and the spanwise directions, the radial one being slightly less important. The crossed components are negligible compared to the normal ones. Whatever the components, the more elevated values are located in the near-wall location where the turbulence fluctuating rates are the more pronounced.

* The fluctuating rate of the turbulence is much more important in the axial direction than the two other directions for the $R_B$ and $R_C$ cases, while in the $R_A$ case, its contribution is similar to the two others.

* A wall model solving the overall boundary layer is used in the LES algorithm resulting in the logarithmic profiles. The wall law is verified in the viscous sublayer, while the profile in the transient region does not fit with the von Kármán constants established for pipe flow. The thickness of the boundary layer is reduced as the rotational velocity of the inner cylinder is raised compared to the axial flow rate.

* Spectral analysis was applied to the fluctuating streamwise velocity in the near-wall region. Turbulence spectra are observed in the frequency domain with three distinct areas of production, inertia, and dissipative zones. The inertia zone is characterized by a strong anisotropy far from the $-5/3$ slope.

* The turbulence scales were considered by using the Q-criterion and the vorticity in the near-wall region. The flow anisotropy is confirmed by the helix shape of the vortex flow structures. The angle of tilt of the structures is progressively reduced as the axial flow rates are raised at the inner and outer cylinders.

* The skin-friction coefficients were assessed, showing a large amplitude imputed to the vortex flow organization. A comparison with Couette–Taylor flow results is provided, which explains the nature of the amplitude linked to the ejection of the fluid between the vortex cells toward the bulk region. The averaged and fluctuating values of the
skin-friction are reduced as the axial flow rates increase; damping the counter-rotating vortex effect.

The present numerical results prove that the axial flow regime plays an important role in turbulence characteristics, which could explain the scarcity found in the Nusselt numbers correlations ($Nu \sim Re_{eff}^{0.6<a<0.8}$) in no-isothermal conditions for a forced convection regime. At low rotational parameters, the streamwise component of the turbulent kinetic energy modifies the momentum exchange in the near-wall region, and the axial flow rate disrupts the vortex generation and gives rise to eddies of different sizes.

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

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

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

## References

1. Viklund, S.B.; Johansson, M.T. Technologies for utilization of industrial excess heat: Potentials for energy recovery and CO₂ emission reduction. *Energy Convers. Manag.* **2014**, *77*, 369–379. [CrossRef]

2. Kongtragool, B.; Wongwises, S. A review of solar-powered Stirling engines and low temperature differential Stirling engines. *Renew. Sustain. Energy Rev.* **2003**, *7*, 131–154. [CrossRef]

3. Gonzales, A.; Riba, J.R.; Puig, R.; Navarro, P. Review of micro and small-scale technologies to produce electricity and heat from Mediterranean forest' wood chips. *Renew. Sustain. Energy Rev.* **2015**, *43*, 143–155. [CrossRef]

4. Nakhchi, M.E.; Hatami, M.; Rahmati, M. Experimental evaluation of performance intensification of double-pipe heat exchangers with rotary elliptical inserts. *Chem. Eng. Process.-Process Intensif.* **2021**, *169*, 108615. [CrossRef]

5. Mungyeko Bisulandu, B.-J.R.; Huchet, F. Rotary kiln process: An overview of physical mechanisms, models and applications. *Appl. Therm. Eng.* **2022**. Submitted.

6. Swann, P.B.; Russel, H.; Jahn, I.H. Taylor-Couette-Poiseuille flow heat transfer in a high Taylor number test rig. *J. Glob. Power Propuls. Soc.* **2021**, *5*, 126–147. [CrossRef]

7. Sögüt, Z.; Oktay, Z.; Karakoç, H. Mathematical modeling of heat recovery from a rotary kiln. *Appl. Therm. Eng.* **2010**, *30*, 817–825. [CrossRef]

8. Engin, T.; Ari, V. Energy auditing and recovery for dry type cement rotary kiln systems-A case study. *Energy Conserv. Manag.* **2005**, *46*, 551–562. [CrossRef]

9. Mujumdar, A.; Ranade, A. Modeling of Rotary Cement Kiln: Applications to Reduction in Energy Consumption. *Ind. Eng. Chem. Res.* **2006**, *45*, 2315–2330. [CrossRef]

10. Caputo, A.C.; Pelegagge, P.M.; Salini, P. Performance modeling of radiant heat recovery exchangers for rotary kilns. *Appl. Therm. Eng.* **2011**, *31*, 2578–2589. [CrossRef]

11. Mastorakos, E.; Massias, A.; Tsakiroglou, C.D.; Goussi, D.A.; Burganos, V.N.; Payatakes, A.C. CFD prediction for cement kilns including flame modelling heat transfer and clinker chemistry. *Appl. Math. Model.* **1999**, *23*, 55–76. [CrossRef]

12. Huchet, F.; Piton, M.; Del Barrio, A.; Le Corre, O.; Cazacliu, B. Air-cooling heat exchanger applied to external rotary kiln wall in forced and natural draft. *Energy Convers. Manag.* **2017**, *154*, 517–525. [CrossRef]

13. Leguen, L.; Huchet, F.; Dumoulin, J.; Baudru, Y.; Tamagny, P. Convective heat transfer analysis in aggregates rotary drum reactor. *Appl. Therm. Eng.* **2013**, *54*, 131–139. [CrossRef]

14. Le Guen, L.; Huchet, F. Thermal imaging as a tool for process modelling: Application to a flight rotary kiln. *Quant. InfraRed Thermogr. J.* **2020**, *17*, 79–95. [CrossRef]

15. Taylor, G.I. Stability of viscous liquid contained between two rotating cylinders. *Philos. Trans. R. Soc. Lond.* **1923**, *223*, 289–343.

16. Nguyen, D.B.; Monico, R.O.; Chen, G. A Visualization Framework for Multi-scale Coherent Structures in Taylor-Couette Turbulence. *IEEE Trans. Vis. Comput. Graph.* **2021**, *27*, 902–912. [CrossRef] [PubMed]

17. Razzak, A.; Akter, S.; Hannan, M.A. The influence of wavy vortex in reducing wall shear stress. *IOP Conf. Ser. Mater. Sci. Eng.* **2020**, *788*, 012093. [CrossRef]

18. Mehrez, I.; Gheith, R.; Aloui, F.; Nasrallah, S.B. Theoretical and numerical study of Couette-Taylor flow with an axial flow using lattice Boltzmann method. *Int. J. Numer. Meth. Fluids* **2019**, *90*, 427–441. [CrossRef]

19. Nourri, J.M.; Whitelaw, J.H. Flow of Newtonian and non-Newtonian fluids in a concentric annulus with rotation of inner cylinder. *J. Fluid Eng.* **1994**, *116*, 821–827. [CrossRef]
20. Escudier, M.P.; Gouldson, I.W. Concentric annular flow with centerbody rotation of a Newtonian and a shear-thinning liquid. *Int. J. Heat Fluid Flow* **1995**, *16*, 156–162. [CrossRef]

21. Naser, A. Prediction of Newtonian and non-Newtonian flow through concentric annulus with centerbody rotation. In *Proceedings of the 1st International Conference on CFD in Mineral and Metal Processing and Power Generation*, Melbourne, Australia, 13–15 December 1996.

22. Poncet, S.; Haddadi, S.; Viazzo, S. Numerical modeling of fluid flow and heat transfer in a narrow Taylor-Couette-Poiseuille system. *Int. J. Heat Fluid Flow* **2011**, *32*, 128–144. [CrossRef]

23. Chung, S.Y.; et Sung, H.J. Large-eddy simulation of turbulent flow in concentric annulus with rotation of an inner cylinder. *Int. J. Heat Fluid Flow* **2005**, *26*, 191–203. [CrossRef]

24. Poncet, S.; Viazzo, S.; Oguic, R. Large eddy simulation of Taylor-Couette-Poiseuille flows in a narrow-gap system. *Phys. Fluids* **2014**, *26*, 105108. [CrossRef]

25. Havlica, J.; Kramolis, D.; Huchet, F. A revisit of the electro-diffusional theory for wall shear stress measurements. *Int. J. Heat Mass Transf.* **2021**, *165*, 120610. [CrossRef]

26. Chung, S.Y.; Rhee, G.H.; Sung, H.J. Direct numerical simulation of turbulent concentric annular pipe flow. Part. 1: Flow field. *Int. J. Heat Fluid Flow* **2002**, *23*, 426–440. [CrossRef]

27. Hadziabdic, M.; Hanjalic, K.; Mullyadzhanov, R. LES of turbulent flow in a concentric annulus with rotating outer wall. *Int. J. Heat Fluid Flow* **2013**, *43*, 74–84. [CrossRef]

28. Lee, M.J.; Kim, J.; Moin, P. Structure of turbulence at high shear rate. *J. Fluid Mech.* **1990**, *216*, 561–583. [CrossRef]

29. Singh, H.; Suazo, C.A.; Liné, A. Log law of the wall revisited in Taylor-Couette flows at intermediate Reynolds numbers. *Phys. Rev. E* **2016**, *94*, 053120. [CrossRef]

30. Belkadi, H.; Laghouati, Y.; Sobolik, V.; Oualli, H.; Bouabdallah, A. Axial Flow Effect on the Stability of Circular Couette Flow. *J. Appl. Fluid Mech.* **2022**, *15*, 271–281.

31. Huchet, F.; Legentilhomme, P.; Legrand, J.; Montillet, A.; Comiti, J. Unsteady flows in milli- and microsystems: Analysis of wall shear rate fluctuations. *Exp. Fluids* **2011**, *51*, 597–610. [CrossRef]

32. Kolmogorov, A.N. A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number. *J. Fluid Mech.* **1962**, *13*, 82–85. [CrossRef]

33. Hinze, J.O. Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes. *AIChE J.* **1955**, *1*, 289–295. [CrossRef]

34. Hunt, J.C.R.; Wray, A.A.; Moin, P. *Eddies, Streams, and Convergence Zones in Turbulent Flows*; Center for Turbulence Research Report CTR-S88; NASA: Washington, DC, USA, 1988; p. 193.