- Email: [email protected]
Direct numerical simulation of the effect of compression on the flow, temperature and composition under engine-like conditions
Available online at www.sciencedirect.com
Proceedings
ScienceDirect
of the
Combustion Institute
Proceedings of the Combustion Institute xxx (2014) xxx–xxx
www.elsevier.com/locate/proci
Direct numerical simulation of the effect of compression on the flow, temperature and composition under engine-like conditions Martin Schmitt a, Christos E. Frouzakis a,⇑, Ananias G. Tomboulides b, Yuri M. Wright a, Konstantinos Boulouchos a a
Aerothermochemistry and Combustion Systems Laboratory, Swiss Federal Institute of Technology, Zurich, Switzerland b Department of Mechanical Engineering, University of Western Macedonia, Kozani, Greece
Abstract The effect of compression on the flow, temperature and composition inside a cylinder is investigated using direct numerical simulation (DNS). The initial conditions are obtained by a separate DNS of the intake stroke in an open-valve setup which includes thermal and species mixing. The results show significant changes of the turbulence and temperature fields during compression: The decrease of kinematic viscosity resulting from the increasing pressure results in smaller turbulent length scales and higher dissipation rates. Temperature fluctuations away from the walls decrease slightly during the first half but increase strongly during the second half of the compression stroke towards the Top Dead Center (TDC) due to heat transfer to and from the walls and turbulent transport. At TDC the turbulent flow field is anisotropic, and the axial fluctuation velocity is approximately 30% smaller than the fluctuation velocities in the radial and azimuthal directions. The integral length scale of temperature is approximately 25% higher than the integral length scale of turbulent kinetic energy. The stratification in the species concentration is found to be practically negligible. Ó 2014 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Direct numerical simulation; Internal combustion engine; HCCI; Compression; Thermal stratification
1. Introduction Homogeneous charge-compression ignition (HCCI) engine concepts offer the potential of increased efficiencies with very low NOx and soot emissions [1]. The higher efficiency can be ⇑ Corresponding author. Fax: +41 44 632 12 55.
E-mail address: (C.E. Frouzakis).
[email protected]
achieved by higher compression ratios, absence of throttling losses and thermodynamically-favorable nearly-isochoric combustion [2]. However, the rapid heat release also results in increased stresses on the engine structure and noisy operation [3]. Reduced NOx and soot emissions result from the relatively low peak combustion temperature attained by lean mixtures or high External Gas Recirculation (EGR) rates and well-mixed operation. In contrast to soot and NOx, higher CO and unburnt hydrocarbon emissions are due
http://dx.doi.org/10.1016/j.proci.2014.06.097 1540-7489/Ó 2014 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
2
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
to temperature and/or mixture inhomogeneities close to the cylinder walls [4]. Inhomogeneities have a strong influence on the temporal development of the heat release, since hotter and fuel rich regions tend to ignite before zones with less favorable ignition conditions. Temperature stratification has a strong influence on the autoignition behavior through the exponential dependence of the reaction rate on temperature, and can be beneficial for reducing the knocking tendency by lowering the rate or pressure rise. Sjo¨berg et al. [5] assessed experimentally the effect of increased thermal stratification on HCCI combustion and found it to be beneficial in significantly extending the high-load operating point. Recently, several studies employed 2-D DNS to describe the influence of weakly thermally stratified initial conditions on the HCCI autoignition chemistry (e.g. [2,6–10]). In Ref. [2,6], the autoignition behavior was investigated for different artificial turbulent and temperature initial conditions of a homogeneous H2/air mixture at 41 atm. The fluctuations of a zero-mean velocity and constant-mean temperature fields were initialized using a Passot–Pouquet spectrum [11], which requires the most energetic wavenumber and the fluctuation velocities as input parameters. Temperature fluctuations were found to have a strong effect on the combustion mode, the timing and the duration of the heat release rate leading to advanced ignition, longer combustion and to a significantly higher fuel consumption rate due to deflagration instead of homogeneous ignition. Under real engine conditions, turbulent quantities like the integral length scales or fluctuation velocities are difficult to measure or estimate and are expected to vary significantly within the domain. Insight into the flow and temperature fields in engines before ignition can be obtained using high-speed Particle Image Velocimetry (PIV) and/or Laser-Induced Fluorescence (LIF) measurements. One of the first applications of 3D PIV measurements in engines [12] recently investigated the flow field in a cube of dimensions 47 35 4 mm3, with a spatial resolution of 0.4 mm in each direction. Recent LIF measurements are reported in Refs. [13–15]. In these works the temperature distributions were obtained on vertical slices with a spatial resolution of approximately 0.4 mm. In contrast to [13], the results of [14,15] showed increasing thermal stratification towards Top Dead Center (TDC). The PIV and LIF measurements were performed on planes or small volumes within the domain, which provide a valuable but limited insight on the flow and temperature fields inside the cylinder. In addition, measurements close to the walls are challenging and it is difficult to resolve the thin transient boundary layer which has a strong influence on the thermal stratification during compression.
DNS offers the advantage of high spatial and temporal resolution. However, due to the high computational cost, few simulations of complex geometries with moving boundaries can be found in the literature. Gu¨ntsch [16] studied the effect of compression on artificial turbulence in 3D cylindrical geometries with a compression ratio of 4.2. The effect of different piston speed profiles and the influence of swirl were investigated by statistical analysis (mean and fluctuating quantities, PDFs, two point correlations, turbulent length scales, statistical moments of third and fourth order and energy density spectra). In Ref. [25], the incompressible flow and cyclic variability in the valve/piston assembly studied experimentally by Morse et al. [24] was investigated using DNS. Starting from the numerical data of the latter, a precursor simulation was performed to compute the mixing of a cold mixture with the hot gases inside the cylinder and obtain engine-relevant initial conditions for the work presented in this paper. Despite its simplicity, the setup considered here preserves many features of the dynamics of real engine flows (e.g. jet breakup during the intake stroke, interaction of the jet flow and shear layers with the walls and the remaining turbulence at TDC of the previous cycle). In addition, the axisymmetric geometry allows for spatial averaging, which heavily reduces the simulation time for comparisons with RANS or LES simulations. The relatively low engine speed of 560 rpm which is dictated by the computational cost is representative of idle conditions in a variety of ICEs. The aim of this work is to obtain 3-D spatio-temporal data at a resolution that cannot be provided even by sophisticated experiments in the cylinder of a laboratory-scale setup and to study the effect of compression on the flow, temperature and species distributions. 2. Computational method The low Mach number form of the conservation equation are integrated in time with a highly-efficient parallel code based on nek5000 [17]. The solution is expressed in terms of tensor products of nth-order Lagrange polynomials based on the Gauss–Lobbato–Legendre quadrature points [18] within hexahedral conforming elements. Time integration employs a high-order time splitting scheme for low Mach number reactive flows [19], whereby the continuity and momentum equations are integrated with a semiimplicit scheme and the species and energy equations are integrated implicitly using CVODE [20]. The Arbitrary Lagrangian/Eulerian (ALE) approach is employed to account for the mesh variation resulting from the piston movement. The ALE implementation described in Ref. [21]
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
was extended in this work to the so-called P n P n formulation, where pressure and velocity are solved on the same mesh. An additional term, to account for the temporally-varying thermodynamic pressure (dp=dt) in the closed domain was also introduced in the energy equation. A mixture-average transport model is adopted for the species diffusion velocity and thermodynamic/ transport properties are evaluated using CHEMKIN libraries [22,23], while reactivity is suppressed. The simulations were performed on a Cray XE6 and required approximately 353,000 CPUh for the simulation of one compression cycle using polynomial orders higher than 9 and a 3rd-order temporal integration. 3. Simulation setup This work is based on a previous DNS study [25], which investigated the incompressible flow past the open valve induced by the moving piston in the geometry experimentally studied by Morse et al. [24] and shown schematically in Fig. 1(a). Eight consecutive cycles were simulated showing significant cyclic variations, which were found to result from the instability of the hollow jet created at the valve caused by the interaction of the jet with the turbulence remaining from the previous cycle and by the orientation and stability of a large central vortex ring formed under the valve. A detailed description of the grid, the resolution requirements and the flow dynamics can be found in Ref. [25]. In this work, a constant temperature of 500 K is imposed on all boundaries while for the species zero flux boundary conditions are specified.
3
3.1. Geometry and computational grid For the calculation of the compression stroke the valve is closed at BDC, i.e. at the time instant of zero mass flux through the valve. The computational domain reduces to only the cylinder part of Fig. 1(a), and the initial conditions are obtained from the precursor mixing simulation at BDC discussed in the next subsection. At the chosen compression ratio of 12, the piston moves from 90 mm (BDC) to 7.5 mm (TDC). Two grids were employed constructed by extruding horizontal slices in the axial direction. The first grid containing 2436 spectral elements on the x y plane (Fig. 1) and 60 layers of elements in the axial direction clustered towards the piston and cylinder head is used for crank angles between 180° and 306°. For a polynomial order of n ¼ 9, an average axial resolution of 180 lm at BDC and 50 lm at 306° is obtained with a total of approximately 90 million nodes. In the radial and azimuthal directions the resolution is approximately 150 lm, while refinement close to the walls results in a resolution of better than 13.8 lm. In order to satisfy the increased resolution requirements during the last part of the compression stroke (>306°CA), a second grid with 3410 elements on the horizontal slice and 30 layers of elements in the axial direction was used and the fields were interpolated using high order spectral interpolation. By increasing the polynomial order to n ¼ 11 (corresponding to 135 million nodes) the average axial resolution is 110 lm at 306° and 30 lm at TDC. In the radial and azimuthal directions the resolution is 85 lm in the bulk of the domain
y
r φ
x
90
z
(a)
(b)
Fig. 1. (a) Schematic view of the domain of the open valve setup experimentally studied in Ref. [24], (b) axial slice of the second mesh used in the compression simulation.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
4
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
and better than 9.75 lm close to the walls. The pffiffiffiffiffiffiffiffiffiffi corresponding y þ ¼ ym sw =q value for the first mesh point at the wall is below 0.35 until 306°CA and below 0.9 until TDC, where y is the dimensional distance from the wall, sw the wall shear stress, and m; q the mixture viscosity and density. The resolution requirements from BDC up to 310°CA are comparable to those of the incompressible flow simulation [25] where they were assessed by a mesh resolution study and by comparison with the experimental data of [24]. As described below, with increasing pressure the size of the smallest scales decrease, resulting in significantly increased resolution requirements close to TDC. The Kolmogorov length scale gK ¼ 1=4 ðm3=Þ at TDC is gK ¼ 12:8 lm, where ¼ 2m sij sij is the mean dissipation rate computed from the strain tensor sij ¼ @ui =@xj þ @uj =@xi . At TDC, the radial and azimuthal resolution is approximately 6:6gK , while in the axial direction dz 2:3gK , which for the employed high order spectral element discretization method is enough to capture accurately most of the dissipation, and thus provide reliable first and second order statistics as pointed out in Ref. [26]. The resolution is also fine enough to resolve the mixing scales (i.e. 1=4 the Batchelor scale gB ¼ ðmD2 =Þ and the Obuk1=4 3 hov–Corrsin scale gOC ¼ ðD =Þ ), which are larger than gK since the Schmidt number is Sc ¼ m=D < 1. 3.2. Initial conditions In order to obtain realistic initial conditions for the temperature and species fields at the BDC of the compression stroke, the intake stroke was simulated as follows: starting with the velocity field at TDC of the third cycle of the cold flow simulation (so that the effect of the quiescent initial field is minimized) of the valve/cylinder assembly [25], the mixture inside the cylinder was assumed to be a homogeneous mixture with Y O2 ¼ 0:115; Y N2 ¼ 0:756 and Y H2 O ¼ 0:129 at 900 K and atmospheric pressure; this mixture corresponds to complete combustion of a / ¼ 0:5 H2/ air mixture and emulates EGR conditions. At the inflow, a homogeneous / ¼ 0:5 H2/air mixture at 500 K was introduced and a constant temperature of 500 K was imposed on all boundaries; for the species, zero flux boundary conditions were specified. The piston speed was 560 rpm and the Reynolds number based on the maximum mean jet speed, the cylinder radius and the viscosity of the intake mixture (m = 4.52 105 m2/s) is Re = 13,643. The downward moving piston draws the colder H2/air mixture into the cylinder. During the intake stroke, a large central vortex ring is formed below the open valve (indicated at Bottom
Dead Center (BDC) by a pressure isosurface in the velocity magnitude plot of Fig. 2(a)), which entrains hot EGR gases, while regions of lower temperature form close to the cylinder liner and the piston head. The initial conditions for temperature and species were found to have a weak effect on the conditions at the end of the compression stroke. As described below, at TDC the species are nearly homogeneously distributed whereas the temperature is primarily affected by heat transfer to and from the walls and not by its initial distribution. The mean velocity field is obtained by averaging in the azimuthal direction. For temperature, statistics are also collected within an inner cylinder of radius 30 mm and axial extend (from the cylinder to the piston) of 7.5 mm until 270°CA and 2 mm after 270°CA. This was done because averages in the whole domain are strongly influenced by the cooler near wall region. In the following, the volume average of any variable w in the whole cylinder will be denoted as hwiV and in the inner cylinder as hwiV ;i , whereas its magnitude is wmag ¼ jwj. The integral length scale Lw is derived by first extracting w along a circle at a fixed radial and axial position and then by integrating its autocorrelation to the point where 10% of its maximum value q is ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi reached [30]. The velocity magnitude is jvj ¼
v2z þ v2r þ v2u , while
the turbulent kinetic energy is defined as 02 02 with v0r ; v0u and v0z being the þ v þ v k ¼ 12 v02 r u z fluctuation velocities and T 0 the fluctuation temperature. 4. Results 4.1. Compression stroke The instantaneous distributions of velocity magnitude, temperature and Y H2 O at BDC are plotted in Fig. 2(a). The dominant flow feature is a large central vortex ring as illustrated in the figure by the pressure isosurface at p = 3.55 Pa; the white arrow indicates the clockwise flow direction around the vortex ring. Entrained inside the vortex ring core is fluid with higher temperature and water vapor mass fraction. In contrast, Fig. 2(b) shows the same fields at the end of compression (TDC) on an axial slice at 3.75 mm below the cylinder head. The temperature field shows high fluctuations in the whole domain with steep gradients close to the cylinder wall. In wall normal units the temperature gradient in the viscous sublayer (y þ < 5) follows the relation T þ ¼ Pr y þ , as reported in Ref. [27] for non-isothermal channel flows. In the inner cylinder region temperatures below 1050 K and above
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
| | [m/s]
| | [m/s]
9.8
3.36
0
0
T [K]
T [K]
600
1225
500
900
Y H2O
Y H2O
0.045
0.0325
0.015
0.0304
(a)
5
(b)
Fig. 2. Velocity magnitude, temperature and water mass fraction on (a) a vertical slice through the axis at BDC, and (b) on a axial slice 3.75 mm below the cylinder head at TDC.
30
2000
Isentropic compression DNS
20
1500
15 10
1000
Temperature [K]
25
pressure [atm]
1200 K can be observed, while the average temperature is close to 1110 K. Water vapor is practically homogeneously distributed inside the cylinder, as indicated by the narrow range of the color scale. The volume averaged pressure and temperature variation in time is shown in Fig. 3 (solid lines) and compared with the values computed for an isentropic compression using variable specific heat capacity ratio c. From BDC to TDC the mean temperature and pressure increase by a factor of 1.93 and 23.09, respectively. This results in a reduction of the kinematic viscosity by a factor of approximately 6.5 which strongly affects the evolution of turbulence. The differences between DNS and the isentropic curves are very small in the first half and increase significantly in the second half of the compression, reaching 3.9 atm and 148 K at TDC. This increase is due to the strongly increased heat transfer towards TDC, related to the decreasing thermal boundary layer thickness and the increasing temperature during compression. The pressure peaks 1.54°CA before TDC (thermodynamic loss angle), which is in
5 180
225
270
315
500 360
piston position [°CA] Fig. 3. Comparison of in-cylinder pressure and mean temperature with the values computed for isentropic compression.
the range of values reported in Ref. [28]. The relative high heat losses are related to the low engine speed (560 rpm) and the relatively cool wall temperatures compared to the average mixture temperature at BDC.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
6
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx 0
The variations of k, and T in the whole and the inner cylinder during the compression stroke are reported in Table 1. The decreasing kinematic viscosity results in higher velocity gradients around the smaller flow structures, and thus increased dissipation during the compression stroke. The peak is reached at 346° CA, i.e. before TDC, since close to TDC the turbulent kinetic energy is strongly reduced by the high dissipation level and zero piston speed. As shown in Ref. [2], temperature fluctuations can significantly affect the timing of autoignition in HCCI combustion and the corresponding peak of the pressure curve. Temperature fluctuations drop from BDC to 225°CA but increase significantly after 270°CA, mainly as a result of
the increasing difference between the average temperature in the cylinder and the fixed wall temperatures; this behavior is in agreement with the LIF measurements reported in Refs. [14,15]. Since autoignition is more likely inside the hotter inner cylinder zone, a more useful indicator for the description of the thermal stratification is the average fluctuations in the inner temperature 0:5 cylinder T 02 V ;i . This quantity decreases from 12.20 to 6.81 K in the first half of the compression stroke, possibly due to the homogenization of the temperature field in the inner cylinder. This trend reverses 0:5 in the second half of the stroke and T 02 V ;i increases significantly. As TDC is approached, the steepening of the velocity and temperature gradients near the walls results in increased heat transfer. Moreover localized ejections transporting cold fluid from the walls towards the cylinder center were found to be the main mechanism for increasing the rate of convective heat transfer during compression. In Fig. 4, the probability density functions (PDFs) of the turbulent kinetic energy and of temperature within the whole cylinder are plotted for several piston positions during the compression stroke. At 180°CA, PDF(k) shows a relatively flat profile with a peak at approximately 0.73 m2/s2. During compression the peak becomes more pronounced, but the position of the maximum remains at values close to 0.7 m2/s2. The relatively
Table 1 Variation of the volume averaged k, and the rms T during compression (in m2 =s2 ; m2 =s3 ; K). 02 0:5 02 0:5 °CA hk iV hiV T V T V ;i 180 225 270 306 330 346 360
5.70 3.78 2.50 1.83 1.79 1.58 0.98
820.00 536.50 973.52 3502.73 19281.7 36641.3 25510.0
16.80 14.84 29.69 58.04 92.99 110.7 119.9
12.20 6.73 6.81 13.41 32.29 39.15 47.31
0.8 0.7
PDF(k)
(a)
180°CA 225°CA 270°CA 306°CA 330°CA 346°CA 360°CA
0.6 0.5 0.4 0.3 0.2 0.1 0
0
1
2
3
4
6
5
7
8
2 2
k [m /s ] 0.06
(b)
PDF(T)
0.05 0.04 0.03 0.02 0.01 0 500
600
700
800
900
1000
1100
1200
T [K] Fig. 4. Probability density functions of (a) turbulent kinetic energy, and (b) temperature at several piston positions during compression.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
large difference in the peak values at 346° and 360°CA can be attributed to the increased dissipation rate towards TDC. According to [9] the direct effect of velocity fluctuations on autoignition is small compared to the effect of thermal stratification. Hence, turbulence is indirectly influencing autoignition by enhancing thermal stratification due to convective heat transport close to the walls. The influence of the initial conditions on temperature can be seen by comparing the PDF ðT Þ at 180° and 225°CA. The minimum and maximum temperatures are similar, but the peak becomes more pronounced and shifts to a higher value, as a result of the temperature homogenization during the early compression stage discussed above. The higher peaks at 225° and 270°CA are due to the low temperature fluctuations in the inner cylinder (Table 1). The PDF shape differs significantly from the almost Gaussian distribution reported in the experimental study of [15], possibly because of the difficulty in resolving the thermal boundary layers in the two-dimensional measurements. Between 306° and 330°CA, due to the effect of wall heat transfer, the maximum PDF ðT Þ values decrease strongly and close to TDC (346° and 360°CA), the PDF ðT Þ becomes almost flat. This can be explained by the high wall heat losses during this period. The maximum temperature 1234.5 K is reached at 360°CA although the average temperature is nearly the same as at 346°CA (Table 1). The temporal development of PDF ðT Þ indicates that the temperature initial conditions are washed out after the first half of the compression stroke, implying that, at least for this setup, the thermal stratification at TDC cannot be controlled by changing the temperature field at BDC. Stronger effects can be expected by modifying the wall temperature, as also found in Ref. [29], where it was shown that lower coolant temperature had a strong effect on the heat release rate for HCCI operation.
7
4.2. Behavior at Top Dead Center The distribution of v0z on an axial slice 3.75 mm below the cylinder head is shown in Fig. 5(a), indicating an inhomogeneous distribution within the domain. The average values and the PDFs of the three fluctuation velocity components are plotted in Fig. 5(b), where their volume averages values are also reported. The PDFs of the radial and azimuthal fluctuation velocities as well as their average values in the radial and azimuthal directions are nearlyDidentical, E while the average in the axial 0 direction vz;mag is approximately 25% lower. V
Accordingly, the peak PDF value is higher, possibly due to the dumping of axial fluctuations in the narrow gap between the piston and the cylinder head. The spatial distribution of the integral length scales of k and temperature are plotted in Fig. 6. Values close to the cylinder center are not shown, since the number of statistically independent points is too low to obtain meaningful values. With the exception of two zones close to the cylinder head, Lk lies in the range between 1 and 2.5 mm, while LT varies between 3 and 4 mm in the cylinder center and decreases towards the walls due to the increased temperature fluctuations in the thermal boundary layer. The average length scale of k is in the whole domain approximately 25% smaller compared to the average scale of the temperature. The PDFs of Lk and LT plotted in Fig. 7 have a similar shape and their peaks are located at approximately 1.1 mm. In contrast, the PDF ðgk Þ is more symmetric and gK ranges from 0.005 to 0.025 with an average value of 12.8 lm. The smaller gk values are mainly located in the relatively cooler thermal boundary layers since m / T 1:5 . The ratio between hLk iV and hgk iV at TDC is approximately 161.4. The low value of hgk iV close to TDC is mainly due to the increase in pressure followed by a reduced kinematic viscosity and a
1.5 1.2
2.8
PDF
z
[m/s]
0.9
<
z,mag
<
r,mag
>V = 0.518 m/s
<
φ,mag
>V = 0.691m/s >V = 0.671 m/s
0.6 0.3
0
0
1
2
3
fluctuation velocities [m/s]
(a)
(b)
Fig. 5. (a) Axial fluctuation velocity on a horizontal slice 3.75 mm below the cylinder head; (b) PDF of the axial, radial and azimuthal fluctuations velocity components at TDC.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
8
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
Lk
= 2.07 mm 5
cylinder wall
cylinder axis
cylinder head
piston
LT
4 3
= 2.57 mm
cylinder axis
z=0 2 1
r=0
r = 30
r = 15
z = -7.5 r = 37.5
Fig. 6. Integral length scales in mm of the turbulent kinetic energy and the temperature at TDC.
LT Lk ηK 100 < ηK>V = 12.8 [μm]
PDF (L T, k )
0.4 0.3 0.2
50
0.1 0
2
4
6
PDF (ηK )
150 0.5
0
8
L [mm] 0
0.005
0.01
0.015
0.02
0.025
ηK [mm]
Fig. 7. PDF of the integral length scales of turbulent kinetic energy and temperature and the Kolmogorov scale at TDC. 15000
YH2
V = 0.01084
PDF
YH2
YO2
10000
0.0109
V = 0.2018
YH2O
5000V = 0.0315 0.0107
-0.001
0
averages. The higher maximum of PDF ðY H2 Þ can be explained by the high diffusivity of hydrogen, which has a Lewis number approximately three times smaller than the other two species. The PDFs of Y O2 and Y H2 O show approximately equal deviations and the local Y H2 O maxima correlate with the local Y O2 minima. However the deviations in the species mass fractions are very small compared to those of temperature, in agreement with the results in Ref. [14].
0.001
deviation from < Y >V
Fig. 8. PDFs of H2, O2 and H2O mass fractions at TDC. The zero value on the x-axis corresponds to the average species mass fraction. The horizontal slice are located 3.75 mm below the cylinder head.
higher dissipation rate. On the other hand, the influence of the increasing temperature on hgk iV is less pronounced, since temperature is less sensitive to volume changes compared to pressure. Thus, the compression ratio has a strong impact on the size distribution of the turbulent scales, which in the case of a propagating flame front (as in an Otto engine) is expected to have a significant impact on the combustion process and the engine efficiency. The PDFs of Y H2 ; Y O2 and Y H2 O reported in Fig. 8 show the deviations from the volumetric
5. Conclusions Direct numerical simulation was performed to study the effect of compression on the flow, temperature and composition fields in the closed cylinder of the valve-piston assembly investigated experimentally [24] and numerically [25]. The results show significant changes of the flow and temperature during compression: The reduction of the kinematic viscosity due to the increasing pressure leads to smaller turbulent length scales and increased dissipation rates close to TDC. Temperature fluctuations decrease during the first half of the compression stroke, and, in agreement with the experimental findings of [14,15], increase significantly during the second half, mainly as a result of the heat transfer to the wall and turbulent convective transport of cooler gas from the walls towards the inner part of the cylinder. The influence of the temperature initial conditions at BDC on the temperature fluctuations at TDC is found to be very small, while, in agreement with [14], the species distributions at TDC are almost homogeneous due to the very efficient mixing during compression. The flow field at TDC is anisotropic, and the average fluctuation velocities in the axial direction are approximately 30% smaller compared to those in the radial and azimuthal direction, possibly due to the smaller axial dimension at TDC. Due to thermal diffusion, the integral length scale hLT i is approximately 25% higher compared to hLk i. The average Kolmogorov scale is 12.8 lm and around 160 times smaller than hLk i. While only one cycle is investigated here, ongoing simulations of the compression stroke using the flow fields from different cycles of the multicycle simulation in Ref. [25] showed that although large flow structures have an influence on the flow and temperature field, the observed trends for the evolution of mean and rms temperatures, k, and the species mixing are very similar for all calculated cycles. We are currently analyzing wall heat transfer in this setup and the phenomena associated with velocity and thermal boundary layers. Future work will consider reactive conditions and investigate the autoignition behavior of the mixture.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
Acknowledgments Financial support of the Swiss National Science Foundation (SNF) under project number 200021_135514 is gratefully acknowledged. The work was also supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID 189. References [1] J.E. Dec, Proc. Combust. Inst. 32 (2009) 2727–2742. [2] E.R. Hawkes, R. Sankaran, H.G. Im, P.P. Pe´bay, J.H. Chen, Combust. Flame 30 (2005) 875–882. [3] M. Yao, Z. Zheng, H. Liu, Prog. Energy Combust. Sci. 35 (2009) 398–437. [4] F. Zhao, T. Asmus, D. Assanis, J. Dec, J. Eng, P. Najt, Homogeneous Charge Compression Ignition (HCCI) Engines, Key Research and Development Issues, SAE Inc., Warrendale, PA, 2002, 1–14 and 531–541. [5] M. Sjo¨berg, J.E. Dec, N.P. Cernansky, SAE Technical paper, 2005. 2005-01-0113. [6] J.H. Chen, E.R. Hawkes, R. Sankaran, S.D. Mason, H.G. Im, Combust. Flame 145 (2006) 128–144. [7] R. Sankaran, H.G. Im, E.R. Hawkes, J.H. Chen, Proc. Combust. Inst. 145 (2006) 145–159. [8] C.S. Yoo, T. Lu, J.H. Chen, C.K. Law, Combust. Flame 158 (2011) 1727–1741. [9] M.B. Luong, Z. Luo, T. Lu, S.H. Chung, C.S. Yoo, Combust. Flame 160 (2013) 2038–2047. [10] G. Bansal, H.G. Im, S.R. Lee, Proc. Combust. Inst. 32 (2009) 1083–1090. [11] J.O. Hinze, Turbulence, McGraw-Hill, New York, 1975. [12] E. Baum, B. Peterson, C. Surmann, D. Michaelis, B. Bo¨hm, A. Dreizler, Proc. Combust. Inst. 34 (2013) 2903–2910. [13] B. Peterson, E. Baum, B. Bo¨hm, V. Sick, A. Dreizler, Proc. Combust. Inst. 34 (2013) 3653–3660. [14] J. Snyder, N. Dronniou, J. Dec, R. Hanson, SAE Technical paper, 2011. -01-1291.
9
[15] S.A. Kaiser, M. Schild, C. Schulz, Proc. Combust. Inst. 34 (2013) 2911–2919. [16] E. Gu¨ntsch, R. Friedrich, Flow Simulation with High-Performance Computers II, Springer, 1996, pp. 213–226. [17] P.F. Fischer, J.W. Lottes, S.G. Kerkemeier, http:// nek5000.mcs.anl.gov, (2008). [18] M.O. Deville, P.F. Fischer, E.H. Mund, Cambridge U. Press, New York, 2002. [19] A.G. Tomboulides, J. Lee, S.A. Orszag, J. Sci. Comput. 12 (1997) 139–167. [20] G.D. Byrne, A.C. Hindmarsh, Int. J. High Perform. Comput. Appl. 13 (1999) 354–365. [21] L.W. Ho, A.T. Patera, Comput. Methods Appl. Mech. Eng. 80 (1990) 355–366. [22] R.J. Kee, F.M. Rupley, J.A. Miller, Chemkin II: A Fortran Chemical Kinetics Package for the Analysis of Gas-phase Chemical Kinetics, Report No. SAND89-8009B, Sandia National Laboratories, 1996. [23] R.J. Kee, G. Dixon-Lewis, J. Warnatz. M.E. Coltrin, J.A. Miller, A Fortran Computer Code Package for the Evaluation of Gas-Phase Multicomponent Transport Properties, Report No. SAND868246, Sandia National Laboratories, 1996. [24] A.P. Morse, J.H. Whitelaw, M. Yianneskis, J. Fluids Eng. 102 (1979) 208–2016. [25] M. Schmitt, C.E. Frouzakis, A.G. Tomboulides, Y.M. Wright, K. Boulouchos, Phys. Fluids 26 (2014) 035105. [26] P. Moin, K. Mahesh, Ann. Rev. Fluid Mech. 30 (1998) 539–578. [27] H. Kawamura, K. Ohsaka, H. Abe, K. Yamamoto, Int. J. Heat Fluid Flow 19 (5) (1998) 482– 491. [28] R. Pischinger, M. Klell, T. Sams, Thermodynamik der Verbrennungskraftmaschine: Der Fahrzeugantrieb, Springer-Verlag, Vienna, 2010. [29] M. Sjo¨berg, J.E. Dec, A. Babajimopoulos, D. Assanis, SAE paper, 2004-01-2994. [30] W.J. Feiereisen, W.C. Reynolds, J.H. Ferziger, Numerical Simulation of a Compressible Homogeneous Turbulent Shear Flow, Tech. Rep., Stanford University, 1989.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
Proceedings
ScienceDirect
of the
Combustion Institute
Proceedings of the Combustion Institute xxx (2014) xxx–xxx
www.elsevier.com/locate/proci
Direct numerical simulation of the effect of compression on the flow, temperature and composition under engine-like conditions Martin Schmitt a, Christos E. Frouzakis a,⇑, Ananias G. Tomboulides b, Yuri M. Wright a, Konstantinos Boulouchos a a
Aerothermochemistry and Combustion Systems Laboratory, Swiss Federal Institute of Technology, Zurich, Switzerland b Department of Mechanical Engineering, University of Western Macedonia, Kozani, Greece
Abstract The effect of compression on the flow, temperature and composition inside a cylinder is investigated using direct numerical simulation (DNS). The initial conditions are obtained by a separate DNS of the intake stroke in an open-valve setup which includes thermal and species mixing. The results show significant changes of the turbulence and temperature fields during compression: The decrease of kinematic viscosity resulting from the increasing pressure results in smaller turbulent length scales and higher dissipation rates. Temperature fluctuations away from the walls decrease slightly during the first half but increase strongly during the second half of the compression stroke towards the Top Dead Center (TDC) due to heat transfer to and from the walls and turbulent transport. At TDC the turbulent flow field is anisotropic, and the axial fluctuation velocity is approximately 30% smaller than the fluctuation velocities in the radial and azimuthal directions. The integral length scale of temperature is approximately 25% higher than the integral length scale of turbulent kinetic energy. The stratification in the species concentration is found to be practically negligible. Ó 2014 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Direct numerical simulation; Internal combustion engine; HCCI; Compression; Thermal stratification
1. Introduction Homogeneous charge-compression ignition (HCCI) engine concepts offer the potential of increased efficiencies with very low NOx and soot emissions [1]. The higher efficiency can be ⇑ Corresponding author. Fax: +41 44 632 12 55.
E-mail address: (C.E. Frouzakis).
[email protected]
achieved by higher compression ratios, absence of throttling losses and thermodynamically-favorable nearly-isochoric combustion [2]. However, the rapid heat release also results in increased stresses on the engine structure and noisy operation [3]. Reduced NOx and soot emissions result from the relatively low peak combustion temperature attained by lean mixtures or high External Gas Recirculation (EGR) rates and well-mixed operation. In contrast to soot and NOx, higher CO and unburnt hydrocarbon emissions are due
http://dx.doi.org/10.1016/j.proci.2014.06.097 1540-7489/Ó 2014 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
2
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
to temperature and/or mixture inhomogeneities close to the cylinder walls [4]. Inhomogeneities have a strong influence on the temporal development of the heat release, since hotter and fuel rich regions tend to ignite before zones with less favorable ignition conditions. Temperature stratification has a strong influence on the autoignition behavior through the exponential dependence of the reaction rate on temperature, and can be beneficial for reducing the knocking tendency by lowering the rate or pressure rise. Sjo¨berg et al. [5] assessed experimentally the effect of increased thermal stratification on HCCI combustion and found it to be beneficial in significantly extending the high-load operating point. Recently, several studies employed 2-D DNS to describe the influence of weakly thermally stratified initial conditions on the HCCI autoignition chemistry (e.g. [2,6–10]). In Ref. [2,6], the autoignition behavior was investigated for different artificial turbulent and temperature initial conditions of a homogeneous H2/air mixture at 41 atm. The fluctuations of a zero-mean velocity and constant-mean temperature fields were initialized using a Passot–Pouquet spectrum [11], which requires the most energetic wavenumber and the fluctuation velocities as input parameters. Temperature fluctuations were found to have a strong effect on the combustion mode, the timing and the duration of the heat release rate leading to advanced ignition, longer combustion and to a significantly higher fuel consumption rate due to deflagration instead of homogeneous ignition. Under real engine conditions, turbulent quantities like the integral length scales or fluctuation velocities are difficult to measure or estimate and are expected to vary significantly within the domain. Insight into the flow and temperature fields in engines before ignition can be obtained using high-speed Particle Image Velocimetry (PIV) and/or Laser-Induced Fluorescence (LIF) measurements. One of the first applications of 3D PIV measurements in engines [12] recently investigated the flow field in a cube of dimensions 47 35 4 mm3, with a spatial resolution of 0.4 mm in each direction. Recent LIF measurements are reported in Refs. [13–15]. In these works the temperature distributions were obtained on vertical slices with a spatial resolution of approximately 0.4 mm. In contrast to [13], the results of [14,15] showed increasing thermal stratification towards Top Dead Center (TDC). The PIV and LIF measurements were performed on planes or small volumes within the domain, which provide a valuable but limited insight on the flow and temperature fields inside the cylinder. In addition, measurements close to the walls are challenging and it is difficult to resolve the thin transient boundary layer which has a strong influence on the thermal stratification during compression.
DNS offers the advantage of high spatial and temporal resolution. However, due to the high computational cost, few simulations of complex geometries with moving boundaries can be found in the literature. Gu¨ntsch [16] studied the effect of compression on artificial turbulence in 3D cylindrical geometries with a compression ratio of 4.2. The effect of different piston speed profiles and the influence of swirl were investigated by statistical analysis (mean and fluctuating quantities, PDFs, two point correlations, turbulent length scales, statistical moments of third and fourth order and energy density spectra). In Ref. [25], the incompressible flow and cyclic variability in the valve/piston assembly studied experimentally by Morse et al. [24] was investigated using DNS. Starting from the numerical data of the latter, a precursor simulation was performed to compute the mixing of a cold mixture with the hot gases inside the cylinder and obtain engine-relevant initial conditions for the work presented in this paper. Despite its simplicity, the setup considered here preserves many features of the dynamics of real engine flows (e.g. jet breakup during the intake stroke, interaction of the jet flow and shear layers with the walls and the remaining turbulence at TDC of the previous cycle). In addition, the axisymmetric geometry allows for spatial averaging, which heavily reduces the simulation time for comparisons with RANS or LES simulations. The relatively low engine speed of 560 rpm which is dictated by the computational cost is representative of idle conditions in a variety of ICEs. The aim of this work is to obtain 3-D spatio-temporal data at a resolution that cannot be provided even by sophisticated experiments in the cylinder of a laboratory-scale setup and to study the effect of compression on the flow, temperature and species distributions. 2. Computational method The low Mach number form of the conservation equation are integrated in time with a highly-efficient parallel code based on nek5000 [17]. The solution is expressed in terms of tensor products of nth-order Lagrange polynomials based on the Gauss–Lobbato–Legendre quadrature points [18] within hexahedral conforming elements. Time integration employs a high-order time splitting scheme for low Mach number reactive flows [19], whereby the continuity and momentum equations are integrated with a semiimplicit scheme and the species and energy equations are integrated implicitly using CVODE [20]. The Arbitrary Lagrangian/Eulerian (ALE) approach is employed to account for the mesh variation resulting from the piston movement. The ALE implementation described in Ref. [21]
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
was extended in this work to the so-called P n P n formulation, where pressure and velocity are solved on the same mesh. An additional term, to account for the temporally-varying thermodynamic pressure (dp=dt) in the closed domain was also introduced in the energy equation. A mixture-average transport model is adopted for the species diffusion velocity and thermodynamic/ transport properties are evaluated using CHEMKIN libraries [22,23], while reactivity is suppressed. The simulations were performed on a Cray XE6 and required approximately 353,000 CPUh for the simulation of one compression cycle using polynomial orders higher than 9 and a 3rd-order temporal integration. 3. Simulation setup This work is based on a previous DNS study [25], which investigated the incompressible flow past the open valve induced by the moving piston in the geometry experimentally studied by Morse et al. [24] and shown schematically in Fig. 1(a). Eight consecutive cycles were simulated showing significant cyclic variations, which were found to result from the instability of the hollow jet created at the valve caused by the interaction of the jet with the turbulence remaining from the previous cycle and by the orientation and stability of a large central vortex ring formed under the valve. A detailed description of the grid, the resolution requirements and the flow dynamics can be found in Ref. [25]. In this work, a constant temperature of 500 K is imposed on all boundaries while for the species zero flux boundary conditions are specified.
3
3.1. Geometry and computational grid For the calculation of the compression stroke the valve is closed at BDC, i.e. at the time instant of zero mass flux through the valve. The computational domain reduces to only the cylinder part of Fig. 1(a), and the initial conditions are obtained from the precursor mixing simulation at BDC discussed in the next subsection. At the chosen compression ratio of 12, the piston moves from 90 mm (BDC) to 7.5 mm (TDC). Two grids were employed constructed by extruding horizontal slices in the axial direction. The first grid containing 2436 spectral elements on the x y plane (Fig. 1) and 60 layers of elements in the axial direction clustered towards the piston and cylinder head is used for crank angles between 180° and 306°. For a polynomial order of n ¼ 9, an average axial resolution of 180 lm at BDC and 50 lm at 306° is obtained with a total of approximately 90 million nodes. In the radial and azimuthal directions the resolution is approximately 150 lm, while refinement close to the walls results in a resolution of better than 13.8 lm. In order to satisfy the increased resolution requirements during the last part of the compression stroke (>306°CA), a second grid with 3410 elements on the horizontal slice and 30 layers of elements in the axial direction was used and the fields were interpolated using high order spectral interpolation. By increasing the polynomial order to n ¼ 11 (corresponding to 135 million nodes) the average axial resolution is 110 lm at 306° and 30 lm at TDC. In the radial and azimuthal directions the resolution is 85 lm in the bulk of the domain
y
r φ
x
90
z
(a)
(b)
Fig. 1. (a) Schematic view of the domain of the open valve setup experimentally studied in Ref. [24], (b) axial slice of the second mesh used in the compression simulation.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
4
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
and better than 9.75 lm close to the walls. The pffiffiffiffiffiffiffiffiffiffi corresponding y þ ¼ ym sw =q value for the first mesh point at the wall is below 0.35 until 306°CA and below 0.9 until TDC, where y is the dimensional distance from the wall, sw the wall shear stress, and m; q the mixture viscosity and density. The resolution requirements from BDC up to 310°CA are comparable to those of the incompressible flow simulation [25] where they were assessed by a mesh resolution study and by comparison with the experimental data of [24]. As described below, with increasing pressure the size of the smallest scales decrease, resulting in significantly increased resolution requirements close to TDC. The Kolmogorov length scale gK ¼ 1=4 ðm3=Þ at TDC is gK ¼ 12:8 lm, where ¼ 2m sij sij is the mean dissipation rate computed from the strain tensor sij ¼ @ui =@xj þ @uj =@xi . At TDC, the radial and azimuthal resolution is approximately 6:6gK , while in the axial direction dz 2:3gK , which for the employed high order spectral element discretization method is enough to capture accurately most of the dissipation, and thus provide reliable first and second order statistics as pointed out in Ref. [26]. The resolution is also fine enough to resolve the mixing scales (i.e. 1=4 the Batchelor scale gB ¼ ðmD2 =Þ and the Obuk1=4 3 hov–Corrsin scale gOC ¼ ðD =Þ ), which are larger than gK since the Schmidt number is Sc ¼ m=D < 1. 3.2. Initial conditions In order to obtain realistic initial conditions for the temperature and species fields at the BDC of the compression stroke, the intake stroke was simulated as follows: starting with the velocity field at TDC of the third cycle of the cold flow simulation (so that the effect of the quiescent initial field is minimized) of the valve/cylinder assembly [25], the mixture inside the cylinder was assumed to be a homogeneous mixture with Y O2 ¼ 0:115; Y N2 ¼ 0:756 and Y H2 O ¼ 0:129 at 900 K and atmospheric pressure; this mixture corresponds to complete combustion of a / ¼ 0:5 H2/ air mixture and emulates EGR conditions. At the inflow, a homogeneous / ¼ 0:5 H2/air mixture at 500 K was introduced and a constant temperature of 500 K was imposed on all boundaries; for the species, zero flux boundary conditions were specified. The piston speed was 560 rpm and the Reynolds number based on the maximum mean jet speed, the cylinder radius and the viscosity of the intake mixture (m = 4.52 105 m2/s) is Re = 13,643. The downward moving piston draws the colder H2/air mixture into the cylinder. During the intake stroke, a large central vortex ring is formed below the open valve (indicated at Bottom
Dead Center (BDC) by a pressure isosurface in the velocity magnitude plot of Fig. 2(a)), which entrains hot EGR gases, while regions of lower temperature form close to the cylinder liner and the piston head. The initial conditions for temperature and species were found to have a weak effect on the conditions at the end of the compression stroke. As described below, at TDC the species are nearly homogeneously distributed whereas the temperature is primarily affected by heat transfer to and from the walls and not by its initial distribution. The mean velocity field is obtained by averaging in the azimuthal direction. For temperature, statistics are also collected within an inner cylinder of radius 30 mm and axial extend (from the cylinder to the piston) of 7.5 mm until 270°CA and 2 mm after 270°CA. This was done because averages in the whole domain are strongly influenced by the cooler near wall region. In the following, the volume average of any variable w in the whole cylinder will be denoted as hwiV and in the inner cylinder as hwiV ;i , whereas its magnitude is wmag ¼ jwj. The integral length scale Lw is derived by first extracting w along a circle at a fixed radial and axial position and then by integrating its autocorrelation to the point where 10% of its maximum value q is ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi reached [30]. The velocity magnitude is jvj ¼
v2z þ v2r þ v2u , while
the turbulent kinetic energy is defined as 02 02 with v0r ; v0u and v0z being the þ v þ v k ¼ 12 v02 r u z fluctuation velocities and T 0 the fluctuation temperature. 4. Results 4.1. Compression stroke The instantaneous distributions of velocity magnitude, temperature and Y H2 O at BDC are plotted in Fig. 2(a). The dominant flow feature is a large central vortex ring as illustrated in the figure by the pressure isosurface at p = 3.55 Pa; the white arrow indicates the clockwise flow direction around the vortex ring. Entrained inside the vortex ring core is fluid with higher temperature and water vapor mass fraction. In contrast, Fig. 2(b) shows the same fields at the end of compression (TDC) on an axial slice at 3.75 mm below the cylinder head. The temperature field shows high fluctuations in the whole domain with steep gradients close to the cylinder wall. In wall normal units the temperature gradient in the viscous sublayer (y þ < 5) follows the relation T þ ¼ Pr y þ , as reported in Ref. [27] for non-isothermal channel flows. In the inner cylinder region temperatures below 1050 K and above
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
| | [m/s]
| | [m/s]
9.8
3.36
0
0
T [K]
T [K]
600
1225
500
900
Y H2O
Y H2O
0.045
0.0325
0.015
0.0304
(a)
5
(b)
Fig. 2. Velocity magnitude, temperature and water mass fraction on (a) a vertical slice through the axis at BDC, and (b) on a axial slice 3.75 mm below the cylinder head at TDC.
30
2000
Isentropic compression DNS
20
1500
15 10
1000
Temperature [K]
25
pressure [atm]
1200 K can be observed, while the average temperature is close to 1110 K. Water vapor is practically homogeneously distributed inside the cylinder, as indicated by the narrow range of the color scale. The volume averaged pressure and temperature variation in time is shown in Fig. 3 (solid lines) and compared with the values computed for an isentropic compression using variable specific heat capacity ratio c. From BDC to TDC the mean temperature and pressure increase by a factor of 1.93 and 23.09, respectively. This results in a reduction of the kinematic viscosity by a factor of approximately 6.5 which strongly affects the evolution of turbulence. The differences between DNS and the isentropic curves are very small in the first half and increase significantly in the second half of the compression, reaching 3.9 atm and 148 K at TDC. This increase is due to the strongly increased heat transfer towards TDC, related to the decreasing thermal boundary layer thickness and the increasing temperature during compression. The pressure peaks 1.54°CA before TDC (thermodynamic loss angle), which is in
5 180
225
270
315
500 360
piston position [°CA] Fig. 3. Comparison of in-cylinder pressure and mean temperature with the values computed for isentropic compression.
the range of values reported in Ref. [28]. The relative high heat losses are related to the low engine speed (560 rpm) and the relatively cool wall temperatures compared to the average mixture temperature at BDC.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
6
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx 0
The variations of k, and T in the whole and the inner cylinder during the compression stroke are reported in Table 1. The decreasing kinematic viscosity results in higher velocity gradients around the smaller flow structures, and thus increased dissipation during the compression stroke. The peak is reached at 346° CA, i.e. before TDC, since close to TDC the turbulent kinetic energy is strongly reduced by the high dissipation level and zero piston speed. As shown in Ref. [2], temperature fluctuations can significantly affect the timing of autoignition in HCCI combustion and the corresponding peak of the pressure curve. Temperature fluctuations drop from BDC to 225°CA but increase significantly after 270°CA, mainly as a result of
the increasing difference between the average temperature in the cylinder and the fixed wall temperatures; this behavior is in agreement with the LIF measurements reported in Refs. [14,15]. Since autoignition is more likely inside the hotter inner cylinder zone, a more useful indicator for the description of the thermal stratification is the average fluctuations in the inner temperature 0:5 cylinder T 02 V ;i . This quantity decreases from 12.20 to 6.81 K in the first half of the compression stroke, possibly due to the homogenization of the temperature field in the inner cylinder. This trend reverses 0:5 in the second half of the stroke and T 02 V ;i increases significantly. As TDC is approached, the steepening of the velocity and temperature gradients near the walls results in increased heat transfer. Moreover localized ejections transporting cold fluid from the walls towards the cylinder center were found to be the main mechanism for increasing the rate of convective heat transfer during compression. In Fig. 4, the probability density functions (PDFs) of the turbulent kinetic energy and of temperature within the whole cylinder are plotted for several piston positions during the compression stroke. At 180°CA, PDF(k) shows a relatively flat profile with a peak at approximately 0.73 m2/s2. During compression the peak becomes more pronounced, but the position of the maximum remains at values close to 0.7 m2/s2. The relatively
Table 1 Variation of the volume averaged k, and the rms T during compression (in m2 =s2 ; m2 =s3 ; K). 02 0:5 02 0:5 °CA hk iV hiV T V T V ;i 180 225 270 306 330 346 360
5.70 3.78 2.50 1.83 1.79 1.58 0.98
820.00 536.50 973.52 3502.73 19281.7 36641.3 25510.0
16.80 14.84 29.69 58.04 92.99 110.7 119.9
12.20 6.73 6.81 13.41 32.29 39.15 47.31
0.8 0.7
PDF(k)
(a)
180°CA 225°CA 270°CA 306°CA 330°CA 346°CA 360°CA
0.6 0.5 0.4 0.3 0.2 0.1 0
0
1
2
3
4
6
5
7
8
2 2
k [m /s ] 0.06
(b)
PDF(T)
0.05 0.04 0.03 0.02 0.01 0 500
600
700
800
900
1000
1100
1200
T [K] Fig. 4. Probability density functions of (a) turbulent kinetic energy, and (b) temperature at several piston positions during compression.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
large difference in the peak values at 346° and 360°CA can be attributed to the increased dissipation rate towards TDC. According to [9] the direct effect of velocity fluctuations on autoignition is small compared to the effect of thermal stratification. Hence, turbulence is indirectly influencing autoignition by enhancing thermal stratification due to convective heat transport close to the walls. The influence of the initial conditions on temperature can be seen by comparing the PDF ðT Þ at 180° and 225°CA. The minimum and maximum temperatures are similar, but the peak becomes more pronounced and shifts to a higher value, as a result of the temperature homogenization during the early compression stage discussed above. The higher peaks at 225° and 270°CA are due to the low temperature fluctuations in the inner cylinder (Table 1). The PDF shape differs significantly from the almost Gaussian distribution reported in the experimental study of [15], possibly because of the difficulty in resolving the thermal boundary layers in the two-dimensional measurements. Between 306° and 330°CA, due to the effect of wall heat transfer, the maximum PDF ðT Þ values decrease strongly and close to TDC (346° and 360°CA), the PDF ðT Þ becomes almost flat. This can be explained by the high wall heat losses during this period. The maximum temperature 1234.5 K is reached at 360°CA although the average temperature is nearly the same as at 346°CA (Table 1). The temporal development of PDF ðT Þ indicates that the temperature initial conditions are washed out after the first half of the compression stroke, implying that, at least for this setup, the thermal stratification at TDC cannot be controlled by changing the temperature field at BDC. Stronger effects can be expected by modifying the wall temperature, as also found in Ref. [29], where it was shown that lower coolant temperature had a strong effect on the heat release rate for HCCI operation.
7
4.2. Behavior at Top Dead Center The distribution of v0z on an axial slice 3.75 mm below the cylinder head is shown in Fig. 5(a), indicating an inhomogeneous distribution within the domain. The average values and the PDFs of the three fluctuation velocity components are plotted in Fig. 5(b), where their volume averages values are also reported. The PDFs of the radial and azimuthal fluctuation velocities as well as their average values in the radial and azimuthal directions are nearlyDidentical, E while the average in the axial 0 direction vz;mag is approximately 25% lower. V
Accordingly, the peak PDF value is higher, possibly due to the dumping of axial fluctuations in the narrow gap between the piston and the cylinder head. The spatial distribution of the integral length scales of k and temperature are plotted in Fig. 6. Values close to the cylinder center are not shown, since the number of statistically independent points is too low to obtain meaningful values. With the exception of two zones close to the cylinder head, Lk lies in the range between 1 and 2.5 mm, while LT varies between 3 and 4 mm in the cylinder center and decreases towards the walls due to the increased temperature fluctuations in the thermal boundary layer. The average length scale of k is in the whole domain approximately 25% smaller compared to the average scale of the temperature. The PDFs of Lk and LT plotted in Fig. 7 have a similar shape and their peaks are located at approximately 1.1 mm. In contrast, the PDF ðgk Þ is more symmetric and gK ranges from 0.005 to 0.025 with an average value of 12.8 lm. The smaller gk values are mainly located in the relatively cooler thermal boundary layers since m / T 1:5 . The ratio between hLk iV and hgk iV at TDC is approximately 161.4. The low value of hgk iV close to TDC is mainly due to the increase in pressure followed by a reduced kinematic viscosity and a
1.5 1.2
2.8
z
[m/s]
0.9
<
z,mag
<
r,mag
>V = 0.518 m/s
<
φ,mag
>V = 0.691m/s >V = 0.671 m/s
0.6 0.3
0
0
1
2
3
fluctuation velocities [m/s]
(a)
(b)
Fig. 5. (a) Axial fluctuation velocity on a horizontal slice 3.75 mm below the cylinder head; (b) PDF of the axial, radial and azimuthal fluctuations velocity components at TDC.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
8
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
Lk
cylinder wall
cylinder axis
cylinder head
piston
LT
4 3
cylinder axis
z=0 2 1
r=0
r = 30
r = 15
z = -7.5 r = 37.5
Fig. 6. Integral length scales in mm of the turbulent kinetic energy and the temperature at TDC.
LT Lk ηK 100 < ηK>V = 12.8 [μm]
PDF (L T, k )
0.4 0.3 0.2
50
0.1 0
2
4
6
PDF (ηK )
150 0.5
0
8
L [mm] 0
0.005
0.01
0.015
0.02
0.025
ηK [mm]
Fig. 7. PDF of the integral length scales of turbulent kinetic energy and temperature and the Kolmogorov scale at TDC. 15000
YH2
YH2
YO2
10000
0.0109
YH2O
5000
-0.001
0
averages. The higher maximum of PDF ðY H2 Þ can be explained by the high diffusivity of hydrogen, which has a Lewis number approximately three times smaller than the other two species. The PDFs of Y O2 and Y H2 O show approximately equal deviations and the local Y H2 O maxima correlate with the local Y O2 minima. However the deviations in the species mass fractions are very small compared to those of temperature, in agreement with the results in Ref. [14].
0.001
deviation from < Y >V
Fig. 8. PDFs of H2, O2 and H2O mass fractions at TDC. The zero value on the x-axis corresponds to the average species mass fraction. The horizontal slice are located 3.75 mm below the cylinder head.
higher dissipation rate. On the other hand, the influence of the increasing temperature on hgk iV is less pronounced, since temperature is less sensitive to volume changes compared to pressure. Thus, the compression ratio has a strong impact on the size distribution of the turbulent scales, which in the case of a propagating flame front (as in an Otto engine) is expected to have a significant impact on the combustion process and the engine efficiency. The PDFs of Y H2 ; Y O2 and Y H2 O reported in Fig. 8 show the deviations from the volumetric
5. Conclusions Direct numerical simulation was performed to study the effect of compression on the flow, temperature and composition fields in the closed cylinder of the valve-piston assembly investigated experimentally [24] and numerically [25]. The results show significant changes of the flow and temperature during compression: The reduction of the kinematic viscosity due to the increasing pressure leads to smaller turbulent length scales and increased dissipation rates close to TDC. Temperature fluctuations decrease during the first half of the compression stroke, and, in agreement with the experimental findings of [14,15], increase significantly during the second half, mainly as a result of the heat transfer to the wall and turbulent convective transport of cooler gas from the walls towards the inner part of the cylinder. The influence of the temperature initial conditions at BDC on the temperature fluctuations at TDC is found to be very small, while, in agreement with [14], the species distributions at TDC are almost homogeneous due to the very efficient mixing during compression. The flow field at TDC is anisotropic, and the average fluctuation velocities in the axial direction are approximately 30% smaller compared to those in the radial and azimuthal direction, possibly due to the smaller axial dimension at TDC. Due to thermal diffusion, the integral length scale hLT i is approximately 25% higher compared to hLk i. The average Kolmogorov scale is 12.8 lm and around 160 times smaller than hLk i. While only one cycle is investigated here, ongoing simulations of the compression stroke using the flow fields from different cycles of the multicycle simulation in Ref. [25] showed that although large flow structures have an influence on the flow and temperature field, the observed trends for the evolution of mean and rms temperatures, k, and the species mixing are very similar for all calculated cycles. We are currently analyzing wall heat transfer in this setup and the phenomena associated with velocity and thermal boundary layers. Future work will consider reactive conditions and investigate the autoignition behavior of the mixture.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097
M. Schmitt et al. / Proceedings of the Combustion Institute xxx (2014) xxx–xxx
Acknowledgments Financial support of the Swiss National Science Foundation (SNF) under project number 200021_135514 is gratefully acknowledged. The work was also supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID 189. References [1] J.E. Dec, Proc. Combust. Inst. 32 (2009) 2727–2742. [2] E.R. Hawkes, R. Sankaran, H.G. Im, P.P. Pe´bay, J.H. Chen, Combust. Flame 30 (2005) 875–882. [3] M. Yao, Z. Zheng, H. Liu, Prog. Energy Combust. Sci. 35 (2009) 398–437. [4] F. Zhao, T. Asmus, D. Assanis, J. Dec, J. Eng, P. Najt, Homogeneous Charge Compression Ignition (HCCI) Engines, Key Research and Development Issues, SAE Inc., Warrendale, PA, 2002, 1–14 and 531–541. [5] M. Sjo¨berg, J.E. Dec, N.P. Cernansky, SAE Technical paper, 2005. 2005-01-0113. [6] J.H. Chen, E.R. Hawkes, R. Sankaran, S.D. Mason, H.G. Im, Combust. Flame 145 (2006) 128–144. [7] R. Sankaran, H.G. Im, E.R. Hawkes, J.H. Chen, Proc. Combust. Inst. 145 (2006) 145–159. [8] C.S. Yoo, T. Lu, J.H. Chen, C.K. Law, Combust. Flame 158 (2011) 1727–1741. [9] M.B. Luong, Z. Luo, T. Lu, S.H. Chung, C.S. Yoo, Combust. Flame 160 (2013) 2038–2047. [10] G. Bansal, H.G. Im, S.R. Lee, Proc. Combust. Inst. 32 (2009) 1083–1090. [11] J.O. Hinze, Turbulence, McGraw-Hill, New York, 1975. [12] E. Baum, B. Peterson, C. Surmann, D. Michaelis, B. Bo¨hm, A. Dreizler, Proc. Combust. Inst. 34 (2013) 2903–2910. [13] B. Peterson, E. Baum, B. Bo¨hm, V. Sick, A. Dreizler, Proc. Combust. Inst. 34 (2013) 3653–3660. [14] J. Snyder, N. Dronniou, J. Dec, R. Hanson, SAE Technical paper, 2011. -01-1291.
9
[15] S.A. Kaiser, M. Schild, C. Schulz, Proc. Combust. Inst. 34 (2013) 2911–2919. [16] E. Gu¨ntsch, R. Friedrich, Flow Simulation with High-Performance Computers II, Springer, 1996, pp. 213–226. [17] P.F. Fischer, J.W. Lottes, S.G. Kerkemeier, http:// nek5000.mcs.anl.gov, (2008). [18] M.O. Deville, P.F. Fischer, E.H. Mund, Cambridge U. Press, New York, 2002. [19] A.G. Tomboulides, J. Lee, S.A. Orszag, J. Sci. Comput. 12 (1997) 139–167. [20] G.D. Byrne, A.C. Hindmarsh, Int. J. High Perform. Comput. Appl. 13 (1999) 354–365. [21] L.W. Ho, A.T. Patera, Comput. Methods Appl. Mech. Eng. 80 (1990) 355–366. [22] R.J. Kee, F.M. Rupley, J.A. Miller, Chemkin II: A Fortran Chemical Kinetics Package for the Analysis of Gas-phase Chemical Kinetics, Report No. SAND89-8009B, Sandia National Laboratories, 1996. [23] R.J. Kee, G. Dixon-Lewis, J. Warnatz. M.E. Coltrin, J.A. Miller, A Fortran Computer Code Package for the Evaluation of Gas-Phase Multicomponent Transport Properties, Report No. SAND868246, Sandia National Laboratories, 1996. [24] A.P. Morse, J.H. Whitelaw, M. Yianneskis, J. Fluids Eng. 102 (1979) 208–2016. [25] M. Schmitt, C.E. Frouzakis, A.G. Tomboulides, Y.M. Wright, K. Boulouchos, Phys. Fluids 26 (2014) 035105. [26] P. Moin, K. Mahesh, Ann. Rev. Fluid Mech. 30 (1998) 539–578. [27] H. Kawamura, K. Ohsaka, H. Abe, K. Yamamoto, Int. J. Heat Fluid Flow 19 (5) (1998) 482– 491. [28] R. Pischinger, M. Klell, T. Sams, Thermodynamik der Verbrennungskraftmaschine: Der Fahrzeugantrieb, Springer-Verlag, Vienna, 2010. [29] M. Sjo¨berg, J.E. Dec, A. Babajimopoulos, D. Assanis, SAE paper, 2004-01-2994. [30] W.J. Feiereisen, W.C. Reynolds, J.H. Ferziger, Numerical Simulation of a Compressible Homogeneous Turbulent Shear Flow, Tech. Rep., Stanford University, 1989.
Please cite this article in press as: M. Schmitt et al., Proc. Combust. Inst. (2014), http://dx.doi.org/ 10.1016/j.proci.2014.06.097