We conduct non-isothermal large eddy simulations (LESs) of flow past a heated cylinder (Re = 3900) to investigate flow physics throughout the wake region and develop a foundation upon which future heat flux wall models can be built (both for wall-modeled LES and other lower fidelity models) for mathematical closure of the energy equation. A rigorous validation of the mesh is made under isothermal conditions with results showing a closer match to experimental data than any other LES studies to date. The insights gained into the mesh design and approach are discussed. Simulation of non-isothermal flow is performed on the validated mesh for temperature differences between the cylinder surface and the freestream of 25 K and 300 K. The mesh design and realistic (temperature-dependent) thermodynamic property variations play key roles in predicting delayed separation, larger re-circulation zones, and enhanced turbulence intensity for the higher temperature difference case. The effect of both temperature differences on the flow is analyzed, and a new scaling of the flow domain is proposed to gain further insight into non-isothermal flow physics. Key scaling variables, friction temperature and friction velocity, are able to reduce nearly all of the temperature dependence of first and second order flow statistics, including turbulent heat fluxes. This leads to the finding that the turbulent heat flux in the wake region scales with the wall heat flux irrespective of the temperature difference in the flow.

## I. INTRODUCTION

The canonical flow of a cylinder in crossflow has been the topic of thousands of scientific inquires, both experimental and computational. The collective benefits from these studies cannot be understated, providing trusted correlations for fluid forces and pressure fluctuations^{1} as well as quantitative measurements of higher order metrics,^{2} thereby enabling robust closure to the turbulent momentum equations. Efforts such as these continue even today with a drive toward higher Reynolds number experimental data^{3} and simulations employing both wall-modeled^{4} and high-fidelity wall-resolved large eddy simulations (LESs).^{5} Similar impacts can be identified for flow past a heated cylinder, greatly increasing trust in design across a wide range of applications (e.g., industrial heat exchangers, pin fin heat sinks, and turbine blade cooling to name a few). However, most of these efforts have focused on average heat transfer or at most the angular dependence of time-averaged Nusselt numbers (e.g., the earlier works of Squire^{1} and Zukauskas^{6}). Higher order terms such as the turbulent heat flux must be provided either from experiments or direct numerical simulations (DNSs) in order to continue making advances in the vein of thermal wall models. The primary purpose of the current work is to provide rigorously validated large eddy simulations (LESs) of flow past a heated cylinder in order to form a solid foundation upon which thermal wall models can be based. To achieve this, we first consider an isothermal flow condition at Re = 3900 using multiple mesh sizes and compare to trusted data from experiments and DNS studies. This yielded LES results in tighter agreement with the underlying data than any other study to date, and findings from this effort (e.g., design of mesh) represent the secondary contribution of our work. We then present results from non-isothermal flow conditions with the cylinder at two different temperatures. There have been some recent efforts to investigate mixed^{7} and forced convection^{8} for flow past a heated cylinder at much lower Reynolds numbers in the range Re = 1–200. Non-isothermal simulations for similar flow configurations at higher Reynolds numbers (Re = 3000–30 000) have been limited to water jets^{9} and wall-bounded pin arrays.^{10} In all of these heat transfer studies, the primary focus has been on the heat transfer characteristics at the cylinder surface. Similarly, in the present study, detailed comparisons are made to the experimental data from the work of Nakamura and Igarashi^{11} where the angular dependence of the Nusselt number is reported in terms of average and root-mean-square (rms) values. Furthermore, the effects of non-isothermal flow on the near wake statistics and heat transfer characteristics are analyzed and contrasted with the recent LES paper of Jogee *et al.*,^{12} which was benchmarked with the same experimental data. The confidence we place in our results is rooted in isothermal and non-isothermal agreement with trusted data. This enables us to formulate physics-based reasoning for the differences in LES results, consider additional higher order metrics of importance, and draw conclusions of the dimensionless form of those metrics that are equally applicable for varying temperature boundary conditions. One of the long term aims of this research is to develop thermal wall models^{13} to be used in high Reynolds number flows by employing wall-modeled LES (WMLES), thus greatly reducing the computational cost while still providing high-fidelity results^{14} important in applications such as gas turbine blade cooling. While the current work is focused on much lower Reynolds numbers where the WMLES shows poor performance,^{15} the wall-resolved LES (WRLES) is used to analyze the flow with the aim of applying the concepts developed here to model heat transfer at conditions more relevant to gas turbine flows in future studies.

The first LES work on flow past a circular cylinder at a Reynolds number of 3900 was performed by Beaudan and Moin.^{16} Since then, there have been multiple efforts to improve the experimental and numerical results at this Reynolds number to better understand the flow dynamics in the wake of the cylinder. The more recent work of Lysenko *et al.*^{17} summarizes these endeavors and discusses the shortcomings of the numerical simulations when compared with the careful experiments performed by Parnaudeau *et al.*^{18} It is noted from these two studies that while the distribution of pressure, vorticity, and skin friction around the cylinder is well understood, there is considerable disparity between predictions of the flow features from numerical simulations and data from experiments in the wake region of the flow. Significant efforts to understand these differences have been undertaken by Parnaudeau *et al.*^{18} and Dong *et al.*^{19} The first of these focuses on resolving differences between LES results and experimental data obtained from Hot Wire Anemometry (HWA) and Particle Image Velocimetry (PIV) for Re = 3900, while the latter analyzes the differences between DNS simulations and PIV data for Re = 3900 as well as Re = 10 000. Both studies focus on the importance of capturing the inertial range of turbulence in the near wake region. This is an indicator of a non-dissipative numerical scheme that is required to capture the flow dynamics effectively. The forces on the cylinder (pressure, drag, and lift) and the near wake flow statistics depend greatly on this factor. Furthermore, it was shown that the main flow feature that governs the downstream velocity profiles and turbulence statistics is the re-circulation bubble length,^{18} i.e., the length of the re-circulation region of the mean flow behind the cylinder. Hence, it is important that the LES models used for the simulations capture these accurately. It should be noted that while LES is the model of choice for studying flows with separation and transitional shear layers, considerable research is being done to improve RANS (Reynolds averaged Navier Stokes) based models for such flows. Pereira *et al.*^{20} have performed a detailed analysis of a number of RANS models, PANS (partially averaged Navier Stokes) models, and a number of Detached Eddy Simulation (DES) models for studying the flow past a cylinder problem. It has been shown that Reynolds Stress Models (RSMs) and PANS models^{21} can predict flow separation and transition. However, the accuracy of those models is poorer than that of LES, especially for the re-circulation length prediction.

The work of Kravchenko and Moin^{22} and Lysenko *et al.*^{17} shows that dynamic (dyn.) sub-grid scale (SGS) models are required for LES to capture the inertial range of turbulence in the near wake region of the cylinder. Hence, the same is adopted for the present study. However, it is noted that the re-circulation length was over-predicted by Lysenko *et al.* and underpredicted by Kravchenko and Moin as compared to experimental data. Hence, although the near-wall quantities such as the coefficient of pressure (C_{P}), coefficient of drag (C_{D}), coefficient of lift (C_{L}), skin friction coefficient (C_{f}), and mean vorticity (Ω) lie within the range of experimental data, the mean velocities and Reynolds stresses downstream of the cylinder deviate significantly from experimental data. Lysenko *et al.*^{17} identified one of the main problems for these deviations as the mesh used by different researchers. Similarly, Kravchenko and Moin^{22} also showed that the grid design plays an important role in the prediction of the shear layer break-up process and hence the near wake flow characteristics that govern the re-circulation zone formation. A wide range of grid types and domain extents have been employed by researchers (see Table I) to ascertain the combination best suited for simulating flow past a cylinder. Although the most popular mesh type is O-type or C-type,^{17} an H-type domain following the work of Wissink and Rodi^{23} is used here since the long-term goal of this project is to move toward complicated geometries and hence unstructured (uns) grids. Recent research into unstructured grids for simulation of flow past a cylinder shows great promise. For example, a variational multiscale dynamic SGS model has been shown to yield very accurate predictions of the separated flow region,^{24} and adaptive mesh refinement techniques on unstructured grids can reduce computational effort significantly.^{25} Currently, however, a structured grid design is adopted for simulating the flow past a cylinder at Re = 3900, with unstructured grid efforts reserved for future work at higher Re.

Numerical . | Turbulence . | Mesh . | Domain extent . | Mesh size . | |||
---|---|---|---|---|---|---|---|

contributors . | model . | type . | L_{x}/D
. | L_{y}/D
. | L_{z}/D
. | N (M) . | N_{z}
. |

Ma et al.^{38} | DNS | H (uns) | −15 to 25 | 18 | 1–2 π | … | … |

Dong et al.^{19} | DNS | H (uns) | −15 to 25 | 18 | 1–2 π | … | … |

Wissink and Rodi^{23} | DNS | O–H | −10 to 15 | 20 | 4–8 | >10 | 512 |

Kravchenko and Moin^{22} | LES (B-splines) | O | 40 | 40 | π | 2.41 | 48 |

Franke and Frank^{40} | LES (comp. Smag.) | O–H | −10 to 20 | 20 | π | 1.14 | 33 |

Parnaudeau et al.^{18} | LES (IBM) | O | 20 | 20 | π | 9 | 48 |

Mani et al.^{41} | LES (comp. dyn. Smag.) | O | 35 | 35 | π | 5.50 | 48 |

Lysenko et al.^{17} | LES (dyn. k-eqn.) | O | 50 | 50 | π | 5.76 | 64 |

Jogee et al.^{12} | LES (comp. dyn. k-eqn.) | H | −5 to 20 | 20 | π | 2.06 | 38 |

Current | LES (dyn. k-eqn.) | O–H | −5 to 20 | 20 | π | 4.81 | 40 |

Numerical . | Turbulence . | Mesh . | Domain extent . | Mesh size . | |||
---|---|---|---|---|---|---|---|

contributors . | model . | type . | L_{x}/D
. | L_{y}/D
. | L_{z}/D
. | N (M) . | N_{z}
. |

Ma et al.^{38} | DNS | H (uns) | −15 to 25 | 18 | 1–2 π | … | … |

Dong et al.^{19} | DNS | H (uns) | −15 to 25 | 18 | 1–2 π | … | … |

Wissink and Rodi^{23} | DNS | O–H | −10 to 15 | 20 | 4–8 | >10 | 512 |

Kravchenko and Moin^{22} | LES (B-splines) | O | 40 | 40 | π | 2.41 | 48 |

Franke and Frank^{40} | LES (comp. Smag.) | O–H | −10 to 20 | 20 | π | 1.14 | 33 |

Parnaudeau et al.^{18} | LES (IBM) | O | 20 | 20 | π | 9 | 48 |

Mani et al.^{41} | LES (comp. dyn. Smag.) | O | 35 | 35 | π | 5.50 | 48 |

Lysenko et al.^{17} | LES (dyn. k-eqn.) | O | 50 | 50 | π | 5.76 | 64 |

Jogee et al.^{12} | LES (comp. dyn. k-eqn.) | H | −5 to 20 | 20 | π | 2.06 | 38 |

Current | LES (dyn. k-eqn.) | O–H | −5 to 20 | 20 | π | 4.81 | 40 |

Significant efforts have been made by researchers to understand the differences between numerical predictions and experimental data in the near wake region. Earlier experiments were not able to provide reliable data both in the near wake and in the far-wake locations. The data from the experiments of Lourenco and Shih^{26} were more reliable in the near wake, while those of Ong and Wallace^{2} were more accurate in the far-wake. These experimental limitations were overcome by Parnaudeau *et al.*^{18} to provide extremely reliable experimental data for this canonical flow configuration. It was identified that the accuracy of prediction of the re-circulation bubble length from simulations was crucial for achieving good agreement with experimental data in the wake region. In fact, a rigorous numerical procedure was laid out by Parnaudeau *et al.* to obtain this value from simulations. Subsequent simulation efforts by Lysenko *et al.*,^{17} Jogee *et al.*,^{12} and others have shown reasonably good agreement with the data of Parnaudeau. The differences in the mean velocity and Reynolds stress predictions can be attributed to the difference in re-circulation bubble length predictions. Attempts have been made to resolve these discrepancies in the near wake in the present study. This is done in order to isolate the effects of the temperature variations and turbulent transport of heat of the non-isothermal flow from inaccuracies in the underlying (isothermal) velocity field. Furthermore, the lack of exhaustive experimental data in the near wake of such flows drives the need for the highest possible accuracy in the isothermal flow field simulation. This will enable proper analysis of the effects of heat transfer on the flow by isolating any numerical simulation errors. The only numerical works in this field are those by Kim and Nakamura,^{27} Salkhordeh *et al.*,^{28} and Jogee *et al.*^{12} The first two works only compare the Nusselt number profiles around the cylinder with the experimental data of Nakamura and Igarashi.^{11} To the best knowledge of the authors, the work of Jogee *et al.*^{12} is found to be the only simulation effort where the effect of the cylinder temperature on the near wake flow is studied in detail. However, it is observed that the re-circulation bubble length predicted by the isothermal LES of Jogee *et al.*^{12} is larger than the PIV data of Parnaudeau *et al.*^{18} and Dong *et al.*^{19} by about 11% and 14%, respectively. This causes inaccurate predictions of the isothermal flow features and significantly affects the non-isothermal flow results. Hence, careful mesh design to enhance the prediction of the re-circulation bubble length and hence the near wake flow statistics is of utmost importance for the present study.

The first part of this paper focuses on understanding how to improve the current state of LES results using better grid design using well-established SGS models. After carefully setting up the mesh to obtain good agreement with experimental data, the same mesh is then used to simulate non-isothermal flow past a cylinder with an elevated wall temperature boundary condition. Two non-isothermal cases are studied with the wall temperature at 25 K and 300 K above the freestream temperature corresponding to the smallest and largest temperature differences studied by Jogee *et al.*^{12} The second part of this paper is dedicated to understanding the effect of the temperature field on the flow and vice versa. As can be seen from the preliminary work by Jogee *et al.*,^{12} the re-circulation bubble length is altered due to the heat transfer that impacts the velocity and Reynolds stress profiles in the near wake region. This, in turn, changes the heat transfer characteristics and the mechanisms through which heat is removed from the cylinder during the vortex shedding process. It is noted that the turbulent heat flux increases significantly as the wall temperature of the cylinder is increased. Efforts have been made to better understand these differences at higher temperatures using scaling laws. Proper scaling of the relevant variables leads to the emergence of new physical insights, which help provide a framework to establish more rigorous wall models to represent the turbulent heat flux. One could then implement such models for higher Reynolds number simulations where the computational cost might preclude a higher fidelity solution such as wall-resolved LES.

## II. NUMERICAL METHODOLOGY

For the simulations to be performed as described above, the open-source Computational Fluid Dynamics (CFD) software OpenFOAM v19.12 was employed. The Navier–Stokes equations in their compressible form along with the energy equation were solved to consider the effects of the density change due to the non-isothermal flow field. The Mach number was kept very low (less than 0.05) for all the cases studied so that the effects of compressibility did not require additional modeling. Most of the models, discretization schemes, and the convergence criterion used by Lysenko *et al.*^{17} are followed since good accuracy of results has been demonstrated with these combinations. Second order schemes have been used for both space and time discretization although further research into developing more advanced discretization schemes for separated flows continues.^{29} Specifically, the *rhoPimpleFoam* solver in OFv19.12 has been employed with some modifications to incorporate the equations described hereafter, along with linear schemes for spatial discretization and backward Euler for time. Variable time stepping is employed with the maximum Courant number constrained to a value of 0.8. This results in the time step varying between 1.6 × 10^{−3} and 3.0 × 10^{−3} *tU*_{∞}/*D* over all meshes tested. As described before, the sensitivity of the mesh is studied in greater detail here while keeping the underlying models fixed. A detailed description of the grid design is provided after introducing the governing equations.

### A. Governing equations, models, and averaging procedures

First, the governing equations in their compressible, density-weighted (Favre-filtered) form are presented in Eqs. (1)–(4). These equations are derived from the continuity [Eq. (1)], momentum conservation [Eq. (2), Navier Stokes], energy conservation [Eq. (3)], and ideal gas equation of state [Eq. (4)] using the LES filtering operation,

Details of the filtering procedure are discussed later. As a first approximation, the CFD mesh serves as a filter for the continuous conservation equations to yield filtered equations, i.e., all scales of flow larger than the local grid scale are resolved by the equations, while the effect of the flow scales smaller than the grid scale is modeled using an appropriate SGS model. The resolved or filtered quantities are represented by over-bars and are directly solved for. Namely, they are the filtered pressure ($p\xaf$) and filtered density ($\rho \xaf$) along with the Favre-filtered or density-weighted velocity ($ui\u0303$ for *i* = 1, 2, 3) and temperature ($T\u0303$). The Favre-filtered temperature is defined as $T\u0303=\rho T\xaf/\rho \xaf$ and similarly for the three velocity components. The energy equation [Eq. (3)] is written in terms of the sensible enthalpy that is related to the temperature field through the isobaric specific heat as $hs\u0303=cPT\u0303$. The additional terms in these equations represent the flux terms that need to be modeled. While $\sigma ij\xaf$ and $qi\xaf$ represent the molecular or viscous momentum and heat fluxes, respectively, $\tau ij\xaf$ and $\theta i\xaf$ represent the effect of the SGS fluctuations of the flow on the filtered quantities, namely, the SGS momentum stress and the SGS heat flux. The first two terms are closed using Newton’s law of viscosity [Eq. (5)] and Fourier’s law of heat conduction [Eq. (6)],

where the filtered strain-rate tensor is given by $Sij\u0303=12\u2202ui\u0303\u2202xj+\u2202uj\u0303\u2202xi$ and

It is noted here that three additional terms requiring SGS modeling are neglected presently. These are the SGS pressure dilatation $DpDt\xaf\u2212Dp\xafDt$ and the SGS viscous dissipation $\sigma ij\u2202uj\u2202xi\xaf\u2212\sigma ij\xaf\u2202uj\u0303\u2202xi$, which appears in the energy equation [Eq. (3)], and the SGS interactions between fluctuating thermal conductivity and temperature gradients $\mu Pr\u2202hS\u2202xi\u0303\u2212\mu \xafPr\u2202hS\u0303\u2202xi$, which appears in the expression for conductive heat flux in Eq. (6). Garnier^{30} discussed recent efforts that show that SGS viscous diffusion is an order of magnitude smaller than the SGS heat flux. Given the low Mach number regime of the current flow, the pressure dilatation is also negligible. In the present study, the effect of the SGS fluctuation in thermophysical properties has also been neglected since high-fidelity LES has been performed, which is expected to resolve most of the fluctuations in temperature and other thermodynamic properties of the fluid.

These equations along with the relationship between enthalpy and temperature contain three thermodynamic variables: the specific heat capacity at constant pressure $cP\xaf$, dynamic viscosity $\mu \xaf$, and thermal heat conductivity $k\xaf$, which can be combined into a single non-dimensional number called the Prandtl number, $Pr=\mu \xafcP\xaf/k\xaf=\nu \xaf/\alpha \xaf$. This represents the ratio of the fluid kinematic viscosity ($\nu \xaf$) to the thermal diffusivity ($\alpha \xaf$). For the present study, the working fluid is air with a constant Prandtl number of 0.7071. The dynamic viscosity and specific heat are considered as functions of temperature only, which subsequently determine the value of the thermal conductivity. The SGS momentum stress tensor is given as $\tau ij\xaf=\rho \xafuiuj\u0303\u2212ui\u0303uj\u0303$. This requires modeling since it involves the non-linear interaction of the velocity field, which is not captured as a part of the filtering process of LES. The eddy-viscosity hypothesis [Eq. (7)] adopted by Smagorinsky^{31} is extensively used for SGS modeling since it follows a procedure similar to Newton’s law of viscosity and, hence, is easy to implement in CFD schemes. According to this hypothesis, the sub-grid stress can be expressed as

The main difference between the molecular stress and the sub-grid stress is that while the coefficient of the molecular strain-rate tensor (*μ*) is a property of the fluid, for the sub-grid stress, the coefficient is a property of the flow and hence requires a turbulence closure model. A wide range of SGS models for LES exist depending on how the term $\mu SGS\xaf$ is modeled. Presently, the dynamic k-equation model similar to that adopted by Lysenko *et al.*^{17} is used since that has been shown to perform well for the case of flow past a cylinder. The details of this modeling approach were developed by Kim and Menon.^{32} They used the solution to an additional transport equation for the sub-grid scale kinetic energy, $kSGS=12\tau kk$, as developed by Yoshizawa,^{33} given by

and to close the SGS eddy-viscosity, the following expression is employed:

The constants *C*_{ε} and *C*_{k} are determined as described by Kim and Menon.^{32} The turbulent heat flux term $\theta i\xaf$ is also modeled using a similar eddy-dissipation hypothesis as given by

where a constant value of sub-grid Prandtl number *Pr*_{SGS} = 0.85 is used.^{34}

This sub-grid model considers non-local and history effects through the direct computation of the subgrid-scale kinetic energy that removes local equilibrium (i.e., equilibrium between SGS energy production and dissipation rate) approximations. This is done by considering a test filter (typically twice the size of the grid-scale filter) such that the sub test-filter scale stresses can also be modeled using an equation similar to Eq. (7) relating the sub test-filter scale stress to the test-filtered strain rate tensor. The approximation that the same $\mu SGS\xaf$ expression can be used for these two cases leads to an equation for *C*_{k}, which is evaluated using a least squares method. Similarly, equivalent expressions for the dissipation rate (see Ref. 32 for equations) for the grid-scale filter and the test-filter are used to evaluate *C*_{ε}. These dynamical evaluations along with the solution of Eq. (8) make this model dynamically adaptive to local flow variations. Hence, it is ideal for use in flows involving transition to turbulence.

Since the flow past a cylinder is highly transient, care needs to be taken in averaging the results when they are compared with experimental data and other numerical simulations. In the grid design, all meshes tested have the same number of equally spaced cells in the spanwise direction with cyclic boundary conditions in that direction. Hence, the fluxes leaving the front face of the domain re-enter through the back face and vice versa. All results presented henceforth are averaged in this homogeneous direction over all spanwise cells. Besides this spatial averaging, temporal averaging is also performed to analyze the results. The flow is initialized with a flat velocity profile with a magnitude *U*_{∞} at the inlet and throughout the entire domain. The value of the velocity is chosen such that *Re* = *ρ*_{∞}*U*_{∞}*D*/*μ*_{∞} = 3900, where *ρ*_{∞} = *P*_{∞}/*R*_{air}*T*_{∞} (from ideal-gas law) with air as the working fluid with *R*_{air} = 287.7 J/(kg K) at standard pressure (*P*_{∞} = 1 atm) and temperature (*T*_{∞} = 300 K) in the freestream and *D* = 0.01 m. The freestream dynamic viscosity is evaluated at the freestream temperature using Sutherland’s law,^{35} given by

where *A*_{S} = 1.460 553 × 10^{−6} kg/(ms K) and *T*_{S} = 110.4 K.

In order to remove the time history of the inlet conditions, the simulation is run for some time before statistics are collected to analyze the flow. Parnaudeau *et al.*^{18} suggested a duration of 150 time-units (*t* = 150 *D*/*U*_{∞}) to obtain a fully-established turbulent wake with no influence of the initial conditions. After these initial cycles, the simulation is re-started to gather time-averaged data. The vortex shedding cycle frequency (*f*_{vs}) is calculated from the fast Fourier transform (FFT) of the lift coefficient. This is then used to obtain statistics over 250 vortex shedding cycles (corresponding to 1097 *tU*_{∞}/*D*) as suggested by Parnaudeau *et al.*^{18} to obtain a converged re-circulation bubble length and at every tenth computational time step. Recent studies also suggest strong sensitivity of the time averaging window to reach statistical convergence for the re-circulating length.^{36} The results presented henceforth denoted as ⟨*X*⟩ represent the spanwise averaging of the temporal averaging (over 250 cycles) of the quantity *X*.

Two non-isothermal cases are studied with a fixed cylinder wall temperature of 325 K and 600 K, respectively. The first is chosen to benchmark the solution with the available Nusselt number data of Nakamura and Igarashi,^{11} and the second is chosen to test the most extreme case simulated by Jogee *et al.*^{12} While the same turbulence models are used for the non-isothermal flow studies, the temperature dependence of key thermodynamic quantities is achieved through trusted correlations. The dynamic viscosity is considered strictly as a function of temperature following Sutherland’s law, given by Eq. (11). The specific heat capacity is a polynomial function of temperature,^{37} given by

where *T* is the temperature of the fluid in Kelvin.

The temperature field is used to evaluate these quantities during each simulation run time. The constant Prandtl number value of *Pr* = 0.7071 is then used to determine the thermal conductivity of the fluid as *k*(*T*) = *μ*(*T*)*c*_{P}(*T*)/*Pr*.

Using these equations, it can be calculated that the thermodynamic fluid properties at the wall deviate significantly from the freestream properties. The dynamic viscosity, specific heat capacity, and thermal conductivity of the fluid change by 6.28%, 0.17%, and 6.07% at T = 325 K and by 63.4%, 4.54%, and 41.46% at T = 600 K, respectively. Since the viscosity and the conductivity show a significant variation from the inlet of the flow to the near-wall (boundary layer) flow, especially at the high temperature condition, it is expected that the flow and heat transfer characteristics will be significantly altered due to these changes in properties.

### B. Computational grid and boundary conditions

The equations described above are discretized onto a finite-volume grid. The mesh design parameters used in recent studies, shown in Table I, are used to fix the domain extent used in the present study. Researchers adopting unstructured grids usually employ an H-type grid for simplicity and generality. While different streamwise (*L*_{x}) and transverse (*L*_{y}) extents of the domain have been used by researchers, presently the domain used by Jogee *et al.*,^{12} shown in Fig. 1, has been invoked to make direct comparisons between the non-isothermal cases. They employed an H-type domain following the work of Wissink and Rodi.^{23} A conventional spanwise extent of *L*_{z} = *πD* is used, consistent with all LES studies. Although Ma *et al.*^{38} observed the sensitivity of DNS results to changing spanwise extent, it is noted that for LES, this effect is reduced and an extent of *πD* is sufficient to capture the large eddies in that direction as long as the spanwise cell size does not exceed 0.1D.^{39} The number of grid points used in that direction is 40, which is slightly larger than that used by Jogee *et al.*^{12} but smaller than other researchers. The choice of boundary conditions for this study is the same as that used by most researchers for a rectangular grid: cyclic boundary conditions in the spanwise direction, symmetric conditions in the transverse direction, and inlet and outlet (wave transmissive condition to prevent reflection of pressure waves) conditions at the inlet and outlet planes, respectively, in the streamwise direction.

Presently, an O–H type grid is overlaid on the H-type domain, as shown in Fig. 2. An O-grid is used close to the cylinder with n_{o} grid points and an expansion ratio (ER) denoted as ER (O-grid). The square region of dimension 2x_{n2} × 2y_{n2}, excluding the circular O-grid region, is the transition region in the mesh where the O-grid transitions into a rectangular H-grid used in the remaining part of the domain. The number of cells used here in the transition region is n_{t} with an expansion ratio ER (transit). Whereas the O-grid (analogous to prism layers at solid surfaces for unstructured grids) is essential for resolving the boundary layer, the transition region captures the separated shear layers. It will be shown how both of these are important to resolve all important flow characteristics. In the O-grid, the cells in the circumferential direction are equally spaced so as to remove issues with improper aspect ratios (ARs).^{23} The near-wall cell aspect ratio (AR), defined as the ratio of the wall-normal cell height to the wall-parallel (circumferential) cell width, plays an important role in capturing the near-wall flow features. At the present Reynolds number of 3900, the oncoming flow is laminar, resulting in a laminar boundary layer at the front stagnation point up to about **θ** ∼ 80°. However, the flow proceeds with a loss of pressure causing flow separation to occur. The separated flow causes shear layers in the wake of the cylinder, resulting in flow fluctuations due to shear layer instabilities.^{39} These fluctuations lead to turbulence and cause a re-circulation zone in the wake of the flow that further affects the near-wall flow development, separation, and transition to turbulence. Hence, the important features of flow past a circular cylinder can be divided into two main groups: the forces on the cylinder due to the flow and the turbulent flow profile behind the cylinder due to the obstruction in the flow. It is found that the cell sizes near the cylinder and the near-wall AR control the resolution of the near-cylinder flow including the boundary layer and the ER in the O-grid and the mesh-transition region control the resolution of the wake region of the flow. However, since the flow separation (which occurs at the near-wall) and the shear layer fluctuations in the wake are dependent on one another (see Ref. 19), there is mutual interdependence as well. In addition to the AR and the ER of the mesh, the cell size at the cylinder wall in the radial direction, denoted as *Δr*_{w}, is also varied in these meshes. *Δr*_{w} governs the friction Reynolds number at the first grid point in the cell center of the cell adjacent to the wall, *r*^{+}, where *r*^{+} = *Δr*_{w}*u*_{τ}/*ν*, $u\tau =\tau w/\rho $ is the friction velocity, *τ*_{w} is the wall shear stress, and *ρ* and *ν* are the fluid density and kinematic viscosity at the wall cell, respectively. For a wall-resolved LES of the flow, the *r*^{+} should remain within 1 to ensure that the boundary layer is well resolved. Although Parnaudeau *et al.* shows good wake region flow predictions even with the removal of this restriction, for the present study, *r*^{+} is maintained below unity since it is found to impact the Nusselt number {*Nu* = *Dq*_{w}/[*k*_{∞}(*T*_{w} − *T*_{∞})], where *q*_{w} is the wall conductive heat flux calculated using fluid properties at *T*_{w}} prediction significantly for the non-isothermal cases.^{28}

### C. Available experimental data and simulation results

Some of the most reliable and widely used experimental data for validation of simulation results are listed in Table II. While hot-wire anemometry (HWA) experimental techniques have mainly been used traditionally, recent particle image velocimetry (PIV) data, especially those of Parnaudeau *et al.*,^{18} provide much more robust validation data. Hence, most LES results from the present research have been compared to PIV data. It should be noted that new experimental techniques such as smoke image velocimetry (SIV) continue to be developed and used for studying fluid flows with separation and transition.^{42} The SIV data are in very close agreement with the PIV data from the work of Parnaudeau *et al.*^{18} The important parameters that govern the near-wall flow and the forces on the cylinder due to the flow past it are controlled by the location on the cylinder where the flow separates from the wall. This is characterized by the angle between the front stagnation point and the location on the cylinder where separation occurs, denoted as the angle of separation (**θ**_{sep}). It is determined as the location on the cylinder (measured from the front stagnation point) where the local mean shear stress goes to 0. The other characterizations include the base pressure coefficient (C_{pb}) that is the coefficient of pressure at the rear stagnation point, the root-mean-square (rms) of the lift coefficient (C_{L}), the mean drag coefficient (C_{D}), and the mean vorticity (Ω). Two parameters governing the transient flow behind the cylinder are commonly noted: the Strouhal number (St), given by *St* = *f*_{vs}*D*/*U*_{∞}, where *f*_{vs} is the vortex shedding frequency, *D* is the diameter of the cylinder, and *U*_{∞} is the freestream velocity at the inlet, and the re-circulation bubble length (L_{rc}) defined as the distance from the rear stagnation point to the point along the centerline where the mean streamwise velocity changes sign from negative to positive.

Experimentalists . | Expt. Tech. . | θ_{sep}
. | −C_{pb}
. | rms C_{L}
. | Mean C_{D}
. | St . | L_{rc}
. |
---|---|---|---|---|---|---|---|

Norberg^{43,44} | (Multiple) | … | 0.85–0.95 | 0.04-0.15 | 0.93–1.05 | 0.21 | 1.4–1.6 |

Cardell^{45} | HWA | … | 0.93 | … | 1.0 | 0.21–0.22 | 1.3–1.5 |

Lourenco and Shih^{26} | PIV | 86 | … | … | 0.97 | … | … |

Ong and Wallace^{2} | HWA | … | … | … | … | 0.21 | … |

Dong et al.^{19} | PIV | … | … | … | … | 0.212 | 1.47 |

Parnaudeau et al.^{18} | PIV | 88 | … | … | … | 0.208 | 1.51 |

Experimentalists . | Expt. Tech. . | θ_{sep}
. | −C_{pb}
. | rms C_{L}
. | Mean C_{D}
. | St . | L_{rc}
. |
---|---|---|---|---|---|---|---|

Norberg^{43,44} | (Multiple) | … | 0.85–0.95 | 0.04-0.15 | 0.93–1.05 | 0.21 | 1.4–1.6 |

Cardell^{45} | HWA | … | 0.93 | … | 1.0 | 0.21–0.22 | 1.3–1.5 |

Lourenco and Shih^{26} | PIV | 86 | … | … | 0.97 | … | … |

Ong and Wallace^{2} | HWA | … | … | … | … | 0.21 | … |

Dong et al.^{19} | PIV | … | … | … | … | 0.212 | 1.47 |

Parnaudeau et al.^{18} | PIV | 88 | … | … | … | 0.208 | 1.51 |

A detailed analysis of five meshes is performed to understand the effect of mesh design on these parameters. Table III shows the mesh design used by researchers (same list of researchers as in Table I) in the past and their corresponding results. Both lift and drag coefficients predicted by LES studies lie within experimental ranges noted in Table II. In addition, the base pressure coefficient and the Strouhal number predicted by LES agree well with the experimental data. The DNS performed by Ma *et al.*^{39} and Wissink and Rodi^{23} predict Strouhal numbers close to 0.22, which is at the higher end of the experimental data range. It is also seen that the higher Strouhal number prediction of Ma *et al.*^{38} compared to LES studies leads to predictions of a slightly lower back-pressure coefficient. However, the greatest discrepancy with experimental results is found in the over-prediction of the re-circulation length, as was noted before. The notable exceptions are the LES of Kravchenko and Moin^{22} and the DNS of Dong *et al.*^{19} who predict much smaller values. Their predicted re-circulation length matches closely with the data obtained by the experiments of Cardell^{45} but under-predicts more recent data of Dong *et al.*^{19} and Parnaudeau *et al.*^{18} It was identified by Kravchenko and Moin^{22} that mesh refinements near the cylinder wall and in the wake region were responsible for these changes in predictions. A similar mesh refinement strategy with the implementation of the expansion ratios (ERs) suggested by Mittal and Moin^{46} leads to considerable improvement in current simulations in the prediction of L_{rc} as compared to PIV data.^{18,19} The values of expansion ratios range from 1.03 to 1.08 for other LES studies; however, reducing these values to the range of 1.01–1.03 in current simulations gives much better prediction of the re-circulation length. In addition, the base pressure coefficient from current simulations falls well within experimental ranges, while the mean drag and the rms lift coefficient are higher. It was noted from the work of Ouvrard *et al.*^{47} that lower L_{rc} predictions lead to larger fluctuations in the lift coefficient. This explains the higher C_{L} rms value as compared to Lysenko’s results. A more rigorous analysis presented in Sec. III explains the reason for the high mean drag coefficient prediction from present simulations. It is also seen that the predicted Strouhal number is closer to the DNS results of Khan *et al.*^{39} and Wissink and Rodi.^{23} The increase in frequency and drag could be due to the shorter re-circulation zone as also noted by Kravchenko and Moin^{22} in comparison to other LES studies.

. | Wall cell height . | . | ER . | ER . | . | . | . | . | . | |
---|---|---|---|---|---|---|---|---|---|---|

Contributors . | Δr_{w}/D
. | Maximum r^{+}
. | AR . | At wall . | In wake . | −Cpb . | CL (rms) . | CD (mean) . | St . | Lrc/D . |

Ma et al.^{38} | … | … | … | … | … | 0.96 | … | 0.84 | 0.219 | 1.59 |

Dong et al.^{19} | … | 0.32 | 0.5 | … | … | 0.93 | … | … | 0.203 | 1.36 |

Wissink and Rodi^{23} | … | 0.66 | 0.38 | … | … | … | … | … | 0.216 | 1.59 |

Kravchenko and Moin^{22} | … | … | … | … | … | 0.94 | … | 1.04 | 0.21 | 1.35 |

Franke and Frank ^{40} | 0.0035 | 1 | … | 1.032 | 1.05 | 0.85 | … | 0.978 | 0.209 | 1.64 |

Parnaudeau^{18} | 0.021 | … | 1.0 | … | … | … | … | … | 0.208 | 1.56 |

Mani and Moin^{42} | 0.000 62 | … | … | 1.04 | … | 0.86 | … | 0.99 | 0.206 | … |

Lysenko et al.^{17} | … | … | … | 1.085 | … | 0.91 | 0.089 | 0.97 | 0.209 | 1.67 |

Jogee et al.^{12} | … | 1 | … | … | … | … | … | … | … | 1.68 |

Current | 0.002 | 0.8 | 0.25 | 1.012 | 1.03 | 0.91 | 0.121 | 1.048 | 0.221 | 1.48 |

. | Wall cell height . | . | ER . | ER . | . | . | . | . | . | |
---|---|---|---|---|---|---|---|---|---|---|

Contributors . | Δr_{w}/D
. | Maximum r^{+}
. | AR . | At wall . | In wake . | −Cpb . | CL (rms) . | CD (mean) . | St . | Lrc/D . |

Ma et al.^{38} | … | … | … | … | … | 0.96 | … | 0.84 | 0.219 | 1.59 |

Dong et al.^{19} | … | 0.32 | 0.5 | … | … | 0.93 | … | … | 0.203 | 1.36 |

Wissink and Rodi^{23} | … | 0.66 | 0.38 | … | … | … | … | … | 0.216 | 1.59 |

Kravchenko and Moin^{22} | … | … | … | … | … | 0.94 | … | 1.04 | 0.21 | 1.35 |

Franke and Frank ^{40} | 0.0035 | 1 | … | 1.032 | 1.05 | 0.85 | … | 0.978 | 0.209 | 1.64 |

Parnaudeau^{18} | 0.021 | … | 1.0 | … | … | … | … | … | 0.208 | 1.56 |

Mani and Moin^{42} | 0.000 62 | … | … | 1.04 | … | 0.86 | … | 0.99 | 0.206 | … |

Lysenko et al.^{17} | … | … | … | 1.085 | … | 0.91 | 0.089 | 0.97 | 0.209 | 1.67 |

Jogee et al.^{12} | … | 1 | … | … | … | … | … | … | … | 1.68 |

Current | 0.002 | 0.8 | 0.25 | 1.012 | 1.03 | 0.91 | 0.121 | 1.048 | 0.221 | 1.48 |

## III. MESH VALIDATION STUDIES

Five meshes were used to test the sensitivity of the results to different refinement parameters. The AR and ER for the five meshes are reported in Table IV consistent with the definitions previously given. All expansion ratios are kept lower than 3% as suggested by Mittal and Moin^{46} for the stability of LES schemes. Although this criterion was developed for high order (4 and 5) upwind schemes for O-grids, it is tested in the present study with second order central schemes following the suggestions of Kravchenko and Moin^{22} and Franke and Frank^{40} who noted high truncation errors of central schemes on coarse grids. As discussed by Franke and Frank, such errors lead to an over-prediction of the re-circulation length due to inadequate resolution of the shear layer. The 3% criterion was hence invoked to remove such over-predictions. In meshes I and II, the grid is highly refined near the wall so that the maximum value of *r*^{+} does not go beyond 0.4 and 0.5, respectively, throughout the simulation run time. This was done to test the capability of refined near-wall cells in capturing the boundary layer formation and separation with high accuracy as this becomes increasingly important for non-isothermal flows during separation.^{28} Next, the *Δr*_{w} is increased for meshes III–V to improve the AR. Different ERs in the O-grid and in the transition region of the mesh are used to test their sensitivity. Mesh V invokes smaller expansion ratios only in the wake region of the flow. For meshes with higher ARs, the ERs in the O-grid and transition regions of the mesh are reduced to ascertain whether better results can be obtained. The ER in the near wake (O-grid and transition regions) is kept the same between meshes II and III and then is progressively reduced for meshes IV and V. As is presented next, all of the tested meshes give quite good results when compared with experimental data, especially in the near wake region of the flow. They show significant improvement from previous LES results. This proves that the current mesh design adopted resolves the flow structures in the wake adequately. In general, it is found that the ER in the near-wall and transition region needs to be refined to about 1.01, i.e., the cell growth should be about 1% in the near wake region (up to 3 D downstream). The AR of the cells near the cylinder wall need to be at least 0.2 or larger in order to resolve the near-wall flow development, which is essential in capturing the force coefficients accurately. Further downstream, the ER of the cells in the streamwise direction is kept at about 3% following previous studies.

Mesh . | Total . | . | Maximum . | AR . | ER . | ER . | . |
---|---|---|---|---|---|---|---|

No. . | cells (M) . | Δr_{w}/D
. | r^{+}
. | (O-grid) . | (transit) . | (wake) . | ER . |

I | 3.17 | 0.001 | 0.4 | 0.10 | 1.02 | 1.024 | 1.028 |

II | 3.84 | 0.001 25 | 0.5 | 0.15 | 1.018 | 1.024 | 1.028 |

III | 3.77 | 0.002 | 0.8 | 0.20 | 1.018 | 1.024 | 1.03 |

IV | 4.02 | 0.002 | 0.8 | 0.25 | 1.015 | 1.015 | 1.03 |

V | 4.81 | 0.002 | 0.8 | 0.25 | 1.012 | 1.010 | 1.03 |

Mesh . | Total . | . | Maximum . | AR . | ER . | ER . | . |
---|---|---|---|---|---|---|---|

No. . | cells (M) . | Δr_{w}/D
. | r^{+}
. | (O-grid) . | (transit) . | (wake) . | ER . |

I | 3.17 | 0.001 | 0.4 | 0.10 | 1.02 | 1.024 | 1.028 |

II | 3.84 | 0.001 25 | 0.5 | 0.15 | 1.018 | 1.024 | 1.028 |

III | 3.77 | 0.002 | 0.8 | 0.20 | 1.018 | 1.024 | 1.03 |

IV | 4.02 | 0.002 | 0.8 | 0.25 | 1.015 | 1.015 | 1.03 |

V | 4.81 | 0.002 | 0.8 | 0.25 | 1.012 | 1.010 | 1.03 |

It should be noted that for all present studies, the O-grid has a diameter of 2D and the transition region is a square around the cylinder of length 6D. Usually, the near-wall grid refinement is maintained over a very small region around the cylinder (e.g., Lysenko had a grid refinement in the region of 1.006D^{17}). However, a larger O-grid region is maintained in present simulations to circumvent a particular limitation of the current O–H grid. The transitioning grid from the O-grid to the rectangular grid in the bulk flow is well behaved only along the primary axes (centerline *x* axis and transverse *y* axis); however, as one moves away from the axes, the discrepancy between the cell sizes across the interface of the O-grid and the transition-grid increases. This discrepancy yields a poor mesh quality, especially at 45° and 135° to the incoming flow. The cells in the transition region, at the interface, are 1.5 times larger in the radial direction as compared to the cells along the streamwise and transverse axes. In order to keep this mesh change away from the wall and the shear layer separation region, the O-grid has a larger diameter of 2D. The size of the transition region was chosen so as to encompass all the major flow features in the transverse direction up to the outlet face of the domain.

### A. Forces on the cylinder

The sensitivity of the current mesh design for predicting the pressure forces and shear forces on the cylinder is studied first. Figure 3 shows that all meshes can capture the mean skin friction coefficient [Fig. 3(a)] and the mean vorticity [Fig. 3(b)] quite well when compared to previous LES efforts of Breuer^{48} and experimental data by Son and Hanratty,^{49} respectively, thus providing confidence on the near-wall flow structure resolution. However, the mean pressure coefficient [Fig. 3(c)], when compared to the experimental data of Cardell,^{45} is not accurately reproduced with meshes I and II. The low AR of the cells near the cylinder walls in those meshes means that the near-wall flow is not adequately resolved. This leads to predictions of lower pressure in the separated region of the flow, although the separation point does not change significantly with different meshes. Meshes with low ARs near the wall depend more on the sub-grid model for the prediction of the streamwise velocities since the grids are stretched more in that direction. This effect is removed as the AR is improved with meshes III–V and more scales of the flow are resolved rather than being modeled, leading to better matching of the pressure data with experiments (see Table V for the values of the back-pressure coefficient, −C_{pb}, compared to the experimental value of 0.9 from Ref. 45). It is also noted that although meshes III and IV improve the C_{P} prediction, further refinement of the mesh in the wake region [lower ER (transit)] is required to get accurate matching of experimental data. This implies that the near wake flow provides feedback to the near-wall flow structures through the re-circulation bubble. Hence, mesh V performs the best in predicting the pressure coefficient in the separated region of the flow. Slight differences between the C_{P} profile from mesh V and the experimental data are observed at the leading edge and at the minimum value of C_{P} for all meshes tested. These differences can be attributed to the compressible solver, and comparable observations were made by Wornom *et al.*,^{50} who also employed a similar compressible solver.

Mesh . | Total . | . | % error from . | . | % error from . | rms . | Mean . | . | . |
---|---|---|---|---|---|---|---|---|---|

No. . | Cells (M) . | θ_{sep}
. | Breuer (θ_{sep} = 87)
. | −C_{pb}
. | Norberg (−C_{pb} = 0.9)
. | C_{L}
. | C_{D}
. | St . | L_{rc}
. |

I | 3.17 | 86.9 | 0.11 | 0.93 | 3.33 | 0.157 | 1.084 | 0.224 | 1.428 |

II | 3.84 | 87.5 | 0.57 | 0.97 | 7.78 | 0.128 | 1.065 | 0.226 | 1.509 |

III | 3.77 | 87.1 | 0.11 | 0.93 | 3.33 | 0.133 | 1.067 | 0.224 | 1.517 |

IV | 4.0 | 86.8 | 0.23 | 0.93 | 3.33 | 0.117 | 1.058 | 0.227 | 1.571 |

V | 4.8 | 86.5 | 0.57 | 0.91 | 1.11 | 0.127 | 1.063 | 0.228 | 1.480 |

Mesh . | Total . | . | % error from . | . | % error from . | rms . | Mean . | . | . |
---|---|---|---|---|---|---|---|---|---|

No. . | Cells (M) . | θ_{sep}
. | Breuer (θ_{sep} = 87)
. | −C_{pb}
. | Norberg (−C_{pb} = 0.9)
. | C_{L}
. | C_{D}
. | St . | L_{rc}
. |

I | 3.17 | 86.9 | 0.11 | 0.93 | 3.33 | 0.157 | 1.084 | 0.224 | 1.428 |

II | 3.84 | 87.5 | 0.57 | 0.97 | 7.78 | 0.128 | 1.065 | 0.226 | 1.509 |

III | 3.77 | 87.1 | 0.11 | 0.93 | 3.33 | 0.133 | 1.067 | 0.224 | 1.517 |

IV | 4.0 | 86.8 | 0.23 | 0.93 | 3.33 | 0.117 | 1.058 | 0.227 | 1.571 |

V | 4.8 | 86.5 | 0.57 | 0.91 | 1.11 | 0.127 | 1.063 | 0.228 | 1.480 |

Table V summarizes the predictions of near-wall characteristics and important parameters of the mean wake flow (Strouhal number and re-circulation length) from the five meshes tested. The near-wall flow is characterized by the separation angle (angle at which the mean skin friction coefficient ⟨*c*_{f}⟩ goes to 0 along the cylinder surface), the base pressure coefficient (C_{pb}), the rms of the lift coefficient, and the mean drag coefficient. As discussed before, with low AR, meshes I and II yield high values of −C_{pb}. In addition, the predicted mean drag and rms lift coefficients are larger than those suggested by experiments (see Table II). As the AR is improved keeping the ER in the near wake region constant from mesh II to mesh III, the back-pressure coefficient prediction significantly improves from −0.97 to −0.93 (experiments suggest a range of −0.85 to −0.95^{44,45}) for meshes III and IV. It is found that a further reduction in the ER in the near wake region is required to get improved predictions of the flow. As the C_{pb} improves for mesh IV, so do the over-predictions of mean drag, rms lift, and St. However, the re-circulation length shows an increased value. Hence, the ER in the transition region of the grid is reduced to test its sensitivity on the L_{rc} prediction. These changes significantly improve the near-wall predictions (especially the C_{P} prediction) and the flow predictions in the near wake region for the finest mesh, mesh V. This demonstrates that the flow separation at the cylinder wall and the shear layer instabilities in the wake of the cylinder are coupled through the re-circulation zone formation. Thus, it is imperative for the mesh to resolve both these flow phenomena with equal accuracy to yield accurate predictions. The successful resolution of this coupling using the current grid design is demonstrated through the improved prediction of the re-circulation length (while accurately capturing the near-wall flow parameters). The value of 1.48 is obtained from 250 vortex shedding cycles as suggested by Parnaudeau *et al.*^{18} The re-circulation length calculated from time-averaged and spanwise averaged streamwise velocity also varies with time until enough vortex shedding cycles have passed through to yield converged statistics. Consistent with that suggested by Parnaudeau *et al.*, 250 such cycles are required to achieve convergence in our present simulations. The re-circulation length obtained presently lies between the two values reported by current PIV experiments of Dong *et al.* (L_{rc} = 1.47)^{19} and Parnaudeau *et al.* (L_{rc} = 1.51).^{18} However, the Strouhal number, mean drag, and rms lift of the fine mesh (mesh V) are larger than other LES results but lie just within the experimental ranges of Norberg. Since Lysenko *et al.*^{17} provides all these numbers from their simulations, comparisons are made with their results. While the Strouhal number and the mean drag are about 5% and 8% larger than our present simulations, as compared to Lysenko’s LES results, the rms of the lift is predicted to be almost 36% larger. These differences are mainly due to the reduced re-circulation length prediction (11% smaller than the predicted value from Ref. 17). It is also noted that the low L_{rc} values also cause high Strouhal number predictions compared to other LES results that are more consistent with DNS results and slightly larger than those obtained from experiments. It has been shown by Ouvrard *et al.*^{47} that a shorter L_{rc} prediction corresponds to a higher rms C_{L}. This explains the higher rms of the lift coefficient compared to that of Lysenko *et al.* since all current meshes predict lower values of L_{rc}.

The reason for the increase in mean drag from current simulations is addressed next. The variation of C_{D} and C_{L} with time from current simulations on mesh V is shown in Fig. 4. It is noted that in the time intervals where the C_{L} is found to fluctuate more (e.g., *tU*_{∞}/*D* of 850–920 and 1100–1140), the mean C_{D} also shows a corresponding increase. Similar observations can be made from Fig. 2 of Lysenko *et al.*^{17} in the *tU*_{∞}/*D* intervals of 900–920 and 1160–1200. This is suggestive of the fact that these two coefficients are strongly connected to each other. In order to test this hypothesis, the entire time history of the coefficients is divided into equal intervals and the corresponding rms C_{L} and mean C_{D} are calculated over those intervals. Figure 5(a) shows the distribution of these for an averaging width interval of *tU*_{∞}/*D* = 10. It is found that these quantities are positively correlated in this interval. In order to assess the influence of the chosen interval, it is varied and the slope of the linear correlation curve is determined for each interval. The dependence of the correlation coefficient on the averaging width is shown in Fig. 5(b). It is found that the slope remains fairly constant up to about *tU*_{∞}/*D* = 20, suggesting a temporal correlation within that time period. Thus, it can be concluded that larger fluctuations in the lift coefficient prediction lead to increased mean drag predictions. Furthermore, the increase in mean drag prediction is about 10% of that of the increase in the rms of lift. Thus, following the argument of smaller L_{rc} leading to larger rms C_{L},^{48} this in turn leads to larger mean C_{D}, as is consistently observed in the higher values of mean C_{D} in all present meshes, as compared to previous numerical simulations. This also explains the high value of mean C_{D} observed by Kravchenko and Moin^{22} (Table III) in their simulation corresponding to their low L_{rc} value, when compared to other researchers.

### B. Flow statistics in the wake

In order to make comparisons of the downstream velocity and Reynolds stress profiles, it is required to test the scale-resolving power of the LES models employed. The power spectrum of the velocity field is a good indicator of this.^{18} The velocity field is probed at four downstream locations of x/D = 3, 5, 7, 10 along the centerline (y/D = 0) following the work of Parnaudeau *et al.*^{18} Figure 6 shows these spectra for the streamwise and transverse velocity components at the four locations for the fine mesh (mesh V). The spectra are calculated using MATLAB’s fast Fourier Transform algorithm using data from the simulation collected over 100 vortex shedding cycles (440 tU_{∞}/D). All grid points in the spanwise direction are treated as independent and used to average the spectra. It is found that the transverse velocity field shows two peaks at *f*_{vs} and 3*f*_{vs}, while the streamwise velocity shows a single peak at 2*f*_{vs}. These are consistent with the observations made by Parnaudeau *et al.*^{18} Furthermore, at almost all the locations, the spectra cover at least a decade of the inertial scaling that further gives confidence that the numerical simulations are not too dissipative and do a good job in resolving the inertial range and capturing the large eddies of the flow. Similar spectra-resolving profiles were observed for all the other meshes tested. Two other metrics important for re-circulation zone formation are flow separation and shear layer instability. Flow separation indicates the separation of the attached boundary layer from the cylinder surface. The predicted values of the mean separation angle from the different meshes are tabulated in Table V, and the available experimental data are tabulated in Table II. It is seen that values from all meshes lie within the experimental ranges. Hence, the near-wall LES resolution is sufficient to capture the effect of the adverse pressure gradient that causes flow separation. The shear layer instability is of importance because this leads to transition to turbulence in the wake and hence controls the formation of the re-circulation zone. Figure 7 shows the profile of instantaneous vorticity at two stages of the vortex shedding cycle: the separation of the shear layer from the cylinder at the top and bottom surfaces in Figs. 7(a) and 7(b), respectively. These continue to develop in the laminar form before the instabilities become too large and the shear layers start rolling up to form the primary vortex. The first image shows that the lower vortex has already been shed, thereby forming the wake in the rear of the cylinder. The top shear layer has rolled up into a primary vortex that is shed in the second image. This alternate vortex shedding process dictates the flow development in the rear of the cylinder, and accurate simulation of this is required to predict the re-circulation bubble length of the mean flow. The detailed dynamics of the shear layer instabilities, frequencies, and transition process leading to vortex shedding can be found in Ref. 39. Before making comparisons of the mean near wake flow field, the structure of the mean flow field is studied using streamlines in Fig. 8. The re-circulation bubble extends just beyond x/D = 2 for the simulation using mesh V. In the transverse direction, it extends to y/D = ±0.5 due to the two large re-circulation zones on either side of the centerline. This zone is the main region of turbulence production through shear in the mean velocity field.

In order to demonstrate the estimation of the re-circulation bubble length, the mean streamwise velocity profile along the centerline and the mean streamwise Reynolds stress variations are shown in Fig. 9. The re-circulation bubble length is defined as the distance from the rear stagnation point to the point where the mean streamwise velocity becomes zero (i.e., changes sign from negative to positive) along the centerline. It is seen that for all meshes tested here, the mean streamwise velocity prediction follows the experimental results of Parnaudeau *et al.*^{18} more accurately than the results of Lysenko *et al.*^{17} and Jogee *et al.*^{12} although they both use the same SGS model with a larger mesh (5.76M cells) by Lysenko *et al.* and a smaller one (2.06M cells) by Jogee *et al.* Hence, these differences can be attributed to the mesh design adopted here. Mesh I predicts the flow very close to the cylinder quite accurately but over-predicts the increase in velocity in the region 1.5 D–2.5 D. Hence, the re-circulation length is under-predicted. Meshes II and III show very similar predictions indicating that the AR and the ER in the far-wake region does not considerably alter the flow predictions within the re-circulation zone. Mesh IV over-predicts the re-circulation length by under-predicting the increase in velocity away from the cylinder. The fine mesh, mesh V, accurately predicts the increase in mean velocity toward the end of the re-circulation zone, thereby producing very close agreement with Parnaudeau’s experimental data. The mean streamwise Reynolds stress variation also shows good agreement with experimental data. It is seen that when the re-circulation zone length is over-predicted grossly (e.g., LES results of Jogee *et al.*^{12}), the peak in the Reynolds stress profile is also significantly larger than experimental predictions. This is because the re-circulation zone is a zone of high strain in the mean flow and hence of large Reynolds stresses. Thus, larger re-circulation zone lengths lead to larger turbulence production. Although very good agreement with the experimental Reynolds stress profile is observed with all meshes up to about x/D = 2.5, the Reynolds stress is under-predicted in the far-wake region, indicating that there is still room for improvement in LES models. Presently, the major focus is on the prediction of the re-circulation zone and the near wake; hence, further advancement in mesh design was not explored. However, mesh V produces very accurate results when compared to experimental data, not only for the near-wall force coefficients but also for the near wake velocity and Reynolds stresses. Figure 10 shows the variation of the re-circulation length with averaging time for all five meshes. The 0 on the *x* axis indicates the end of the flow initialization simulation and the start of the flow statistics collection simulation. The flow-field averaging is started at this point, and the longitudinal velocity is averaged in the spanwise direction at the end of every cycle to produce profiles similar to Fig. 9(a). It is observed from Fig. 10 that for all the meshes, the variation in the prediction of this length settles down after about 200 cycles. Hence, the values after 250 cycles are used as converged values and reported in Table V. These observations are similar to those by Parnaudeau.

In order to analyze the performance of the fine mesh (mesh V), the mean streamwise and transverse velocity and Reynolds stress profiles at three downstream locations, x/D = 1.06, 1.54, 2.02, were compared with the experimental data of Parnaudeau, DNS of Khan *et al.*^{39} and the LES of Lysenko *et al.*^{17} These three sets of data represent the most accurate results for flow past a cylinder at Re = 3900 from the literature. Figure 11 shows the results of current simulations compared with these three datasets. For the mean velocities, the current results show excellent agreement with both experimental and DNS data, suggesting that the rigorous grid design efforts enable significant improvements over previous LES results. The streamwise Reynolds stress at x/D = 1.06 shows some deviations from experimental and DNS data, which is attributed to the inability of the current LES to fully resolve the negative velocity due to re-circulation from the rear stagnation point to the point of maximum negative velocity in the flow (see Fig. 9). However, the results at the other two locations and the transverse Reynolds stress profiles show very high level of accuracy as compared to both DNS and experimental data. These results hence constitute one of the most accurate LES results for this flow situation to date. It should be noted that excellent agreement with the data of Parnaudeau *et al.* was also obtained by Moussaed *et al.*^{24} using a dynamic variational multiscale SGS model. However, only mean streamwise velocity profiles were reported while presently both streamwise and transverse velocity and their corresponding Reynolds stresses in the near wake have been validated.

## IV. RESULTS AND DISCUSSIONS OF NON-ISOTHERMAL CASES

Experiments on the flow past cylinders with a temperature difference, defined as the difference in cylinder wall temperature and freestream fluid temperature (Δ*T* = *T*_{w} − *T*_{∞}), were first performed by Schmidt and Wenner in 1941. They measured the variation of the Nusselt number around the cylinder for Δ*T* = 25 K. While a vast number of studies have yielded valuable predictive correlations across a large Re range, a very limited number of studies have included higher order statistics of the heat transfer behavior. One of the first to report fluctuating Nusselt numbers in a frequency spectrum important to turbulent crossflows was Scholten and Murray.^{51} This work was extended by Nakamura and Igarashi in 2004^{11} and represents the highest fidelity experimental data known to the authors for Nusselt number distributions in addition to the averaged quantities. Therefore, we employ these data to benchmark the numerical simulation performed presently at Δ*T* = 25 K. The most recent numerical study that tried to reproduce the work of Nakamura and Igarashi was by Jogee *et al.*^{12} They performed LES of non-isothermal flow past cylinders at four different Δ*T* of 25 K, 100 K, 200 K, and 300 K. It was found that up to 100 K, the numerical results were not significantly affected by Δ*T* and its influence was most pronounced at the highest temperature difference of 300 K. Hence, for the present study, the cases of Δ*T* of 25 K and 300 K were chosen to benchmark the simulations with the experimental data and more importantly provide insight into the flow and heat transfer physics across this temperature range. Only after a mature understanding is achieved can the ultimate goal be addressed, namely, the development of robust heat flux wall models to simulate non-isothermal flow past bluff bodies at higher Reynolds numbers. The main differences between the simulations of Jogee *et al.*^{12} and the current LES are a finer mesh (4.81M cells compared to 2.06M cells used by Jogee *et al.*) and the use of temperature-dependent thermophysical properties (Jogee *et al.* used constant properties). The same fine mesh used for the detailed isothermal flow simulation (mesh V) as previously described was employed for the non-isothermal cases. In addition to the Navier Stokes equations, the energy equation was also solved in a coupled manner with temperature-dependent fluid properties, as described in Sec. II. The color scheme employed to present the results for the three cases studied (one isothermal and two non-isothermal cases with Δ*T* = 25 K and 300 K) is kept the same all subsequent figures, i.e., Figures 12–27 with the exception of Figure 14. The current isothermal results are shown in green dashed-dotted lines with circle symbols, the non-isothermal case with Δ*T* = 25 K employs blue solid lines with square symbols, and the Δ*T* = 300 K case employs red solid lines with diamond shaped symbols. All experimental data are shown with black symbols, isothermal simulation results of other researchers are provided as black dashed lines, and the non-isothermal LES results of Jogee *et al.*^{12} (the only other LES of non-isothermal flow past a cylinder) are shown in the same color scheme as used by Jogee *et al.*, namely, magenta for Δ*T* = 25 K results and yellow for Δ*T* = 300 K results.

### A. Effect on cylinder forces and heat fluxes

The effect of non-isothermal flow on the friction coefficient and the pressure coefficient at the cylinder surface is illustrated in Fig. 12. The peak skin friction coefficient increases slightly from the isothermal case to the Δ*T* = 25 K case but shows a significant increase for Δ*T* = 300 K in the attached flow region. This increase is due to the increased viscosity of the fluid at higher temperatures. Since the boundary layer flow near the wall is dominated by viscous effects, the increase in the viscosity of air due to temperature has a pronounced effect in this region. After separation, the friction coefficient goes to 0 for all three cases. It is noted that the increase in the skin friction coefficient for Δ*T* = 300 K leads to a delayed separation of the boundary layer from the cylinder wall. This is inferred from the displacement of the skin friction profile toward the rear of the cylinder in the region where the friction coefficient goes to 0. Although the location of separation does not change appreciably from the isothermal flow (86.71°) to the Δ*T* = 25 K case (86.75°), there is a significant shift of this point toward the rear stagnation point for the Δ*T* = 300 K case (88.89°). Thus, a higher cylinder surface temperature causes a delayed separation of the boundary layer. The pressure coefficient also shows very weak sensitivity to the temperature field when the cylinder temperature is only 25 K higher than the freestream temperature. However, a significant increase in pressure is observed for the larger Δ*T* case. The back-pressure coefficient (pressure coefficient at the rear stagnation point) increases from the isothermal case (−0.908) by 1.2% (−0.897) for Δ*T* = 25 K and by 14.9% (−0.773) for Δ*T* = 300 K. Thus, the pressure at the cylinder surface is significantly higher when the cylinder is at a higher temperature. It is also noted that the pressure rise from the minimum value of the pressure coefficient to the back-pressure value is less steep for the higher Δ*T* case (−1.016 to −0.773) as compared to the lower Δ*T* case (−1.210 to −0.897). This lower adverse pressure gradient along the cylinder surface also explains the delayed separation for the higher Δ*T* case.

In order to benchmark present non-isothermal simulations, the variation of the mean Nusselt number (Nu) and rms of the fluctuating Nusselt number, normalized by the square-root of the flow Reynolds number, around the cylinder is compared with the experimental data of Nakamura and Igarashi.^{11} The current simulation is performed at Re = 3900 with Δ*T* = 25 K and 300 K, while Nakamura and Igarashi had conducted experiments at Re = 3000, 5900, 7400, 8900, 12 000, and 18 900, all of them at Δ*T* = 25 K. Figure 13 shows the current LES results, LES results of Jogee *et al.* (also performed at Re = 3900), and experimental data at Re = 3000 and Re = 5900 from Nakamura and Igarashi. Also shown in Fig. 13(a) is the Nusselt number at the rear stagnation point (Nu_{r}) based on correlation developed from experimental data as a function of Reynolds number from a separate experimental effort by Nakamura and Igarashi;^{52} values for Nu_{r} at three Reynolds numbers (3000, 3900, and 5900) are shown as dotted lines. It is seen that the experimental data show slightly lower values of Nu_{r} compared to the correlation, while the current LES slightly over-predicts this value. Since experimental data are not available at the flow Reynolds number of 3900, the close agreement with the correlation at the rear stagnation point and the fact that the mean Nusselt number trend lies within the experimental data of Re = 3000 and Re = 5900 in the separated region of the flow give good confidence in present simulations.

The experimental data show that the mean Nusselt number ⟨*Nu*⟩ decreases from a value of unity at the front stagnation point to about 90°. After reaching a minimum value, ⟨*Nu*⟩ increases toward the rear stagnation point. The minimum value of ⟨*Nu*⟩ is reached slightly after flow separation ensues. The decrease in ⟨*Nu*⟩ from the front stagnation point to the point of separation is due to the heating up of the fluid near the wall by conduction in the laminar boundary layer. This decreases the temperature gradient between the cylinder wall and the adjacent fluid, thereby leading to a reduction in ⟨*Nu*⟩ along the cylinder surface in the direction of flow. Nakamura and Igarashi mentioned that their data are slightly over-predicted near the front stagnation point since they do not correct for blockage effects, which were not insignificant in their studies.^{53} The experimental work of Schmidt and Wenner^{53} does not suffer from these effects and predict normalized ⟨*Nu*⟩ values slightly less than unity at the front stagnation point. This is also validated by the numerical work of Kim and Nakamura.^{27} After separation, the flow near the wall transitions to turbulent flow, which increases fluid mixing and transports heat away from the wall and into the bulk of the flow more effectively. This manifests as a rise in ⟨*Nu*⟩ toward the rear stagnation point. It is also noted that at higher Re, ⟨*Nu*⟩ reaches its minimum value earlier and the heat transfer recovery is enhanced due to greater mixing, leading to a higher value at the rear stagnation point. The rms of the fluctuating Nusselt number shows a very low constant value before separation and a monotonically increasing trend after separation takes place. It is observed that the location of minimum ⟨*Nu*⟩ lies between the locations of minimum ⟨*Nu*⟩ for Re = 3000 and Re = 5900 observed by Nakamura and Igarashi. However, it is noticed that both the minimum value of ⟨*Nu*⟩ itself and the increase in ⟨*Nu*⟩ just after the minimum value are closer to the data for Re = 5900. This suggests that the predicted flow transition from laminar to turbulent at the wall is slightly accelerated, i.e., the transition occurs faster than that observed in experiments, thereby preventing ⟨*Nu*⟩ from reaching its minimum value and causing an over-prediction of ⟨*Nu*⟩ just after. This suggests that sub-grid models for LES still require some fine tuning in predicting transition to turbulence where most approximations break down. The increasing trend of ⟨*Nu*⟩ toward the rear stagnation point is well captured by the current LES. It is noted that the value of Nu_{r} is predicted 4.2% larger than the value from the correlation expression shown in Fig. 13. It is observed that the current LES results match quite closely with the simulation by Jogee *et al.* up to the location of minimum ⟨*Nu*⟩. However, the LES of Jogee *et al.* predicts a steeper increase in ⟨*Nu*⟩, which then remains constant in the separated region of the flow. This over-prediction could be due to the coarser mesh used by those authors, thereby relying more on the sub-grid model and not resolving the essential flow structures during transition. For the higher Δ*T* = 300 K case, the normalized ⟨*Nu*⟩ has a value slightly larger than unity at the front stagnation point due to the increased thermal conductivity of the fluid at a higher temperature. A decreasing trend of ⟨*Nu*⟩ is observed up to flow separation, similar to the lower Δ*T* = 25 K case. Since the flow separation is delayed due to reasons already noted, a larger departure from experimental data is observed in the transition region. Once the flow becomes turbulent near the wall, the increasing trend of ⟨*Nu*⟩ is observed again until the rear stagnation point. Although in the 120°–150° range the ⟨*Nu*⟩ for both Δ*T* cases are very similar, the value of the predicted Nu_{r} for the higher Δ*T* case is only 0.6% larger than the correlation value, suggesting a decrease in Nu_{r} with higher Δ*T*. Both the shift in the transition point toward the rear and the decrease in Nu_{r} suggest that the higher Δ*T* case behaves like a flow with a lower Reynolds number. This is mainly due to the increased viscosity at elevated temperatures.

The fluctuating Nusselt number profiles in Fig. 13(b) show values very close to zero at the front stagnation point from both LES simulations (Jogee and current) as one would expect for the known laminar flow conditions in that region. However, the non-zero values from experimental data tell a different story, one that is much more difficult to explain with flow physics. Fluctuation signatures in the experimental data are due to both temperature fluctuations and unavoidable random noise in the sensor signal, and it is difficult to isolate those effects without a rigorous calibration of the temperature sensors. For turbulent flow physics, the electrical noise is small compared to the fluctuations derived from the flow condition, but the same is not true for laminar flows. Therefore, the near constant fluctuation intensity of the experimental data up to about 90° could likely be attributed to sensor noise. In both simulations, a slight increase in fluctuating Nusselt number is observed from about 45°, which could be due to fluctuations in the separated shear layer.^{27} For current simulations, the sharp increase in the fluctuating Nusselt number starts from the point of separation due to transition to turbulent flow. However, the simulations by Jogee *et al.*^{12} show that the sharp increase occurs from a point corresponding to the minimum ⟨*Nu*⟩ of their results. This discrepancy could also be due to the differences in the mesh. Both sets of LES results show a much steeper rise in the fluctuating Nusselt number, indicating a quicker transition to turbulence than is evidenced by experimental data. As before, this could be attributed to the inaccuracy of SGS models of LES to accurately capture transition flow. It is noted that both simulations predict higher fluctuations for the lower Δ*T* case. The reduction in fluctuations for the higher Δ*T* case can be explained by the reduction in effective Reynolds number due to increased viscosity of the flow at elevated temperatures. This justifies the use of temperature-dependent thermodynamic properties for our current simulations. However, it is to be noted that the simulations of Jogee *et al.* do not consider this effect and hence could result in some of the differences.

It is found that the fluid properties play an important role in the mean and fluctuating Nusselt number prediction. In order to understand their impact on the flow and the differences in predictions with the LES data of Jogee *et al.*, additional simulations were performed with constant fluid properties (*μ*, *c*_{P}, and *k*) for the higher Δ*T* case where the range of thermophysical property values in the flow is greatest. Three new simulations were performed only allowing fluid density *ρ* to change with temperature according to the ideal-gas equation of state, consistent with the simulations by Jogee *et al.* The other properties were separately evaluated at T = *T*_{∞}, *T*_{w}, and *T*_{f} [note that the film temperature is defined as *T*_{f} = 0.5(*T*_{∞} + *T*_{w})] for the three unique simulations, respectively. These were compared to the results using variable fluid properties after each simulation was run for ten vortex shedding cycles only. Figure 14 shows the mean Nusselt number for these four simulations. It is observed that when the properties are fixed at *T*_{∞} and *T*_{w}, the ⟨*Nu*⟩ in the attached region of the flow is significantly underpredicted and overpredicted, respectively. However, when the properties are fixed at the film temperature, the results are quite agreeable with the simulation using variable fluid properties. ⟨*Nu*⟩ at the front stagnation point is slightly smaller than that predicted from the simulation with variable properties since the conductivity is evaluated at T = *T*_{f} = 450 K for the former and at T = *T*_{w} = 600 K for the latter. As mentioned before, all thermophysical properties except density increases with increasing temperature for air. Thus, the simulation with properties evaluated at *T*_{f} predicts lower heat flux due to the lower thermal conductivity and the graph moves closer to the case of lower *ΔT*. It is also observed that in the separated region of the flow, the ⟨*Nu*⟩ prediction does not follow a monotonic trend with the temperature at which the properties are evaluated. This is because the changes in properties also impact the flow in the wake region, which, in turn, affects the predictions at the cylinder surface as well as the re-circulation bubble formation.

In order to gain further insight into the effect of constant thermophysical properties, these simulations need to be performed for more vortex shedding cycles to gather meaningful statistics. As this was not the main focus of the current paper, further cycles were not gathered. However, it can be concluded that considering realistic fluid properties is of utmost importance for accurate simulation of the non-isothermal flow past a cylinder. Further, since the cases with constant fluid properties can predict the increasing trend of ⟨*Nu*⟩ in the separated region of the flow, flattening of ⟨*Nu*⟩ in the results of Jogee *et al.* can be attributed mostly to the coarser mesh used, thereby underpinning the importance of mesh design.

### B. Effect on flow field

The influence of the non-isothermal flow on the flow development itself is studied next. It is observed that the effect of temperature on viscosity dictates the near-wall flow development and hence the forces on the cylinder and the heat fluxes. However, as we move away from the cylinder and into the wake region, the effect of temperature on other properties needs to be considered as well. As the temperature of the flow is increased, one could have expectations based on two different trains of thought. First, we recognize that the fluid density decreases as the temperature is raised (following the ideal-gas law with no pressure change due to a low Mach number approximation). This, coupled with the increase in fluid dynamic viscosity, yields a significant increase in the kinematic viscosity (*υ* = *μ*/*ρ*) of the fluid, thereby reducing the local flow Reynolds number (*Re*_{l} = *ul*/*υ*) and suppressing turbulence. Alternatively, one could recognize that the decrease in density would bring about an increase in velocity in order to maintain continuity in the flow, which would then lead to increased flow turbulence. In order to test how reality balances these competing effects, we analyze the mean streamwise velocity and the turbulent kinetic energy [$TKE=0.5u\u2032u\u2032\u0303+v\u2032v\u2032\u0303+w\u2032w\u2032\u0303$] at the wake centerline (results are plotted in Fig. 15). Looking at the streamwise velocity data and assessing L_{rc}, there is a slight increase in this value for the current LES results between isothermal and *ΔT* = 25 K (6.6%) but very little shift between 25 K and 300 K (increase by 1.8%). See Table VI for the exact values of L_{rc}. The increase in L_{rc} with increasing *ΔT* is also suggestive of the apparent low Re behavior of the higher *ΔT* case similar to the observations of the near-wall phenomena. This indicates that the changes in properties due to increasing *ΔT* dictate not only the near-wall flow development but also the flow in the near wake. The decrease in density and increase in viscosity of the fluid at regions of high temperatures lead to a reduction in the local flow Reynolds number, thereby causing a larger re-circulation zone. However, the same argument does not explain the smaller increase in L_{rc} from the low *ΔT* case to the high *ΔT* case as compared to the increase from the isothermal case to the low *ΔT* case. In order to explain this difference, we need to resort to the second train of thought. For the higher *ΔT* case, the wall temperature is double that of the freestream temperature; hence, the fluid density is almost halved. This results in local flow acceleration in order to maintain continuity and hence an increase in turbulence. At higher temperatures, this effect starts to become more important and opposes the effect of turbulence suppression due to increased viscosity. By comparison, those same values of *ΔT* suggest a substantial increase in L_{rc} when analyzing the results of Jogee *et al.* Two main attributes of the current work help explain these differences: careful mesh design including rigorous isothermal validation in the present simulations and the use of temperature-dependent fluid properties. For TKE, a peak is observed just beyond x = L_{rc}, and the inconsistencies seen for the horizontal shift and its dependence on *ΔT* are largely explained by the difference in predictability of L_{rc} as noted. The TKE profile moves downstream with a reduction in the peak value from the isothermal case to the lower *ΔT* case because of the increase in L_{rc} due to turbulence suppression. The increase in TKE for the higher *ΔT* case further strengthens the arguments provided before of enhanced turbulence production at *ΔT* = 300 K. It is noted that the TKE further downstream (beyond x/D = 4) shows an increased value for the higher *ΔT* case. Hence, turbulence production at elevated temperatures is not only constrained in the vicinity of the re-circulation zone but also further downstream. There is a significant difference in the peak TKE values themselves with the current simulations, suggesting an increase greater than 50% compared to those of Jogee *et al.* This large difference in magnitude is because the TKE is not scaled by the freestream velocity and represents dimensional quantities that depend on the value of D and the freestream fluid properties since *U*_{∞} = *μ*_{∞}*Re*/*ρ*_{∞}*D*. Thus, although the same Re is used in both simulations, the local freestream velocity could be quite different, thus leading to different values of TKE (m^{2}/s^{2}).

. | L_{rc}/D
. | X/L_{rc} (for x/D = 1.06)
. | Current x/D . | X/L_{rc} (for x/D = 1.54)
. | Current x/D . | |||
---|---|---|---|---|---|---|---|---|

. | Jogee . | Current . | Jogee . | Current . | at Jogee x/L_{rc}
. | Jogee . | Current . | at Jogee x/L_{rc}
. |

Isothermal | 1.68 | 1.480 | 0.486 | 0.535 | 0.963 | 0.706 | 0.778 | 1.399 |

∆T = 25 K | 1.74 | 1.577 | 0.473 | 0.510 | 0.983 | 0.688 | 0.741 | 1.428 |

∆T = 300 K | 2.04 | 1.606 | 0.417 | 0.503 | 0.879 | 0.606 | 0.731 | 1.277 |

. | L_{rc}/D
. | X/L_{rc} (for x/D = 1.06)
. | Current x/D . | X/L_{rc} (for x/D = 1.54)
. | Current x/D . | |||
---|---|---|---|---|---|---|---|---|

. | Jogee . | Current . | Jogee . | Current . | at Jogee x/L_{rc}
. | Jogee . | Current . | at Jogee x/L_{rc}
. |

Isothermal | 1.68 | 1.480 | 0.486 | 0.535 | 0.963 | 0.706 | 0.778 | 1.399 |

∆T = 25 K | 1.74 | 1.577 | 0.473 | 0.510 | 0.983 | 0.688 | 0.741 | 1.428 |

∆T = 300 K | 2.04 | 1.606 | 0.417 | 0.503 | 0.879 | 0.606 | 0.731 | 1.277 |

The streamwise and transverse Reynolds stress terms in Fig. 16 show some trends similar to the TKE profile while bringing out some differences in the contributions of these terms to the overall TKE. The predicted streamwise Reynolds stress shows a shift of the profile away from the cylinder for the lower *ΔT* due to an increased re-circulation length. The profile for the higher *ΔT* shows a broadening effect with a larger (almost 20%) peak, indicating increased turbulence production. This 20% increase in the streamwise stress at the edge of the re-circulation zone explains the increase in TKE at higher *ΔT*. As discussed before, this is likely due to enhanced turbulence as a result of lower fluid density. It is also observed that after a distance of 4D from the center of the cylinder, all three profiles merge. Thus, the enhancement of the streamwise Reynolds stress at high *ΔT* is restricted to the near wake region. The transverse Reynolds stress profiles show trends very similar to the TKE. The increased stress beyond 4D for the higher *ΔT* demonstrates that the increase in TKE far downstream is due to the enhanced transverse turbulent fluctuations. This should lead to increased turbulent mixing further away from the cylinder due to a strongly non-isothermal flow.

From the discussion of results for the isothermal flow simulation, it was identified that the larger second peak in the isothermal simulations of Jogee and Parnaudeau was tied to the over-prediction of the re-circulation length. Since the re-circulation zone is a region of large shear stresses (see Fig. 8), most of the turbulence production due to shear in the mean flow occurs here. This is further evidenced by the increasing trend of TKE (shown in Fig. 15) along the centerline, within the re-circulation zone. It is noted that the TKE reaches a peak at the edge of the re-circulation zone where the two counter-rotating pairs interact with the outer flow. In the far-wake region, the TKE starts falling off but does not go to zero due to the vortex shedding phenomenon. Figure 16 shows that the peak in TKE at the edge of the re-circulation zone is due to the peak in the streamwise Reynolds stress. Hence, large re-circulation length simulations lead to a continued rise in the streamwise fluctuations until the end of the zone, thereby causing over-predictions in the Reynolds stresses further away from the cylinder. However, for cases of non-isothermal flow, the increase in the re-circulation length is due to changes in the fluid thermophysical properties. Although this increase is due to local turbulence suppression, the effect of an increased re-circulation length is an increase in the total TKE as observed above. Since it is known that increased levels of turbulence lead to shorter re-circulation zones,^{19} this enhanced TKE prevents further growth of the re-circulation length. Thus, giving further evidence for the smaller difference of re-circulation lengths between the lower *ΔT* and higher *ΔT* cases.

It was identified by Parnaudeau *et al.*^{18} and supported by present simulations that the re-circulation zone prediction can be considered one of the most important flow features for the flow past a cylinder. Predictions of mean velocities and Reynolds stresses depend greatly on this prediction. After having analyzed the flow behavior in the streamwise direction, the mean velocity and Reynolds stress variations in the transverse direction are presented in Figs. 17 and 18, respectively. It is observed that the results of Jogee *et al.* indicate that the flow is much more sensitive to *ΔT* than current results. However, this is mainly due to the discrepancies in the re-circulation bubble length prediction. The current results hardly show any differences in the streamwise velocity profiles for the three cases, although some differences, especially near the centerline, are observed for the transverse velocity profiles. At the location x/D = 1.06, the transverse velocity in the region |y/D| < 0.5 decreases when a small *ΔT* is introduced in the flow but increases when the *ΔT* is increased. Hence, the hypothesis of local flow acceleration due to increased *ΔT* is validated with these results. However, the flow in the |y/D| < 0.5 region for x/D = 1.54 is quite complex as the flow changes sign three times for the non-isothermal flow cases. The current results also show a non-monotonic behavior of the Reynolds stress variation with increasing *ΔT*. Higher Reynolds stresses for the larger *ΔT* as compared to the smaller *ΔT* case at all downstream locations further indicate the effect of local enhancement of turbulence.

### C. Re-circulation length scaling

In order to understand whether the differences in results are governed mainly by the considerations of variable properties or due to the mesh design, an analysis of the flow inside the re-circulation zone is performed. Table VI shows the re-circulation length predictions for the isothermal and two non-isothermal cases by the LES of Jogee *et al.* and current LES. The results presented before are at two downstream locations, x/D = 1.06 and 1.54. However, these correspond to different locations with respect to the position within the re-circulation bubble of each case, which affects the flow profiles extensively.^{54} Hence, a new parameter x/L_{rc} is introduced. It is recalled that the domain origin is at the center of the cylinder. Hence, the distance x is measured from the center of the cylinder. However, the re-circulation length is measured from the rear stagnation point conventionally. Thus, this new parameter is evaluated by considering the re-circulation length from the center of the cylinder as x/L_{rc} = (x/D)/(L_{rc}/D + 0.5). It is seen that the same x/D = 1.06 corresponds to different x/L_{rc} for the three cases due to the different re-circulation bubble lengths (see Table VI for values). In order to compare the results of Jogee with present simulations at the same relative locations within the re-circulation zone, the x/D locations that correspond to the x/L_{rc} of Jogee *et al.* at x/D = 1.06 are calculated. Similar operation is performed for the other downstream location x/D = 1.54. The calculations for both these locations using the relations between x/D and x/L_{rc} discussed above are shown in Table VI. Figure 19 shows the mean velocities and mean Reynolds stresses at these x/D locations. It is seen that the mean streamwise and transverse velocity predictions agree very closely between the two LES results (that of Jogee *et al.* and current). This implies that the relative velocity distribution within the re-circulation zone is quite accurately captured by both simulations. Notably, the large sensitivity to *ΔT* observed by Jogee *et al.* simulations is reproduced by present simulations when the flow field is probed at matching locations with respect to the re-circulation zone lengths of each case. Whereas the transverse velocity profiles for the high *ΔT* case show very good agreement at both locations, the isothermal and low *ΔT* case exhibit some differences near the centerline, in the region |y/D| < 0.5. Although the re-circulation length scaling yields similar velocity profiles between the two simulation studies, the same cannot be said of the Reynolds stress profiles. While similar trends are observed for both simulations, i.e., a decrease in turbulent stress with increasing *ΔT*, the magnitudes and variations predicted are quite disparate. The Reynolds stresses are more sensitive to the correct prediction of the re-circulation length itself. This is because the size of the re-circulation zone controls the shear production of turbulence, as shown in Figs. 15 and 16. Thus, although the velocity profiles within the re-circulation zone are relatively self-similar, different predicted zone lengths lead to changes in the prediction of turbulent stresses. Furthermore, the transverse Reynolds stress predicted by Jogee shows a profile very different from the current simulation, while the profiles for the streamwise Reynolds stress are similar in shape but varying in magnitude. Thus, the transverse Reynolds stress that is mainly involved in lateral mixing (as opposed to streamwise transport) depends more strongly on the accurate prediction of the re-circulation length.

The previous analysis shows that a better physical understanding of the flow can be obtained by looking at the velocity and Reynolds stress profiles at the same relative location in reference to the predicted re-circulation zone. Figure 20 shows the variation of the streamwise velocity and the turbulent stresses along the centerline for each of the three cases, with the distance from the center of the cylinder normalized by the re-circulation length of each case. Hence, the streamwise velocity for each case passes through (abscissa, ordinate) = (1, 0). It is noted that the streamwise velocity for all three cases agrees very closely after the point of minimum velocity, which occurs at the same location in the x/L_{rc} space. Between the rear stagnation point and the minimum velocity location, the flow toward the cylinder is decelerated for the lower *ΔT* case and accelerated for the high *ΔT* case. This acceleration due to density reduction is also tied to the increase in turbulent stresses noted on the right-hand side of Fig. 20. It is noted that the turbulent stress variation for the isothermal case and the lower *ΔT* case are quite identical in the x/L_{rc} space. However, at higher *ΔT*, there is a significant increase in the turbulent stresses. Furthermore, it is observed that the peak TKE for the high *ΔT* case moves closer to the cylinder and lies within the re-circulation zone (at x/L_{rc} < 1), while the peak for the isothermal and low *ΔT* case lies at the edge of the zone (at x/L_{rc} ∼ 1). The increase in TKE is mainly due to the increase in the transverse Reynolds stress, and this trend continues beyond the re-circulation zone as noted earlier in Fig. 16. However, the increase in TKE at the higher *ΔT* case is larger within the re-circulation zone than outside it. This is because the fluid is heated up by the hot cylinder more in this region, thus causing a greater reduction in fluid density.

In Table VI, the locations x/D = 1.06 and 1.54 refer to the locations x/L_{rc} = 0.535 and 0.778, respectively, for the current isothermal flow simulations. Using the L_{rc}/D values of the non-isothermal cases from Table VI, the data for the non-isothermal cases are extracted at these two x/L_{rc} and compared with the isothermal case in Fig. 21. The streamwise velocity shows very small differences for the three cases, while the transverse velocity is more sensitive to *ΔT*. This is because of the larger increase in transverse Reynolds stresses and the complex flow behavior noted earlier in the region |y/D| < 0.5. The Reynolds stresses show that the isothermal stresses and the turbulent stresses at the lower *ΔT* case follow similar profiles with a slight increase for the lower *ΔT* case, especially for the transverse Reynolds stress. However, the stresses for the higher *ΔT* case show a significantly larger magnitude. It is also noted that the stress profiles are broader for the higher *ΔT* case. This indicates that the re-circulation zone is larger in the transverse direction as compared to the re-circulation transverse width of the isothermal and lower *ΔT* case. Thus, the increased transverse stresses and velocities for the case of greater temperature differences in the flow also result in broadening of the re-circulation zone in the transverse direction. It can be concluded that the new scaling x/L_{rc} implicitly considers the first order effects of changes in the flow due to changes in thermophysical properties through non-dimensionalization using the re-circulation length itself. The results analyzed using this scaling provides greater insight into the differences in results for different non-isothermal flow cases, which arise due to the non-linear coupling of flow and heat transfer. At elevated temperatures, the transverse Reynolds stresses increase significantly and hence aid in heat transfer away from the centerline. This causes an expansion of the width of the re-circulation zone and hence enhanced turbulence.

### D. Effect on thermal field

Next, the temperature field is investigated in more detail for the two non-isothermal cases. The increase in the temperature of the flow is quite low for *ΔT* = 25 K; however, this still increases the re-circulation length due to the changes in the thermophysical fluid properties and hence the flow field. Figure 22 shows the increase in the temperature of the fluid for both the non-isothermal cases. It is seen that at a distance of x/D = 1.06, the maximum increase in the temperature of the fluid for the higher *ΔT* case is about 25%. Similarly, the rms of the fluctuating temperature field is also significantly higher for the larger *ΔT* case. The peaks of both of these profiles occur before |y/D| = 1 at both downstream locations and correspond to the edge of the re-circulation zone. The temperature profile for the larger *ΔT* case shows a steep gradient in temperature due to mixing of the cold fluid (almost at freestream temperature) outside and the hot fluid inside the re-circulation zone. This gradient causes strong heat fluxes in this region (see Fig. 24) and is accompanied by large streamwise Reynolds stresses (see Fig. 16). These combine to give rise to the peaks in the temperature fluctuations. Differences in the results with Jogee *et al.* profiles arise due to similar reasons as described before, i.e., mainly due to the differences in the re-circulation length prediction. When the current results are extracted at the same x/L_{rc} as those of Jogee *et al.*, similar mean and rms temperature profiles are observed, although the magnitudes remain higher in present computations due to the higher stresses in the shorter re-circulation zone. Figure 23 shows that there is a sharp rise in the rms of temperature very close to the cylinder but falls off rapidly as well. This peak corresponds to the sharp rise in TKE near the cylinder, as shown in Fig. 16. In order to further understand the temperature fluctuations along the centerline, the fluctuation intensities are studied. A temperature fluctuation intensity, $IT=T\u2032/T\u221e=T\u2032T\u2032/T\u221e$, is defined analogous to the definition of turbulence intensity, $I=u\u2032/U\u221e=2TKE/3/U\u221e$. The variation of these quantities for each case is shown along the centerline on the right-hand side of Fig. 23. It is seen that the peak temperature fluctuation intensity near the rear stagnation point of the cylinder corresponds to the sharp rise in turbulence intensity there. The peak intensity for the higher *ΔT* case is almost 12 times larger than the peak for the lower *ΔT* case. This is due to the enhanced turbulent heat transfer (see Fig. 24) from a larger temperature difference and is intimately tied to the increase in turbulence intensity for the high *ΔT* case near the rear stagnation point. After this initial peak, the temperature fluctuation intensity reduces as the temperature gradients in the flow weaken, as evidenced by the flattening of the temperature profile further away from the cylinder in Fig. 22. However, a secondary, much smaller peak is observed near the edge of the re-circulation zone where the turbulence intensity is maximum and heat transfer is enhanced due to the mixing of the colder fluid outside the re-circulation zone with the hotter fluid inside. This is further evidenced by the negative peak in the streamwise turbulent heat flux in Fig. 26. It is also noted that the second peak moves upstream toward the cylinder for the higher *ΔT* case similar to the shift in the turbulence intensity peak. Thus, the turbulent fluctuations and the temperature fluctuations are affected by one another.

The increased temperature fluctuation and turbulent heat flux near the cylinder are responsible for heating up of the fluid downstream of the cylinder. Figure 24 shows the turbulent heat fluxes for the two non-isothermal cases at the two downstream locations. The profiles show peaks in locations similar to the rms of temperature, indicating that the peaks arise at the edge of the re-circulation zone. The negative peaks denote the absorption of heat by the cold fluid, while the positive peaks denote the transfer of heat from the hot fluid inside the re-circulation zone. The transverse heat flux has a value of zero at the centerline as the transverse velocity is also zero there. The opposite signs of the transverse turbulent heat flux denote the transfer of heat in opposite directions. Turbulent heat flux profiles follow the results of Jogee *et al.* quite closely. It is noted that there is a great difference in the turbulent heat flux magnitudes between the two cases of *ΔT* = 25 K and 300 K. The higher *ΔT* case shows much larger values of the turbulent heat flux as the temperature gradients are larger, thereby causing increased transfer of heat due to turbulent fluctuations in the flow.

### E. Friction scaling for turbulent heat flux

A more in depth understanding of the physics is required in order to contribute to the ultimate goal of the current efforts, namely to provide a framework for the development of thermal wall models. In particular, we aim to develop a thermal analog to hydrodynamic wall models gaining traction in WMLES efforts. An excellent review of recent progress on those efforts is from Bose and Park.^{55} For the thermal wall models, this represents an important step in achieving significant reductions in computational cost while maintaining the high degree of fidelity and accuracy offered by LES in a wide range of important applications. Implicit in a useful wall model is the ability to quantitatively express flow metrics across as wide of operating conditions as possible. While the current study is conducted at a single Reynolds number, the range of *ΔT* considered here, and more importantly, the worth that is placed on our validated results, enables a first step toward the ultimate goal. The variation noted in the turbulent heat flux profiles for different *ΔT* suggests that finding self-similar behavior would be challenging. The key to overcoming that challenge is to normalize the profiles by turbulence-based quantities. This leads to looking at friction quantities at the boundary rather than the freestream conditions. The friction velocity *u*_{τ} and friction temperature *T*_{τ} [see Eq. (17.45) of Ref. 56) are given by

where *τ*_{w} is the wall shear stress and *q*_{w} is the wall heat flux, with *ρ* and *c*_{P} representing the density and specific heat capacity of the fluid, respectively. These are used to normalize the turbulent heat fluxes in place of the freestream quantities to ascertain whether better scaling can be obtained.

Figure 25 shows the distribution of these two quantities around the cylinder for the two non-isothermal cases. The differences in the normalized *T*_{τ} profiles are greater than the normalized *u*_{τ} profiles since the two cases differ in *ΔT* but not in *U*_{∞}. However, it is noteworthy that the differences in *T*_{τ} are mainly after separation, while the differences in *u*_{τ} are in the attached region of the flow. The former is because of different heat transfer characteristics in the separated region, while the latter is mainly due to thermodynamic property changes with temperature. It is expected that the *u*_{τ} profiles would show greater differences for cases at different Reynolds numbers of the flow, which is the basis of ongoing and future work. Both profiles show a shift in the separation point toward the rear of the cylinder for the higher *ΔT* case. At this point, there are local minima in *u*_{τ} and local maxima in *T*_{τ} due to flow detachment from the surface and transition from a laminar boundary layer to a turbulent one. Further away from the point of separation, the *u*_{τ} decreases due to the reduction in skin friction and the *T*_{τ} increases as turbulent fluxes aid the wall heat transfer. It is seen that *u*_{τ} reaches a constant value at about 135° while the variations in *T*_{τ} also decrease but show small variations after this point. These variations are due to the high Nusselt number fluctuations shown in Fig. 13. Hence, the values at the rear stagnation point are used to normalize the turbulent heat flux variables in the wake region instead of *U*_{∞} and *T*_{∞}. Figure 26 shows the variation of the streamwise component of turbulent heat flux along the centerline. When normalized by freestream quantities, the higher *ΔT* case shows very large turbulent heat fluxes due to the larger temperature gradients and the increased turbulence intensity of the flow. The turbulent heat flux increases sharply from the rear stagnation point due to the increased turbulence intensity and temperature fluctuation intensity in the flow (see Fig. 23). Within the re-circulation zone, the heat flux remains positive, indicating heating up of the fluid due to heat transfer from the cylinder to the flow. This heat transfer decreases away from the cylinder due to the decrease in temperature fluctuations and goes to zero outside the re-circulation zone, indicating the transfer of heat away from the fluid to the colder fluid in the neighboring regions outside the zone. The differences with the results of Jogee *et al.* are mainly due to the different re-circulation length predictions and because the increased turbulence at the high *ΔT* case is not adequately captured by the LES of Jogee *et al.* When scaled by the friction quantities at the rear stagnation point, there is a much better collapse of the turbulent heat flux profiles from the two non-isothermal cases. The peak near the rear stagnation point reaches a value close to 1 for both cases: 0.93 for the lower *ΔT* case and 1.08 for the higher *ΔT*. Further downstream, there is a very good agreement of these profiles for both non-isothermal cases with slightly larger values for the high *ΔT* case throughout. This gives great confidence in the scalings adopted presently.

The same scaling is used to normalize the heat flux quantities in the near wake. Figure 27 shows that even downstream of the flow, these profiles collapse for both non-isothermal cases. All profiles are of order one, thereby indicating the correct choice of the scaling variables. The higher *ΔT* case shows slightly more broadening of the profiles. This is due to the broadening of the re-circulation zone from enhanced turbulence intensity. The scaling adopted presently implies that the turbulent heat flux in the near wake flow is well correlated with the wall heat flux since Eq. (13) can be re-arranged as *u*_{τ}*T*_{τ} = *q*_{w}/(*ρc*_{P}). Thus, the new scaling *u*_{τ}*T*_{τ} differs from the wall heat flux only through constant (although temperature dependent) fluid properties. This shows that the turbulent heat flux behavior in the near wake region is tightly coupled to the wall heat flux, irrespective of the *ΔT* in the flow. In order to test the applicability of the friction scaling in the near wake region of the flow, the contours of the scaled streamwise and transverse turbulent heat flux for both cases are shown in Fig. 28. It is seen that at all points in the wake region, the values stay within the range of ±2. The turbulent heat flux characteristics flip signs at the edge of the re-circulation zone for all cases shown here. In addition, the major turbulent transport of heat in the transverse direction occurs outside the re-circulation zone.

## V. CONCLUSIONS

The first part of the work presented here focuses on the improvement of large-eddy simulations of flow past a circular cylinder using rigorous grid design. Grid aspect ratios (ratio of circumferential cell width to wall normal cell height) of 0.2 and higher are required near the cylinder wall to capture the boundary layer and its separation from the surface in order to yield accurate predictions of pressure and skin friction distributions. Expansion ratios of about 1% and 3%, respectively, are required in the near wake (up to three times the diameter of the cylinder downstream) and far wake. These restrictions yield predictions very close to the well-validated experimental results of Parnaudeau *et al.*^{18} With these grid considerations, the velocity and shear stresses within the re-circulation zone are found to match experimental data extremely well, much better than any other LES results known to the authors. This is achieved through the rigorous grid design treatment, which enabled a more accurate prediction of the re-circulation length compared to other researchers. This is also because averaging of the results was performed over 250 vortex shedding cycles following the recommendations of Parnaudeau *et al.* to obtain a converged value of the re-circulation length. It has been demonstrated that this reduction in re-circulation length also causes an increase in predictions of mean drag and rms of lift.

In non-isothermal flow simulations, while the improved near-wall grid design has important implications for the accuracy of heat transfer predictions from the cylinder, the near wake grid design governs the re-circulation length predictions through resolution of shear layer instabilities. Additionally, it is also important to account for the temperature dependence of thermophysical quantities, especially in highly non-isothermal flow situations. It has been found that both grid design and variable properties play an important role in establishing transition to turbulence due to the fact that flow separation is delayed by at least 2° when the cylinder temperature is double that of the freestream temperature of 300 K as compared to the isothermal flow case. In the wake region of the flow, the re-circulation zone is broadened for non-isothermal flow as the fluid is heated up, causing increased viscosity and lower densities. However, as the *ΔT* (difference in cylinder temperature and freestream temperature) is increased, the effect of broadening reduces. Current simulations indicate that this is due to broadening of the re-circulation zone in the transverse direction and enhanced turbulent heat transfer, which leads to increased turbulence intensity. Higher turbulence intensity implies shorter re-circulation zone lengths similar to the reduction of the re-circulation zone at higher Reynolds numbers. It has been identified that due to these competing effects of non-isothermal flows, comparisons of flows with different *ΔT* are more insightful when the streamwise direction is non-dimensionalized by the re-circulation length rather than the more conventional choice of the cylinder diameter. This allows the conclusion that while the mean flow characteristics are reasonably self-similar within the re-circulation zone, the Reynolds stresses depend on the length of the zone itself. These stresses are found to increase with increasing *ΔT* at a particular x/L_{rc} location due to the increase in turbulence intensity at higher *ΔT*. Increased turbulence intensities and increased temperature gradients in the flow at high *ΔT* cause increased turbulent heat fluxes in the wake region of the flow. Furthermore, driven by the broader goal of developing heat-flux wall models for use in WMLES at higher Re, an attempt was made to relate these turbulent heat fluxes to the wall heat transfer. It was found that the introduction of the new friction scaling *u*_{τ}*T*_{τ} to normalize the heat fluxes for the two studied cases of *ΔT* = 25 K and 300 K not only brings both the streamwise and normal turbulent heat fluxes to values of order one but also collapses the profiles for both cases. The new friction scaling term *u*_{τ}*T*_{τ} depends directly on the wall heat flux. This presents the most novel finding of the current paper, i.e., the turbulent heat fluxes in the wake region of the flow scale directly with the wall heat flux (and thermodynamic quantities *ρc*_{P}) irrespective of the temperature of the cylinder wall.

The correlation of turbulent heat fluxes in the flow with the wall heat transfer over a large range of *ΔT* gives great confidence in the ability to develop heat-flux wall models for flow past bluff bodies. Currently, wall-resolved simulations at higher Reynolds numbers using the grid design procedure outlined here are ongoing to test the sensitivity of these scalings to the Reynolds number. A wall-modeled framework based on the implementations by Mukha *et al.*^{57} has also been developed that probes the near-wall flow to model the wall shear stress at high Re. These will be benchmarked against the recent high Re experimental data from the work of Desai *et al.*^{58} and compared to other simulation methods such as detached eddy simulations^{59} and techniques such as adaptive mesh refinement^{25} to assess the computational savings from using wall models. After these new friction scalings have been tested over a range of Re and a range of *ΔT*, efforts will be made to come up with a heat-flux wall model in the same framework as the wall shear stress model. This is expected to significantly decrease the computational cost of wall-resolved simulations at Reynolds numbers relevant to gas turbines while maintaining accuracy.

## ACKNOWLEDGMENTS

We gratefully acknowledge financial support through the Small Business Technology Transfer (STTR) Phase I program from the U.S. Navy under the project titled “A Framework for Modeling Turbulent Flow with Conjugated Heat Transfer,” Contract No. N6833519C0775, awarded to Advanced Cooling Technologies and Texas A&M University. This research was driven by the need to address the gap in the available data observed during the Phase I STTR program. Interactions with Dr. John Willard and Dr. Allan Aubert, the technical points of contact at the Naval Air Warfare Center—Aircraft Division, Propulsion and Power, helped to enhance the long-term relevance of our efforts. We also recognize the contributions of Texas A&M University’s High-Performance Research Computing (HPRC) who supplied the requisite resources for this computational investigation. The early stages of this project also partially benefited from the high-performance computing infrastructure of the Extreme Science and Engineering Discovery Environment (XSEDE), Allocation Grant No. DMR180017. The XSEDE is supported by the National Science Foundation (NSF), Grant No. ACI-1053575. NAVAIR Public Release 2020-623. Distribution Statement A - Approved for public release; distribution is unlimited.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

^{3}to 10

^{5}