A Wall-Adapted Anisotropic Heat Flux Model for Large Eddy Simulations of Complex Turbulent Thermal Flows

In this paper, a wall-adapted anisotropic heat flux model for large eddy simulations of complex engineering applications is proposed. First, the accuracy and physical consistency of the novel heat flux model are testified for turbulent heated channel flows with different fluid properties by comparing with conventional isotropic models. Then, the performance of the model is evaluated in case of more complex heat and fluid flow situations that are in particular relevant for internal combustion engines and engine exhaust systems. For this purpose large eddy simulations of a strongly heated pipe flow, a turbulent inclined jet impinging on a heated solid surface and a backward-facing step flow with heated walls were carried out. It turned out that the proposed heat flux model has the following advantages over existing model formulations: (1) it accounts for variable fluid properties and anisotropic effects in the unresolved temperature scales, (2) no ad-hoc treatments or dynamic procedure are required to obtain the correct near-wall behavior, (3) the formulation is consistent with the second law of thermodynamics, and (4) the model has a similar prediction accuracy and computational effort than conventional isotropic models. In particular, it is shown that the proposed heat flux model is the only model under consideration that is able to predict the direction of subgrid-scale heat fluxes correctly, also under realistic heat and fluid flow conditions in complex engineering applications.


Introduction
Many energy systems, such as internal combustion engines or exhaust after treatment devices, are confined by solid walls. Thereby, combustion and other energy conversion processes lead to an intensive heat transfer through the solid walls along with complex unsteady mixing dynamics that determine the lifetime and overall thermodynamic performance of thermo-fluid systems. These complex heat and fluid flow phenomena make a description and optimization of such energy systems very challenging and require accurate and viable tools for their engineering design.
On the numerical side, large eddy simulation (LES) technique has proved to be a promising approach to predict complex heat and fluid flow phenomena in many thermo-fluid systems, likewise in automotive applications Hasse et al. 2010;Nishad et al. 2012;Goryntsev et al. 2014;Rutland 2011;Goryntsev et al. 2010). In LES, the large three-dimensional unsteady turbulent motions are explicitly computed, whilst a turbulence closure model accounts for the influence of the unresolved more universal scales (Klein 2008;Sagaut 2006;Pope 2009). The benefit of such an approach is quite obvious. First, the computational expense of LES is significantly lower than in fully resolved direct numerical simulations (DNS). Secondly, only small scale turbulent structures with a small amount of turbulent energy have to be modeled, which are believed to be universal, homogeneous and isotropic. This simplifies the turbulence modeling, improves the predictive capability compared to approaches based on the solution of the Reynolds-averaged Navier-Stokes equations (RANS) and makes LES valid for a broad range of flow situations with complex physics.
In principle, most strategies that are used to close the LES momentum equation can be also applied to deal with the unresolved subgrid-scale heat flux vector in the energy equation (Sagaut 2006). However, based on turbulence theory, the rationale behind LES is less obvious when dealing with turbulent heat transport. The main reason is that the temperature variance spectrum as well as the dynamics of small temperature scales have a less universal character than the kinetic energy spectrum and velocity scales. As a consequence, the spectral scalar transfer across the LES cutoff can strongly depend on the physical regime in which it is located (Sagaut 2006). Moreover, small temperature scales are also influenced by the interaction of the velocity gradient and the scalar fluctuations, which causes anisotropic behavior even at smallest temperature scales. Consequently, more advanced subgrid-scale models are required in the case of turbulent heat transport in order to justify the cut off and modeling of the unresolved temperature scales. However, despite of the anisotropic behavior of small turbulent thermal structures, isotropic models are most often employed in LES of turbulent heat transport. Thereby, the subgrid-scale thermal diffusivity is traditionally represented based on the Reynolds analogy and the concept of turbulent Prandtl number. Many researchers intended to improve this simple approach by using a dynamic procedure to calculate the turbulent Prandtl number (Lilly 1992;Moin et al. 1991), including buoyancy effects (Eidson 1985), using a definition of the thermal diffusivity based on Kolmogorov scaling (Wong and Lilly 1995) or including the effects of local fluid properties in the subgrid-scale model (Otic 2010).
As implied by the discussion above, a better representation of the subgrid-scale heat flux vector for complex heat and fluid flow situations can be obtained by accounting for the anisotropic behavior of small temperature scales and introducing a tensor subgrid-scale thermal diffusivity. Some of these models were derived in analogy to the general gradient diffusion hypothesis (Daly and Harlow 1970) as it is often applied in the RANS context. In the approach of Peng and Davidson (2002), the authors developed a non-linear subgrid-scale heat flux model based on considerations of the transport equation of the subgrid-scale heat flux. Wang et al. (2006Wang et al. ( , 2007Wang et al. ( , 2008 proposed a series of models which includes the resolved strain-rate tensor S ij , rotation rate tensor Ω ij and the temperature gradient T∕ x i . Similar, Rasam et al. (2017) proposed an anisotropic scalar flux model that is based on the modeled transport equation of the subgrid-scale scalar flux and depends also on S ij , Ω ij and T∕ x i . In contrast to these model formulations, the anisotropic model by Huai and Sadiki Huai (2006);Pantangi et al. (2010) is based on the second law of thermodynamics in conjunction with the invariant theory. In this way, the irreversibility requirements of the second law of thermodynamics are automatically fulfilled by the model formulation. Besides, scale similarity and mixed models were also proposed in the literature, e.g. Jaberi and Colucci (2003); Porté-Agel et al. (2001); Salvetti and Banerjee (1995).
From this brief literature review it appears that only a few advanced anisotropic heat flux models exist in the literature and they are almost never used in LES practice. Most of these models are relatively complex, do not account for the effects of local fluid properties and do not provide the correct asymptotic near-wall behavior. To overcome these limitations, a novel wall-adapted anisotropic heat flux model for LES is proposed in this paper. The main features of the novel model are that (1) it accounts for variable fluid properties and anisotropic effects in the unresolved temperature scales, (2) no ad-hoc treatments or dynamic procedure are required to obtain the correct near-wall behavior, and (3) the formulation is consistent with the second law of thermodynamics. First, the prediction accuracy and computational cost of the novel heat flux model are evaluated in the present work for a turbulent heated channel flow with different Prandtl numbers. Then, the performance of the proposed heat flux model is demonstrated for complex heat and fluid flow situations.
The present paper is organized as follows: First, the applied LES approach and the novel anisotropic heat flux model are introduced (Sect. 2). Then, the accuracy and consistency of the proposed model is evaluated for turbulent heated channel flows at different Prandtl numbers (Sect. 3). Thereby comparison with available DNS data and with achievements obtained by means of conventional isotropic models are presented and discussed. Subsequently, the model is applied to more complex heat and fluid flow situations that are in particular relevant to internal combustion engines and exhaust gas systems, namely, a strongly heated pipe flow, a turbulent inclined jet impinging on a heated surface and a backward-facing step flow with heated walls (Sect. 4). Finally, some concluding remarks are provided at the end (Sect. 5).

Model Description and Numerical Treatment
In the case of LES with implicit filtering of an incompressible fluid flow with variable physical properties, the transport equations of mass, momentum and enthalpy can be formulated as Sagaut (2006) where (.) are spatially filtered quantities, ( .) Favre-filtered quantities and (.) sgs denotes the subgrid-scale quantities. is the mass density, U i the velocity, h the sensible enthalpy, ij the stress tensor, sgs ij the subgrid-scale stress tensor, and q sgs j the unresolved heat flux. In the case of a Navier-Stokes-Fourier fluid, the stress tensor is described as where is the molecular viscosity and p is the pressure. The resolved heat flux is expressed by means of Fourier's law for incompressible flow as where is the thermal conductivity and c p the isobaric heat capacity. The gravitational force and the radiation are not considered in Eqs. 2 and 3. Furthermore, the unsteady pressure term in the enthalpy equation is usually assumed to be small in incompressible flows (Poinsot and Veynante 2017). It is therefore neglected in this work. In order to close the filtered momentum and energy transport equations, the subgridscale stress tensor sgs ij = Ũ i U j −Ũ iŨj and the subgrid-scale heat flux vector q sgs i = Ũ i h −Ũ ih have to be postulated by being related to the resolved velocity and temperature fields, respectively, in terms of subgrid-scale models.

Modeling of the Subgrid-Scale Momentum Transport
In this work, the eddy viscosity approach is applied to model the subgrid-scale transport of momentum. Thereby, the isotropic part of sgs ij is included into the modified filtered pressure and the remaining deviatoric part is expressed using the eddy viscosity sgs and the Boussinesq approximation as The eddy viscosity is represented by means of the SIGMA-model (Nicoud et al. 2011) with standard model coefficient C = 1.5.

Modeling of the Subgrid-Scale Heat Transport
Because turbulent flows are thermodynamic processes and the directions of all such processes are restricted by the second law of thermodynamics, it is recommendable to take account of this fact at every level and kind of closure formulations. In analogy to the anisotropic scalar flux model of Huai (2006) and Pantangi et al. (2010), the explicit algebraic anisotropic heat flux model proposed here is based on the second law of thermodynamics in conjunction with the invariant theory. From this formalism, a general expression for the subgrid-scale heat flux vector in non-rotating observer system can be written as (see Huai 2006;Pantangi et al. 2010): where the first term on the right-hand side represents the contribution by linear diffusion, the second term expresses the influence of the subgrid-scale flow fluctuations acting on the resolved enthalpy gradient and the last term is associated with production/dissipation of subgrid-scale variance. For more details on the derivation of Eq. 8 the reader is referred to Ahmadi et al. (2006) and Sadiki and Hutter (2000).
Restricting ourselves to linear subgrid-scale viscosity models for the velocity field in Eq. 8 to represent sgs ij and by means of dimensional analysis, Eq. 8 can be reformulated as where k sgs is the unresolved turbulent kinetic energy, k sgs the dissipation rate of k sgs , sgs the subgrid-scale temperature variance and sgs its dissipation rate. c is a characteristic subgrid-scale time scale that can be either represented by means of a mechanical c = k sgs ∕ k sgs , a thermal c = sgs ∕ sgs or a mixed subgrid-scale time scale c = √ sgs k sgs k sgs sgs . The latter one was used for example in the model formulation proposed by the first author in Ries (2019). As shown by DNS (Kawamura et al. 1998), the time scale ratio between thermal and mechanical time scales behaves proportionally to ∼ √ Pr for turbulent fluid flows with a wide range of molecular Prandtl numbers. In order to account for fluid flows with variable molecular Prandtl numbers, the same behavior is assumed for the characteristic subgridscale time scale in the present model formulation leading to c ∼ √ Pr k sgs k sgs .
At this stage, it is obvious that closure relation for k sgs , k sgs , sgs and sgs are required to close the model formulation in Eq. 9. However, by applying the inertial-convective subrange theory (Obukhov 1968;Corrsin 1951;Schmidt and Schumann 1989) and combining all isotropic terms on the right-hand side, Eq. 9 can be rewritten in a fully algebraic form as where the coefficients C I and C II are determined in analogy to the generalized gradient diffusion hypothesis (Daly and Harlow 1970) and in accordance with the inertial-convective subrange theory (Obukhov 1968;Corrsin 1951;Schmidt and Schumann 1989) as C I = 2.381 and C II = 0.081.
Finally, in order to obtain the correct asymptotic behavior of q sgs i ∼ O y 3 in the vicinity of the wall, a damping factor f q is added to the second term on the right-hand side of Eq. 10. This leads to the final formulation of the present explicit algebraic anisotropic heat flux model as where f q is defined as which accounts for shear damping effects in the near-wall region. From scaling analysis, it appears that f q ∼ O y 3 leads to the correct asymptotic behavior of q sgs i ∼ O y 3 for y → 0 in the case that sgs ∼ O y 3 . Notice that the asymptotic behavior of sgs switches from cubic to quadratic when density fluctuations are present at the wall, which is correctly represented by the SIGMA-model (Nicoud et al. 2011). This holds also for q sgs i in case the SIGMA-model is used as subgrid-scale closure for the flow field. The coefficients of the novel model are summarized as C I = 2.381 , C II = 0.081 , C k = 0.094 and C = 1.048 . A detailed derivation of the model can be found in Ries (2019). An analysis of the asymptotic behavior of the proposed model for different spatial resolutions is provided in the Appendix.

Numerical Treatment
All simulations were carried out using OpenFOAM 2.4.0 (http://foam.sourc eforg e.net/ docs/Guide s-a4/Progr ammer sGuid e.pdf). Thereby, a low Mach-number approach was employed that is suitable for flows under incompressible conditions ( Ma < 0.3 ). This consists of a merged PISO (Issa 1985)-SIMPLE (Patankar and Spalding 1972) algorithm for the pressure-velocity coupling along with a second-order implicit backward-differencing scheme for time integration. The solution procedure was applied with a low-dissipative second-order flux-limiting differencing scheme for the convection terms and a conservative scheme for the Laplacian and gradient terms. A detailed description, verification and validation of the code employed in this study can be found in Ries et al. (2017); Ries (2019).

Model Evaluation
The evaluation study of the proposed anisotropic heat flux model (see Eq. 11) is divided into three parts. First, the physical consistency of the model is testified for turbulent heated channel flow at Re = 395 and Pr = 0.71. For comparison purpose, the DNS dataset of Kawamura et al. (1999) is utilized. Additionally, results of the linear thermal diffusivity model with a constant subgrid-scale Prandtl number of Pr sgs = 0.7 and with a dynamic procedure to calculate Pr sgs are also provided. In the second part, the influence of the fluid properties on the prediction accuracy of the heat flux models is analyzed. For this purpose LES of turbulent heated channel flow at different molecular Prandtl numbers were carried out and predictions are compared and assessed with the reference DNS dataset of Abe et al. (2004). Finally, in the last part of the model evaluation study, the prediction accuracy and computational cost of the proposed heat flux model are quantified and compared with those of isotropic heat flux models.
An illustration of the computational domain employed in the present evaluation studies is depicted in Fig. 1, where denotes half the height of the channel and N 1,2,3 is the number of grid points in x-, y-and z-direction.
Three numerical grids with (N 1 × N 2 × N 3 ) = (81 × 91 × 81) , (97 × 111 × 97) , and (121 × 137 × 121) cells, denoted here as grid no. 1, 2, 3, respectively, are used in the evaluation study. Thereby, all grids are refined towards the near wall region to ensure a non-dimensional wall distance much smaller than one. Figure 2 shows non-dimensional mean and rms temperature profiles predicted by the different heat flux models in comparison with DNS data. Thereby, the non-dimensional temperature is defined as

Physical Consistency
is the friction temperature taken from the DNS. As focus is put on thermal quantities only, velocity fields are not shown.
As it is expected, predictions of mean and rms temperature profiles become more accurate with increasing spatial resolution. This holds true for all heat flux models under consideration. Thereby, LES results of the different models are very close to each other which indicates a similar prediction accuracy of the models.
Next the physical consistency of predicted subgrid-scale heat fluxes are analyzed. For this purpose, Fig. 3 shows predicted wall-normal (a) and axial (b) turbulent heat fluxes in comparison with the DNS data. Solid lines denote resolved heat fluxes, while dashed lines represent subgrid-scale heat fluxes. Thereby, in case of turbulent heated channel flow, axial velocity u and temperature are always greater than or equal to zero. From this reason and by using the Steiner translation theorem, it follows directly that resolved and subgrid-scale axial heat fluxes should be also greater than or equal to zero.
Regarding wall-normal heat fluxes as depicted in Fig. 3a, it can be clearly seen that predictions of all subgrid-scale heat flux models are quite similar. In contrast, resolved and subgrid-scale axial heat fluxes have only same positive sign in case of the anisotropic heat flux model (see Fig. 3b). This reflects the physical consistency of the proposed anisotropic heat flux model. Considering the isotropic models, both, the standard and dynamic thermal diffusivity models are unable to reproduce the direction of axial subgrid-scale heat flux correctly. Nevertheless, although the resolved and subgrid-scale axial heat fluxes are not correctly predicted by means of the isotropic models, the additional anisotropic contribution appears quite small, which might explain that the overall prediction accuracy of thermal statistics with such isotropic models is often comparable to anisotropic model. This will be quantified in Sect. 3.3.

Influence of the Molecular Prandtl Number
The influence of the fluid properties on the prediction accuracy of the subgrid-scale heat flux models is analyzed next. For this purpose, Fig. 4 shows LES results of a turbulent heated channel flow at Re = 395 and molecular Prandtl numbers of Pr = 0.025, 0.71 and Fig. 1 Computational domain for the LES study of heated channel flow. N 1,2,3 represents the number of grid points 10. Thereby, predictions: The Sigma model has been implemented in a dedicated LES/DNS code and is tested against Direct Numerical Simulation (DNS) data from channel flow and grid turbulence data obtained from DNS and from measurements in a TOJ configuration. The flow through the turbulent of the proposed anisotropic heat flux model, the isotropic heat flux model with Pr sgs = 0.7 and the isotropic heat flux model using a dynamic procedure to calculate Pr sgs are compared to DNS data of Kawamura et al. (1999). In addition, results of the anisotropic heat flux model without the anisotropic term ( C II = 0 in Eq. 11) are provided in Fig. 4 in order to assess the advantage of the proposed √ Pr dependency of the model. Figure 4a depicts predicted mean temperature profiles and Fig. 4b shows rms temperature profiles as a function of non-dimensional wall distance. Thereby, only LES results of grid no. 1 are presented. Similar results are obtained for grid no. 2 and 3, and are therefore simply omitted for the sake of clarity.
As it can be clearly seen in Fig. 4, the physics of turbulent heat transfer in channel flow differs significantly for different Prandtl numbers. Non-dimensional mean and rms temperatures ( Θ + and Θ + rms ) increase with increasing molecular Prandtl number and peak values of Θ + rms are shifted towards the wall. This tendency is well retrieved by all tested subgridscale heat flux models, and, just as in the case of Pr = 0.71, predictions of the different models are also quite similar. However, it can be seen that predicted mean temperature profiles and peak values of Θ + rms at Pr = 0.025 and Pr = 10 are closer to the reference DNS data when the proposed anisotropic model is applied in comparison to cases conventional isotropic models are used. Similar results are obtained by using the proposed heat flux model without the anisotropic contribution term ( C II = 0 ). This suggests that mainly the √ Pr dependency, that reduces the residual contribution in case of fluid flows with Prandtl numbers smaller than one and increases the modeling contribution regarding fluid flows with Prandtl numbers higher than one, leads to a better prediction of flows with variable fluid properties.

Prediction Accuracy and Computational Cost
It appears that LES predictions of first and second order thermal statistics in turbulent heated channel flows are very similar for anisotropic and isotropic heat flux models. This is quantified next by means of an error analysis. Thereby, a turbulent heated channel flow test case at Re = 395 and Pr = 0.71 is selected for the error analysis and the DNS dataset of Kawamura et al. (1999) is utilized as reference. The normalized relative error of the mean  Kawamura et al. (1999) ( e + ) and rms ( e + rms ) temperatures with respect to the non-dimensional wall distance y + are shown in Fig. 5. Results are shown for grid no. 1. Similar error characteristics are obtained for the other spatial resolutions. Notice that the relative error is normalized in this study by means of the difference between the maximal and minimal value of the reference data, following the error analysis procedure for LES described in Ries et al. (2018c).
It is visible in Fig. 5 that errors are small in the near-wall region, increase rapidly in the buffer layer and finally decrease in the outer region. This trend is similar for all tested LES heat flux models and also for both thermal statistics, mean and rms temperatures. Especially in regard to e + rms , it is interesting to observe that the isotropic model with Pr sgs = 0.7 and in particular the proposed anisotropic heat flux model are more accurate in the nearwall region than the isotropic heat flux model with dynamic procedure. Further away from the wall, in the outer region, the error contribution is similar for all models. As shown by the authors in Ries et al. (2018c), localized dynamic procedures can produce a non-physical amount of residual contribution in the near-wall region when the dynamic procedure is not applied over homogeneous planes parallel to the walls. Such a procedure is generally not feasible in complex geometries and therefore not applied in this study. Thus, it seems to be likely, that this physical inconsistency of the dynamic procedure is responsible for the higher errors of the dynamic model in the near-wall region.
After analyzing the error characteristics of the different heat flux models for a given spatial resolution, the overall prediction accuracy of the models is now examined with respect to spatial resolution. Following the procedure described in Ries et al. (2018c), the normalized mean absolute error (nMAE) is employed as global error metric to quantify the overall prediction accuracy of each model. Thereby, the locations at which the nMAEs are computed are logarithmically distributed along the channel height in order to obtain an approximately equal number of sampling points in each flow regime. Results with respect to the spatial averaged ratio of Kolmogorov length scale K and grid width Δ grid are depicted in Fig. 6.
As it is expected the nMAEs decrease with increasing spatial resolution ( K ∕Δ grid ↑ ). This holds true for all tested heat flux models and also for both statistics, which confirms the consistency of all these modeling approaches in terms of LES. Thereby, in particular the anisotropic heat flux model and the isotropic model with constant Pr sgs have lowest values of nMAEs, reflecting a smaller modeling error and therefore best prediction accuracy. However, deviations are small, and it can be concluded that all the tested models have a Finally, the required computational cost of the subgrid-scale heat flux models is analyzed. Following the procedure described in Ries et al. (2018c), the relative computational cost of a subgrid-scale model CPUh * can be defined as the ratio of the CPU time spent for the calculation of the subgrid-scale model and the total computation time of the simulation. Regarding the selected subgrid-scale heat flux models in this study, it is found that the eddy diffusivity model with constant Pr sgs has the lowest relative computational cost ( CPUh * ∼ 0.25% ), whereas the eddy diffusivity model using a dynamic procedure to calculate Pr sgs is the most expensive one ( CPUh * ∼ 3% ). The proposed anisotropic heat flux model is also quite inexpensive ( CPUh * ∼ 0.5% ) since it is fully algebraic and does not use any dynamic procedure. Nevertheless, the CPU time spent for the calculation of all subgrid-scale heat flux models under consideration is fairly small compared to the total computation time of the simulations.
Considering the present evaluation study for turbulent heated channel flow, it turned out that the proposed anisotropic subgrid-scale heat flux model as well as isotropic models are able to predict first and second order thermal statistics accurately for this test case, regardless a dynamic procedure is used or not. However, only the proposed anisotropic model is able to reproduce the correct direction of the axial subgrid-scale heat flux, accounts for variable fluid properties and exhibits the proper near-wall behavior. This reflects the physical consistency of the proposed model. Furthermore, it is shown that the anisotropic heat flux model has not major impact on computational cost. Its prediction capability in complex flows is therefore demonstrated in the following section.

Application to Flow Configurations Relevant to Internal Combustion Engines and Exhaust Gas Systems
Among various energy systems, internal combustion (IC) engines features very complex heat and fluid flow situations. Besides being confined by solid walls, the IC engine motor is connected to an exhaust gas system. All together includes processes like (1) thermo-viscous boundary layer flows, (2) impinging cooling/heating, (3) recirculation, (4) flow separation and many more. Thus, given the complexity of heat and fluid flows in IC-engines (a) (b) Fig. 6 Normalized mean absolute error (nMAE) of predicted mean and rms temperatures with respect to spatial resolutions (legend see Fig. 5) and exhaust gas systems, it is useful to divide the evolving flow and mixing phenomena into different canonical flow situations that represents the most of the physical processes relevant to such applications. This allows to evaluate a modeling approach for such complex engineering systems in generality by considering only process relevant unit problems.
After analyzing the performance of the proposed anisotropic heat flux model for turbulent heated channel flow and different fluid properties, the novel approach is now applied to more complex heat and fluid flow situations, namely a strongly heated turbulent air flow in a pipe, a turbulent inclined jet impinging on a heated surface, and a backward-facing step flow with heated walls. These test cases are selected since they feature essential heat and fluid flow situations that are in particular relevant for internal combustion engines and exhaust gas systems. An illustration of the generic test cases and the location where such flow situations can be found in automotive technologies is shown in Fig. 7. A short description of each test case and the obtained LES results are presented and discussed in the following.

Strongly Heated Turbulent Air Flow in a Pipe
In exhaust gas systems and many other engineering applications, large temperature differences occur that leads to strongly varying thermo-fluid properties. In order to establish the validity of the present anisotropic heat flux model for such extreme operating conditions, LES of a strongly heated air flow in a vertical pipe with constant heat flux have been carried out and simulation results are compared with measurements of Shehata and McEligot (1998) and DNS data of Bae et al. (2006). In addition, LES results using the isotropic linear thermal diffusivity model with Pr sgs = 0.7 are also provided for comparison purpose. An illustration of the strongly heated pipe flow configuration is shown in Fig. 8, where D denotes the inner diameter of the pipe.
In the test section, a fully developed turbulent flow of dry air ( Re = 6000 , T 0 = 298.15 K, p = 0.1 MPa) enters a DN-25 pipe ( D = 0.0272 m, L = 30 D) and is heated up after an entrance length of 5D. The heated pipe region has a length of 25D with a constant wall heat flux of q w = 4.11kW∕m 2 . In line with the DNS study of Bae et al. (2006), air is treated in the current LES study as an ideal gas using the ideal gas equation. Other thermo-physical properties are obtained by means of power laws in the temperature as described in Bae et al. (2006). A block-structured numerical grid with 649536 control volumes is used to discretize the pipe flow domain. Thereby, the near-wall region is refined in order to fully resolve the small turbulence scales in the vicinity of the wall. At the pipe wall, a no-slip condition is set for the velocity and a zero Neumann conditions for the pressure. A constant wall heat flux of q w = c p h r | |r=R = 4.11kW∕m 2 is imposed at the heated wall while a zero temperature gradient condition is set at the adiabatic wall. In order to obtain realistic inflow turbulence, the velocity field is extracted for each time step at the x = 5D plane downstream of the inlet and used to prescribe the velocity field at the inflow plane. At the outlet, a convective boundary condition is used for the velocity to maintain the overall mass conservation, while the pressure is set to a constant value. Figure 9a shows predicted mean wall temperature and Nusselt number as a function of axial distance, where heating starts at a axial position of x∕D = 5.
As it can be seen in Fig. 9, both, the anisotropic heat flux model as well as the isotropic model with Pr sgs = 0.7 show excellent agreement with the experiment (Shehata and McEligot 1998) and also with the reference DNS (Bae et al. 2006) in case of streamwise distribution of the wall temperature and Nusselt number. Furthermore, LES results are very similar to each other, which suggests that both models are well suited to predict such a strongly heated turbulent air flow in a pipe with variable thermo-physical properties, at least in case of mean wall temperatures and Nusselt number.

Turbulent Inclined Jet Impinging on a Heated Solid Surface
Several canonical mixing and fluid flow situations that occur in internal combustion engines and exhaust gas systems can be also found in impinging jet flows. These complex phenomena include (1) thermo-viscous boundary layers, (2) impinging heating/cooling, (3) wall-jets, (4) recirculation and (5) separation. In order to establish the validity of the novel anisotropic heat flux model under such flow conditions, the second application test case consists of a turbulent square jet impinging on a heated solid surface. The heat and fluid flows within this configuration were investigated numerically using DNS technique (see Ries et al. 2018a, b). A schematic of the impinging jet configuration used in the LES study is provided in Fig. 10.
In accordance with the reference DNS, a turbulent jet of dry air ( T inlet = 290 K, p = 1 atm) leaves a square nozzle ( D = 40 mm) and impinges on a heated flat plate. The heated wall has a constant wall temperature of T w = 330 K, a jet-to-plate distance of H∕D = 1 , and an inclination angle of = 45 • . At the impinged wall, the jet is divided into two opposed wall-jets directed outward along the solid wall and gets heated up.
A block-structured numerical grid with 1,699,375 control volumes is employed in the LES study, that is refined in the near-wall region to ensure a non-dimensional wall distance smaller than one. Regarding the inflow, synthetic turbulent inlet condition is employed at the nozzle exit section. Thereby, realistic turbulence is generated using the digital filter approach proposed by Klein et al. (2003), while the mean velocity profile is taken from the DNS study. Figure 11a shows LES predictions of the distribution of the local Nusselt number along the -axis at x = = 0 in comparison with the DNS data. Thereby, the local Nusselt number is defined as Nu = h t D∕ , where h t denotes the local heat transfer coefficient and the thermal conductivity. Figure 11b depicts the turbulent wall-parallel heat flux as a function at ∕D = − 0.15 . At this location, the Nusselt number is maximal, heat is transported counter to the temperature gradient and heat fluxes appear highly anisotropic (see Ries et al. 2018a). In Fig. 11b solid lines denote resolved heat fluxes and dashed lines represent modeled subgrid-scale heat fluxes. LES results of the proposed anisotropic heat flux model and the isotropic heat flux model with Pr sgs = 0.7 are presented.
Examining Fig. 11, a clear peak in the Nusselt number can be observed at ∕D = −0.15 . Additionally, the higher values perceived around the stagnation point rapidly decrease away from this location. This tendency is well reproduced by both  Fig. 11b, it can be clearly seen that the direction of < U � Θ � > changes close to the wall in the DNS data. This trend is only reproduced correctly by means of the proposed anisotropic heat flux model. Furthermore, predictions of < U � Θ � > using the proposed anisotropic heat flux model compare much better to the reference DNS than heat fluxes obtained by using the standard isotropic model. This confirms that only the proposed heat flux model is able to predict resolved and subgridscale heat fluxes in a physically consistent way for such a complex heat and fluid flow situation, while standard isotropic models fail.

Backward-Facing Step Flow with Heated Walls
The last test case in the evaluation study deals with a backward-facing step flow with a constant heated surface behind a sudden expansion. This generic test case features complex flow situations such as recirculation and flow separation and is therefore an excellent benchmark flow for exhaust gas systems, in particular to mimic heat and fluid flow phenomena within exhaust silencer devices. The backward-facing step flow with heated walls was investigated experimentally by Vogel and Eaton (1985). A representation of the computational domain used in the LES study is shown in Fig. 12.
In the backward-facing step test case, a turbulent stream of dry air ( T = 298 K, Pr = 0.71) enters a wind tunnel, expands suddenly after 2 h and is finally heated up with a constant heat flux of q w = 270W∕m 2 behind the sudden expansion. The channel expansion ratio is 1.25 with a Reynolds number of Re = 28, 000 (based on the freestream velocity and step height, h). A block-structured numerical grid with 2,745,504 control volumes is employed in the LES study, that is refined in the near-wall region to ensure a non-dimensional wall distance smaller than one. Realistic inflow turbulence is generated using the digital filter approach proposed in Klein et al. (2003), while the mean velocity profile equals a boundary layer flow profile with a boundary layer thickness of 99 = 1.07 h. Figure 8 depicts (a) temperature profiles at different axial positions and (b) the distribution of the Stanton number along the axial direction at the heated lower wall behind the expansion. Thereby, the Stanton number is defined as St = q w ∕(U ∞ c p (T − T w )) , where U ∞ is the freestream velocity, c p the specific heat capacity of the fluid, the fluid density and T w the wall temperature.
As it can be observed in Fig. 13, there is excellent agreement between LES predictions and the experiment. Mean temperature profiles are very close to the experimental data and peak values in the computed profiles of St compare qualitatively and quantitatively very well with the experiment. This holds true for the anisotropic heat flux model as well as for the isotropic model with Pr sgs = 0.7 . A significant improvement in using an anisotropic heat flux model cannot be determined.

Conclusion
A novel anisotropic heat flux model for large eddy simulations of complex engineering applications have been proposed and evaluated. The prominent features of the proposed model are that (1) it accounts for variable fluid properties and anisotropic effects in the unresolved temperature scales, (2) no ad-hoc treatments or dynamic procedure are required to obtain the correct nearwall behavior, and (3) the formulation is consistent with the second law of thermodynamics. It is shown in this work that only the proposed anisotropic model from the tested ones is able to predict subgrid-scale heat fluxes in a physically consistent way, while both, the standard and dynamic thermal diffusivity models are unable to reproduce the direction of subgrid-scale heat flux correctly. However, the proposed anisotropic heat flux model has a similar prediction accuracy and computational expense than conventional isotropic models. This was confirmed by comparison with DNS and experimental data from the literature for several test cases that are  Vogel and Eaton (1985) relevant for internal combustion engines and exhaust gas systems, namely, a turbulent heated channel flow, a strongly heated air flow in a vertical pipe, a turbulent inclined jet impinging on a heated solid surface and a backward facing step flow with heated walls.