Predictions of Conjugate Heat Transfer in Turbulent Channel Flow Using Advanced Wall-Modeled Large Eddy Simulation Techniques

In this paper, advanced wall-modeled large eddy simulation (LES) techniques are used to predict conjugate heat transfer processes in turbulent channel flow. Thereby, the thermal energy transfer process involves an interaction of conduction within a solid body and convection from the solid surface by fluid motion. The approaches comprise a two-layer RANS–LES approach (zonal LES), a hybrid RANS–LES representative, the so-called improved delayed detached eddy simulation method (IDDES) and a non-equilibrium wall function model (WFLES), respectively. The results obtained are evaluated in comparison with direct numerical simulation (DNS) data and wall-resolved LES including thermal cases of large Reynolds numbers where DNS data are not available in the literature. It turns out that zonal LES, IDDES and WFLES are able to predict heat and fluid flow statistics along with wall shear stresses and Nusselt numbers accurately and that are physically consistent. Furthermore, it is found that IDDES, WFLES and zonal LES exhibit significantly lower computational costs than wall-resolved LES. Since IDDES and especially zonal LES require considerable extra work to generate numerical grids, this study indicates in particular that WFLES offers a promising near-wall modeling strategy for LES of conjugated heat transfer problems. Finally, an entropy generation analysis using the various models showed that the viscous entropy production is zero inside the solid region, peaks at the solid–fluid interface and decreases rapidly with increasing wall distance within the fluid region. Except inside the solid region, where steep temperature gradients lead to high (thermal) entropy generation rates, a similar behavior is monitored for the entropy generation by heat transfer process.


Introduction
As classified in [1], heat transfer processes can be investigated using conjugate, coupled or adjoint formulations corresponding to problems that contain two or more subdomains with phenomena that are described by different differential equations. In many energy systems, such as internal combustion engines, gas turbines, heat exchangers and many more, like converter monolith channels in after-treatment devices, conjugate heat transfer features a thermal energy transfer process that involves the interaction of conduction within a solid body and convection from the solid surface by fluid motion [2]. Note that catalytic converters, which are usually multiple-channel reactors with a honeycomb structure, are additionally characterized by relatively low pressure drop, enhanced mass transfer and large geometric surface area and thickness of the catalyst film on the substrate wall in which heterogeneous chemical reactions take place [3]. The structural field of a solid wall in terms of structural stresses and deformation as well as chemical reactions is not considered in the present paper. Nevertheless, a realistic prediction of conjugate heat transfer problems is very challenging as it requires a coupling of the conduction in the solid part and the convection in the fluid region.
Focusing on the numerical simulation of turbulent conjugate heat transfer, the large eddy simulation (LES) technique has been proven to be an accurate approach to predict such turbulent thermal processes. This was shown in the literature for both generic flow configurations like heated cavity flows [4] and also for complex engineering applications such as cooling of gas turbine blades [5]. However, it is well established that the computational cost of LES with conjugated heat transfer is very high. This is mainly because of the thin momentum and thermal boundary layers at the solid surface that have to be fully resolved in classical LES which requires very fine spatial resolution, in particular for turbulent flows with high Reynolds and Prandtl numbers. Therefore, in order to overcome this issue, it is common practice in LES to use a near-wall modeling approach to reduce the required computational effort of the simulation. In general, such near-wall modeling strategies can be divided in the context of LES into approaches based on wall functions (WFLES), two-layer RANS-LES (zonal LES) and hybrid RANS-LES methods [6].
Regarding WFLES, the momentum and thermal boundary layers are not explicitly resolved by the numerical grid. Instead, they are bridged with a single cell while suitable assumptions are made about the near-wall velocity and temperature profiles [7]. Thereby, in the case of classical wall functions, a linear variation of the near-wall velocity is assumed very close to the wall and a semi-logarithmic variation away from it (e.g., [8][9][10]). Based on Reynolds analogy assumptions, similar variation is also found for the near-wall temperature profile (e.g., [11][12][13][14]). However, it is well known that such simplified formulations for the near-wall velocity and temperature do not apply to complex boundary layers because they are based on simple equilibrium flow assumptions. Therefore, advanced wall functions for the velocity (e.g., [7,[15][16][17][18]) and also for the temperature [7,18] were proposed in the literature that account for additional non-equilibrium effects like time rate change, convection or pressure gradients. This allows accurate predictions of realistic heat and fluid flow applications with a reasonable computational cost, as was recently shown by the authors in [7].
In the case of the zonal LES approach, a numerical grid with a fine spatial resolution is embedded between the matching location of the outer mesh and the solid surface within the fluid region. A simplified set of RANS-based turbulent boundary-layer equations are solved at the embedded mesh region. By means of this, the required wall shear stress is calculated and employed as a wall boundary condition for the LES calculation on the overlapped outer mesh [6,19]. Thereby, pressure gradient and convection effects are taken into account from the solution of the RANS-based turbulent boundary-layer equations. Consequently, two-layer models are also able to capture, to some extent, nonequilibrium heat and fluid flow effects. Nevertheless, it was observed in many numerical studies [20][21][22] that two-layer models tend to overpredict the wall shear stress. In addition to this, the generation of two separate numerical grids can be very challenging, in particular for complex geometries, which impedes the use of zonal LES for practical engineering applications with conjugate heat transfer.
In hybrid RANS-LES modeling approaches, a RANS model is applied in the vicinity of the solid surface, while LES equations with a subgrid-scale model are solved away from it. In this framework, different strategies can be used for the transition from a RANS behavior to a LES behavior, based on criteria updated during the computation [23]. For instance, the turbulent length scale can be changed from a RANS mixing length scale to a grid size-related length scale, or a blending function can be used to merge the RANS and subgrid-scale eddy viscosities [6]. In contrast to wall-resolved LES of conjugate heat transfer, where the grid has to be refined isotropically in all three directions in the vicinity of the solid surface in the fluid region, hybrid RANS-LES requires only grid refinement in the wall-normal direction, leading to a significant reduction in the computational cost [24]. Prominent examples of hybrid RANS-LES models are detached eddy simulations (DESs) [25], delayed detached eddy simulations (DDESs) [26], improved delayed detached eddy simulations (IDDESs) [27], very large eddy simulations (VLESs) [28] or scale-adaptive simulations (SASs) [29].
From this short literature review, it appears that numerous LES near-wall modeling approaches have been proposed in the literature. However, it is worth mentioning that an assessment of the prediction accuracy and computational cost of these wall models regarding turbulent flow with conjugate heat transfer are rarely reported. This motivates the present work that reports on comparative predictions of conjugate heat transfer achieved by means of advanced near-wall modeling approaches in the context of LES for turbulent flows with conjugate heat transfer. For this purpose, wall-modeled LES of turbulent channel flow with conjugate heat transfer was conducted and near-wall statistics were compared with DNS and wall-resolved LES results [30]. A two-layer RANS-LES approach based on the Spalart-Allmaras model [31], an improved delayed detached eddy simulation method (IDDES) [27] and a non-equilibrium wall function model (WFLES) [7] were assessed in terms of prediction accuracy, physical consistency and computational cost. To the authors' knowledge, this is the first study in the literature that presents (i) a comprehensive comparison study of zonal LES, IDDES and WFLES approaches in terms of conjugated heat transfer, (ii) a non-equilibrium WFLES approach applied to simulate conjugate heat transfer and (iii) predictions of entropy production rates based on wall-modeled LES in turbulent conjugate heated channel flow. Furthermore, wall-resolved LES data of conjugate heated channel flow at Re τ = 640 are provided for evaluation purposes that to date have not been available in the literature. This paper is organized as follows. The different near-wall modeling approaches are introduced in Section 2. Subsequently, the turbulent channel flow configuration with conjugate heat transfer is described in Section 3 and the numerical procedure employed for the simulation is outlined. Then, in Section 4, results of the wall-modeled LES are analyzed and the approaches are subsequently assessed. Afterward, an entropy generation analysis using the various models employed is reported. Finally, some concluding remarks are summarized in Section 5.

Wall-Modeled LES Approaches with Conjugate Heat Transfer
Since the structural field of the solid wall in terms of structural stresses and deformation and chemical reactions are not considered, the conjugate heat transfer problem under investigation in the present paper can be divided into three main regions, namely a nonisothermal fluid flow region, a transient heat conduction through the solid and a thermal solid/fluid interface. In the case of turbulent incompressible Navier-Stokes-Fourier fluid flow with conjugate heat transfer and constant physical properties, the transport equations for mass, momentum and energy with respect to RANS and LES can be formulated for the fluid region as: where the concept of eddy viscosity is applied in order to close the LES and RANS equations. In Equations (1)-(3), U i denotes the flow velocity, T f the fluid temperature, p the kinematic pressure, ν the kinematic viscosity, α the molecular thermal diffusivity and f T f , f U i are additional source terms. Note that, in the present study, purely forced convection is investigated. Therefore, no additional source terms for buoyancy effects appear in the balance equations and the temperature is treated as a passive scalar.
Regarding LES, the operator (•) in Equations (1)-(3) represents spatially filtered quantities, ν t is the subgrid-scale viscosity and Pr t the subgrid-scale Prandtl number. Thereby, the wall-adapting local eddy viscosity (WALE) model [7,19,32] is employed in this study to calculate ν t and the subgrid-scale Prandtl number is set to Pr t = 0.5 in accordance with [30,33]. In the context of RANS, the operator (•) denotes time-averaged quantities, ν t is the turbulent eddy viscosity and Pr t the turbulent Prandtl number. In this work, the Spalart-Allmaras turbulence model [31] is used to close the RANS equations and the turbulent Prandtl number is set to Pr t = 1. This turbulent Prandtl number is selected in accordance with the DNS study of turbulent heated channel flow of Kawamura et al. [34].
In the solid region, the velocity is zero in all the balance Equations (1)-(3) and only the energy equation has to be solved, which further simplifies to the classical heat equation: Here, T s represents the solid temperature and α s the thermal diffusivity of the solid. Finally, the solid and fluid regions are coupled via a thermal interface. Here, the temperature and the heat flux of both phases have to be equal, which leads to the following boundary conditions at the fluid-solid interface: where ρ f is the fluid density, c f p the specific heat capacity of the fluid, λ s the thermal conductivity of the solid, α t the turbulent thermal diffusivity and n represents the direction normal to the solid surface.
Next, the different LES near-wall modeling approaches applied in this study to model the momentum and thermal boundary layers in the vicinity of the solid surface are described.

LES with Non-Equilibrium Wall Functions (WFLES)
The basic idea of wall function models is to bridge the momentum and/or thermal boundary layers with a single grid cell and make suitable assumptions about the near-wall velocity and temperature profiles in order to obtain the required wall shear stress and wall heat flux, respectively. In this work, the non-equilibrium wall function approach as proposed by the authors in [7] is employed. In contrast to classical formulations, the nonequilibrium wall functions are continuously valid over the whole range of dimensionless wall distance y + and include transient as well as local non-equilibrium effects like time rate change, adverse pressure gradients, convection and additional source terms. In the following section, only the formulation of the non-equilibrium wall function for the momentum boundary layer is described. A similar procedure applies for the thermal boundary layer and is therefore not shown here for the sake of clarity. A detailed description as well as computational details on the non-equilibrium wall function approach can be found in [7].
Based on the momentum boundary layer equation, an analytical solution for the near-wall velocity profile can be formulated as [7]: with where κ ≈ 0.41 and C 0 = 9.6 · 10 −4 . Here, U + = U/u τ and y + = yu τ /ν are the dimensionless mean velocity and wall distance, respectively, and u τ is the friction velocity. In this formulation, all local non-equilibrium effects are combined into C + U = νC u /ρ f u 3 τ , where C U is calculated at the cell centroid P of the first cell at the wall as: with ζ as the flow direction and η the wall-normal direction. Based on Equation (6), the required wall shear stress can be determined using an iterative procedure (e.g., Newton-Raphson or regula falsi methods) and applied as a boundary condition for the LES with the wall function approach. A similar procedure was used in this work to bridge the thermal boundary layer. Thereby, the thermal diffusivity α t was calculated in an iterative procedure from the temperature wall function and employed as a boundary condition.

Two-Layer RANS-LES Approach (Zonal LES)
As mentioned above, in the two-layer RANS-LES approach (zonal LES), a numerical grid with a fine spatial resolution, is embedded between the matching location of the outer mesh and the solid surface within the fluid region. Thereby, in the present study, the Spalart-Allmaras [31] model was applied to solve the RANS equations at the inner layer on the embedded mesh. The Spalart-Allmaras RANS model reads [31]: where the turbulent viscosity is calculated as The coefficients in the Spalart-Allmaras RANS model are defined as where Ω ij = 1/2 ∂U i /∂x j − ∂U j /∂x i is the rotation tensor and d a characteristic length scale defined as the distance to the wall. The model constants are given as Note that in contrast to the original Spalart-Allmaras RANS model, the trip term in Equation (11) is not considered in the present model formulation.
At the outer layer, the wall-adapting local eddy viscosity model (WALE) [32] is employed to solve the LES equations. In the WALE model, the subgrid-scale viscosity is expressed as [32] where C W = 0.325 is the model coefficient, ∆ = (∆ x ∆ y ∆ z ) 1/3 the grid filter and S d ij is the traceless symmetric part of the square of the velocity gradient tensor.
In the zonal LES approach, the resolved velocity and the pressure gradient from the LES calculation serve as boundary conditions for the inner-layer RANS simulation. Thereby, the wall stress from the RANS is returned as wall boundary condition for the LES calculation. Regarding the treatment of the thermal boundary layer, the turbulent thermal diffusivity α t is calculated in the zonal LES approach based on Reynolds analogy assumptions as α t = ν t /Pr t .

Improved Delayed Detached Eddy Simulation (IDDES)
Similar to [27], the Spalart-Allmaras eddy viscosity transport equation (see Equation (11)) is used in the present IDDES approach in order to achieve an eddy viscosity. Thereby, in contrast to the classical Spalart-Allmaras RANS model, a hybrid turbulent lengthscale formulation, that blends between a RANS and a LES length scale, is used for the approximation of d as [35]: Here, l RANS is a RANS-based turbulent length scale and l DES a length scale that depends on the grid width. The RANS-based turbulent length scale equals the distance to the wall l RANS = d. The grid-based length scale is calculated as l DES = ΨC DES ∆, where Ψ is the low Reynolds number correction function (see [27]), ∆ = (∆ x ∆ y ∆ z ) 1/3 the grid filter and C DES = 0.65 a model constant. The blending function f d in Equation (15) is defined in such a way that l RANS is predominantly used in regions with low mesh resolution and l DES in regions where the grid resolution is sufficient for LES. The elevation function f e aims at preventing an excessive reduction in the Reynolds stresses in the vicinity of the RANS-LES interface [35]. A detailed description of the IDDES model and the blending functions can be found in [27,35]. Similar to the zonal LES approach, the turbulent thermal diffusivity α t is calculated in IDDES based on Reynolds analogy assumptions as α t = ν t /Pr t .

Configuration and Numerical Procedure
In line with the numerical study of Flageul et al. [30], a turbulent heated channel flow test case with conjugate heat transfer was selected in this work. The heated channel flow was simulated for a fluid with a molecular Prandtl number of Pr = ν/α f = 0.71 and at Re τ = 395, 640, 1020, where Re τ = u τ δ/ν is the Reynolds number based on the friction velocity. A fluid-to-solid thermal diffusivity ratio of G 1 = α f /α s = 1 and a solid-to-fluid thermal conductivity ratio of G 2 = λ s /λ f = 1 were selected, leading to a thermal activity ratio of K = 1/G 2 √ G 1 = 1 [30]. These values of G 1 = G 2 = K = 1 were selected in accordance with the reference DNS of [30] and represent the case of a coupled scalar with the same thermal properties in the fluid and solid region. A sketch of the computational domain used for the simulations is shown in Figure 1, where x, y and z are the spanwise, wall-normal and streamwise directions, respectively. The entire computational domain has a length of 6.4δ, a width of 3.2δ and a height of 4δ, where δ is the channel half-height. Thereby, similar to [30], the fluid domain is bounded at −δ < y < δ and the solid domains are located at y > δ and y < −δ, respectively. Both solid domains have an height of δ, which ensures that the boundary condition used at the outer wall has no significant impact on the statistics at the fluid-solid interface [30].
Periodic boundary conditions were applied for the velocity and temperature in streamwise and spanwise directions. At the solid surface, a no-slip condition was employed for the velocity and a coupled thermal boundary condition was used for the temperature (see Equation (5)). The pressure and temperature gradients that drive the heat and fluid flow in the fluid region are adjusted dynamically to maintain a constant mass flux and mean mixed temperature. Therefore, the pressure and temperature were split into a periodic and a non-periodic part. Source terms for the non-periodic part, f U x and f T f , were added to the momentum and temperature equation, respectively (see [34]). Synthetic turbulence was used to initialize the channel flow simulations in the fluid region. Thereby, to avoid uncertainties caused by the initial solution, the start-up phase of the simulations was chosen to be long enough to ensure that the channel flow was fully developed before sampling started. A detailed description of the initialization method can be found in [36].
Three-dimensional block-structure numerical grids with different spatial resolutions were used for each Re τ and also for each LES near-wall modeling approach. To complete the DNS dataset, additional wall-resolved LES (WRLES) were conducted for each Re τ , including the thermal cases for large Re-numbers for which DNS data were not available (Re τ = 640, 1020). In the case of IDDES and WRLES, the numerical grid was refined towards the wall in order to ensure a non-dimensional wall distance y + w smaller than one. A representation of each grid used for the WFLES, WRLES, IDDES and zonal LES approaches is shown in Figure 2. Thereby, the coarsest spatial resolutions for each of the test cases at Re τ = 395 are shown (see Table 1   The balance equations for turbulent flow with conjugate heat transfer were solved numerically using an incompressible version of chtMultiRegionFoam from the open-source software OpenFOAM v1912 [37]. Thereby, the temperature transport equation, the LES near-wall modeling approaches and the source terms that drive the channel flow were added to the source code. Regarding the fluid region, a merged PISO-SIMPLE ( [38,39]) algorithm was employed for the pressure-velocity coupling. The solution procedure was applied with a second-order implicit backward-differencing scheme for the time integration, a low-dissipative second-order flux-limiting differencing scheme for the convection terms and a conservative scheme for the Laplacian and gradient terms. The time step of the simulations was chosen to be small enough to ensure that the Courant-Friedrichs-Lewy number remained smaller than one. Convergence optimization and acceleration techniques were incorporated in order to speed up the calculations. In particular, a geometric agglomerated algebraic multigrid solver was considered for the resolution of the pressure Poisson equation and for the momentum predictor. Convergence of the iterative procedure is obtained if the normalized residuals of all governing equations are reduced by more than three orders of magnitude within each time step.
Several wall-resolved and wall-modeled LESs have been carried out in the present work. A summary of all the test cases is given in Table 1.

Results
The comparative predictions of the conjugate heat transfer process due to conduction and convection as achieved by the different near-wall treatments for LES are divided into five parts. First, predictions of instantaneous temperature and velocity fields are presented and compared for the different modeling approaches. Then, fluid flow statistics from the different wall-modeled LES approaches are reported together with their effect on the heat transfer. The obtained predictions are compared with DNS and wall-resolved LES in order to analyze the accuracy of the different near-wall treatments. The comparison includes mean and rms values in the fluid and solid region as well as friction coefficients and Nusselt numbers. Next, results of a systematic grid-dependency study are presented that allow us to assess the influence of the spatial resolution on predictions obtained by the different near-wall treatments. Subsequently, the physical consistency of the modeling approach is investigated. In particular, predicted entropy generation rates obtained by the different wall-modeled LES approaches are compared with results of wall-resolved LES. Finally, the computational cost of the wall-modeled LES approaches is quantified and compared to classical wall-resolved LES methods in order to highlight the benefit of each near-wall treatment for practical LES.  Table 1). Note that the temperature color scale is subdivided into a range for the fluid region and the solid region in order to better visualize the wide range of temperature scales.   Table 1).

Instantaneous Temperature and Velocity Fields
It can be clearly seen in Figure 3 that the velocity and temperature fields are highly turbulent in the fluid region with steep velocity/temperature gradients close to the wall. Thereby, due to the finer grid resolution, more of the small-scale turbulent structures are resolved in the WRLES than in case of IDDES, zonal LES and WFLES. In contrast, the temperature field in the solid region is homogeneous distributed with a steep gradient in the wall-normal direction. Thereby, it appears that the predicted temperature fields obtained by the different wall-modeled approaches (IDDES, zonal LES, WFLES) are quite similar and compare well with the WRLES. Figure 4 presents non-dimensional mean velocity U + and turbulent kinetic energy k + profiles as a function of dimensionless wall distance y + for the turbulent heated channel flow configuration at Re τ = 395, 640, 1020. Predictions of WFLES, IDDES and zonal LES are compared with DNS data from Abe et al. [40] as well as results of wall-resolved LES that have been computed in the present study (cases 1, 11, 21 of Table 1). Note that results of the wall-modeled LES are only shown for the medium spatial resolution (cases 3, 6, 9, 13, 16, 19, 23, 26, 29 of Table 1). Similar trends are found for the other spatial resolutions and are therefore not shown for the sake of clarity.  In the case of WRLES, it is apparent in Figure 4 that the agreement with DNS is very satisfactory for both U + and k + . This confirmed the validity of the present WRLES results and allowed us to use this dataset as a reference for further assessment of the LES near-wall modeling strategies in this study. Regarding the predictions of the wall-modeled LES, it appears that mean velocity profiles agree very well with the DNS and WRLES data. This holds true for all near-wall treatments and also for all Re τ under consideration. Less good agreement can be observed for the predicted turbulent kinetic energy k + . In particular, the peak value of k + is underestimated in the case of WFLES and zonal LES, while it is predominantly overestimated in the case of IDDES. Nevertheless, the overall agreement of the wall-modeled LES approaches is still satisfactory, which confirms the physical consistency of all LES near-wall modeling approaches under consideration in terms of turbulent channel flow. Figure 5 shows predicted non-dimensional mean temperature Θ + and rms temperature Θ + rms profiles in the fluid region based on WFLES, IDDES and zonal LES treatment. Thereby, the non-dimensional temperature is defined as Θ + = (T w − T)/T τ , where T τ = q w / ρ f c f p u τ is the friction temperature with q w the wall heat flux which is defined in the case of wall-modeled LES as q w = α f + α t [∂T/∂y] y=0 . For comparison purposes, DNS data from [30] of a turbulent channel flow with conjugate heat transfer at Re τ = 395 were employed. For higher Re τ , the results of wall-resolved LES (cases 1, 11, 21 of Table 1) were used for comparison because thermal statistics from DNS of this specific configuration at higher Re τ are not available in the literature. Note that in Figure 5, the value of Θ rms is not zero at the wall. As indicated in [30], the temperature variance at the wall depends strongly on the value of the thermal activity ratio K. Thereby, lower values of K correspond to conjugate cases similar to isothermal boundary conditions, while higher values of K correspond to isoflux conditions at the thermal interface. . Dimensionless mean and rms temperature Θ + , Θ + rms as a function of non-dimensional wall distance y + for Re τ = 395, 640, 1020. Comparison of wall-modeled with wall-resolved LES (cases 1, 3,6,9,11,13,16,19,21,23,26,29). : WFLES; : IDDES; : zonal LES; : WRLES; : DNS.

Fluid Flow Statistics and Impact on Heat Transfer
According to the fluid flow statistics in Figure 4, predictions of mean and rms temperature profiles based on WFLES, IDDES and zonal LES compare well with the WRLES reference data. Significant differences in the predictions from the individual wall-modeled LES approaches cannot be observed. Therefore, it can be concluded that WFLES, IDDES and zonal LES are able to reproduce thermal statistics in the fluid region properly in the case of turbulent channel flow with conjugate heat transfer. This holds true for all Re τ under consideration.
Next, predictions of skin friction coefficients C f = 2τ w /(ρU 2 b ) and Nusselt numbers Nu = ((α + α t )/α) · δ[∂T/∂y] y=δ /(T w − T δ ) based on WFLES, IDDES and zonal LES are compared in Figure 6 with respect to WRLES for different Re τ . Here, τ w is the wall shear stress, U b the bulk velocity, T w the wall temperature and T δ the temperature at y = δ. The integral quantities C f and Nu are particularly relevant for the development and design of industrial/engineering applications. Accurate predictions of such quantities are therefore important in the context of practical LES with near-wall modeling. Regarding C f , it can be clearly seen in Figure 6 that WFLES, IDDES and zonal LES are able to reproduce the skin friction coefficient accurately for all Re τ under consideration. In contrast, discrepancies in the prediction of Nusselt numbers are more significant. Thereby, in particular, the zonal LES approach slightly overestimates values of Nu, while results of WFLES and IDDES are still very close to the reference WRLES. Nevertheless, predictions of Nu obtained by zonal LES are still acceptable for wall-modeled LES, which leads to the conclusion that all wall-modeled LES approaches are suitable to predict heat and fluid flow statistics including skin friction coefficients and Nusselt number at the fluid region of turbulent channel flow with conjugate heat transfer.
Finally, mean temperature Θ + and rms temperature Θ + rms profiles in the solid region are presented in Figure 7. Predictions of WFLES, IDDES and zonal LES are compared with WRLES. Results are solely shown for the turbulent channel flow at Re τ = 1020. Similar trends are found for Re τ = 395 and 640 and are therefore not shown for the sake of clarity. As can be clearly seen in Figure 7, predictions of Θ + obtained by the individual wall-modeled LES approaches agree very well with the reference wall-resolved LES. This holds true for the interface as well as for the rest of the solid region. Similarly, profiles of Θ + rms are reproduced well by the wall-modeled LES, except in the case of zonal LES. Here, values of Θ + rms decrease too rapidly towards the outer wall. This unphysical behavior of the thermal fluctuation penetration may be caused by overestimated values of Nu at the thermal interface (see Figure 6), resulting in too intensive heat transfer in the context of zonal LES.
By examining predictions of heat and fluid flow statistics within the turbulent channel flow configuration with conjugate heat transfer, it turned out that WFLES, IDDES and zonal LES are able to reproduce the physics of such heated flows properly. The influence of the spatial resolution on the prediction accuracy is analyzed next. Figure 8 shows predictions of mean velocity U + and rms velocity U + rms by IDDES, zonal LES and WFLES. Three different spatial resolutions are shown, denoted here as coarse, medium and fine. For comparison, DNS data of Abe et al. [40] and results of WRLES are depicted.  Table 1). Comparison with DNS data of [40] and wall-resolved LES (cases 1, 11, 21 of Table 1). : coarse grid; : medium grid; : fine grid; : WRLES; : DNS.

Grid Dependency
From Figure 8, it appears that predictions of U + are very close to the reference data and nearly independent of the spatial resolution. This holds true for IDDES, zonal LES and WFLES over the entire range of y + . In contrast, predicted profiles of k + are more affected by the spatial resolution. In particular, the zonal LES and WFLES approaches predict a non-physical peak of k + close to the wall. Such a non-physical peak of k + might be caused by grid resolution requirements [41], the numerical schemes employed [42] and effects of the subgrid-scale model, among other numerical or modeling errors [43]. However, it can be seen that results are still reasonably close to the reference WRLES and DNS. Obviously, the grid dependency of mean and rms velocities obtained by wall-modeled LES is not very significant, at least for turbulent channel flow at Re τ = 1020.
Next, the grid dependency of mean temperature Θ + and rms temperature Θ + rms is analyzed in Figure 9. Regarding mean temperature profiles, it is visible in Figure 9 that predictions based on IDDES compare very well with the reference DNS. This holds true for all spatial resolutions that are considered in the present study. In contrast, WFLES and zonal LES are more affected by the spatial resolution. Thereby, values of Θ + are slightly underestimated, in particular for the fine grid resolution. As indicated in [41], such a log-layer mismatch might be caused by numerics, grid resolution requirements, effects of the subgrid model, etc. Regarding Θ + rms , the influence of the spatial resolution is less significant. However, it can be observed that the peak value of Θ + rms is slightly shifted in the case of the IDDES approach and coarse grid. Nevertheless, similar to the fluid flow statistics, the grid dependency of mean and rms temperatures in the fluid region obtained by wall-modeled LES seems to be not very significant in the case of turbulent channel flow with conjugate heat transfer at Re τ = 1020.
Finally, the grid dependency of thermal statistics in the solid region were studied. For this purpose, Figure 10 presents mean temperature Θ + and rms temperature Θ + rms profiles in the the vicinity of the wall. Results are shown for IDDES, zonal LES and WFLES for different spatial resolutions. For comparison, WRLES data were employed. Just as is the case for the fluid region, predictions of mean temperature profiles are nearly independent of the spatial resolution in the solid region. Similarly, values of Θ + rms are more or less independent of the spatial resolution, except in the case of zonal LES. Here, Θ + rms decreases too rapidly towards the outer wall. This is particularly visible for the coarse and medium grid resolutions.
In summary, the grid dependency of heat and fluid flow statistics obtained by wallmodeled LES is not very significant. This holds true for IDDES, WFLES and zonal LES in the case of turbulent channel flow with conjugate heat transfer at Re τ = 1020.

Physical Consistency of the Modeling
In addition to precise predictions of thermal and fluid flow statistics, it is also important that a wall-modeled LES approach is consistent with the second law of thermodynamics to ensure that thermodynamic processes, as they occur in heat and flow systems, are correctly described [7]. In the case of LES of turbulent fluid flows with conjugate heat transfer, the second law of thermodynamics can be written in the form of the filtered local entropy inequality at the continuum level [44][45][46][47] as: where (•) denotes spatial filtering. The terms on the left-hand side are the local change, convection and flux of entropy density s [7]. The source terms on the right-hand side represent the local entropy production rates by viscous dissipation Π v and by heat transfer Π q .
In the context of Navier-Stokes-Fourier fluid and LES with conjugate heat transfer, the mean of the filtered source terms, Π v and Π q , can be formulated as [46,47]: where • denotes temporal averaging, (•) res is the resolved contribution and (•) sgs represents the subgrid-scale contribution. In this work, C OC = 1.34 is the Obukhov-Corrsin constant [48] and C S the Smagorinsky coefficient [49]. The latter can be directly related to the WALE model as C S = C W / √ 11.27 [32]. Figure 11 shows predicted entropy production rates (a) by viscous dissipation Π v and (b) by heat transfer Π q in the vicinity of the wall. Thereby, Π v is normalized as where T w is the temperature at the fluid-solid interface. Results are shown for the turbulent channel flow configuration at Re τ = 1020 and medium grid resolution (cases 21,23,26,29). Note that Π + v is zero in the solid region due to the absence of velocity gradients. It is visible in Figure 11 that Π + v is zero inside the solid region, peaks at the solid-fluid interface and decreases rapidly with increasing wall distance within the fluid region. Obviously, irreversibilities in such flows occur predominantly at the solid-fluid interface where velocity gradients are high. A similar conclusion can be drawn for entropy generation by heat transfer Π + q , except inside the solid region, where steep temperature gradients lead to high entropy generation rates. Both characteristic trends are well captured by the different wall-modeled LES approaches, which confirms their physical consistency in terms of turbulent channel flow with conjugate heat transfer. It turns out that the behavior predictions close to the wall using zonal LES or WFLES are more affected by the dimensionless wall distance of the first cell at the solid surface than when using IDDES.

Computational Cost
One of the key objectives of using near-wall modeling in the context of LES is to reduce the computational effort in order to allow the calculation of high Reynolds number flows [19]. Therefore, it is of practical interest to address the required computational cost to achieve an acceptable prediction accuracy. For this purpose, Figure 12 shows the required relative computational cost of the different wall-modeled LES approaches with respect to wall-resolved LES. Thereby, the relative computational cost of a wall-modeled LES approach is defined in this work as the ratio of the CPU time spent for the calculation of a wall-modeled LES and the CPU time that is required for a wall-resolved LES of the same configuration. The computational cost was estimated in the present study on a Linux 3.10.0-514.26.1.el7.x86 64 Red Hat 4.8.5-11 (x86 64) system using an Intel(R) Core(TM) i5-6600K CPU @ 3.50GHz and 32GiB RAM. Only one CPU core was used to quantify the computational cost and the maximal Courant-Friedrichs-Lewy number of the simulations was set to CFL = 0.3. The use of only one CPU core for estimating the runtime performance was chosen because the parallel scalability is not in the scope of this study. Note that the relative computational cost depends generally not only on the near-wall treatment model, but also on the selected test case, the particular code implementation and the solution procedure applied. Regarding the turbulent channel flow test case with conjugate heat transfer, it was found in this work that the computational cost of WFLES and zonal LES is about 50-100 times lower than in the case of WRLES. The computational cost of IDDES is significantly higher than that of WFLES and zonal LES but still 3-10 lower compared to WRLES. Furthermore, it can be observed that the relative computational cost of all tested wallmodeled LES approaches decreases with increasing Reynolds number. However, due to the considerable extra work to generate numerical grids in the case of IDDES and especially for zonal LES, this work suggests that WFLES offers, in particular, a promising near-wall modeling strategy for LES of conjugate heat transfer in realistic engineering applications.

Conclusions
The impact of different wall-modeled LES approaches, namely, LES with non-equilibrium wall functions (WFLES), two-layer RANS-LES (zonal LES) and improved delayed detached eddy simulation (IDDES), on the heat transfer predictions has been investigated in a turbulent channel flow with conjugate heat transfer at Re τ = 395, 640, 1020. Thereby, the heat transfer process involved the interaction of conduction within the channel solid wall and convection from the solid wall surface by the fluid flow. In addition, the physical consistency, accuracy and computational cost of the different near-wall modeling strategies have been assessed.
Using available DNS data as a reference, the fluid flow statistics obtained were first Regarding the turbulent channel flow test case with conjugate heat transfer, it was found in this work that the computational cost of WFLES and zonal LES is about 50-100 times lower than in the case of WRLES. The computational cost of IDDES is significantly higher than that of WFLES and zonal LES but still 3-10 lower compared to WRLES. Furthermore, it can be observed that the relative computational cost of all tested wallmodeled LES approaches decreases with increasing Reynolds number. However, due to the considerable extra work to generate numerical grids in the case of IDDES and especially for zonal LES, this work suggests that WFLES offers, in particular, a promising near-wall modeling strategy for LES of conjugate heat transfer in realistic engineering applications.

Conclusions
The impact of different wall-modeled LES approaches, namely, LES with non-equilibrium wall functions (WFLES), two-layer RANS-LES (zonal LES) and improved delayed detached eddy simulation (IDDES), on the heat transfer predictions has been investigated in a turbulent channel flow with conjugate heat transfer at Re τ = 395, 640, 1020. Thereby, the heat transfer process involved the interaction of conduction within the channel solid wall and convection from the solid wall surface by the fluid flow. In addition, the physical consistency, accuracy and computational cost of the different near-wall modeling strategies have been assessed.
Using available DNS data as a reference, the fluid flow statistics obtained were first evaluated and compared to wall-resolved LES (WRLES). Then, the effect of these predictions on the heat transfer was assessed. The following important findings are worth mentioning: (i) WFLES, IDDES and zonal LES are able to reproduce the physics of turbulent channel flow with conjugate heat transfer properly. (ii) The grid dependency of fluid flow statistics obtained by IDDES, WFLES and zonal LES is not very significant. The same is valid for the thermal field. (iii) The computational cost of IDDES, WFLES and zonal LES of turbulent channel flow with conjugate heat transfer is considerably lower than in the case of WRLES. In particular, WFLES and zonal LES allow accurate predictions with a reasonable computational cost. (iv) The relative computational cost of the wall-modeled LES decreases with increasing Reynolds number.
Finally, it was found that IDDES and especially zonal LES require considerable extra work to generate numerical grids. Therefore, in the authors' opinion, WFLES also offers an especially promising near-wall modeling approach for LES of conjugate heat transfer in realistic engineering applications. An entropy generation analysis using the various models was carried out. From this analysis, it turned out that the viscous entropy production is zero inside the solid region due to the prevailing zero velocity gradient, peaks at the solid-fluid interface and decreases rapidly with increasing wall distance within the fluid region. Except inside the solid region, where steep temperature gradients lead to high (thermal) entropy generation rates, a similar conclusion can be drawn for the entropy generation by the heat transfer process.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

LES
Large eddy simulation DNS Direct numerical simulation RANS Reynolds-averaged Navier-Stokes CFD Computational fluid dynamics WRLES Wall-resolved LES WFLES LES with wall functions Zonal LES Two-layer RANS-LES IDDES Improved delayed detached eddy simulation