Numerical Simulation of Turbulent Flow in Eccentric Co-Rotating Heat Transfer

: Heat transfer engineering is signiﬁcant in many applications, especially in buoyancy natural convection in concentric and eccentric cavities. The biggest practical challenges, in this context, are capturing the self-natural ﬂow, estimating the mixing performance, and determining what parameters affect the temperature distribution in the cavity. In this paper, we focus on the improvement of a mathematical model, in order to enhance the accuracy of the solution, by investigating a new source term in the SST k − ω turbulence model based on the ﬁnite volume technique. The commercial numerical simulation software ANSYS Fluent 2021R1 is implemented to validate the accuracy. A concentric cavity was chosen for validation, the obtained temperature proﬁles at θ = 0 ◦ , θ = 30 ◦ , θ = 60 ◦ , θ = 90 ◦ , θ = 120 ◦ , θ = 150 ◦ , and θ = 180 ◦ were compared with previous experimental data. We applied this model to four eccentric rotating scenarios, including inner counterclockwise rotation, outer counterclockwise rotation, inner–outer clockwise rotation, and inner clockwise–outer counterclockwise rotation. The numerical simulation results reveal that the new source term in the momentum equation can produce superior results in the concentric test-case. The proposed mathematical model can describe the heat transfer under the eccentric co-rotation scenario well. Furthermore, the results for eccentric cases conﬁrm that the rotational direction affects the mixing temperature by generating a large vortex in the cavity, which increases the temperature mixing performance.


Introduction
When considering the motion of air, the fluid movement is caused by differences in fluid density called natural circulation. Natural convection has attracted a great deal of attention from researchers, due to its significance in nature and in engineering applications. In nature, convection cells formed by air rising above sun-warmed land or water lead to significant weather systems [1][2][3]. In engineering applications, for example, convection is commonly visualized in terms of the formation of microstructures during the cooling of nuclear power plants (i.e., when heat dissipation occurs in the main atomic reactor) [4,5]. A widespread industrial application of natural convection is free-air cooling without the aid of fans; this can happen on small-scale (e.g., computer chips) [6,7] to large-scale process equipment (e.g., natural mixed with forced circulation in a huge nuclear cooling tower) [8][9][10].
The understanding and control of heat transfer and fluid flow are essential for the success of the process. To this end, many researchers have focused on the behaviors of rotating cylinders. Lioa and Lin [11] investigated heat transfer under the conditions of different rotating speeds between a concentric circular cylinder and a square enclosure. The effect of rational velocity between two cylinders, where the inner cylinder is rotating and the outer cylinder is stationary, has been considered by Abou-Ziyan et al. [12]. Many researchers have focused on the ratio between concentric and eccentric tubes. Bagheri and Wang [13] have investigated the effect of the radius ratio on turbulent heat transfer in a concentric annular pipe flow. The results show that the radius ratio significantly impacts the Nusselt numbers and friction coefficients of the inner and outer cylinder walls. Abed et al. [14] have simulated flow behaviors and convective heat transfers in the gaps of concentric and eccentric cylinders. The rotation of the inner cylinder has a significant impact on the flow behavior. Alnakeeb et al. [15] conducted a numerical study of the eccentricity of the inner flat tube and improved the performance of the malting process. The results were compared with those of previous research, in order to confirm that the ratio of the inner flat tube has a significant effect on the melting process. Tang et al. [16] established a simplified model to predict the pressure gradient for laminar flow. Numerical simulation showed that the model was in agreement with the experimental data and the pressure gradient decreased with the increasing inner pipe rotation speed.
Solving a heat transfer problem using classical analytical methods is difficult if the situation involves complicated boundary conditions [17]. A recent development in the numerical simulation of heat transfer in eccentric and concentric pipes has led to physically interesting and mathematically challenging problems related to numerical solutions. Hassan et al. [18] developed a new finite volume-immersed boundary method to simulate fluid flow and heat transfer when considering complex geometries. Goodarzi et al. [19] used the finite volume method to simulate internal natural convection heat transfer inside cavities and enclosures. Over the past decade, many finite volume method-related techniques for solving the turbulence model have been proposed and investigated. Svacek et al. [20] investigated numerical solutions for incompressible turbulent flow over the backward facing step with various inclinations, based on the finite volume and finite element methods. Their results showed that the two methods were in good agreement, compared with previous data. Asadi et al. [21] investigated the effects of the heat transfer characteristics of a double-pipe heat exchanger equipped with turbulence-inducing elements. Numerical simulations were analyzed by ANSYS Fluent based on the finite volume method. They concluded that the double-pipe heat exchanger with the spherical elements had the best thermal performance.
The remainder of this paper is structured as follows. Explanations of the geometry, mesh details, and mathematical model are presented in Section 2. Then, the temperaturemixing performance of the proposed method is evaluated in Section 3. The results and conclusions are presented in Sections 4 and 5, respectively.

Methodology
In this study, we aim to evaluate a mathematical model implemented considering a non-isothermal turbulence system. The critical primary variable is the temperature distribution in the test-case cavity. Additionally, we include the other relevant and feasible behaviors for future research applications and other industries, such as rotation and corotation. A numerical simulation technique based on the finite volume method is applied in the data analysis. To better understand the flow locations in the cavity gap, with different areas between the small gap of the cavity, we chose seven locations to monitor the temperature. Furthermore, rotating effects were included, in order to investigate the temperature mixing performance.
The geometry of the eccentric annular cylinder was modeled following the standard academic test-case. Then, the problem description of the annular horizontal cylinders were calculated using the ANSYS Fluent 2021R1 software, as shown in Figure 1. A hexahedral structure grid was applied to obtain the best numerical accuracy. The mathematical modeling, grid-independence analysis, and boundary conditions are discussed in the following section.

Geometry and Mesh Details
A three-dimensional eccentric cylindrical tube was modeled in the numerical investigation following the previous experiment. The computational domain, specific boundary conditions, and evaluation data with different angles are shown in Figure 2. We emphasize that we studied the temperature variation between the cavity gap and the mixing performance during wall rotation under different conditions. Initially, the inner wall temperature (T i ) was set to 300 K and the outer wall temperature (T o ) was equal to 328 K, without a rotating velocity. The radii of the outer (R o ) and inner (R i ) cylinders were 46.3 mm and 17.8 mm, respectively. In the eccentric annulus case, the eccentricity was λ = 0.6245. This is close to the value of 0.623 reported in the experiment in [22,23]. The eccentricity is the measure of the distance that the inner cylinder is moved from the concentric position [24,25], defined as The mathematical examination of the spatial convergence of a simulation is a straightforward method for determining the ordered discretization error in a numerical simulation. This procedure involves performing the simulation on progressively finer grids, which is commonly known as grid-independence analysis. In order to obtain the most accurate velocity profile in the near-wall region, computational equations and intensive mesh with high quality are required to be applied on the wall during the simulation. We successively assessed four cases with different grid numbers [12,14,26], as shown in Figure 3. By con-ducting a grid-independence evaluation based on a critical parameter and monitoring the temperature profiles at θ = 60 degrees along with the radial distance (R), we suspected that it presented complex flow characteristics. The examination results suggested that case 4 (which represents 27,669,660 nodes and 2,668,800 elements) had consistent grid convergence, independent of the solution.

Mathematical Model
The difference of density in a fluid is the key driving mechanism. In fluid dynamics, the Boussinesq approximation is used in the field of buoyancy-driven flow, but is inaccurate when considering the dimensionless density difference. As previously mentioned, there has been no research on the efficiency of new source terms for natural convection heat transfer based on finite volume discretization. As a result, the development of mathematical models with new source terms has remained an open question. To solve this problem, we investigate how the mathematical process acts when applied to the simulation of non-concentric co-rotations in certain stages.
We aim to propose a new source term for the numerical investigation of the heat transfer of a turbulent natural convection heat transfer in an eccentric channel. The flow field is simulated using a Reynolds-averaged Navier-Stoke equation (RANS) with a two-equation eddy-viscosity model called the SST k − ω turbulence model; the naturally buoyant force is taken into account [27].

Turbulent Flow of an Incompressible Ideal Fluid
In terms of the mathematical equation, the SST k − ω model includes of all the refinements of the BSL k − ω model and accounts for the turbulence shear stress transport. The BSL k − ω model is similar to the STD k − model, but includes the following refinements [28,29]: • The STD k − ω model and the transformed k − model are multiplied by a blending function and added together. The blending function is designed to take a value of one in the near-wall region, which activates the STD k − ω model, and zero away from the surface, which activates the transformed k − model. • A damped cross-diffusion derivative term is included in the BSL model equation.

•
The model includes a new source term in order to capture the natural convection turbulent heat transfer. • The modeling constants are different.
These features make the SST k − ω model more accurate and reliable for a more comprehensive class of flows than the standard and BSL k − ω models. The conservation equation system is written as follows [30,31]: On the right-hand side of Equation (4), representing the production of turbulent kinetic energy, the dissipation of turbulent kinetic energy and the last term account for buoyancy due to turbulent natural convection flow, which has the model relation σ k , which is transformed by the first blending function (F 1 ) and β * . These are defined as follows: β * = 0.09 0.2667 + (ρk/8µω) 4 1 + (ρk/8µω) 4 (3) where ρ is fluid density, ρ 0 is reference fluid density, U is velocity, p is pressure, µ t is turbulent viscosity, k is turbulent kinetic energy, ω is specific dissipation rate, µ is dynamics viscosity, E is total enthalpy, T re f is reference temperature, y is the distance to the next surface, and D + ω is the positive portion of the cross-diffusion term. On the right-hand side of Equation (5), representing the generation of ω and the dissipation of ω due to turbulence, the cross-diffusion term with the blending function and the buoyancy terms are related by σ ω , α, α * , and β i , which are also transformed by the first blending function, and are defined by Additionally, the model relation parameter defined in Equation (7) has the sub-model constant, which needs to be blended between the inner layer and outer layer using the first blending function, and is defined by The mathematical system described above combines the advantages of the Wilcox and k − models, but still fails to properly predict the onset and amount of flow separation from smooth surfaces. The main reason is that these models do not account for the transport of the turbulent shear stress. This results in an over-prediction of the eddy-viscosity. A limiter can obtain the proper transport behavior for the formulation of the eddy-viscosity, which is modified as shown in Equation (11). This is a function of the strain rate magnitude, which is also used for the second blending function (F 2 ).
where S is the modulus of the mean strain rate tensor and S ij is strain rate tensor, defined as Equation (6) represents the energy equation, where k e f f is the effective conductivity (k + k t ), and k t is the turbulent thermal conductivity, which is defined according to the fluctuation of turbulent flow behavior. In this study, when heat is added to a fluid and the fluid density varies with temperature, a flow may be induced due to gravity acting on the density variations. Such buoyancy-driven flows are termed natural convection or mixed convection flows, which can be modeled by the new source term in Equation (3).

Initial and Boundary Conditions
A heated cylinder is placed inside another cylinder, trapping nitrogen in the resulting annular cavity. The inner cylinder is set in two configurations: one in which the cylinders are concentric for the purposes of the validation study and another in which the inner cylinder is displaced downwards. As the inner cylinder is hotter than the outer, as shown in Figure 2, a buoyancy-induced flow results, and turbulent natural convection occurs.
The assumptions, boundary conditions, and governing equations were employed to solve the problem on the studied computational domain. The commercial software, ANSYS Fluent 2021R1, based on the finite volume method (FVM), was adopted to solve the gas-flow turbulent natural convection heat transfer with and without a rotating inner/outer cylinder. The numerical simulation settings are detailed in Table 1, and all of the temperature results for each angle were compared with those in previously published papers [22,23,32].

Validation Study
Before we investigate the effect of natural convection under the rotating wall, in the next section, we discuss how we need to ensure that the SST k − ω with additional buoyancy source terms provides precise predictions that are in excellent agreement with experimental measurements [22,23]. The numerical simulation was performed using the commercial software ANSYS Fluent 2021R2, based on the finite volume technique. Simultaneously, the obtained temperature profiles at θ = 0 • , θ = 30 • , θ = 60 • , θ = 90 • , θ = 120 • , θ = 150 • , and θ = 180 • were compared with previous experimental data. Figure 4a expresses the contour of temperature distribution due to density variation, resulting in a significant mixing temperature between the tubes. Figure 4b-h compare the temperatures at each angle. The results suggest that the proposed mathematical model can capture the temperature in good agreement with experimental data. In the following section, we extend the model to evaluate eccentric cases with buoyancy turbulence-driven flow.

Results
In this section, we discuss the numerical results after performing the validation study and confirming that the SST k − ω turbulence model with the new buoyancy source in the momentum equation is in good agreement with the experimental test-case. Hence, the next numerical experiment included the rotating wall effect; in order to evaluate and investigate the flow behavior, it was divided into four cases, as discussed below.

The Effect of the Inner Counterclockwise Rotating Velocity
We aim to extend the capability of the proposed model to enhance the basic knowledge regarding the fundamental heat transfer, by including rotating cylinder walls in four scenarios, as shown in Figure 5. First, we start by assessing the effect of the inner counterclockwise rotating velocity, for which it is suggested that, when the inner rotating wall is between 2 and 10 rad/s, the high-temperature tail appears on the top of the inner wall but changes position, following the direction of rotation. As a result, it induced more mixing of temperatures between the top and bottom side regions at low rotation velocity, as shown in Figures 5a-c. Additionally, when the rotation velocity increased to 15-25 rad/s, the high-temperature tail broke out, causing more mixing between low and high temperature areas and leading to a lower temperature. If our system requires a low temperature, we determined that the rotation velocity must be more than 15 rad/s to achieve maximum mixing performance. On the contrary, if our system requires a high temperature, we must control the rotation velocity to be less than 15 rad/s. Finally, we can conclude that the critical rotation velocity in this study, based on forced convection, is 15 rad/s.

The Effect of Outer Counterclockwise Rotating Velocity
This section focuses on the outer wall rotating counterclockwise with velocity between 2-25 rad/s, as shown in Figure 6. The high-temperature tail remains at the top of the inner wall when the outer wall's rotation velocity does not meet 5 rad/s; however, when the rotation speed increases to more than 10 rad/s, the high-temperature tail is destroyed, which is suspected to be due to the fluid momentum induced by the outer wall. In addition, Figures 6d-f clearly show that the rotation of the outer wall produces a flow that changes its direction until a separation flow is formed on the right side of the inner wall. From the numerical results, it can be concluded that the external wall's rotation can improve the heat dissipation efficiency, although it is less than the internal counter-clockwise rotating wall. In practical engineering design, once the high-temperature area moves from the left-hand side to the right-hand side, it will likely affect the environment outside the eccentric cavity.
For safety, we should evaluate and consider whether any flammable material is located outside on the right-hand side and re-design or conduct prevention measures based on the findings.

The Effect of Inner-Outer Counterclockwise Rotating Velocity
The temperature distribution under the effect of inner-outer clockwise rotation is shown in Figure 7. We set the rotation velocity to 2 rad/s and observed that the hightemperature tail appeared at the top of the inner cylinder. When we increased the rotation velocity, the high-temperature tail began to be destroyed at 5 rad/s. Interestingly, the increasing counterclockwise rotation velocity shown in Figures 7c-f led to faster mixing than the observation test-case for the outer wall rotation. From this point of view, our research confirms that adding counterclockwise rotation into the system can improve the heat transfer performance and facilitate the development of more effective engineering applications in future equipment design in order to achieve high-accuracy eco-efficiency.

The Effect of Inner Clockwise, Outer Counterclockwise Velocity
In this section, we discuss the effect of co-rotation in the counterclockwise direction of the inner-outer eccentric cylinder wall. The numerical simulation results suggest that the co-rotation of the counterclockwise inner-outer wall affects the temperature mixing and temperature-distribution plume. The full results are shown in Figure 8. In the first scenario, we set both cylinder walls to rotate at 2 rad/s. The high-temperature tail followed the outer rotation direction. Meanwhile, when we increased the rotation velocity to 5 rad/s, the high-temperature tail started to curve following the direction of the outer rotation.
Next, the high-temperature tail disappeared in the flow system when the rotation velocity increased to 10 rad/s. Finally, a rotation velocity of 15 rad/s changed the direction and shape of the high-temperature tail in the inner right-hand side wall region, which decreased dramatically until the rotating velocity reached 25 rad/s. From these findings, we can conclude that the co-rotation of the inner and outer walls can increase the temperaturemixing performance by destroying the high-temperature tail, according to the momentum transfer from the wall motion to the fluid motion, which led to a decrease in mixing temperature in the computational domain.

Temperature Mixing Performance
Next, we discuss the mixing performance when hot and cold temperatures are considered between the eccentric cylinder walls. Figure 9 describes the static temperature after approaching the steady-state under a variety of wall rotation configurations.
First, when increasing the wall velocity, no matter the configuration, the mixing temperature homogeneously increased. However, the mixing temperature decreased when the rotational wall speed was higher than 10 rad/s. We know that, with this speed, the high-temperature tail will be destroyed due to the instability flow behavior shown in Figures 5-8. Thus, from this finding, we can conclude that the temperature tail depends on the mixing temperature and is affected by the wall rotation speed. Hence, in terms of engineering designs or applications, we can obtain two ideas from this finding. The first idea, which is applicable when using high mixing temperatures, is that we must rotate only the outer wall in the same direction or in the opposite direction; however, if we require a lower mixing temperature, we can rotate the inner wall more than 10 rad/s.

Velocity Distribution
The previous section considered the mixing performance due to the rotational speed of the walls. The conclusion was that an increase in rotational speed, where either configuration has the same behavior, led to a decrease in the range of the final temperature. In this section, we discuss how the flow structure is affected by the rotational wall and temperature-mixing characteristics when using 5 rad/s and 25 rad/s, in consideration of the above. In Figure 10a, the inner wall rotates under the right-hand side rule with 5 rad/s, creating two vortexes that are not symmetrical. The rotation induces a larger vortex on the right side than on the left side, compared with Figure 10b, where the inner wall rotates with 25 rad/s and induces a large single vortex between the eccentric cylinder. This is in agreement with the fact that the mixing performance of large vortexes leads to a higher capacity to increase heat transfer between the inner and outer walls. Next, Figures 10c,d consider only the outer wall rotating between 5 rad/s and 25 rad/s. The flow structure reveals that, with a low rotational velocity speed, there are multiple small vortexes on the left-hand side of the system, but, once the vortex is larger on the left, it creates an additional small vortex on the right side of the eccentric cavity. This behavior causes the final temperature to be lower than in the beginning, but still higher than in the inner rotating scenario. Furthermore, the case of inner-outer rotation with the same velocity is depicted in Figure 10e, which represents 5 rad/s, and Figure 10f, which represents 25 rad/s. This configuration interpretation confirms that a lower rotation velocity does not lead to more momentum transfer from the wall to the fluid, which generates the circulation vortex inside the cavity but can only produce a vortex on the top of the left side. This is the same as in other cases, and the temperature-mixing performance was the same as in the previous case. Finally, we consider the counterclockwise configuration, and the results of the vortex structure clearly differed from other cases: a medium vortex region was produced on the right side and a long shape vortex on the bottom left side. When we increase the speed, this flow structure converged to form only one large vortex on the top right side of the cavity, and the mixing performance was the same as that for the other cases, except for the inner wall rotation case. From the analysis, our findings, in terms of the flow structure based on the vortex circulation area, were in compliance with the flow physics discussed above. Moreover, we ensured that the new source term included in the momentum equation led to superior performance for investigating the turbulence buoyancy flow between concentric and eccentric cylinders.

Conclusions
In this paper, we focused on the numerical simulation of buoyancy natural heat transfer by including a new source term in the SST k − ω turbulence model. In this context, the new source term significantly improved the modeling of the flow inside the concentric cavity, extending the capability of the model to an eccentric rotating cavity. Numerical solutions based on the finite volume method considering the temperature distribution, mixing temperature, and vortex region described the situation well when the results were compared to previous validation test-cases. The results indicated that the inner-outer rotation velocity and the direction of the rotation of the inner and outer walls affected the mixing temperature. The vortex region that occurs in the eccentric cavity depends on the speed of the rotating wall. This new scientific knowledge can provide elaboration for engineering designs that require higher accuracy, such as those related to oil drilling, nuclear power cooling, and aerospace structures. Future research will consider the transient mix between natural convection and radiation heat transfer under the variation of parameters.