M Gelain

O Gicquel

A Couilleaux

R Vicquelin
email: ronan.vicquelin@centralesupelec.com

Physical study of the non-equilibrium development of a turbulent thermal boundary layer

The direct numerical simulation of a non-equilibrium turbulent heat transfer case is performed in a channel flow, where non-equilibrium is induced by a step change in surface temperature. The domain is thus made of two parts in the streamwise direction. Upstream, the flow is turbulent, homogeneous in temperature, and the channel walls are adiabatic. The inflow conditions are extracted from a recycling plane located further downstream, so that a fully developed turbulent adiabatic flow reaches the second part. In the domain located downstream, isothermal boundary conditions are prescribed at the walls. The boundary layer, initially at equilibrium, is perturbed by the abrupt change of boundary conditions, and a non-equilibrium transient phase is observed until, further downstream, the flow reaches a new equilibrium state, presenting a fully developed thermal boundary layer. The work aims at identifying the non-equilibrium effects that are expected to be encountered in comparable flows, while providing the means to understand them. In particular, the study allows for the identification of an inner region of the developing boundary layer where several quantities are at equilibrium. Other quantities, instead, exhibit a behaviour of their own, especially in proximity to the leading edge. The analysis is supported by mean and root-mean-square profiles of temperature and velocity, as well as by budgets of first-and second-order moment balance equations for the enthalpy and momentum turbulent fields.

Introduction

Turbulent heat transfer is encountered in numerous industrial applications, and in most cases the turbulent heat exchange takes place in non-equilibrium flows. It is the case, for instance, of the surface air-oil heat exchangers implemented in modern aircraft engines by-pass ducts, where an essentially temperature-homogeneous flow makes contact with a heated surface, and is therefore abruptly subjected to a temperature gradient.

It is the interest in the behaviour of the atmospheric boundary layer which drew first attention to turbulent heat transfer in non-equilibrium flows, both theoretically and experimentally (see [START_REF] Antonia | Response of a turbulent boundary layer to a step change in surface heat flux[END_REF]. The case where non-equilibrium is induced by a step change in surface temperature is the simplest example of non-equilibrium turbulent heat transfer, and numerous experimental works can be found in the literature. [START_REF] Johnson | Velocity, temperature and heat transfer measurements in a turbulent boundary layer downstream of a stepwise discontinuity in wall temperature[END_REF] studied the development of a thermal boundary layer on a smooth flat plate, and presented mean temperature and velocity profiles while focusing on fluctuating profiles in a later work [START_REF] Johnson | Velocity and temperature fluctuation measurements in a turbulent boundary layer downstream of a stepwise discontinuity in wall temperature[END_REF]. [START_REF] Blom | Experimental determination of the turbulent prandtl number in a developing temperature boundary layer. 4th int[END_REF] carried out a similar study, comparing mean temperature profiles to the theoretical predictions of [START_REF] Spalding | Heat transfer to a turbulent stream from a surface with a step-wise discontinuity in wall temperature[END_REF], and presenting the evolution of the turbulent Prandtl number. Similar experimental works can be found in [START_REF] Fulachier | Contribution à l'étude des analogies des champs dynamique et thermique dans une couche limite turbulente: effet de l'aspiration[END_REF], [START_REF] Hoffmann | The development of turbulent thermal layers on flat plates[END_REF], [START_REF] Ng | The turbulent boundary layer over a flat plate with strong stepwise heating[END_REF], [START_REF] Taylor | Heat transfer measurements in incompressible turbulent flat plate boundary layers with step wall temperature boundary conditions[END_REF] and, more recently, [START_REF] Biles | The design and validation of a thermal boundary layer wall plate[END_REF], while [START_REF] Antonia | Response of a turbulent boundary layer to a step change in surface heat flux[END_REF] and [START_REF] Teitel | A step change in wall heat flux in a turbulent channel flow[END_REF] studied the case of a step change in wall heat flux in a flat plate and turbulent channel flow configuration, respectively.

In spite of its Reynolds number limitation, direct numerical simulation (DNS) is undoubtedly the most accurate tool available for the investigation of turbulent heat transfer, since all the turbulent scales are solved. The first DNS addressing turbulent heat transfer is that of [START_REF] Kim | Transport of passive scalars in a turbulent channel flow[END_REF], where the transport of three passive scalars at different molecular Prandtl numbers is considered in a fully developed channel flow at a friction Reynolds number of Re τ = 180. Several works followed, aiming at understanding the influence on the flow statistics of the Reynolds number and molecular Prandtl number (e.g. [START_REF] Papavassiliou | Transport of a passive scalar in a turbulent channel flow[END_REF][START_REF] Kawamura | Dns of turbulent heat transfer in channel flow with respect to reynolds and prandtl number effects[END_REF][START_REF] Abe | Direct numerical simulation of a fully developed turbulent channel flow with respect to the reynolds number dependence[END_REF], as well as of the isothermal or iso-flux boundary conditions (e.g. [START_REF] Kasagi | Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy[END_REF]Kawamura et al. 2000). Other authors focused on the coupling between turbulence and temperature gradient, an analysis which is possible only if temperature is handled as an active scalar. Most of these DNS were performed in supersonic channel flows (see [START_REF] Coleman | A numerical study of turbulent supersonic isothermal-wall channel flow[END_REF][START_REF] Huang | Compressible turbulent channel flows: Dns results and modelling[END_REF][START_REF] Morinishi | Direct numerical simulation of compressible turbulent channel flow between adiabatic and isothermal walls[END_REF][START_REF] Tamano | Effect of different thermal wall boundary conditions on compressible turbulent channel flow at m= 1.5[END_REF] while only a few focused on low speed flows with high temperature gradients (see [START_REF] Nicoud | Numerical study of a channel flow with variable properties[END_REF][START_REF] Toutant | Turbulence statistics in a fully developed channel flow submitted to a high temperature gradient[END_REF].

All the numerical works mentioned so far represent cases of equilibrium flows, and are undoubtedly the reference for understanding turbulent heat transfer in all its aspects (flow statistics, turbulent Prandtl number, wall scaling and much more). Nevertheless, questions arise concerning the validity of these findings in non-equilibrium configurations. Once more, DNS can be a powerful tool for the analysis of such flows, yet the literature is extremely less abundant on the matter. [START_REF] Seki | Dns of turbulent heat transfer in a channel flow with streamwisely varying thermal boundary condition[END_REF] performed the DNS of a fully developed channel flow with temperature as a passive scalar, where the wall temperature is constant everywhere but in a small fraction of the bottom wall, where it evolves along the streamwise direction attaining a peak. The step change of wall temperature significantly perturbs the mean and fluctuating temperature, as well as turbulent Prandtl number. [START_REF] Hattori | Direct numerical simulation of stable and unstable turbulent thermal boundary layers[END_REF], [START_REF] Hattori | Dns study of effects of suddenly-vanishing wall heating in turbulent boundary layer[END_REF] and [START_REF] Hattori | Dns, les and rans of turbulent heat transfer in boundary layer with suddenly changing wall thermal conditions[END_REF] performed several DNS of non-equilibrium thermal boundary layers in the flat plate configuration. [START_REF] Hattori | Direct numerical simulation of stable and unstable turbulent thermal boundary layers[END_REF] analysed the effects of buoyancy on mean and fluctuating properties, as well as on the turbulent budgets for a turbulent boundary layer at equilibrium facing a step change in wall temperature. In [START_REF] Hattori | Dns study of effects of suddenly-vanishing wall heating in turbulent boundary layer[END_REF], non-equilibrium is induced through the sudden vanishing of wall-heat flux and, in one of the two cases presented, with the addition of a forward facing step. Similar analyses are carried out in [START_REF] Hattori | Dns, les and rans of turbulent heat transfer in boundary layer with suddenly changing wall thermal conditions[END_REF], yet temperature, in this case, is not handled as a passive scalar. Large-eddy simulation (LES) is certainly another valid tool for studying non-equilibrium turbulent heat transfer yet, despite the large turbulent scales being solved, a certain degree of modelling is needed for the smaller scales. [START_REF] Sanchez | Large eddy simulations of thermal boundary layer spatial development in a turbulent channel flow[END_REF] and [START_REF] Bellec | Large eddy simulations of thermal boundary layer developments in a turbulent channel flow under asymmetrical heating[END_REF] performed the LES of a temperature-homogeneous fully developed channel flow making contact with anisothermal walls for Re τ = 180 and Re τ = 395, respectively. They showed the evolution of the mean and fluctuating temperature profiles, as well as the impact on the mean and fluctuating velocity components.

The objective of this paper is to deepen the knowledge of turbulent non-equilibrium heat transfer, which, as the literature review shows, has been investigated in a limited number of studies, especially in the case of flows with temperature-dependent properties.

The intention is not only to describe and characterise a non-equilibrium flow, but also its gradual evolution towards a new equilibrium state.

To achieve the scope of the study, we perform the direct numerical simulation of a channel flow, where a turbulent, fully developed, temperature-homogeneous flow at Re τ = 395 makes contact with an isothermal wall. The step change in surface temperature leads to the development of a thermal layer, and, since temperature is an active scalar, the velocity boundary layer is also modified. The length of the isothermal wall allows the boundary layer to reach a new equilibrium state, characterised by a fully developed thermal boundary layer. The simulation is performed under one single flow regime, in terms of Reynolds number, Mach number and heat flux parameter. Nevertheless, the case is representative of a larger class of non-equilibrium flows where the thermal boundary layer responses to a sudden perturbation. Therefore, beyond the case-specific conditions, the study primarily aims at identifying the physical phenomena that are expected to be encountered in comparable non-equilibrium flows, while providing the theoretical means to understand them. In particular, we will see that the sudden perturbation is associated with (i) a local effect near the hot-wall leading-edge specific to the change in wall temperature, and (ii) a more general subsequent longitudinal and transverse development of the boundary layer towards a new equilibrium.

The work is organised as follows. In section 2, the problem is described in detail; the governing equations are presented as well as details about the study configuration are given. In section 3, the numerical set-up and the proposed configuration are validated with respect to reference results found in the literature for equilibrium flows. In section 4, the results are presented. In §4.1, the behaviour of the flow is analysed in what appears to be the most perturbed region, i.e., the leading edge of the isothermal wall. In §4.2, an equilibrium region which characterises the development of the thermal boundary layer is identified and quantitatively defined. In §4.3, the findings of §4.2 are contrasted to the evolution of several turbulent quantities. This is followed by the conclusions, in §5.

Problem description

In this section, details about governing equations and numerical scheme ( §2.1) as well as geometry, mesh and boundary conditions ( §2.2) are given.

Governing equations and numerical scheme

The full compressible Navier-Stokes equations are solved in the fluid without any turbulence model. The set of equations is the following:

∂ρ ∂t + ∂ρu i ∂x i = 0 , (2.1
)

∂ρu i ∂t + ∂ρu i u j ∂x j = - ∂p ∂x i + ∂τ ij ∂x j , (2.2
)

∂ρh ∂t + ∂ρu j h ∂x j = Dp Dt - ∂q cd j ∂x j + τ ij ∂u i ∂x j + S ener , (2.3)
where ρ, u i , p, h are respectively the mass density, velocity components, static pressure and enthalpy per mass unit of the fluid; the fluid is considered to be an ideal gas with temperature-tabulated thermodynamic coefficients and mass-specific gas constant r = 288.18 J/(kg.K) (assuming a mixture of 79% N 2 and 21% O 2 ) and state equation p = ρrT , where T is the static temperature; gravity is neglected and not included in momentum 2); the fluid is Newtonian and its dynamic viscosity µ is computed through the Sutherland law and is, therefore, temperature-dependent; the conductive heat flux q cd i follows the classic Fourier's law and the fluid thermal conductivity λ is computed via the Prandtl number Pr = 0.71; finally, S ener is a source term added to the energy equation ( 2.3), the role of which will be clarified in §2.2.

The set of equations is solved by the parallel code AVBP (Schonfeld & Rudgyard 1999;Moureau et al. 2005) using a time-explicit finite-element two-step Taylor-Galerkin scheme (Colin & Rudgyard 2000) which provides third-order accuracy in space and time.

Geometry, mesh and boundary conditions

The geometry of the configuration studied is shown in figure 1. The computational domain is a channel flow of size 22πδ × 2δ × πδ (where δ = 0.002 m) which is made of two parts in the streamwise direction.

The upstream part has size 4πδ × 2δ × πδ, the upper and lower walls (with respect to the Y direction) are adiabatic, and a no-slip boundary condition is prescribed, while periodic boundary conditions are applied in the spanwise direction Z. Differently from what happens in bi-periodic channel flows, in this case the flow is driven by a streamwise pressure gradient compensating the head losses generated along the whole domain shown in figure 1. In the Z direction the mesh is uniform; in the wall-normal direction the mesh size, expressed in wall units, varies from ∆Y + = 0.75 at the upper and lower walls, to ∆Y + = 7 at the centre of the channel; along the streamwise direction, the mesh size decreases with a constant ratio 1.05 from ∆X + = 12 at the inlet to ∆X + = 1 at the interface with the downstream part of the channel flow (see figure 1). The role of this part of the domain is to generate a temperature-homogeneous boundary layer at equilibrium, with the regime conditions specified in table 1. In order to do so at a moderate cost, a recycling strategy is used: the three velocity components and the temperature imposed at the domain inlet via the Navier-Stokes Characteristic Boundary Condition (NSCBC) formalism [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]Moureau et al. 2005) are extracted from a recycling plane situated downstream at a distance of 2πδ. The energy source term S ener added to (2.3), is zero in this portion of the channel flow.

The downstream part of the domain has size [18πδ, 2δ, πδ], the upper and lower walls are isothermal, both at temperature T w = 400 K, and a no-slip boundary condition is prescribed, while, analogously to the upstream part, the domain is periodic in the spanwise direction. At the outlet, pressure is imposed via the NSCBC formalism. The mesh has the same point distribution of the upstream sub-domain along the Y and Z directions, while the mesh size increases with a constant ratio 1.05 from ∆X + = 1 at the

Size nX × nY × nz ∆X + ∆Y + ∆Z + Reτ M T b (K) [4πδ, 2δ, πδ] 420 × 179 × 200 [1.0 -12.0] [0.75 -7.0]
6.2 395 0.16 304.5

Table 1: Size, mesh, spatial resolution and regime conditions of upstream (adiabatic) sub-domain: n i is the number of grid points and ∆X + i the spatial resolution in wallunits; Re τ is the friction Reynolds number, M the average Mach number and T b the bulk temperature.

Size

nX × nY × nz ∆X + ∆Y + ∆Z + [18πδ, 2δ, πδ] 1771 × 179 × 200 [1.0 -12.0] [0.75 -7.0] 6.2
Table 2: Size, mesh and spatial resolution of downstream (isothermal) sub-domain.

adiabatic-isothermal interface until ∆X + = 12 further downstream, where wall units are still referred to the adiabatic conditions of the upstream sub-domain, which are the most restrictive in terms of spatial resolution. This refinement along the streamwise direction serves to well capture the very first development region of the thermal boundary layer at the leading edge of the isothermal wall, where streamwise gradients can be important. The role of the downstream part is to allow the thermal boundary layer to be fully developed and to attain a new equilibrium state before the exit. The fully developed regime is characterized by a scaled mean temperature profile which is homogeneous in the streamwise direction. The role of the source term S ener added to the energy equation (2.3), is that of guaranteeing that a temperature gradient is preserved between the wall and the centre of the channel, allowing for a better appreciation of the fully developed regime. The mean temperature profile, therefore, results from the balance between the prescribed source term and the wall heat fluxes of the established flow q w eq : q w eq ≈ -S ener δ .

(2.4)

The source term S ener takes a constant value specifically designed to obtain at equilibrium a bulk temperature of T eq b = 304.5 K, i.e., the same of the upstream part of the domain (see table 1), which yields a temperature of T eq c ≈ 290 K at the centre of the channel. Thanks to (2.4), the value at equilibrium of the heat flux parameter B q and of the Nusselt number Nu can be deduced:

B eq q = - q w eq ρ w c pw u τ T w Nu eq = -2δ q w eq λ w (T w -T c ) , (2.5)
where u τ = τ w /ρ w , with τ w the wall shear stress, is the friction velocity, c p is the constant-pressure specific heat and the subscript (•) w denotes a fluid quantity taken at the wall. Table 2 summarises the size and the mesh resolution, while table 3 details the regime conditions of the downstream part of the domain.

Equilibrium states and validation

This study is characterised by two distinct equilibrium states, between which the flow evolves along the streamwise direction. The former is an adiabatic, fully developed Table 3: Regime conditions of downstream (isothermal) sub-domain. T w /T in is the temperature ratio between the wall and the inlet.

turbulent flow in the upstream part of the domain; the latter is a fully developed turbulent thermal boundary layer at the outlet. The objective of this section is twofold: first, to validate the numerical setup introduced in §2.1 with respect to reference results of the literature; second, to verify that the configuration described in §2.2 allows us to to meet the scope of the study, i.e., to observe the evolution of the boundary layer between the aforementioned equilibrium states.

In the following, (•) and (•) denote Reynolds and Favre averages, while (•) and (•) denote their respective fluctuating parts. Spatially, quantities are only averaged along the Z axis, which is the only homogeneous direction for this configuration. Given the streamwise evolution of the boundary layer, unless specified, scaling of quantities (dimensionless numbers included) has to be intended locally with respect to the X coordinate. The origin of the axes is placed at the interface between the adiabatic and isothermal walls, so that the adimensionalised X coordinate x/δ is defined in the range [-4π, 18π] between the inlet and the outlet. The subscripts (•) w and (•) c indicate that the given property is evaluated at the wall and at the centre of channel, respectively. The superscript (•) + denotes classic wall scaling, that is:

y + = ρ w yu τ µ w and u + i = u i u τ with u τ = τ w ρ w , (3.1)
T + = |T -T w | T τ with T τ = |q w | ρ w c pw u τ . (3.2)
The superscript (•) * , instead, denotes semi-local scaling (see [START_REF] Huang | Compressible turbulent channel flows: Dns results and modelling[END_REF][START_REF] Patel | Semilocal scaling and turbulence modulation in variable property turbulent channel flows[END_REF], which uses local values (with respect to the wall-normal direction Y ) of ρ, µ and c p , instead of their values at the wall to account for variable properties In the following sections, the behaviour of the flow is analysed at two streamwise locations: at x/δ ≈ 0 in §3.1, where an adiabatic, fully developed turbulent flow is expected; at x/δ ≈ 18π in §3.2, where the development of the thermal boundary layer should be concluded.

Upstream equilibrium state

Figure 2 shows the mean streamwise velocity and the three root-mean square (r.m.s.) velocity profiles at x/δ = -0.1π, i.e., slightly downstream of the leading edge of the isothermal wall. Results at x/δ ≈ 0, represented by the black solid lines, are compared with two reference data. The first reference data, represented by the gray lines, have been obtained by simulating a classic bi-periodic channel flow of size [2πδ, 2δ, πδ] with the regime conditions of table 1; the second reference data, represented by the symbols, are taken from Kawamura et al. (2000).

The collapsing between the small channel flow results and Kawamura et al. (2000), is required to validate the numerical setup of §2.1; as can be seen, good agreement is obtained for both the mean streamwise velocity and the r.m.s. velocities. -present results with 2πδ-long bi-periodic channel flow (gray) and at x/δ ≈ 0 (black); results from [START_REF] Kawamura | Dns of turbulent heat transfer in channel flow with respect to reynolds and prandtl number effects[END_REF]. Profiles of r.m.s. streamwise, wall-normal and spanwise velocity respectively (b): --, -• -• -and ----present results (colors as above); , and results from [START_REF] Kawamura | Dns of turbulent heat transfer in channel flow with respect to reynolds and prandtl number effects[END_REF].

The accord between the small channel flow and the present results at x/δ ≈ 0, instead, is required to validate the study configuration described in §2.2; excellent agreement is obtained on the mean streamwise velocity, which shows that the flow reaching the downstream part of the domain is a fully developed turbulent temperature-homogeneous boundary layer as desired; on the other hand, even though good agreement is obtained for the wall-normal and spanwise r.m.s. velocities, a disparity of around 10% is observed for the streamwise component in the log-layer; a possible explanation for this behaviour is given in the following section and, more in detail, in Supplementary Materials A.

Downstream equilibrium state

At x/δ ≈ 18π, the flow is expected to have evolved towards a new equilibrium state, characterised by a fully developed thermal boundary layer.

Figure 3 shows the mean streamwise and r.m.s. velocities at x/δ = 287/16π compared to both the results of [START_REF] Kawamura | Dns of turbulent heat transfer in channel flow with respect to reynolds and prandtl number effects[END_REF] and, similarly to what is done in §3.1, to the results obtained with a classic bi-periodic channel flow at the regime conditions of table 3. Note that, in this case, semi-local scaling is used. As can be seen, good agreement is obtained for all the profiles, and even the streamwise r.m.s. velocity is considerably closer to the reference results. The disparity illustrated in figure 2b, can thus be attributed to some local perturbation of the upstream flow, most probably induced by the recycling method described in §2.2, the effects of which disappear further downstream.

The profiles of mean and r.m.s. temperature are instead shown in figure 4. Good agreement is obtained among all the curves for the mean temperature, while a slight discrepancy is observed on the r.m.s. profile. Therefore, if on the one hand the average thermal boundary layer seems to be fully developed at the outlet, on the other hand a longer isothermal wall would have been necessary to observe the same level of convergence for the temperature fluctuations.

Both figures 3 & 4 also illustrate how, thanks to semi-local scaling, fluid property variations in the wall-normal direction are properly taken into account, with all velocity and temperature profiles collapsing to those obtained by [START_REF] Kawamura | Dns of turbulent heat transfer in channel flow with respect to reynolds and prandtl number effects[END_REF] with a passive scalar. Fluid mean property variations, at least at equilibrium, seem thus to --present results with 2πδ-long bi-periodic channel flow (gray) and at x/δ ≈ 18π; results from [START_REF] Kawamura | Dns of turbulent heat transfer in channel flow with respect to reynolds and prandtl number effects[END_REF].

be the only remarkable impact of the heated wall on the flow, as suggested by [START_REF] Morkovin | Effects of compressibility on turbulent flows[END_REF] and seen in several studies involving compressible flows (e.g. [START_REF] Huang | Van driest transformation and compressible wallbounded flows[END_REF][START_REF] Huang | Compressible turbulent channel flows: Dns results and modelling[END_REF][START_REF] Nicoud | Numerical study of a channel flow with variable properties[END_REF]). However, the effects of the temperature gradient on the flow can be much more important in other portions of the channel, namely in proximity of the leading edge of the isothermal wall. This can be seen in figure 5, which plots the evolution of the Nusselt number and of the heat flux parameter along the channel; both are very important near the leading edge, and the behaviour of the flow needs to be investigated in that region. The later convergence of Nu and B q towards the values of table 3, instead, provides further confirmation of the evolution towards equilibrium of the thermal boundary layer.

As a whole, thanks to the results of §3.1 and of this section, it is possible to conclude that the given configuration and the numerical setup allow us to meet the aim of the study and study the boundary layer development between two equilibrated conditions.

Results

This work aims at analysing the non-equilibrium evolution of an adiabatic, equilibrium boundary layer towards a new equilibrium state, which is characterised by a fully developed thermal boundary layer. These two equilibrium states have been presented and analysed in §3. The objective of this section, instead, is to investigate the most relevant phenomena which can be identified in the non-equilibrium phase of the flow. Although the DNS is performed at one single flow regime, certain physical events and flow characteristics are expected to be distinguishable in phenomenologically similar nonequilibrium flows with a streamwise development, regardless of the particular regime. We aim at identifying these phenomena and at providing the theoretical tools to investigate them.

In §4.1, the region in proximity of the leading edge is analysed, where due to the abrupt change of boundary conditions, non-equilibrium effects are the strongest. Particular attention is devoted to the velocity profiles and to the momentum balance, which, at x/δ ≈ 0, are significantly perturbed by the temperature gradient. Figure 6, for instance, shows some instantaneous contours of the wall shear stress near the leading edge. Observe how the contours are condensed around x/δ = 0 to the point that the interface becomes clearly visible, signalling a discontinuity.

In §4.2, the gradual development of the thermal boundary layer is investigated. Along the isothermal wall, the edge of the growing thermal layer (of which an instantaneous visualisation is given in figure 7) is at a non-equilibrium state. Yet, in a more inner region, a sub-layer can be identified, where several quantities are equilibrated. The notion of equilibrium and non-equilibrium layers are clarified in this section, and a definition based on the mean energy balance is proposed.

However, the analysis of section 4.2 is effectuated in a mean sense and, therefore, the findings concerning the equilibrium sub-layer need to be contrasted to the evolution of several turbulent quantities. This is done in §4.3, where the behaviour of the fluctuating temperature (of which an instantaneous field is provided in figure 8), turbulent heat flux and turbulent Prandtl number is investigated.

In the following, the same notation of section 3 is used.

Analysis of the leading edge

The abrupt change of boundary conditions taking place at x/δ = 0 generates enduring

• • • • • • equilibrium profile of §3.1; -• -• -equilibrium profile of §3.2.
non-equilibrium effects which, as will be seen more in detail in §4.2, do not vanish completely even at the distance of x/δ = 18π. However, there is a portion of the channel where these effects are so strong that the flow has a behaviour of its own that needs to be carefully investigated. This region is the leading edge, where the temperature gradient is the most elevated, as a consequence of the sudden change of thermal wall conditions (see the evolution of Nu and B q in figure 5). Since the temperature is treated as an active scalar in our simulation, an impact on the velocity fields and momentum balance is expected. On the other hand, in this region the thermal boundary layer is too thin for any appreciable phenomenon to be observed on the temperature.

Figure 9 shows the mean streamwise velocity for x/δ ranging between [-0.06, 0.18] compared to the canonical profiles seen in §3. Even if semi-local scaling is used, none of the profiles agrees with the equilibrium ones, showing that, at these close distances from the leading edge, the perturbation cannot be ascribed to the only variation of mean fluid properties. This is particularly true for the profile at x/δ = 0.015, for which the velocity boundary layer seems to be at equilibrium only for y * 2, i.e., in the viscous sub-layer. The discrepancy of the profile at x/δ = -0.06 with respect to the equilibrium profile of §3.1, instead, shows how the perturbation slightly propagates upstream. A similar impact on the streamwise velocity has been observed by [START_REF] Sanchez | Large eddy simulations of thermal boundary layer spatial development in a turbulent channel flow[END_REF] with a higher B q , arguing that the perturbation should be ascribed to the increase of the wall-normal velocity induced by the temperature gradient. It is certainly true that the wall-normal velocity considerably increases near the leading edge, as it can be seen in figure 10. Yet, the wall-normal velocity is also important at x/δ = -0.06, and, above all, always greater than zero. If the wall-normal velocity were the only responsible for the deviation of the streamwise velocity, the impact would be expected to be similar to the one observed downstream of the leading edge, where the normal velocity is also always positive and of the same order of magnitude. Instead, as can be seen in figure 9, it is not the case, since an upper and lower shift is respectively observed at x/δ = -0.06 and, for example, x/δ = 0.18.

Hence, we propose a different interpretation, i.e., that the destabilisation of the boundary layer is due to to the abrupt variation of the wall shear stress. Figure 11a shows the evolution of the skin friction coefficient

C f = τ w / 1 2 ρ b u 2 b
, where ρ b and u b are respectively the bulk density and the bulk velocity. At x/δ = 0 + , the skin friction coefficient is at its maximum, which can be explained by the abrupt increase of the fluid dynamic viscosity at the leading edge of the isothermal wall. Then, the following relaxation and adaptation of the velocity gradient at the wall to the increased viscosity, leads to a quickly decreasing C f until it attains a plateau. At x/δ = 0 -, instead, the skin friction coefficient is at its minimum. This is due to the fact that the flow, being subsonic, perceives the upcoming isothermal wall, and adapts by decreasing the gradient of streamwise velocity at the wall, while the dynamic viscosity remains constant since the wall is adiabatic. This perturbation, limited to a very small portion of the channel flow between x/δ ∈ [-1/2π, 1/2π], results in a local perturbation of the wall streamwise pressure gradient, as it can be seen in figure 11b. While ∂p ∂x | w / τw δ ≈ -1 far from the leading edge, a positive (and thus adverse) and strongly negative (and thus favourable) pressure gradient is observed at x/δ = 0 -and x/δ = 0 + , respectively.

In order to take the effect of the streamwise pressure gradient into account, we introduce, following [START_REF] Simpson | A model for the backflow mean velocity profile[END_REF], the velocity scale u p defined as:

u p = µ ρ 2 ∂p ∂x w 1 3 (4.1)
u τ p = u 2 τ + u 2 p . (4.2)
We propose here a slight modification of this definition, which distinguishes positive and negative pressure gradients:

u τ p = u 2 τ + sign ∂p ∂x w u 2 p . (4.3)
In this manner, the decrease of u τ in case of adverse pressure gradient, and its increase in case of favourable pressure gradient, are compensated by u p , correcting the upward and downward deviation of the streamwise velocity at x/δ = 0 -and x/δ = 0 + , respectively. We can thus introduce the following wall scaling:

y * τ p = yu τ p ν u * τ p = u u τ p , (4.4)
where, analogously to semi-local scaling, fluid properties depend on the wall-normal direction in order to take the effect of the temperature gradient into account.

Figure 12a shows the same velocity profiles of figure 9 in the newly introduced wall scaling, while figure 12b illustrates the r.m.s. velocity profiles. Very good agreement is obtained between the different profiles in every region of the boundary layer, with the only exception of the r.m.s. streamwise velocity, for which the upward shift in the log layer has in any case to be ascribed to what is explained in §3.2.

The fact that the velocity profiles are corrected by adopting a different wall-scaling, suggests that, in our case, the perturbation is confined in the near-wall region and does not propagate towards the core of the flow. In the most general case, however, the leading edge effects are expected to become more important as B q increases, with, in particular, larger wall-normal velocities and stronger variations of the pressure gradient. In these conditions, a modified wall-scaling might not suffice to correct these effects.

In order to corroborate our interpretation of the leading edge perturbations, the where on the right-hand side of the equation there are two mean convective terms (I x and I y ), the pressure gradient (II), two viscous terms (III x and III y ) and two terms associated with turbulent transport (IV x and IV y ). Figure 13 shows the different terms non-dimensionalised with respect to the average local wall shear stress τ w as a function of y/δ for two values of x/δ, slightly upstream and downstream of x/δ = 0. Figure 13a is relative to x/δ = -0.06. At this close distance from the leading edge, the terms of (4.5) are strongly perturbed. The pressure gradient flux does not follow a straight line, showing that ∂p/∂x is not uniform along the wall-normal direction and that the pressure field is bidimensional; it is also evident how the streamwise pressure gradient is positive (and thus adverse) until y/δ ≈ 0.05. Both mean convective terms are important; note the positive slope of the streamwise convective term near the wall, indicating that ∂ u/∂x is negative and resulting in the aforementioned decrease of τ w . The remaining terms seem not to be significantly modified, indicating that the perturbation of the pressure gradient is fully compensated by the two convective terms. This is in agreement with the hypothesis that the non-equilibrium contributions to the momentum balance tend to self-compensate (see [START_REF] Larsson | Large eddy simulation with modeled wall-stress: recent progress and future directions[END_REF].

Figure 13b is relative to x/δ = 0.18. The pressure gradient flux still does not follow a straight line and the pressure field is thus two-dimensional; the positive slope of the pressure gradient term indicates that ∂p/∂x is now negative for every y/δ. Both mean convective terms are still important; near the wall, like at x/δ = -0.06, ∂ u/∂x < 0, indicating that τ w is decreasing. The perturbation of the momentum fluxes becomes negligible at x/δ ≈ π/2, and completely disappears by x/δ ≈ 3π, where the canonic momentum fluxes are finally retrieved, and the effects of the leading edge can be considered fully dissipated (not shown here).

Note that the results described in this section can only be observed if the temperature is handled as an active scalar, i.e., if the momentum and energy equations are coupled.

Identification of an equilibrium sub-layer

Observe the behaviour of the instantaneous temperature iso-surface illustrated in figure 7. Near the leading edge of the isothermal wall, the thermal layer is extremely thin, and the heated pockets are entirely located in the purely shear region of the flow. As the distance increases, these pockets thicken, breach through the viscous region and gradually become turbulent, until a fully developed state is attained at the outlet. The goal of this section is to describe the development of the thermal boundary layer from the leading edge to the outlet, and to identify, thanks to the mean energy balance, the different non-equilibrium regions of the flow.

The progressive development of the thermal boundary layer is shown in figure 14, where several mean temperature profiles are plotted at different x/δ, and compared to the equilibrium profile of section 3.2. For every x/δ, three distinct parts of the developing thermal boundary layer can be clearly identified. The first is the near wall region where, for y * ∈ [0, h * eq ], where h * eq depends on and increases with x/δ, the thermal boundary layer is fully developed and the mean temperature profile agrees with the one at equilibrium; note that this equilibrated behaviour is visible even at x/δ = 0.015, where the leading edge effects described in §4.1 are very much present. The second, for y * ∈ [h * n-eq , δ * ], where h * n-eq also increases with x/δ, is a region that is still not affected by the isothermal wall and the temperature profile is flat. The third, for y * ∈ [h * eq , h * n-eq ], is instead the actual non equilibrium developing portion of the thermal boundary layer, where the mean temperature profile is neither flat nor agreeing with the equilibrium one. As x/δ increases, both h * eq and h * n-eq tend towards the mid-height of the channel δ * , and the non equilibrium region disappears.

The existence of an equilibrium region in the development of the thermal boundary layer has been shown experimentally, for example by [START_REF] Blom | Experimental determination of the turbulent prandtl number in a developing temperature boundary layer. 4th int[END_REF] (who called this zone adapted region) or [START_REF] Teitel | A step change in wall heat flux in a turbulent channel flow[END_REF]. The objective is to analyse the evolution of the equilibrium layer (delimited by h * eq ) and of the non-equilibrium region (delimited by h * n-eq ), yet it is first necessary to quantitatively define the notions of equilibrium and non-equilibrium.

We propose to do so through the analysis of the different terms of the energy balance. The integration along the wall-normal direction of the local averaged energy equation gives a local heat budget for every crosswise section: where on the right-hand side of the equation there are two mean convective terms (I x and I y ), two mean conductive terms (II x and II y ), two terms associated with turbulent heat transport (III x and III y ), one flux combining compressibility effects (IV ), and the source term contribution (V ). The evolution of these terms, normalised by q w , is shown for several x/δ in figure 15. Figures 15a & 15b are relative to x/δ = 0.73 and x/δ = 2π, respectively. At these short distances from the leading edge, many contributions which do not appear in equilibrium channel flows, are important. It is the case, for example, of the wall-normal convective flux since, as illustrated in figure 10, the mean vertical velocity is considerably higher compared to equilibrium. It is also the case of the streamwise turbulent flux, indicating that the correlation u h is evolving and increasing (in module). Yet, the preponderant contribution is that of the streamwise convective flux. Its positive slope near the wall indicates that the thermal boundary layer is developing and, as can be seen, the wallnormal conductive and turbulent fluxes are active. Farther away from the wall, its slope is negative and linear; since all the other contributions are constant, it is evident that the only active phenomenon is the uniform cooling caused by the source term. It is thus a portion of the boundary layer which still has not perceived the presence of the heated wall.

q w (x) =
Further away from the leading edge, the canonic heat fluxes are gradually retrieved. At x/δ = 5.7π (see figure 15c), ρ v h > 0 everywhere but at y/δ = 1, suggesting that the impact of the heated wall has by now reached the centre of the channel. At x/δ = 16.5π (see figure 15d), the thermal boundary layer seems to be at equilibrium. However, notice how at y/δ ≈ 1 the streamwise convective flux is negative and decreasing, indicating that the temperature at the centre of the channel is still cooling under the action of the source term. This signals that even at the outlet the thermal boundary layer is not perfectly developed, given the limited extent of the simulated channel.

The evolution of the different contributions to the energy balance shows that the fluxes can be divided into two categories. The former is that of the equilibrium terms, which tend to become more and more important with x/δ until stabilisation. The contributions that belong to this category are the wall-normal conductive flux and the wall-normal turbulent flux, which, summed, form the mean equilibrium flux q eq . The latter is that of the nonequilibrium terms which, even if potentially preponderant near the leading edge, tend to become negligible and disappear with x/δ. The contributions belonging to this category are the streamwise and wall-normal convective terms as well as the streamwise turbulent flux, which form the mean non-equilibrium flux q n-eq . The remaining contributions, at least in the present case, are always negligible. Now, defining the average total heat flux q tot as:

q tot = q w (x) + yS ener , (4.7)
(4.6) can be rearranged in the following form:

q tot = q eq + q n-eq , (4.8)

where the compressibility effects and the streamwise conductive flux have been neglected. For every x/δ and y/δ, it is thus possible to define the ratio:

R n-eq tot = |q n-eq | |q eq | + |q n-eq | , (4.9)

quantifying the importance of the non-equilibrium terms in the energy balance.

Figure 16a shows the evolution of R n-eq tot along the wall-normal direction for different x/δ. The ratio R n-eq tot allows to define a quantitative criterion for distinguishing the three aforementioned regions of the developing thermal boundary layer. Indeed, as it can Figure 16: Evolution of R n-eq tot with respect to y/δ for x/δ = 0.18, 1.5π, 2.5π and 16.5π (from light gray to black) (a). Evolution of h * eq , h * n-eq and δ * 99% along the channel flow:

• • • • • • h * eq ; -• -• -h * n-eq ; --δ * 99% (b).
be seen in figure 16a, for every x/δ there is a near-wall portion of the boundary layer where the equilibrium fluxes are predominant with respect to the non-equilibrium ones; consequently R n-eq tot ≈ 0 and the energy equation can be simply approximated as

q w + yS ener ≈ q cd y + ρ v h .
The equilibrium region of the boundary layer, delimited by h * eq (x), can thus be defined as the region laying beneath an isoline of R n-eq tot in the x-y plane, for example R n-eq tot = 0.1. For every x/δ, there is also a portion of the boundary layer adjacent to the centre of the channel, where the non-equilibrium terms are predominant; consequently R n-eq tot ≈ 1, the only active physical phenomenon is the cooling caused by the source term and the mean temperature profile is flat along the wall-normal direction. This region of the boundary layer, delimited by h * n-eq (x), can thus be defined as the region laying above an isoline of R n-eq tot , for example R n-eq tot = 0.9. Finally, for every x/δ and for y * ∈ [h * eq , h * n-eq ], the equilibrium and non-equilibrium fluxes are of the same order magnitude, and the thermal boundary layer is in development.

Figure 16b shows the evolution of h * eq (defined as the isoline R n-eq eq = 0.1) and h * n-eq (defined as the isoline R n-eq eq = 0.9) compared to the wall-scaled thermal boundary layer thickness δ * 99% := h * : T w -T (h * ) / T w -T c (x) = 0.99 . The iso-line h * n-eq is defined by a very loose criterion on R n-eq tot . Therefore, the advance of h * n-eq towards the centre of the channel can benefit from the slightest expansion of the wall-normal heat flux, which takes place in regions of the flow where heat transport is efficiently carried out by turbulence. Consequently, the region beyond h * n-eq quickly decays and disappears before x/δ = 4π. Downstream of this crosswise section, the effect of the heated wall has reached the centre of the channel and, as shown in figure 14, the mean temperature profile is nowhere flat any more.

The iso-line h * eq , instead, is defined by a very stringent criterion on the development of the equilibrium terms. Therefore, its growth is initially driven by the slow development of the wall-normal conductive flux, which takes place in regions of the flow where heat transport is essentially carried out by diffusion. As h * eq grows and reaches the most turbulent parts of the boundary layer, its expansion accelerates, and, as can be qualitatively seen in figure 16b, at x/δ ≈ 10π its growth rate is similar to that of h * n-eq at x/δ ≈ 3π, where h * eq ≈ h * n-eq ≈ 300. However, h * eq does not reach the centre of the channel, and attains a peak at y * ≈ 350 before slowly decreasing. The criterion R n-eq tot = 0.1 is thus so stringent that the slight cooling of the source term at the centre of the channel (recall figure 15d) becomes appreciable at a certain x/δ. With a higher criterion for h * eq , the same behaviour of h * n-eq would be observed. However, note that h * eq > δ * 99% , indicating that the equilibrium sub-layer attains in any case an essentially temperature-homogeneous region.

The evolution of h * eq and h * n-eq , is understandably case-dependent. However, note that their definition, based on R n-eq tot , is general and can be used in all cases. By computing or estimating the size of h * eq and comparing it with respect to a characteristic length of a given flow configuration, one can thus quantify the importance of the non-equilibrium effects in the thermal boundary layer.

Finally, figure 17a the evolution of R n-eq tot as a function of x/δ for several wall distances. The intersection of every iso-y * with the isoline R n-eq tot = 0.1 can be seen as the non-equilibrium distance associated to that wall distance, i.e., the distance from the leading edge which is necessary for the thermal boundary layer to be at equilibrium between the wall and the specified y * . Some contours of T * are given in figure 17b at y * = 12.5 (above) and y * = 40 (below). By comparing figures 17a and 17b, one can appreciate the correspondence between the values of R n-eq tot at these wall distances and the development of the heated turbulent structures as x/δ increases.

Development of turbulent heat transfer

The previous section has allowed us to appreciate and quantitatively define an equilibrium sub-layer, where both the mean temperature and the mean energy balance are equilibrated. In this section, we investigate its relevance with respect to the evolution along the isothermal wall of several turbulent quantities.

Figure 18 shows the evolution of the temperature fluctuations and of the wall-normal turbulent heat flux for different crosswise sections compared to equilibrium. Observe how both quantities progressively grow with x/δ, and how the peak shifts towards the higher y * , until it stabilises at y * ≈ 20 for the r.m.s temperature and y * ≈ 45 for the turbulent heat flux.

The wall-normal turbulent heat flux (see figure 18b) clearly exhibits an equilibrated behaviour in the near-wall region for every x/δ. Interestingly, the same cannot be said about the fluctuating temperature (see figure 18a), which starts agreeing with equilibrium only at a certain distance from the leading edge (x/δ ≈ 3/2π). This illustrates how the r.m.s. temperature has a stronger inertia, and needs a larger distance to settle, an aspect which has also been observed experimentally (see [START_REF] Teitel | A step change in wall heat flux in a turbulent channel flow[END_REF]. Also observe, in figure 8, the scant fluctuating structures near the leading edge. Figure 18b also outlines the fast wall-normal expansion of the turbulent heat flux during the early boundary layer development. As suggested in §4.2, this explains the fast decay of the region y * > h * n-eq seen in figure 16b.

The turbulent budgets of the enthalpy variance are now analysed. The time-averaged transport equation of h 2 reads (see, for a similar derivation, [START_REF] Vicquelin | Effects of radiation in turbulent channel flow: analysis of coupled direct numerical simulations[END_REF]

): 0 = - ∂ ∂x ρ u h 2 2 Ix - ∂ ∂y ρ v h 2 2 Iy - ∂ ∂x q cd x h IIx - ∂ ∂y q cd y h IIy + + q cd x ∂h ∂x + q cd y ∂h ∂y + q cd z ∂h ∂z III + h N IV -ρ u h ∂ h ∂x Vx -ρ v h ∂ h ∂y Vy + - 1 2 ∂ ∂x ρ u h h V Ix - 1 2 ∂ ∂y ρ v h h V Iy + h - ∂q cd x ∂x - ∂q cd y ∂y + N V II
, (4.10) composed of the streamwise and wall-normal mean convective terms (I x and I y ) and molecular diffusion (II x and II y ), molecular dissipation (III), a term of correlation between the enthalpy and N (IV ), where N = Dp Dt + τ ij ∂ui ∂xj + S ener regroups the compressibility effects and the source term, the streamwise and wall-normal production (V x and V y ) and turbulent diffusion terms (V I x and V I y ), and, finally, a term related to the enthalpy-density correlation (V II). Figure 19 shows the evolution of the different terms (scaled with respect to q w (x) 2 /µ w (x)) for different crosswise sections. For each x/δ, the values of h * eq and h * n-eq defined in §4.2 are put into evidence. Some terms (IV and V II), being always small and negligible, are not plotted. Finally, we have observed non-negligible levels of spurious numerical dissipation in our simulation, which has been added to the molecular dissipation (III) to yield a total dissipation term and a closed balance. The justification of this approach is given in Supplementary Materials B.

Figure 19a is relative to x/δ = 0.73. At this close distance from the leading edge, the different terms of (4.10) are active only in a small portion (approximately 10%) of the boundary layer, delimited by h * n-eq , beyond which the flow is still not affected by the heated wall. Observe how it is between h * eq and h * n-eq that the non-equilibrium effects are located, with non-negligible streamwise convective and production terms. The classic contributions (wall-normal production, molecular dissipation and wall-normal turbulent and molecular diffusion), instead, appear to be under-developed. This is true even between the wall and h * eq , which is in agreement with the low level of r.m.s. temperature observed near the leading edge.

Figure 19b is relative to x/δ = 2π. It is interesting to notice how all the streamwise non-equilibrium terms in the temperature variance budget have vanished. In the region between h * eq and h * n-eq , while the mean energy balance presents strong non-equilibrium terms as seen in §4.2, the enthalpy variance budget seems to merely feature the gradual development of the classic terms. Observe how in the equilibrium sub-layer, instead, the different contributions controlling the r.m.s. temperature have by now settled, as can be seen by comparing figures 19b, 19c and 19d for y * < h * eq . Accordingly, the r.m.s. temperature exhibits an equilibrated behaviour for y * < h * eq at these crosswise sections (see figure 18a).

Concerning the wall-normal turbulent heat flux, the time-averaged transport equation reads (see, for a similar derivation, [START_REF] Vicquelin | Effects of radiation in turbulent channel flow: analysis of coupled direct numerical simulations[END_REF]

): 0 = - ∂ ∂x ρ u v h Ix - ∂ ∂y ρ v v h Iy + ∂ ∂x τ 2,1 h -q cd x u IIx + ∂ ∂y τ 2,2 h -q cd y v IIy + -τ 2,1 ∂h ∂x + τ 2,2 ∂h ∂y + τ 2,3 ∂h ∂z + q cd x ∂v ∂x + q cd y ∂v ∂y + q cd z ∂v ∂z III -h ∂p ∂y IV + v N V + -ρ u h ∂ v ∂x -ρ u v ∂ h ∂x V Ix -ρ v h ∂ v ∂y -ρ v v ∂ h ∂y V Iy - ∂ ∂x ρ u v h V IIx - ∂ ∂y ρ v v h V IIy + + h - ∂p ∂y + ∂τ 2,1 ∂x + ∂τ 2,2 ∂y + v N - ∂q cd x ∂x - ∂q cd y ∂y V III , (4.11)
composed of the mean streamwise and wall-normal convective terms (I x and I y ) and molecular diffusion (II x and II y ), the molecular dissipation (III), the enthalpy-pressuregradient correlation (IV ), the correlation between the wall-normal velocity and N (V ), the streamwise and wall-normal production (V I x and V I y ) and turbulent diffusion (V II x and V II y ), and, finally, two terms related to the enthalpy-density and wall-normalvelocity-density correlations (V III). Figure 20 shows the evolution of the different terms (scaled with respect to q w (x)τ w (x)/µ w (x)) for different crosswise sections with the respective h * eq and h * n-eq defined in §4.2. The terms V and V III are always small and negligible, and are not plotted.

Most of the observations made concerning the budgets of enthalpy variance are still valid. The non-equilibrium region h * eq < y * < h * n-eq manifests itself with a gradual development of the classic terms (including, in this case, the enthalpy-pressure-gradient), rather than with the presence of strong streamwise terms as seen in §4.2 for the mean energy balance. In this case, the only streamwise terms that are important near the leading edge are the streamwise production and convection, as illustrated in figure 20a. Yet, even these terms become negligible at x/δ = 2π (see figure 20b).

The only difference with the enthalpy variance budget, is the behaviour for y * < h * eq . For the wall-normal turbulent heat flux, the classic budgets have settled in the equilibrium sub-layer even for the smallest x/δ, which is in agreement with the equilibrated profiles observed for y * < h * eq in figure 18b. Figure 21 shows the details of the budgets of enthalpy variance and wall-normal turbulent heat flux for y * ≈ h * n-eq at x/δ = 0.73, i.e., the same crosswise section of figures 19a and 20a. In both figures 21a and 21b, it can be seen how the turbulent diffusion plays an important role in the development of both h 2 and v h . In the case of the enthalpy variance, turbulent diffusion is actually the only term inducing a gain of The presence of non-negligible turbulent diffusion around y * ≈ h * n-eq corroborates the hypothesis formulated in §4.2, i.e., that the fast decay of the temperature-homogeneous region beyond y * > h * n-eq might be due to the efficient heat transport carried out by turbulence.

Finally, the behaviour of the turbulent Prandtl number Pr t is investigated. Pr t offers a simple way of relating the turbulent kinematic viscosity ν t and conductivity λ t in Reynolds-averaged Navier-Stokes (RANS), LES or wall-modelled LES (WMLES) turbulence modelling, and, analogously to the molecular Prandtl number, Pr t is defined as:

Pr t = ρν t c p λ t = -ρ u v -ρ v h /c p
dT /dy du/dy . (4.12)

Experimentally, [START_REF] Blom | Experimental determination of the turbulent prandtl number in a developing temperature boundary layer. 4th int[END_REF] analysed the evolution of the turbulent Prandtl number in a developing thermal boundary layer, and found that it varies in both wall-normal and streamwise directions, while suggesting that a universal distribution of Pr t can only be expected in the near-wall region. [START_REF] Antonia | Response of a turbulent boundary layer to a step change in surface heat flux[END_REF] carried out similar investigations, and found that the turbulent Prandtl number is generally in the order of magnitude of one, thus showing the overall validity of the Reynolds analogy; according to their experimental results, values of Pr t strongly exceeding unity can only be found near the leading edge of the thermal boundary layer. Figure 22 displays the evolution of Pr t for different x/δ, compared to its equilibrium profile. The latter, is represented by the dashed-dotted line. It can be seen that Pr t is slightly greater than unity at the wall, before decreasing and attaining Pr t ≈ 0.7 in the channel core flow. Its mild gradient along the wall-normal direction can justify the use of a constant average turbulent Prandtl number (usually Pr t ≈ 0.9), at least in equilibrium boundary layers (see, for instance, Kawai & Larsson 2012).

Figure 22a shows the evolution of Pr t for four crosswise sections which are relatively close to the leading edge. As can be seen, the turbulent Prandtl number does not exhibit an equilibrated behaviour in the near-wall region, where its value largely exceeds unity. Since the turbulent viscosity and conductivity, on their part, are at equilibrium in the wall vicinity (not shown here), this behaviour is explained by the abrupt increase of c p astride x/δ = 0, and represents another leading-edge effect. Further away from x/δ = 0, Pr t remains for the most part greater than unity and, as x/δ increases, its perturbation expands in the wall-normal direction, while its intensity decreases. This tendency is particularly evident in figure 22b, for x/δ ranging between [π, 1.75π]. The profiles of figure 22c, instead, illustrate how the turbulent Prandtl number essentially returns to equilibrium at x/δ ≈ 5π.

The larger values of Pr t during the development of the thermal boundary layer, are explained by the fact that its numerator, ν t , is already well-established with the incoming fully developed adiabatic flow; its denominator λ t , instead, is initially small, since heat transfer is essentially carried out by mean convection, as seen in §4.2. The results are, therefore, in qualitative agreement with those of [START_REF] Antonia | Response of a turbulent boundary layer to a step change in surface heat flux[END_REF]. On the other hand, there is no evidence of the existence of a universal distribution in the near-wall region, as indicated in [START_REF] Blom | Experimental determination of the turbulent prandtl number in a developing temperature boundary layer. 4th int[END_REF]. It appears to be quite the opposite since the highest values of Pr t can be found near the leading edge in the very proximity of the wall. Figure 22: Evolution of turbulent Prandtl number in the wall-normal direction for different x/δ: --present results at (from light to dark gray) x/δ = 0.18, 0.73 0.5π and 0.75π (a), at x/δ = π, 1.25π, 1.5π and 1.75π (b); at x/δ = 5.67π, 8.25π, 10.8π and 16.5π (c); -• -• -equilibrium profile.

Conclusion

The direct numerical simulation of a non-equilibrium turbulent thermal boundary layer is performed in a channel flow. Non-equilibrium is induced by a step-change in wall temperature, which perturbs a turbulent adiabatic flow initially at equilibrium. The flow later evolves towards a new equilibrium state, characterised by a fully developed thermal boundary layer. The simulation is performed under one single flow regime, in terms of Reynolds number, Mach number and heat flux parameter. The study throws light on several physical phenomena which are expected to appear in phenomenologically similar non-equilibrium flows where a streamwise and transverse development the thermal boundary layer follows a sudden perturbation.

The leading edge of the isothermal wall appears to be the most perturbed region of the flow, with effects that also slightly propagate upstream. The temperature stepchange causes a discontinuity of the wall shear stress and of the pressure gradient, which modify the velocity profiles. In our case, these modifications seem to be fully corrected by the adoption of a new wall-scaling, which takes into account both the mean fluid property variations and the pressure gradient. This suggests that the disturbance is confined in the very near-wall region, and does not propagate to the core of the flow. In the most general case, these effects are expected to appear at any discontinuity of thermal wall condition and, being regulated by the heat flux parameter B q , to expand and aggravate as B q increases, i.e., for stronger temperature gradients or lower friction velocities. In these conditions, one should anticipate larger wall-normal velocities and stronger pressure gradient discontinuities astride the leading edge, the effects of which might not be corrected by only adopting a different wall-scaling.

The study has allowed for the identification of an equilibrium sub-layer, where several quantities are equilibrated. These quantities are the mean temperature, the energy fluxes, the wall-normal turbulent heat flux and its budgets. Defined as the part of the flow where the equilibrium contributions to the energy balance are preponderant, the equilibrium region is initially confined in the conductive sub-layer and slowly driven by heat diffusion, before it reaches the most turbulent part of the flow. The existence of such an equilibrium layer is an important assumption for turbulence modelling. It may a priori justify, for instance, the use of equilibrium wall-models in the wall-modelled large-eddy simulations of non-canonic flows, as long as the modelled portion of the flow lays within the aforesaid equilibrium layer. Experimentally observed by several authors, there is reason to believe that it might be identifiable in any comparable non-equilibrium flow, where the same criterion on the energy balance contributions can be adopted to analyse it. However, stronger leading edge effects might initially disrupt its structure and delay its development.

Other quantities do not have the expected behaviour in the equilibrium sub-layer. Low levels of temperature fluctuations, for instance, have been observed in the early development of the thermal boundary layer, accompanied by under-developed budgets of enthalpy variance. This illustrates the stronger inertia of the r.m.s. temperature, and it is an aspect that has also been observed experimentally. Finally, the turbulent Prandtl number exceeds unity in a relatively large portion of the flow, which includes part of the equilibrium sub-layer. The perturbation is due to the fact that heat transfer is initially carried out by mean convection, and only to a minor extent by turbulence. An important role is also played by the abrupt variation of fluid properties at the beginning of the isothermal wall, thus representing another leading edge effect. Even in the present case, characterised by a mild temperature gradient, different models for the turbulent Prandtl number might be required to take into account its modified behaviour in RANS or LES modelling of turbulent heat transfer.

In perspective, further investigations are required to deepen our knowledge of heat transfer taking place in non-equilibrium turbulent flows. Different flow regimes, in terms of Reynolds number, Mach number and heat flux parameter, should allow for wider generalisations of the effects observed in this study, with a particular allusion to the leading-edge perturbation and the evolution of the equilibrium sub-layer. Entirely different non-equilibrium heat transfer cases, instead, could bring to light other physical phenomena that our flow configuration has not allowed us to observe. The predominant terms identified outside the equilibrium sub-layer can be different in another kind of flow topology.

Supplementary Materials A: shift of streamwise root-mean-square velocity near the leading edge of the isothermal wall

As seen in Section 3 of the article, an upward shift of streamwise root-mean-square (r.m.s.) velocity u rms is observed near the leading edge of the isothermal wall, i.e., for x/δ ≈ 0. This perturbation vanishes as x/δ increases, and we have attributed this anomaly to the recycling method described in §2.2. In fact, between the inlet (at x/δ = -4π) and the recycling plane (at x/δ = -2π), there is not a real periodicity in the sense of absence of boundary conditions, like in a bi-periodic channel flow or like along the spanwise direction in our simulation. In our case, we extract the fields at x/δ = -2π and re-inject them at the inlet using relaxed characteristic boundary conditions, and although we use an extremely high relaxation coefficient (∼ 10 +5 ), at x/δ = -4π we do not have the exact velocity values of x/δ = -2π.

Other aspects might contribute to the shift we observe on u rms , namely the discontinuity of surface temperature at x/δ, and potential auto-correlation issues between the inlet and the recycling plane. The objective of this document is to provide further elements to support the explanation given in Section 3.

We have performed two additional direct numerical simulations. The first simulation's domain is 6πδ long and the recycling is prescribed at x = 2πδ (i.e., exactly like the non-equilibrium simulation of our paper); the second simulation's domain is still 6πδ long, yet the recycling is prescribed at x = 4πδ, so that we are sure to avoid any autocorrelation problem. In both cases, the walls are adiabatic and the point distribution is uniform along the streamwise direction, thus avoiding any discontinuity in the streamwise direction. Table 1 summarises the main numerical parameters of the two simulations (S1 and S2).

In the following, (•) + denotes classic wall-scaling, as in the paper, and x/δ denotes the non-dimensional distance from the inlet (therefore, x/δ ∈ [0, 6π]); results are compared to the equilibrium adiabatic channel flow of Section 3.1. Let us commence by analysing the evolution of the mean streamwise velocity. Figure 1 shows the profiles obtained at x/δ = 0.5π (a), x/δ = 2.5π (b), x/δ = 4.5π (c) and x/δ = 5.5π (d) with both S1 and S2. First of all, notice that no appreciable difference can be observed between S1 and S2 at any x/δ, indicating that the recycling location has no influence on the mean velocity. As can be seen, the profiles at x/δ = 0.5π are strongly perturbed, with an upward shift of around 10%; however, the perturbation quickly decreases and, even if it is still noticeable at x/δ = 2.5π, it appears to have fully disappeared by x/δ = 4.5π. Note that, in the paper, the leading edge of the isothermal wall is at a distance of x/δ = 4π from the inlet, which explains why no perturbation of the mean streamwise velocity has been reported in Section 3.1. Now, let us focus on the r.m.s. velocity profiles. Figure 2 shows the streamwise, wallnormal and spanwise r.m.s. velocity profiles obtained at x/δ = 0.5π (a), x/δ = 2.5π (b), x/δ = 4.5π (c) and x/δ = 5.5π (d) with both S1 and S2. In this case, larger differences between S1 and S2 are observed. However, these deviations become remarkably smaller as x/δ increases. Therefore, also in this case we can conclude that the location of the recycling plane does not seem to have a strong influence on the flow statistics. The most perturbed profiles appear to be the ones at x/δ = 0.5π, i.e., the closest to the inlet, where none of the profiles agree with the equilibrium ones; this is true for the streamwise component (as observed in the paper), yet also for the wall-normal and spanwise components. As x/δ increases, v rms and w rms quickly return to equilibrium, which explains why no shifts have been reported for these components in Section 3.1. Concerning u rms , instead, notice how the perturbation persists even until x/δ = 5.5π, consistently with what has been observed in Section 3.1 of the paper. In this case, however, the channel flow's length is limited, and, therefore, the return to equilibrium of u rms does not occur before the outlet.

Our results can be summarised as follows:

-The mean and r.m.s. velocity profiles are perturbed by the recycling method; -The perturbation appears near the inlet and vanishes as x/δ increases; -The most persisting impact, as observed in the non-equilibrium simulation of the paper, seems to be the one on the streamwise r.m.s. velocity, which exhibits a shift even at x/δ = 5.5. These results allow us to conclude that:

-The perturbation cannot be attributed to the discontinuity of wall temperature since, in this case, the walls of the channel flow are entirely adiabatic;

-The perturbation cannot be ascribed to potential auto-correlation problems since the results of S1 and S2, obtained with two different recycling plane locations, are extremely similar;

-The perturbation is due to the recycling method itself and by the use of relaxed characteristic boundary conditions at the inlet towards target values determined at the recycling plane. Despite the very short response time, this introduces a non-ideal recycling, all the more so as negative streamwise velocity can occasionally be encountered, which the inlet characteristic boundary cannot handle properly. This conclusion, in particular, is corroborated by the fact that the most affected profiles are the ones at x/δ = 0.5π and the deviation decreases as the distance from the inlet is increased.

Figure 1 :

 1 Figure 1: Computational domain of the present study.

Figure 2 :

 2 Figure 2: Upstream equilibrium state. Mean profile of streamwise velocity (a): -present results with 2πδ-long bi-periodic channel flow (gray) and at x/δ ≈ 0 (black);results from[START_REF] Kawamura | Dns of turbulent heat transfer in channel flow with respect to reynolds and prandtl number effects[END_REF]. Profiles of r.m.s. streamwise, wall-normal and spanwise velocity respectively (b): --, -• -• -and ----present results (colors as above); , and results from[START_REF] Kawamura | Dns of turbulent heat transfer in channel flow with respect to reynolds and prandtl number effects[END_REF].

Figure 3 :Figure 4 :

 34 Figure 3: Downstream equilibrium state. Profiles of mean streamwise (a) and r.m.s. streamwise, wall-normal and spanwise velocities (b) at x/δ ≈ 18π. See figure 2 for a reference to colours and symbols.

Figure 5 :

 5 Figure 5: Nusselt number Nu (a) and heat flux parameter B q (b) as a function of x/δ: --present results; -• -• -equilibrium values of table 3.

Figure 6 :

 6 Figure 6: Leading edge effects. Contours of wall shear stress (τ w ) near the leading edge of the isothermal wall. Values are non-dimensionalised with respect to τ w , averaged over the surface of the wall.

Figure 7 :

 7 Figure 7: Instantaneous iso-temperature θ = T -T in / T w -T in = 0.9 along the channel, coloured by the friction velocity u * .

Figure 8 :Figure 9 :

 89 Figure 8: Instantaneous iso-surface of fluctuating temperature (T * ) = T * -T * /T τ = 2, coloured by the friction velocity u * .

Figure 11 :

 11 Figure 10: Wall-normal velocity profile for different x/δ in semi-local scaling (a) and normalised with the mean streamwise velocity (b): --present results at x/δ = 0.015, 0.18, 0.73 (from lightest to darkest); -• -• -equilibrium profile.

Figure 12 :

 12 Figure 12: Mean streamwise velocity profile for different x/δ in the newly introduced wall scaling (a). See figure 9 for a reference to lines and colours. Profiles of r.m.s. streamwise, wall-normal and spanwise velocity, respectively (b): --, -• -• -and ----present results at x/δ = -0.06, 0.015, 0.09, 0.18 and equilibrium (see §3).

Figure 13 :

 13 Figure 13: Momentum flux balance at x/δ = -0.06 (a) and 0.18 (b): --pressure gradient; ----streamwise (gray) and wall-normal (black) viscous terms; • • • • • • streamwise (gray) and wall-normal (black) mean convective terms; -• -• -streamwise (gray) and wall-normal (black) turbulent terms.

Figure 14 :

 14 Figure 14: Profiles of mean temperature for different x/δ: --present results at x/δ = 0.015 -0.21 -1 2 π -3 2 π -3.6π -7.75π -17.9π (from lighter to darker); -• -• -equilibrium profile of §3.2.

Figure 15 :

 15 Figure 15: Energy flux balance at x/δ = 0.73 (a), 2π (b), 5.7π (c) and 16.5π (d): -source term; ----streamwise (gray) and wall-normal (black) conductive fluxes; • • • • • • streamwise (gray) and wall-normal (black) convective fluxes; -• -• -streamwise (gray) and wall-normal (black) turbulent fluxes; -• • -compressibility effects.

Figure 17 :

 17 Figure 17: Evolution of R n-eq tot as a function of x/δ for several wall distances: --R n-eq tot for y * = 5, 12.5, 30 and 40 (from light gray to black respectively); -• -• -isoline R n-eq tot = 0.1 (a). Instantaneous contours of temperature at y * = 12.5 (above) and y * = 40 (below) (b).

Figure 18 :

 18 Figure 18: Profiles of r.m.s temperature (a) and wall-normal turbulent heat flux (b) for different x/δ: --present results at x/δ = 0.015 -0.21 -1 2 π -3 2 π -3.6π -7.75π -17.9π (from lighter to darker); -• -• -equilibrium profiles of §3.2.

Figure 19 :

 19 Figure 19: Budgets of enthalpy variance at x/δ = 0.73 (a), 2π (b), 5.7π (c) and 16.5π : --streamwise (V x : gray) and wall-normal (V y : black) production; -• • -• • -molecular dissipation III; -----streamwise (II x : gray) and wall-normal (II y : black) molecular diffusion; -• -• -streamwise (V I x : gray) and wall-normal (V I y : black) turbulent diffusion; • • • • • • streamwise (I x : gray) and wall-normal (I y : black) mean convective terms.

Figure 20 :

 20 Figure 20: Budgets of wall-normal turbulent heat flux at x/δ = 0.73 (a), 2π (b), 5.7π (c) and 16.5π: --streamwise (V I x : gray) and wall-normal (V I y : black) production; -• • -• • -molecular dissipation III; -----streamwise (II x : gray) and wall-normal (II y : black) molecular diffusion; -• -• -streamwise (V II x : gray) and wall-normal (V II y : black) turbulent diffusion; • • • • • • streamwise (I x : gray) and wall-normal (I y : black) convective terms; ----enthalpy-pressure-gradient correlation IV .

Figure 21 :

 21 Figure 21: of enthalpy variance (a) and wall-normal turbulent heat flux (b) at x/δ = 0.73 zoomed at y * ∼ h * n-eq . See figures 19 and 20 for the notation.

Figure 1 :

 1 Figure 1: Mean streamwise velocity profiles at x/δ = 0.5π (a), x/δ = 2.5π (b), x/δ = 4.5π (c) and x/δ = 5.5π (d). Black solid line, S1; Gray solid line, S2; bi-periodic adiabatic channel flow from Section 3.1.

Figure 2 :

 2 Figure 2: Profiles of r.m.s. streamwise, wall-normal and spanwise velocity, respectively, at x/δ = 0.5π (a), x/δ = 2.5π (b), x/δ = 4.5π (c) and x/δ = 5.5π (d): --, -• -•and ----present results (black lines for S1, gray lines for S2); , and results from bi-periodic adiabatic channel flow of Section 3.1.

Figure 1 :

 1 Figure1: Budgets of k. Production (red), dissipation (blue), molecular diffusion (green), turbulent diffusion (yellow), sum of terms (black). M1 (dashed-dotted lines), M2 (solid lines).

Table 1 :

 1 Size, number of nodes and resolutions of S1 and S2.

	Size	Number of nodes Recycling location
	S1 6πδ, 2δ, πδ	598, 179, 200	x = 2πδ
	S2 6πδ, 2δ, πδ	598, 179, 200	x = 4πδ

Table 1 :

 1 Size, number of nodes and resolutions of M1 and M2.

	Size	Number of nodes ∆X + ∆Z +	∆Y +
	M1 100, 179, 100	1.79M	11.3 5.65 0.34-3.1
	M2 256, 192, 256	12.5M	4.42 2.2 0.34-2.23

Acknowledgements

The numerical study was carried out thanks to the HPC resources of TTGC under the allocation 2018-A0042B10159 made by GENCI, and the HPC resources from the "Mésocentre" computing centre of Centrale-Supélec and École Normale Supérieure Paris-

Saclay supported by CNRS and grants from Ministry of Research via the Contrat Plan État-Région.

Declaration of interests

The authors report no conflict of interest.

Supplementary Materials B: molecular and numerical dissipation in root-mean-square budgets

As mentioned in Section 4.3, we have observed non-negligible levels of spurious numerical dissipation with the CFD solver that we have retained for the study, i.e., the parallel code AVBP (Schonfeld & Rudgyard 1999;Moureau et al. 2005), which employs a time-explicit finite-element two-step Taylor-Galerkin scheme (Colin & Rudgyard 2000) which provides third-order accuracy in space and time.

AVBP is commonly used for unstructured meshes on very complex industrial applications. The numerical scheme relies on a distributed-residual scheme in a cell-vertex approach, which is less-widespread and makes delicate the computation of the different terms in the various root-mean-square budgets. Furthermore, since we have preferred to post-treat the different terms with standard spatial discretisation on the structured grid, our post-processing includes a possible discretisation error with the actual terms seen by the numerical scheme. These errors lead to an underestimated level of molecular dissipation, and yield unclosed root-mean-square balances. We have observed this issue on both budgets of enthalpy h h and turbulent kinetic energy k, and we have further investigated. As we explain hereinafter, we have identified that the molecular dissipation could be safely estimated as the complementary term to all the others, so that the sum of all terms gives exactly zero, yielding a closed balance.

We have performed the DNS of an adiabatic 2πδ channel flow at Re τ = 180 with two meshes. The former has similar resolution to the one used in our work, and in agreement with the usual resolutions of the literature, in particular Kawamura et al. (2000); we will call it M1. The latter, instead, has been extremely (and excessively) refined: 256 × 192 × 256 points in the streamwise, wall-normal and spanwise directions, respectively, for a total of ∼ 12.5M nodes; we will call it M2. More information is given in table 1.

We have extracted the budgets of k in both cases, and obtained the results shown in figure 1. As one can see, M1 does not allow us to close the balance when each term is computed according to its definition. With M2, the balance is, on the other hand, perfectly closed. Observe that the only term which strongly varies between the two is the computed dissipation; other than that, there are negligible discrepancies near the wall for the molecular diffusion and near the peak for the production. Hence, the imbalance seen with M1 is largely retrieved within the molecular dissipation computed on M2. This proves that our problem is due to the computed molecular dissipation which (i) can be degraded in our post-processing of the cell-vertex scheme and (ii) should be augmented with the residual numerical dissipation.

Since the resolution of M2 would be infeasible in our simulation's domain (which, recall, is 22πδ-long in the streamwise direction), we have retained a pragmatic determination of dissipation in the present study: we have added the imbalance to the computed physical dissipation to yield a total dissipation. Concerning the wall-normal turbulent heat flux v h , we have encountered much fewer issues, since molecular dissipation plays a less significant role and the dominating loss term is h ∂p /∂y, which is well predicted.