# A Two-Zone Multigrid Model for SI Engine Combustion Simulation Using Detailed Chemistry

Authors: Hai-Wen Ge; Harmit Juneja; Yu Shi; Shiyou Yang; Rolf D. Reitz
Journal: Journal of Combustion
 (2010)

---

## Abstract

An efficient multigrid (MG) model was implemented for spark-ignited (SI) engine combustion modeling using detailed chemistry. The model is designed to be coupled with a level-set-G-equation model for flame propagation (GAMUT combustion model) for highly efficient engine simulation. The model was explored for a gasoline direct-injection SI engine with knocking combustion. The numerical results using the MG model were compared with the results of the original GAMUT combustion model. A simpler one-zone MG model was found to be unable to reproduce the results of the original GAMUT model. However, a two-zone MG model, which treats the burned and unburned regions separately, was found to provide much better accuracy and efficiency than the one-zone MG model. Without loss in accuracy, an order of magnitude speedup was achieved in terms of CPU and wall times. To reproduce the results of the original GAMUT combustion model, either a low searching level or a procedure to exclude high-temperature computational cells from the grouping should be applied to the unburned region, which was found to be more sensitive to the combustion model details.


---

## Body

## 1. Introduction

CFD-based engine development is widely used in the automotive industry due to its much lower cost compared to experiment-based engine development [1]. Additionally, computer simulation provides more insight about the complex engine combustion process which is inaccessible or costly using current measurement techniques. Improved accuracy and efficiency in numerical models and rapidly advancing computer capabilities makes CFD-based engine development more feasible than ever. Thus, traditional empirical hardware-based optimization is being replaced with CFD-based optimization. As a result, the development time of a vehicle in the automotive industry has been reduced from average 60 months twenty years ago to 18 months today [2].Detailed chemistry is necessary to reproduce kinetics-controlled processes, for instance, the processes of auto-ignition and pollutant emissions [3]. Understanding such processes is essential for improving engine performance, especially when considering strict emissions regulations. Advanced engine design concepts, including homogeneous charge compression ignition (HCCI), premixed charge compression ignition (PCCI), and modulated-kinetics (MK), are all chemical kinetics-controlled, and simulation of these combustion processes requires detailed chemistry. Detailed chemistry is often implemented into CFD software through coupling with the CHEMKIN package. For example, Kong and Reitz [4] developed a KIVA-CHEMKIN model which implemented the CHEMKIN II solver [5] into the KIVA code [6]. The KIVA-CHEMKIN model assumes that each cell is a perfectly stirred reactor (PSR) and that species conversion rates can be modeled using the CHEMKIN solver. The KIVA-CHEMKIN model has been proven to be a very reliable tool for diesel engine simulation and thus widely used in industry and academia. However, even with a relatively simple detailed mechanism, such as the ERC reduced n-heptane mechanism (34 species, 74 reactions) [7], computational time is increased by 1 to 2 order(s) of magnitude compared with conventional simplified combustion models, such as the Characteristic Time Combustion (CTC) model [8]. Therefore, efficiency is of vital inherent when applying detailed chemistry to engine simulations, especially when used in engine optimization studies. Depending on the total number of design parameters considered in the optimization, hundreds or thousands of individual simulations are required for one optimization task [9, 10]. Thus, many efforts have been devoted to reduce the computational time spent in solving the high-dimensional and stiff ordinary differential equations (ODEs) that arise from the inclusion of detailed chemistry. Strategies for this purpose can be categorized into the following classes:(1)
reduce the number of species and reactions in the mechanism: for example, [3];(2)
dynamically reduce the size of the Jacobian matrix in the CHEMKIN solver by lumping or assuming equilibrium species and/or reactions that have little impact on the global reaction system:  for example, QSSA, ILDM, CSP, partial equilibrium [3, 11], and dynamic adaptive chemistry (DAC) scheme [12];(3)
reduce the frequency of calling the CHEMKIN solver: for example, multizone models [13], dynamic cell clustering [14], and adaptive multigrid chemistry (AMC) [15];(4)
solve the ODEs using alternative faster numerical schemes: for example, [16].In the present paper, a one-zone and a two-zone multigrid (MG) models were implemented into the level-set-G-equation model for SI engine simulation. The models were validated using comparisons with the full chemistry solver for a gasoline direct-injection (GDI) engine case. The computational efficiency and accuracy of the models are discussed, and it is concluded that speedup techniques that have previously been developed for diesel engine simulations are successful also for modeling spark-ignition engines.


## 2. Mathematical Models

### 2.1. Multigrid Model

The multigrid (MG) model [13] is designed to reduce the frequency of calling the CHEMKIN solver, which is the most time-consuming part of simulation when detailed chemistry is used. The MG model groups thermodynamically similar cells via a mapping procedure into one “grouped cell” for the CHEMKIN calculation, followed by a procedure to redistribute the group information back onto the individual cells using remapping. The grouping region (searching level) is either fixed [13] or determined adaptively [15] accounting for the local inhomogeneities in the thermodynamic conditions to deal with the nonuniform spatial distribution of thermochemistry properties. The concept of an adaptive neighbor search is further illustrated in Figure 1 using a 2D schematic mesh. It can be seen that the first level search only covers the four adjacent cells (six if 3D mesh) of the reference cell. If the in-cylinder temperature inhomogeneity is below prespecified values, the search can then be expanded to a second level or higher. In this study, the maximum search level was limited to four, where a maximum of 129 cells can be reached using the fourth level search in a 3D block structural mesh. The fifth searching level is a global search, which means that grouping is not limited to neighbor cells and all the computational cells are considered for grouping. The model offers in total seven options: no MG, fixed first level, fixed second level; fixed third level, fixed fourth level, fixed fifth level, and finally the adaptive level. When adaptive searching is used, the searching level at each time step is determined considering the inhomogeneity of all the computational cells. The grouping criteria considered in the present work include the cell temperature and progress equivalence ratio [13, 15]. In the present work, a difference in temperature of 20 K and a difference in equivalence ratio of 10% were used to determine if cells are thermodynamically “similar” [15]. Temperature inhomogeneities of 400 K, 350 K, 300 K, 250 K, and 50 K were used to determine the first through fifth (global) searching levels, respectively [15].Figure 1
Schematic of adaptive multigrid grouping procedure.After the calculation of chemical reactions of all the grouped cells, the mass fractions of each species were remapped back onto the original cells. A gradient-preserving remapping method [13, 15] was used to preserve temperature and species composition gradients. Firstly, before the mapping procedure, a new quantity, ch, is defined using the number of C and H atoms of all participating species, except the combustion products, CO2 and H2O, where(1)ch=2C-CO2#+H-H2O#2.The ch number is an indicator of reactivity of a mixture and defined based on the equivalent number of O atoms that are required to form final products, CO2 and H2O. Thus, CO2 and H2O in the mixture are excluded in the definition of the ch number. The ch number of a group is the sum of the ch number of all cells in that group. After the chemistry calculation, all species, except CO2, H2O, O2, and N2, are assigned back to the group’s cells based on ch. In this case, the mass of species k in an individual cell is obtained from the ratio:(2)mk,cell=mk,groupchcellchgroup.In this way the mass of each species in the group is also conserved. Consequently, some cells can have more or fewer C or H atoms than before the mapping process, and thus the rest of the cells in the group need to be adjusted to maintain the total number of C and H atoms in that group. After that O2 is distributed to maintain the total number of O atoms in each cell and, finally, adjustment of N2 is used to conserve the mass of each cell if necessary. Since after the remapping process the mass fraction of each species in each cell is known, the change of the specific internal energy of each cell can be obtained from the difference between the internal energy of formation of the species present in the cell before the grouping process and that after the remapping process. The cell temperature can be computed from the updated specific internal energy and the mass fraction of species.As mentioned earlier, the multigrid chemistry model has been successfully used for DI engine simulations [15] and for optimization of diesel engines operating under different engine loads [1, 17]. It has also been successfully coupled with the DAC scheme and parallelized using a load balancing scheme [18]. It was found in these studies that significant amounts of computational time were saved with the adaptive multigrid chemistry model.


### 2.2. Level-Set-G-Equation Model

The present level-setG-equation model is a turbulent combustion model which tracks flame propagation. The flame front propagation is represented by an iso-scale surface (usually the scalar G=0), and the flame propagation process is modeled by tracking the location of this surface.  The transport equation of the iso-G surface is [19](3)ρ∂G∂t+ρU⃗·∇G=ρsL|∇G|,
and it is solved using a level-set method. The second term, on the left-hand side, is the convective term. The term on the right-hand side is the propagation term which accounts for the chemical properties of the mixture. Stretch effects on the flame are reflected in the curvature of the G surface. Therefore, turbulence-combustion interactions are also described by the G-equation. Considering the moving meshes used in engine codes, the modeled equations for the Favre-averaged G are written as [20](4)∂G̃∂t+(U⃗̃bf-U⃗vertex)·∇G̃=ρ̅uρ̅sT0|∇G̃|-DTκ̃M|∇G̃|.
The computational domain is separated into three parts: the region ahead of the flame front (unburned region, G<0), the region on the flame front (G=0), and the regions behind the flame front (burned region, G>0). To consider autoignition ahead of the flame front (e.g., leading to “knocking” combustion) and stratification behind the flame front, the cells ahead and behind the flame front are treated as PSRs and their thermochemical states are computed using the CHEMKIN solver. With this consideration, the model is then capable of modeling of both turbulent premixed and nonpremixed combustion. Thus, the model is also referred to as the G-equation for All Mixtures, a Universal Turbulent (GAMUT) combustion model [20]. Figure 2 shows the schematic flame structure of a partially premixed flame and the corresponding modeling strategies for the different regimes. Computational cells in the unburned region with temperatures greater than 1100 K and a fraction of combustion products greater than certain amount are regarded as autoignition sites. As soon as an ignition kernel is formed by either the spark plug model or autoignition, the flame propagation model is activated.Figure 2
GAMUT combustion model for partially premixed combustion.In the present, the discrete particle ignition kernel (DPIK) model by Fan and Reitz [21] was used for modeling of spark ignition. Once the characteristic length of the ignition kernel region is larger than the averaged turbulent integral length scale, a G surface is initialized. The cells on the flame front are assumed to be in chemical equilibrium; that is, complete combustion is assumed to occur across the flame. The change of composition in a computational cell containing the flame front is determined from the swept volume of the propagating flame:(5)dYidt=(Yiu-Yib)AVsT0,
where A and V are the flame surface area and volume of the computational cell, respectively. Yi is the mass fraction of the ith species. The superscripts “u” and “b” refer to the values in the unburned (G<0) and burned (G>0) regions, respectively. The turbulent burning velocity in (4) and (5) is related to the laminar burning velocity and turbulence intensity as well as the Damköhler number [19]:(6)st0=sl0+v′{-a4b322b1Da+[(a4b322b1Da)2+a4b32Da]1/2}.Here,sl0 is the laminar burning velocity, and v' is the turbulence intensity, equal to (2/3)k, where k is the turbulent kinetic energy. The model constants a4=0.78, b1=2.0, and b3=1.0 are derived from the turbulence model [19].  The Damköhler number is defined as(7)Da=sl0ltv'lF,
where lt is the turbulence integral length scale and lF is the thickness of the flame:(8)lF=(λ/Cp)0(ρsl0)u.The heat conductivity,λ, and heat capacity, Cp, are evaluated at the inner layer temperature. The laminar burning velocity, sl0, is determined from an experimental correlation, for example, that of Metghalchi and Keck [22] is widely used:(9)sl0=sl,ref0(TuTu,ref)γ(pupu,ref)β(1-2.1Ydil),
where Tu,ref=298K, pu,ref=1atm, and Ydil is the mass fraction of diluents present in the fuel-air mixture. The subscript “u” indicates the unburned gas, and the reference laminar burning velocity is(10)sl,ref0=BM+B2(ϕ2-ϕM)2
with BM=26.32 cm/s, B2=-84.72 cm/s, and ϕM=1.13. The exponents γ and β are given as(11)γ=2.18-0.8(ϕ-1),β=-0.16+0.22(ϕ-1),
where ϕ is the equivalence ratio.The flame surface from theG-equation describes the mean position of the flame brush, which has a thickness on the order of a computational cell.  The chemical reaction rate for each species is calculated in each time step as(12)∂Yk∂t=(Yku-Ykb)AVst0,
where Yku and Ykb are the mass fractions of species i in the unburnt and burnt gases on either side of the flame surface. Ykb are determined by assuming chemical equilibrium in the burnt gases. V is the volume of the computational cell and A is the total area of the flame front in the cell.The GAMUT combustion model has been developed and implemented into the KIVA-3V code by Tan and Reitz [20, 23]. It has been extensively used to simulate combustion in spark-ignited gasoline engines, in which the combustion process is a stratified premixed combustion [20, 23–28]. It has also successfully been applied to simulate diesel/natural-gas dual-fuel CI engines [29, 30] and conventional diesel engines [30–33]. Yang and Reitz [28] have made improvements to the combustion submodels for modeling SI combustion with the level set G-equation method and detailed chemical kinetics, including a transport equation residual model, modeling flame front quenching in highly stratified mixtures, and the inclusion of a recently developed PRF mechanism [34].


### 2.3. Two-Zone Multigrid Model

When the GAMUT combustion model is used, chemical reactions in the cells ahead and behind the flame front are computed using the CHEMKIN solver. Thus, similar to the KIVA-CHEMKIN model, the GAMUT combustion model spends the most CPU time in the CHEMKIN solver. Since additional CPU time is needed to solve theG-equation, the GAMUT combustion model usually requires more CPU time than the KIVA-CHEMKIN model. This motivated integration of the MG model into the level-set-G-equation model for higher computational efficiency. If the original MG model is directly applied to the GAMUT combustion model, the process of searching and grouping thermodynamically similar cells covers all the computational cells, the “one-zone MG model”. However, the burnt and unburned regions have very different thermochemical conditions in both temperature and composition. The burnt region usually has much higher temperatures, more combustion products, and less reactants than the unburned region. These two regions are separated by the flame front, which is numerically represented by the iso-surface of G=0. Thus, it is logical to apply the MG model separately to the burnt and unburned regions, the “two-zone MG model”. The mapping and remapping procedures in the burnt and unburned regions are independent, and thus different searching levels can be used for the two regions according to their local conditions. Both the one-zone and two-zone MG models were implemented into the KIVA3v code for comparison.


### 2.4. Adaptive Load-Balancing Algorithm for Parallelization

In the present work a load balancing scheme was also used to parallelize the CHEMKIN computations, which are the most time-consuming part of the code. The parallelization is based on the message passing interface (MPI). Before the CHEMKIN solver is called, all cells that need to be computed using CHEMKIN are listed. At the initial step, those cells are evenly assigned to the available processors. Then, the speed of each processor can be defined as the ratio of the number of cells computed by this processor and the CPU time taken in the previous time step:(13)Pin+1=Nckcell,intin,
where the superscript n denotes the time step and subscript i indicates the ith processor. Starting from the next time step, the cells are assigned based on the speed of the processor Pin+1:(14)Nckcell,in+1=Pin+1∑i=1NcpuPin+1Nck,
in which Nck is the total number of cells required in the CHEMKIN computation. In this way, a dynamic balance between the processors is achieved. To smooth the process, a relaxation factor, α, is introduced in (12) as (15)Pin+1=αPin+(1-α)Nckcell,intin.This algorithm has been validated on diesel engine simulations and has proven to be very efficient [18, 34].Other models used in the KIVA code are described in detail in [28]. A recently developed PRF mechanism [34], which consists of 41 species and 130 elementary reactions, was used in the present paper.


## 2.1. Multigrid Model

The multigrid (MG) model [13] is designed to reduce the frequency of calling the CHEMKIN solver, which is the most time-consuming part of simulation when detailed chemistry is used. The MG model groups thermodynamically similar cells via a mapping procedure into one “grouped cell” for the CHEMKIN calculation, followed by a procedure to redistribute the group information back onto the individual cells using remapping. The grouping region (searching level) is either fixed [13] or determined adaptively [15] accounting for the local inhomogeneities in the thermodynamic conditions to deal with the nonuniform spatial distribution of thermochemistry properties. The concept of an adaptive neighbor search is further illustrated in Figure 1 using a 2D schematic mesh. It can be seen that the first level search only covers the four adjacent cells (six if 3D mesh) of the reference cell. If the in-cylinder temperature inhomogeneity is below prespecified values, the search can then be expanded to a second level or higher. In this study, the maximum search level was limited to four, where a maximum of 129 cells can be reached using the fourth level search in a 3D block structural mesh. The fifth searching level is a global search, which means that grouping is not limited to neighbor cells and all the computational cells are considered for grouping. The model offers in total seven options: no MG, fixed first level, fixed second level; fixed third level, fixed fourth level, fixed fifth level, and finally the adaptive level. When adaptive searching is used, the searching level at each time step is determined considering the inhomogeneity of all the computational cells. The grouping criteria considered in the present work include the cell temperature and progress equivalence ratio [13, 15]. In the present work, a difference in temperature of 20 K and a difference in equivalence ratio of 10% were used to determine if cells are thermodynamically “similar” [15]. Temperature inhomogeneities of 400 K, 350 K, 300 K, 250 K, and 50 K were used to determine the first through fifth (global) searching levels, respectively [15].Figure 1
Schematic of adaptive multigrid grouping procedure.After the calculation of chemical reactions of all the grouped cells, the mass fractions of each species were remapped back onto the original cells. A gradient-preserving remapping method [13, 15] was used to preserve temperature and species composition gradients. Firstly, before the mapping procedure, a new quantity, ch, is defined using the number of C and H atoms of all participating species, except the combustion products, CO2 and H2O, where(1)ch=2C-CO2#+H-H2O#2.The ch number is an indicator of reactivity of a mixture and defined based on the equivalent number of O atoms that are required to form final products, CO2 and H2O. Thus, CO2 and H2O in the mixture are excluded in the definition of the ch number. The ch number of a group is the sum of the ch number of all cells in that group. After the chemistry calculation, all species, except CO2, H2O, O2, and N2, are assigned back to the group’s cells based on ch. In this case, the mass of species k in an individual cell is obtained from the ratio:(2)mk,cell=mk,groupchcellchgroup.In this way the mass of each species in the group is also conserved. Consequently, some cells can have more or fewer C or H atoms than before the mapping process, and thus the rest of the cells in the group need to be adjusted to maintain the total number of C and H atoms in that group. After that O2 is distributed to maintain the total number of O atoms in each cell and, finally, adjustment of N2 is used to conserve the mass of each cell if necessary. Since after the remapping process the mass fraction of each species in each cell is known, the change of the specific internal energy of each cell can be obtained from the difference between the internal energy of formation of the species present in the cell before the grouping process and that after the remapping process. The cell temperature can be computed from the updated specific internal energy and the mass fraction of species.As mentioned earlier, the multigrid chemistry model has been successfully used for DI engine simulations [15] and for optimization of diesel engines operating under different engine loads [1, 17]. It has also been successfully coupled with the DAC scheme and parallelized using a load balancing scheme [18]. It was found in these studies that significant amounts of computational time were saved with the adaptive multigrid chemistry model.


## 2.2. Level-Set-G-Equation Model

The present level-setG-equation model is a turbulent combustion model which tracks flame propagation. The flame front propagation is represented by an iso-scale surface (usually the scalar G=0), and the flame propagation process is modeled by tracking the location of this surface.  The transport equation of the iso-G surface is [19](3)ρ∂G∂t+ρU⃗·∇G=ρsL|∇G|,
and it is solved using a level-set method. The second term, on the left-hand side, is the convective term. The term on the right-hand side is the propagation term which accounts for the chemical properties of the mixture. Stretch effects on the flame are reflected in the curvature of the G surface. Therefore, turbulence-combustion interactions are also described by the G-equation. Considering the moving meshes used in engine codes, the modeled equations for the Favre-averaged G are written as [20](4)∂G̃∂t+(U⃗̃bf-U⃗vertex)·∇G̃=ρ̅uρ̅sT0|∇G̃|-DTκ̃M|∇G̃|.
The computational domain is separated into three parts: the region ahead of the flame front (unburned region, G<0), the region on the flame front (G=0), and the regions behind the flame front (burned region, G>0). To consider autoignition ahead of the flame front (e.g., leading to “knocking” combustion) and stratification behind the flame front, the cells ahead and behind the flame front are treated as PSRs and their thermochemical states are computed using the CHEMKIN solver. With this consideration, the model is then capable of modeling of both turbulent premixed and nonpremixed combustion. Thus, the model is also referred to as the G-equation for All Mixtures, a Universal Turbulent (GAMUT) combustion model [20]. Figure 2 shows the schematic flame structure of a partially premixed flame and the corresponding modeling strategies for the different regimes. Computational cells in the unburned region with temperatures greater than 1100 K and a fraction of combustion products greater than certain amount are regarded as autoignition sites. As soon as an ignition kernel is formed by either the spark plug model or autoignition, the flame propagation model is activated.Figure 2
GAMUT combustion model for partially premixed combustion.In the present, the discrete particle ignition kernel (DPIK) model by Fan and Reitz [21] was used for modeling of spark ignition. Once the characteristic length of the ignition kernel region is larger than the averaged turbulent integral length scale, a G surface is initialized. The cells on the flame front are assumed to be in chemical equilibrium; that is, complete combustion is assumed to occur across the flame. The change of composition in a computational cell containing the flame front is determined from the swept volume of the propagating flame:(5)dYidt=(Yiu-Yib)AVsT0,
where A and V are the flame surface area and volume of the computational cell, respectively. Yi is the mass fraction of the ith species. The superscripts “u” and “b” refer to the values in the unburned (G<0) and burned (G>0) regions, respectively. The turbulent burning velocity in (4) and (5) is related to the laminar burning velocity and turbulence intensity as well as the Damköhler number [19]:(6)st0=sl0+v′{-a4b322b1Da+[(a4b322b1Da)2+a4b32Da]1/2}.Here,sl0 is the laminar burning velocity, and v' is the turbulence intensity, equal to (2/3)k, where k is the turbulent kinetic energy. The model constants a4=0.78, b1=2.0, and b3=1.0 are derived from the turbulence model [19].  The Damköhler number is defined as(7)Da=sl0ltv'lF,
where lt is the turbulence integral length scale and lF is the thickness of the flame:(8)lF=(λ/Cp)0(ρsl0)u.The heat conductivity,λ, and heat capacity, Cp, are evaluated at the inner layer temperature. The laminar burning velocity, sl0, is determined from an experimental correlation, for example, that of Metghalchi and Keck [22] is widely used:(9)sl0=sl,ref0(TuTu,ref)γ(pupu,ref)β(1-2.1Ydil),
where Tu,ref=298K, pu,ref=1atm, and Ydil is the mass fraction of diluents present in the fuel-air mixture. The subscript “u” indicates the unburned gas, and the reference laminar burning velocity is(10)sl,ref0=BM+B2(ϕ2-ϕM)2
with BM=26.32 cm/s, B2=-84.72 cm/s, and ϕM=1.13. The exponents γ and β are given as(11)γ=2.18-0.8(ϕ-1),β=-0.16+0.22(ϕ-1),
where ϕ is the equivalence ratio.The flame surface from theG-equation describes the mean position of the flame brush, which has a thickness on the order of a computational cell.  The chemical reaction rate for each species is calculated in each time step as(12)∂Yk∂t=(Yku-Ykb)AVst0,
where Yku and Ykb are the mass fractions of species i in the unburnt and burnt gases on either side of the flame surface. Ykb are determined by assuming chemical equilibrium in the burnt gases. V is the volume of the computational cell and A is the total area of the flame front in the cell.The GAMUT combustion model has been developed and implemented into the KIVA-3V code by Tan and Reitz [20, 23]. It has been extensively used to simulate combustion in spark-ignited gasoline engines, in which the combustion process is a stratified premixed combustion [20, 23–28]. It has also successfully been applied to simulate diesel/natural-gas dual-fuel CI engines [29, 30] and conventional diesel engines [30–33]. Yang and Reitz [28] have made improvements to the combustion submodels for modeling SI combustion with the level set G-equation method and detailed chemical kinetics, including a transport equation residual model, modeling flame front quenching in highly stratified mixtures, and the inclusion of a recently developed PRF mechanism [34].


## 2.3. Two-Zone Multigrid Model

When the GAMUT combustion model is used, chemical reactions in the cells ahead and behind the flame front are computed using the CHEMKIN solver. Thus, similar to the KIVA-CHEMKIN model, the GAMUT combustion model spends the most CPU time in the CHEMKIN solver. Since additional CPU time is needed to solve theG-equation, the GAMUT combustion model usually requires more CPU time than the KIVA-CHEMKIN model. This motivated integration of the MG model into the level-set-G-equation model for higher computational efficiency. If the original MG model is directly applied to the GAMUT combustion model, the process of searching and grouping thermodynamically similar cells covers all the computational cells, the “one-zone MG model”. However, the burnt and unburned regions have very different thermochemical conditions in both temperature and composition. The burnt region usually has much higher temperatures, more combustion products, and less reactants than the unburned region. These two regions are separated by the flame front, which is numerically represented by the iso-surface of G=0. Thus, it is logical to apply the MG model separately to the burnt and unburned regions, the “two-zone MG model”. The mapping and remapping procedures in the burnt and unburned regions are independent, and thus different searching levels can be used for the two regions according to their local conditions. Both the one-zone and two-zone MG models were implemented into the KIVA3v code for comparison.


## 2.4. Adaptive Load-Balancing Algorithm for Parallelization

In the present work a load balancing scheme was also used to parallelize the CHEMKIN computations, which are the most time-consuming part of the code. The parallelization is based on the message passing interface (MPI). Before the CHEMKIN solver is called, all cells that need to be computed using CHEMKIN are listed. At the initial step, those cells are evenly assigned to the available processors. Then, the speed of each processor can be defined as the ratio of the number of cells computed by this processor and the CPU time taken in the previous time step:(13)Pin+1=Nckcell,intin,
where the superscript n denotes the time step and subscript i indicates the ith processor. Starting from the next time step, the cells are assigned based on the speed of the processor Pin+1:(14)Nckcell,in+1=Pin+1∑i=1NcpuPin+1Nck,
in which Nck is the total number of cells required in the CHEMKIN computation. In this way, a dynamic balance between the processors is achieved. To smooth the process, a relaxation factor, α, is introduced in (12) as (15)Pin+1=αPin+(1-α)Nckcell,intin.This algorithm has been validated on diesel engine simulations and has proven to be very efficient [18, 34].Other models used in the KIVA code are described in detail in [28]. A recently developed PRF mechanism [34], which consists of 41 species and 130 elementary reactions, was used in the present paper.


## 3. Results and Discussion

The models were applied to an experimental single cylinder homogeneous gasoline direct injection (GDI) engine [35]. This engine features a wedge type combustion chamber and a shallow piston bowl. A full mesh was created using WERC’s in-house grid generation tool capable of producing structured grids. Figure 3 shows a snapshot of the computational mesh with a maximum grid size of 238,000 cells. The direct-injection (DI) system employs an outboard side spark location and an intake-side-mounted multihole injector. The specifications of the engine and the selected operating conditions are listed in Table 1. Experiments showed that autoignition occurs in this engine case [35]. Thus, this case is a critical case for validation of the MG and GAMUT combustion models. The GAMUT combustion model has been validated in [35] by comparing with available experimental data, including pressure trace and knock intensity. Excellent agreement with the measured data was obtained. The present paper focuses on the efficiency and accuracy of the new MG model. Thus, the numerical results from the original GAMUT combustion model are taken as a reference for model validation. The calculation in the present paper starts at -6.3 deg ATDC and ends at 70 deg ATDC. At this moment, all the liquid fuel has vaporized. All cases were calculated using the parallel code on 4 nodes. Each node is equipped with an Intel Pentium 3.0 G Hz CPU and 2.0 GB RAM memory.Table 1
Engine specifications and selected operating conditions.
Bore×Stroke (mm)99.6×92.0Connecting Rod length (mm)154.9Compression Ratio11 : 1Engine Speed (rev/min)1600MAPwide-open throttle (WOT)Air-Fuel Ratio12 : 1Injector Type6-hole InjectorFuel Injected (mg)60Start of Injection (deg ATDC)-243Spark Timing (deg ATDC)-6.0Computational mesh for the single cylinder homogeneous GDI engine.
(a)(b)The GDI case was first simulated using the one-zone MG model. Figure4 shows the in-cylinder pressure and heat release rate (HRR) of the GDI engine predicted using the different models, including the original GAMUT combustion model (without MG model) and the GAMUT combustion model with the MG model. Different searching levels were tested: fixed three level (MG3), fixed four level (MG4), fixed global level (MG-global), and the adaptive search (MG-adaptive). The corresponding computational times of the different solution methods are listed in Table 2. The speedup data in the table refers to the CPU time of the original combustion model and wall clock time of each case. It can be seen that the parallel code works well. Threefold or more speedup was achieved for all cases, except for the fixed global adaptive search case, in which the CPU time spent on the chemistry solver is very low and thus the efficiency of the parallel scheme is low. Combined with the MG model, the computational time was reduced by at least one order of magnitude. With the fixed global search, the total computational time was reduced from 173 hours to 5 hours. With other searching levels, around a 10-fold speedup was achieved. The computational cost of the adaptive searching level is similar to the one of the third searching level. However, the predicted pressure and HRR using the one-zone MG model are not good. As seen in Figure 4, none of them can accurately reproduce the HRR curve. All of them predicted an earlier main combustion event than the original GAMUT combustion model.Table 2
Computational time of one-zone MG model with different searching levels (CA=-6.3 to 70).
Casew/o MGMG3MG4MG-globalMG-adaptiveCPU time [hrs]173.355.345.65.552.6Wall time [hrs]48.618.416.55.517.9Speedup [-]3.69.410.531.59.7Pressure and HRR of SI engine simulated using one-zone MG model with different searching levels.
(a)(b)The two-zone MG model was then used to simulate the same case. The two-zone MG model enables the use of different searching levels in the burnt and unburned regions. The unburned region has higher chemical potential and therefore is more sensitive to the combustion model details. Applying a higher searching level can lead to failure in accurately predicting autoignition in the unburned region.The mixture in the burnt region is usually close to the state of chemical equilibrium. Thus, it is safe to apply global search, which results in the most aggressive grouping and significantly reduces computational cost. Simulations with global search in the burnt region (5) and different searching level (0 to 2, labeled as “MG0+5”, “MG1+5”, and “MG2+5”, resp.) in the unburned region were conducted on the same case. Table 3 lists the computational times of these simulations. The predicted pressure and HRR using the two-zone MG model were compared with the original GAMUT combustion model, as illustrated in Figure 5. It can be seen from Figure 5 that the two-zone MG model gives much better predictions than the one-zone MG model. Only when the searching level in the unburned region is greater than or equal to 2, do the curves of pressure and HRR deviate from the results predicted using the original GAMUT combustion model. Figure 6 shows total mole of in-cylinder NOx (NO2, NO, and N2O) computed using the original GAMUT combustion model and the multigrid model. The numerical results of the multigrid model are very close to the one of the original GAMUT combustion model before crank angle of 30 deg ATDC. Due to earlier predicted main combustion event when searching level is set to 2, NOx formation is slightly earlier than others. After crank angle of 30 deg ATDC, the multigrid model consistently predicts slightly higher NOx emission. Figure 7 shows iso-surface of the propagating flame front with contour plot of temperature on a cut plane right below the fire deck at CA=3 and 18 deg ATDC. Figure 8 shows iso-surface of the flame front on the same cut plane with contour plot of NOx at CA=18 deg ATDC. Autoignition at CA=18 deg ATDC is predicted by all the models. Before autoignition event (e.g., CA=3 deg ATDC), the multigrid models give very similar results comparing to the one from original GAMUT combustion model. It can be seen that using searching level of 2 overpredicted the hot spots in the unburnt region and consequently overpredicted the NOx in the corresponding regions. When searching level is set to 1, the predicted temperature and NOx distributions are very close to the results of the original GAMUT combustion model. In overall, applying a global search in the burnt region is very reliable, and the prediction accuracy mainly depends on the searching level in the unburned region.Table 3
Computational time of two-zone MG model with different combination of searching levels (CA=-6.3 to 70).
Casew/o MGMG0+5MG1+5MG2+5CPU time [hrs]173.335.727.413.7Wall time [hrs]48.613.111.17.2Speedup [-]3.613.215.624.1Pressure and HRR of SI engine simulated using two-zone MG model with different searching levels.
(a)(b)Figure 6
Total mole of NOx emission simulated using two-zone MG model with different searching levels.Iso-surface of flame front and temperature distribution atCA=3 (a) and 18 (b) deg ATDC.
(a)(b)Figure 8
Iso-surface of flame front  (black line) and NOx distribution atCA=18 deg ATDC.Although more conservative searching levels were applied to the unburned region, it can be seen from Table3 that more savings in computational time were obtained. Even when the MG model is not used in the unburned region (searching level=0), about 5 times speedup in computational time was achieved, which is better than all the one-zone MG models except for the global search in terms of its computational efficiency. If the searching level in the unburned region is one, the model is able to reproduce the results of the original GAMUT combustion model satisfactorily. In this case, more than 6 times speedup in computational time was achieved, corresponding to more than 15 times speedup in terms of wall time.As shown above, the unburned gas region is very sensitive to the details of the combustion solution technique. The risk of using the MG model is that it may artificially enhance species and temperature mixing. When the one-zone MG model is employed, all the computational cells near the propagating flame front may be grouped. The cells in the unburned region have lower temperature than those in the burnt region that will then have higher temperature. Consequently, this can lead to a faster flame propagation speed and even autoignition. The two-zone MG model avoids such artificially enhanced mixing in temperature between the burnt and unburned regions. Thus, it gives much better prediction in terms of pressure and HRR. However, as shown in Figure5, the use of too aggressive searching level in the unburned region will result in faster flame propagation speed during the main combustion stage, at least for the current stratified GDI engine case. Using a low searching level is one option. Another option is to skip the high-temperature computational cells during the mapping procedure, that is, to avoid allowing ignition-sensitive cells to be grouped.As a test of this concept, a temperature criterion was set to 800 K. In this case, cells are skipped if the temperature exceeds 800 K. Again, different searching levels for the unburned region were tested. The corresponding predicted pressure and HRR are plotted in Figure9 and the computational times are listed in Table 4. It can be seen from Figure 9 that by excluding the high-temperature cells, the MG model with almost all the searching levels gives close enough predictions to the original GAMUT combustion model. Only when global searching was applied in the unburned region, does the predicted main combustion occur slightly earlier than the one predicted by GAMUT model. As a compromise, fewer cells were grouped since the high-temperature cells are excluded. The speedup in computational times is similar to the ones in Table 3, even though more aggressive searching was used. Figure 10 shows the corresponding NOx emissions computed using the original GAMUT combustion model and the multigrid model. Similar to Figure 6, the NOx emission predicted by original GAMUT combustion model and multigrid models very close before crank angle of 30 deg ATDC. Slight differences between them show up after crank angle of 30 deg ATDC. Figure 11 shows iso-surface of the propagating flame front with contour plot of temperature on a cut plane right below the fire deck at CA=3 and 18 deg ATDC. Figure 12 shows iso-surface of the flame front on the same cut plane with contour plot of NOx at CA=18 deg ATDC. When global searching is used, the hot spots and NOx emission in the unburnt region are slightly overpredicted.Table 4
Computational time of two-zone MG model with different combination of searching levels and skipping high-temperature cells (T>800 K) in the unburned region (CA=-6.3 to 70).
CasewoMGMG2+5MG3+5MG4+5MG5+5CPU time [hrs]173.329.329.228.327.4Wall time [hrs]48.611.911.61110.5Speedup3.614.614.915.816.5Pressure and HRR of SI engine simulated using two-zone MG model with different searching levels and skipping high-temperature cells (T>800 K) in the unburned region.
(a)(b)Figure 10
Total mole of NOx emission of SI engine simulated using two-zone MG model with different searching levels and skipping high-temperature cells (T>800 K) in the unburned region.Iso-surface of flame front and temperature distribution atCA=3 (a) and 18 (b) deg ATDC simulated using two-zone MG model with different searching levels and skipping high-temperature cells (T>800 K) in the unburned region.
(a)(b)Figure 12
Iso-surface of flame front  (black line) and NOx distribution atCA=18 deg ATDC simulated using two-zone MG model with different searching levels and skipping high-temperature cells (T>800 K) in the unburned region.A higher criterion temperature, 900 K, was also tested. However, the performance results were worse than those of 800 K and therefore were not pursued.


## 4. Conclusions

In the present paper a multigrid model was implemented together with a level-set-G-equation model for spark ignition engine combustion modeling using detailed chemistry (GAMUT combustion model). The model was validated against the numerical results of the original GAMUT combustion model. The following conclusions can be drawn.(i)
A one-zone MG model was unable to accurately reproduce the main combustion event predicted by the original GAMUT combustion model.(ii)
A two-zone MG model, which treats the unburned and burnt regions separately, was able to reproduce the results of the original GAMUT combustion model and large reduction in computer time was seen (1 to 2 order(s) of magnitude).(iii)
The burnt region is insensitive to the combustion model because the mixture in the burnt region is close to chemical equilibrium. Thus, global searching can be applied to the burnt region to maximize the computational efficiency without loss in accuracy.(iv)
The unburned region is very sensitive to the details of the combustion model implementation. Either using a low searching level or excluding high-temperature computational cells from the grouping in the unburned region was found to provide accurate and efficient predictions.



---

*Source: 201780-2010-07-28.xml*