Numerical Simulation of Combustion in Porous Media

1. Introduction

The demand of fossil fuel resources and air pollution are the major problems that caused by using of fuels. In recent years, many researchers have developed new methods for efficient combustion of fuels. Porous media combustion, also known as filtration combustion in a packed bed, due to the interaction between two different phases, solid and gas or liquid. The theory of filtration combustion involves a new type of flame with exothermic chemical reactions during fluid flow in a porous medium. The term 'filtration combustion' was introduced by Russian scientists for combustion of gas flow through porous media. This term does not correspond to western scientific terminology, still it can be found in special literature as a synonym to combustion within porous media (PM). This process facilitates a combustion process with stability in a wide range of reactant fluid velocities, air-fuel ratios, and power density. PM combustion has some unique characteristics. It gives rise to high radiant output, low NO_x (Oxides of Nitrogen) and CO (Carbon Monoxide) emissions, high flame speed and higher power density. A porous material means a material with connected voids that flow can easily penetrate through its structure. This technology is different from conventional combustion, with free flame, thin reaction zone and high temperature gradients. Compared to conventional combustion devices, the combustion efficiency of PM burner is reasonably high and better heat transfer from burned gases to unburned mixture, takes place. On the other hand, in PM combustion three modes of heat transfer conduction, convection, and radiation are significant. In addition, there is a better homogenization of temperature across the PM and the significant amount of radiation helps to preheat the incoming air-fuel mixture at upstream. The technique of premixed combustion within PM has been studied and applied to steady combustion with great success. The porous media combustion has proved to be one of the applicable options to solve the problems to a remarkable extent in both technical and economic perspectives. This technique has been used for both gaseous and liquid fuels in steady or unsteady combustion. Flame stability in PM with lean and rich mixtures, significant reduction in pollutants and increasing combustion efficiency, was proven.

In recent years many researchers investigated the PM combustion technology both experimentally and theoretically. The most of researches are in field of steady combustion in PM and few of them are about transient flame propagation, both approaches are employed in PM combustion. Steady combustion is widely used in radiant burners and surface combustor-heaters due to its high radiant emissivity of the solid. The combustion zone is stabilized by its solid. The other, transient leads to an unsteady reaction zone freely propagate as a filtration combustion wave in the downstream direction. Combustion in PM differs considerably from the homogeneous flames flame front. Considerable features of PM for application of combustion technology are: large specific surface area, excellent heat transfer properties, heat capacity, transparency for fluid flow, thermal resistance, mechanical resistance, recuperation of energy and electrical properties.

2. Background of combustion in porous media

2.1. Stability of flame

The flame stabilization and propagation in a PM are governed by the modified Peclet number:

$$Pe=\mathrm{ }\frac{\left(S_{L\mathrm{ }\mathrm{ }}\mathrm{ }{d}_m\mathrm{ }{C}_p\mathrm{ }\rho \right)}{k}$$E1

where S_L is the laminar flame speed, d_m is the equivalent diameter of the average hollow space of the porous material, C_p is the specific heat of the gas mixture, 𝜌 is the density of the gas mixture and k is the thermal conductivity of the gas mixture. For flame propagation through a porous material, the critical Peclet number of 65 has been found. Thus, Pe < 65 for quenching, and Pe > 65 for flame propagation.

2.2. Premixed and non-premixed mixture

PM may work with a premixed flow or with a non-premixed fuel flow. Premixed porous burners consist of two zones: the premixed fuel–air mixture first enters a hot solid matrix, where it is heated until it enters to the second hot solid matrix. Depending on its application, a third section, a compact heat exchanger may be added to the burner. The schematic of a two-layer premixed combustor is shown in Fig. 1.

A premixed mixture of methane-air and hydrogen-air were used in different burners with two section PM. Measurements show considerable reduction in the concentrations of NO_x, CO in the fluid gases. Also the effects of hydrogen addition to methane, were investigated. For the porous burners hydrogen was found to lower the NO_x emissions slightly, while for the other burners an increase, or no obvious effect, was found. The enhancement of the radiation flux from PM burners operating with non-premixed flames by using a vane-rotary burner, in which the swirling fuel flow was confined by an air duct, is necessary. They also studied emission characteristics of ceramic foam burners operating with non-premixed flame.

Figure 1.

Schematic of a premixed combustor (M. A. Mujeebu et al., 2009)

2.3. Reciprocating flow

For the excess enthalpy combustion, a combustion system using reciprocating flow in PM was introduced. By the reciprocating flow, the combustion gas enthalpy is regenerated into increase in enthalpy of the combustible gas through the PM, which store heat. For this technique a new arrangement of the PM that stabilized flame for a wide operating condition, was used. The mixture first flow in, and the gas and solid temperatures reaches to a maximum at the exit side. Then the flow direction is reversed by means of valves. On the reverse flow half-cycle, the fresh mixture encounters much higher solid temperatures at the entering side. Therefore, the amount of heat recycled becomes larger than that with the single flow direction. Hence in the reciprocating flow system, the heat transfer from the combustion gases raises the solid temperatures from both directions.

2.4. Hydrogen production

Production of hydrogen from gases such as methane and hydrogen-sulphide is another potential application of PM combustion. These reactants convert into products such as hydrogen, syngas (H₂ and‏ CO), and sulfur. For both methane and hydrogen-sulphide combustion, upstream propagation corresponded to the range of equivalence ratios from stoichiometry to 1.7, and downstream wave propagation was observed for ultra-rich (1.7- 4) mixtures.

2.5. Materials for porous media combustion

Aluminum oxide (Al₂O₃), silicon carbide (SiC), and zirconium dioxide (ZrO₂) proposed as suitable materials for application. Al₂O₃ and ZrO₂ were recognized as high temperature resistant materials. SiC shows good thermal shock resistance, mechanical strength, and conductive heat transport. SiC also has high melting point (3260 K), against cyclic thermal stress and strength retention at the peak regenerator temperature (1673 K), and excellent oxidation resistance. Metallic materials were found less suitable for PM because of their inadequate thermal stability and high thermal inertia. Fe–Cr–Al-alloys and nickel-base alloys were found suitable for some applications but they were said to be comparatively less heat resistant. Structures of ceramic foams with different base materials were observed to possess high porosity, good conduction heat transport, low thermal inertia, low radiation heat transport properties and relatively high pressure drop. The effective thermal conductivity of anisotropic porous composite medium could vary largely with the component fractions.

3. Numerical modeling

The thermophysical properties of the air such as density, thermal conductivity and specific heat are assumed to be functions of the temperature. Usually, the pressure drop through the porous burner is not that high (with high porosity of PM) and its effect on the thermophysical properties can be neglected. In general, the properties of the solid phase may be assumed to be constant and assumed that there is thermal non-equilibrium between the gas and solid phases. Therefore, there are two energy equations to model the energy transport in the system. The porous material can be assumed as a scattering, emitting and absorbing medium. Gaseous radiation is assumed to be negligible compared to the solid radiation.

3.1. One-dimensional modeling

Many researchers studied the combustion and heat transfer in porous radiant burners. The PM was assumed to emit, absorb, and scatter radiant energy. Non-local thermal equilibrium between the solid and the gas is assumed for and combustion was described by a one-step, multi-steps or kinetic reactions. The effect of the optical depth, scattering albedo, solid thermal conductivity, upstream environment reflectivity, and interphase heat transfer coupling on the burner performance can be considered. Also, low solid thermal conductivity, low scattering albedo, and high inlet environment reflectivity produced a high radiant efficiency. The system consisted of a packed bed or foams in which a natural gas–air mixture combusts inside it. Radiative heat transfer in the packed bed or foams was modeled as a diffusion process, and the flow and temperature distribution in the packed bed or foams can be determined. The numerically results usually were compared with available experimental data for a similar system. Subsequently, the one-dimensional predictions of methane/air combustion in inert PM can include full mechanism (49 species and 227 elemental reactions), skeletal mechanism (26 species and 77 elemental reactions), 4-step reduced mechanism (9 species and 1-step global mechanism). The effects of these models on temperature, species, burning speeds and pollutant emissions were examined by researchers. Experimental and numerical investigations found that the flammability limits of the gaseous mixture in PM were more sensitive to their geometric properties than the physical properties.

3.1.1. One-dimensional governing equations

The following assumptions are made to simplify the problem:

Gas radiation is neglected, gas flow and heat transfer are one-dimensional.

PM is considered to be non-catalytic, homogeneous and optically thick.

The radiation of solid phase is treated using Rosseland approximation.

The PM consists of solid dispersed homogeneously, and the porosity variation near the tube wall is neglected. Under the above assumptions, a set of differential equations can be obtained in the following form.

Continuity equation:

$$\frac{\partial }{\partial t}\left({\rho }_g\epsilon \right)+\frac{\partial }{\partial x}\left({\rho }_{g}u\epsilon \right)=0.$$E2

Gas phase energy equation:

$${\rho }_g{c}_g\epsilon \frac{\partial T_g}{\partial t}+{\rho }_g{c}_g\epsilon u\frac{\partial T_g}{\partial x}+{\displaystyle \sum _i{\rho }_g}\epsilon Y_i{V}_i{c}_{gi}\frac{\partial T_g}{\partial x}+\epsilon {\displaystyle \sum _i{\dot{\omega }}_i{h}_i{W}_i}=\epsilon \frac{\partial }{\partial x}\left((k_g+{\rho }_g{c}_g{D}_{\left|\right|}{}^D)\frac{\partial T_g}{\partial x}\right)-h_v(T_g-T_s).$$E3

Solid phase energy equation:

$${\rho }_S{c}_S(1-\epsilon )\frac{\partial T_s}{\partial t}=k_s(1-\epsilon )\frac{{\partial }^2{T}_s}{\partial x^2}+h_v(T_g-T_s)-\frac{d{q}_r}{dx}=0.$$E4

Species transport equation:

$${\rho }_g\epsilon \frac{\partial Y_i}{\partial t}+{\rho }_g\epsilon u\frac{\partial Y_i}{\partial x}+\frac{\partial }{\partial x}\left({\rho }_g\epsilon Y_i{V}_i\right)-\epsilon {\dot{\omega }}_i{W}_i=0.$$E5

$$V_i=-\left(D+D_{im}^d\right)\frac{1}{X_i}\frac{\partial X_i}{\partial x}.$$E6

density:

$$\rho =\text{ constant density}.$$E7

Radiation:

$$q_r(x)=-\frac{16}{3}\frac{\sigma T_s{}^3}{\beta }\frac{d{T}_s}{dx}.$$E8

3.1.2. Boundary conditions for one-dimensional model

The following boundary conditions are considered in the computations:

Inlet	$$T_g=T_{f,inlet}$$	$$\frac{\partial T_s}{\partial x}=0$$	$$\frac{\partial T_s}{\partial x}=0$$
Outlet	$$\frac{\partial T_g}{\partial x}=0$$	$$\frac{\partial T_s}{\partial x}=0$$	$$\frac{\partial Y_i}{\partial x}=0$$

3.1.3. Solution procedure

The time integration is performed using combination of implicit finite difference and finite volume method on a uniform- adaptive or non uniform mesh. The convective terms are upwinded and the diffusive terms are discretized using a second-order technique. The initial condition was also given to the system of discretized equations. For more information interested reader can see (Henneke et al, 1999)

3.2. Two-dimensional modeling

In the two dimensional model a comparison is made between the local thermal equilibrium and thermal non-equilibrium approaches by different researchers. The volume-averaged treatments unable to predict the pore-level, local high temperature region in the gas phase and the pore-level variation in the flame speed with respect to the flame location in the pore. The gas and the solid phases consider in non-local thermal equilibrium, and separate energy equations use for the two phases. The solid phase can be assumed absorbing, emitting and scattering, while the gas phase was considered transparent to radiation. The alternating direction implicit (ADI) scheme can be used to solve the transient two dimensional energy equations. Methane–air combustion with one-step, multi step and detailed chemical kinetics are used to model the combustion, like one dimensional combustion modeling. The radiative part of the energy equation can be modeled using simple Rosseland method until complex discrete ordinate method. In the flat plate burner, air flows axially through a constant area duct filled with a porous layer of thickness L. In the cylindrical and spherical burners, the air flows radially through an annular porous matrix.

3.2.1. Two dimensional governing equations

The energy equation for the gas phase is as follows:

$$\frac{\partial }{\partial t}\left(\phi {\rho }_g{C}_{pg}{T}_g\right)+\frac{1}{r^n}\frac{\partial }{\partial r}\left(\phi {\rho }_g{C}_{pg}{r}^{n}v{T}_g\right)=\frac{1}{r^n}\frac{\partial }{\partial r}\left(\phi k_g{r}^n\frac{\partial T_g}{\partial r}\right)-\left(1-\phi \right)h_v\left(T_g-T_s\right)+\phi {∆H}_c{S}_{fg}$$E9

Where 𝜑, 𝜌_g, r, C_pg, T_g, v, k_g, h_v, H_c and S_fg are the porosity, density, radial position, specific heat, temperature, velocity, thermal conductivity, volumetric heat transfer coefficient, enthalpy of combustion and rate of fuel consumption per unit volume, respectively. Subscripts g and s refer to gaseous and solid phases, respectively.

The energy equation for the solid phase is:

$$\frac{\partial }{\partial t}\left({\rho }_s{C}_s{T}_s\right)+\frac{1}{r^n}\frac{\partial }{\partial r}\left(r^n{k}_s\frac{\partial T_s}{\partial r}\right)-h_v\left(T_s-T_g\right)-\nabla .F.$$E10

The term F represents the radiative transport equation and is given by:

$$\nabla .F=-\left(1-\omega \right)(G-4E_b),$$E11

where ω is the single scattering albedo and the irradiance G is governed by

$${\nabla }^{2}G=\mathrm{ }{\eta }^2\left(G-4E_b\right),\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }{\eta }^2=3{\beta }^2\left(1-\omega \right)\left(1-\mathrm{g}\mathrm{\omega }\right),$$E12

where E_b is the Planck black body emitted flux and G is the radiative flux.

The conservation equation for the mass fraction of the fuel is given as follows:

$$\frac{\partial }{\partial t}\left({\rho }_g{m}_f\right)+\frac{1}{r^n}\frac{\partial }{\partial r}\left({\rho }_g{r}^{n}v{m}_f\right)=\frac{1}{r^n}\frac{\partial }{\partial r}\left(r^n{D}_{AB}{\rho }_g\frac{\partial m_f}{\partial r}\right)-S_{fg}\mathrm{ },$$E13

where m_f is the fuel mass fraction and D_AB is the diffusion coefficient. The n value is set to 0, 1 and 2 for spherical, radial and axial flow burners, respectively.

A single-step Arrhenius type chemical kinetic equation as given below is normally adopted in modeling the combustion:

$$S_{fg}=f{\rho }_{g\mathrm{ }}^2{m}_f{m}_{o_2}\exp \left(-\frac{E}{R{T}_g}\right),$$E14

where f, m_O2, E and R refer to pre-exponential factor, oxygen mass fraction, activation energy and gas constant, respectively.

3.2.2. Boundary conditions for two-dimensional model

The following boundary conditions are adopted for the gas, solid and species:

Gas:$T_{g|r={\mathrm{ }r}_{in}}=\mathrm{ }\mathrm{ }0\mathrm{ }\mathrm{ }\mathrm{a}\mathrm{t}\mathrm{ }\mathrm{ }\mathrm{ }r=r_{in}$

$${\left.\frac{\partial T_g}{\partial r}\right|}_{r\mathrm{ }={\mathrm{ }r}_{out}}=\mathrm{ }\mathrm{ }{T}_{in}\mathrm{ }\mathrm{ }\mathrm{a}\mathrm{t}\mathrm{ }\mathrm{ }\mathrm{ }r=r_{out}.$$E15

Solid:

$h_{in}\left[(T_{g\mathrm{ },\mathrm{ }in}-T_{s|r={\mathrm{ }r}_{in}}\mathrm{ })\right]+\sigma {ϵ}_{in}\left[({T_{\mathrm{ }}}_{in\mathrm{ },\mathrm{ }amb}^4-{\left.T_s^4\right|}_{s|r={\mathrm{ }r}_{in}}\mathrm{ })\right]=\mathrm{ }-k_s{\left.\frac{\partial T_s}{\partial r}\right|}_{r\mathrm{ }={\mathrm{ }r}_{in}}\mathrm{ }\mathrm{ }\mathrm{a}\mathrm{t}\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }r=r_{in}.$

$$h_{out}\left[(T_{out\mathrm{ },\mathrm{ }in}-T_{s|r={\mathrm{ }r}_{out}}\mathrm{ })\right]+\sigma {ϵ}_{out}\left[({T_{\mathrm{ }}}_{out\mathrm{ },\mathrm{ }amb}^4-{\left.T_s^4\right|}_{s|r=out}\mathrm{ })\right]=\mathrm{ }-k_s{\left.\frac{\partial T_s}{\partial r}\right|}_{r\mathrm{ }={\mathrm{ }r}_{out}}\mathrm{ }\mathrm{ }\mathrm{a}\mathrm{t}\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }r=r_{out}.$$E16

Species:

$m_f=\mathrm{ }{m}_{f\mathrm{ },\mathrm{ }in}\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{a}\mathrm{t}\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }r=r_{in}.$

$$\frac{{\partial m}_f}{\partial r}=\mathrm{ }0\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{a}\mathrm{t}\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }r=r_{out}.$$E17

The control volume approach, or the finite-difference method, can be used to solve the governing equations. The solution is advanced in time by using a fully implicit technique and this was necessary due to the stiffness of the governing matrix of the problem. Also, it is necessary to use an adaptive grid, or a very fine grid, to insure the accuracy of the solution.

3.2.3. Solution procedure

The based code can be solved using alternative direction implicit (ADI) finite volume formulation. The pressure field is solved using the SIMPLE method in steady state and PISO in transient state. The gridding system should prove to be sufficient by testing several grids sizes. For more information interested reader can see (Mohammadi, 2010, Hackert et al, 1998).

4. Three-dimensional modeling

Although there exist many simulation of one-dimensional and two dimensional model in literature but only several papers discussed three-dimensional simulation in PM. Three-dimensional simulations of combustion inside the PM burner, are complex and time consuming. A finite-volume calculation for three-dimensional reacting flow in a porous burner can be used. The Navier–Stokes equations, energy for both phase of PM and species transport equations were solved, and radiative heat transfer under local thermal non equilibrium between the solid and gas phases was considered.

4.1. Three-dimensional governing equations

Following assumptions were used in modeling and simulation of the PM in modified KIVA-3V code:

There is thermal non-equilibrium between gas and solid phases.

Solid is homogeneous, isotropic, variable property with temperature and has no catalyst effects.

Only radiation heat transfer from the solid phase is considered and the gas phase is transparent.

Based on the above assumptions, the general governing equations are simplified as follows:Continuity equation for species i:

$$\frac{\partial \left(\phi {\rho }_{i\mathrm{ }}\right)}{\partial t}\mathrm{ }+\mathrm{ }\nabla .\left(\phi {\rho }_i\mathrm{ }u\mathrm{ }\right)=\mathrm{ }\mathrm{ }\nabla .⌈\rho \mathrm{ }{\phi \mathrm{ }D}_{im\mathrm{ }}\nabla \left(\frac{{\rho }_i}{\rho }\right)⌉+\phi \mathrm{ }{\dot{\rho }}_i^c+\mathrm{ }{\dot{\rho }}^s\mathrm{ }\mathrm{ }{\delta }_i\mathrm{ }\mathrm{ }$$E18

where diffusion coefficient modified according to kinetic theory of gases. ρ is mixture density, φ is porosity (void fraction) of PM, ρ_i is density of species i, D_im is diffusion coefficient of species i in the mixture and u is the velocity vector.

Gas phase momentum equation:

$$\frac{\partial \left(\rho u\right)}{\partial t}+\nabla .\left(\rho u u\right)=-\frac{1}{a^2}\nabla P-\nabla \left(\frac{2}{3} \rho k\right)+\nabla .\sigma +F^s-\left(\frac{∆P}{∆L}\right)\mathrm{ }.$$E19

The last term on the right-hand side of Eq. (19) is due to pressure drop caused by PM according to Ergan equation and σ is stress tensor. The dimensionless quantity a is used in conjunction with the Pressure Gradient Scaling (PGS) method. This is a method for enhancing computational efficiency in low Mach number flows, where the pressure is nearly uniform.

$$\left(\frac{∆P}{∆L}\right)= \left(\frac{\mu }{\alpha } u+ c_2 \frac{1}{2} \rho \left|u\right| u\right).$$E20

where 𝛼 is the permeability and c₂ is the inertial resistance factor of PM, which are determined according to Eq. 21.

$$\alpha =\frac{d^2}{150}\frac{{ϵ}^3}{{(1-ϵ)}^2}\mathrm{ }\mathrm{ }\mathrm{ },\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }{c}_2=\frac{3.5}{d}\frac{{(1-ϵ)}^{}}{{ϵ}^3}\mathrm{ }.\mathrm{ }$$E21

Gas phase energy equation:

$\frac{\partial }{\partial t}\left(\phi \rho c_p{T}_g\right)+\nabla .\left(\phi \rho {c_{p}T}_{g}u\right)+\phi \sum _i{\dot{\omega }}_i{H}_i{W}_i=-\phi P\mathrm{ }\nabla .u+\phi A_0\rho ϵ+\left(1-A_0\right)\sigma :\nabla u$

$$\phi \mathrm{ }\nabla .\left((k_g+{\rho }_g{c}_g{D}_{\parallel }^d)\nabla T_g\right)-h_v\left(T_g-T_s\right)+{\dot{Q}}^s\mathrm{ }.$$E22

where the fourth term on the right hand side is the conduction heat transfer due to thermal conductivity of fluid and dispersion term due to existence of PM and the fifth term represent the convective heat transfer between gas and solid phase of PM which are determined according to Eqs. (23 – 25):

$$D_{\parallel }^d=0.5\mathrm{ }{{\alpha }_{g}Pe}_{}$$E23

$${\text{Nu}}_{\text{v}}=\text{ 2 }+\text{ 1}.{\text{1 Re}}^{0.\text{6}}{\text{Pr}}^{0.\text{33}},$$E24

$$h_v=\frac{6\phi }{d^2}\mathrm{ }{k}_g{Nu}_v$$E25

c_p is specific heat of the mixture, T_g is gas temperature, Y_i is mass fraction species i, is rate of reaction i, H_i is enthalpy of species i, W_i is molecular weight of species i, k_g is thermal conductivity of the fluid, is thermal dispersion coefficient along the length of the porous medium, h_v is volumetric heat transfer coefficient between solid and gas phase of PM. Correlation (23 - 25) was estimated from experimental data by Wakao and Kaguei for heat transfer between packed beds and fluid.

Solid phase energy equation:

$$\frac{\partial }{\partial t}\left(\left(1-\phi \right)\mathrm{ }{\rho }_s\mathrm{ }{c}_s\mathrm{ }{T}_s\right)=\mathrm{ }\nabla \bullet \left[k_s\mathrm{ }\left(1-\phi \right)\nabla T_s\right]+h_{v\mathrm{ }}\left(T_g-\mathrm{ }{T}_s\right)-\mathrm{ }\nabla .q_r.$$E26

T_s is solid temperature, ρ_s is solid density, c_s is specific heat of solid phase, k_s is thermal conductivity of solid phase, q_r is the radiation heat transfer in solid in r direction.

Chemical species continuity equation:

$$\frac{\partial }{\partial t}\left(\phi {\rho Y}_i\mathrm{ }\right)+\mathrm{ }\nabla .\left(\phi \rho Y_i\mathrm{ }u\mathrm{ }\right)+\nabla .\left(\phi \rho Y_i\mathrm{ }{v}_i\mathrm{ }\right)-\mathrm{ }\phi \mathrm{ }{\dot{\omega }}_i\mathrm{ }{W}_i=0,\mathrm{ }$$E27

$$v_i=\mathrm{ }-\left(D+D_{m\parallel }^d\right)\frac{1}{X_i}\mathrm{ }\nabla X_i,$$E28

$$Pe=\frac{{\rho }_{\mathrm{ }{C}_{p\mathrm{ }\left|u\right|\mathrm{ }d}}}{k_g}\mathrm{ }$$E29

$$D_{m\parallel }^d=0.5\mathrm{ }{D}_{im}Pe$$E30

X_i is molar fraction of species i, Pe is Peclet number, d is diameter of sphere in packed bed, is species dispersion coefficient, Y_i is mass fraction of species i and v_i is diffusion velocity of species i in the mixture.

Turbulence model

Since there is no model presented for simulation of turbulent compressible-flow in PM by any researcher. Hence the basic κ- ε equations were used without any modification.

The transport equation for κ turbulent kinetic energy:

$$\frac{\partial \left(\rho k\right)}{\partial t}+\mathrm{ }\nabla .\left(\rho u\mathrm{ }k\right)=\mathrm{ }-\frac{2}{3}\mathrm{ }\mathrm{\rho }\mathrm{ }\mathrm{\kappa }\mathrm{ }\nabla .u+\mathrm{ }\mathrm{ }\sigma \mathrm{ }:\nabla u+\nabla .\left[\left(\frac{\mu }{{Pr}_k}\right)\nabla k\right]-\rho \epsilon +{\dot{W}\mathrm{ }}^s\mathrm{ },$$E31

where σ is stress tensor. With a similar consideration dissipation rate, ε:

$$\frac{\partial \left(\rho \epsilon \right)}{\partial t}+\nabla .\left(\rho u\epsilon \right)-\left(\frac{2}{3}C\in {}_1{}^{}-C\in {}_3{}^{}\right)\rho \epsilon \nabla .u+\nabla .\left[\left(\frac{\mu }{{\mathrm{Pr}}_{\epsilon }}\right)\nabla \in \right]+\frac{\in }{k}\left[C{\in }_1\sigma :\nabla u-C{\in }_2\rho \epsilon +c_s{\dot{W}}^s\right]$$E32

The parameters,,, and are constant whose values are determined from experiments and some theoretical considerations and is viscous stress tensor.

Equation of state:

$$P=\mathrm{ }\rho R\mathrm{ }{T}_g\mathrm{ }/\mathrm{ }\stackrel{-}{W}$$E33

R is universal gas constant, average molecular weight of mixture, P is pressure inside the combustion chamber and the PM.

Radiation model

Due to extreme temperature of combustion zone and solid phase, radiation heat transfer is very important. Gas phase radiation in comparison with solid phase radiation that has a high absorption coefficient, is negligible. Several relations for modeling of radiation heat transfer and derived radiation intensity are presented. The heat source term, due to radiation in solid phase that appears in Eq. 26, can be calculated from Rosseland model.

$$q_r=\mathrm{ }-\mathrm{ }\frac{16}{3}\mathrm{ }\frac{\sigma \mathrm{ }{T_s}^3}{\beta }\mathrm{ }\nabla T_s$$E34

where σ is Boltzmann constant and β is extinction coefficient.

Combustion model

Chemical mechanism for oxidation of methane fuel is considered chemical production rate:

$${\omega }_i={\displaystyle \sum _{i=1}^{NR}(v_{k,i}^{\text{'}\text{'}}}-v_{k,i}^{\text{'}})R_{i,}$$E35

and are stoichiometric coefficients. Combustion process includes ten equations and twelve species. These equations are presented in Table.1, which includes one-step reaction for methane fuel and equations 2-4, are Zeldovich mechanism for NO formation. Rate of reactions are computed with Arrhenius method. But for six other equations that reaction rate is very quick relative to last four equations hence, equilibrium reactions are considered. In order to consider effects of turbulence on combustion the common Eddy-Dissipation model can be used.

Number	Equation
1	CH4 + 2 O2 → CO2 + 2 H2O
2	O2 + 2 N2 → 2 N + 2 NO
3	2 O2 + N2 → 2 O + 2 NO
4	N2 + 2 OH → 2 H + 2 NO
5	${\mathrm{H}}_2\leftrightarrows 2\mathrm{ }\mathrm{H}$
6	${\mathrm{O}}_2\leftrightarrows 2\mathrm{ }\mathrm{O}$
7	${\mathrm{N}}_2\leftrightarrows 2\mathrm{N}$
8	$\mathrm{O}+\mathrm{H}\leftrightarrows \mathrm{O}\mathrm{ }\mathrm{H}$
9	${\mathrm{O}}_2+\mathrm{ }2\mathrm{ }{\mathrm{H}}_2\mathrm{O}\leftrightarrows 4\mathrm{ }\mathrm{O}\mathrm{H}$
10	${\mathrm{O}}_2+\mathrm{ }2\mathrm{ }\mathrm{C}\mathrm{O}\mathrm{ }\leftrightarrows 2\mathrm{ }{\mathrm{C}\mathrm{O}}_2$

Table 2.

Kinetic and equilibrium reactions

Spalding suggested that combustion processes are best described by focusing attention on coherent bodies of gas, which squeezed and stretched during their travel through the flame. This model relates the local and instantaneous turbulent combustion rate to the fuel mass fraction and the characteristic time scale of turbulence. The application of this model requires adjustment of a specific coefficient and ignition time to match the experimental combustion rate with the computational combustion rate. In combustion chamber of engine, high turbulence intensity exists and hence combustion for such device lies in the flamelets-in-eddies regime. The intrinsic idea behind the model is that the rate of combustion is determined by the rate at which parcel of unburned gas are broken down into the smaller ones, such that there is sufficient interfacial area between the unburned mixture and hot gases to permit reaction and also the turbulence length scale which is quite important can determine the turbulent burning rates. In this case, the coefficients of model were determined from experimental analysis of conventional engine.

5. Foundation of reaction

5.1. Chemical equilibrium

A general chemical equilibrium reaction with v′_i,s and v′'_i,s representing the stoichiometric coefficient of reaction and products for the chemical species M_i

$$\sum _{i=1}^N{v^{\mathrm{\text{'}}}}_{i,s}\mathrm{ }{M}_i\mathrm{ }\rightleftarrows \mathrm{ }\sum _{i=1}^N{v^{\mathrm{\text{'}}\mathrm{\text{'}}}}_{i,s}\mathrm{ }{M}_i\mathrm{ }.$$E36

State of equilibrium can be interpreted as a situation, in which both the forward as well as reverse reactions progresses with identical speed.

$$K_p\left(T\right)=\mathrm{ }\prod _i{\left(\frac{pi}{p0}\right)}^{\left(v^{\mathrm{\text{'}}\mathrm{\text{'}}}i,s\mathrm{ }–\mathrm{ }{v}^{\mathrm{\text{'}}}i,s\right)}.$$E37

The equilibrium constant K_p now contains the information about the equilibrium material composition in term of partial pressure p_i of the various species i. For more information interested reader can see (Mohammadi 2010).

5.2. Reaction kinetics

A one step chemical reaction of arbitrary complexity can be represented by the following stoichiometric equation:

$$\sum _{i=1}^N{v^{\mathrm{\text{'}}}}_i\mathrm{ }{M}_i\mathrm{ }\to \mathrm{ }\sum _{i=1}^N{v^{\mathrm{\text{'}}\mathrm{\text{'}}}}_i\mathrm{ }{M}_i,\mathrm{ }\mathrm{ }$$E38

where v′_s are the stoichiometric coefficient of reactions and v′'_i representing the stoichiometric coefficient of products, M_i the specification of molecule of ith chemical species, and N total number of component involved. Usually, they are represented with an Arrhenious formulation form:

$$K=A{T}^b\exp \left(-\frac{E_A}{RT}\right).$$E39

The constant A and the exponent b as well as the so-called activation energy E_A are summarized for many chemical reactions in extensive table.

6. The finite volume method

Customarily, CFD codes work with the finite volume method. This approach guaranties the numerical preservation of conservative quantities for the incompressible flows. The finite volume (FV) method uses the integral form of the conservation of equations. The solution domain is subdivided into a finite number of control volumes, and the conservation equations are applied to each control volume. As result, an algebraic equation for each CV is obtained. The FV method accommodates any type of grid, so it is suitable for complex geometries. However, the computational mesh ideally, be built hexahedratically. The conservation law for transport of a scalar in an unsteady flow has the general form:

$$\frac{\partial }{\partial t}\left(\rho \mathrm{ }\Phi \right)+\nabla .\left(\rho u\Phi \right)=\mathrm{ }\nabla .\left(\Gamma \nabla \Phi \right)+\mathrm{ }{S}_{\Phi }$$E40

(𝜌uΦ) designates convection, (𝛤𝛻Φ) diffusion flows of and S_Φ the corresponding local source. For more information interested reader can see (Mohammadi 2010).

6.1. The finite volume equations formulation

Finite volume equations are derived by the integration of above differential equations over finite control volumes that taken together fully cover the entire domain of interest (Fig. 2). These control volumes are called ''cells'' P, for which the fluid-property value, are regarded as representative of the whole cell. It is surrounded by neighboring nodes which we shall denote by N, S, E, W, B and T. Cells and nodes for velocity components are ''staggered'' relative to those for all other variables.

Figure 2.

Computational molecule in 3D domain (Patankar, 2002)

For more information interested reader can see (Mohammadi, 2010).

6.2. Discretisation and numerical solution of the momentum equation

Finally, the momentum equation for the calculation of velocity and pressure by use of continuity equation should be considered. For numerical reasons, it is recommendable to resort to so called staggered grid, i.e. pressure and velocity are calculated on computational grids shifted to each other, the pressure for example in the cells and the velocity on the nodes. The calculations of velocity commonly take place iteratively, for which several algorithms are known (e.g. SIMPLE, PISO, SIMPLER…). In final analysis, all have the fact in common that is first step the momentum equation is solved for the velocities of momentums kept constant. In the second step, pressure corrections are then calculated with the help of a Poisson equation. For pressure with these pressure corrections, new velocities are then calculated again, and that again, until a pre-given break off threshold for the convergence is reached.

6.2.1. Discretisation of transient convection diffusion equation

Transient three dimensional convection diffusion of a general property Φ in a velocity field that govern by equation (40). The fully implicit discretisation equation is:

$${\text{a}}_{\text{p}}{\Phi }_{\text{p}}={\text{ a}}_{\text{w}}{\Phi }_{\text{w}}+{\text{ a}}_{\text{E}}{\Phi }_{\text{E}}+{\text{ a}}_{\text{s}}{\Phi }_{\text{s}}+{\text{ a}}_{\text{N}}{\Phi }_{\text{N}}+{\text{ a}}_{\text{B}}{\Phi }_{\text{B}}+{\text{ a}}_{\text{T}}{\Phi }_{\text{T}}+{\text{ a}}^{º}{}_{\text{p}}{\Phi }^{º}{}_{\text{p}}+{\text{ S}}_{\text{u}}{,}_{}$$E41

where:

$${\text{a}}_{\text{p}}={\text{ a}}_{\text{w}}+{\text{ a}}_{\text{E}}+{\text{ a}}_{\text{s}}+{\text{ a}}_{\text{N}}+{\text{ a}}_{\text{B}}+{\text{ a}}_{\text{T}}+{\text{ a}}^{º}{}_{\text{p}}+\Delta \text{F}-{\text{ S}}_{\text{p} }$$E42

with $a_p^{°}=\mathrm{ }\frac{{\rho }_{p\mathrm{ }\Delta V}^{\mathrm{º}}}{\Delta }$ and $\overline{S}\Delta V=S_u+S_p{\Phi }_p$

For more information interested reader can see (Mohammadi, 2010).

7. Simple algorithm

As discussed in the preceding section, the governing equation for the flow may be solved in terms of derived variables, or in term of primitive variables consisting of the velocity components and the pressure.

Figure 3.

Staggered location for the velocity components in a two dimensional flow (Patankar, 1980)

However, in the advent of Simple (Semi Implicit Method for Pressure Linked Equations) algorithm, along with its revised version Simpler and the enhancement such as Simplec, the solution of the equations using primitive variable approach has become very attractive. In fact, Simple and Simple like algorithms are extremely popular for the solution of problems involving convective flow and transport. The basic approach involves the control volume formulation, with the staggered grid, as outlined in the proceeding section. This avoids the appearance of physically unrealistic wavy velocity fields in the solution to equations. The pressure at a chosen point is taken at arbitrary value and the pressures at other points are calculated as differences from the chosen pressure value.

Following (Patankar, 1980), if a guessed pressure field p^* is taken, the corresponding velocity field can be calculated from the discretised equations for the control volume shown in Fig. 4 These equations are of the form:

$$a_e{u}_{e\mathrm{ }\mathrm{ }}=\mathrm{ }\sum a_{nb}\mathrm{ }{u}_{nb\mathrm{ }\mathrm{ }}^{\mathrm{*}}+\mathrm{ }b+\left(p_p^{\mathrm{*}}-p_E^{\mathrm{*}}\right)\mathrm{ }{A}_e,$$E43

where the asterisk on the velocity indicates the erroneous velocity field based on guessed pressure field. Here, a_nb is a coefficient that accounts for the combined convection-diffusion at the faces of the control volume, with nb referring to the neighbors e to the control volume, b includes the source terms except the pressure gradient, and A_e is the area on which pressure acts, being Δy*Δz for 3D. The numbers of neighbor terms are 6 for three dimensional ones. Similar equations can be written for v_n^* and w_t^*, where t lies on the z-direction grid line between grid points P and T. if p is the correct pressure and p is the correct pressure and p' the pressure correction, we may write:

$$\text{p }={\text{ p}}^{*}+\text{ p}\prime ,\text{ u }={\text{ u}}^{*}+\text{ u}\prime ,\text{ v }={\text{ v}}^{*}+\text{ v}\prime ,\text{ w }={\text{ w}}^{*}+\text{ w}\prime ,$$E44

where the prime indicate corrections needed to reach the correct values that satisfy the continuity equation. Omitting the correction terms due to the neighbors, an iterative solution may be developed to solve for the pressure and the velocity field. Then, the velocity correction formula becomes:

$u_{e\mathrm{ }\mathrm{ }}=u_e^{\mathrm{*}}+\mathrm{ }\frac{A_e}{a_e}\left(p_p^{\mathrm{\text{'}}}-p_E^{\mathrm{\text{'}}}\right)$

$$v_{n\mathrm{ }\mathrm{ }}=v_e^{\mathrm{*}}+\mathrm{ }\frac{A_n}{a_e}\left(p_p^{\mathrm{\text{'}}}-p_N^{\mathrm{\text{'}}}\right),\mathrm{ }$$E45

And similarly for w_t. From the time dependent continuity equation, the pressure correction equation in then developed as:

$${\text{a}}_{\text{p}}\text{p}{\text{'}}_{\text{p} }={\text{ a}}_{\text{E}}\text{p}{\text{'}}_{\text{e}}+{\text{ a}}_{\text{w} }\text{p}{\text{'}}_{\text{w}}+{\text{ a}}_{\text{N}}\text{p}{\text{'}}_{\text{N}}+{\text{ a}}_{\text{s}}\text{p}{\text{'}}_{\text{s}}+{\text{ a}}_{\text{T}}\text{p}{\text{'}}_{\text{T}}+ {\text{a}}_{\text{B}}\text{p}{\text{'}}_{\text{B}}+\text{b},$$E46

where b is a mass source which must be eliminated through pressure correction so that continuity is satisfied. Here, T and B are neighboring grid points on the z direction grid line.

Figure 4.

Control volume for driving the pressure correction equation (Patankar, 1980)

The simple algorithm has the following main steps:

Guess the pressure field p*.

Solve the momentum equation to obtain u*,v*, and w*.

Solve the pressure correction equation to obtain p'.

Add p' to p* to obtain the corrected pressure p.

Calculate u, v and w from u*, v* and w* using velocity correction equations.

Treat the corrected pressure p as the new guess p* and iterate the preceding procedure to convergence.

The revised version Simpler is quite similar to preceding algorithm and was developed mainly to improve the rate of convergence. In this case, the mail steps are:

Guess the velocity field

Solve the pressure equation, which is similar to pressure correction equation, Eq. 46, to obtain the pressure distribution. In this equation p' is replaced by p and a different expression arise for b.

Treating the pressure field as p*, solve the momentum equations to obtain u*, v* and w*.

Solve the pressure correction equation to obtain p'.

Correct the velocity field but not the pressure.

Use the velocity field as the guessed distribution and iterate the preceding procedure to convergence.

The pressure at any arbitrary point in the computational domain is specified and pressure differentials from this value are computed. The boundary condition may be a given pressure, which makes p' = 0, or a given normal velocity which makes the velocity a known quantity at the boundary and not a quality to be corrected so that p' at the boundary is not needed. For further details, (Patankar, 1980) may be consulted.

8. Multi-step reaction models

Models for premixed PM combustion are complicated by the highly nonlinear radiative exchange terms in the energy equation for the solid matrix in addition to the stiffness of the set of gas phase equations. Therefore researchers have simulated the gas phase reactions using single-step chemistry. However, few researchers had taken up this issue and presented multi-step reaction models. It was concluded that use multi-step kinetics is essential for accurate predictions of the temperature distributions, energy release rates, and emissions. Single-step kinetics was shown to be adequate for predicting all the flame characteristics except the emissions for the very lean conditions under which equilibrium favors the more complete combustion process dictated by global chemistry. full mechanism (49 species and 227 elemental reactions), skeletal mechanism (26 species and 77 elemental reactions), 4-step reduced mechanism 9 species and 1-step global mechanism.

In the open literature, there are few articles concerning the interaction between a fuel spray and a PM. The PM under study was of high porosity with uniformly distributed spheres and with uniform distribution of cavities with equal mean pore size in the porous medium. The interaction between droplets during evaporation was neglected. Physical properties of fuel such as latent heat of evaporation and healing value were constant. The temperature distribution inside of the fuel droplet was uniform, but time-varying. Laminar isobaric flow consits of air, fuel droplets, gaseous combustible mixture and hot products of combustion. After evaporation the fuel mixes immediately with air to form a homogeneous combustible mixture. Mass fraction of the liquid was negligible, and the gas was optically transparent. Also, radiative heat transfer between the skeleton surfaces of the porous medium can be considered.

9. Examples

9.1. Mesh preparation

Prior to CFD simulation, computational mesh of cylinder was generated with Kiva-Prep (preprocessor for mesh generation for KIVA-3V main code). The geometry of a mesh is composed of one block. Fig. 5 shows the grid configuration of porous tube, About 300000 grids were generated for computational studies.

Figure 5.

Computational mesh for the CFD calculation

9.2. Initial and boundary conditions

The test section was a vertical quartz glass-tube with 1.3 m in length and 0.076 m in diameter, that was isolated from the environment. The test section was filled with a packed bed of 0.0056 m solid alumina spheres. For simulation a cylinder with 0.076 m in diameter and 0.60 m in length that filled with PM, was considered (Fig. 5). The boundary condition applied to the momentum and energy equation with the assumption of zero gradients for temperature of both phase of PM and for species transport through the downstream boundary. At the upstream boundary, the gas temperature is 300K, composition is premixed methane-air with equivalence ratio 0.15, and velocity is 0.43 m/s of the premixed reactants and zero gradient for solid phase, were specified. For initial temperature for both phase of PM experimental measured data was used. Fuel is methane, porosity of PM is 0.4. The laminar flow considered for simulation. For validation of numerical simulation, modified KIVA code was used for simulation of unsteady combustion is a cylindrical tube with the experiments of Zhdanok et al. Fig. 6 plots a comparison of computation results to the experimental results of Zhdanok at the time of 147 s, which shows that the computed speed of combustion wave agrees well with the same condition of the experimental results. It is seen that methane is completely consumed in flame front that has maximum temperature.

Figure 6.

Comparison of combustion wave propagation between CFD and experiment in time 147 s

9.3. Discussion

In Figs. 7-10, contours of methane, carbon dioxide, temperature in gas and solid phases of PM, in cross section of tube for several times (10, 50, and 100 second) are shown. In Fig. 7 mass fraction of methane is shown. With entering of methane-air to porous tube that has high initial temperature, approximately 1800 K, in a narrow zone near inlet location, the temperature of gas increases until it reaches to self-ignition temperature. Methane consumes in a narrow region that is thicker relative to conventional flame front. In Figs. 8a, b, c, it is seen that after 10, 50 and 100 second combustion is started in respectively at x = 0.8, 2, 6 cm from entrance of mixture. The value of mass fraction of methane in PM tube is between 0.002 and this Fig indicates that the flame front has arc-shape.

Fig. 8 shows mass fraction of carbon dioxide in different time after start of simulation. Mixture flows through the porous tube that has initially heated and combustion in narrow zone of high temperature takes place. It is seen that reactions occur around the flame front in PM and its thickness is about 0.4 cm, which is very thick in compare with flame front in normal combustion. Carbon-dioxide disperses in pre heat region of entering mixture by diffusion of CO₂ and disperses in post flame by flow motion. In Fig. 9 temperature distribution in gas phase of PM in different time is shown. Flame front is recognizable from its high temperature region. Maximum fluid temperature is about 1600 K. Also, because effect of solid phase of PM, part of heat release of combustion is absorbed by it and prevents from high temperature gradient in fluid. Energy is re-circulated to the unburned gas mixture through the heat combustion and radiation of the solid. Fig. 10 shows temperature distribution in solid phase. At initial condition, mixture temperature is 300 K. The inlet temperature of solid phase is higher than gas temperature and heat is transferred from solid to gas, so allow it to reach the ignition temperature. Maximum temperature in solid phase is about 1600 K. Then the gas delivers its energy to the solid. Also, due to high heat capacity of solid phase of PM, low-temperature gradient occurs in it.

Figure 7.

Mass fraction of methane in cross section x = 0 after a) t = 10 s b) t = 50 s c) t = 100 s

Figure 8.

Mass fraction of Carbon-dioxide in cross section x = 0 after a) t = 10 s b) t = 50 s c) t = 100 s

Figure 9.

Gas phase temperature in cross section x = 0 after a) t = 10 s b) t = 50 s c) t = 100 s

Figure 10.

Temperature of solid phase in cross section x = 0 after a) t = 10 s b) t = 50s c) t= 100 s

Figs. 11, 12 show results for distribution of temperature in center line of PM tube in both phase of PM (solid and fluid) for different times versus axial direction. Heat transport is related to thermal properties of the solid material and fluid property. Flame core transports heat to incoming methane-air mixture with conduction and radiation, that its temperature is 300 K. At time t = 10 s maximum temperature in gas phase is about 1700 K. After 10 s flame moves to right with constant speed and maximum temperature of about 1600 K. After this time equilibrium between conduction, convection and radiation, causes to no change in maximum value of combustion. At the end of the tube temperature in all cases is about 325 K. In solid phase due to preheating of inlet mixture, at t = 10 s in inlet solid phase temperature is about 1450 K, in t = 50 s and t = 100 s, upstream temperature is about 850 K and 580 K, respectively and this temperature finally reach to 300 K approximately. Fig. 13 shows distribution of methane mass fraction from 10 to 100 s. Inlet mass fraction of methane is 0.016. Decrease in mass fraction of methane shows location of flame front. Flame location after 10 s is in location 3.8 cm and its thickness is thickness of 4 cm, after 50 s, is about 1.6 cm with the thickness 4.3 cm, and after 100 s is about 3.3 cm with the thickness 2.4 cm. From this Fig. inferred in 50 s and 100 s after simulation variation in flame thickness is very low value and the CH₄ is almost completely consumed in this zone. Fig. 14 shows mass fraction of CO₂ value in axial direction. Carbon dioxide is produced during combustion and its mass fraction reaches to highest value in x = 2.4, 5.5, 7.1 cm respectively to 0.036, 0.035, 0.024 mass fraction after 10, 50 and 100 s from simulation. Fig. 15 shows mass fraction of CO value in axial direction. Carbon monoxide is produced during combustion with entering of mixture and as an intermediate species produced and consumes gradually in axial direction. CO concentration reaches its highest value near the flame front at x = 1.8, 4.1, 6 cm respectively after 10, 50 and 100 s and gradually decreases at x = 4.8, 7.2 and 9.1 cm to approximately zero. With completeness of combustion and it oxidizes slowly and converted to CO_2, because of enough accessible oxygen for converting CO to CO₂. Fig. 16 shows mass fraction of CH₃ value in axial direction. Methyl concentration reaches its highest value in x = 1.4, 3.6, 9.2 cm respectively 10, 50 and 100 s after simulation. But after 100s due to fluid flow, it disperses in porous tube and oscillation occurs in value of it.

10. Two applications of PM technology

10.1. Internal combustion engines

The major target for further development of the current IC engines is to reduce their harmful emissions to environment. The most important difficulty with existing IC engines that currently exists is non-homogeneity of mixture formation within the combustion chamber which is the cause heterogeneous heat release and high temperature gradient in combustion chamber which is the main source of excess emissions such as NO_x, unburned hydrocarbons (HC), carbon monoxide (CO), soot and suspended particles. At present, the IC engine exhaust gas emission could be reduced by catalyst, but these are costly, sensitive to fuel and with low efficiency. Another strategy has been initiated to avoid the temperature gradient in IC engines that is using homogeneous charge compression ignition (HCCI) engines but there still exist some challenges including, higher HC and CO emissions and control of ignition time and rate of heat release under variable engine operating condition. In such engines, lean mixture with high amount of exhaust gas recycled (EGR), could be used that increases amount of CO and soot. Also by increasing the load of HCCI engine, NO_x formation and fuel consumption increases. It means that low compression ratio should be used for these engines, while in reality, the compression ratio must be high enough that temperature near the end of compression process, lead to self-ignition of mixture with reasonable time delay. Therefore, these engines are suitable for low and medium loads and it better for high loads the engine can change mode of operation to compression ignition. Also direct fuel injection engines generally have some unresolved problems due to lack of homogeneity of mixture formation and combustion. Several other technologies have been used to reduce emissions in engines, such as electronically controlled high pressure fuel injection systems, variable valve timing, EGR but still in these methods still could not solve the problem completely under all engine operating conditions. Could there be any other homogeneous combustion in IC engines to meet all operational conditions (various load and speed)? The demand target may be possible with homogeneous mixture formation and a 3D-ignition of a homogeneous charge to prevent formation of flame front that lead to temperature gradient in the entire combustion chamber which is ensuring a homogeneous temperature field. In conventional direct injection engines mechanisms also there is a lack of homogenization of combustion process. PM-engine is defined as an engine with homogeneous combustion process. The following distinct process of PM-engine is realized in PM volume: energy recirculation in cycle, fuel injection in PM, fuel vaporization for liquid fuels, perfect mixing with air, homogeneity of charge, 3D-thermal self-ignition, and homogeneous combustion. PM-engine may be classified as heat recuperation timing in an engine as: Engine with periodic contact between PM and cylinder is called closed PM-chamber and Engine with permanent contact between PM and cylinder which is called open PM-chamber. In this paper an open PM-chamber is studied. Permanent contact between working gas and PM-volume is shown schematically in Fig. 17. The PM is placed in cylinder head. During the intake stroke there is a not much influence from PM-heat capacitor with in-cylinder air thermodynamic conditions. Also during early stages of compression stroke only a small amount of air is in contact with hot porous medium. The heat transfer process (non-isentropic compression) increases during compression, and at TDC the air penetration is cut to the PM volume. At final stages of compression stroke the fuel is injected into PM volume and with liquid fuels rapid fuel vaporization and mixing with air occurs in 3D-structure of PM-volume. A 3D-thermal self-ignition in PM-volume together with a volumetric combustion is characterized by a homogeneous temperature distribution. Therefore, all essential conditions exist for having homogeneous combustion in the PM engine. The initial idea to use PM in IC engines was proposed by Weclas. Their investigation was performed in a single-cylinder air-cooled PM Diesel engine without any catalyst. A Silicon Carbide (SiC) PM was mounted in the cylinder head between the intake and exhaust valves and fuel was injected through the PM volume. The implementation has improved engine thermal efficiency, reduced emissions and noise in comparison to the base engine. The mean cylinder temperature was about 2200 K for base engine without having any PM. The temperature reduces to about 1500 K when PM is used which is significantly even lower during combustion.

Figure 11.

Mean temperature distribution in gas phase of PM versus axial direction

Figure 12.

Mean mass fraction of methane versus axial direction

Figure 13.

Mass fraction of methane versus axial direction

Figure 14.

Mass fraction of Carbone dioxide versus axial direction

Figure 15.

Mass fraction of CO versus axial direction

Figure 16.

Mass fraction of CH₃ value versus axial direction

Figure 17.

A permanent contact PM-engine in operation (Weclas, 2001)

The effect of SiC PM as a regenerator was simulated by Park and Kaviany. In their study a PM disk like shape was connected to a rod and was moving near piston within the cylinder of diesel engine. A two-zone thermodynamic model with single-step reaction for methane-air combustion is carried out. It is shown that the maximum cylinder pressure during combustion increases and more work is done during a full cycle, also engine efficiency increases, but due to high temperature of PM that its temperature is higher than adiabatic flame temperature of methane-air, the production of NO_x is rather higher while as its soot decreases. Macek and Polasek simulated and studied a PM engine with methane and hydrogen respectively and its potential for practical application was shown. Weclas and Faltermeier investigated penetration of liquid-fuel injection into a PM (as arrangement of cylinders which were mounted on a flat plate with different diameters). The arrangement was changed to obtain optimum geometry which produces the best mixture formation.

Porous regenerator as shown in Fig. 18 has the potential to improve fuel–air mixing and combustion. The porous insert is attached to a rod and moves in the cylinder, synchronized, but out of phase with the piston. During the regenerative heating stroke, the regenerator remains just beneath the cylinder head for most of the period and moves down to the piston (as it approaches the TDC position). During the regenerative cooling stroke, the regenerator moves up and remains in the original position until the next regenerative heating stroke. Following the combustion and expansion, the products of combustion (exhaust gases) retain an appreciable sensible heat. During the regenerative cooling stroke, the hot exhaust gas flows through the insert and stores part of this sensible heat by surface-convection heat transfer in the porous insert (with large surface area). For the proposed engine, a thermal efficiency of 53% was claimed, compared to 43% of the conventional Diesel engines. Macek and Polasek presented a finite volume based simulation of porous medium combustion for reducing emissions from reciprocating internal combustion engines.

Figure 18.

Sequence of motion of the regenerator and physical of fuel injection and air blowing during the regenerative heating stroke (Park and Kaviany, 2002)

The application of a highly porous open cell structures to internal combustion engines for supporting mixture formation and combustion processes was introduced by Weclas. Novel concepts for internal combustion engines based on the application of PMC technology were presented and discussed. His study proved that gas flow, fuel injection and its spatial distribution, vaporization, mixture homogenization; ignition and combustion could be controlled or positively influenced with the use of porous media reactors. The key features of the highly porous medium for supporting the mixture formation, ignition and combustion in IC engines are illustrated in Fig. 19. A study on the use of PMC in direct injection (diesel or gasoline) IC engines was performed by Durst and Weclas. Polasek and Macek presented the simulation of properties of IC engine equipped with a PM to homogenize and stabilize the combustion of CI engines. The purpose of the PIM matrix use was to ensure reliable ignition of lean mixture and to limit maximum in-cylinder temperature during combustion.

Figure 19.

Main feature of porous structure to be utilized to support engine process (Durst and Weclas, 2001)

10.2. Gas turbines and propulsion

A porous burner with matrix stabilized combustion for gas turbine applications. A numerical model for analyzing the evaporation processes in PM for gas turbine applications had been developed. Evaporation of a point wise-injected kerosene spray in a carbon-carbon porous medium was considered. The effects of porous medium temperature, fuel flow rate, air inlet temperature and porous medium geometry on the evaporation of spray can be analyzed. Evaporation characteristics were not found to vary much with porous medium geometry, as the porous medium was modeled as a momentum sink. But thermal effects of PM were found to be more dominant. The characteristics of combustion within porous media which are attractive in a propulsion context are the ability to burn leaner and hotter than a free flame with low emissions, there no cooling requirement for the combustor itself and the potential to operate free from combustion-induced noise. The performance of a PM combustor is applicable for gas turbines, at elevated pressures and inlet temperatures. The combustor was formed of reticulated porous ceramics, untreated to augment or sustain chemical reaction. The results showed that the combustor could operate in a ‘‘super-adiabatic” mode, with low emissions.