Structure of the thermal boundary layer in turbulent channel flows at transcritical conditions

Abstract Enhanced fluctuations, steep gradients, and intensified heat transfer are characteristics of wall-bounded turbulence at transcritical conditions. Although such conditions are prevalent in numerous technical applications, the structure of the thermal boundary layer under realistic density gradients and heating conditions remains poorly understood. Specifically, statistical descriptions of the temperature field in such flows are provided inconsistently using existing models. To address this issue, direct numerical simulations are performed by considering fully developed transcritical turbulent channel flow at pressure and temperature conditions that cause density changes of a factor of up to $O(20)$ between the hot and cold walls. As a consequence of the proximity of the Widom line to the hot wall, significant asymmetries are observed when comparing regions near the cold wall and near the hot wall. Previous transformations that attempt to collapse the near-wall mean temperature profiles among different cases to a single curve are examined. By addressing model deficiencies of these transformations, a formulation for an improved mean temperature transformation is proposed, with appropriate considerations for real fluid effects that involve strong variations in thermodynamic quantities. Our proposed transformation is shown to perform well in collapsing the slope of the logarithmic region to a single universal value with reduced uncertainty. Coupled with a predictive framework to estimate the non-universal shift parameter of the logarithmic region using a priori information, our transformation provides an analytic profile to model the near-wall mean temperature. These results thus provide a framework to guide the development of models for wall-bounded transcritical turbulence.


Introduction
Fluids exceeding the critical pressure, p c , are ubiquitous in many environmental and engineering settings. Well-known natural occurrences exist in the atmosphere of Venus (Lebonnois & Schubert 2017) as well as subsurface flows near terrestrial undersea volcanoes (Parisio et al. 2019). Technical applications include fluids used in chemical impregnation and extraction processes, chromatography, and polymer processing (Knez et al. 2014). Especially in many industrial applications, high-pressure fluids are also subject to significant heat transfer, with convective heat transfer being a dominant mechanism for energy transport. Examples of such applications include power generation systems and heat engines Yu et al. 2019). Additionally, the operational designs of many energy conversion systems under development are predicated on efficient heat transfer with working fluids at supercritical pressures, including regenerative cooling in rocket engines (Ruan & Meng 2012), nuclear reactor coolant systems (Yoo 2013), and supercritical water oxidation (Pizzarelli 2018).
For the large temperature ranges found in many technical and environmental applications, many such high-pressure fluids exist at transcritical conditions. At these conditions, the pressure p exceeds p c but the temperature T straddles the critical temperature T c . In transcritical fluids, a number of non-ideal effects have been observed when the reduced pressure p r ≡ p/p c and the reduced temperature T r ≡ T/T c approach unity. Transitioning from the liquid-like to gas-like regions, the density, viscosity and diffusivity display sharp gradients, while the isothermal compressibility, specific heat capacity and thermal conductivity exhibit pronounced peaks (Clifford & Williams 2000;Kiran, Debenedetti & Peters 2012). These behaviours are illustrated in figure 1, which shows thermoviscous properties -density ρ, constant-pressure specific heat capacity c p , thermal conductivity λ, and dynamic viscosity μ -of nitrogen for different supercritical pressures.
At turbulent conditions in the transcritical regime, the presence of both liquid-like and gas-like thermodynamic property dependencies yields two spatial regions in the fluid domain, each with distinctive behaviour as influenced by the thermodynamics. Compared with classical constant-property trends, observations of turbulence in the gas-like region have shown -among other features -increased overall streamwise and spanwise anisotropy, enlarged spanwise spacing of large-scale near-wall structures , and intensified relative turbulent fluctuation levels of thermodynamic transport properties (Kim, Hickey & Scalo 2019). Reverse trends have been observed in the liquid-like regime. In addition, at the transition between liquid-like and gas-like conditions, particularly abrupt variations in thermodynamic properties induce sharp increases in turbulent fluctuations to further alter the flow dynamics. Overall, the coupling between the transcritical property variations and the turbulence dynamics yields behaviours that contrast with those of the more extensively-studied compressible turbulent wall-bounded flows that are well-described as ideal gases.
Additionally, in realistic transcritical flows, these coupled effects of modulated turbulence and property variations can be combined with a number of other physical phenomena, including thermoacoustic oscillations and buoyancy, all of which contribute towards further modification of the resultant turbulence dynamics and heat transfer, as discussed in Kim et al. (2019) and Nemati et al. (2015). The presence of buoyancy effects can cause transcritical fluids to undergo significant reduction in the magnitude of heat transfer. Known as heat transfer deterioration, this is a phenomenon that is known to arise from reduced turbulence production and decreased near-wall fluid density (Pizzarelli 2018), but which current analytical models and correlations have not been able Thermal field in transcritical turbulent channel flows to successfully predict (Yoo 2013). In recent decades, while the momentum behaviour of turbulence has been examined thoroughly for both variable-property and constant-property wall-bounded flows, the detailed characterization of the thermal transport is not nearly as well understood (Patel, Boersma & Pecnik 2017). Further physical study and understanding are necessary to predict accurately the temperature statistics and associated heat transfer in transcritical flows.
The complex interplay of the physical processes that are present in turbulent transcritical flows poses a number of challenges when investigating the heat transfer. Experimental measurements at such conditions lack spatial and temporal resolution; the available data are mostly limited to averaged heat transfer quantities with minimal information on turbulence statistics (Ma, Yang & Ihme 2018). Presently, models used in Reynolds-averaged Navier-Stokes (RANS) methods for transcritical fluids are unable to predict accurately flow statistics of interest, particularly for heat transfer values (Yoo 2013). Current wall-modelled large eddy simulation (WMLES) techniques have been developed largely to simulate ideal gas flows (Ma et al. 2018) and thus are not suited for calculations in transcritical conditions. With such constraints, current and ongoing investigations of turbulence in transcritical fluids are predominantly direct numerical simulation (DNS) studies that resolve the full range of turbulent scales in the flow.
Efforts to characterize the thermal boundary layer have led to the development of the temperature law of the wall as a thermal analogy of the well-known velocity law of the wall. The temperature law of the wall, first proposed by von Kárman (1939), is a functional relation for the mean temperature differenceθ = |T − T w |, where T is the fluid temperature, T w denotes the temperature at the wall, and· indicates averaging over homogeneous and steady-state conditions. A general analytical form of the temperature law of the wall, which can be obtained through scaling arguments, is expressed as (Kader 1981) where Pr = c p μ/λ is the averaged Prandtl number, κ T is a constant, and β(Pr) is a function of Pr; ζ is a modified wall-normal coordinate that provides the distance from the pertinent wall. The superscript '+' on ζ andθ indicates quantities measured in wall units, via normalization by friction quantities as where τ w ≡ |μ ∂u/∂y| w is the wall shear stress, u τ ≡ √ τ w /ρ w is the friction velocity, q w ≡ |λ ∂T/∂y| w is the wall heat flux, ν w ≡ μ w /ρ w is the kinematic viscosity evaluated at the wall, and T τ ≡ q w /(ρ w c p,w u τ ) is the friction temperature. The subscript w indicates variables evaluated at the wall, or calculated using averaged quantities that are then evaluated at the wall when applicable. Here, y is a global wall-normal coordinate. For historical context, the approach of using self-similarity methods to obtain the temperature law of the wall closely resembles the development of the Monin-Obukhov (MO) similarity theory (Monin & Obukhov 1954). In the surface layer of the atmospheric boundary layer, the MO similarity theory utilizes the near-constant vertical fluxes to describe successfully the wind speed and temperature profiles in stable boundary layers. Even with considerations for buoyancy-driven physics, the MO similarity theory also yields a log-linear profile that remains similar in form to (1.1) (Stull 1988).
The viscous sublayer solution of (1.1) is well characterized for a broad range of flows. In contrast, there is currently no consensus regarding a consistent representation of the logarithmic region. The most successful and widely-used correlation was developed by Kader (1981), who employed an assumption of constant turbulent Prandtl number in the logarithmic region and arrived at the expressions κ T = 0.4717 and β(Pr) = (3.85Pr 1/3 − 1.3) 2 + ln(Pr)/κ T . In flows with spatially homogeneous Pr, Kader's correlation agrees well with mean temperature profiles from experimental data for 0.7 ≤ Pr ≤ 170, for data from fully developed turbulent flows in flat-plate boundary layers (Simonich & Bradshaw 1978;Blair 1983) and in channels (Kader 1981). Agreement with numerical data has also been observed in DNS (Kim & Moin 1989;Pirozzoli, Bernardini & Orlandi 2016;Abe & Antonia 2017) as well as in RANS and large eddy simulation (LES) settings (Duponcheel et al. 2014;van Cauwenberge et al. 2015). However, characterizing the mean temperature profile using Kader's correlation in other flow set-ups has shown a number of deficiencies, specifically in characterizing the logarithmic region quantitatively. This poor agreement is seen in flows with strong temperature gradients (Toutant & Bataille 2013), in chemically reacting flows (Artal & Nicoud 2006), and in flows with strongly varying thermodynamic properties (Patel et al. 2017;Guo, Yang & Ihme 2018).
In the following, we review the relevant body of literature for variable-property wall-bounded DNS studies that attempt to characterize the mean temperature profile. We introduce here the density ratio, Ω, defined as Ω ≡ ρ cold /ρ hot , with ρ cold being the density at the cold wall and ρ hot being the density at the hot wall. A number of early numerical simulations of variable-property flows have considered relatively small Ω values of O(1) and simplified thermodynamics wherein one or more transport properties are either held 934 A45-4 constant artificially or prescribed by a simplified algebraic dependence on temperature. While general insight into the behaviour of the temperature profile in variable-property turbulence was accessible from these early studies, limited understanding is available into the more extreme conditions that are representative of typical transcritical turbulence in real-world applications (with Ω values of O(10-100)). Nicoud (1998) and Nicoud & Poinsot (1999) performed simulations of variable-property isothermal wall channels with Ω = 2 and showed significant deviation from Kader's mean temperature correlation. Later, Patel et al. (2017) analysed a set of nine variable-density turbulent channel flows with Ω 2.5, and showed that both Kader's correlation and a proposed modification developed through the analysis by Lee et al. (2014) for flat-plate turbulent boundary layers have limited success in representing the mean temperature accurately. Instead, Patel et al. (2017) proposed an 'extended van Driest transformation' for the mean temperature profile, demonstrating good collapse across cases with different behaviours of the semi-local friction Reynolds number Re * τ =ρ √ τ w /ρ L y /μ but constant Pr. Here, L y is the channel half-height. Recent computational investigations have utilized more realistic thermodynamic models in an effort to examine conditions that are more representative of real-world transcritical turbulence. Toki, Teramoto & Okamoto (2020) presented results for turbulent transcritical channel flow calculations at density ratios Ω = 4 and Ω = 8 before proposing a mixing-length model to transform the mean temperature profile. Comparisons of results from their transformation with incompressible DNS data by Morinishi, Tamano & Nakamura (2007) appeared successful, conditional on allowing the amount of shift in the logarithmic region to be fitted empirically a posteriori. Wan et al. (2020) presented channel flow simulations of transcritical carbon dioxide with a maximum density ratio Ω ≈ 3. They then extended the capabilities of Patel et al. (2017)'s 'extended van Driest transformation' to include the effects of variable c p . These more recent studies demonstrate relative success in the characterization of the mean temperature profile and other temperature statistics, but still at limited ranges of density ratios and heating conditions examined. At the transcritical conditions relevant to practical applications, current modelling efforts have not yet demonstrated confidence in the ability to characterize or predict temperature statistics.
To address this issue, this investigation provides a detailed analysis of DNS calculations of fully developed transcritical turbulent channel flow, with a focus on the analysis of the thermal boundary layer. Our calculations involve (1) a compressible Navier-Stokes formulation, (2) realistic thermodynamic considerations to represent accurately the behaviour of state variables and transport properties in the transcritical regime, and (3) density ratios Ω of up to O(20) to better represent the behaviour of transcritical turbulence in real-world settings. The main objective of this study is to provide a theoretical framework to best characterize the mean temperature profile in the logarithmic region. As a secondary objective, we provide observations of the turbulent statistics and structures of the thermal flow field, and how they are affected by the operating conditions.
The remainder of this paper is organized as follows. Section 2 presents the problem formulation, which includes details of the governing equations, model relations utilized, domain definition and overall computational set-up. Section 3 presents results from the simulations along with associated interpretations, discussions, and modelling efforts. Finally, § 4 offers conclusions.

Governing equations, thermodynamic models, and domain definition
The governing equations that are solved are the conservation equations of mass, momentum and total energy: where u i is the ith component of the velocity vector, p is the pressure, and e t is the total specific energy, combining the specific internal energy e and the specific kinetic energy (u i u i )/2. Here, f i is the ith component of the body force vector, which imposes a prescribed bulk streamwise momentum. This body force vector serves as a streamwise-homogeneous substitute for a pressure gradient. For the Reynolds numbers in the current investigation, the choice of a body force to replace a conventional pressure gradient has a negligible effect on the pressure drop from skin friction. Previous authors (Huang, Coleman & Bradshaw 1995;He, He & Seddighi 2016;Quadrio, Frohnapfel & Hasegawa 2016) have shown additionally that characteristic turbulence statistics are largely robust to the choice of forcing scheme. Following Huang et al. (1995), f i is written as where the direction i = 1 indicates the streamwise direction, the subscript '0' denotes a volume-averaged quantity, and δ ij is the Kronecker delta operator. Here, τ ij is the viscous stress tensor defined as and q i is the heat flux vector given as where λ is the thermal conductivity. The working fluid considered is nitrogen (N 2 ) with critical pressure p c = 3.3958 MPa, critical temperature T c = 126.19 K, molecular weight W = 28.0134 g mol −1 , and acentric factor ω = 0.0372. We close the system of equations with the Peng-Robinson (PR) equation of state (EoS) (Peng & Robinson 1976;Poling, Prausnitz & O'Connell 2001), which has been used widely in supercritical and transcritical settings because of its accuracy in predicting thermodynamic state variables in the vicinity of the critical point (Miller, Harstad & Bellan 2001). At transcritical conditions comparable to those used in our investigation, Peng & Robinson (1976) reports root-mean-square (r.m.s.) errors below 1.25 % in enthalpy departure prediction when compared to experimental values, while Congiunti & Bruno Thermal field in transcritical turbulent channel flows (2003) and Hickey et al. (2013) show results with similar error levels when representing ρ and c p . The PR EoS is expressed as where R is the gas constant. The parameters a and b account for real fluid effects and are given as and with ω being the previously introduced acentric factor. For N 2 , b = 8.58 × 10 −4 and c = 0.432. Chung's model for high-pressure fluids (Chung, Lee & Starling 1984;Chung et al. 1988) is used to evaluate molecular transport properties. For N 2 at temperatures and pressures comparable to those used in our current investigation, Chung's model demonstrates average absolute deviations below 1.24% for μ and 7.30% for λ when compared to experimental data; both values lie within typical experimental uncertainty ranges. With its relative accuracy and computational efficiency, this model has also been used by a number of previous studies of transcritical turbulence (Ruiz 2012;Toki et al. 2020).
A schematic of the computational domain is given in figure 2. The channel is periodic in streamwise and spanwise directions, while the wall-normal coordinate y extends from −L y to L y . The cold wall for each case is thus at y/L y = −1, and the hot wall is at y/L y = 1. As introduced in § 1 and also depicted in figure 2, the coordinate ζ provides the wall-normal distance from each pertinent wall and satisfies ζ = 0 at the wall being considered. Gravity is not considered in our calculations, so the relative orientation of the hot and cold walls is arbitrary. The domain dimensions are L x × 2L y × L z , with L x /L y = 2π, L z /L y = 4π/3, and the channel height measuring 2L y = 9.0132 × 10 −5 m. In their comprehensive study, Lozano-Durán & Jiménez (2014) showed that dimensions of L x = 2πL y and L z = πL y are sufficient for capturing one-point statistics in channel flows with Re τ up to 2009 for isothermal ideal gas flows. Further investigation is needed to extend this conclusion to real-fluid turbulent flows. However, with the grid and domain validation results presented in Appendix A and in Ma et al. (2018), we conclude that the domain dimensions in our study are adequate for capturing the desired statistics of our study.

Computational set-up
In our simulations, we used a compressible finite-volume solver (Khalighi et al. 2011;Hickey et al. 2013;Larsson et al. 2015;Ma et al. 2018). The governing equations are solved in dimensional units using a strong stability-preserving Runge-Kutta scheme with third-order accuracy for time stepping, and a flux reconstruction central scheme with fourth-order accuracy on uniform grids and third-order accuracy on non-uniform grids for spatial discretization. Ensuring that the effects of numerical errors are minimal for the current choice of numerical methods at the high density ratios studied is important for 934 A45-7  Figure 2. Schematic for the turbulent channel at transcritical conditions. Note that x denotes the streamwise coordinate (with velocity component u 1 = u), y denotes the wall-normal coordinate (with velocity component u 2 = v), and z denotes the spanwise coordinate (with velocity component u 3 = w). The wall-distance coordinate ζ is also depicted for each wall. The hot and cold walls are each independently assigned isothermal temperature values at T hot and T cold , respectively. confidence in the results. To this end, for the currently-used numerical procedure, Ma, Lv & Ihme (2017) and Ma et al. (2019) demonstrated minimal dispersion errors via solution profiles that are robust against the formation of spurious numerical oscillations, as well as evidence for the negligible effect of dissipation errors through convergence studies. In our current study, we also provide evidence towards the robustness of the numerical procedure against dissipation errors using the well-resolved energy spectra provided in Appendix A.
For our simulations, a relative solution (RS) sensor has been applied (Ma et al. 2017); in regions where the local density ratio, defined using the reconstructed face value and the control volume centre value, exceeds 50 %, a second-order hybrid essentially non-oscillatory (ENO) scheme is employed along with a Harten-Lax-van Leer contact (HLLC) Riemann flux computation. The same framework is also applied for the local pressure ratio. An entropy-stable double-flux model, developed as part of the procedure to simulate transcritical flows with physically-realizable solutions, is used (Ma et al. 2017).
We use a structured Cartesian grid with mesh size N x × N y × N z = 384 × 256 × 384. The streamwise and spanwise spacings are uniform, while the wall-normal spacing is stretched using a hyperbolic tangent profile to resolve the viscous sublayer. Validation of the current numerical procedures for stretched grids is provided in Ma et al. (2017Ma et al. ( , 2018. Details of the grid resolution are provided in Appendix A. Details of individual cases considered are provided in table 1, and a thermodynamic state diagram is provided in figure 3. The operating conditions were chosen to sample density ratios Ω for values from 1 to ∼20. To study different heating conditions, we vary the temperature of the top wall T hot among different cases. Cases are named based on the ratio of the hot wall to the cold wall temperatures, which we denote as the temperature ratio (TR). Of note, cases TR3, TR1.9, TR1.4 and TR1.3 are transcritical, whereas cases TR1.25 and TR1 are strictly subcritical in temperature. The wall temperatures of case TR1 are the same, and the case is included as a reference. The maximum Mach number, calculated usingū and the mean speed of sound, is less than 0.16 for all cases. For all cases, the bulk pressure is set to a constant value 3.87 MPa, which corresponds to a reduced pressure p r = 1.14, and the bulk Reynolds number is set to a constant value Re 0 = 2ρ 0 u 0 L y /μ 0 = 3.5 × 10 4 . At the given bulk pressure, the transition temperature of the Widom line (WL), T WL , is found at a reduced temperature T WL /T c ≈ 1.022; this  (see table 1 for parameters), with the circle denoting the condition at the cold wall and ×-marks denoting conditions at the hot wall. All cases are at p r,0 = 1.14.
Cases  Table 1. Summary of cases and conditions, where T r,cold = T cold /T c , T r,hot = T hot /T c , and ρ r,0 = ρ 0 /ρ c is the reduced volume-averaged density, ρ c is the critical density.
is the temperature corresponding to maximum c p at the given supercritical pressure. For each case, two different friction Reynolds numbers are reported at each of the two walls, calculated using property values evaluated at the conditions of the appropriate wall. For the length scale, the friction Reynolds number Re τ uses the channel half-height L y . In contrast, the modified friction Reynolds number Re τ,max(ū) uses the distance of maximum u from the wall, L y,max(ū) . Because the bulk pressure p 0 is held constant while the bulk temperature T 0 changes among cases as a result of the adjusted temperature boundary conditions, the bulk density ρ 0 of the system needed to be adjusted accordingly for the flow to remain statistically stationary. For each case, the values for the bulk momentum (ρu) 0 and bulk density ρ 0 were found by iterating. After reaching a statistically stationary state, each case is run for more than six flow-through times, where each flow-through time is defined as L x /u 0 .
Values of the Eckert number, calculated as are provided in table 2. Here, h is the specific enthalpy, defined as h = e + p/ρ, h w denotes the mean enthalpy at the wall, andh CL is the mean enthalpy at the channel centreline (CL). As observed, viscous heating and dissipation effects become more significant as the temperature ratio decreases. Cases TR3, TR1.9, TR1.4 and TR1.3, being transcritical, all cross the Widom line within the fluid domain and thus contain both liquid-like and gas-like fluid behaviours. From past studies regarding the coupled effects of property variations and turbulence modulation (as discussed in § 1), we expect many of the previously-observed features to be present in these transcritical cases; these details will be examined and discussed in § 3.
Although case TR1.25 is subcritical in temperature, the relative proximity of the hot wall temperature (125 K) to the Widom line transition temperature (∼129 K) means many of the discussed effects that are intrinsic to transcritical flows are still observable. Case TR1, on the other hand, has wall temperatures (both at 100 K) markedly far from the Widom line transition temperature; we expect the behaviour in this case to be similar to that of constant-property turbulence.

Results
In this section, we present DNS results for all cases, along with associated analysis and modelling. Note that investigation of the momentum characteristics of case TR3 were presented in Ma et al. (2018). Subsection 3.1 provides observations regarding instantaneous contours and fluctuations profiles, and § 3.2 discusses mean profiles. These two subsections provide the foundation for characterization of the near-wall mean temperature behaviour. Consequently, § 3.3 discusses the performance of previous mean temperature transformations and models, before we provide the mathematical formulation and results for our improved mean temperature transformation and model in § 3.4.

Instantaneous contours and fluctuation profiles
As a means of visualizing the instantaneous thermal field, figure 4 shows contours of c p normalized by the cold wall value, which we denote c p,cold , for the four transcritical cases. As the temperature ratio decreases, we observe a shift in location of the transition point of the Widom line towards the hot wall. Increased intensity of fluctuations in c p is visible for cases with larger temperature ratio.
For observations of how variations in temperature boundary conditions affect the resultant temperature field, figures 5 and 6 show temperature fluctuations near the cold and hot walls. Two different wall-normal distances are shown to capture the behaviour in both the viscous sublayer (with planes at ζ * = 5) and the logarithmic region (with planes 934 A45-10 Downloaded from https://www.cambridge.org/core. Thermal field in transcritical turbulent channel flows  at ζ * = 100). Here, ζ * is the semi-local coordinate and is defined as ζ * ≡ ζ √ τ wρ /μ. As a modification of ζ + , ζ * has been used extensively in variable-property wall-bounded turbulence in order to account for the effect of local property variations on the magnitude of near-wall characteristic turbulent scales (Huang et al. 1995).
In the viscous sublayer, the structure of the streaks is visibly distinct between the cold and hot walls. Similar qualitative characteristics in temperature fluctuations were also reported by Sengupta et al. (2017). When values are normalized by the friction temperature T τ , we observe an increase in the relative intensity of temperature fluctuations near the hot wall at both reported ζ * distances for case TR1.3, due to the proximity of the hot wall temperature to T WL .
There exists a sizeable body of literature studying near-wall velocity structures (especially for u ), from which observations and associated theoretical predictions have 934 A45-11 Downloaded from https://www.cambridge.org/core. been made. A key conclusion from previous studies is the spanwise spacing of 100 wall plus units for near-wall streaks in u (Kline et al. 1967;Kim, Moin & Moser 1987).
Observations of the streak structures in temperature in figure 5 show similar spacing values of ∼100 when measured using semi-local units; these visual observations are confirmed using two-point spanwise correlations of temperature fluctuations, as plotted in figure 7. Departing from this trend is the hot wall of case TR3, which is characterized by a larger spacing approaching 400 semi-local units. Deviations from the value of 100 are consistent with the conclusions reached by Patel et al. (2015); the liquid-like structures at the cold wall tend to exhibit decreased spacing, while the gas-like structures at the hot wall tend to have increased spacing. We note also that with the high degree of correlation observed between u and T in the near-wall regime (Guo et al. 2018), many conclusions regarding structures in T can be compared directly with those from the near-wall u literature. Further details of transcritical velocity structures and statistics, including evidence supporting the presence of attached eddies and other comparisons with results from classical constant-property wall-bounded turbulence, are provided in Ma et al. (2018). For a more quantitative view of the apparent asymmetry in temperature fields, figure 8 displays probability density functions (p.d.f.s) of temperature in the viscous sublayer and in the logarithmic region. In the cold wall profiles, we observe pronounced skewness towards the hot wall temperature in the viscous sublayer. This behaviour is an indication of significant mixing of thermal energy between the two walls; similar behaviour in the density p.d.f. distributions was reported by Ma et al. (2018). Profiles in the logarithmic 934 A45-12 Downloaded from https://www.cambridge.org/core. region are less skewed (consistent with increasing isotropy in flow structures closer to the channel centre) and also display decreased overall variance compared to the viscous sublayer profiles. For the hot wall profiles, a skewness is again observed in the viscous sublayer towards the temperature of the opposite (cold) wall. However, distinct from the behaviour of the cold wall profiles, many of the hot wall profiles are observed to cluster near T r = 1 as a consequence of the proximity of T c to T WL .
934 A45-13 To examine the importance in fluctuations of thermodynamic quantities, figure 9 shows profiles for relative magnitudes of r.m.s. fluctuations in ρ, c p , μ and λ. Just as was observed for the instantaneous snapshots and the temperature p.d.f.s, a pronounced asymmetry between cold and hot wall values is observed in all cases except for case TR1, with much larger fluctuation magnitudes near the hot wall. All asymmetric cases with temperature ratio > 1 are also characterized by r.m.s. fluctuation values that are significantly larger than for the symmetric case TR1.
Especially for cases with large temperature ratio (cases TR3 and TR1.9), the fluctuations in c p are most significant and consistently exceed 50 % of the local mean value. This observation suggests the importance of capturing the behaviour of c p on characterizing the thermal boundary layer, a conclusion that was also noted by previous authors (Wan et al. 2020) and which we will utilize in § 3.4. We note that fluctuations in ρ, μ and λ exceed 20 % of the mean value near the hot wall -with highly localized peak values for cases TR1.3 and TR1.4 -and are by no means negligible. Notably, cases TR3 and TR1.9 have fluctuation magnitudes 10% of the corresponding mean value across the majority of the channel for each of the four plotted quantities. Observations of these fluctuation profiles have important implications for the temperature characterization. With sufficiently large temperature ratio, fluctuations of all thermodynamic properties in transcritical turbulence cannot be neglected presumptively and require consideration in boundary layer modelling. This conclusion is in direct contrast to classical results. Constant-property turbulence by its nature has negligible thermodynamic fluctuation levels, and supersonic turbulent boundary layers (with significant density variation) generally follow Morkovin's hypothesis and are characterized by fluctuation magnitudes typically not exceeding 10-15 % of the mean value (Coleman, Kim & Moser 1995;Huang et al. 1995). Favre-averaged mean. The molecular heat fluxes exhibit similarity in behaviour among all cases. However, the same similarity is not observed for the turbulent heat fluxes, including in the near-wall thermal boundary layer (regions with large gradients in mean temperature, as can be observed in profiles of the mean temperature provided in figure 25a). The significant dissimilarity in turbulent heat flux profiles seen here was not observed in cases with lower temperature ratio values (Toki et al. 2020). For the momentum transport, figure 11 displays the decomposition of the total mean stress into the viscous stress μ ∂u/∂y and the Reynolds stress −ρ u v . In all cases, the profiles of viscous stress collapse well near the cold wall, while the asymmetry among different cases is evident moving towards the hot wall. It can be seen that as the temperature ratio decreases, the stress profile becomes more and more symmetric, with the Reynolds stress profile equalling zero closer to the channel centreline, and τ w,hot increasing in magnitude to approach that of τ w,cold . Overall, similarity in the behaviour among Reynolds stress profiles is evident -a feature not observed in the turbulent heat flux profiles in figure 10.

Mean profiles
Despite the differences between the behaviours of the turbulent heat fluxes and Reynolds stresses, similarity is observed for the ratio of the turbulent and molecular components, as done in figure 12. The profiles provide a metric for the relative importance of each component for the momentum and energy transport as a function of ζ * . As observed, the dynamics of both the momentum and thermal fields are governed primarily by the molecular component for ζ * 5 and by the turbulent component for ζ * 30. The presence of distinct spatial regimes and transport mechanisms supports early theory that led to the recognition of the logarithmic scaling of the velocity distribution, thus providing key insights towards the development of Townsend's hypothesis of equilibrium layers (Townsend 1961) and the Monin-Obukhov similarity theory (Monin & Obukhov 1954). The observations made here, in addition to justifying the choice of ζ * used in studying the near-wall temperature fluctuations in § 3.1, will also be utilized for defining the extent of the temperature logarithmic region in § 3.3. Before considering the near-wall mean temperature behaviour, it is informative to gain insight into the scales that govern the near-wall dynamics. The friction velocity and temperature are plotted in figure 13(a). From observations of u τ across both cold and hot walls, and in all conditions considered, variations by a factor of less than 2 are seen. The variability of the cold wall T τ is somewhat larger, with overall variation of a factor of approximately 6, while the hot wall T τ displays much larger variations of more than two orders of magnitude. We observe that for the hot walls, trends in both u τ and T τ have a noticeable non-monotonicity for cases with T hot ≈ T c . This observation is rationalized by 934 A45-16  Figure 13. (a) Plots of friction velocities u τ and friction temperatures T τ , as a function of reduced hot wall temperature T r,hot . Normalization in each line uses the values for case TR1. (b) Plots of magnitude of mean heat flux q w = |λ ∂T/∂y| w at hot and cold walls, plotted as a function of the reduced hot wall temperature T r,hot . Normalization by the value found in case TR1 is used. In both plots, T WL is also plotted as a vertical dotted line. the significant gradients in thermodynamic properties near T WL . Figure 13(b) compares the magnitude of the mean wall heat flux q w = |λ ∂T/∂y| w , for each case and for the cold and hot walls. The overall trend of decreasing heat fluxes at both walls for smaller values of T hot is expected from observations of the mean temperature profiles in figure 25. Deviating from this trend, we observe a considerable increase in q w,hot for T hot /T c = 1.03 (corresponding to case TR1.3) that can be explained by the peak in λ that occurs at T WL for weakly supercritical pressures (as observed in figure 1). Knowledge of the variation in characteristic scales given by figure 13 provides a foundation to investigate the temperature behaviour. Figure 14(a) shows the near-wall mean temperature profiles, normalized using the friction temperature T τ . The cases each appear to display a region with linear slope as indication of the logarithmic region; this observation is also substantiated using the diagnostic function plotted in figure 14(b). If friction quantities by themselves were sufficient to describe the mean temperature dynamics, then all profiles would collapse into a single curve and thus satisfy a universal temperature law of the wall, in the form as presented by (1.1) but with ζ * replacing ζ + . While this is not the case, the observed profiles provide insight into the ability of T τ 934 A45-17 Downloaded from https://www.cambridge.org/core.  Figure 14. Profiles for (a) normalized mean temperatureθ + , and (b) the log law diagnostic function ζ * dθ + /dζ * . For each profile, the diagnostic function approaches a constant value in the logarithmic region, corresponding to 1/κ T in (1.1). Dashed lines indicate cold wall profiles, and solid lines indicate hot wall profiles. All 12 curves (6 cases with 2 walls each) are displayed. Profiles are plotted versus the semi-local coordinate ζ * .
to describe accurately the near-wall scales for the temperature dynamics. The cold wall profiles collapse more tightly -with similar log region slopes and amount of vertical shift -than the hot wall profiles do; this indicates that T τ at the cold walls performs relatively well as a normalization quantity. This good collapse of the cold wall profiles is expected; profiles at the cold walls are more similar in nature to classical constant-property results, with comparatively less effect of near-wall property gradients on the mean flow. In contrast, the hot wall profiles of cases TR3 and TR1.9 display noticeably smaller slopes and are shifted vertically downwards compared to the cold wall profiles and to all other hot wall profiles. This suggests that near the hot wall, T τ,TR3 and T τ,TR1.9 overestimate the true near-wall temperature scales. The opposite conclusion can be made for the hot wall profiles of TR1.3 and TR1.25; the log regions in these cases have slope values that are too large and profiles that are shifted vertically upwards. These profiles of near-wall mean temperature serve as a baseline for the characterization efforts that follow. Additional mean profiles to illustrate the flow are provided in Appendix B.

Assessment of previous mean temperature theoretical work
As discussed in § 1, the mean temperature behaviour of the viscous sublayer is well characterized. For completeness, we show these profiles in § 3.2. In this subsection, we focus on the quantitative assessment of the performance of previous transformations from literature in characterizing the logarithmic region of the mean temperature profile for our simulations. Note that while we included case TR1 in plots and discussions up to this point as a reference case, we will not consider it for purposes of characterization and modelling.
Before beginning the assessment, the boundaries of the temperature log region must be defined properly. For the velocity log law, the log region is defined as being between ζ * > 30 and ζ /L y < 0.3; these bounds have been well-validated using data from both computational and experimental investigations (Pope 2000).
The lower bound value of ζ * = 30 for the velocity log law is chosen based on observations from past studies that the direct effect of the molecular viscosity on the total viscous stress is negligible in the region further than ∼30 viscous length scales from the wall (Pope 2000). From the observations and discussion associated with figure 12, it is seen that both the viscous stress and the molecular heat flux are negligible for our cases 934 A45-18  Table 3. Distance from wall to location of maximum mean streamwise velocity, L y,max(ū) /L y , and the corresponding location of the upper bound of the temperature log region, ζ = 0.3 L y,max(ū) , in semi-local units.
for ζ * 30. Thus we choose to use the same value of ζ * = 30 for the lower bound of the temperature log region. Observations of figure 14 support this choice for the cases under investigation, with a region of linear slope for all cases beginning at ζ * ≈ 30 and extending away from the wall.
For the upper bound of the temperature log region, modification of the velocity log law criterion (ζ /L y < 0.3) is necessary. The criterion from the velocity log law stems from the fact that ζ = L y is the location of the maximum ofū in a symmetric channel. In our current investigation, it can be seen in figure 25(c) that the peak inū deviates significantly from ζ = L y for several cases, making ζ = L y an ill-suited choice in our current investigation. Instead, we choose L y,max(ū) -the distance from the wall to the location of maximumūas the relevant length scale for each case, as was also used to calculate Re τ,max(ū) in table 1. The modified upper bound for the temperature log region is thus ζ /L y,max(ū) < 0.3. Values for the location of this temperature log region upper bound are provided in table 3.
With our clarification of the boundaries for the temperature log region, we proceed by reviewing previous transformations from literature. We quantify the performance of a transformation by (1) how the slopes of the temperature log region differ among cases, and (2) how well the shift parameter of the temperature log region is characterized. Here, the shift parameter of the temperature log region is defined to be the value of the transformed temperature at ζ * = 30, and can be interpreted as a measure of the vertical variation in the log region of the transformed temperature.
Note that the presented integrals are of the form   of the temperature profiles, and each study attributes this result to the influence of near-wall property variations that are unaccounted for by the transformation. We plot the results of this transformation for our data in figure 15(a) and note that, in addition to the variation in the log region shift parameter, the slopes among different cases vary across a wide range. It is evident that accounting only for variations inρ and c p is not sufficient to fully describe the mean temperature profile. In contrast, Ma et al. (2018) demonstrated good agreement with the established velocity log law correlation when using the velocity van Driest transformation on theū profile.

A45-20
Thermal field in transcritical turbulent channel flows A second transformation proposed by Toki et al. (2020) is derived using the mixing length hypothesis for both momentum and energy and also by assuming linear shear stress and heat flux. This transformation is written as The results for this transformation are provided in figure 15(b). The Toki transformation shows a slight improvement over the van Driest transformation when observing the collapse of the slope. However, significant deviations from the reported slope value of 2.8 (provided in Toki et al. 2020) are seen, especially for the hot wall of case TR3. We now analyse the behaviour of the extended van Driest (eVD) transformation, first proposed by Patel et al. (2017) and later extended to variable c p conditions by Wan et al. (2020). The derivation of the transformation begins from the low-Mach energy equation, with a final expression written as where the turbulent eddy conductivity is defined to be α t = (−ρ v θ )/(dθ/dy). The local Prandtl number is defined as Pr * = Pr w (c p /c p,w )(μ/μ w )/(λ/λ w ) = c pμ /λ and is plotted and discussed in Appendix B. Notably, while (3.2) and (3.3) are direct integral transformations for temperature, (3.4) is a spatial integral. The results of this transformation for our data are provided in figure 15(c). Though significant improvement in the collapse of the slope is achieved when compared to previous models, we observe a large spread in slope values of ∼ ± 0.8 (corresponding to an uncertainty percentage of approximately ±25 % from an average value of ∼3.3). Increasingly large slope values are found for the hot wall of cases with larger temperature ratios, with an especially large value for the hot wall of case TR3. This difficulty in collapse indicates that the governing physical processes of the near-wall temperature dynamics are not sufficiently accounted for at higher temperature ratios. Our current study seeks to address this issue in § 3.4. An additional feature of the extended van Driest transformation is a mathematical description of the log region shift parameter by decomposing (3.4) into two portions: One portion,θ * T ,eVD , sets Pr * = 1 to remove variations in the viscous sublayer (and thus remove variation in the log region shift parameter) to isolate the log region slope. A second portion,θ * P,eVD , captures the different log region shift parameters of each case by removing the effect of the log region slope. Figure 16 shows this decomposition. With this formulation, the log region shift parameters given in figure 15(c) can be found by adding the asymptotic values ofθ * P,eVD to the log region shift parameter ofθ * T ,eVD . We expand upon this decomposition in our model development in § 3.4.
In these presented transformations, we observe significant non-universality in the value of log region slopes obtained. Difficulties in collapsing the slope values are especially pronounced near the hot wall for cases TR3 and TR1.9 (with T r,hot = 2.58 and 1.51, respectively); these coincide with large values of T τ and q w , as shown in § 3.2 through the 934 A45-21 Downloaded from https://www.cambridge.org/core. discussion associated with figure 13. We observe deviations from the general trend of hot wall slopes for cases TR1.3 and TR1.25 as a result of the proximity of the transition point of the Widom line; however, we do not observe the same degree of difficulty in collapse of slopes as in the hot walls of cases TR3 and TR1.9. Specifically for the log region shift parameter, a common theme that we observe is the larger shift parameter values found near the hot walls for cases TR1.3 and TR1.4 (with T r,hot = 1.03 and 1.11, respectively). This common observation can be rationalized by the proximity of these cases to T WL and the corresponding peak Pr * that occurs at this temperature. We note that in the transformations presented in this section, the log region shift parameter is either not predicted by the proposed mathematical transformation (for the van Driest and Toki transformations) or characterized a posteriori using the transformed DNS profiles (for the extended van Driest transformation).

Mean temperature characterization and modelling
From observations of the performance of currently available transformations in § 3.3, an improved transformation for the near-wall mean temperature profile would provide a characterization of the log region slope with higher accuracy. An improved transformation would also contain a predictive model for the log region shift parameter. If the log region shift parameter could be determined using only values known a priori, then the details of the temperature profile in the log region can be quantitatively predicted beforehand, increasing the applicability of the transformation to turbulence modelling. These are the goals that motivate our proposed transformation, which we derive next.
To consider the entire set of relevant physical processes, we begin with the energy equation (2.1c); after subtracting the transport of kinetic energy and then averaging over time and over homogeneous directions, we obtain where q represents the total heat flux. Note that in this equation, the turbulent heat flux term −ρ v h and the molecular heat flux term λ dθ/dζ are present, as well as two additional terms: the pressure work term ζ 0 τ ij ∂u i /∂x j dϕ. These integrals are evaluated from the wall to the desired ζ location. Notably, these last two terms are neglected by the energy equation of the low-Mach limit of the Navier-Stokes equations.
In simplifying the energy equation, we make two assumptions. The first is that the fluctuating molecular heat flux is negligible, |λ dθ/dζ | |λ dθ /dζ |. Although observations of the r.m.s. fluctuations as seen in figure 9 have shown that such fluctuations of thermodynamic properties (including λ) can be significant, figure 17(a) shows that the contributions of the fluctuating heat flux are small when normalized byq. Additionally and more importantly, the profiles decay to 0 for ζ * > 30. Because our primary focus is the characterization of the temperature log region, this first assumption is justified, and we approximate the molecular heat flux term as As a second simplifying approximation, we neglect the pressure work term while keeping the viscous dissipation of kinetic energy. Figure 17(b) shows the pressure work, and figure 17(c) shows the viscous dissipation of kinetic energy. When assessing the absolute value of the relative magnitude across all cases, the pressure work is two orders of magnitude smaller than the other terms in consideration in figure 17, and thus does not provide a significant contribution to the overall heat flux. In contrast, the viscous dissipation contributes 5-10 % of the overall heat flux in certain cases. These contributions 934 A45-23 are most significant in the temperature log region for ζ * > 30. We thus neglect the fluctuating molecular heat flux and the pressure work but keep the viscous dissipation. With this simplification, we rewrite (3.6) as Normalizing each term with the wall heat flux q w = ρ w u τ c p,w T τ and rearranging leads to where the turbulent eddy conductivity α t , accounting for enthalpy fluctuations, is defined as (3.10) The potential of (3.9) in collapsing theθ profile rests on how well the quantity collapses among different cases as a function of ζ * . Here, we introduce the dθ * /dζ * label following the notation used by Patel et al. (2017). Following the definition in (3.11), we can rearrange the terms in (3.9) to obtain as an approximate representation of (3.11), which will be used to develop a mathematical transformation to model the mean temperature. Here, C 1 is given as (3.13) The profiles of dθ * /dζ * , as defined in (3.11), are plotted in figure 18(a). We observe significant spread in magnitudes of profiles among cases; in the temperature log region (defined by the discussion in § 3.3), relative magnitudes of cases span a factor of up to 3. Notably, the magnitude of profiles for the hot walls of cases TR3 and TR1.9 are consistent outliers. To correct for these observed discrepancies among profiles, we introduce two 934 A45-24 Thermal field in transcritical turbulent channel flows correction factors to be multiplied by dθ * /dζ * as given in (3.11): (3.14) with C 1 as given by (3.13) and C 2 given as As seen in figure 18(b), these two correction terms achieve significant improvement in the collapse of dθ * /dζ * , when compared with the uncorrected values in figure 18(a). The term C 1 corrects for the observed variation in total heat flux across cases. The term C 2 corrects for an implicit approximation made by the eddy conductivity definition in (3.10) -a more appropriate eddy conductivity would use the enthalpy gradient dh/dy, instead of c p dθ/dy, in the denominator of the definition of α t . Note that the specific enthalpy differential for a real fluid is a function of two independent state variables -such as T and p -and can be written as (3.16) Employing dh = T ds + (1/ρ) dp leads to where the Maxwell relation has been used. Substituting (3.17) into (3.16) leads to (3.19) thus showing that dh / = c p dθ for a real fluid. The combination of C 1 and C 2 has demonstrated their ability to provide a physics-motivated collapse of the profiles of (3.11). Our goal is to achieve an accurate transformation to represent (3.11), and in turn represent the mean temperature profile. Given the success of the correction factor procedure, we apply the same procedure to the approximate representation for dθ * /dζ * as given by (3.12) and obtain  Figure 18. Profiles of (a) dθ * /dζ * , and (b) C 1 C 2 dθ * /dζ * . Dashed lines indicate cold wall profiles, and solid lines indicate hot wall profiles. With these corrections in place, we can now define a temperature transformation in the near-wall region by integrating (3.20) to obtain Profiles of the full transformationθ * as given in (3.21) are plotted in figure 19, along with the slope and shift parameter that best fit the logarithmic region of the profile. The slopes are all clustered nicely around a mean value of ∼3.2 (with a reduced range of ∼ ± 0.56 among all cases, corresponding to an improved uncertainty percentage range of approximately ±18 %).
We observe a shift in the log region shift parameter that is seen in previous models and can be explained by Pr * effects. Borrowing the procedure of Patel et al. (2017) and Wan et al. (2020), the final transformationθ * can be decomposed into two separate terms, as Thermal field in transcritical turbulent channel flows One term,θ * T , characterizes the slope of the log region and assumes unity Pr * : The second term,θ * P , contains effects of variable Pr * on the shift parameter of the log region, and is written asθ * Figure 20(a) shows profiles ofθ * T ; with the removal of variations in log region shift parameter, we observe good collapse among profiles. To provide a modelled profile using an empirical fit, we integrate an ODE-based wall model equation for near-wall mean temperature (as an analogy of the near-wall mean velocity model equation of Kawai & Larsson 2012): where the mixing length parameter l * T is given by (3.26) Parameter choices used in figure 20(a) to best describe the data are κ * T = 2.90, A = 23 and b = 1. Figure 20(b) shows profiles ofθ * P . The observation thatθ * P reaches an asymptotic value for ζ * 30 indicates that the variation in the shift parameter of the log region forθ * comes nearly exclusively from variations in dθ * P /dζ * for ζ * 30 (in the viscous sublayer and in the buffer layer). The log region shift parameter forθ * can now be described as the sum of (1) the log region shift parameter value forθ * T , which is ∼14 for all cases, and (2) θ * P (ζ * = 30). As mentioned previously, the ability to predict the shift parameter of the log region for θ * has direct utility in turbulence modelling. In analysing the distribution of values for θ * P (ζ * = 30), we observe a significant increase for wall temperatures near T WL . We make the following observations and approximations.
(i) From the expression forθ * P given in (3.24), the value of dθ * P /dζ * at the wall (ζ * = 0) can be estimated for all cases to be Pr w − 1.
(ii) Given that Pr reaches a maximum at T = T WL , an estimated theoretical maximum value of dθ * P /dζ * for ζ * = 0 is (Pr WL − 1), where Pr WL is the value of Pr evaluated at T = T WL . A reasonable estimate for the integral of dθ * P /dζ * from ζ * = 0 to ζ * = 30 to obtainθ * P (ζ * = 30) is 3(Pr WL − 1). (iii) At values of |T w − T WL |/T WL exceeding unity, Pr will approach the same value as found in an ideal gas, which we denote as Pr ig . In this regime, a reasonable estimate for the integral of dθ * P /dζ * from ζ * = 0 to ζ * = 30 to obtainθ * P (ζ * = 30) is 5(Pr ig − 1). Figure 21, developed using these observations, shows our prediction of the values of θ * P (ζ * = 30), along with an exponential function that best fits the data. Although the exact value of the log region shift parameter would still require the evaluation of the expression in (3.24), this approximation provides a framework to estimate the value quantitatively with reasonable accuracy.
This concludes the derivation of our mean temperature transformation. Previous formulations began their model development by starting off with one or more simplifying assumptions -namely the low-Mach equations and hypotheses of linear heat flux/stressthat restrict the applicability and performance of the transformations. Instead, we consider, from the outset, all possible physical processes allowed by the energy equation.

Thermal field in transcritical turbulent channel flows
In assessment of performance, our transformation offers reduced uncertainty in the value of log region slope. Note that the current transformation improves the collapse of the slope near the hot walls for cases TR3 and TR1.9 -cases that previous transformations have struggled with, likely due to the presence of steeper gradients in mean temperature (figure 25a), larger relative values of fluctuations (figure 9), real fluid considerations for the enthalpy, and larger heat flux and friction temperature values ( figure 13). Additionally, we have presented a methodology to estimate the log region shift parameter a priori, thus providing a quantitative, predictive view for the near-wall mean temperature profile.

Conclusions
The results of DNS calculations for a fully developed transcritical turbulent channel flow at a reduced pressure p r = 1.14 are presented. The ratio of wall temperatures of the cases is varied up to a value of 3, resulting in four transcritical and two non-transcritical cases with density ratios of up to O(20). For the transcritical cases, the asymmetric thermodynamic behaviour of the liquid-like and gas-like regimes, as well as the direct influence of the Widom line, strongly interact with the turbulence dynamics and heat transfer in the near-wall regime. The impact of the transcritical thermodynamics is also evident through a number of modifications to flow structures and statistics.
(i) Compared to classical boundary layer trends, the spanwise spacing of temperature structures is widened in the gas-like regime and contracted in the liquid-like regime. (ii) Probability distributions of temperature fluctuations exhibit increased skewness, which is indicative of significant mixing between the hot and cold walls. Although profiles of the momentum and thermal transport appear disparate in nature, the molecular and turbulent transport were each found to be the dominant mechanism in the same wall-normal region of both fields.
The previous characterization efforts for the near-wall mean temperature profile are examined, along with a quantitative assessment of their strengths and weaknesses. Difficulty in consistently characterizing the near-wall temperature dynamics is seen in conditions with large density ratios and the associated steep gradients in mean temperature, larger relative fluctuation magnitudes, real fluid considerations for the enthalpy, and increased heat flux and friction quantity magnitudes. In certain conditions, the viscous dissipation of kinetic energy is also found to be a non-negligible portion of the total heat flux. To address these findings, a new characterization is introduced and derived systematically. Results demonstrate improved collapse in the log region slope to within ±18 % uncertainty. Following the procedure proposed by Patel et al. (2017), the full transformation is split into two component terms. The first term,θ * T , characterizes the log region slope and collapses well among all cases into a single, near-universal profile. The second term,θ * P , characterizes the non-universal log region shift parameter that deviates significantly for different operating conditions. After a predictive framework is developed that models the log region shift parameter using the results ofθ * P , the near-wall mean temperature profile in the logarithmic region can be described fully as a modelled analytical profile. This theoretical framework, by capturing the impact of select physical 934 A45-29 processes on characterizing the turbulence dynamics, thus guides the development of more accurate turbulence models for wall-bounded transcritical flows. dissipation rate¯ = τ ij ∂u i /∂x j /ρ, and η T = η u / Pr for the thermal Kolmogorov length scale. Across the majority of the channel, x ≈ (2.5-7.0)η u ≈ (3.5-11)η T and y ≈ (0.4-3.2)η u ≈ (0.5-8.3)η T -these resolutions are similar to those used in previous studies (Patel et al. 2017;Wan et al. 2020). The larger outlying values of x 12.7η 37.3η T are localized to near the hot wall, where η u and η T are not the relevant scales to characterize the flow. Even so, everywhere in the channel domain, the grid spacing in the current study is in line with previous estimations that state that the majority of the turbulent dissipation occurs at scales larger than 15η u in channel flows .
To ensure that the relevant small scales are sufficiently resolved, figure 23 shows the streamwise momentum energy spectrum, and figure 24 shows the enthalpy energy spectrum; each spectrum has been premultiplied by the wavenumber to show the resolution of the near-wall inner peak. Plots for case TR3 are shown to assess the grid resolution for the largest values of heat transfer simulated across all cases. Plots for case TR1.3 are shown to assess the grid resolution when the Widom line occurs very close to the wall, resulting in steep gradients in thermodynamic properties in the near-wall regime. Contour plots of the other four cases (TR1.9, TR1.4, TR1.25 and TR1) have been omitted for brevity, as these cases are less stringent versions of cases TR3 and TR1.3 from a grid resolution perspective. For both momentum and enthalpy, at least three decades of premultiplied energy are resolved. As another important observation attesting to the well-resolvedness 934 A45-30 of the simulations, for all cases the pile-up of under-resolved dissipative energy at smaller wavelengths (where under-resolution of dissipative scales is expected to occur, if present) appears minimal for all ζ * distances. Spanwise spectra yield similar observations and thus have been omitted. From these findings, we conclude that the computational grid used for the current simulations satisfies the requirements to be of DNS resolution.
Appendix B. Additional momentum and thermal field mean profiles:ū,T ,ρ and Pr * We here present mean profiles for temperature (figure 25a), density (figure 25b), and streamwise velocity (figure 25c). Results are computed using Reynolds averaging, as Reynolds-and Favre-averaged mean profiles for temperature and velocity are nearly identical for all cases (not shown). Increasing asymmetry between the cold and hot walls is observed with increasing temperature ratio.
To provide a quantitative evaluation of the relative importance of energy transport and storage mechanisms, the local Prandtl number Pr * is plotted for each case in figure 26. The behaviour of Pr * among different cases is quite distinct, especially in the region near the hot wall. Because variations inμ/λ are negligible compared to variations in c p , Pr * captures closely the mean behaviour of c p . As a consequence, the locations of the peaks in Pr * essentially coincide with the locations of the transition points of the Widom line; for all cases whose temperature ranges straddle T WL , the location of the transition point of the Widom line lies in the upper half of the channel (y/L y ≥ 0). As noted by Toki et al. (2020), spatial property variations in transcritical fluids are associated with volume 934 A45-31  Figure 23. Contour plots of pre-multiplied energy spectra for streamwise momentum fluctuations, k x Φ ρuu , where k x is the streamwise wavenumber, so that its inverse 1/k x is the streamwise wavelength. The energy spectra values have been normalized using characteristic density scaleρ and characteristic velocity scale u * τ = √ τ w /ρ. The wavelength on the x-axis is normalized with friction length μ w / √ τ w ρ w . The y-axis gives the semi-local wall-normal distance ζ * . Panel (a) shows contours for case TR3, and (b) shows contours for TR1.3. For each case, (i) displays contour plots near the cold wall, and (ii) displays contour plots near the hot wall.  Figure 24. Contour plots of pre-multiplied energy spectra for enthalpy fluctuations, k x Φ hh , where k x is the streamwise wavenumber, so that its inverse 1/k x is the streamwise wavelength. The energy spectra values have been normalized using a characteristic specific energy scale u * 2 τ = τ w /ρ. The wavelength on the x-axis is normalized with friction length μ w / √ τ w ρ w . The y-axis gives the semi-local wall-normal distance ζ * . Panel expansions that lead to differences in turbulent structures and transport mechanisms, when compared with the analogous features of constant-property turbulence. These differences result in observable discrepancies in velocity and temperature profiles. The quantity Pr * incorporates the effects of property variations on the temperature behaviour, shown through mathematical manipulation of the mean energy equation in § 3.4. As discussed in § 1, the behaviour of the mean temperature profile in the viscous sublayer is well-characterized as a linear function of wall-normal distance, with slope equal to Pr * . An observation of (3.21) and figure 20 shows that within the viscous sublayer, where we can approximateθ * ≈θ + and (1/c p )(dh/dθ) ≈ 1, the transformation indeed collapses to the relationθ + =Pr * ζ * . Thermal field in transcritical turbulent channel flows We zoom in on this part of the near-wall region in figure 27. For most of the ζ * < 5 region, good collapse with (B1) is observed, in agreement with previous studies in both variable-property and constant-property turbulent flows.

A45-33
Downloaded from https://www.cambridge.org/core.   Figure 28. Profiles of (a)θ * T and (b) the full transformationθ * , along with comparison with reference data by Kim & Moin (1989) (shown with dotted lines). For our current calculation, dashed lines indicate cold wall profiles, and solid lines indicate hot wall profiles. Note that because Kim & Moin (1989) used an internal heating source to generate the passive scalar, the numerator in the integrands of bothθ * T andθ * were multiplied by a factor 1 − y/L y to account properly for the heat flux profile. The relatively low friction Reynolds number Re τ = 180 for the reference data limits the ζ * extent of the profiles.