Three-Dimensional Computational Fluid Dynamics Analysis of Convective Heat Transfer From a Heated Horizontal Cylinder Rotating in Air: From Laminar to Turbulent Flow Open Access

A 3D computational fluid dynamics model is developed to reproduce the results of previous experiments and to investigate the correlation between Nusselt numbers and convection heat transfer phenomena surrounding an isothermal rotating cylinder. The simulation is conducted in a quiescent air domain and a fixed Grashof number of 2.32 × 10⁸ for a horizontal cylinder placed in air with rotational speeds ranging from 2.43 to 103.22 RPM. The effects of buoyancy-induced flows and the rotational Reynolds number Re_r on convective heat transfer characteristics are investigated. At low Re_r, buoyancy-driven Rayleigh–Bénard convection dominates, forming vertically extended thermal plumes obstructing heat convection on the upper side of the cylinder, leading to lower Nusselt number in these regions. As Re_r increases, rotational effects intensify, flow plumes merge with the cylinder surface and thicken the thermal boundary layers, on the other hand enhancing turbulent mixing, thus ultimately improving heat transfer. The circumferential Nusselt number distribution further highlights that plume formation lowers Nusselt numbers on the descending side, while heat transfer is enhanced along the axial direction toward the cylinder ends, where the thermal boundary layer thickness gradually decreases.

For centuries, rotating cylinders have been a significant topic of interest in the engineering community. Most machinery with rotating components eventually encounters convective heat transfer problems unique to heat convection from static surfaces. These problems often appear in tandem with thermal design considerations of flywheels, turbine rotors, rotary heat exchangers, rotary kilns and various other industrial applications.

An ample amount of experimental research has been conducted over the last several decades regarding the heat transfer phenomena surrounding heated rotating cylinders. Anderson and Saunders [1] carried out experiments to investigate the heat transfer characteristics of heated rotating cylinders subjected to varying rotational speeds. Becker [2] experimentally measured the heat transfer of a horizontal cylinder rotating in water. Ball and Farouk [3] experimentally investigated mixed convection around an isothermal rotating cylinder using a Schlieren method and recorded the deflection angle of the trailing plume under different rotating speeds. It was observed that when the Richardson number $γ=Gr/Re2$ reaches just below unity, the plume from the rotating cylinder becomes unstable and eventually breaks down. Jones et al. [4] experimentally investigated the combined forced and natural convection from a rotating horizontal heated cylinder in a low-speed wind tunnel. $O¨zerdem$ [5] experimentally investigated conductive heat transfer from a horizontal cylinder rotating in quiescent air and proposed the correlation $Nu=0.318Rer0.571$ for rotational Reynolds numbers ranging from 2000 to 40,000. Gschwendtner [6] investigated the heat transfer of a rotating cylinder in a crossflow using optical measuring techniques based on light deflection. Cheng et al. [7] investigated the convective heat transfer of a rotating heated cylinder and observed that the rotation-induced crossflow has a significant influence on the heat transfer coefficient on the cylinder surface. Ma et al. [8–11] conducted several experimental studies with heated rotating cylinders and observed critical rotational Reynolds numbers as the average Nusselt number varies with the rotational Reynolds number (Re_r). In another study by Ma et al. [12], rotation effects on the cylinder temperature and concentration boundary layer were investigated experimentally.

Apart from experimental observations, an abundant amount of research has also been conducted using computational fluid dynamics (CFD) to investigate convective flow and heat transfer over a rotating cylinder in free streams or confined spaces. Badr and Dennis [13] numerically investigated the laminar forced convective heat transfer from an isothermal rotating cylinder placed in a uniform stream, where the direction of the forced flow is normal to the cylinder axis. Makkonen [14] proposed a boundary-layer model for the heat transfer from the front half of a rough cylinder based on the integral equation of the boundary layer. Abu-Hijleh and Heilen [15] calculated entropy generation due to laminar mixed convection flow from an isothermal rotating cylinder. Yan and Zu [16] simulated a viscous liquid flowing past an isothermal rotating cylinder using the lattice Boltzmann method and reported the impact of the peripheral-to-translating speed ratio, Reynolds number, and Prandtl number on the flow and heat transfer characteristics. Pramane and Sharma [17] studied the two-dimensional (2D) freestream flow and forced convective heat transfer across a rotating cylinder, where vorticity dynamics involves vortex shedding at a critical rotational velocity. Sharma and Dhiman [18] investigated the effects of the Prandtl number on the heat transfer characteristics of an unconfined rotating cylinder. Bouakkaz et al. [19] studied the 2D flow pattern, time-averaged lift and drag coefficients, and Nusselt number around a heated rotating cylinder by constructing a 2D model using a finite volume-based commercial solver. Liao and Lin [20] investigated natural and mixed convection around a heated rotating cylinder placed within a square enclosure using an immersed boundary method. Elghnam [21] performed experimental investigations and 2D CFD simulations on a horizontal cylinder with various rotational speeds. Selimefendigil and $O¨ztop$ [22] investigated the mixed convection phenomenon of a rotating cylinder near a backward step and immersed in a nanofluid. Salimipour and Anbarsooz [23] revealed how the surface temperature of a rotating cylinder affects a compressible flow passing around a cylinder. Sasmal and Chhabra [24] studied laminar natural convection heat transfer in a power-law fluid from an isothermal rotating cylinder placed coaxially in a square duct. Bhagat and Ranjan [25] performed a computational study on the flow and heat transfer characteristics of a rotating cylinder subjected to a crossflow. Hassanzadeh et al. [26] investigated natural convective heat transfer around a rough rotating cylinder inside a square cavity.

Some studies were proposed to investigate the application of rotating cylinders or objects for thermal management. Fatla et al. [27] used CFD to investigate gas circulation generated by rotating cylinders inside high-temperature coil annealing furnaces during the annealing treatment of grain-oriented electrical steel. Loksupapaiboon and Suvanjumrat [28] experimentally and computationally investigated the forced convective heat transfer around a rotating hand-shaped former to improve rubber glove curing processes. Teamah and Hamed [29] conducted a numerical and experimental study of a multiphase flow inside a self-contained drum motor drive system (SCDMDS) which is used in food and pharmaceutical industries due to its contamination-free operations. Abd Al-Hasan et al. [30] used CFD to investigate heat convective flow in a vessel-tube array with a rotating baffle that represents the energy exchange in nuclear and chemical reactors.

Previous experimental studies [8–12] reported the heat convection characteristics along the radial and tangential direction of a cylinder. The question of how the space-averaged Nusselt number of a rotating cylinder is influenced by heat and fluid flow along the axial direction of the cylinder is little known. This knowledge is certainly critical as it is highly correlated with the effective dimensions in design of industrial applications. Prompted by this consideration, the present study aims to use a 3D CFD model based on the experimental setup of Ma et al. [8–12] to reveal the heat and flow patterns caused by a heated horizontal rotating cylinder.

The computational model in this study is constructed following the experimental setup in the study of Ma. et al. [10]. As shown in Fig. 1, a 3D cylinder with a length of 0.9 m and a diameter of 0.5 m is placed in a cylindrical air domain that contains a length of 4.9 m and a diameter of 10 m. The diameter of computational domain is 20-time larger than that of the cylinder to ensure a sufficient development of the flow and convection zone [31]. The rotational axis is parallel with the z-axis and rotates counterclockwise, and the origin of the Cartesian coordinate system is placed at an edge the rotating cylinder. The no-slip and no-penetration boundary condition is applied to the entire cylinder surface. The isothermal boundary conditions are imposed to the circumferential cylinder surface while the adiabatic boundary condition is set for both lateral walls of the rotating cylinder. The pressure outlet boundary condition is given to the circumference of the air domain using a gauge pressure of 0 Pa (absolute 101,325 Pa) to match ambient conditions and prevent unintended pressure gradients. With a low Mach number (Ma = 8 × 10^–3), the flow is considered incompressible. Since only pressure gradients matter, the choice of gauge pressure over absolute pressure does not affect the results. The symmetry boundary condition is applied for both lateral sides of the air domain. The surface temperature of the rotating cylinder and the temperature of the pressure outlet are maintained at 307.15 K and 293.15 K, respectively. Fluid motion in the current CFD setup is driven by the buoyancy force and rotating cylinder surface without imposed inlets.

Under our simulation conditions, the maximum rotational speed is 103.22 RPM, corresponding to a Mach number of 8 × 10⁻³, at which the assumption of incompressible flow remains valid [32]. Therefore, the CFD model describing the physical model is based on the following assumptions: (i) thermal and fluid properties of the working fluid are all constant except that the viscosity value is based on the kinetic theory of gases, and the density value is set with the air modeled as “incompressible ideal gas.” (ii) the radiation heat transfer between the cylinder and air, according to the experimental study, is negligible compared with the values of convection heat transfer [8].

We employ ansysfluent v2023-R1 for solving the dimensional governing equations of laminar heat convection [33] of air in a Cartesian coordinate system.

where $V$⁠, $T$⁠, g, and $p$ are the velocity vector [V = (u, v, w)], temperature, gravitational acceleration and pressure of air, respectively; $ρ$⁠, $t$⁠, $I$⁠, and $∇$ are density, time, unit tensor, and gradient operator [$∇ =(∂/∂x$⁠,$∂/∂y$⁠,$∂/∂z)$], respectively. The parameters $μ$⁠, $cp$⁠, $κ$⁠, $pop$⁠, $R$⁠, and $MW$ are the gas dynamic viscosity, specific heat capacity, thermal conductivity, operating pressure, gas constant and molecular weight, respectively. In this study, the operating pressure in the ideal gas state equation is assumed to be constant in accordance with the low Mach number condition applied in our simulation, as previously demonstrated by Chenoweth and Paolucci [34]. The operating pressure is taken to be the standard atmospheric pressure at sea level and 0 °C, with a value of 1 atm.

Our preliminary assessment indicates that convective flow solutions above the critical rotational Reynolds number obtained from the transition SST model [36] and shear stress transport (SST) k–ω turbulence model show no significant differences from each other. For saving the computational costs, the present study uses the SST k–ω model that combines the k–ω model in the near-cylinder region and the far-field calculations of the standard k-ε model without the inclusion of the transport equation for the intermittency and transition momentum thickness Reynolds number.

where κ_t is the turbulent thermal conductivity.

where the values of the model constants are: β* = 0.09, β = 0.075, σ_k = 0.5, σ_ω = 0.5, σ_ω,2 = 0.856, a₁ = 5.0, and α* = 0.555.

where G_k represents the production term of turbulent kinetic energy.

The first-order upwind scheme is applied to discretize the equations of turbulent kinetic energy and specific dissipation rate. A second-order upwind scheme is used for discretizing the 3D transient conservation equations of momentum and energy equations with the least squares cell-based gradient method. The time derivatives are discretized by the first-order implicit scheme. The discretized equations are solved by means of an implicit finite volume scheme based on the iterative SIMPLE algorithm employed for pressure–velocity coupling.

where $ν$, β,$h$⁠, and $r′$ are the gas kinematic viscosity, thermal expansion coefficient, convective heat transfer coefficient and dimensionless radial coordinate (⁠$r′$= r/D), respectively.

The temperature variation ratio, $ΔT/Tref$⁠, where $ΔT$ is $(Tcs−Ta)$ and $Tref$ is $0.5(Tcs+Ta)$⁠, shows the value of 0.0467 obtained from the fixed cylinder surface temperature at 307.15 K and ambient air temperature at 293.15 K. Since the present $ΔT/Tref$ is not extremely small, it is more reasonable to use an incompressible ideal gas rather than the Boussinesq approximation to account for density variations with temperature.

To investigate the effect of rotational speed on the heat convection, we fix Pr = 0.74 and Gr = 2.32 × 10⁸ while varying the rotational Reynolds number from 2 × 10³ to 8.5 × 10⁴. For Re_r = 2 × 10³–1.69 × 10⁴, corresponding to RPM of 2.43 to 20.57 in the experiment [10], the laminar flow equations are used to simulate the convective flow. For Re_r = 2.3 × 10⁴–8.5 × 10⁴, corresponding to the RPM of 27.93–103.22, the turbulent flow equations are employed for the fluid flow simulation. Although the cylinder is made of steel sheet covered with a thin layer of chromium [10], the turbulent flow model considers roughness on the cylinder surface due to minor eccentricity of a rotating cylinder. We perform the sensitivity tests of roughness constants (C_s) on the dimensionless temperature along the boundary layer radial coordinate at Re_r = 8.5 × 10⁴ (103.22 RPM) as shown in Fig. 2. Three roughness constants adopted for relatively uniform surface show that the temperature profiles are not much altered. In this study, we choose C_s = 0.5 without roughness height that properly describes the effect of eccentricity on the heat convection along a rotating cylinder wall [33]. In the present study, the convergence criterion of scaled residuals for continuity, momentum, energy, k, and ω equations are set to 10⁻⁵, 10⁻⁵, 10⁻⁶, 10⁻³, and 10⁻³, respectively. Further reduction of the scaled residuals is observed to yield no significant changes in the solutions. Figure 3 shows the sensitivity test for the time-step independence using errors of temperature profiles along the boundary-layer radial coordinate in the turbulent flow regime relative to the results with a dimensionless time-step size of Δt = 4.296 × 10⁻⁸. The selected Δt = 4.296 × 10⁻⁵, 4.296 × 10⁻⁶, 4.296 × 10⁻⁷, and 4.296 × 10⁻⁸ correspond to 0.5, 0.05, 0.005 and 0.0005 s in dimensional time, respectively. The comparison indicates that solutions computed with Δt = 4.296 × 10⁻⁷ and Δt = 4.296 × 10⁻⁸ have relatively small errors in temperature profiles. Accordingly, for both laminar and RANS models, the dimensionless time-step size of 4.296 × 10⁻⁸ (0.0005 s) is applied throughout the study. Each case is conducted over varying periods to ensure that the steady-state or periodically steady-state Nu numbers are achieved. Table 2 lists the RPM and corresponding Re_r numbers investigated in this study. In the experiments of Ma et al. [10], for conditions matching ours in terms of Prandtl and Grashof numbers, the reported Re_r ranges from 5 × 10³ to 5.1 × 10⁴. In this study, the rotational speeds selected for this investigation are to encompass both the laminar and turbulent regimes, while ensuring the consistency with the experimental conditions of Ma et al. [10] for further comparisons. In particular, a rotational speed of 20.57 RPM (Re_r = 1.69 × 10³) identified in the present simulation agrees well with their finding.

where t_p is the last 25 s in each numerical simulation which corresponds to the time taken for one revolution of the rotating cylinder operated at 2.43 RPM, the lowest one in the present study. Accordingly, 25 s is sufficiently long to capture the steady-state or periodically steady-state Nusselt numbers for each RPM investigated. L in Eq. (8) is the axial length of the cylinder.

where a is the order of the numerical method, f₁ and f₂ are the solutions for the fine and coarse grids, respectively. N₁ and N₂ refer to the number of elements on the fine and coarse grids, respectively. The safety factor f_s is chosen to be three, as recommended by Roache [37,38].

We perform the numerical simulation to reproduce the existing experimental data including $Nu¯¯cs$ on the cylinder surface and trailing vortex deflection angles. Figure 5 shows that the computed $Nu¯¯cs$ obtained at a wide range of the rotational Reynolds number generally agree with the experimental data and empirical correlations [10]. From Re_r = $2.0×$10³ to 1.69 × 10⁴, $Nu¯¯cs$ values obtained from the solutions of the laminar flow equations overlap completely the data estimated with the empirical correlations but are slightly higher than the experimental $Nu¯¯cs$ values. It is seen that the laminar flow model very much underpredicts $Nu¯¯cs$ value obtained at Re_r = 2.6 × 10⁴, implying that the laminar flow modeling is valid up to the critical Re_r (∼1.7 × 10⁴). This outcome confirms that turbulent flow equations need to be used at the rotational speed above critical Re_r. In the forced convection-dominated region, where the rotational speed is above the critical Re_r, the computed $Nu¯¯cs$ values overall are proportional to the rotational speed but with fluctuation. The data of empirical correlations lie in between the fluctuation of the computed $Nu¯¯cs$ values within the forced convection-dominated region. The computed $Nu¯¯cs$ values are close to the experimental measurements between Re_r of 2.3 × 10⁴ and 3.45 × 10⁴ but deviate from the experimental results at the higher rotational speed. In the region beyond the critical Reynolds number (∼1.7 × 10³), the present simulation demonstrates that the convective flow above the critical Re_r needs to be treated with turbulent flow modeling. The normalized root-mean-square errors (NRMSE) of the Nusselt numbers obtained from the present simulations versus experimental measurements and empirical relations are shown in Table 3. We evaluate the accuracy between RANS and large eddy simulation (LES) in terms of $Nu¯¯cs$ estimations. The SST k–ω model employed in our simulations achieves a remarkably low NRMSE (0.035) compared with the LES model, ensuring good agreement with LES solutions while significantly reducing computational cost. A discrepancy is observed between the predicted and measured $Nu¯¯cs$ values by Ma et al. [10], with an NRMSE of 0.156. However, the predicted $Nucs$ values exhibit strong agreement with those obtained from empirical models [10], as indicated by an NRMSE of less than 0.1. The NRMSE at 0 deg and 180 deg highlights the sensitivity of $Nu¯cs,ϕ$ in regions with pronounced turbulence and thermal effects, where discrepancies may stem from experimental uncertainties and measurement limitations. With the local Nusselt number difference normalized by $Nu¯¯cs$⁠, NRMSE of 0.095 further reflects reasonable overall consistency between the computation and measurements.

The comparison between the computed and measured Nusselt number at varied Re_r is further extended to the local $Nu¯cs,0o,180o$⁠, as shown in Fig. 6. It is seen that the computed $Nu¯cs,180o$ profile well captures the experimental data at Re_r < 3 × 10⁴ [9]. Figure 6 indicates that computed $Nu¯cs,180o$ at the descending side exhibits a local minimum time-averaged $Nu¯¯cs$ at the critical Reynolds number, which is consistent with the experimental measurements. The simulations that underpredict the computed $Nu¯cs,0o$ and $Nu¯cs,180o$ at relatively high Re_r are likely due to the RANS model that has insufficient resolutions to capture eddies correlated with heat convection. However, the difference between $Nu¯cs,0°$ and $Nu¯cs,180o$ that reveals how the heat transfer mechanism changes from natural to forced convection is in accordance with the experimental observation as shown in Fig. 7. The computed profile of (⁠$Nu¯cs,0o−Nu¯cs,180o$⁠)/$Nu¯¯cs$ surging at Re_r lower than the critical Re_r number supports the experimental finding that the mixed convection exists before the cylinder rotation reaches the critical Re_r number.

The deflection angle defined by the angle between the trailing vortex centerline and cylinder perpendicular bisector [8] is analyzed in the present simulation. As seen in Fig. 8, the computed deflection angles varied with Re_r/Gr^0.5 show very well agreement with the experimental measurement and the empirical correlation. The computed plume tilting toward the descending side with increasing Reynolds number aligns well with the experimental measurements obtained for Richardson numbers up to 2.69 in the study by Ma et al. [8]. Figure 9 further shows the computed temperature fields averaged in the last 25 s of each rotational speed with the Schlieren graph data obtained from Ma et al. [8]. Figs. 9(a,1)–9(a,4) and 9(b,1)–9(b,4) illustrate that as Re_r increases, the computed plumes exhibit a qualitative agreement with the features of plume dynamics observed in Figs. 9(c,1)–9(c,4). The NRMSE of the computed deflection angles is within 4.6% when compared with the Schlieren images by Ma et al. [8], and less than 5.6% compared with the empirical equation reported in Ref. [8] (see Table 3).

We present the contours of velocity magnitude and temperature averaged in the last 25 s in each numerical simulation. Figures 10–12 show the contours captured at the center plane along the axial direction. The range of the dimensionless velocity, $V′$⁠, is scaled to be equivalent to the dimensional velocity of 0.02–0.24 m/s for investigating the effect of rotational speed on convective flows. Within the laminar flow regime, the plume formed at the low Re_r of 2000 is mainly driven by the buoyancy force that moves the gas flow upward from the top of the cylinder and the plume is accelerated by the approach flow moving inward to fill the space yielded by the lifted gas (Fig. 10(a,1)). The temperature field shown in Fig. 10(b1) indicates that the thermal boundary layer is thin and heat transfer is primarily dominated by the natural convection that occurs on the top of the cylinder. As Re_r in the laminar flow regime is increased to the critical Re_r (Fig. 10(a,2)), the influence of the cylinder rotation on the air velocity field and thermal plume is considerable. At the critical Re_r, the velocity magnitude at the center of the plume is weakened but the structure of the plume is expanded in the air flow disturbed by the cylinder rotation. Although $Nu¯¯cs$ is not altered by Re_r increased to 1.7 × 10⁴, the surface heat convected away to the surrounding at the critical Re_r migrates toward the lateral side of the cylinder (Fig. 10(b,2)).

At Re_r = 2.3 × 10⁴, where the convected flow is simulated by unsteady Reynolds-averaged Navier-Stokes (RANS), the fluid plume is detached from the separation point located at around the middle of the deflected angle (Fig. 11(a,1)) and the size of the thermal plume is greater than that in the laminar flow regime (Fig. 11(b,1)). When Re_r is increased to 5.6 × 10⁴, the flow separation point is retarded toward the descending side, we observe that the fluid plume is merged with the rotating cylinder while the thermal plume leads to thickened thermal boundary layer (Figs. 11(a,2) and 11(b,2)). As Re_r is further increased to 5.7 × 10⁴, the stagnation point of the fluid flow continues to shift downward on the descending side (Fig. 11(a,3). At this moment, the thermal boundary layer continues to increase (Fig. 11(b,3)). Figure 12 shows that continuously increasing Re_r from 5.8 × 10⁴ to 8.5 × 10⁴ may lead to the merge of the flow plume and cylinder surface (Figs. 12(a,1)–12(a,3)) that prohibits the heat convection by thickening the thermal boundary layer thickness (Figs. 12(b,1)–12(b,3)). The simulation indicates that the center-plane flow plume detaching from or attaching to the cylinder surface is not able to be classified in terms of the Re_r range. More importantly, the results analyzed herein imply that the heat transfer enhancement proportional to Re_r depends on the heat and fluid flow behaviors along the axial direction.

Since the convective flow varies also along the axial direction of the rotating cylinder, we present the yz-plane contours of velocity magnitudes and temperature captured at deflection angles of their corresponding Re_r, as shown in Figs. 13–15. In the laminar flow regime, the Rayleigh–Bénard convection is observed on the top of the cylinder surface along the axial direction. With the dominance of natural convection, the flow plumes are primarily formed by the buoyant force (Figs. 13(a,1) and 13(a,2)), and the thermal plumes appear to extend along the axial direction (Figs. 13(b,1) and 13(b,2)).

As the rotational speed is at Re_r = 2.3 × 10⁴ (Fig. 14(a,1)), high enough to drive gas flow to overcome natural convection, the thickness of the gas flow stream circumferentially circulating above the cylinder surface near the center is enlarged, and explains the separation of the flow plume from the cylinder observed in Fig. 11(a,1). When the flow changes from laminar to turbulence in the yz-plane, the surrounding flow direction is entirely altered and dominated by the cylinder rotation that drives gas flow moving upwards together with the plumes. When Re_r is increased to 5.6 × 10⁴, we observe that the thumb-like flow plume is merged with the cylinder surface (Fig. 14(a,2)). The corresponding thermal plume shown in Fig. 14(b,2) is formed around the center section of the cylinder and the thermal boundary layer thickness near the edge is reduced. With the RANS simulation, multiple small mushroom-like flow and thermal plumes formed at the bottom of the cylinder at Re_r = 5.7 × 10⁴ begin to present at the cylinder bottom along the yz plane (Figs. 14(a,3) and 14(b,3)). As Re_r continues to increase, Figs. 15(a,1)–15(a,3) shows that the cylinder-driven gas velocity magnitude around the cylinder becomes profound. By collating Figs. 14(b,2), 14(b3), and 15(b1)–15(b3), it is evident that the thermal boundary layer thickness decreases outward from the cylinder center along the axial direction, primarily due to the thermal plumes emerging at the center.

To elucidate results in Secs. 3.1 and 3.2, we analyze the time-averaged $Nu¯(ϕ,z)$ along the circumferential angle of the cylinder surface at three different axial locations, i.e., the center (⁠$z′$ = 0), the end of axial direction (⁠$z′$ = 0.88) and the midway which locates between the center and end of the cylinder (⁠$z′$ = 0.5). Figure 16 presents the circumferential distributions of $Nu¯(ϕ,z)$ in the laminar flow regime at Re_r = 2.0 × 10³ and 1.7 × 10⁴. Notably, at both Reynolds numbers, $Nu¯(ϕ,z)$ exhibits significantly lower values along the upper side of the cylinder compared with the bottom side. This trend strongly correlates with the presence of the primary thermal plume above the cylinder, as observed in Fig. 10, which obstructs heat convection. Furthermore, the azimuthal locations where $Nu¯(ϕ,z)$ appears to reach its minimum coincide with the regions of plume formation, reinforcing the impact of plume dynamics on heat transfer. At Re_r = 2.0 × 10³, the variation of $Nu¯(ϕ,z)$ across different axial positions is insignificant, reflecting the uniform spread of the thermal plume, as illustrated in Figs. 13(a1)–13(b1).

The natural convection-driven flow and the rotationally induced flow are aligned on the ascending side and become increasingly significant with rising Re_r. In contrast, on the descending side, the natural convection flow opposes the cylinder-driven flow, resulting in a lower $Nu¯(ϕ,z)$⁠. Figures 17–19 illustrate the circumferential distribution of $Nu¯(ϕ,z)$ in the turbulent flow regime at Re_r = 2.3 × 10⁴, 5.8 × 10⁴ and 8.5 × 10⁴. As shown in these figures, $Nu¯(ϕ,z)$ profiles along the axial direction at Re_r = 8.5 × 10⁴ show higher values compared with those at Re_r = 5.8 × 10⁴ and 2.3 × 10⁴, indicating intensified convective heat transfer. This occurrence is further reflected in the average Nusselt number, where $Nu¯¯cs$ at Re_r = 8.5 × 10⁴ is greater than that at Re_r = 5.8 × 10⁴, which in turn is greater than that at Re_r = 2.3 × 10⁴ (see Fig. 5). This progressive enhancement in heat transfer corresponds to the increased turbulent mixing with increasing rotational speed, as observed in Figs. 11(a1)–11(b1) (Re_r = 2.3 × 10⁴) and Figs. 12(a1)–12(b1) (Re_r = 5.8 × 10⁴), 12(a,3)–12(b,3) (Re_r = 8.5 × 10⁴). In this turbulent regime, the cylinder-induced gas flow fully dominates heat convection, leading to enhanced heat advection.

The visualizations of plume formations along the axial direction of the cylinder at Re_r = 5.8 × 10⁴ and 8.5 × 10⁴, as observed in Figs. 15(a1)–15(b1) and 15(a3)–15(b3), reveal a thumb-like plume structure originating from the central region. These visualizations suggest that the influence of the plume on $Nu¯(ϕ,z)$ is most pronounced at the axial center for both Reynolds numbers. Consequently, $Nu¯(ϕ,z)$ at the center exhibits significantly lower values on the upper side, where the plume is present, compared with the bottom side of the cylinder. Moving away from the center, the difference between the upper and bottom sides of $Nu¯(ϕ,z)$ profiles becomes less pronounced at the midway and end location for Re_r of 5.8 × 10⁴ and 8.5 × 10⁴. This reduction in disparity indicates that the impact of the plume weakens at these locations compared with the central region.

Additionally, along the center of axial direction in the turbulent regime (Re_r = 2.3 × 10⁴, 5.8 × 10⁴ and 8.5 × 10⁴), $Nu¯(ϕ,z)$ profiles exhibit considerably lower values on the descending side, where the plume forms, which is consistently illustrated in Figs. 11(a1)–11(b1) and 12(a1)–12(b1), 12(a,3)–12(b,3) and suggests that plume formation hinders heat transfer in this region due to the increase in thermal boundary layer thickness. This observation implies a strong correlation between thermal boundary layer thickness and $Nu¯(ϕ,z)$⁠, as further evidenced in Fig. 20. In the laminar flow regime, at Re_r = 2.0 × 10³, the thermal boundary layer on the descending side is observed to be thinner than that on the ascending side (Fig. 20(a)), which corresponds to the slightly higher $Nu¯(ϕ,z)$ found on the descending side (see Fig. 16(a)). At this Re_r, the plume locates on the cylinder top, thus exhibiting less profound influence on the thermal boundary layer thickness difference between the two sides compared with that in turbulent flow regime. As Re_r increases to 2.3 × 10⁴ in the turbulent flow regime, a larger disparity in the thickness of the thermal boundary layers on the two sides of the cylinder can be observed (Fig. 20(b)), which aligns well with the observation that the difference in $Nu¯(ϕ,z)$between the two sides becomes more pronounced (Fig. 17). At$$Re_r = 8.5 × 10⁴$,$ this thickness difference becomes highly significant (Fig. 20(c)), further reinforcing relatively higher $Nu¯(ϕ,z)$on the ascending side compared with that on the descending side. As shown in Figs. 20(b) and 20(c), the thermal boundary layer at Re_r = 8.5 × 10⁴ is considerably thicker than that at Re_r = 2.3 × 10⁴. While the thicker boundary layer at Re_r = 8.5 × 10⁴ implies more limited heat transfer compared with that at Re_r = 2.3 × 10⁴, the increase in Re_r also enhances turbulent mixing where momentum exchanges are intensified, thereby improving heat transfer. Consequently, the overall heat transfer is significantly increased at Re_r = 8.5 × 10⁴, as consistently illustrated by the much higher $Nu¯¯cs$ at Re_r= 8.5 × 10⁴ compared with that at Re_r = 2.3 × 10⁴ (Fig. 5). As observed from the cylinder center toward its end section, heat transfer appears to be enhanced, as evidenced by the overall increase in $Nu¯(ϕ,z)$ at the midway and end sections of the cylinder. This observation aligns with the gradual reduction in thermal boundary layer thickness toward the cylinder ends at Re_r = 2.3 × 10⁴, 5.8 × 10⁴ and 8.5 × 10⁴, as depicted in Figs. 14(a,1)–14(b,1) and 15(a,1)–15(b,1), 15(a,3)–15(b,3), respectively.

This study presents a 3D CFD model for heat convection around an isothermal rotating cylinder of 307.15 K to complement experimental measurements of Nusselt numbers, deflection angles and Schlieren data. The model is investigated with fixed Pr = 0.74 and Gr = 2.32 × 10⁸, while varying Re_r from 2 × 10³ to 8.5 × 10⁴. The results successfully reproduce experimental data both quantitatively and qualitatively, providing a clear distinction between laminar and turbulent flow regimes.

Computed results reveal how cylinder rotation enhances heat transfer to the surrounding air. $Nu¯¯cs$ predictions confirm that laminar flow assumptions hold only up to the critical rotational Reynolds number, beyond which the SST k–ω turbulence model is required. In the laminar flow regime, the computed $Nu¯cs,180o$ profiles on the descending side align with experimental observations of inferior heat transfer performance compared with ascending side. However, at higher Reynolds numbers, the unsteady RANS model underpredicts heat transfer, likely due to limited near-wall resolution.

The plumes behaviors are governed by buoyancy forces in the laminar flow regime and cylinder-driven flow in the turbulent flow regime, as revealed by convective flow contours in both the circumferential and axial planes. Predicted $Nu¯(ϕ,z)$ profiles indicate that heat transfer enhancement in the turbulent regime is primarily driven by thermal boundary layer thinning at the cylinder edge. This 3D CFD study clarifies the heat transfer characteristics of a heated rotating cylinder, complementing previous experimental findings. Furthermore, the observed axial heat transfer variations highlight the necessity of a 3D model for accurately capturing convective heat transfer and fluid dynamics.

These findings provide insights into heat transfer in rotating machinery. In rotary kilns, plume detachment and boundary layer thickening affect efficiency and stability. In journal bearings, thickened boundary layers impact lubrication, risking failure. In heat exchangers and reactors, controlling boundary layer development optimizes performance. Future work will refine predictive models to capture 3D instabilities in flow transition and assess surface roughness effects under varying Grashof numbers. These advances will enhance heat transfer understanding and improve rotating thermal system design.

• National Science and Technology Council (NSTC) in Taiwan (Contract No. 110-2221-E-007-062-MY3; Funder ID: 10.13039/501100004663).

c_p =

gas specific heat capacity at constant pressure (Jkg⁻¹K⁻¹)

C_s =

roughness constant

D =

diameter of the cylinder (m)

f_s =

safety factor

Gr =

Grashof Number

h =

convective heat transfer coefficient (Wm⁻²K⁻¹)

I =

unit tensor

k =

kinetic turbulent energy (m²s⁻²)

Ma =

Mach number

n =

rotational speed of the cylinder (rpm)

$Nucs$ =

Nusselt number at the cylinder surface

$Nu¯¯cs$ =

time-space averaged Nusselt number

$Nu¯cs,z$ =

time- and local circumferentially averaged Nusselt number

$Nu¯cs,ϕ$ =

time- and local axially averaged Nusselt number

Pr =

Prandtl number

Pr_t =

turbulent Prandtl number

r =

radial coordinate (m)

$r′$ =

dimensionless radial coordinate, $r′=$r/D

R_c =

radius of the cylinder (m)

Re_r =

rotational Reynolds number

Re_r,cri =

critical rotational Reynolds number

t =

time (s)

$t′$ =

dimensionless time, $t′=$αt/D²

$T′$ =

dimensionless temperature, $T′$$= (T−Ta)/(Tcs−Ta)$

T_a =

temperature of ambient air (K)

T_cs =

temperature of the cylinder surface (K)

T_ref= =

qualitative temperature, T_ref = (T_a+ T_cs)/2 (K)

V = (u, v, w) =

velocity vector (ms^-1)

$V′=(u′, v′,w′ )$ =

dimensionless velocity vector, $V′=$VD/α

X = (x, y, z) =

Cartesian coordinate (m)

$X′=(x′, y′,z′ )$ =

dimensionless Cartesian coordinate, $X′=X/D$

ΔT =

difference between T_a and T_cs (K)

α =

gas thermal diffusivity (m²s^–1)

β =

gas thermal expansion coefficient (K^–1)

γ =

Richardson number

κ =

thermal conductivity of air (Wm⁻¹K⁻¹)

κ_t =

turbulent thermal conductivity of air (Wm⁻¹K⁻¹)

ν =

gas kinematic viscosity (m²s^–1)

ϕ =

angle along circumferential direction (deg)

φ =

general mean variable

ω =

specific dissipation rate (s^–1)

cri =

critical

cs =

cylinder surface