**3. PM modeling**

### **3.1 Gas phase energy equation**

Combustion inside a chamber is considered as thermodynamically modeling. Hence, gas phase energy equation and reaction rate are solved simultaneously. The closed system first law of thermodynamics (energy equation) for a differential crank angle, *dθ*, is [10, 11]:

$$\frac{dQ}{d\theta} - \frac{dV}{d\theta} = m\frac{du}{dt} \tag{9}$$

The definition of enthalpy *h = u + pv* and differentiating in constant pressure process gives:

$$\frac{dQ}{d\theta} = m \left. \frac{dh}{d\theta} \right| \tag{10}$$

The enthalpy in terms of chemical compositions of gases is:

$$h = \frac{\sum\_{i=1}^{N} N\_i \overline{h\_i}}{m} \tag{11}$$

where *Ni* and *hi* are the moles and molar enthalpy of species i, so:

$$\frac{dh}{d\theta} = \frac{1}{m} \left[ \sum\_{i} \left( \overline{h}\_{i} \frac{dN\_{i}}{d\theta} \right) + \sum\_{i} \left( N\_{i} \frac{d\overline{h}\_{i}}{d\theta} \right) \right] \tag{12}$$

Assuming ideal gas relation:

$$\frac{d\overline{h}\_i}{d\theta} = \frac{\partial \overline{u}\_i}{\partial T} \bullet \frac{dT}{d\theta} = \overline{c}\_{p,i} \frac{dT}{d\theta} \tag{13}$$

where *cp*,*<sup>i</sup>* is the constant molar pressure specific heat of species i. Concentration of species i is obtained as:

$$\left[\mathbf{X}\_{i}\right] = \frac{\mathbf{x}\_{i}}{R} \frac{P}{T\_{\mathcal{g}}} \tag{14}$$

The porosity of PM is φ (0 <φ<1), so the volume of fluid in PM is φ*V*, and the volume of solid is (1-φ) *V.* So *Ni* ¼ *φV Xi* ½ �, with differentiating of this equation:

$$\frac{dN\_i}{d\theta} = \rho V \dot{\alpha}\_i \tag{15}$$

*Mathematical Modeling of a Porous Medium in Diesel Engines DOI: http://dx.doi.org/10.5772/intechopen.108626*

The *ω*\_ *<sup>i</sup>* is the rate of reaction species i.

Gas-phase energy equation by applying PM is determined based on a relation of convective heat transfer between the gas and solid phases of PM [10, 11]:

$$\frac{dT\_{\rm g}}{d\theta} = \frac{-h\_{\rm v}(T\_{\rm g} - T\_{\rm s}) - \varrho \sum\_{i} (\dot{\alpha}\_{i} h\_{i}) + \varrho \mathop{\rm R\,}\nolimits\_{\rm g} \sum\_{i} \dot{\alpha}\_{i}}{\sum\_{i} \Big[ [X\_{i}] (c\_{p,i} - R) \big]} \tag{16}$$

*Q*\_ *=V* � � is heat loss to the environment.

The empirical relation Eq. (17) is applied to compute the volumetric heat transfer coefficient [1]:

$$h\_v = \frac{k\_\text{g} N u\_v}{d^2} \tag{17}$$

Volumetric Nusselt number is calculated by relation Eq. (18) [10, 11]:

$$Nu\_v = 2 + 1.11 \ Re^{0.60} Pr^{0.33} \tag{18}$$

Reynolds and Prandtl's numbers are computed by time-varying conditions of the fluid inside the combustion chamber. Reynolds was calculated according to injection velocity. Eq. (19) is used for the calculation of the Reynolds number [10, 11]:

$$Re = \frac{\rho\_g u\_c \, d}{\mu} \tag{19}$$

The density of gases from the ideal gas equation:

$$\rho = \frac{P}{T \left( \mathbb{X}\_{\!\!\!\!\!M \!\!W\_{\text{mix}}} \right)} \tag{20}$$
 
$$\mathbf{M} \mathbf{W}\_{\text{mix}} = \sum\_{i}^{N} \mathbf{x}\_{i} \mathbf{M} \mathbf{W}\_{i}$$

Relation Eq. (21) applies for calculating velocity where ρ<sup>f</sup> is the fuel density, the chamber pressure is pc, and injection pressure is clarified by pi [10, 11].

$$
\mu\_c = \sqrt{\frac{2\left(p\_i - p\_c\right)}{\rho\_f}}\tag{21}
$$

Eq. (22) shows the calculation of chamber pressure [10, 11]:

$$P = \sum\_{i} [X\_i] \, R \, T\_{\text{g}} \tag{22}$$

The modified Arrhenius equation is used for modeling the fuel oxidation as shown in Eq. (23) [1].

$$\dot{m}\_{\text{Fuel}} = \frac{d[\mathbf{X}\_{\text{Fuel}}]}{d\theta} = -\mathbf{A} \cdot \exp\left(-\mathbb{E}\_{\text{V}}\boldsymbol{\zeta}\_{\text{RT}}\right) [\mathbf{X}\_{\text{Fuel}}]^m \left[\mathbf{X}\_{\text{Oxygen}}\right]^n \tag{23}$$

A, m, and n are constant coefficients of single-step fuel oxidation Eq. (23), selected as **Table 1**. After solving coupling Eqs. (16) and (23), the gas temperature, species, and fuel consumption rates are determined. After updating the temperature and concentration of species, the pressure was calculated in Eq. (24) [1].

$$P = \sum\_{i} [\mathbf{X}\_i] \, R \, T \tag{24}$$

NASA seven-term polynomials are applied to compute the thermodynamic properties. For the calculation of constant pressure of specific heat, NASA polynomial can be used, which depends on temperature:

$$\mathcal{L}\_p = \mathbb{R}\left(a\_1 + a\_2T + a\_3T^2 + a\_4T^3 + a\_5T^4\right) \tag{25}$$

Enthalpy of gas can be determined from [10, 11]:

$$h = R\left(a\_1 + \frac{a\_2T}{2} + \frac{a\_3T^2}{3} + \frac{a\_4T^3}{4} + \frac{a\_5T^4}{5} + \frac{a\_6}{T}\right)\tag{26}$$

$$h(T) = \Delta l \mathfrak{y}(\mathbf{298}) + [h(T) - h(\mathbf{298})]$$

#### **3.2 Solid phase energy equation**

Energy equation for the solid phase of PM is:

$$
\dot{Q} = m \frac{dh}{d\theta} \tag{27}
$$

By chain rule differentiable:

$$\frac{dh}{d\theta} = \frac{\partial h}{\partial T\_s} \frac{dT\_s}{dt} = c\_s \frac{dT}{d\theta} \tag{28}$$

Solid phase volume is (1-φ) V, so:

$$\frac{\dot{Q}}{(1-\rho)\ V\rho\_s} = c\_s \frac{dT\_s}{d\theta} \tag{29}$$

where ρ<sup>s</sup> is the density of solid phase, and *cs* is the specific heat of solid phase.


**Table 1.**

*Single-step oxidation of several fuels [10, 11].*

Due to heat transfer between energy gas and solid phase, *hv Tg* � *Ts* is the heat transfer from the gas phase to solid phase. Also, *q*\_*<sup>r</sup>* is radiation from solid phase to environment. Therefore, Eq. (30) calculates solid temperature.

$$\frac{dT\_s}{d\theta} = \frac{h\_v \left(T\_g - T\_s\right) - \dot{q}\_r}{(1 - \rho)\rho\_s \mathbf{C}\_s} \tag{30}$$

The solid phase energy is computed based on Eq. (31) [10, 11]:

$$\frac{dT\_s}{d\theta} = \frac{h\_v \left(T\_g - T\_s\right) - \varepsilon \sigma \frac{A}{V} \left(T\_s^4 - T\_0^4\right)}{(1 - \rho)\rho\_s \mathbf{C}\_s} \tag{31}$$

where *<sup>A</sup> V <sup>s</sup>* is the surface area to volume ratio of the PM, and *T0* is the environment's temperature where radiation loss occurs.

#### **3.3 Solution method of the equations**

Three coupled Eqs. (32)–(34) that are relevant to solid and gaseous phase temperature and species concentrations are solved.

$$\frac{dT\_{\rm g}}{d\theta} = f(\theta, [X\_i], T\_{\rm g}, T\_s) \tag{32}$$

$$\frac{dT\_s}{d\theta} = \mathbf{g}\left(\theta, [\mathbf{X}\_i], T\_\mathbf{g}, T\_s\right) \tag{33}$$

$$\frac{d[X\_i]}{d\theta} = h\left(\theta, [X\_i], T\_{\rm g}\right) \tag{34}$$

Runge–Kutta fourth-order method is applied to solve the fluid and solid phase energy equations in PM and compute the species concentrations. Due to the high rate of combustion, the small crank angle (Δθ) should be considered. The step size for solving Runge–Kutta method was assumed 10�<sup>6</sup> s or 0.1 crank angle that depend on engine speed [10, 11].

$$k\_{1\xi} = hf\left(\theta\_n, T\_{\xi}, T\_s\right)$$

$$k\_{1,\varepsilon} = h\,\mathcal{g}\left(\theta\_n, T\_{\mathcal{g}}, T\_s\right)$$

$$k\_{2\mathcal{g}} = hf\left(\theta\_n + \frac{h}{2}, T\_{\mathcal{g}} + \frac{k\_{1\mathcal{g}}}{2}, T\_s + \frac{k\_{1\mathcal{g}}}{2}\right)$$

$$k\_{2,\varepsilon} = h\,\mathcal{g}\left(\theta\_n + \frac{h}{2}, T\_{\mathcal{g}} + \frac{k\_{1\mathcal{g}}}{2}, T\_s + \frac{k\_{1\mathcal{s}}}{2}\right)$$

$$k\_{3\mathcal{g}} = hf\left(\theta\_n + \frac{h}{2}, T\_{\mathcal{g}} + \frac{k\_{2\mathcal{g}}}{2}, T\_s + \frac{k\_{2\mathcal{s}}}{2}\right)$$

$$k\_{3,\varepsilon} = h\,\mathcal{g}\left(\theta\_n + \frac{h}{2}, T\_{\mathcal{g}} + \frac{k\_{2\mathcal{g}}}{2}, T\_s + \frac{k\_{2\mathcal{s}}}{2}\right)$$

$$k\_{4,\mathcal{g}} = hf\left(\theta\_n + h, T\_{\mathcal{g}} + k\_{3\mathcal{g}}, T\_s + k\_{3\mathcal{s}}\right)$$

$$k\_{4,\mathcal{s}} = h\,\mathcal{g}\left(\theta\_n + h, T\_{\mathcal{g}} + k\_{3\mathcal{g}}, T\_s + k\_{3\mathcal{s}}\right)$$

**189**

Updated gas and solid phase temperatures of PM can be obtained from Eqs. (35) and (36). An update of species concentration is computed from Eq. (37):

$$\theta\_{n+1} = \theta\_n + \Delta\theta$$

$$T\_{\mathfrak{g}, n+1} = T\_{\mathfrak{g}, n} + \frac{1}{6} \left(k\_{1\mathfrak{g}} + 2k\_{2\mathfrak{g}} + 2k\_{3\mathfrak{g}} + k\_{4\mathfrak{g}}\right) \tag{35}$$

$$T\_{s,n+1} = T\_{s,n} + \frac{1}{6} \left(k\_{1\circ} + 2k\_{2\circ} + 2k\_{3\circ} + k\_{4\circ}\right) \tag{36}$$

$$\left[\left[X\_{i}\right]\_{n+1} = \left[X\_{i}\right]\_{n} + h\left(\theta\_{n}, \left[X\_{i}\right]\_{n}, T\_{\lg n+1}\right) \tag{37}$$
