© A. Wu et al., published by IFP Energies nouvelles, 2019

Symbols and/or notations

A_n, B_n: Fourier series expansion coefficients

C_p: Heat capacity

f: Filter width

G: Gaussian filter

h: Heat transfer coefficient

k: Thermal conductivity

k_resolved: Resolved turbulent kinetic energy

k_sgs: Subgrid turbulent kinetic energy

M: LES quality index

: Spatially averaged LES quality index

N: Number of cycles

n: Harmonic number

q_s: Surface heat flux

T_melt: Melting temperature

: Cycle-averaged surface temperature

: Cycle-averaged in-depth temperature

U_i: Velocity

: Fluctuating velocity

U_i,std: Standard deviation of the velocity

x_i: Position

Greek symbols

α: Thermal diffusivity

δ: Depth

θ: Crank angle degree

ρ: Density

ω: Angular velocity

Operators

〈∙〉: Ensemble averaging

: Spatial averaging

: Spatially filtered

Abbreviations

aTDCc: After Top Dead Center compression

CAD: Crank Angle Degree

CARS: Coherent Anti-Stokes Raman Scattering

CCV: Cycle-to-Cycle Variation

CFD: Computational Fluid Dynamics

CHT: Conjugate Heat Transfer

EVC: Exhaust Valve Closing

EVO: Exhaust Valve Opening

EXP: Experiment

FoV: Field of View

GDI: Gasoline Direct Injection

ICE: Internal Combustion Engine

IVC: Intake Valve Closing

IVO: Intake Valve Opening

LDV: Laser Doppler Velocimetry

LES: Large-Eddy Simulation

NWR: Near-Wall Region

NOx: Nitric Oxides

PIV: Particle Image Velocimetry

PTV: Particle Tracking Velocimetry

RANS: Reynolds-Averaged Navier Stokes

TDC: Top Dead Center

uniT: Uniform Temperature

1 Introduction

Health and environmental awareness have been driving more stringent Internal Combustion Engine (ICE) regulations on harmful emissions and engine performance. To meet such requirements, it is important to increase the ICE efficiency by reducing loss mechanisms, such as thermal loss across the cylinder wall. The boundary layer experiences large thermal gradients, and therefore most of the engine heat transfer occurs over this layer. This wall heat transfer affects the surface temperature and temperature profile in the Near-Wall Region (NWR), and is one of the major contributors to engine efficiency losses. In addition, the wall heat transfer affects combustion phasing [1, 2], flame quenching [3], engine knock, [4], and the formation of harmful nitric oxides (NOx) emissions [5]. Accurate simulations of ICEs require that the surface temperature and temperature profile in the NWR are predicted correctly.

Experimental studies have been performed since the 1980s to provide a better understanding of the boundary layer development in ICEs. Optical measurement techniques, such as Laser Doppler Velocimetry (LDV) [6, 7], Particle Image Velocimetry (PIV) [8, 9], hybrid micro-PIV and Particle Tracking Velocimetry (PTV) [10], and micro-PIV [11], have been developed and applied to ICE momentum boundary flows. Schlieren photographs [12] and Coherent Anti-Stokes Raman Scattering (CARS) [13] both revealed thickening of the thermal boundary layer during the expansion stroke.

These experimental studies provide guidance to how the in-cylinder flow and heat transfer in ICEs need to be modeled. The momentum boundary layer is often modeled by the law of the wall [14], or the Werner and Wengle model [15], but improvements have been suggested [16, 17]. Ma et al. [16] found several shortcomings of the law of the wall model applied to a motored ICE flow. From these shortcomings, they developed a non-equilibrium wall model which includes the transient, convective, pressure gradient, and pressure work terms. Their model improved the velocity predictions in the NWR compared to the law of the wall. They also predicted non-monotonic behavior of the temperature profile in the NWR, but did not have experimental data to validate with. In a follow-up study, Ma et al. [17] extended their non-equilibrium wall model to a fired engine flow. The model includes the heat release rate in the energy equation. With their model, they improved velocity and surface heat flux predictions compared to the law of the wall model and empirical heat transfer correlations. However, the non-equilibrium wall model requires solving a set of partial differential equations, and therefore requires resolving the boundary layer thickness. This might be impractical in ICE simulations due to the increased mesh and computational requirements and to date, this improved wall model has not yet been implemented into Computational Fluid Dynamics (CFD), especially in the Large-Eddy Simulation (LES) framework.

On the contrary, modeling of the thermal boundary layer has been investigated and improved over decades to incorporate mainly the variable density effect on heat transfer [18–23]. Due to the strong temperature gradient in the NWR, Angelberger et al. [18] modified the thermal law of the wall model to include non-isothermal effects on the variable fluid properties. Han and Reitz [19] included effects of variable density and turbulent Prandtl number in the wall model formulation, and their model was optimized in the Reynolds-Averaged Navier-Stokes (RANS) framework. In addition to the effects of variable density and variable turbulent Prandtl number, Keum et al. [20] also included variable viscosity in their wall model. Berni et al. [21] used mean thermodynamic properties from the near-wall cells in their thermal wall model with time-variable Prandtl number, which performed better for engines with smaller temperature gradients at the wall, compared with Angelberger or Han and Reitz models, which worked well for engines with higher wall temperature gradients. Šarić et al. [22, 23] developed a hybrid model using the Han and Reitz model, and used a more advanced RANS k–ζ–f turbulence model to improve wall heat flux predictions.

To complete the wall modeling, the thermal boundary condition at the engine wall is an important consideration. Typical wall modeling uses uniform and constant wall temperature as the thermal boundary condition, which affects the accuracy of the predicted heat flux [24]. However, experiments have revealed that there is significant spatial and temporal variation of the wall temperature, as well as Cycle-to-Cycle Variation (CCV), in ICE flows [25]. Negative surface heat flux from the wall to the hot gases has also been seen to occur during the engine cycle due to pressure work in the boundary layer, even with bulk gas temperatures being much higher than the wall temperature [26]. Thus, using a constant and uniform wall temperature as the thermal boundary condition is inaccurate and can be detrimental to engine heat transfer predictions.

Conjugate Heat Transfer (CHT) solves the coupled heat transfer between fluid and solid domains. Recently, CHT has been applied to engine modeling in conjunction with RANS turbulence models [27–33]. Three typical CHT modeling approaches are as follows. CHT can be used to model the coupled heat transfer between the engine solid components and liquid cooling jacket, while the in-cylinder flow is simulated separately using CFD [27, 29, 30, 32]. Some studies only considered the coupling of the in-cylinder gas and solid components [34, 35]. Or, CHT can be used to simulate the heat transfer between the engine solid components and the in-cylinder flow, while the liquid cooling jacket is simulated using CFD [28, 31, 33]. The CFD and CHT domains are coupled by the resulting wall temperature from the CHT simulation, and the heat transfer coefficient and fluid temperature from the CFD calculation.

Xin et al. [27] applied CHT to the solid and coolant components, excluding heat transfer in the valves, and simulated the in-cylinder flow using CFD. Their results show that the predicted heat flux from CHT can differ more than 20% from the heat flux obtained from uniform wall temperature boundary condition. Urip et al. [34] implemented the CHT model in KIVA-3V, but excluded the CHT calculation in the valves, and only simulated the engine cycle from Intake Valve Closing (IVC) to Exhaust Valve Opening (EVO). Significant heat transfer can occur during valve events, and the entire engine cycle should be simulated. Li and Kong [28] modeled the solid engine components and the in-cylinder flow separately, but enforced continuity between the fluid and solid domains. Their results show that CHT influenced the NOx, soot, and heat release rate predictions, due to CHT being able to predict spatially and temporally varying temperatures. Fontanesi et al. [29] examined the influence of uniform versus spatially varying heat flux obtained from CHT calculations, on the heat transfer predictions of a multi-cylinder engine. They found that only the CHT simulation was able to predict the correct engine temperature field. In addition, using a temperature-dependent thermal conductivity leads to improved temperature field predictions. Cicalese et al. [30] applied CHT to the cylinder head and valves of a Diesel engine. The valves were fixed due to their CHT model not being able to simulate moving parts. To account for valve motions, a variable thermal resistance was used between the valves and the valve seats. A separate piston model was created to account for the effects its motion on the heat transfer of the cylinder liner. They were able to predict heat rejection to the coolant within 3% error. Iqbal et al. [31] were able to predict temperatures within 10% of the experimental measurement at several locations using a CHT model for three engine load conditions. Wu et al. [32] found that the resulting spatially varying temperatures significantly influence the predicted velocity field of a Gasoline Direct Injection (GDI) engine. Leguille et al. [33] preserved the valve motion in their CHT simulations to better account for the cooling effects of the valves when they are in contact with the valve seats.

Misdariis et al. [35] coupled an LES with a solid simulation of the cylinder head with fixed valves. They were able to predict CCV in the instantaneous heat flux and the cycle-averaged heat flux. However, the impact of the cyclic variation of the heat flux on the surface temperature was low. The nonuniformity of the surface temperature affected the flow field and therefore the combustion process. The use of fixed valves is questionable since the influence of the moving valves on the predicted flow field cannot be neglected.

These efforts show that more accurate heat transfer predictions can be obtained with CHT, since it can provide better thermal boundary conditions at the gas-solid interfaces, compared to applying constant wall temperature boundary conditions in conventional in-cylinder flow modeling. Validation of the CHT approach with measured heat flux and surface temperature data has been done extensively. However, evaluation of the effects of CHT on the flow field has not yet been performed.

In this study, the effects of the CHT approach on the bulk flow field, as well as the near-wall flow field, and heat transfer predictions are evaluated against previously published experimental bulk flow, near-wall flow, and heat transfer data [9, 36]. Two LES computations are performed. The first case uses a fixed, uniform temperature boundary condition. The second case uses CHT to couple the in-cylinder flow and solid domains. For the second case, an additional liquid coolant simulation is performed. The CHT engine modeling approach is improved upon by including the solid cylinder head with its coolant circuit, moving valves, cylinder liner, moving piston, and spark plug components. First, the experimental and numerical methodology are discussed, followed by the method of analysis of the PIV and LES data. Next, a discussion of the statistical convergence, bulk flow, near-wall flow, and temperature fields, the surface temperature predictions from the CHT simulation, and heat flux profiles is provided.

2 Experimental methodology

The Transparent Combustion Chamber (TCC-III) engine was utilized in this study. It is a two-valve, four-stroke, single cylinder engine with a pancake chamber and centrally-installed spark plug in the cylinder head. The engine provides optical access through a quartz cylinder liner and a quartz flat piston window. Detailed description of the engine can be found in [37], and a schematic of the engine is shown in Figure 1.

Bulk flow measurement data from the central plane at y = 0 mm [36, 37] were used to start the validation of the simulation results of this study, as shown in Figure 2. The dataset analyzed was S_2013_10_24_01. The engine was motored at 1300 RPM, 40 and 98 kPa intake and exhaust MAP, with intake air heated to 45 °C. Every five Crank Angle Degree (CAD) from −360 to +330 CAD after Top Dead Center compression (aTDCc) was imaged, and 240 cycles were recorded. PIV processing was done in LaVision DaVis. The images were interrogated with a decreasing window size from 128 × 128 pixels, 50% overlap, to a final window size of 32 × 32 pixels at 50% overlap, resulting in a spatial resolution of 2.4 mm, oversampled to a final vector spacing of 1.2 mm.

In addition to bulk flow measurements, near-wall experimental flow and heat transfer data were also used to validate the simulated results [9, 36]. The dataset analyzed was S_2016_03_13_04. PIV was used to measure the in-plane velocity components in a vertical plane at 28 mm offset from the central tumble plane. The location of the near-wall PIV plane and heat flux probe are shown in Figure 3. The nominal Field of View (FoV) of the PIV plane was 5 × 6 mm, which allowed sampling of the core flow. The engine was motored at 1300 RPM, 40 and 98 kPa intake and exhaust MAP, with intake air heated to 80 °C. Every 10 CAD of the entire engine cycle was imaged, and 141 cycles were recorded. PIV processing was done in LaVision DaVis. The images were interrogated with a final window size of 32 × 32 pixels at 50% overlap, resulting in a spatial resolution of 250 μm, oversampled to a final vector spacing of 125 μm. The first vector from the wall was located at 250 μm.

The surface and in-depth temperatures were measured at every 0.5 CAD using a Medtherm heat flux probe, located at 35.5 mm from the cylinder axis, mounted flush with the cylinder head. The instantaneous heat flux was calculated from the measured temperatures using the method by Nijeweme et al. [26]:

(1)Here, q_s is the instantaneous surface heat flux at the cylinder head surface, k and α are the thermal conductivity and thermal diffusivity of the probe material, δ is the distance between the thermocouples of the in-depth and surface temperature measurement, is the cycle-averaged surface temperature, is the cycle-averaged in-depth temperature, A_n and B_n are the coefficients from the Fourier series expansion of the surface temperature trace, n is the harmonic number, and ω is the angular velocity.

3 Numerical methodology

Two simulations were conducted based on the thermal boundary condition applied. The first case uses uniform, constant temperature as the thermal boundary condition (hereafter called the uniT model). For the second case, CHT is used to calculate the surface temperature of the engine cylinder (hereafter called the CHT model). Simulations were performed using a commercial CFD software package (CONVERGE v2.4.18 [38]). Turbulence at the subgrid scales was modeled using the dynamic structure model [39]. The momentum boundary layer was modeled using the Werner and Wengle model [15], and the thermal boundary layer using the Han and Reitz model [19]. The computational domain for both simulations includes the intake and exhaust plenums, intake and exhaust ports, and the cylinder. The CHT model includes additional solid components of the engine, i.e., the cylinder head, the cylinder liner, exhaust and intake valves, valve seats, the full spark plug geometry, and the piston. The material properties of the solid components are listed in Table 1. The computational domain for the uniT and CHT models are shown in Figures 4 and 5. Eleven cycles were simulated using the CHT model, while 31 cycles were simulated using the uniT model. The first cycle was discarded from both models to remove initial condition bias. Cycles 2–11 of each model is analyzed in this study.

Orthogonal cubic meshes are used during run-time. In this study, a base mesh of 8 × 8 × 8 mm is used for the uniT and CHT models, with additional embedding of finer grids within the combustion chamber. The mesh gradually reduces from 8 mm to 1 mm from the plenums to the cylinder, with 0.5 mm mesh in the spark plug vicinity, at the valve top and bottom, and head, liner, and piston surfaces. The in-cylinder mesh topology of the uniT model is shown in Figure 6. In the CHT model, 0.5 mm mesh is used between the solid-solid and fluid-solid interfaces.

Inlet and outlet boundary conditions at the plenums were obtained from GT-Power following the procedure defined by [40], with intake temperature controlled at 80 °C. For the uniT model, a constant wall temperature of 353.16 K is applied to all surfaces. For the CHT model, the temperature of the top of the cylinder head, valve stems, and head groove is set to 353.16 K. Convective boundary condition is applied to the exposed sides and bottom of the cylinder head, with heat transfer coefficient h set at 10 W/m² K and ambient temperature at 298.16 K. The piston skirt, liner bottom, and inner liner surface are set to adiabatic. The piston bottom and the liner outer surface are exposed to a cooling air jet, so a convective boundary condition is applied, with h set at 200 W/m² K, and ambient temperature of 298.16 K.

In the CHT model, the motion of the piston and valves are preserved by coupling the heat transfer of the in-cylinder flow and the solid domains. A separate liquid coolant simulation was performed for the coolant jacket in the cylinder head. The turbulence model for the coolant jacket is the RANS RNG k-ε model [41]. Adaptive Mesh Refinement (AMR) was used to increase the grid resolution. The mesh topology is shown in Figure 7. The scalable wall function is used for the wall model, which ensures that proper log-law behavior is reproduced if grid refinement due to AMR puts the first cell in the buffer region [38]. A mixture of 50:50 by volume of ethylene glycol and water was fed through the inlet at 3.15e⁻⁴ m/s, with an inlet temperature of 353.16 K. The coolant wall temperature is initially assumed to be uniform, and was set at 353.16 K. This coolant model provides a spatially varying convective boundary condition for the CHT model. The CHT model was first run with the RANS RNG k-ε model. This CHT model provides an updated spatially varying wall temperature for the second iteration of the coolant model. This iterative process was performed three times until a converged boundary condition at the coolant wall was obtained. The RANS CHT simulation results provided the initial and boundary conditions for the LES CHT model.

In CHT, the time-scale for the heat transfer in the solid is several orders of magnitudes larger than the fluid time-scale [33]. This discrepancy would lead to an increase in the number of engine cycles needed before the solid components reached thermal equilibrium. To address this with less computational effort, CONVERGE uses a super-cycling method [38]. This method iterates between a fully-coupled transient fluid-solid solver, and a steady-state solid solver which calculates and updates the surface temperature at the fluid-solid interface.

4 Analysis

A velocity vector from PIV and LES is defined as the resolved velocity U_i. Each component of the velocity vector at each CAD θ and position vector x_i was ensemble averaged as

(2)where N is the number of cycles. The velocity can be decomposed into a fluctuating velocity component using Reynold’s decomposition:

(3)The standard deviation of the velocity U_i,std represents the CCV of the flow field, and was obtained using

(4)The ensemble and standard deviation values of the temperature and heat flux were obtained in the same manner as above.

Due to the difference in spatial resolution of the near-wall PIV and LES results, the near-wall PIV data is smoothed using a Gaussian spatial filter at each point (x_o, y_o). The Gaussian filter is defined as

(5)

The filtered velocity at point (x_o, y_o) can then be obtained using

(6)A filter width f of 1 mm is used to match the in-cylinder mesh size of 1 mm in the LES models. Figure 8 shows the effects of applying the Gaussian spatial filter on an individual cycle at −180 CAD aTDCc. The Gaussian filter averages out the small-scale features while the large-scale features are conserved. The magnitude of the velocity vectors are reduced after filtering since they are averaged with neighboring vectors. In addition, slight changes in the direction of the velocity vectors can be seen.

5 Results and discussion

The first step in validating the LES models is to compare ensemble averaged cylinder pressures of the simulations to the experimental ensemble averaged pressures as shown in Figure 9. The peak pressures are underestimated by both models, but are within 3% of the experimental peak pressures.

Multiple consecutive cycles have been simulated extensively in LES framework by others, mainly for CCV studies [42–44]. Statistical convergence of the simulations is therefore an important consideration. The ensemble average and standard deviation of the velocity field was obtained for increasing number of cycles (i.e., cycles 2–3, cycles 2–4, cycles 2–5, etc.). These velocities were then spatially averaged over the entire cylinder domain, and the magnitude of the spatially averaged velocity vector was obtained. Figure 10 shows the ensemble average and standard deviation of the spatially averaged velocity for increasing number of cycles in the statistical analysis. For both models, satisfactory statistical convergence of ensemble-averaged velocity is acquired after including five cycles (cycles 2–6) in the ensemble averaging process. At least seven cycles (cycles 2–8) are needed for the standard deviation of the velocity to reach statistical convergence. Figure 11 shows the ensemble averaged gas temperature for an increasing number of cycles, which was then spatially averaged over the cylinder domain. The ensemble average temperature of both models has reached statistical convergence after averaging cycles 2–3, while the standard deviation of the temperature becomes steady after cycles 2–8. It is also evident that the ensemble averaged gas temperature of the CHT model is higher than for the uniT model.

In LES, large scales are resolved while the subgrid scales are modeled. A quality index M can be used to measure the turbulence resolution of the LES model [45], and it is an indicator of the grid resolution [46]. It is defined as,

(7)where k_resolved is the turbulent kinetic energy of the resolved velocity field, obtained using

(8)k_sgs is the subgrid-scale kinetic energy, which is obtained from the dynamic structure model. M can take values between 0 and 1, where a value of 0 corresponds with direct numerical simulation (all turbulence scales are resolved), and a value of 1 corresponds with RANS (all turbulence is modeled). Typically, a value of 0.2 or less is recommended, which indicates that 80% or more of the turbulent kinetic energy is resolved by the LES model [45].

The quality index is calculated for both the uniT and CHT models over the engine cycle with increasing number of cycles in the ensemble averaging process, and spatially averaged over the entire cylinder domain. The spatially averaged quality index is shown in Figure 12. Evident from the figure is that as the number of cycles in the ensemble averaging process increases, decreases. The exception is cycles 2–3 in the CHT model which has low value. This could be due to how the first cycle of the LES CHT model was initialized using the RANS CHT results. This initialization could still influence cycle 2. Once more cycles were added to the averaging process, follows the same trend as the uniT cycle. Since this quality index has the influence of both the ensemble average and the fluctuating velocity , it can be used as another indicator of statistical convergence of the simulations. During the compression stroke, the integral length scale decreased due to the piston compressing the flow into a decreasing engine volume, and the valve opening events introduced additional small-scale structures in the flow [47]. The quality index reflects such a decrease in the integral length scale as it increases to a maximum during the compression stroke and during the valve opening events. This index indicates that during these events, the computational mesh is unable to resolve the small turbulent scales. At times of large M values, grid refinement is needed to resolve more of the turbulent kinetic energy.

Ensemble averaged bulk flow fields are shown in Figures 13–15 from the PIV measurement and the LES models at three CAD. Figure 13 shows the velocity field at −260 CAD aTDCc during the intake stroke. At this CAD, the spatially average MI values are equal to 0.10 for both models. This figure shows that the bulk flow is well captured by both models. The direction and magnitude of the large-scale turbulent intake jet are well captured, and the entrainment of the flow on both sides of the intake jet is also seen in both simulations. However, differences in the bulk flow field between the simulations are evident. The flow entrainment on the left of the intake jet (region A) is much stronger in the CHT model compared to the uniT model, and compares better with the magnitude of the flow in the PIV measurement. In addition, the area of the high-magnitude flow at the bottom-right corner (region B) is larger in the CHT model compared to the uniT model. The location of the vortex center on the right of the intake jet (region C) from both models compare qualitatively well with the PIV measurement, although slight differences can be seen between the models. On the other hand, the flow at the bottom-left corner near the piston (region D) compares much better between the uniT model and the experiment.

For quantitative analysis of the velocity field, the local Relevance Index (RI) and the local Kinetic Energy Index (KEI) have been calculated. The local RI compares the vector alignment of two velocity vectors [48], and is defined by

(9)Here, denotes the magnitude of the velocity vector. If the velocity vectors are perfectly aligned, RI will be equal to 1. If the velocity vectors are orthogonal to each other, RI will be equal to 0. And if the velocity vectors are opposite of each other, RI will be −1.

RI only compares the angle between two velocity vectors, while the kinetic energy is disregarded. To compare the kinetic energy between two velocity vectors, the local KEI can be used [49], and is calculated using

The local KEI can take on values equal or larger than 0. A KEI of 1 means that the kinetic energy is perfectly captured by the simulation, while a KEI of less than or larger than 1 means that the kinetic energy is either too low or too high. This index also indicates whether the velocity magnitude is well-captured by the simulation. The local RI and KEI values can be spatially averaged over the plane to obtain global RI and KEI.

In Figure 13, regions A and C are marked by low RI values in both models. These regions contain the vortical structures near the intake jet. The low RI values indicate that the flow direction differs from measurement due to differences in the vortex centers. It can be seen that in the uniT model, both regions encompass a larger area of low RI values compared to the CHT model. The flow direction in region D is better captured by the uniT model. This is due to the flow entrainment on the left side of the intake jet, which is much weaker in the uniT model. A larger part of the flow is entrained in the CHT model, and this causes flow in region D to move upwards. The RI values in region B is close to 1, indicating that the velocity direction is well captured by both models. The KEI values in this figure indicate that the CHT tends to underpredict the velocity magnitude in the bulk flow domain, when compared to the uniT model. The kinetic energy in the intake jet near the Intake Valve Opening (IVO) is overpredicted by both models. The region of large KEI values in the intake jet is much larger in the uniT model than the CHT model.

The ensemble averaged bulk flow fields at −30 CAD aTDCc are shown in Figure 14. The RI fields show that the velocity field is very similar, except for the regions close to the spark plug (regions A and B). The kinetic energy is also overestimated in these two regions, with higher kinetic energy values in the uniT model. This will have implications for combustion simulations, since the kinetic energy and the flow velocity are variables used in combustion modeling. At this CAD, the spatially averaged MI values are equal to 0.14 for the CHT model and 0.17 for the uniT model. The grid resolution may not be fine enough in regions A and B. Increasing the grid resolution in these regions will allow smaller flow structures to be resolved, leading to possible reduction in the discrepancy of the flow direction and kinetic energy.

The ensemble averaged bulk flow fields at 120 CAD aTDCc is shown in Figure 15. At this CAD, the spatially average MI values are equal to 0.07 for the CHT model and 0.09 for the uniT model. This velocity field corresponds to right before EVO, when the LES quality deteriorates quickly. The flow direction is well captured by both models with better match between the CHT and the bulk flow measurements. On the other hand, kinetic energy is overestimated by the uniT model in much of the FoV, while the CHT predicts smaller kinetic energy.

Figure 16 shows the temperature fields of the uniT and CHT models. Experimentally determined temperature fields are not yet available. It is evident that differences in the thermal boundary condition lead to differences in the bulk gas temperature predictions. At −260 CAD aTDCc, near the exhaust valve (region A), a hotter region exists in the CHT case compared to the uniT model. Region B in the intake jet of the CHT case is also at a higher temperature compared to the uniT model. Near the piston surface (region C), differences in the surface temperature also lead to different thermal behavior of the nearby gases in both models. Not shown here is the surface temperature of the spark plug. The CHT model predicts a higher surface temperature than the prescribed temperature of the uniT model. This leads to higher gas temperature in close proximity to the spark plug, which will lead to differences in the flame kernel development due to changes in heat transfer. The differences in bulk gas temperature predictions lead to differences in gas density, hence the velocity fields are also affected by the thermal boundary condition, as shown in Figure 13. The temperature fields at −30 CAD aTDCc are comparable to each other, although the CHT model predicts slightly higher temperatures compared to the uniT model. During the expansion stroke, at 120 CAD aTDCc, the effects of the surface temperature lead to temperature field differences, for example in regions D where a cooler region exists under the intake valve, and a hotter region E, in the CHT model.

The ensemble averaged near-wall flow fields from the near-wall PIV and both simulations are shown in Figure 17 at −260 CAD aTDCc. Note that the first cells at the fluid-solid interfaces (at y_head = 0 mm) are not available for display from the CHT model due to data export routine limitations. Both models underpredict the velocity magnitude, and differences in the flow direction can be seen. From the RI values, the uniT model captures the near-wall flow direction better than the CHT model. The kinetic energy is underpredicted by both models, with a larger region of small kinetic energy in the CHT model. With increasing distance from the cylinder head, at large y_head values, the influence of the wall model on the flow predictions is reduced, and the direction of the flow in the CHT model aligns better with the measurement.

Figure 18 shows the near-wall flow field at −30 CAD aTDCc. The CHT model provides better predictions in the flow and kinetic energy with increasing wall distance. In contrast, the flow direction of the uniT model is opposite that of the experiment, leading to negative RI values. Figure 19 shows the near-wall flow field at 120 CAD aTDCc. The velocity field directions are well captured, while kinetic energy is overpredicted in both models. Compared to the experimental flow field, both models do not predict the near-wall flow well. This can be explained by the wall model used in the simulations, which may not describe the near-wall flow well, but it also has to be kept in mind that only 10 cycles have been averaged.

The near-wall temperature fields are shown in Figure 20. At −260 CAD aTDCc, the CHT model predicts higher temperatures near the wall compared to the uniT model. In contrast, far away from the wall, the gas temperature is lower in the CHT model. At −30 CAD aTDCc, similar temperatures are obtained near the wall, but with increasing wall distance, the uniT model predicts higher gas temperatures than the CHT model. At 120 CAD aTDCc, again similar temperatures are obtained near the wall, while the gas temperature of the uniT model is lower far away from the wall. These differences arise from the surface temperature, which is used in the heat transfer model for the heat flux calculation at the first cell.

In Figure 21, the spatially averaged, ensemble averaged y ⁺ value is plotted for both simulations in the near-wall PIV FoV. The y ⁺ of 11.05 indicates the switching point in the Werner and Wengle model between viscous and outer layer for the calculation of u ⁺. The CHT model has y ⁺ ≤ 11.05 for most of the engine cycle except for a brief period after EVO. In contrast, the uniT model has y ⁺ > 11.05 after EVO for the duration of the exhaust stroke, for brief periods after IVO during the intake stroke, and late compression and early expansion strokes. No discernable trends are observed in the y ⁺ and predicted velocity vectors in the first two near-wall cells. At −260 CAD aTDCc, both models fall within y ⁺ ≤ 11.05 and the uniT model predicts velocity vector direction better than the CHT model in the first two near-wall cells. At −30 CAD aTDCc, y ⁺ = 11.9 in the uniT model, and y ⁺ = 5.6 in the CHT model. Both models fail to predict the velocity direction in the first two near-wall cells. At 120 CAD, y ⁺ = 6 in the uniT model, while y ⁺ = 1.3 in the CHT model. The velocity direction is well captured, while the velocity magnitude is over-estimated in the first two near-wall cells. On the other hand, from Figure 20, it can be seen that the temperature gradient in the first two near-wall cells is large and this leads to conclusions that the LES may be under-resolved. It may be necessary to increase the grid resolution near the wall to below y ⁺ < 1, which may not be feasible in practical engine simulations.

Another concern is the Han and Reitz heat transfer model, which uses a single equation for T ⁺ for all y ⁺, but the turbulent Prandtl number in their formulation depends on y ⁺ [19]. In fact, the transition point from viscous sublayer to turbulent boundary layer for the calculation of the turbulent Prandtl number was not chosen based on experimental observation but purely mathematical. Therefore, it may be that the turbulent Prandtl number is calculated using the wrong formula or the wrong transition point , which could lead to incorrectly predicted heat transfer values.

The heat flux probe provides surface temperature measurements over the entire engine cycle at the cylinder head. The measured and predicted surface temperatures from the CHT model are plotted in Figure 22a. Note that a low-pass filter of 100 Hz was applied to the measured surface temperature to filter out the noise in the measurement. The shaded area represents one standard deviation of the ensemble average temperature, and indicates CCV of the surface temperature. The CHT model overestimates the surface temperature by 5 K, or an over-estimation of approximately 1.4%. This can be due to approximations in the thermal boundary conditions applied at the outer surfaces of the solid geometries, which were not all measured. In addition, the material properties could be different due to simplifications made to the solid geometries in the CHT simulation. However, the temporal behavior of the surface temperature matches very closely with the measurement during compression. The increase in the surface temperature is well captured, but the drop in the temperature happens earlier in the compression stroke compared to the measurements. The peak temperature occurs at −5 CAD aTDCc in the simulation, while in the experiment the low-pass filtered surface temperature peaks at 14 CAD aTDCc. The unfiltered peak surface temperature occurs at −0.5 CAD aTDCc. The discrepancy in the CAD location of the peak temperature is likely due to the difference in the crank-angle resolution between the simulations (5 CAD) and the measurements (0.5 CAD) in the compression stroke. The standard deviation is also smaller compared to the experiment.

For easier assessment of the temperature prediction, Figure 22b shows the normalized temperature profile . The surface temperature is quasi-steady from −260 CAD to −90 CAD and the mean value from this interval is used to normalize the surface temperature. From that, it is noticed that the surface temperature is underpredicted in the intake, expansion, and exhaust stroke. The normalized temperature reveals that the increase in the surface temperature in the compression stroke is well captured, but the surface temperature decrease happens earlier in the simulation compared to the experiment.

The temperatures at the inner engine surfaces are difficult to measure and thus the following section only presents simulation results from the CHT model that can predict the temporal, spatial, and cyclic variation of the surface temperature. Figure 23 shows the surface temperature of the piston for cycles 2–11 for three different CAD. The piston consists of a metal crown that surrounds a cylindrical quartz window. The thermal conductivity of the quartz material is low compared to the metal crown, which results in a much higher surface temperature in the quartz window compared to the metal crown. The surface temperature is also nonuniform due to the asymmetry of the in-cylinder flow. As the piston comes to TDC, the surface temperature increases due to compression heating of the surrounding gas temperatures. Evident from Figure 23 is the CCV in the surface temperature. These results are in stark contrast with the conventional uniform wall temperature model, where a constant, uniform temperature of 353.16 K was applied to all engine surfaces.

The surface temperatures from the piston, inner liner surface, and cylinder head are first ensemble averaged, then spatially averaged for the entire engine cycle in Figure 24 for each surface. The shaded area indicates one standard deviation of the spatially averaged temperature, i.e., how stratified the temperature distribution is over the surfaces. The temporal variation of the surface temperature in the piston and cylinder liner over the engine cycle is much larger than that of the cylinder head. This is due to the cooling circuit in the cylinder head stabilizing the cylinder head temperature. The cylinder head surface also has the smallest temperature stratification compared to the inner liner and the piston surfaces. The temperature stratification of the inner liner surface is large due to the scrubbing motion of the piston, which exposes the inner liner surfaces to different thermal boundary conditions during the cycle. The piston surface has the largest temperature stratification since it is composed of two different materials. This figure shows that the temperature at each surface is different and has significant spatial and temporal variation, which is in contrast with the uniform temperature assumption. Due to the inherent difficulty in measuring the temperatures at these engine surfaces, the CHT can provide more accurate thermal boundary conditions for engine simulations than only involving the gas-phase.

The surface heat flux measured at the cylinder head is shown in Figure 25 along with the results from the simulations. A 100 Hz low-pass filter was applied to the measured heat flux to filter out the noise in the measurement. Both models predict the correct temporal behavior of the surface heat flux. The peak locations and peak values correspond well with the experiment. The experimental CCV is larger than both simulations, but the CHT model predicts larger CCV than the uniT model with the same number of cycles in the ensemble averaging process. The CHT model also predicts a larger negative heat flux during the expansion stroke.

6 Conclusion

LES of the TCC-III engine was performed using CONVERGE v2.4.18 to study the effects of the thermal boundary condition on the predicted flow and temperature fields, surface temperatures, and surface heat flux. A uniform surface temperature model and a CHT model were utilized for this purpose.

The LES quality index was calculated for both models over the entire engine cycle. It was shown that the LES quality index deteriorates during valve opening events and in the compression stroke mainly due to decreasing integral length scales. The LES quality index is also affected by the number of cycles in the ensemble averaging process, and it can be used as another indicator of statistical convergence of the simulations.

With the CHT model, substantial spatial, temporal, and cycle-to-cycle variability of the wall heat transfer was predicted. The temporal and cyclic behavior of the wall heat transfer compared well to measurements, whereas the spatial variation was not validated due to lack of measured heat transfer data. The surface temperature dynamics of the CHT model at the cylinder head was validated against measurement, with the absolute magnitude about 5 K off. The differences in the predicted temperature fields between the uniform surface temperature and the CHT models arose out of differences in the thermal boundary conditions at the fluid-solid interfaces in each simulation. This leads to differences in the gas density and therefore the predicted flow field. The impact of the surface temperature boundary condition could potentially lead to changes in the temperature and flow fields, e.g., near the spark plug, which will lead to differences in the flame kernel development due to changes in heat transfer demonstrating the need for dynamic modeling of heat transfer in engine simulations.

The effects of super-cycling on LES-CHT simulations have not been evaluated in this study. In addition, the thermal boundary condition from the RANS-CHT could be used in LES to reduce the computational requirements for the more expensive LES-CHT.

Acknowledgments

This material is based upon work supported by the National Science Foundation under the Grant CBET 1402707, “Volumetrically resolved single-shot single-access-point imaging of translucent objects.” Angela Wu acknowledge the support through NSF INTERN program, supplement to CBET 1402707. We greatly appreciate Convergent Science for their technical support and for providing CONVERGE licenses to the University of Michigan.

References

All Tables

All Figures