PDF

Homogeneous Charge Compression Ignition of Methane
Introduction
Homogeneous Charge Compression Ignition (HCCI) engines are being considered as an alternative to traditional spark- and compression-ignition engines. As the name implies, a homogeneous fuel/oxidant mixture is auto-ignited by compression with simultaneous combustion occurring throughout the cylinder volume. Combustion temperatures under lean burn operation are relatively low, resulting in low levels of NOx emission. Furthermore, the fuel’s homogeneous nature, as well as the combustion process itself, lead to low levels of particulate matter being produced.
Although HCCI combustion shows promise, the method has several recurring problems: an important one to be addressed is ignition timing. This example examines the HCCI of methane, investigating ignition trends as a function of initial temperature, initial pressure, and fuel additives.
This example solves the mass and energy balances describing the detailed combustion of methane in a variable-volume system. The large amount of kinetic and thermodynamic data required to set up the problem is easily available by importing the relevant files into the Reaction Engineering interface.
Model Definition
It is difficult to form the uniform mixtures required for HCCI with conventional diesel fuel. Natural gas fuels, on the other hand, readily produce homogeneous mixtures and have the potential to serve as HCCI fuels. This example considers the combustion of methane, as described by the GRI-3.0 mechanism, incorporating a detailed reaction mechanism of 53 species taking part in 325 reactions. The files describing the reaction kinetics and thermodynamics of the GRI-3.0 mechanism are available on the Internet (Ref. 1), and you can import the files into the Reaction Engineering interface.
Variable Volume Reactor
This model represents the combustion cylinder with a perfectly mixed batch system of variable volume, a predefined reactor type available with the Reaction Engineering interface. Figure 1 shows an engine cylinder and it includes the relevant parameters to calculate the instantaneous cylinder volume.
Figure 1: The volume of a combustion cylinder can be expressed as a function of time with the slider-crank relationship. The diagram shows the key geometric parameters. La is the length of the crank arm, Lc is the length of the connecting rod, D equals the cylinder diameter, and α is the crank angle.
The volume change as a function of time is described by the slider-crank equation:
where V is the cylinder volume (SI unit: m3), Vc is the clearance volume (SI unit: m3), CR equals the compression ratio, and R denotes the ratio of the connecting rod to the crank arm (Lc/La). Further, α is the crank angle (SI unit: rad), which is also a function of time
where N is the engine speed in rpm, and t is the time (SI unit: s).
The engine specifications are:
Equation 1 includes the clearance volume Vc which is calculated from
(1)
Vs is the volume swept by the piston during a cycle from the equation
Figure 2 shows the calculated cylinder volume as a function of the crank angle. The piston is initially at bottom dead center (BDC), corresponding to a crank angle of 180 degrees.
Figure 2: Cylinder volume as a function of crank angle. The crank angle is defined as being zero at top dead center (TDC).
Methane Combustion reaction
The kinetic and thermodynamic data for methane combustion is available in the form of CHEMKIN data input files. The files are imported into the Reaction Engineering interface within the Reversible Reaction Group and Species Thermodynamics (belongs to Species Group) features. This automatically sets up the mass and energy balances for a batch reactor of constant volume.
In this example methane is combusted under lean conditions, that is, supplying more than the stoichiometric amount of oxidizer. The stoichiometric requirement of the oxidizer (air) to combust methane is found from the overall reaction:
Assuming that the composition of air is 21% oxygen and 79% nitrogen, the stoichiometric air-fuel ratio is
(2)
The equivalence ratio relates the actual air-fuel ratio to the stoichiometric requirements
(3)
This model sets the equivalence ratio to Φ = 0.5.
From Equation 2 and Equation 3 it is possible to calculate the molar fraction of fuel in the reacting mixture as
and subsequently the initial concentration is
Results and Discussion
Figure 3 shows the cylinder pressure as a function of time when a methane-air mixture is compressed and ignites. The piston starts at bottom dead center (BDC) and reaches top dead center (TDC) after 0.02 s. At BDC the pressure is set to 1.5 × 105 Pa, Φ is 0.5, and the compression ratio is CR = 15. The initial temperature is varied from 400 K to 800 K.
Figure 3: Pressure traces illustrating the compression and ignition of fuel in an engine cylinder. The initial temperature varies between 400 K and 800 K.
Consistent with literature results, methane does not ignite at an initial temperature of 400 K (Ref. 2). Furthermore, the induction delay decreases with increasing initial temperature. The induction delay time can be evaluated from the pressure gradient. For instance, the induction delay is 0.0193 s when Tinit = 500 K.
Figure 4 illustrates the pressure traces as the initial pressure varies from × 105 Pa to × 105 Pa. The initial temperature is 500 K. An increase in pressure means an increase in the species concentrations in the fuel-air mixture, resulting in the expected advance in ignition times.
Figure 4: Increased initial gas pressure advances ignition times.
As mentioned, a significant challenge to the realization of HCCI engines is ignition control. In this regard, combustion at TDC is suggested as the optimum timing (Ref. 3). These results show that the inlet temperature of the fuel-air mixture is a potential tuning parameter for ignition. However, relatively high inlet temperatures are often required for proper timing. This adversely affects engine performance because the trapped mass as well as the volumetric efficiency decreases. An alternative that facilitates ignition is to mix small amounts of additives into the fuel-air mixture (Ref. 4). These additives chemically activate the reaction mixture even at relatively low temperatures. This approach alleviates the requirements of high intake temperatures. Figure 5 shows how small amounts of formaldehyde (CH2O) cause ignition at an initial temperature of 400 K, which is a temperature insufficient to induce combustion with a pure methane fuel.
Figure 5: Small amounts of formaldehyde stimulate ignition of the fuel-air mixture.
The increased reactivity observed in the presence of CH2O is explained by the opening of a new chemical pathway leading to the formation of hydroxyl radicals. Specifically, CH2O reacts with O2 to produce H2O2:
H2O2, in turn, decomposes to reactive OH radicals, which subsequently react violently with the fuel molecules to cause ignition:
The results in the following figures show the species molar fractions of CH2O, HO2, H2O2, and OH during the combustion of methane. Figure 6 shows molar fraction plots for the case when 0.13% CH2O is added to the fuel; Figure 7 is the equivalent species plots for the case when pure methane is combusted. In each case conditions are tuned to produce ignition near TDC, thus providing a reference point for comparing the species concentrations.
Figure 6: Selected species molar fractions as a function of crank angle. 0.13 molar percent CH2O is added to the reacting mixture, which is initially at 400 K and 1.5 bar.
Figure 7: Selected species molar fraction as a function of crank angle. Only methane is combusted. The initial temperature is 469 K and the initial pressure is1.5 bar.
The implications of the CH2O reaction path are apparent by comparing Figure 6 and Figure 7: CH2O stimulates the production of HO2 and H2O2, which in turn produce OH radicals in amounts critical to fuel ignition.
References
1. G.P. Smith, D.M. Golden, M. Frenklach, N.W. Moriarty, B. Eiteneer, M. Goldenberg, C.T. Bowman, R.K. Hanson, S. Song, W.C. Gardiner, Jr., V. V. Lissianski, and Z. Qin, GRI-Mech 3.0 home page, http://www.me.berkeley.edu/gri_mech/.
2. S.B. Fiveland and D.N. Assanis, “A four-stroke homogeneous charge compression ignition engine simulation for combustion and performance studies, “SAE Paper 2000-01-0332, 2000.
3. D.L. Flowers, S.M. Aceves, C.K. Westbrook, J.R. Smith, and R.W. Dibble, “Detailed Chemical Kinetic Simulation of Natural Gas HCCI Combustion: Gas Composition Effects and Investigation of Control Strategies,” J. Eng. Gas Turbine Power, vol. 123, no. 2, pp. 433–439, 2001.
4. M.H. Morsy, “Ignition control of methane fueled homogeneous charge compression ignition engines using additives,” Fuel, vol. 86, no. 4, pp. 533–540, 2007.
Application Library path: Chemical_Reaction_Engineering_Module/Ideal_Tank_Reactors/compression_ignition
Notes about the COMSOL Implementation
The kinetic and thermodynamic data required for this model are available on the Internet. Find the GRI-Mech 3.0 input files at (Ref. 1):
http://www.me.berkeley.edu/gri_mech/version30/text30.html.
Download the reaction mechanism and rate coefficient file (grimech30.dat), as well as thermodynamic data file (thermo30.dat) and store them on your computer so you can import these into the Reaction Engineering interface.
Modeling Instructions
From the File menu, choose New.
New
In the New window, click Model Wizard.
Model Wizard
1
In the Model Wizard window, click 0D.
2
In the Select Physics tree, select Chemical Species Transport>Reaction Engineering (re).
3
Click Add.
4
Click Study.
5
In the Select Study tree, select Preset Studies>Time Dependent.
6
Click Done.
Reaction Engineering (re)
In the Model Builder window’s toolbar, click the Show button and select Discretization in the menu.
Global Definitions
Import the model parameters from a text file.
Parameters
1
On the Home toolbar, click Parameters.
2
In the Settings window for Parameters, locate the Parameters section.
3
Click Load from File.
4
Import also some necessary variables, amongst these the cylinder volume function, from a text file.
Definitions
Variables 1
1
In the Model Builder window, under Component 1 (comp1) right-click Definitions and choose Variables.
2
In the Settings window for Variables, locate the Variables section.
3
Click Load from File.
4
Reaction Engineering (re)
1
In the Model Builder window, under Component 1 (comp1) click Reaction Engineering (re).
2
In the Settings window for Reaction Engineering, locate the Reactor section.
3
From the Reactor type list, choose Batch.
4
Locate the Energy Balance section. From the list, choose Include.
Use uniform scaling of the concentration variables to improve the computational performance.
5
Click to expand the Discretization section. Select the Uniform scaling of concentration variables check box.
6
Locate the Mass Balance section. In the Vr text field, type Vol.
Import the reaction kinetics data, available as a CHEMKIN file (grimech30.dat).
Reversible Reaction Group 1
1
Right-click Component 1 (comp1)>Reaction Engineering (re) and choose Reversible Reaction Group.
2
In the Settings window for Reversible Reaction Group, click to expand the CHEMKIN import for kinetics section.
3
Locate the CHEMKIN Import for Kinetics section. Select the Import CHEMKIN data check box.
4
Click Browse.
5
6
Click Import.
Species Thermodynamics 1
Import also the thermodynamic data, available as a CHEMKIN file (thermo30.dat).
1
In the Model Builder window, expand the Component 1 (comp1)>Reaction Engineering (re)>Species Group 1 node, then click Species Thermodynamics 1.
2
In the Settings window for Species Thermodynamics, click to expand the CHEMKIN import for thermodynamic data section.
3
Locate the CHEMKIN Import for Thermodynamic Data section. Click Browse.
4
5
Click Import.
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1)>Reaction Engineering (re) click Initial Values 1.
2
In the Settings window for Initial Values, locate the General Parameters section.
3
In the T0 text field, type T_init.
4
Locate the Volumetric Species Initial Value section. In the table, enter the following settings:
Study 1
Step 1: Time Dependent
Set up the time dependent study, modify the default tolerance settings to improve the accuracy of the solution.
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Times text field, type 0 0.026.
4
From the Tolerance list, choose User controlled.
5
In the Relative tolerance text field, type 1e-6.
Solution 1 (sol1)
1
On the Study toolbar, click Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node, then click Time-Dependent Solver 1.
3
In the Settings window for Time-Dependent Solver, click to expand the Absolute tolerance section.
4
Locate the Absolute Tolerance section. From the Tolerance method list, choose Manual.
5
In the Absolute tolerance text field, type 1.0E-7.
6
Clear the Update scaled absolute tolerance check box.
7
On the Study toolbar, click Compute.
Results
Global 1
1
In the Model Builder window, expand the Concentration (re) node, then click Global 1.
2
In the Settings window for Global, click Replace Expression in the upper-right corner of the y-axis data section. From the menu, choose Model>Component 1>Reaction Engineering>comp1.re.c_N2 - Concentration.
3
Click Add Expression in the upper-right corner of the y-axis data section. From the menu, choose Model>Component 1>Reaction Engineering>comp1.re.c_CH4 - Concentration.
4
Click Add Expression in the upper-right corner of the y-axis data section. From the menu, choose Model>Component 1>Reaction Engineering>comp1.re.c_O2 - Concentration.
5
Click Add Expression in the upper-right corner of the y-axis data section. From the menu, choose Model>Component 1>Reaction Engineering>comp1.re.c_CH2O - Concentration.
6
Click to expand the Coloring and style section. Locate the Coloring and Style section. In the Width text field, type 2.
7
Click to expand the Legends section. From the Legends list, choose Manual.
8
9
On the Concentration (re) toolbar, click Plot.
10
Click the Zoom Extents button on the Graphics toolbar.
The following steps create a plot of the pressure versus time.
1D Plot Group 3
1
On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, click to expand the Title section.
3
From the Title type list, choose None.
Global 1
1
Right-click 1D Plot Group 3 and choose Global.
2
In the Settings window for Global, click Replace Expression in the upper-right corner of the y-axis data section. From the menu, choose Model>Component 1>Reaction Engineering>comp1.re.p - Pressure.
3
On the 1D Plot Group 3 toolbar, click Plot.
The following steps reproduce Figure 6.
1D Plot Group 4
1
On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Mole fraction in the Label text field.
3
Locate the Title section. From the Title type list, choose None.
4
Locate the Plot Settings section. Select the x-axis label check box.
5
6
Select the y-axis label check box.
7
Global 1
1
Right-click Mole fraction and choose Global.
2
In the Settings window for Global, click Replace Expression in the upper-right corner of the y-axis data section. From the menu, choose Model>Component 1>Reaction Engineering>comp1.re.c_OH - Concentration.
3
Locate the y-Axis Data section. In the table, enter the following settings:
4
Locate the x-Axis Data section. From the Parameter list, choose Expression.
5
In the Expression text field, type comp1.crank_angle.
6
Click to expand the Coloring and style section. Locate the Coloring and Style section. In the Width text field, type 2.
7
Click to expand the Legends section. From the Legends list, choose Manual.
8
Mole fraction
1
In the Model Builder window, under Results click Mole fraction.
2
In the Settings window for 1D Plot Group, locate the Axis section.
3
Select the y-axis log scale check box.
4
On the Mole fraction toolbar, click Plot.
Global 2
1
In the Model Builder window, under Results>Mole fraction right-click Global 1 and choose Duplicate.
2
In the Settings window for Global, locate the y-Axis Data section.
3
4
Global 3
1
Right-click Results>Mole fraction>Global 2 and choose Duplicate.
2
In the Settings window for Global, locate the y-Axis Data section.
3
4
Global 4
1
Right-click Results>Mole fraction>Global 3 and choose Duplicate.
2
In the Settings window for Global, locate the y-Axis Data section.
3
4
5
On the Mole fraction toolbar, click Plot.
Mole fraction
1
In the Model Builder window, under Results click Mole fraction.
2
In the Settings window for 1D Plot Group, locate the Axis section.
3
Select the Manual axis limits check box.
4
In the x minimum text field, type -30.
5
In the x maximum text field, type 30.
6
In the y minimum text field, type 1e-8.
7
In the y maximum text field, type 1e-2.
8
Click to expand the Legend section. From the Position list, choose Lower right.
9
On the Mole fraction toolbar, click Plot.
The following steps reproduce Figure 7. First change the temperature and the initial CH2O concentration, then resolve.
Global Definitions
Parameters
1
In the Model Builder window, expand the Global Definitions node, then click Parameters.
2
In the Settings window for Parameters, locate the Parameters section.
3
Study 1
On the Home toolbar, click Compute.
Results
Mole fraction
1
In the Model Builder window, under Results click Mole fraction.
2
On the Mole fraction toolbar, click Plot.
To reproduce Figure 3, create a parametric sweep over the intial temperature parameter.
Study 1
Parametric Sweep
1
On the Study toolbar, click Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
Click Add.
4
5
In the Model Builder window, click Study 1.
6
In the Settings window for Study, locate the Study Settings section.
7
Clear the Generate default plots check box.
8
On the Study toolbar, click Compute.
Results
1D Plot Group 3
1
In the Model Builder window, under Results click 1D Plot Group 3.
2
In the Settings window for 1D Plot Group, locate the Data section.
3
From the Data set list, choose Study 1/Parametric Solutions 1 (sol2).
4
On the 1D Plot Group 3 toolbar, click Plot.
5
Click to expand the Legend section. Locate the Axis section. Select the Manual axis limits check box.
6
In the x minimum text field, type 0.012.
7
In the x maximum text field, type 0.026.
8
In the y minimum text field, type 0.
9
In the y maximum text field, type 1.2e7.
10
Locate the Legend section. From the Position list, choose Upper left.
Global 1
1
In the Model Builder window, under Results>1D Plot Group 3 click Global 1.
2
In the Settings window for Global, click to expand the Coloring and style section.
3
Locate the Coloring and Style section. In the Width text field, type 2.
4
Click to expand the Legends section. From the Legends list, choose Manual.
5
6
On the 1D Plot Group 3 toolbar, click Plot.
Global Definitions
Parameters
To reproduce Figure 4, sweep instead over the initial pressure parameter. Set the initial temperature to 500 K first.
1
In the Model Builder window, under Global Definitions click Parameters.
2
In the Settings window for Parameters, locate the Parameters section.
3
Study 1
Parametric Sweep
1
In the Model Builder window, under Study 1 click Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
4
On the Home toolbar, click Compute.
Results
1D Plot Group 3
1
In the Model Builder window, under Results click 1D Plot Group 3.
2
In the Settings window for 1D Plot Group, locate the Axis section.
3
In the y maximum text field, type 2.5e7.
Global 1
1
In the Model Builder window, under Results>1D Plot Group 3 click Global 1.
2
In the Settings window for Global, locate the Legends section.
3
4
On the 1D Plot Group 3 toolbar, click Plot.
Global Definitions
Parameters
To reproduce Figure 5, sweep instead over the initial CH2O mole fraction. Set the initial temperature to 400 K first.
1
In the Model Builder window, under Global Definitions click Parameters.
2
In the Settings window for Parameters, locate the Parameters section.
3
Study 1
Parametric Sweep
1
In the Model Builder window, under Study 1 click Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
4
On the Home toolbar, click Compute.
Results
1D Plot Group 3
1
In the Model Builder window, under Results click 1D Plot Group 3.
2
In the Settings window for 1D Plot Group, type Reactor pressure in the Label text field.
3
Locate the Axis section. In the y maximum text field, type 1.8e7.
Global 1
1
In the Model Builder window, under Results>Reactor pressure click Global 1.
2
In the Settings window for Global, locate the Legends section.
3
4
On the Reactor pressure toolbar, click Plot.