Numerical investigation of boundary layer flow and wall heat transfer in a gasoline direct-injection engine

International Journal of Heat and Mass Transfer 120 (2018) 1189–1199

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Numerical investigation of boundary layer flow and wall heat transfer in a gasoline direct-injection engine Xueqi Fan, Zhizhao Che, Tianyou Wang ⇑, Zhen Lu State Key Laboratory of Engines, Tianjin University, Tianjin 300072, China

a r t i c l e

i n f o

Article history: Received 19 May 2017 Received in revised form 16 August 2017 Accepted 23 September 2017

Keywords: Turbulent boundary layer Heat transfer DDES GDI engine

a b s t r a c t Near-wall turbulent boundary-layer has a substantial impact on the wall heat transfer process in internal combustion engines, and ultimately affects the engine performance. However, the heat transfer processes are not well understood because of a lack of information on gas side velocity and temperature distribution. In the present study, the transient velocity and thermal boundary layers in a motored spark-ignition direct-injection engine are predicted using a Delay Detached Eddy Simulation Shear-stress Transport (DDES-SST) model. The near-wall ensemble-averaged and fluctuation velocity and temperature fields at the cylinder head are investigated throughout the compression stroke. The numerical results of the flow field in engine core region agree well with high-speed Two Dimensional Three ComponentsParticle Image Velocimetry (2D3C-PIV) experimental data. The simulated data show that boundarylayer separation occurs around the closing phase of intake valves and the turbulence in the near-wall region is anisotropic. The thickness of thermal boundary layer shows a non-monotonic variation with time. The dimensionless velocity profiles agree well with the law of the wall in the viscous sublayer but deviate from the log layer. The dimensionless temperature profiles tend to be uniform in the loglayer. The maximum local wall heat flux and heat transfer coefficient are 70 kW/m2 and 200 W/(m2 K) at top dead center (TDC), respectively. The spatial distribution of turbulent heat flux shows a strong link between the velocity and thermal boundary layers, and the peaks in the buffer layer occupy about 70– 85% of the local wall heat flux, indicating that the forced convection heat transfer is the main heat transfer mechanism in internal combustion engines. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction In-cylinder turbulent flow has a substantial impact on air-fuel mixture formation and flame propagation in internal combustion engines (ICEs), and the turbulent boundary layer contributes a lot to the wall heat transfer which ultimately affects the engine efficiency, pollutant formation and component thermal stresses [1– 3]. Therefore, quantitative prediction of local and transient heat transfer is important not only to accurate simulations but also to the optimization of ICEs [4]. Due to experimental limitations in measuring the velocity and temperature distribution in the nearwall region, the heat transfer processes are still not well understood [5]. Besides, the boundary layer flow in ICEs is influenced by the non-equilibrium turbulence due to high-speed reciprocating motion of the piston and the periodic opening and closing of the valves, so the current wall function models are not suitable for the complex conditions in engines [6]. Hence, understanding the ⇑ Corresponding author. E-mail address: [email protected] (T. Wang). https://doi.org/10.1016/j.ijheatmasstransfer.2017.09.089 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

details of ICE boundary layer flows and wall heat transfer is crucial for the design and optimization of ICEs. Despite the importance of understanding wall heat transfer in ICEs, only a few studies have been focused on the near-wall turbulence flow compared with the flow structures in engine core region. Laser Doppler Velocimetry (LDV) was firstly applied to investigate the near-wall velocity filed in ICEs in 1980s [7]. Foster & Witze [8] measured as close as 60 mm from the cylinder head and found that the boundary layer thickness was less than 1 mm. Pierce et al. [9] conducted LDV with Particle Image Velocimetry (PIV) to measure the 3D velocity field on the cylinder head and proved that the boundary layer in ICEs was at least a 2D velocity field. Recently, high-speed micro-PIV and Particle Tracking Velocimetry (PTV) were applied in the near-wall region in ICEs, and greatly improved the measurement accuracy. Alharbi & Sick [10–13] measured as close as 45 mm from the cylinder head on a motored engine and observed reversal flows and millimeter vortex structures. Besides, they found that the dimensionless velocity profiles agreed well with the law of the wall in the viscous sublayer but deviated from the log layer. Jainski et al. [14] extended the

1190

X. Fan et al. / International Journal of Heat and Mass Transfer 120 (2018) 1189–1199

above measurement to three different engine speeds and found that the thickness of the viscous sub-layer decreased as the engine speed increased. Koehler et al. [15] conducted time-resolved PIV (TR-PIV) to measure as close as 1.5 mm from the piston crown, but the experiment could only be carried out during the intake stroke because of the working fluid was water rather than air. The thermal boundary layers in ICEs can be accessed by Schlieren photography, Coherent Anti-Stokes Raman Scattering (CARS) and Laser-Induced Fluorescence (LIF) techniques. Lucht et al. [16,17] measured the thermal boundary layer on the cylinder head using CARS, and found that the thermal boundary layer thickness in fired engine increased during the expansion stroke. Cundy and Sick [12,18,19] conducted LIF to measure the temperature distribution in the combustion chamber, and the results showed that the thermal boundary layer thickness was thinner than the results obtained by Lyford-Pike [20] using Schlieren photograph. They also observed thermal stratification on a cooled curved metal plate from the transient temperature fields. Tran et al. [21] observed the non-uniform temperature distribution in a Homogeneous Charge Compression Ignition (HCCI) engine by using LIF. Furthermore, the instant wall heat flux in ICEs was generally measured by thermocouples. Nijeweme et al. [22] found that the wall heat flux was influenced by the test location, engine speed, throttle setting and ignition time. Gingrich et al. [23] compared the instant heat flux on the piston crown under different combustion regimes. Ma et al. [24] measured the instantaneous heat flux on the cylinder head, the peak heat flux was about 60 kW/m2 under motored condition at the engine speed of 500 rpm. Since those experiments could only provide velocity and temperature distribution in limited areas, multi-dimensional simulations are necessary to solve the transient flow and heat transfer in ICEs. However, the wall models often used in ReynoldsAveraged Navier-Stokes (RANS) and Large-Eddy Simulation (LES) in ICEs are too simple and need to be refined [25]. The most widely used wall model is the classic equilibrium wall-function model derived by Prandtl [26]. Han and Reitz [4] improved this model by considering the density and turbulent Prandtl number variation in the near-wall region, and the improved model has also been widely used in multi-dimensional simulations in ICEs. After that, many researchers further improved the Han-Reitz heat transfer model in ICEs by considering the pressure work term [27], the variation of gas viscosity and turbulent Prandtl number [28,29] and the subgrid-scale (SGS) turbulent viscosity [30]. Ma et al. [31] conducted k-omega turbulence model with pressure gradient to solve the 2D unsteady velocity boundary layer under the experimental condition conducted by Jainski and Sick [11,14], the gas temperature was about 640 K at 330 CAD under motored condition at the engine speed of 800 rpm. In addition to the wall-function models, hybrid LES/RANS method and Detached Eddy Simulation (DES) have been increasing popular to simulate in-cylinder flows. Jhavar et al. [32] solved the in-cylinder flow based on Very Large Eddy Simulation (VLES) and RANS k-e model. Hasse et al. [33] conducted Shear-Stress Transport DES (SST DES) to simulate in-cylinder flow and the results agreed better with experimental data than the RANS approach. Hartmann et al. [34] investigated the velocity boundary layer at the intake valve seat based on SST DES model and observed flow separation in the near-wall region. Buhl et al. [35] investigated the velocity boundary layer on the piston crown during the intake stroke using Scale Adaptive Simulation Shear-stress Transport (SAS-SST) model, and the results showed that the boundary layer thickness varied strongly along the piston surface due to the large-scale tumble flow. Direct numerical simulation (DNS) has been widely used to investigate the dynamics of the coherent structures and the physics of turbulent heat transfer within the turbulent boundary layers

of simple geometry like channels, pipes and plates [36]. Shadloo et al. [37,38] investigated the effect of wall heat transfer on turbulent statistics and near wall behaviors in supersonic turbulent boundary layers under different wall conditions using DNS. The results showed that increasing the disturbance amplitude as well as perturbation frequency moved the laminar-to-turbulent transition upstream. Wu & Moin [39] simulated an incompressible, zeropressure-gradient flat-plate boundary layer and observed that hairpin vortices developed into the downstream hairpin forests in the transitional region (800 < Re < 1900). However, DNS method is seldom used in ICEs because of the high CPU demand. Schmitt et al. [40–42] conducted DNS to simulate the boundary layer flow with wall heat transfer in an engine-relevant condition. The results showed that the velocity and thermal boundary layer thicknesses decreased towards top dead center, and the dimensionless velocity and temperature profiles deviated strongly from the law of the wall. The DNS data showed that the turbulent heat flux occupied almost 60% to 80% of total wall heat flux. These conclusions were significant but the simulation geometry was simplified without combustion chamber structure and moving valves, so it could not reflect the real boundary layer flow in ICEs. In summary, the experiments could only access limited areas on wall surfaces in ICEs and there is a lack of simulation of velocity and thermal boundary layers in realistic engines to the authors’ knowledge. The objective of the present study is to analyze the velocity boundary layer and thermal boundary layer together in a realistic engine, and provide a further understanding of wall heat flux and turbulent heat transfer in ICEs. The simulation setup and experimental validation are presented in Section 2. In Section 3, the evolution of velocity and thermal boundary layers during the compression stroke are described, the dimensionless velocity and temperature profiles are compared against the law of the wall, and local wall heat flux and heat transfer coefficient are predicted. Conclusions are drawn in Section 4.

2. Numerical method 2.1. Governing equations and turbulence modeling Since the characteristic velocity in ICEs is on the order of 10 m/s, the in-cylinder flow can be described by the low Mach number compressible Navier-Stokes equations [31]. The compressibility of the fluid is ruled by the ideal gas law. The governing equations for mass, momentum and energy are as follows:

@q þ r  ðquÞ ¼ 0 @t
 @ui @P q þ rðl  gradui Þ þ Si þ rðui  uÞ ¼  @xi @t


qcp

@T þ u  rT @t

 ¼

@P0 þ rðk  rTÞ @t

P0 ¼ qRT

ð1Þ ð2Þ

ð3Þ ð4Þ

where P is the hydrodynamic pressure, P0 is the thermodynamic pressure, and the physical pressure is the sum of this two terms. cp is the specific heat at constant pressure, l is the dynamic viscosity and k is the thermal conductivity. Among them [43],

l ¼ 1:46  106

T 1:5 T þ 111

k ¼ 2:089  103

T 1:5 T þ 111

ð5Þ

ð6Þ

1191

X. Fan et al. / International Journal of Heat and Mass Transfer 120 (2018) 1189–1199

cp ¼ 1:05  0:365


2
3 T T T þ 0:85  0:39 1000 1000 1000

ð7Þ

In the momentum equation, ui are the velocity components, xi are the coordinates, Si are the generalized sources.

Si ¼

@ @xi




l



 @ui @ @u @ @u @ þ l j þ l k þ ðkruÞ @xj @xk @xi @xi @xi @xi

ð8Þ

where k is the second viscosity coefficient. The turbulence model used in this study is DDES-SST which is a seamless hybrid URANS/LES method. A major advantage of this model is that it only requires grid refinement in the wall-normal direction in the near-wall region to solve wall effect. Therefore, it has higher precision than wall function methods and less CPU demand than LES approach. The SST RANS model is a combination of standard k-e and k-x model. It equals the k-x model within the boundary layer to solve attachment and the k-e model in the outer region to solve separation. The turbulent length scale of SST RANS

Fig. 1. Domain of the simulation.

model is lt ¼ k =ðb xÞ. LES approach would be applied if the grid resolution is sufficiently fine to resolve smaller turbulent structures than lt , so DES length scale should be defined as: 1=2

lDES ¼ minðlt ; C DES DÞ

Table 1 Engine parameters.

ð9Þ

where D ¼ maxðDi Þ is the maximum length scale of the computational cell and C DES is a modeling constant. This idea successfully shifts the standard RANS SST model to a DES formulation. Additionally, this work utilizes a Delayed-DES method to keep LES region out of the boundary layer and ease the modeled-stress depletion problem. The dissipation term could be written as:

 b kx  F DES DkDES ¼ q F DES


 lt ¼ max ð1  F SST Þ; 1 C DES D

ð10Þ ð11Þ

where F SST equals the blending function in SST RANS model. However, the DDES-SST turbulent model requires high quality grids and the first grid point should be located within the viscous sublayer (y+ < 5) in RANS. The LES region in DES also requires a sufficiently fine isotropic grid resolution to support subgrid-scale model to resolve the isotropic small scale turbulent flow. It is challenge to meet all these requirements in the complex geometries especially with moving boundaries [33]. Spalart [44] recommended to locate the first grid point at y+ = 5 and to use a ratio of yj+1/yj = 1.3 for the first attempt. 2.2. Simulation setup The flow field and heat transfer in a motored Direct-Injection Spark-Ignition (DISI) optical engine was simulated using commercial CFD software FLUENT. We focused on the transient velocity and thermal boundary layers on the cylinder head in the tumble plane (x = 0 in Fig. 1) during the compression stroke. The optical engine was simplified from Volkswagen MAGOTAN and the fuel injector and the spark-plug were removed from the pent-roof combustion chamber. The engine specifications and valve timings are summarized in Table 1. We only simulated a half of the symmetric geometry for the sake of computing resources. In this paper, 0 crank angle degree (CAD) refers to top dead center at the beginning of the intake stroke. According to the requirements of DDES-SST model, the computational mesh was constructed in The Integrated Computer Engineering and Manufacturing code for Computational Fluid Dynamics (ICEM CFD) with the hybrid mesh strategy as shown in Fig. 2(a). The ports and the upper part of the combustion chamber utilized unstructured tetrahedral grids and the cylinder area

Parameters

values

Bore Stroke Connecting rod Compression ratio TDC clearance height Intake valve opens Intake valve closes Exhaust valve opens Exhaust valve closes Engine speed

82.5 mm 90 mm 168 mm 9.7:1 1.4 mm 22 CAD 241.9 CAD 491.3 CAD 706.6 CAD 800 rpm

utilized structured hexahedron grids. The grid reference size was selected as 1.3 mm to balance the computational accuracy and efficiency. The elements on the cylinder head were refined with 10 prism layers in the near-wall region. The first grid point spaced 0.1 mm from the wall, and the growth ratio was 1.1. Fig. 2(b) shows the dimensionless height of the first grid point at x = 0, y = 8 mm versus crank angle degree (CAD). The red dashed line represents the mesh requirement of DES approach (y+ < 5), so the current computational mesh satisfied the requirement in most regions. Although y+ of three piston positions close to TDC exceed the dashed line, the current computing resource cannot support further mesh refinement and the calculation accuracy is acceptable. The Reynolds number (Re ¼ us  B=m) of the in-cylinder boundary layer flow is proportional to y+ and it is raised from 1500 to 5000 during the compression stroke. A measure Mðx; tÞ [45] of the turbulence resolution is proposed to quantitatively validate the mesh resolution in LES region, and it is defined as follows:

Mðx; tÞ ¼

K res ðx; tÞ K res ðx; tÞ þ K sgs ðx; tÞ

K res ðx; tÞ represents the turbulent kinetic 1 ðu2rms þ 2rms þ w2rms Þ, K sgs ðx; tÞ represents 2

ð12Þ

energy calculated from the turbulent kinetic energy calculated from LES subgrid-scale model. So Mðx; tÞ shows the ratio of resolved turbulent kinetic energy to the amount of total turbulent kinetic energy. In general, a LES calculation is considered to be good, if the resolved structures contain at least 80% of turbulent kinetic energy (Mðx; tÞ > 0:8) [46]. Fig. 3 shows the Mðx; tÞ phase-averaged distribution on mid-section plane at BDC. Mðx; tÞ is higher than 0.8 in most areas, so the mesh resolution is fine enough in LES region. Overall, the whole computational domain contains approximately 1.4 million cells at the bottom dead center (BDC), and about 0.65 million cells at TDC.

v

1192

X. Fan et al. / International Journal of Heat and Mass Transfer 120 (2018) 1189–1199

Fig. 2. Computational mesh setting (a) mesh at BDC; (b) first grid y+ at x = 0, y = 8 mm.

Fig. 3. M(x, t) phase-averaged on mid-section plane at 180CAD.

The initial temperature within the cylinder is 318 K. The initial temperature of the intake port and the exhaust port are 313 K and of 333 K, respectively. The wall temperature on the cylinder head is fixed at 348 K and all the other region is fixed at 318 K. The initial pressure of the inlet and the outlet are 325 Pa (gauge pressure) and of 1325 Pa, respectively. The values of temperature and pressure in boundary condition are the same as the initial condition. The PISO algorithm was used to solve the velocity-pressure coupling. The time scheme is second order implicit and the time step varies between 13 and 5.2 ls. For momentum equation, the bounded center differencing is used. The discretization scheme applied for energy and turbulent kinetic energy is second order upwind. The interpolation scheme for pressure is PRESTO!, the interpolation scheme for density is second order upwind. The case begins at 380 CAD and runs through a total of 32 continuous engine cycles. The first 2 cycles are discarded to ignore the impact of initial conditions. The simulation was performed on a 20-core workstation and required approximately 96 h for each cycle. 2.3. Comparison with experimental results A set of two-dimensional three-component particle image velocity (2D3C-PIV) experiments were made on an optical DISI engine to validate the simulation results. The experiments were conducted with Lavision high speed 3D-PIV system, the schematic of the experiment is shown in Fig. 4. The three-dimensional velocity

Fig. 4. Schematic of the optical engine set-up and the PIV system for tumble measurements.

field on the tumble plane was measured during the intake and compression strokes (30–330 CAD). A dual oscillator single head 527 nm laser (LDY 300) was operated at a frequency of 2.4 kHz and the laser sheet thickness is 1 mm. 1 lm oil tracker particles were captured by two CCD cameras (HighSpeedStar 6) with the full 1024  1024 pixels chip. The interrogation window in postprocessing is 64  64 pixels. The final view size is 64.8  59.4 mm and the spatial resolution of the velocity vector is 1.03 mm. The calibration error of CCD camera is 0.029 pixel, the maximum theoretical error of three velocity components is less than 0.72 m/s at BDC. The simulation results and experimental data are compared base on the macro scale flow field rather than the near-wall region, because it is challenging in the experiment to capture the unsteady velocity boundary layer in a complex cylinder head structure. The ensemble-averaged 3D velocity vectors of PIV and DDES at 240 CAD and 270 CAD are shown in Fig. 5 and the white point in the figures represents the position of the tumble center. The red point in Fig. 5a represents the origin of the coordinate system. The experimental data were averaged from 100 cycles to reduce the experimental uncertainty due to cycle-to-cycle variation and the simulated data were averaged from 30 cycles. The relevance index (RI) proposed in Ref. [47] is employed to quantitatively analyze the resemblance between PIV and DDES velocity fields, and it is defined as follows:

X. Fan et al. / International Journal of Heat and Mass Transfer 120 (2018) 1189–1199

1193

Fig. 5. Validation of the DDES velocity fields against PIV results. (a) 240 CAD DDES; (b) 240 CAD PIV, RI = 0.91 at 240 CAD; (c) 270 CAD DDES; (d) 270 CAD PIV, RI = 0.94 at 270 CAD.

RIu;v ¼

ðu; v Þ kuk  kv k

ð13Þ

RIu;v ranges from 1 to 1.RIu;v = 1 means that the two velocity fields are identical; RIu;v = -1 indicates the two velocity fields have the same magnitude but the opposite flow direction; and RIu;v = 0 means the two velocity field are completely unrelated. The DDES and PIV results show good agreement in the aspects of the position of the tumble center, the velocity magnitude distribution and the flow direction. The RI is 0.91 at 240 CAD and 0.94 at 270 CAD.

3. Results and discussion 3.1. Velocity boundary layer Wall-bounded flow is influenced by the fluid viscosity in the near-wall region, which is defined as a 3 mm-thick layer from the wall in this paper. The typical point is selected at x = 0 mm, y = 8 mm on the cylinder head surface as shown in Fig. 5a. The ensemble-averaged velocity fields and the velocity fluctuation fields are analyzed below. 3.1.1. Temporal evolution of the velocity boundary layer Fig. 6 shows the evolution of the ensemble-averaged velocity profiles of stream direction during the compression stroke. The velocity decreases from 180 to 240 CAD because of the intake valve closing. After the intake valve fully closed at 241.9 CAD, the velocity rapidly increases to around 3 m/s and the momentum is due to

the upward piston. After 320 CAD, the velocity gradually decreases as the piston slows down. At the same time, the gas momentum dissipates because of the large-scale tumble flow breaks down to small-scale flow structures at the end of the compression stroke. The velocity magnitude is about 2 m/s at the TDC. Affected by the large-scale tumble flow, the flow direction near the cylinder head surface should be negative as shown in Fig. 5(a). However, a positive velocity appears at 220 CAD as shown in Fig. 6 (a). Analyzing the reason, there would form a small gap between the intake valve and the valve seat before the complete close of the intake valve, and this small gap will lead to a high-speed low-pressure region which results in a local reversal flow. Besides, boundary layer separation occurs within 0.5 mm spacing from the wall at 210 CAD and 230 CAD, the near-wall transient velocity field at 210 CAD in the 24th cycle is shown in the inset of Fig. 6(a). The top boundary refers to the cylinder head surface and the white dot on the wall indicates the separation point of the boundary layer. Therefore, the stability of the boundary layer is significantly influenced by the moving valves. In addition to the flow separation, sub-millimeter vortex structures randomly occur in the near-wall region. Fig. 7 shows the movement of a near-wall vortex from 195 to 205 CAD in the 28th cycle. The top boundary refers to the cylinder head surface. The counterclockwise vortex generated in the boundary layer around 195 CAD and moved towards the wall. It finally disappeared after hitting the wall at 205 CAD so the lifetime lasted for more than 10 CAD. Similar vortex structures were observed in the micro-PIV experiment in ICEs [11,14] and in plate boundary layer flow [48,49]. These vortex structures can contribute substantially to heat and mass transfer in the boundary layer.

1194

X. Fan et al. / International Journal of Heat and Mass Transfer 120 (2018) 1189–1199

Fig. 6. Ensemble-averaged velocity profiles of stream direction at selected CADs. (a) Velocity profiles from 180CAD to 240 CAD; (b) Velocity profiles from 260 CAD to 360 CAD.

Fig. 7. Near-wall vortex movement from 195 to 205 CAD in the 28th cycle. (a) 205 CAD; (b) 210 CAD; (c) 215 CAD.

3.1.2. Anisotropic turbulence in the near wall region The flow in the near-wall region of ICEs is a complex wall turbulent flow and the turbulence intensity is quantified by the fluctuation of the ensemble-averaged velocity. Fig. 8 shows the evolution of the stream and normal fluctuation intensities during the compression stroke. Generally, the stream and normal fluctuation intensities gradually increase as the piston moves upward, but there is another peak at 240 CAD. The reason is that the closing of the intake valve leads to local reversal flow and momentum dissipation. In the normal direction, the velocity fluctuation increases almost linearly as a function of the distance from the wall and the maximum is about 1 m/s at TDC. In the stream direction, the tendency of the velocity fluctuation profiles is similar to that of the mean velocity profiles, which grows rapidly within 0.5 mm spacing from the wall and grows gently after that. The stream

fluctuation intensities are on the same order as the ensembleaveraged velocities and almost twice as that in the normal direction. Reynolds stresses represent the mean momentum fluxes by turbulence and generally are larger than viscous stress. Fig. 9 shows the Reynolds stress in the near-wall region at 180 and 300 CAD. The two normal stresses gradually increase with the increase of the distance from the wall while the shear stress has no significant change. More specifically, Reyy is much larger than Rezz and Reyz throughout the compression stroke. Reyz is always close to zero and could be negative or positive. The results reflect that the turbulence in the near-wall region is anisotropic, and this is a major difference from the isotropic turbulence in the engine core region. Hence the mesh strategies should be different in the near-wall region and engine core region, also known as zonal approach [50].

Fig. 8. Velocity fluctuation profiles at selected CADs. (a) Fluctuation intensity of the stream direction; (b) Fluctuation intensity of wall-normal direction.

X. Fan et al. / International Journal of Heat and Mass Transfer 120 (2018) 1189–1199

1195

Fig. 9. Reynolds stresses distribution in near-wall region: the normal stress in the stream direction (y): Reyy ¼ v 02 , the normal stress in the normal direction (z): Rezz ¼ w02 , the shear stress on the yoz plane: Reyz ¼ v 0 w0 . (a) Reynolds stresses at 180 CAD; (b) Reynolds stresses at 300 CAD.

where j ¼ 0:41 is the von Kármán constant and B = 5.2 is the loglaw constant. Fig. 10 shows the dimensionless velocity profiles at selected CADs during the compression stroke. In this study, the length of the near-wall velocity profiles is fixed at 3 mm, so the starting point and the length of the dimensionless curves are different at different piston positions. Overall, the simulated results agree well with the standard wall function in the viscous sublayer (y+ < 5) and

deviate in the log layer (y+ < 30). The buffer layer (5 < y+ < 30) was not taken into account because of irregularity. The thickness of the viscous sublayer is defined as y+ = 5 according to the law of the wall, but the thickness of the viscous sublayer length increases as the CAD increases. The viscous sublayer thickness is about y+ = 5 at 220 CAD, but it exceeds y+ = 30 close to the TDC (360 CAD). In the log layer, the dimensionless velocity increases as the piston moves upward. So there is a negative offset between the simulated curves and the standard wall function at the beginning of the compression stroke, and the offset changes to positive later. In addition, the trend of the u+ profile is similar to the log law, but there is a decline tendency on the outside at several crank angle degrees due to the complex in-cylinder flow. The standard wall function was derived under ideal condition such as wall-parallel flow, quasi-steady flow, zero-pressure gradient. The in-cylinder flow is three-dimensional unsteady with a complex combustion chamber structure and a moving piston and valves, and these differences results in the inconsistency between the in-cylinder boundary layer velocity profiles with the law of the wall. Moreover, the in-cylinder velocity boundary layer is formed by the shear effect of the large-scale tumble flow, hence, it is affected by the evolution of tumble flow including the formulation, moving and breaking process. Fig. 11 shows the evolution of the viscous scale and the friction velocity. The viscous scale increases and the friction velocity decreases at the beginning of the compression stroke, which is mainly caused by the decrease of the velocity gradient at the wall before flow reversing. Both parameters show strong irregularity at the instant of the intake valve closing. After the near-wall flow is relatively stabilized, the in-cylinder gas is compressed as the piston moves upward and

Fig. 10. The dimensionless velocity profiles at selected CADs. The dashed line refers to the widely used standard wall function based on Eq. (4).

Fig. 11. Evolution of the viscous scale and the friction velocity.

3.2.1. Dimensionless velocity profiles The dimensionless ensemble-averaged velocity profiles are compared with the commonly used velocity law of the wall, which was derived under the condition of two-dimensional steady and fully developed flow over a plate. The law of the wall is usually non-dimensionalized using the viscous scale dv and the friction velocity us :

dv ¼

mw

ð14Þ

us

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  lw @ u us ¼ q @y w

ð15Þ

w

where qw and lw are the gas density and dynamic viscosity at the wall, mw ¼ lw =qw : Then, the distance from the wall and velocity can be expressed according to the parameters above yþ ¼ y=dv , uþ ¼ u=us . The law of the wall consists of two layers: the viscous sublayer and the logarithm layer:

(

þ

u ¼

yþ 1

j

if yþ < 11 þ

lnðy Þ þ B if yþ P 11

ð16Þ

1196

X. Fan et al. / International Journal of Heat and Mass Transfer 120 (2018) 1189–1199

leads to a rapid decrease in the kinematic viscosity, which finally results in a gradual decrease of the viscous scale and the friction velocity at the end of the compression stroke. 3.2. Thermal boundary layer The following discussion focuses on the thermal boundary layer. The complex local wall heat transfer process is analyzed based on the ensemble-averaged and fluctuation temperature and velocity distribution in the near-wall region. The same region as the previous discuss for velocity boundary layer is used. 3.2.1. Temporal evolution of the thermal boundary layer Fig. 12 shows the ensemble-averaged temperature profiles during the compression stroke. The intake temperature is 313 K and the wall temperature on the cylinder head is fixed at 348 K. The gas temperature is lower than the wall temperature at the beginning of the compression stroke, and becomes higher after 270 CAD because the gas is heated up. The temperature gradient at the wall increases towards TDC. The temperature profile varies non-monotonically from the wall at 260 CAD. The gas temperature slightly increases within 0.2 mm from the wall and then decrease. The temperature variation is caused due to the advection by the velocity field. This behavior indicates that the large-scale flow structures have a great impact on the temperature distribution in the near-wall region. The classic standard wall function cannot reflect non-monotonic variation of the temperature boundary layer [31]. For all piston positions, the temperature gradient tends to zero with the increase of the distance from the wall and the temperature stabilized at around 1 mm spacing from the wall. The thermal boundary layer thickness is defined as the distance between the wall and the point where the gas temperature reaches to 99% of the main stream gas temperature (calculated in the range from 1 mm to 3 mm spacing from the wall in this study). Fig. 13 shows the evolution of thermal boundary layer thickness during the compression stroke. The thickness increases as the excess temperature trends to zero and reaches to a peak value of 1.7 mm at 260 CAD. After that, the thermal boundary layer thickness drops sharply around 270 CAD because of the gas temperature and the wall temperature are almost the same. Then the thickness increases again as the temperature gradient reverses and the excess temperature increases. At around 310 CAD, the thickness decrease gradually because of the local wall temperature gradient increases. Although the mean temperature profiles are much more regular than the mean velocity profiles, the temperature fluctuation in the near wall region could not be ignored as shown in Fig. 14. Each

Fig. 12. Ensemble-averaged temperature profiles at selected CADs.

Fig. 13. Evolution of the thermal boundary layer thickness.

curve has a peak value within 0.5 mm spacing from the wall, indicating that the temperature fluctuation firstly increases and then decreases with the increase of the distance from the wall, and the fluctuation intensity is almost constant after 2 mm spacing from the wall. As the piston upwards, the peak value increases and the position of the peak moves towards the wall. The maximum is around 65 K at 0.2 mm spacing from the wall at TDC. The temperature fluctuation is associated with the convective heat flux in the in-cylinder flow, which will be discussed in Section 3.2.3. 3.2.2. Dimensionless temperature profiles The ensemble-averaged temperature profiles were nondimensionalized and compared with the universal thermal law of the wall. The thermal law of the wall was derived under the same ideal conditions as the momentum law of the wall with a constant inlet temperature and wall temperature. The friction temperature is defined as:

Ts ¼

  kw @@yT 

ð17Þ

w

qw cp us

where kw is the thermal conductivity at the wall and cp is the specific heat capacity of the gas. The dimensionless temperature is T þ ¼ ðT  T w Þ=T s . The near-wall temperature profiles can be described by the following formulas:

(

Tþ ¼

Pryþ

if yþ < 11 þ

Pr t ðj lnðy Þ þ BÞ þ P 1

if yþ P 11

ð18Þ

Fig. 14. Temperature fluctuation profiles at selected CADs during the compression stroke.

X. Fan et al. / International Journal of Heat and Mass Transfer 120 (2018) 1189–1199

where j ¼ 0:41 is the von Kármán constant and B = 5.2 is the loglaw constant. Pr is the laminar Prandtl number, Pr t is the turbulent Prandtl number, and the term P is a function of Pr and Prt . Fig. 15 shows the dimensionless temperature profiles, y+ was plotted in logarithmic scale and T+ was divided by the laminar Prandtl number. The dimensionless temperature profiles agree well with the thermal law of the wall in the viscous sublayer, although the curve has an obvious positive offset from the dashed line at 360 CAD because the first grid point exceeds y+ = 5. In addition, all the curves tend to be horizontal in the log layer and gradually increase as the piston moves upward. T+/Pr reaches to a peak value of around 12 at TDC. Fig. 16 shows the temporal evolution of the friction temperature and Prandtl number. The friction temperature is inversely proportional to the friction velocity according to Eq. (17), so it increases at the beginning of the compression stroke mainly due to the decrease of the friction velocity while the gas temperature is almost constant. The temperature gradient direction reverses around 260 CAD because the gas is compressed by the moving piston, so the friction temperature decreases in this stage. After that, the friction temperature gradually increases because of the increase of the local wall heat flux and reaches to a peak of about 42 K at TDC. The Prandtl number reflects the ratio of momentum diffusion and heat diffusion in the near-wall flow, and it is calculated from the dynamic viscosity, thermal conductivity and specific heat capacity of the gas at wall-bound temperature using Pr ¼ lcp =k. It can be seen that the Prandtl number rises from 0.705 to 0.727 during the compression stroke in Fig. 16. The Prandtl number was often set to be a constant in most flow and heat transfer simulations in ICEs and care should be taken in using such simplification in future simulations. 3.2.3. Wall heat transfer A main purpose of studying transient velocity and thermal boundary layers in ICEs is to understand the complex heat transfer process. In this work, the local wall heat flux is calculated based on the Fourier law, thus, thermal conductivity is calculated from Eq. (19) and the temperature gradient at the wall is estimated according to the ensemble-averaged temperature profiles. Fig. 17(a) shows the evolution of local wall heat flux during the compression stroke. A negative value means the in-cylinder gas absorbs heat from the cylinder head and a positive value is just the opposite. The in-cylinder gas release heat after 260 CAD because the excess temperature changes to positive. The local heat flux reaches to a peak around 70 kW/m2 at TDC, and the maximum heat absorption is about-5 kW/m2 at BDC.

Fig. 15. The dimensionless temperature profiles at selected CADs during the compression stroke. The bold dashed line refers to the viscous sublayer according to Eq. (18). The buffer layer and the log layer are not plotted because of a lack of proper experience parameters.

1197

Fig. 16. Temporal evolution of the friction temperature and Prandtl number.

The heat transfer coefficient is defined as

h¼

qw Tv  Tw

ð19Þ

where T v is the cylinder-averaged temperature and T w is the wall temperature of 348 K. As shown in Fig. 17(b), the heat transfer coefficient decreases from 100 W/(m2 K) at BDC to 30 W/(m2 K) around 260 CAD. The reason is that the excess temperature and the local wall heat flux decrease concurrently at this stage. After 270 CAD, the gas temperature is higher than 348 K, the excess temperature and local wall heat flux increase simultaneously which results in an increase of heat transfer coefficient after 280 CAD. The maximum heat transfer coefficient is about 200 W/(m2 K) at TDC. As mentioned earlier (in Section 3.2.1), the direction of temperature gradient reverses at around 260 CAD, so another peak at 270 CAD is caused by a large wall temperature gradient and the small temperature difference between the engine core region and the wall. Heat transfer mechanism in the near wall region of ICEs includes gas-phase convection and fuel film conduction if the influence of high-temperature gas and soot radiation is ignored [4]. Therefore, the one-dimensional heat flux equation is written as [51]

q ¼ qcp hv 0 T 0 i þ k

@T @z

ð20Þ

where the first term in Eq. (20) is the convective term which is also known as the turbulent heat flux and the second term is the conductive term. More specifically, v 0 and T 0 represent the wallnormal direction velocity fluctuation and the temperature fluctuation. v 0 T 0 reflects to the correlation between the velocity and thermal boundary layers. q refers to the gas density and cp refers to the specific heat capacity, these two parameters are treated as a function of the distance from the wall. Fig. 18 shows the spatial distribution of turbulent heat flux at several piston positions after 300 CAD. The turbulent heat flux increases as the piston moves upward and the local wall heat flux increases at the same time. For any piston position, the turbulent heat flux profile has a peak value near the wall, and the peak value increases and moves towards the wall as the piston moves upward. The black dashed line connects the peak point of all the curves and shows the moving tendency of the peak location. According to Eq. (20), the turbulent heat flux distribution depends strongly on the correlation between the wall-normal velocity fluctuation and temperature fluctuation [52]. Compared with the temperature fluctuation peak point in Fig. 14, the peak location of the turbulent heat flux is farther from the wall. Generally, most peaks are located within the buffer layer which reflects that the buffer layer is the main region for turbulence heat transfer. The maximum turbulent heat flux is about 70%85% of the local wall heat flux, and this data proved that forced convection heat transfer is the main heat transfer mechanics in ICEs.

1198

X. Fan et al. / International Journal of Heat and Mass Transfer 120 (2018) 1189–1199

Fig. 17. Evolution of local wall heat flux and heat transfer coefficient. (a) Temporal evolution of heat flux during the compression stroke, (b) Temporal evolution of heat transfer coefficient during the compression stroke.

Fig. 18. Turbulent heat flux profiles at selected CADs.

moves upward, but the evolution of the thermal boundary layer thickness is non-monotonic under the combined impact of excess temperature and temperature gradient, the peak thickness is 1.7 mm at 260 CAD. The dimensionless temperature profiles agree well with the thermal law of the wall in the viscous sublayer and tend to be uniform after that. In addition, the friction temperature and Prandtl number increase concurrently after the near-wall flow is relatively stabilized. The in-cylinder gas firstly absorb heat from the cylinder head and release heat after 260 CAD, the peak local wall heat flux is around 70 kW/m2 at TDC. The heat transfer coefficient decreases before 260 CAD and gradually increases to a peak around 200 W/ (m2 K) at TDC. The spatial distribution of turbulent heat flux shows a strong link between the velocity boundary layer and thermal boundary layer, and the peaks located in the buffer layer occupy about 70%–85% of the local wall heat flux, which indicates that the forced convection heat transfer is the main heat transfer mechanic in internal combustion engines.

4. Conclusions Acknowledgments In this work, the transient velocity and thermal boundary layers with wall heat transfer during the compression stroke in a motored SIDI-engine were predicted using DDES-SST model. The simulated results of the in-cylinder macro-flow field agreed well with the high-speed 2D3C-PIV experimental data in the tumble plane. We firstly analyzed the temporal evaluation of the velocity and thermal boundary layers, compared the dimensionless curves with the standard wall functions, and predicted the evolution of local wall heat flux and the turbulent heat transfer in the near-wall region. From these results, the following conclusions can be drawn: The temporal evolution of the ensemble-averaged stream velocity in the near-wall region is significantly influenced by the moving piston and valves during the compression stroke. Flow separation and sub-millimeter vortex structures randomly occur in the velocity boundary layer. The significant changes of velocity fluctuation and Reynolds stress in different directions proved that the turbulence in the near-wall region is anisotropic. The dimensionless velocity profiles agree well with the law of the wall in the viscous sublayer and deviate from the log layer, these differences come from the temporal complex flow and the variable fluid properties in ICEs. The thickness of the viscous sublayer gradually extends and the dimensionless velocity in loglayer increases as the piston moves upward. In addition, the viscous scale and the friction velocity decrease concurrently after the near-wall flow is relatively stabilized. The ensemble-averaged temperature and the temperature gradient in the near-wall region gradually increase as the piston

The study is financially supported by the National Science Fund for Distinguished Young Scholars (No. 51525603). Conflict of Interest Statement We declare that we have no conflict of interest to this work. The manuscript has been read and approved by all named authors and there are no other persons who satisfied the criteria for authorship but are not listed. We confirm that the order of authors listed in the manuscript has been approved by all of us. References [1] M. Schmitt, C.E. Frouzakis, A.G. Tomboulides, Y.M. Wright, K. Boulouchos, Direct numerical simulation of the effect of compression on the flow, temperature and composition under engine-like conditions, Proc. Combust. Inst. 35 (3) (2015) 3069–3077. [2] J.B. Heywood, Internal combustion engine fundamentals, McGraw-Hill, New York, 1988. [3] J. Chang, O. Güralp, Z. Filipi, D. Assanis, T.-W. Kuo, P. Najt, R. Rask, New heat transfer correlation for an HCCI engine derived from measurements of instantaneous surface heat flux, in: SAE Paper, 2004. [4] Z. Han, R.D. Reitz, A temperature wall function formulation for variabledensity turbulent flows with application to engine convective heat transfer modeling, Int. J. Heat Mass Transf. 40 (3) (1997) 613–625. [5] G. Borman, K. Nishiwaki, Internal-combustion engine heat transfer, Progress in Energy and Combust. Sci. 13 (1) (1987) 1–46. [6] I. Celik, I. Yavuz, A. Smirnov, Large eddy simulat ions of in-cylinder turbulence for internal combustion engines: a review, Int. J. Engine Res. 2 (2) (2001) 119– 148.

X. Fan et al. / International Journal of Heat and Mass Transfer 120 (2018) 1189–1199 [7] M.J. Hall, F.V. Bracco, Cycle-resolved velocity and turbulence measurements near the cylinder wall of a firing S.I. engine, in: SAE paper, 1986. [8] D.E. Foster, P.O. Witze, Velocity measurements in the wall buondary layer of a spark-ignited research engine, in: SAE paper, 1987. [9] P.H. Pierce, J.B. Ghandhi, J.K. Martin, Near-wall velocity characteristics in valved and ported motored engines, in: SAE paper, 1992. [10] A.Y. Alharbi, High-speed high-resolution vector field measurements and analysis of boundary layer flows in an internal combustion engine, University of Michigan, 2010. [11] A.Y. Alharbi, V. Sick, Investigation of boundary layers in internal combustion engines using a hybrid algorithm of high speed micro-PIV and PTV, Exp. Fluids 49 (4) (2010) 949–959. [12] V. Sick, High speed imaging in fundamental and applied combustion research, Proc. Combust. Inst. 34 (2013) 3509–3530. [13] L. Lu, V. Sick, High-speed particle image velocimetry near surfaces, J. Visualized Exp. 76 (2013). [14] C. Jainski, L. Lu, A. Dreizler, V. Sick, High-speed micro particle image velocimetry studies of boundary-layer flows in a direct-injection engine, Int. J. Engine Res. 14 (3) (2012) 247–259. [15] M. Koehler, D. Hess, C. Brücker, Flying PIV measurements in a 4-valve IC engine water analogue to characterize the near-wall flow evolution, Meas. Sci. Technol. 26 (12) (2015) 125302. [16] R. Lucht, M. Maris, Cars measurements of temperature profiles near a wall in an internal combustion engine, in: SAE paper, 1987. [17] R. Lucht, D. Dunn-Rankin, T. Walter, T. Dreier, S. Bopp, Heat transfer in engines: comparision of cars thermal boundary layer measurements and heat flux measurements, in: International Congress & Exposition, 1991. [18] M.E. Cundy, Development of high-speed laser diagnostics for the investigation of scalar heterogeneities in engines, University of Michigan, 2012. [19] M. Cundy, P. Trunk, A. Dreizler, V. Sick, Gas-phase toluene LIF temperature imaging near surfaces at 10 kHz, Exp. Fluids 51 (5) (2011) 1169–1176. [20] E.J. Lyford-Pike, J.B. Heywood, Thermal boundary layer thickness in the cylinder of a spark-ignition engine, Int. J. Heat Mass Transf. 27 (10) (1984) 1873–1878. [21] K.H. Tran, P. Guibert, C. Morin, J. Bonnety, S. Pounkin, G. Legros, Temperature measurements in a rapid compression machine using anisole planar laserinduced fluorescence, Combust. Flame 162 (10) (2015) 3960–3970. [22] O. Nijeweme, J. Kok, C. Stone, L. Wyszynski, Unsteady in-cylinder heat transfer in a spark ignition engine: Experiments and modelling, Proc. Inst. Mech. Engineers, Part D: J. Automobile Eng. 215 (6) (2001) 747–760. [23] E. Gingrich, J. Ghandhi, R.D. Reitz, Experimental investigation of piston heat transfer in a light duty engine under conventional diesel homogeneous charge compression ignition, and reactivity controlled compression ignition combustion regimes, SAE Int. J. Engines 7 (1) (2014) 375–386. [24] P.C. Ma, M. Greene, V. Sick, M. Ihme, Non-equilibrium wall-modeling for internal combustion engine simulations with wall heat transfer, Int. J. Engine Res. (2017) 1–11. [25] C.J. Rutland, Large-eddy simulations for internal combustion engines - a review, Int. J. Engine Res. 12 (2011) 1–31. [26] H. Schlichting, K. Gersten, boundary layer theory, 8th revised and, enlarged edition., Springer, Berlin, 2003. [27] C. Rakopoulos, G. Kosmadakis, E. Pariotis, Critical evaluation of current heat transfer models used in CFD in-cylinder engine simulations and establishment of a comprehensive wall-function formulation, Appl. Energy 87 (5) (2010) 1612–1630. [28] S. Keum, H. Park, A. Babajimopoulos, D.N. Assanis, D. Jung, Modelling of heat transfer in internal combustion engines with variable density effect, Int. J. Engine Res. 12 (6) (2011) 513–526. [29] F. Berni, G. Cicalese, S. Fontanesi, A modified thermal wall function for the estimation of gas-to-wall heat fluxes in CFD in-cylinder simulations of high performance spark-ignition engines, Appl. Therm. Eng. 115 (2017) 1045–1062.

1199

[30] M. Jia, E. Gingrich, H. Wang, Y. Li, J.B. Ghandhi, R.D. Reitz, Effect of combustion regime on in-cylinder heat transfer in internal combustion engines, Int. J. Engine Res. 17 (3) (2015) 331–346. [31] P.C. Ma, T. Ewan, C. Jainski, L. Lu, A. Dreizler, V. Sick, M. Ihme, Development and analysis of wall models for internal combustion engine simulations using high-speed micro-PIV measurements, Flow, Turbul. Combust. (2016) 1–27. [32] R. Jhavar, Using Large Eddy Simulation to Study Diesel DI-HCCI Engine Flow Structure, University of Wisconsin - Madison, Wisconsin, Mixing and Combustion, 2007. [33] C. Hasse, V. Sohm, B. Durst, Numerical investigation of cyclic variations in gasoline engines using a hybrid URANS/LES modeling approach, Comput. Fluids 39 (1) (2010) 25–48. [34] F. Hartmann, S. Buhl, F. Gleiss, P. Barth, M. Schild, S.A. Kaiser, C. Hasse, Spatially resolved experimental and numerical investigation of the flow through the intake port of an internal combustion engine, Oil & Gas Science and Technology – Revue d’IFP Energies nouvelles 71 (1) (2016) 1–129. [35] S. Buhl, F. Gleiss, M. Köhler, F. Hartmann, D. Messig, C. Brücker, C. Hasse, A Combined Numerical and Experimental Study of the 3D Tumble Structure and Piston Boundary Layer Development During the Intake Stroke of a Gasoline Engine, Flow, Turbul. Combust. (2016) 1–22. [36] J.M. Wallace, Highlights from 50 years of turbulent boundary layer research, J. Turbul. 13 (53) (2012) 1–70. [37] M.S. Shadloo, A. Hadjadj, F. Hussain, Statistical behavior of supersonic turbulent boundary layers with heat transfer at M1=2, Int. J. Heat Fluid Flow 53 (2015) 113–134. [38] M.S. Shadloo, A. Hadjadj, D.J. Bodony, F. Hussain, S.K. Lele, Effects of heat transfer on transitional states of supersonic boundary layers, center for turbulence research, Stanford University, Proc. Summer Program (2016) 175–184. [39] X. Wu, P. Moin, Transitional and turbulent boundary layer with heat transfer, Phys. Fluids 22 (8) (2010) 1–8. [40] M. Schmitt, C.E. Frouzakis, Y.M. Wright, A.G. Tomboulides, K. Boulouchos, Direct numerical simulation of the compression stroke under engine-relevant conditions: evolution of the velocity and thermal boundary layers, Int. J. Heat Mass Transf. 91 (2015) 948–960. [41] M. Schmitt, C.E. Frouzakis, Y.M. Wright, A.G. Tomboulides, K. Boulouchos, Investigation of wall heat transfer and thermal stratification under enginerelevant conditions using DNS, Int. J. Engine Res. 17 (1) (2015) 63–75. [42] M. Schmitt, C.E. Frouzakis, Y.M. Wright, A. Tomboulides, K. Boulouchos, Direct numerical simulation of the compression stroke under engine relevant conditions: Local wall heat flux distribution, Int. J. Heat Mass Transf. 92 (2016) 718–731. [43] M. Bellec, A. Toutant, G. Olalde, Large Eddy Simulations of thermal boundary layer developments in a turbulent channel flow under asymmetrical heating, Comput. Fluids (2016) 1–18. [44] P.R. Spalart, young person’s guide to detached-eddy simulation grids, 2001. [45] F.D. Mare, R. Knappstein, M. Baumann, Application of LES-quality criteria to internal combustion engine flows, Comput. Fluids 89 (2) (2014) 200–213. [46] S.B. Pope, Turbulent Flow, Cambridge University Press, UK, Cambridge, 2000. [47] K. Liu, D.C. Haworth, Development and assessment of POD for analysis of turbulent flow in Piston Engines, in: SAE paper, 2011. [48] D. Lengani, D. Simoni, M. Ubaldi, P. Zunino, F. Bertini, Turbulent boundary layer separation control and loss evaluation of low profile vortex generators, Exp. Thermal Fluid Sci. 35 (8) (2011) 1505–1513. [49] R. Dekou, J.M. Foucaut, S. Roux, M. Stanislas, J. Delville, Large scale organization of a near wall turbulent boundary layer, Int. J. Heat Fluid Flow 61 (2016) 12–20. [50] U. Piomelli, Wall-layer models for large-eddy simulations, Prog. Aerosp. Sci. 44 (6) (2008) 437–446. [51] F. Nicoud, numerical study of a channel flow with variable properties, Center for Turbul. Res. (1998) 289–310. [52] H. Zhao, A. Wei, K. Luo, J. Fan, Direct numerical simulation of turbulent boundary layer with heat transfer, Int. J. Heat Mass Transf. 99 (2016) 10–19.