© E. Yildar et al., published by IFP Energies nouvelles, 2017

Introduction

Considering the limited fossil fuel resources and strict limitation of pollutant emissions, the Homogeneous Charge Compression Ignition (HCCI) engine is a promising technology as it combines advantages of both the diesel and spark ignition engines [1–4]. Firstly, the efficiency of HCCI is comparably high for the diesel engine due to the high compression ratios. Secondly, in HCCI due to the highly diluted charge, high local peak temperatures are avoided, thus reducing the NO_x formation [1–4]. Furthermore, since the charge is well-mixed, soot formation is prevented [1,2,4]. However, since high peak pressures and heat releases occur at high loads, the HCCI technology is applied for a limited range of load.

To overcome this issue, strategies are developed to control HCCI combustion which refers to the Controlled Auto Ignition (CAI) concept [2,3]. Exhaust Gas Recirculation (EGR) as a dilutant of the fresh fuel-air mixture is one of the strategies used for controlling HCCI combustion [1–3]. The EGR is either achieved by trapping hot exhaust gas in-cylinder with the Negative Valve Overlap (NVO) concept or burnt gas load during intake, so-called Burnt Gas Re-Breathing (BGRB) [2]. The combustion timing of CAI is controlled by the chemical kinetics and highly depends on the properties of the mixture field [2,5–10]. Inhomogeneities in the mixture and temperature field designate the combustion characteristics of CAI [2,5]. Thus, the prediction of the auto-ignition process requires an accurate description of the chemistry within the whole range of thermodynamic conditions given by thermal and composition inhomogeneities.

For a better understanding of the chemical and physical phenomena in CAI engines CFD is a powerful tool. In industrial applications due to their convenient computational costs Reynolds-Averaged Navier Stokes (RANS) models are used for simulation of turbulent flows. However, RANS models solve the temporal mean value of the relevant quantities of the system. With regard to a more accurate prediction of situations that intrinsically depend on temporal and spatial variations associated with turbulence, the Large Eddy Simulation (LES) has established as a research tool. The large-scale motions are explicitly resolved to accurately predict the conditions prior to ignition. In particular those are the mixing of the injected fuel with the air, the EGR distribution within in the cylinder and the thermal conditions associated with them.

In order to accurately predict the auto-ignition process the reaction kinetics as well as the turbulence–chemistry interaction have to be considered. For the complex reaction chemistry of the auto-ignition process, detailed reaction mechanisms have been developed [11,12]. Coupling these detailed reaction kinetics directly with a CFD code with full spatial and temporal resolution is, in principle, possible. However, such detailed chemistry simulations are computationally extremely expensive, rendering them infeasible as simulation tools for practical systems of technical scale. One way is to use strongly reduced kinetic schemes which consider only the most important species and their reactions. Methods such as Computational Singular Perturbation (CSP) [13] and the intrinsically low dimensional manifold (ILDM) [14] perform the reduction automatically based on a mathematical analysis that considers the time scales of the underlying chemical kinetics. Approaches like the flamelet model [15] explicitly account for the coupling of the chemical kinetics and molecular transport. For the prediction of finite rate chemistry effects, flamelet models have been extended with a progress variable approach [16,17]; variants are the Flamelet Generated Manifolds (FGM) [18], the Flamelet Prolongation of ILDM (FPI) [19], and the Reaction Diffusion Manifold (REDIM) method [14]. For cases where the chemical reaction occurs without a strong coupling to diffusive processes, for instance, auto-ignition in mixtures with weak scalar gradients, other models have been designed [20]. A simple approach is to map pre-computed information about the chemical dynamics in detailed homogeneous reaction systems onto a reduced model, whose source term is then tabulated and used in CFD simulations [21].

Regarding the turbulence–chemistry interaction, additional uncertainties arise from the averaging within RANS simulations or filtering within LES yielding the subgrid variations. Within ICE simulations [22] the models are chosen according to the physics required for the given combustion mode. Accordingly, for spark ignition engines usually the Flame Surface Density (FSD) approach [23,24], the G-equation [25], or the Artificially Thickened Flame (ATF) model [26] are employed to treat the flame propagation. For applications with auto-ignition different mechanisms are dominating [2] and those simulations often use Probability Density Functions (pdf) to account for the unresolved scalar distribution and accordingly obtain a more accurate evaluation of the non-linear relations. However, assumptions have to be made to obtain the pdf. Within RANS, approaches based on presumed shapes [27,28] (e.g. β-pdf) as well as transported pdfs [29] have been used to simulate the ignition in engines. Yu et al. [30–32] conducted several LES studies based on a three-dimensional chemistry table. For the subgrid treatment they employed a two-dimensional pdf build by assuming independent Gaussian distributions for the total enthalpy and the reaction progress variable.

LES has been applied to HCCI in several works. In [30–32] a detailed reaction mechanism is coupled with the LES code using a progress variable approach. Within these consecutive works the authors investigated the sensitivity of the HCCI combustion with respect to an initially imposed turbulence level and temperature inhomogeneity. In [33,34] the knock phenomenon in HCCI engine is studied with detailed chemical kinetics coupled to LES codes. A model called Ignition to Propagation Reduced Scheme (IPRS) is introduced in [35] to predict abnormal combustions in IC engines in the context of LES. This model uses a single reduced chemical scheme (with two reactions) to describe both Auto-Ignition (AI) and propagation.

Particularly, in this work the gasoline fueled CAI engine with the NVO concept built and experimentally investigated at the Institute of Internal Combustion Engines, Karlsruhe Institute of Technology has been simulated. The computations have been performed using the KIVA-4mpi code where the LES methodology [36,37] is adopted. In order to capture the complex chemistry of auto-ignition, a detailed reaction mechanism [12] is pre-tabulated and coupled with the CFD code. A progress variable is defined to reflect the chemical progress and its dynamics in the reacting system. A multi-dimensional chemistry table is generated to cover inhomogeneous engine conditions being the air–fuel ratio (λ), the unburnt gas temperature (T_u), pressure (p), and EGR. The objectives of the work are:

• 1.

  Demonstrate the application of the multi-dimensional chemistry table to a real engine.

• 2.

  Perform a statistical comparison to measured pressure curves.

• 3.

  Characterize the thermal and mixture inhomogeneities to explain the occurrence of ignition kernels.

• 4.

  Provide an understanding of the combustion process and the interaction between cycles.

1 Experimental Setup and Configuration

The experiments are performed at the Institute of Internal Combustion Engines, Karlsruhe Institute of Technology on a single-cylinder four-stroke Direct Injection (DI) gasoline engine (Fig. 1) with four canted valves, derived from a BMW F650 motorcycle manufactured by Rotax with several constructive modifications for numerous analyses. The valves located on the intake side of the engine have a curved bottom contour, whereas on the exhaust side the valves have an even bottom. The exhaust ducts are straight, whereas the intake ducts are straight only for the second half of the duct length. The maximum valve lifts are 0.186 cm and 0.137 cm for the inlet and exhaust valves, respectively. In Table 1 the geometric properties of the engine are summarized.

Two conventional hydrodynamic phase shifters for intake and exhaust are installed to regulate the residual gas fraction. Extensions of the camshafts are provided with an angle encoder brand while the Top Dead Center (TDC) signal is set to the maximum valve lift. In this way, the information about the phase position is available through the indexing system. The cam profiles are calculated with a 1D-simulation for a previously defined operating point at 2000 rpm with an IMEP of 2.0 bar.

A high-pressure multihole injector HDEV5 from Bosch GmbH is placed centrally in the head of the cylinder with an installation angle of 18°. The injected fuel is EuroSuper with overall 84.1 mass% carbon, 13.7 mass% hydrogen and 2.2 mass% oxygen. For the simulations Toluene Reference Fuel (TRF) is used, which consists of 85% iso-octane, 5% n-heptane and 10% toluene (liquid volume%) and has the same Research Octane Number (RON) of 95 as the EuroSuper. The fuel mass is m_inj = 10.5 mg, being injected with a pressure of 100 bar and a temperature of T_inj = 350 K. The processes throughout a cycle will be illustrated in Section 3.1.

A high-pressure sensor was adapted between the exhaust valves. The lateral spark plug visible on the left of Figure 1 was needed just for starting operations. The cooling water and oil are set to 95 °C with an external conditioning system. The intake air is also conditioned, and the temperature is set to 25 °C with a humidity of 35% and a pressure of 1050 mbar absolute. The air mass was measured with a mass flow sensor from Bosch GmbH, the fuel mass is measured with a PLU121 from Pierburg. The air–fuel ratio λ is measured by an LSU4.9 from Bosch GmbH and additionally calculated from the measured exhaust emissions with the Brettschneider formula. The engine was also fully indicated, for intake and exhaust pressure an uncooled pressure sensor 4045 from Kistler is used. The in-cylinder pressure is measured with an uncooled 6054 sensor.

Based on 200 measured cycles the average EGR rate defined by (1) is approximately 65%. This is calculated with an in-house pressure trace analysis tool at the Institute of Internal Combustion Engines. The calculations are done on single cycles, in which a Savitzky-Golay filter is used to reduce signal noise or on averaged pressure curves.

2 Numerical Modeling

The simulations are conducted with the engine code KIVA-4mpi. The code utilizes the ALE approach to solve the governing equations on a moving mesh that follows the geometry. Convective fluxes are approximated by a quasi-second-order upwind discretization and the Crank Nicolson scheme is used for the time integration. Details about the numerical methods can be found in Torres and Trujillo [38].

Within this work, the flow field is predicted by the LES-filtered, density weighted transport equations for mass, momentum and internal energy reading (2) (3) (4)

Herein, the subgrid flux of momentum was closed using the standard Smagorinsky model [39] (5) where C_s = 0.1 denotes the Smagorinsky constant and S̃ is the filtered rate of strain. No special wall treatment is applied. The unresolved transport for any scalar ϕ is treated by the gradient approach (6)where σ_ϕ denoted the Prandtl or Schmidt number being set to 0.7. Likewise the heat flux towards the wall follows from (7)where denotes the wall normal direction. The coupled system is closed by the filtered equation of state (8)

2.1 Chemistry Treatment

To predict the auto-ignition process, the LES code is coupled with a chemistry tabulation approach specifically developed for this purpose. The chemistry reduction is presented in [20] and the coupling to the LES code used here has been addressed in a foregoing work by Yildar et al. [40]. The approach will be summarized in the following. Further details can be found in the references while the focus of this work is the application of the method to a real engine.

The chemistry table for the Toluene Reference Fuel (TRF) is constructed by detailed chemistry simulations employing the reaction mechanisms of Andrae and Head [12] that accounts for 137 species in 633 reactions. According to the physical situation of the auto-ignition, isobaric, adiabatic homogeneous reaction trajectories given by (9) (10) (11) are computed using the HOMREA code [41]. To include all possible states found in the engine, multiple trajectories are calculated for the relevant conditions. At this, the operation of the engine with spray injection and a high amount of EGR yields a four-dimensional space of possible initial conditions for the ignition. This multi-dimensional thermo-chemical space is covered by considering in each dimension:

• 14 values for the air–fuel ratio (0.3 ≤ λ ≤ 2.0)

• 25 unburnt gas temperatures (500 K < T_u < 1250 K )

• 11 ratios of EGR (0% < EGR < 100%)

• 9 pressure levels (10 bar < p < 63 bar)

yielding a total of 34 650 detailed chemistry simulations for the table generation. The chemistry table then accordingly spans into five directions with the table access parameters

• Z the mixture fraction representing the air–fuel ratio in a numerical convenient way

• T_u the unburnt gas temperature

• Y _EGR being the EGR mass fraction as given by equation (1)

• p the pressure from the equation of state

• Y _PV the reaction progress variable

Hence, in addition to equations (2)–(4), it requires to solve for the three access parameters: (12) (13) (14) where describes the effect of droplet evaporation onto the mixture fraction and is the chemical source term of the reaction progress variable. This latter as well as the heat release ( in Eq. (4)), the specific heat, the mean molar mass, and viscosity are taken from the chemistry table during the simulation. In the current state of the model there is no further treatment of the chemistry–turbulence interaction at the subgrid level. We rely on the zero-dimensional auto-ignition being the dominant process and that the corresponding scalar fields representing the thermal and mixing state at ignition are sufficiently resolved. Further development and knowledge in that regard are desirable, but certain approaches like a pdf pre-integration seem currently unfeasible due to the high number of dimensions as well as the corresponding assumptions about their statistical dependency.

For the auto-ignition the progress variable is a non-trivial choice. As detailed in [20], within this work, we use the specific chemical entropy with the production rate: (15)

Integrating this specific chemical entropy production rate the system will deliver a monotonic variable that can be transported with the source term . The integrated entropy offers several advantages as a progress variable, including a guaranteed temporal monotonicity in closed adiabatic systems, and also a good sensitivity to the early stages of the chemical reaction, well before the onset of the main heat release. For a better interpretation the reaction progress can be normalized by: (16) as will be used in the results section.

To comply with the tabulation a further special treatment of the EGR and PV transport is required. During a cycle those are ordinary mass fractions. However, according to its definition, EGR is 100% after the exhaust valves closed. Accordingly the EGR of these gases is set to this value before the ordinary transport continues for the next cycle. In that regard the EGR is a passive scalar that should be viewed as a table access parameter whose connection to actual scalars is reset in between cycles. Consistently, also the progress variable gets reset, i.e. after combustion the burnt gases are defined to be the unburnt EGR of the consecutive one. The chemistry table is build such that this joint transition in access parameters is continuous for the coupling variables meaning that the variables restored from the table are identical before and after this step. For the non-reacting phases of the cycle where the thermo-chemical properties cannot be obtained from detailed chemistry simulations an additional table is employed where the properties simply follow from the thermodynamic state.

In Figure 2, three-dimensional plots of the five-dimensional tabulation, obtained by varying the independent variables λ, T_u, p, and EGR are visualized. For this figure, an engine relevant condition is chosen, namely λ = 1.52, T_u = 1110 K, p = 31 bar, and EGR = 60%. The reaction progress is represented by PV in the subfigures (a)–(d). The figure shows how the source term of progress variable changes during the chemical reaction for CAI engine-relevant ranges of λ, T_u, p and EGR. Comparing the slopes, the plots illustrate that p has a weaker influence on compared to λ, T_u and EGR. The EGR influence is provided with more detail in Figure 3 showing slices of the table. One can see that the source terms value can reduce by one order of magnitude when increasing the EGR from 50% to 70%. Furthermore the range of λ where significant reactivity is found reduces.

2.2 Computational Domain and Simulation Setup

The computational domain and the corresponding boundary conditions are illustrated in Figure 4. The extent to which the intake and exhaust ports are included in the simulations has been chosen to comply with the measurements available. As indicated, within the experiment, temperature and pressure have been monitored at 3.7 D and 4.5 D upstream of the intake and downstream the exhaust port, respectively. These data have been gathered fully temporal resolved for all cycles and the phase-averaged curve as provided in the bottom of Figure 4 has then been applied equally to all cycles in the simulation. Furthermore, the wall temperatures applied to the simulation are visualized in the top of Figure 4. Those were not directly available from the experiment. Accordingly, the temperatures of the intake and exhaust ports have been set to the averaged measured gas phase temperature. Furthermore, the temperature of all in-cylinder walls (head, piston, liner) could only be roughly estimated based on the cooling water.

The injector is not reproduced in the engine geometry. Therefore one point is defined as the injection point with the properties of the injector. For the spray injection, a solid cone spray profile with a square wave injection pulse is specified. The description of the spray follows the discrete droplet model of Dukowicz [42]. Due to the large memory requirement it was not possible to treat each droplet individually and therefore multiple droplets of similar properties are represented by a computational parcel. Accordingly, about 250 000 Lagrangian parcels have been tracked, each of them representing an initial liquid volume of 249 μm in diameter. The spray is treated with the Taylor Analogy Breakup (TAB) model [43] and several models to account for the fluid interactions as described in [44,45].

As depicted in Figure 5 the domain is discretized with a block-structured mesh containing 2.3 million cells at the bottom dead center (BDC) yielding an average cell size of about 0.8 mm. The time step size is dynamically adjusted by several regulations like CFL-criteria or the given maximum mesh movement per time step and is 3 × 10⁻⁸ s ≤ Δt ≤ 1 × 10⁻⁶ s.

To perform a comparison with the phase-averaged measurements as well as to obtain a solid database for our further analysis, a series of cycles have been computed. Especially due to the EGR of this engine, the operation shows a strong history effect spanning over several cycles. To increase the sample number we started two independent sets of simulations with different stochastic initial conditions. Each of them ran for 20 consecutive cycles out of which we removed the first ten for the statistical evaluation to avoid an influence by the startup process. Accordingly 20 valid cycles enter our phase averaged data being sufficient to provide averaged values but rather insufficient for statistics of higher orders (variances). Furthermore, as given by our setup we have ten consecutive cycles available for the analysis of cycle-to-cycle variations.

3 Results

This result presentation is divided into three parts. First, Section 3.1 outlines the processes taking place throughout a typical cycle to provide a basic understanding. Additionally, the computed pressure curve gets compared to measurements and potential uncertainties are discussed. Then, in Section 3.2, a detailed analysis is conducted to identify the thermal and mixing conditions that determine the auto-ignition. Finally, Section 3.3 considers the interaction of consecutive cycles to explain the occurrence of cycle-to-cycle variations.

3.1 Cycle Illustration and Pressure Curves

3.1.1 Cycle Illustration

The illustration of a cycle is done by means of the global quantities depicted in Figure 6 supplemented by the field quantities given in Figures 7 and 8. The top of Figure 6 shows the pressure evolution of a typical cycle spanning from −360 < cad < 360 (cad = crank angle degree) supplemented by the valve lift curves given in blue. For a better understanding of the process its preceding and subsequent cycle have been added with dashed lines. The process will be explained following the green numbered positions.

The first point at −500 cad represents the exhaust valve opened of the preceding cycle as given in the top row of Figure 7. Before this opening, by definition, the cylinder contains only exhaust gas leaving through the valve. However, the valve timing is chosen such that a considerable amount of this exhaust gas remains within the cylinder for the next cycle being the reason why we start the explanation at this point. According to this early closing of the valve being 87 cad before the TDC a significant pressure increases follows. It is worth noting the strong inhomogeneity of the exhaust gas temperature shown in Figure 7. Its causes are the combustion of different mixing conditions, locally different ignition histories as well as heat losses to the cold wall.

After the TDC we arrive at the second point being the fuel injection done at −300 cad (second row in Figure 7). The liquid fuel is injected into the exhaust gases of the preceding cycle by a multi-hole injector. The droplets are indicated on the right of Figure 7. Furthermore, the fuel mass fraction that already evaporated is given by the contours within the two slices. The middle column of Figure 7 shows the ratio of EGR. One can see that the injection direction spreads the fuel mostly towards the positive x-direction. Likewise, the fuel injected at 350 K causes regions of low temperature visible on the right.

After this short injection phase the inlet valves are opened at point three. Pure air enters at 317 K leading to a further reduction of the EGR and the temperature. At this stage, the third row in Figure 7 reveals certain similarities between the EGR and temperature field. Even though there was little initial correlation (see top row), from this point on, the hot EGR is the dominating factor determining the temperature away from walls as will be detailed in the next section.

Going over to point four at −25 cad given in the first row of Figure 8. On the left we introduce the air–fuel ratio λ. From this moment on it provides reasonable values to visualize the mixing state. However, it is reconstructed based on the mixture fraction and since the exhaust gas contains unburnt oxygen it does not represent the actual air-to-fuel ratio. Accordingly, the stoichiometric ratio is shifted to lower values (i.e. below one) as we will see in the next section. With regard to this mixing state it is important to mention that all regions in the cylinder are within the flammability limits by now. This point can be considered as the pre-auto-ignition state. No ignition took place yet but the evolution of the combustion process is mostly determined by the conditions found here. Compared to point three, the EGR became more homogeneous by diffusion but still significant spatial variations exist (colormap has been rescaled compared to point three). The temperature increased significantly by the compression, still, it is insufficient for the auto-ignition.

Finally at point five (−5 cad) the first ignition spots visualized by the normalized reaction progress variable (PV) in Figure 8 become visible. The EGR distribution is very similar and that is why we stated that the process is mostly pre-determined by point four. From this moment on the processes take place in a very short amount of time. Accordingly, the mixing state is mostly frozen and the compression acts in increasing the temperature of these frozen states until sufficiently low ignition delay times are reached at some favorable positions.

Point six being at 3 cad then describes a state where the combustion process is already well developed. The corresponding progress variable depicted on the right in the last row of Figure 7 indicates that some regions are already fully reacted while numerous ignitions spots of intermediate states can be observed.^a Compared to point five now the temperature shows larger variations being dominated by the combustion process and no longer by the EGR distribution. The combustion initially evolving only at the most favorable positions now increases the pressure such that also in the other regions sufficient activation energy is present.

3.1.2 Pressure Curves

Figure 9 compares the phase-averaged pressures predicted by the LES with the measurements. The motored case is also shown here to demonstrate that the numerical setup is suited for this configuration. Specifically, the pressure boundary conditions set in the inlet and exhaust pipes, together with the prescribed valve lift and piston movement lead to the correct trapped mass. These motored curves show the ordinary evolution imposed by the compression at which the LES and the experiment are in close agreement.

Going over to the fired operation, the results are given by the red curves. Compared to the motored case, initially it shows a lower pressure slope during compression. This is caused by the motored case having the larger trapped mass at which the hot exhaust gases of low density take up parts of the in-cylinder volume in the fired case. At about 5 cad the onset of combustion causes a distinct slope increase. Especially at the end of the combustion process differences between the simulation and the experiments are visible as detailed in the following.

As will become more clear in the next sections, the ignition process and the pressure curve are influenced by numerous parameters where uncertainties in the simulation exist. Taking a closer look at Figure 9, the LES shows a slightly lower pressure rise in the compression phase when compared to the experiment. It then crosses the experimental curve due to an earlier onset of ignition after which the slopes during the combustion phase are similar. As detailed in Table 2 the simulated peak pressure is about 2.3 bar higher and delayed by 1.9 cad. Besides the averaged pressure the envelopes have been added to the graph. While the variation of the pressure curve between different cycles was negligible in the motored case (not shown) a scatter is very visible in the fired operation. At this the simulation predicts a larger range being 18% and 14% of the average pressure in the LES and experiment, respectively. Accordingly, the LES overpredicts the average pressure as well as the cycle-to-cycle variation. The potential causes are discussed in the following. However, those can interfere with each other and some relations may require the next sections to become more clear.

The initial compression slope indicates a slightly lower trapped mass in the LES. After that, the curve indicates that a larger amount of work is delivered in the simulation. For the latter, one possible reason is an imperfect reaction in the experiment. Indeed, components such as CO and HC were found in the exhaust gas indicating that the fuel was not fully burned. The simulation showed an average reaction progress of 98% seeming more complete. However, the measurements do only provide a qualitative information such that the actual contribution of this effect remains uncertain.

Within the simulation, the spray treatment is in general prone to errors. Deviations associated with them should be rather local since the global energy balance is hardly affected as long as the fuel fully evaporates. Besides uncertainties in the initial conditions such as the droplet size distribution and spray angle, the breakup and evaporation process require modeling. Since the fuel distribution has a strong impact onto the ignition, the simulation might be sensitive to these errors. Specifically one major concern was related to the grouping of the droplets to parcels. However, as we analyzed for this specific engine, the spray treatment should not be critical: As visible in Figure 6 the spray gets injected at −300 cad into a high temperature environment of about 800 K. This causes a very quick evaporation process taking only 10 cad in the simulation. This means, the evaporation is completed before any further processes take place in the LES and likely also in the experiment. Accordingly, inaccuracies in this very short period should barely affect the overall process.

Also the chemistry treatment is obviously a suspect for uncertainties when considering the complex reaction paths involved in the auto-ignition process. Specifically the detailed chemistry simulations used in the pre-tabulation process do not resemble the actual engine operation. At this, the real process is not isobaric and homogeneous and the corresponding equations (9)–(11) will have contributions on their right hand side. However, in [40] we evaluated the assumptions including compression and enthalpy diffusion and they proved to be accurate for the gradients found in the real engine. Indeed, the ignition process itself shows a reasonable slope comparable to the measurements.

The heat transfer treatment is a further uncertainty impacting the global behavior. This comprises the accuracy of the temperature set at the walls as well as the computation of the heat flux. The latter (see Eq. (7)) obviously inherits the inaccuracies of the turbulence model. Regarding the wall temperatures we outlined our approximations in Section 2.2 and Figure 4. As mentioned therein, a simplification due to the lack of more data was to assume identical temperatures for all in-cylinders walls. It is well known that certain variations exist [46,47], specifically the temperature of the exhaust valves is increased. On one hand, the ignition depends on hot spots where this might be important. On the other hand however, as detailed in the next sections, even a typical exhaust valve temperature of 600 K is significantly below the auto-ignition temperature of about 1000 K where significant reaction takes place. Accordingly, the boundary temperatures may affect the global balance but not the local ignition being dominated by other conditions (see next sections).

Finally, and also related to the boundary conditions, the interaction of heat and mass fluxes does represent a major uncertainty. Especially in interaction with the high rate of EGR we believe those significantly contribute to the deviations observed in the average pressure as well as in the envelope. Regarding the mean values, we observed that the LES generally shows a lower EGR rate than the experiment (approx. 58% to 65%). Throughout combustion its inert contribution reduces the maximum temperature and therewith pressure. The root causes for this are likely related to the boundary conditions in the cylinder and exhaust port. As mentioned, only a rather short section of the port has been included in the domain (see Figure 4). In addition wall temperatures were not measured. It is quite possible that this causes a wrong pressure drop during the exhaust valve opening phase such that a larger amount of gases leaves the cylinder. Obviously this strongly interacts with the gas density influenced by the heat transfer. Regarding the pressure envelope, we observed that the EGR rate of the LES varied between cycles which contributes to these variations. As mentioned in Section 1, within the experiment this is countered by an active valve control to keep the EGR rate identical within each cycle. Such a control was not possible in the LES since, among other problems, it requires a specific grid generation for each valve phase.

According to this discussion, the following can be concluded based on the comparison of the measured and simulated pressure curves:

• The larger pressure envelope observed in the simulation is likely related to fluctuations of the EGR which is actively countered in the real engine.

• The higher average pressure is likely caused by uncertainties in the boundary conditions with minor contributions from the other effects mentioned. Probably initiated by deviations in the heat flux, the mass fluxes become erroneous and the EGR rate is lowered. As will be explained in Section 3.3 it is reasonable to assume the following mechanism: The EGR in the LES has a higher temperature than in the experiment due to deviations in the heat transfer. By that it takes up a similar volume even if it is of lower total mass. This yields a similar pressure when the intake valves open and accordingly the same mass of fresh air enters. As a result, the total in-cylinder mass is lower in the simulation yielding the slight underestimation of the pressure increase during compression. Furthermore, the EGR rate is lower as observed. As will be detailed in Section 3.2, the higher EGR temperature causes the simulation to ignite earlier than the experiment. Likewise the higher EGR mass limits the absolute temperature and pressure in the experiment.

As mentioned, some of these relations may require the next sections to become more clear. According to the interaction and possibly amplification of effects (heat flux, boundary conditions, EGR rate) it is apparent that a very high accuracy is required for the simulation to match the measurements. Besides unavoidable modeling uncertainties this comprises more information about the boundary conditions being not available.

Summary of Section 3.1

The illustration of a cycle demonstrated how the conditions evolve towards the auto-ignition process. The EGR initialization by the valve timing is followed by the liquid fuel injection and intake phase after which certain similarities of the temperature and the EGR evolve. Towards the end of the compression the fuel diffused such that all regions are within the flammability limit, however being still strongly inhomogeneous. Likewise the EGR and temperature field show strong variations. Accordingly, at the onset of the auto-ignition, inhomogeneities in all three parameters (λ, EGR, T) are found in which the reaction evolves favored at certain positions. This latter process will be detailed in the next sections. The comparison of the simulated and measured pressure curves showed a qualitative reasonable behavior of the LES with a certain deviations of the peak pressure predicted. The analysis conducted in the following is quite general with respect to such a CAI-operation and is unaffected by these uncertainties corresponding to a specific configuration.

3.2 Characterization of the Inhomogeneities and their Relation to the Ignition Behavior

The analysis is conducted in three steps. First, Section 3.2.1 shows the ignition behavior of different temperature and mixing states. Section 3.2.2 then characterizes the inhomogeneities found at a certain cad within the engine. Finally section 3.2.3 then connects these findings to outline the combustion process of a cycle.

3.2.1 First View on the Auto-Ignition under Different Conditions

In this section we analyze the influence of different temperature and mixing states onto the ignition process. The significance of inhomogeneities within this engine is demonstrated in Figure 10. We first consider the pressure shown in the top. The dashed curve represents a single (arbitrary) cycle taken from the simulation. Out of this cycle we selected several spatial locations of different ignition quality as given by the mixture and temperature which are summarized in Table 3. The two cases with more and less favorable conditions are denoted by strong and weak, respectively. With these initial conditions we performed individual one-dimensional simulations, i.e. each simulation corresponds to an adiabatic piston with a homogeneous mixture featuring the same compression ratio as the real engine. One further simulation has been performed with the conditions that would be present if the distribution within the real engine would be homogeneous (denoted by hom). All simulations start at −25 cad where the corresponding initial conditions have been taken from the real engine. As one can see in the top of Figure 10 the pressure then increases equally for all simulations until a very early ignition is observed for the strong conditions. As expected it then follows the homogeneous simulation with a strong delay since its initial temperature is significantly lower. After that, the weak simulation reacts. Its initial temperature is equal to the homogeneous case but it has a less favorable mixing state given by λ. Not only the ignition delay is increased but also the time required for the reaction. This is visible in the second plot in Figure 10 showing the reaction progress variable. The time required from the initial increase until the end of combustion is about twice for the weak case in comparison to the strong one. The real cycle evolves like a smoothened curve in between these conditions.

Finally, the lowest plot in Figure 10 shows the ignition delay derived from the corresponding homogeneous reactor simulations, i.e. when freezing a certain state of T, EGR, λ, and p, an isenthalpic (isobaric and adiabatic) reaction would require this time to react. For the visualization we transferred this delay time into the cad it would require at the 2000 rpm of the engine. According to this derivation, this delay time is not the one of the real engine but provides an idea whether the ignition can already evolve. In the initial state at −25 cad the delay is at least tens of degree, representing obviously insufficient conditions for the engine. The temperature increased by the compression and then reduces this delay and when it approaches about 10 cad the delay time is in a comparable order with the engine relevant combustion range such that an increase of the progress variable can evolve. From this one-dimensional study the following can be concluded: First, the local inhomogeneities found within the real cylinder show significant differences in the ignition behavior. Second, due to these inhomogeneities, the combustion evolves very smooth over a certain cad range within the real engine. Third, without these inhomogeneities, the operation mode would be less reliable. Considering that we used a rather strong cycle and adiabatic conditions for the one-dimensional piston, the ignition strength under homogeneous conditions may be insufficient for the real engine.

3.2.2 Spatial Inhomogeneities at a Given Crank Angle

This first view provided the information about the ignition behavior of certain thermo-chemical states. The conditions used therein have been given as granted without explaining how the combinations of T_u, λ and EGR actually form. Therefore, as the next step in analyzing the process, in this section we want to characterize the spatial inhomogeneities. Figure 11 shows the scalar fields of two different cycles at −25 cad on the left and right, respectively. A similarity is the central vertical line of high λ and corresponding low temperature and EGR. It is visible more or less pronounced in all cycles and geometrically caused by the flow issuing through the inlet valve. However, except for this structure, the spatial distribution in between cycles is very different. Likewise, the occurrence of ignition spots is found at different spatial locations (not shown). However, as mentioned in the previous section the ignition is strongly linked to the scalar field. If conditioned on this latter the cycles show identical behavior. Furthermore, as already indicated, the scalars are not independent from each other which will be quantified in the following for representative cycles.

The left of Figure 12 shows the joint pdf of T and EGR found within a slice at −25 cad where a clear correlation can be observed as it was very clear for all cycles. As indicated it is caused by the high temperature of the EGR as illustrated on the left of Figure 13. Here one can see the strong alignment of the temperature field with the EGR lines, i.e. in most of the regions their gradients are oriented towards the same direction. The EGR lines are colored to provide the information that high EGR levels are associated with a high temperature and vice versa. As expected there is one exception of this correlation: At boundaries the temperature is reduced by the heat flux into the wall (see left and right border) which does not affect the EGR. This region has been removed for the samples entering the pdf in Figure 12 but the influence was not significant.

Next we consider the joint pdf of λ and EGR shown on the right of Figure 12 with the corresponding field shown on the right in Figure 13. The pdf shows two distinct branches. Those were not always that distinct as the one shown but existed for all cycles. Their physical significance is explained on the left of Figure 13. First we have the diagonal branch where the EGR drops with increasing λ. It simply results from mixing where the fresh air entering the intake increases λ and reduces EGR. Accordingly, low EGR values are associated with a high λ as visible in most of the top left region of the slice. One extremum is marked by the circle 1 where the spatially lowest EGR rate and highest λ are found. The second branch shows a reduction of the EGR and constant λ. This is caused by mixing of air and EGR as well. However, here also additionally fuel is entered into the mixture such that λ remains constant while EGR significantly reduces. Such a situation is marked by circle 2 in Figure 13. Specifically in this region, the EGR changes strongly from 50% to 60% while λ is literally constant.

Now we have a good understanding on how certain states form so the next step is to link them to the ignition. Figure 14 shows the ignition delay as derived from isenthalpic simulations. For this given pressure corresponding to the one found at −5 cad it is a function of EGR, λ and T. The dependencies are clear and expected: First, the delay decreases strongly with the temperature due to the exponential link of the reaction to the temperature by the Arrhenius law. Second, the delay increases with EGR since it mostly consists of inert products that reduce the reactivity. Third, the excess in oxidizer causes the delay to increase with λ.

This would suggest that for example a high temperature and a low EGR provide good ignition conditions. However as we know from their correlation in Figure 12 such conditions do not exist in the engine. In Figure 15 the pdf of the states found in the engine is superimposed onto the thermochemical states determining the ignition delay. Regarding this delay, the increasing EGR and T have opposed effects but one can see that the slope is such that their correlated increases evolve with a slope leading to an overall decrease of the ignition delay. The pdf shown here is taken from −10 cad to show realistic delay times. Furthermore, wall regions have not been excluded for this pdf and do show that some regions of exceptionally low temperature exist. Excluding these low values influenced by heat losses; throughout the slope of the pdf the ignition quality varies from about 40 cad to 10 cad for this crank angle. This is a first strong indicator that the combustion process gets determined by this distribution.

To supplement this assumption, we make one further step. The top left of Figure 16 shows again the pdf of T and EGR but with normalized coordinates. Here we introduce the Heating-Coordinate given by the red arrow. It represents a joint coordinate linearizing the dependency of EGR and T. It has a certain error (all deviations from this line) but allows to reduce the three-dimensional state description f(T, EGR, λ) to a two-dimensional f(Heating-Coordinate, λ) one to ease further insight. Accordingly the delay given in Figure 14 reduces to a slice being given in the middle left in Figure 16. Here one can see that the major factor is the preheating but also, that for a given value of this preheating the local air–fuel ratio λ can then have a certain influence.

The top right of Figure 16 shows the joint pdf of the Heating-Coordinate and the ignition delay. They are indeed strongly linked. The two plots below show the average values of λ and the EGR superimposed onto this pdf. These contours show that variations orthogonal to the evolution are associated with variations in λ while the EGR is aligned rather parallel to the Heating-Coordinate which confirms its derivation. So far the analysis only showed a certain ignition likeliness. The proof that these conditions are actually determining the combustion process is given in the bottom left of Figure 16. Since we chose a cad where first slight ignition spots already formed it is possible here to demonstrate that a clear dependency of the reaction progress to the Heating-Coordinate exists.

3.2.3 Process throughout the Full Cad Range to Outline the Course of the Combustion Process

From the analysis of the last section we could see that strong spatial inhomogeneities exist. We showed that a preheating by EGR takes place along a trajectory in the thermo-chemical space which we termed the Heating-Coordinate. Due to its slope it is the main vector determining ignition. A certain scatter orthogonal to this coordinate exists which is mostly caused by variations in λ. The analysis was conducted for a given crank angle and holds for all cycles. Within this section we consider the full evolution of the combustion process throughout a cycle to complete the picture.

First we follow the combustion process as monitored at the five points marked in Figure 17, each of them characterizing one of the typical evolutions found. The variables found at these points are provided in Figure 18 as a function of the cad being detailed in the following according to their number:

• 1.

  The first point is the earliest to ignite. One can see that it has a high initial temperature, a corresponding high EGR and a rather low λ at onset of ignition. Accordingly, its ignition delay is very short. The EGR variation seen here corresponds to the convection of a structure spanning over approx. 3 mm. This is about the lower limit expected from the mesh size and accordingly it is quite possible that higher frequencies exist in reality.

• 2.

  This point is close to point one but ignites independent. It has a comparable evolution but all ignition favoring parameter are just a bit lower and accordingly its reaction process is delayed.

• 3.

  The third point is an example where good conditions are countered by heat losses to the walls. It is positioned within the gap at the cylinder walls and from its temperature and progress variable evolution one can see that the reaction gets slowed down despite an initially high temperature, EGR and low λ.

• 4.

  The fourth point has no wall influence but very poor mixing conditions. Its EGR is still average but λ is high. Accordingly, it has a low temperature and ignites rather late.

• 5.

  Finally this point has poor mixing conditions and is also submitted to heat losses. These conditions represent the last to complete the reaction. However, since the volume-to-surface ratio is low at the TDC, those conditions have a significant contribution to the overall process.

To provide a better view on the spatial distribution Figure 19 shows the temporal evolution found along a line connecting two of these points. First, Figure 19a is monitored along the line connecting the points one and two. As one can see, in between point one (at r = 0 mm) and point two the reaction evolves slower. The main cause is that they are separated by a region with a significantly lower EGR and corresponding initial temperature as illustrated in the left plot. Furthermore, as the plot on the right indicates, also λ is higher and accordingly it requires further compression to sufficiently lower the ignition delay.

Still for all positions along this line, λ was still within the range for a good ignition. This is contrasted in Figure 19b showing the same information for a line connecting the points one and four. As one can see in Figure 17 on about half the way one enters regions of high λ representing poor ignition conditions with the correspondingly increased delay.

As the last part of this section we want to supplement the observations at certain positions by statistics gathered within the whole in-cylinder volume. Therefore we performed a conditioned analysis as exemplified in Figure 20. It shows the pdf of the EGR conditioned on the progress variable, i.e. it shows which EGR values are found for a given progress variable. The first instant shown is −12 cad where the first spots approach the fully reacted state. The evolution shows what we already know from the previous section: only at high EGR values reaction progress is found. At −10 cad reaction is already found also for lower EGR values. Furthermore, the distribution became broader since the ignition criteria are now no longer limited to the best spots exclusively. Within the third plot showing −5 cad a completed reaction progress is found for all kinds of EGR values and broadening continued. Finally at −2 cad the reaction mostly completed where now also the lowest EGR states (50%) react.

This representation of the reaction progress by means of the conditioned pdf is summarized in Figure 21. For a more compact presentation it only shows the mean values as indicated by the white line in Figure 20. For λ an identical evolution as outlined for the EGR is found where the lowest values react first and the behavior converges towards a vertical line with increasing cad. The temperature behaves different due to its physical significance. In contrast to the EGR and λ it increases by the compression and determines the ignition for a given mixing state. Even though the curves show a certain statistical inaccuracy an interesting and somehow expected observation can be made. Looking for example at which unburnt gas temperature T_u certain conditions reach half of the reaction progress (PV = 0.5) the following evolution is visible: Up to approx. −4 cad, the temperature required to reach this state increases. This is likely linked to the dependencies shown in Figure 14, i.e. with increasing cad the ignition takes place at less suited mixing conditions which require a higher temperature to ignite. Within the two latest cad it becomes visible that the reaction gets strongly influenced by the wall being that intermediate states of the reaction progress variable are found for rather low temperatures.

Summary of Section 3.2

First Section 3.2.1 quantified how the three parameters (λ, EGR, T_u) influence the ignition quality by showing how the ignition delay evolved for conditions found within different regions of the cylinder. Section 3.2.2 then derived how these conditions actually form. Based on the correlation of EGR and T it was possible to reduce the three-dimensional composition space to a two-dimensional given by the heating coordinate and λ. It was shown that those are the decisive parameters that predetermine the auto-ignition. Section 3.2.3 then illustrated how the reaction process evolves from this predetermined state which was presented based on individual points as well as on statistical analysis of the whole cylinder.

3.3 Global Behavior and Consecutive Cycles

Up to now we considered a single cycle. The analysis provided us with a good understanding of the process, specifically the distribution of scalar quantities, how they interact and accordingly determine the combustion process. Within this last results section, we take a look at the global behavior and compare different cycles to complete the analysis of this engine.

Plot 1 in Figure 22 shows the ignition time for eleven consecutive cycles. For the illustration in this section we chose the cycles 11 to 20 of the first set of simulations. The first ten cycles have been excluded to remove effects from the initialization which persist for a long time due to the EGR. This ignition time is defined to be the cad where the slope of the spatially averaged reaction progress variable has its maximum. This criterion was found to be appropriate to judge on the ignition strength of a cycle and is consistent with other quantities like the instant of pressure increase when comparing the cycles amongst each other. One can observe a total variation of about 7 cad and a maximum difference of about 4 cad between two consecutive cycles. This is quite significant and worth a deeper look. From our previous analysis we know that a high EGR can produce good ignition conditions but as the first plot in Figure 23 demonstrates there is no clear correlation – at least not to the average EGR of the respective cycle. Still it could be that the EGR is more homogeneous for the cycles that show a late ignition. But this is not the case as given in Figure 24. It shows the joint pdf as already introduced in Figure 12 but for cycles of different ignition strength. As the corresponding pressure curves on the right illustrate we chose a rather weak, average, and strong cycle respectively. One can see, that for a given cad before ignition, indeed the EGR distribution spans a similar range indicating equal inhomogeneities. However, the slope is very different such that the strong cycle reaches high temperatures that do not even exist in the weak one and have a lower probability for the average cycle. Accordingly, the inhomogeneities we analyzed in the last section are important for the combustion process for a given cycle. But when comparing individual cycles there seem to be further determining factors. The cycles apparently have a different average temperature superimposed onto a similar inhomogeneous mixture. Indeed, a clear correlation exists between this average temperature before ignition (T_Ign defined to be at −25 cad where no combustion yet occurred) and the ignition time as given by their correlation in plot 2 of Figure 23.

Accordingly, the question is what this increased average temperature caused. To understand this, we first want to see from where it evolved simply from an observing point of view and then follow it up by an explanation. Since the engine has no significant thermal variations and fuel gets injected at identical temperatures for each cycle, it is clear that the difference is in the history of the gases. Therefore we backtrack it through its state at the Inlet-Valve-Close and Open (IVC and IVO), the preceding’s cycles Exhaust Valve Close and Open (EVC and EVO) and its state found after the combustion of the preceding cycle (aComb) defined to be at 80 cad. Starting from the third plot in Figure 23, as expected T_Ign translates from T_IVC which in turn correlates with T_IVO (plot 4). This latter is not necessarily expected but will be explained later. Going further back to plot 5, the temperature T_IVO correlates with the one found after the exhaust valve close of the preceding cycle. As plot 6 indicates also the gases after and before the exhaust phase correlate but with certain deviations. The latter in turn is perfectly correlated to the state after combustion (plot 7). Accordingly, we can bridge these subsequent steps in plot 8 showing the temperature before ignition as a function of the one found after combustion of the preceding cycle. This correlation has a certain scatter arising from the sum of the intermediate processes but still shows a visible trend. Via the correlation in plot 2 we can now directly link these conditions after combustion to the cad of ignition in plot 9.

So we observed that a cycle leaving exhaust gases with a high temperature very likely lowers the ignition delay of its successor. The scatter indicates that there are other parameters but a trend definitely exists. Shedding further light onto this requires to consider many parameters. For example, from the equation of state one guess could be that a high temperature is caused by a large trapped mass, but those are not correlated (not shown). Also the peak pressure is not a good indicator since it does not strongly depend on the overall energetic state but rather on the cad where it occurs. Also, a bit surprising, as plot 10 in Figure 23 indicates, the temperature after combustion is not correlated with the one before combustion. This could be suggested when assuming a given energy release that causes an offset to its initial state. The root causes are not straight-forward and some correlations appear counterintuitive at a first glance. We found it best to outline the mechanisms along some consecutive cycles to provide a good understanding.

We consider the cycles starting with cycle number 15. As given in Figure 22, cycle 17 has a very late ignition and from our above analysis cycle 16 then should have a low T_aComb which holds true (plot 2 in Figure 22), but we will start at cycle 15. As one can see in Figure 22, cycles 15 and 16 have the same mass, the same temperature before ignition and corresponding ignition time (plots 3, 4 and 1, respectively). However, the temperature after combustion of cycle 16 is lower by 44 K being a visible drop in this regard. The difference to cycle 15 is that in cycle 16, the same mass constitutes of a lower EGR mass and a larger amount of air issued through the inlet valve. The causes for this latter will become clear when we consider the processes following the combustion of cycle 16 and first take the low T_aComb of it as granted.

Given the low T_aComb of cycle 16, the pressure is correspondingly low and as a result a lower amount of mass will leave the cylinder when the exhaust valves open. Hence, as one can see in plot 6 of Figure 22, cycle 17 has a high EGR mass. It is important to keep in mind that this exhaust gas is comparatively cold as given by its previous cycle. It is visible in plot 7 of Figure 22 that cycle 17 has a rather low EGR temperature and indeed, as plot 11 in Figure 23 indicates, the EGR mass and its temperature are generally negatively correlated. This is a first property to notice. From our local analysis we know that a locally high EGR favors ignition. However, when considering a globally high EGR it is partially countered by having a low temperature. When we approach the IVO with this EGR mass, two countering effects come into play. The in-cylinder mass is rather high but it is also cold such that the pressure at this instant does not correlate with the in-cylinder mass as indicated in plot 12 of Figure 23. Indeed, as visible in this plot those pressures are almost identical for cycles 16 and 17 yielding the same mass of air issuing during the intake phase (see plot 5 in Figure 22). This means, that cycle 17 has a high EGR mass and the same mass of air, i.e. a larger total mass. Indeed, plot 13 of Figure 23 confirms that those are often correlated. This now also explains why the temperatures before and after the inlet valve phase are correlated (plot 4 in Figure 23): Given that the same amount of air enters in the cycles 16 and 17 with the same air-temperature, the state is obviously determined be the exhaust gas temperatures.

Looking at the properties of cycle 17 in Figure 22 one can see that this history results in one of the latest reactions leaving a high temperature after combustion. At this, the strong ignition delay is caused by the low exhaust gas temperature of its preceding cycle giving little potential for hot spots as discussed in the previous section. However, it also has a high temperature after combustion due to its large total mass yielding high pressure and temperature conditions for a given cylinder volume. Still, we want to remind that it is not the total mass alone but also its composition as outlined for this cycle. According to this high T_aComb, cycle 18 has a low ignition delay again. The procedure outlined here for cycles 15–17 repeats at cycles 18–20. It would require a significant amount of more cycles computed to gather statistics of this processes spanning over several cycles. However, our explanation represents a logical basis to understand such large cycle-to-cycle variations. Also for some of the smaller variations in between consecutive cycles observed this explanation holds true but the physical significance is less dominant such that other parameters like local inhomogeneities can interfere.

Conclusion

This work applied LES to a CAI engine using a five-dimensional pre-tabulation strategy to predict the auto-ignition process. The individual physical mechanisms evolving throughout the cycle have been outlined to provide a basic understanding of the overall process. A comparison to measured pressure curves showed that the simulation is able to reproduce the basic characteristics of the process with certain deviations regarding the peak pressure and envelope width. Explanations for those have been provided being mostly related to uncertainties in the boundary conditions and their unfortunate interaction with the EGR rate. We outlined the dependency of the ignition on scalar quantities being the temperature, air–fuel ratio, EGR, and pressure. The effect has been quantified by means of one-dimensional piston simulations for typical combinations of those as found in the actual engine.

A statistical analysis was performed to characterize the strong spatial inhomogeneities within the cylinder. It revealed a strong correlation between the EGR and the temperature as well as two separate branches along which the probability of EGR and the air–fuel ratio evolved. These evolutions have been explained by outlining and visualizing the respective processes. We showed that a preheating by EGR takes place along a trajectory in the thermo-chemical space which we termed the Heating-Coordinate. By projecting this heating coordinate onto the thermo-chemical space of the ignition table we showed that its slope is the main vector in producing conditions able to ignite. Accordingly the reaction process was strongly correlated with this composition coordinate. A certain scatter orthogonal to this coordinate exists which is mostly caused by variations in the air–fuel ratio. It was demonstrated that the combustion process is pre-determined by these conditions which we supplemented by tracking individual points and lines through the reaction phase.

Finally, the analysis of consecutive cycles revealed a strong influence caused by the EGR. It was shown that the cycles show non-negligible differences in the ignition delay. An explanation was provided for a variation that spans over 3–4 cycles supported by correlations of the decisive quantities. This study was somehow limited by the sample number but still enabled to shed light onto the mechanisms evolving.

Acknowledgments

The authors gratefully acknowledge financial support by the FVV-DFG and FAUDI foundation. The computations have been carried out on the Lichtenberg High Performance Computer of TU Darmstadt (HHLR).

References

---

All Tables

All Figures