Comprehensive Correlation for the Prediction of the Heat Release Characteristics of Diesel/CNG Mixtures in a Single-Zone Combustion Model

Abstract

Fuel combinations with substantial differences in reactivity, such as diesel/CNG, represent one of the most promising alternative combustion strategies these days. In general, the conversion from diesel to dual-fuel operation can be performed in existing in-use heavy-duty compression-ignition engines with minimum modifications, which guarantee very little particles, less nitrogen oxide (NOx), and reduced noise by half compared to diesel. These factors make it feasible to retrofit a CNG fuel system on an existing diesel engine to operate it in dual fuel mode. However, the single-zone combustion models using the traditional single-Wiebe function are exceptionally adopted to assess the dedicated dual fuel engines, whereas the heat loss to the walls is estimated by using the Woschni heat loss formulation. It means that the fast and preliminary analysis of the unmodified engine performance by 1-zone models becomes complicated due to the obvious deterioration of the energy parameters, which, in turn, was predetermined from the deviation in the thermodynamic cycle variables as the calculation outcome. In this study, the main novelty lies in the fact that we propose a novel composition-considered Woschni correlation for the prediction of the heat release duration characteristics of diesel/CNG mixtures for the unmodified diesel engine. The elimination of former deficiencies distinctive to a single-zone thermodynamic model by applying the interim steps described became the core of the research presented in this paper. It led to successful derivation of the necessary correlation for modelling the heat release duration characteristics of an ICE operated in the dual fuel mode.

1. Introduction

The ever-increasing emissions of CO₂ raise massive efforts from industrial experts and researchers around the world who are trying to overcome impending challenges and taking actions toward faster decarbonization of industry and transport sectors [1,2,3]. In recent years, diesel engine emissions were a target for alleviation, to a certain extent, through the research of experts and scholars [4,5,6]. For the road transport sector, a time span between 2022 and 2035 marks an uncertainty on how to ensure lower levels of GHG emissions if (i) in the case of a pessimistic scenario, the electrification of long-distance trucks will fail, and (ii) the large numbers of long-distance trucks driven by ICE will be left abandoned due to the era of electrified transport [7]. For the medium-term scenario, natural gas (CNG/LNG) vehicles have been promoted worldwide to mitigate air pollution [3,8,9,10]. Furthermore, there are several vehicle types that can run on two to five energy carriers (bioethanol E85, 10% hydrogen and 90% methane, biomethane, natural gas (CNG), and gasoline (The Volvo Multi-Fuel)), referred to as multi-fuel vehicles, including bi-fuel vehicles, flexible fuel (or dual fuel) vehicles, and plug-in hybrid electric vehicles [11]. A study by Armaroli et al. (2022) [2] has found that natural gas, biogas, and biomethane are ranked number two (after electricity) among alternative fuels in Italy and could trigger the widespread adoption of alternative fuel vehicles (all motorcycles, passenger cars and light-duty vehicles, short- and long-haul buses, short- and long-haul trucks, short- and long-haul trains, and short- and long-haul ship vessels) to support energy policy. Fuel combinations with substantial differences in reactivity, such as diesel/CNG, represent one of the most promising alternative combustion strategies these days [12]. Dual fuel combustion consists of the preparation of a premixed fuel (CNG) and intake air, whose ignition is triggered by the injection of a more ignitable diesel fuel [13]. However, direct conversion of a diesel engine is not possible due to the very low Cetane number of CNG, which means that the extra cost will be reimbursed by the savings in operating cost due to fuel cost [14]. These minor aspects can be compensated for with the potential of CNG to meet strict engine emission regulations and be cheaper than other fuels in many countries [15]. To accelerate this technology into a real application for the next two decades, especially trying to evaluate a wide variety of engine families in the used truck or second-hand vehicle market, theoretical validation of experimental results through mathematical modeling helps in modeling engine combustion and performance parameters and reduces experimental engine testing cost [16]. The alternatives to using the thermodynamic combustion model (TCM) or computational fluid dynamic (CFD) model are based on the main aim of modeling and simulation. Major efforts are focused on the development of computational fluid dynamic (CFD) technology, which performs the three-dimensional modeling of gas flow and fuel injected into the manifolds and cylinders of IC engines. The most popular softwares are KIVA (Los Alamos) [17], FIRE (AVL) [18], VECTIS (Ricardo) [19], and STAR-CD (Computational Dynamic Ltd., Berwyn, IL, USA) [20]. The programs allow modeling the movements of gas and burnt fuel flows in the combustion chamber, as well as temperature and concentration fields during the combustion process. Parallel to the methodology used in CFD, software is being developed that analyzes the thermodynamic process in one-dimensional or dimensionless space. The programs BOOST (AVL) [19], WAVE (Ricardo) [19], GT-Power (Gamma Technologies) [21], and HIROYASU [22] are considered world leaders here. Among the programs developed in Russia, IMPULS [23] and DIZEL-2/4t software have considerable commercial success. The competing programs BOOST, WAVE, and GT-Power provide the user with a fairly wide range of design options: a convenient control system, the ability to analyze one-dimensional gas exchange models, which allow one to evaluate the flow in the fuel tube system, to evaluate the unevenness of the cylinder filling, and to optimize gas distribution phases. To calculate the formation and combustion of a combustible mixture, these programs use methodologies based on equations derived in 1962 by Ivan Ivanovitch Wiebe (1902–1969) [24] or by subsequent works of scientists who used very similar research principles. The biggest challenge is modeling the formation and combustion processes of the combustible mixture. The later version of the IMPULS software uses the calculation method proposed by Razliaitsev [25] in 1980. This method made it possible to study the dependence of the heat release rate on parameters, such as compression ratio, injection delay, diameter and number of nozzle holes, and injection characteristics. Thermodynamic combustion models, especially the single-zone (1-z) combustion model (the simplest TCM) cannot provide analysis with the same prediction level as CFD. Even when the cylinder charge is assumed to be uniform in density and the space inside the cylinder is regarded as a single zone, the model is capable of predicting the combustion process under various operating conditions of the engine with the required level of precision [26]. According to Stepanenko and Kneba (2019) [26], single-zone models are often used if there is a need for a fast and preliminary analysis of engine performance. In other words, 1-z models could be considered a thermodynamic system in which combustion takes place in the engine, which exchanges heat (energy) and mass transfer with the environment, and thermodynamics laws could be applied to determine heat transfer from the system [16].

Based on the insights of Gautam et al. (2022) [16], there are limited studies on the combustion behavior of diesel engines using thermodynamic modeling, and also there have been no studies available on the detailed statistical analysis between experimental and numerical data. To bridge these gaps, our study presents the experimentally derived Woschni correlations for the assessment of the heat release duration of diesel/CNG blends using a single-zone thermodynamic model. This study will exclusively attempt to find a solution on how to expand the boundaries of a single-zone thermodynamic model to be able to address the challenges faced by today’s long-haul transport, by addressing the following research questions:

(1) What is the actual range of calculation error when a single-zone model is applied for the assessment of the unmodified diesel engine running in dual fuel mode? Are they acceptable or should further strategies be undertaken?

(2) What is the nature of the overall relationship between the indicated thermal efficiency values and the excess air coefficient for a wide range of indicated mean effective pressure values and injection angles for a wide spectrum of D/CNG ratios? If clear trends are observed, how does the experimental outcome suite the Vibe’s predictive staged correlation, which represents the gross energy release of fuel as a function of the crank angle?

(3) Which parameter has a primary influence on the duration of the heat release and can be labelled as an example of a deterministic relationship? Can the generalized equations for all investigated D/CNG ratios be obtained following our methodology?

(4) Is the accuracy of the newly derived Woschni correlation satisfactory if compared with the results obtained by other researchers?

2. Materials and Methods

The rationale and methodology of the study were considered before the investigation was carried out. As a result, the study employs an extended experimental investigation of a diesel engine running on diesel/CNG blends (see Figure A1 and

Table A1. Engine specification.

Parameter	Description
Displacement (cm3)	1896
Bore × stroke (mm)	79.5 × 95.5
Maximum power (kW/rpm)	66/4000
Maximum torque (Nm/rpm)	180/2000–2500
Cooling type	Water-cooled
Fuel supply system	Direct injection
Cylinders	Inline-four
Compression ratio	19.5:1
Aspiration	Turbocharger

in Appendix A) to obtain a data set of parameters required for the theoretical analysis of Woschni correlations applied for combustion and heat transfer, respectively, in a single-zone model, and its further validation.

2.1. An Assessment of the Combustion Characteristics of Compression Ignition Engine Using a Single-Zone Thermodynamic Model

The single-zone mathematical model of basic software IMPULS comprised of 18 sub-models was used for the simulation of the engine energy parameters. This software has previously been successfully used in the development and modification of high-speed transport engines (15/15, 15/18, 16.5/18) [27]. A single-zone thermodynamic model (a closed system that undergoes various changes due to temperature, pressure, and volume; however, its final and initial states are equal) was calibrated to assess the diesel engine equipped with a turbocharger. The quasi-static approximation implies that the set of thermodynamic equations and the equations of gas dynamics were able to describe precisely the design parameters of the exhaust system, variable coefficients of efficiency for gas turbine and compressor unit, heat losses to the engine cooling system, etc. In-cylinder processes were described with the help of differential equations related to energy conservation and mass conservation in their general form [27]:

$$\frac{dU}{d\tau }=\frac{d{Q}_{re}}{d\tau }-\frac{d{Q}_e}{d\tau }-p·\frac{dV}{d\tau }+h_s·\frac{d{m}_s}{d\tau }-h_{ex}·\frac{d{m}_{ex}}{d\tau },[\mathrm{kJ}/\mathrm{s}]$$

(1)

$$\frac{dm}{d\tau }=\frac{d{m}_s}{d\tau }+\frac{d{m}_{inj}}{d\tau }-\frac{d{m}_{ex}}{d\tau },[\mathrm{kg}/\mathrm{s}]$$

(2)

$$\frac{dp}{d\tau }=\frac{m·R}{V}·\frac{dT}{d\tau }+\frac{m·T}{V}·\frac{dR}{d\tau }+\frac{R·T}{V}·\frac{dm}{d\tau }-\frac{p}{V}·\frac{dV}{d\tau },[\mathrm{Pa}/\mathrm{s}]$$

(3)

where U is the internal energy of a thermodynamic system (J); Q_re and Q_ex are heat energy released and heat energy transferred, accordingly (J), respectively; p is the pressure (Pa); V is the volume (m³); h_s the amount of internal energy contained in a compound (enthalpy) (J/kg); h_ex the exhaust enthalpy (J/kg); m is the total mass (kg); m_s is the supply (intake) air mass (kg); m_inj is the mass of fuel sprayed (kg); m_ex is the mass of exhaust gases (kg); τ is the time (s); R is the gas constant (J/kg·K); and T is the temperature (K).

The Vibe function itself is parameterized by studying the heat release [28]; in its current form, the function includes the empirical observation made by Woschni (1967) [29]. The heat release was obtained from cylinder pressure data by methods described in [30,31,32]. In general, the empirical observation of Woschni improves the accuracy of the approximation of the actual heat release characteristics of the Vibe function through the full attribution of m (form factor) and φ_z (conditional duration of combustion) in a thorough evaluation of the partial loads of a diesel engine [27]. The in-cylinder heat transfer is governed by Gerhard Woschni’s equations separately for the piston, cylinder head, and cylinder bushings by taking into account the intensity of the in-cylinder air motion. The cyclic portion of a dual fuel blend (q_cycl) and the lower heating values (LHV) for both fuels were established by employing the Equations (4) and (5) [27]:

$$q_{cycl}=\frac{q_{cyc{l}_D}·LH{V}_D+q_{cyc{l}_{CNG}}·LH{V}_{CNG}}{LH{V}_{mix}}$$

(4)

where $q_{cycl}$ is the overall fuel consumption per cycle (g/cycle); $q_{cyc{l}_D}$ and $q_{cyc{l}_{CNG}}$ are the diesel fuel and CNG consumptions per cycle (g/cycle), respectively; $LH{V}_{mix}$ is the lower heating value of the D/CNG mixture; LHV_D and LHV_CNG are the lower heating values for diesel fuel and CNG (MJ/kg), respectively.

The lower heating values for D/CNG mixtures later were determined from Formula (5), as described in [27]:

$$LH{V}_{mix}=337.5 · C+1025 · H-108.3 · O$$

(5)

where

$$C=C_D · \left(100-CC{R}_{CNG}\right)+C_{CNG} · CC{R}_{CNG}$$

(5a)

$$H=H_D · \left(100-CC{R}_{CNG}\right)+H_{CNG}· CC{R}_{CNG}$$

(5b)

$$O=O_D · \left(100-CC{R}_{CNG}\right)+O_{CNG} · CC{R}_{CNG}$$

(5c)

$$CC{R}_{CNG}=\frac{q_{cyc{l}_{CNG}} · LH{V}_{CNG}}{q_{cyc{l}_{CNG}} · LH{V}_{CNG}+q_{cyc{l}_D} · LH{V}_D} · 100\%$$

(5d)

Here, $CC{R}_{CNG}$ is the percentage of CNG in the combustible mixture (co-combustion ratio) (%), and $C_D$ and $C_{CNG}$ are the numbers of carbon atoms in diesel and CNG, respectively.

The indicated efficiency (${\eta }_i$) was established using the following formula [27]:

$${\eta }_i=\frac{3.6·p_{mi}}{LH{V}_{LD} · G_{fD}+LH{V}_{CNG} · G_{fCNG}}$$

(6)

where $G_{fD} \mathrm{and} G_{fCNG}$ are the hourly fuel consumption for diesel fuel and CNG consumption, respectively (kg/h); p_mi is the mean effective pressure (kPa).

2.2. Fuel Specifications

During the experiment, two fuel types were used: diesel fuel (EN 590) and CNG (ISO 6976:1995) (see

Table 1. Fuel properties.

Fuel Property	CNG	Diesel Fuel
Density (kg/m3)	0.74	829.0
Cetane number	-	49
Lower heating value (MJ/kg)	51.7	42.8
Viscosity (cSt 40 °C)	-	1.485
H/C ratio	-	1.907
Component (% vol.)	Methane: 91.97	Carbon: 86.0
	Ethane: 5.75	Hydrogen: 13.6
	Propane: 1.30	Oxygen: 0.4
	Butane: 0.281
	Nitrogen: 0.562 Carbon dioxide: 0.0

) [27].

2.3. Theoretical Background

Determination of the amount of heat energy generated during the combustion of the fuel charge Q_i, based on the defined pressure inside the cylinder p(φ), which is obtained from experimental data. It is also needed to define the heat transfer through the walls Q_W(φ); if only one pressure is defined then only the amount of useful heat Q can be determined, according to the Equation (7):

dQ = dQ_i − dQ_W

(7)

The amount of useful heat coupled with the amount of heat being transferred through the walls can be considered as two summands separately assessed by the IMPULS software to obtain the heat release characteristics Q_i or (Q_i = f(φ) of an engine:

dQ_i = dQ + dQ_W

(8)

The Q_i input values are determined experimentally from indicator diagrams which are used to assess the performance of an ICE; it can be expressed as:

$$Q_i=\frac{3600·P_{mi}}{LHV·G_f}$$

(9)

Opposite to Q_i values, for the establishment of heat transfer through the walls, the mathematical model uses the installed library to render math equations, which supports a large subset of IMPULS functions. The lack of calculation precision for the Q_W values leads to systematic error in the accurate establishment of the heat release characteristics when D/CNG is used in dual fuel mode. From the earlier study published by the coauthors of this paper [27], it was found that the use of a single-zone model for the detailed assessment of an ICE in dual fuel mode gives a calculation error of 25% for the establishment of the heat balance characteristic Q_i = f(φ), compared to the manual input of experimentally derived Q_W values. Furthermore, the difference between the calculated and experimentally established φ_z (duration of heat release) values may exceed 307%, while the indicators of the diesel engine (p_me, η_e, etc.) may differ more than two times [27,33,34,35]. Eliminating these deficiencies through the applyication of the interim steps described below became the core of the research presented in this paper. As already mentioned, the latest version of IMPULS software uses an approximation that heat transfer through the walls and Q_W is a single constituent of energy losses [27,33,34,35]. A thorough analysis of the use of D/CNG as fuel for ICE in dual fuel mode was extensively published in Refs. [33,34,35], suggesting that the total energy losses of heat energy generated during the combustion (Q_total) should be assessed as a sum of energy losses due to valves and fittings (friction losses), Q_fr and heat transfer through the walls Q_W:

$$Q_{total}=Q_{fr}+Q_W$$

(10)

$$Q_{fr}=\frac{p_{mm}·3600}{Q_i}$$

(11)

where $p_{mm}$ is the average mechanical pressure exerted at a point in the fluid that can be established from the mathematical relationship ${\eta }_m=\frac{p_{me}}{p_{mi}}=\frac{p_{mi}-p_{mm}}{p_{mi}}$) ⇒ $p_{mm}=\frac{p_{me}·\left(1-{\eta }_m\right)}{{\eta }_m}$; ${\eta }_m$ is mechanical efficiency.

The mean indicative pressure $p_{mi}$ that an engine can put out may be derived from a simple formula for the efficiency of a heat engine (${\eta }_e$):

$${\eta }_e=\frac{p_e · \left(1-{\eta }_m\right)}{{\eta }_m}$$

(12)

As a result, an assessment of energy losses due to valves and fittings can be performed using Equation (13):

$$Q_{fr}=\frac{{\eta }_e\left(1-{\eta }_m\right)}{{\eta }_m}$$

(13)

where: ${\eta }_m$ is the coefficient of performance (${\eta }_m=0.727$ for p_mi = 8.2 bar, ${\eta }_m=0.637$ for p_mi = 2.2 bar, and ${\eta }_m=0.474$ for p_mi = 4.2 bar).

Q_total, Q_fr, and Q_W values for an ICE were evaluated when D/CNG is used in dual fuel mode—see

Table 2. Total energy losses of heat energy generated during combustion (Q_total), %.

	D100		D60/CNG40		D40/CNG60		D80/CNG20
	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC
pmi = 4.2 bar	35.3	35.5	45.5	43.2	51.6	53.3	59.7	61.4
pmi = 6.2 bar	31.1	32.3	35.0	35.7	32.9	38.0	46.2	44.6
pmi = 8.2 bar	27.1	28.7	30.0	32.6	32.6	33.6	33.2	33.6

,

Table 3. Energy losses due to valves and fittings (friction losses) (Q_fr), %.

	D100		D60/CNG40		D40/CNG60		D80/CNG20
	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC
pmi = 4.2 bar	27.2	27.9	22.3	24.6	19.2	22.3	14.8	17.2
pmi = 6.2 bar	17.8	18.4	15.5	17.6	14.8	16.9	12.8	15.6
pmi = 8.2 bar	12.7	13.5	12.1	13.3	11.7	13.0	12.1	13.0

and

Table 4. Energy losses due to heat transfer through the walls (Q_W), %.

	D100		D60/CNG40		D40/CNG60		D80/CNG20
	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC
pmi = 4.2 bar	8.1	7.6	23.2	18.6	32.4	31.0	44.9	44.2
pmi = 6.2 bar	13.3	13.9	19.5	18.1	18.1	21.1	33.4	29.0
pmi = 8.2 bar	14.4	15.2	17.9	19.3	20.9	20.6	21.1	20.6

, respectively. Two different injection timings (φ_inj = −1 degCA BTDC and (φ_inj = −13 degCA BTDC) and three different load modes: (i) BEMP = 6 bar (p_mi = 8.2 bar), high load mode, (ii) BEMP = 4 bar (p_mi = 6.2 bar), medium load mode, and (iii) BEMP = 2 bar (p_mi = 4.2 bar), low load mode were a focus of research in this study.

The results revealed that the total energy losses of heat energy generated during the combustion tend to increase drastically when an engine is switched from diesel fuel to run in a dual fuel mode (see Table 2): +9.2%/+3.9%/+1.6% for the D60/CNG40 blend, +16.3%/+1.8%/+5.5% for the D40/CNG60 blend, and +24.4%/+15.1%+6.1% for the D20/CNG80 blend, respectively. The largest energy losses were observed for low loads (p_mi = 4.2 bar).

In other words, an increase in %CNG resulted in a linear growth of Q_total mainly associated with reduced friction losses (see Table 3) and sudden jump of Q_W values at p_mi = 4.2 bar (see Table 4): +186.4% for D60/CNG40 blend, +300.0% for D40/CNG60 blend, and +445.7% for D20/CNG80 blend, respectively. An interesting observation is that the retardation in the injection timing from φ_inj = −1 degCA relative to the TDC to φ_inj = −13 degCA relative to the TDC has a positive influence on the reduction of Q_W, especially at p_mi = 4.2 bar (see Table 4): −19.8% for the D60/CNG40 blend, −4.3% for the D40/CNG60 blend, and −1.6% for the D20/CNG80 blend, accordingly. Thermal balance analysis is a useful method to determine the energy distribution and efficiency of the compression ignition engine which is introduced as a specific piece of IMPULS software. A general mathematical model of software estimates the energy balance using the steady-flow energy equation, as expressed in Equation (14):

$$E_{in}-E_{out}=dE/dt=0$$

(14)

where: $E_{in}$ is the input (fuel) energy; $E_{out}$ is the output energy.

In the steady-flow process, the energy content of a control volume remains constant, so the rate of change in total energy is equivalent to zero. Therefore, the energy balance can be rewritten as Equation (15):

$$E_{in}=E_{out}$$

(15)

The software model is developed based on the underlying assumption that the $E_{in}=E_{out}=1$, not depending on the type of fuel used.

Table 5. Juxtaposition of the Q_total (counter) and Q_{total_corr} (denominator) values.

	D100		D60/CNG40		D40/CNG60		D80/CNG20
	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC	φinj = −1 degCA BTDC	φinj = −13 degCA BTDC
pmi = 4.2 bar	1.10/1.01	0.96/0.99	0.92/0.98	0.94/0.98	0.84/0.99	0.85/0.99	0.75/0.97	0.86/0.99
pmi = 6.2 bar	1.02/1.00	0.99/0.99	0.95/0.98	0.96/0.99	0.83/0.97	0.91/0.98	0.87/0.97	0.87/0.98
pmi = 8.2 bar	1.01/0.99	0.99/1.00	0.87/1.00	0.99/1.01	0.98/1.01	0.99/1.01	0.93/0.99	0.96/1.00

represents the result of the calculation of the amount of heat energy released during the entire cycle of combustion in the form of approximated heat release characteristics before the correcting of heat losses (Q_total) and after the correcting of heat losses (Q_{total_corr}) by employing the data presented in Table 2, Table 3 and Table 4. A relative number greater than 1 indicates that the amount of heat energy released is greater than the energy stored in a fuel (chemical energy), and opposite, if Q_total or Q_{total_corr} < 1, this number suggests that the chemical energy of a fuel was higher compared to the output energy. Both cases show inadequacy in the calculations. What can be concluded from the analysis of Q_total values is that the most accurate results were obtained for the performance of the engine with diesel fuel at the high load mode (p_mi = 8.2 bar) when the compromise in the indicated thermal efficiency was negligible due to the absence of CNG in a fuel blend. The performance of a diesel engine in dual-fuel operation at lower loads has shown a significant deviation between the input energy and the heat energy released. The manual replacement of originally available energy loss values with the ones presented in Table 2, Table 3 and Table 4 was made using the guide on the parameter controls in the Adjustment section of the IMPULS software. After repeated calculations to compute formulas that refer back to the same algorithm, the corrected values (Q_{total_corr}) were derived to characterize the amount of heat energy released during the entire cycle of combustion with a calculation bias not exceeding 2%.

3. Results

3.1. Engine Coherency Analysis for Identification of the Aspects and Impacts of %CNG on a Combustion Cycle

Figure 1 presents the results obtained which reflect the two characteristic features of a reciprocal interaction between the ratio of the mean indicated pressure of a cycle to the indicated efficiency (η_i) and excess air coefficient (α) (it is worth mentioning that the use of mechanical efficiency η_m (work output/heat input) instead of η_i is not recommended for this case due to skewed outcome data). First, clear linear relationships were obtained between the indicated thermal efficiency values and α for a wide range of indicated mean effective pressure p_mi values (or alternatively, brake effective mean pressure p_me values) and injection angles for a wide spectrum of D/CNG ratios, including pure diesel fuel as a reference fuel [34]. Second, since a diesel engine is defined as one “that has operating characteristics significantly similar to those of the theoretical diesel combustion cycle”, it can be clearly seen from Figure 1 that the gradual increase in p_mi (p_me) values is accompanied by the parallel shift of an excess air coefficient to a lower range of distinct zones being characterized by lower η_i. Thus, not only the main influence of η_i on α is being ensured, but also the quantitative (albeit smaller) impact of η_i = f(α) on other dynamic indicators of a combustion cycle, namely the pressure at the end of the compression stroke, p_e, and the maximum pressure of the cycle, p_max ($\alpha =p_{max}/p_e$) is also being successfully secured. Parameters p_max and p_e depend mainly on the size of the diesel fuel portion injected at various retarded angles φ_inj (see Figure 1) in a form of liquid phase.

The findings revealed that the linear increase in %CNG causes the corresponding reduction in the indicated thermal efficiency. It is especially evident in the D20/CNG80 case, where the variation of the excess air coefficient has the shortest range and, in turn, after being shifted to the right side of the x-axis, is limited to the lowest values compared to other blends (see Figure 1). These facts testify to a deterioration of the combustion cycle that was determined by the different physicochemical properties of diesel fuel and CNG. From the process modelling point of view, these differences are attributed to the tangible change in combustion duration and shape factor m of the Vibe function which are coupled to each other by scaling factors. The reduction in η_i/α ratio with an increase of %CNG at p_mi = idem (idem: the same as the actual condition [36]) is predetermined by the disproportionate increase in the mass fraction of natural gas (G_CNG) entering the cylinder, which is much higher than the ratio of lower heating values for both fuels LHV_CNG/LHV_D = 1.18 ratio. This ratio is in good conformity with the relation of stoichiometry constants (see Equation (16)) for the same fuels, L_CNG/L_D = 1.19:

$$\alpha =\frac{G_{air}}{L_D·G_D+L_{CNG}·G_{CNG}}$$

(16)

where: $G_{air}$ is the air flow rate (kg/h).

The phenomenon of disproportionate increase in the mass fraction of CNG entering the cylinder is formulated as an outcome of rejecting the null hypothesis. First, both fuels are characterized by similar LHV values, which should not have an evident impact on a steep reduction in α when diesel fuel is replaced by diesel blends with high concentrations of CNG. Moreover, the minor change then compensates for the higher value of stoichiometry constant of CNG. Second, we rejected the hypothesis that when D/CNG is used in dual fuel mode, the reduction in excess air coefficient was associated with an observed decrease in air pressure after compression, as the alteration of p_k values (at p_mi = idem) was in the range of only 4–5%.

Adding CNG to a diesel fuel significantly changes the ability of the unmodified ICE to release chemical energy through combustion reactions, because natural gas is relatively light and highly diffusive [37,38]. This fact predetermines the different velocities at which a laminar flame propagates into the flammable mixture normal to the reaction zone [37]. The laminar flame speed (LFS) is a fundamental flame property which depends on the type of fuel, the fuel-air ratio, and its initial thermodynamic conditions: pressure, temperature, and the equivalence ratio [38,39]; it increases with the reaction rate and the thermal conductivity of the unburnt gases since the preheated zone gets smaller [38]. A detailed assessment of the hierarchical nature of the key factors that affect combustion led to a posteriori assumption that LFS has the greatest influence on the duration of the combustion cycle when diesel fuel and D/CNG blends are juxtaposed under equal conditions. This difference will evidently be higher, the higher the share of CNG in the fuel blend. The same tendency for the paradigms applied in this study to compare and to contrast in-cylinder processes is illustrated in Figure 1.

3.2. Analytical Approach for Evaluation of Factors Influencing the Heat Release

As demonstrated in the previous chapter describing the differences between the single-zone and two-zone thermodynamic models of combustion, due to the implied averaging of temperature inside the cylinder during the combustion duration distinctive to a former method [40], the dual fuel combustion mode cannot be accurately predicted by the Woschni model in its original formulation. An over-expected amount of heat loss occurred in the internal combustion engine that, in our case, was predicted by the different laminar flame speed velocities for diesel fuel and natural gas. Therefore, it is necessary to modify the model and to verify that the proposed model remains valid under various operating conditions [41]. The single-zone model uses Vibe’s predictive staged correlation, which represents the gross energy release of fuel as a function of the crank angle [16]. The shape parameter (m) governs the PRR (pressure rise rate) and the location of the point of inflection of the burned mass fraction curve [16], while the angle of duration of the simple Vibe function (φ_z) characterizes the duration of heat release [16] at P_mi = 4.4, 6.2, and 8.2 bar. The values of m₀ and φ_z0 correspond to the shape parameter and the duration angle obtained at the indicated mean effective pressure P_mi = 8.2 bar. For a particular engine, the parameter φ_z0 is adjusted until a numerical model produces the pressure-φ curve to tune the profile of the experimental pressure-φ signal [16]. To determine the values of m and φ_z for the test fuels (D100, D60/CNG40, D40/CNG60, D20/CNG80) from the Woschni correlations (see Equations (17) and (18) [42,43], we used the experimentally established parameters φ_i, p_k, T_k, n, and α at P_mi = 4.4, 6.2, and 8.2 bar. The lower index ‘0’ stands for the same parameters φ_i, p_k, T_k, n, and α obtained at P_mi = 8.2 bar (base-pressure mode).

$$\frac{m}{m_0}={\left(\frac{{\phi }_{\tau i0}}{{\phi }_{\tau i}}\right)}^{0.5}·\left(\frac{p_k·T_{k0}}{p_{k0}·T_{k0}}\right)·{\left(\frac{n_0}{n}\right)}^{0.6}$$

(17)

$$\frac{{\phi }_z}{{\phi }_{z0}}={\left(\frac{{\alpha }_0}{\alpha }\right)}^{0.6}·{\left(\frac{n}{n_0}\right)}^{0.5}$$

(18)

where m—shape parameter, φ_z—duration of the heat release, φ_τi—induction period (crank angle degrees between the start of injection and the start of combustion), α—excess air coefficient, p_k—air pressure after compression, T_k—air temperature after compression, and n—crankshaft rotation speed.

The mathematical model of the IMPULS software works as follows: f and φ_z were set-up to one mode (mainly peak power), and for intermediate load modes, they were recalculated according to Woschni’s statistically determined dependencies—φ_i, p₀, T₀, α, n [27,44]. The good correlation between the calculated and experimentally established variables of the combustion cycle of a diesel engine running exceptionally on diesel fuel indicates that proper calibration of a model, which is an imperative condition for the accurate analysis of the ICE working process (see Figure 2,

Table 6. Compliance of the calculated vs. experimentally established variables of the combustion cycle of a diesel engine after the parameters m and φ_z being introduced (n = 2000 rpm) [27].

	Pme = 8.2 bar		Pme = 6.2 bar		Pme = 4.4 bar
Parameters	Experimental Values	Calculated Values	Experimental Values	Calculated Values	Experimental Values	Calculated Values
qcycl, g/cicl	0.0169	0.0169	0.0143	0.0143	0.01255	0.01255
pk, bar	1.35	1.35	1.25	1.25	1.175	1.175
Tk, K	333	333	318	318	315	315
pmax, bar	100.6	100.0	79.0	76	67.3	66.3
ηe	0.345	0.345	0.273	0.276	0.155	0.157
ηi	0.475	0.473	0.433	0.438	0.330	0.333
ηm	0.726	0.730	0.630	0.630	0.470	0.470
pme, bar	5.95	5.96	3.97	4.05	1.98	2.00
pmi, bar	8.2	8.17	6.30	6.40	4.20	4.28
α∑	2.54	2.43	2.95	2.79	3.19	3.08

).

In the next step, the precision of a single-zone model was further evaluated to assess the variables of the combustion cycle of a diesel engine when running on D20/CNG80 (p_mi = 4.4, 6.2., 8.2 bar). A positive correspondence of the results was obtained for the high load mode only (see

Table 7. Comparison of the calculated vs. experimentally established variables of the combustion cycle of a diesel engine running on D20/CNG80 (n = 2000 min⁻¹) [23].

Parameters	${\mathit{P}}_{\mathit{m}\mathit{i}}=4.4 \mathbf{b}\mathbf{a}\mathbf{r}$	${\mathit{P}}_{\mathit{m}\mathit{i}}=6.2 \mathbf{b}\mathbf{a}\mathbf{r}$	${\mathit{P}}_{\mathit{m}\mathit{i}}=8.2 \mathbf{b}\mathbf{a}\mathbf{r}$
${\phi }_{inj},°\mathrm{CA}$	−1	−1	−1	−13
$p_k,\mathrm{bar}$	1.25/1.28	2.38/1.36	1.47/1.45	1.35/1.34
$\alpha$	2.93/3.0	2.57/2.54	2.30/2.28	2.54/2.23
$T_k,\mathrm{K}$	319/325	329/336	343/342	333/335
$q_{cycl},\mathrm{g}/\mathrm{cycle}$	0.01457	0.01746	0.01975/	0.01695
$p_{max},\mathrm{bar}$	50.1/70.5	57.8/74.1	64.2/62.8	99.7/99.5
${\phi }_{comb},°\mathrm{CA}$	11/11	8.0/9.5	7/7	−5/−5
$P_{me},\mathrm{bar}$	1.988/4.1	3.976/5.18	5.964/6.00	5.964/6.08
$P_{mi},\mathrm{bar}$	4.194/6.25	6.242/7.32	8.204/8.14	8.204/8.25
${\eta }_m$	0.474/0.65	0.637/0.71	0.727/0.73	0.727/0.73
${\eta }_e$	0.131/0.27	0.218/0.293	0.288/0.310	0.366/0.35
${\eta }_i$	0.275/0.424	0.342/0.415	0.396/0.42	0.462/0.47

), where the difference between the calculated and experimental values did not exceed 2–6% [23]. When it comes to partial load modes (medium and low load modes), the disparity between the calculated and experimental values was reaching 50%. In particular, the most important discrepancy in the combustion cycle parameters (p_me, p_mi, η_e, η_i) was in range of 15–52% [23].

A thorough assessment of m and φ_z values distinctive to diesel fuel and D20/CNG80 was presented to understand the reasons for the discrepancy in calculated and experimentally established variables of the combustion cycle presented in Table 7.

The data presented in

Table 8. Comparison of the calculated m and φ_z values with experimentally established ones by using Equations (17) and (18) variables of the combustion cycle of a diesel engine after the parameters m and φ_z being introduced (φ_inj = −1 degCA BTDC).

	Pmi = 8.2 bar		Pmi = 6.2 bar		Pmi = 4.4 bar
	D100	D20/CNG80	D100	D20/CNG80	D100	D20/CNG80
m	0.58/0.58	0.6/0.6	0.55/0.57	0.57/0.6	0.52/0.54	0.52/0.54
φz	60/60	65/65	52/52	60/105	41/42	56/172

reveal three important aspects of the hypothesis testing. First, the investigation of trends and relationships using quantitative data (see Figure 1) confirms the direct dependence of α and η_i. Second, there is an evident difference in the combustion duration values for D20/CNG80 obtained for partial load modes if one compares the calculation output data (Equations (17) and (18)) with the experimental output data. Third, a trend that needs to be specified by further investigation was obtained: the switch from P_mi = 8.2 bar to P_mi = 4.4 bar resulted in the reduced combustion duration of a diesel fuel charge by 18 degCA BTDC (φ_z = 60 and φ_z = 42, respectively). Although the D20/CNG80 blend demonstrated a completely opposite trend: the load reduction from 8.2 bar to 4.4 bar resulted in a drastic increase in the φ_z value showing a huge difference of 112 degCA BTDC compared to the high load mode for P_mi = 8.2 bar (see Table 8).

To overcome the calculation deficiencies presented in Table 7 and Table 8, it was necessary to calibrate and configure a single-zone model separately for the D/CNG case. Experimentally derived data served as input parameters in the assessing the heat release characteristics X = f(φ) of an engine when D/CNG is used in dual fuel mode. The results are presented in Figure 3. The maximum inadequacy between the experimental values and the calculated values does not exceed 5–7% [27].

3.3. Exploration of the Indicated Parameters Associated with a Given Heat Release Pattern When D/CNG Is Used in Dual Fuel Mode

In our case, the relative ratios of calculated parameters φ_z/φ₀ and m/m₀ were further assessed by establishing the increase or decrease in their values compared to the values obtained for the ‘base-pressure mode’. This means that at P_mi = 8.2 bar, the relative ratio of φ_z/φ₀ becomes equal to 1. The proposed method allows the research results to be considered generalizable when the findings obtained for various P_mi values and D/CNG ratios formed a basis for the further rearrangement of the universally applicable Woschni equations to better characterize the dual fuel engine. The interpretation of the tendencies presented in Figure 4 and initial screening of the duration of heat release for diesel/CNG blends revealed that the excess air coefficient can be considered as the main factor influencing the angle of duration of the simple Vibe function. To better explain the relationships between the measured variables and latent variables (φ_z/φ_z0) in Equation (18), we employed the built-in function known as the power function, which is used to calculate the power of a given number (see Figure 4). It takes two arguments: the base and the exponent. The following trigonometric functions were obtained for various diesel/natural gas ratios:

D100: (φ_z/φ_z0) = 1.013 (α₀/α)^0.50

(19)

D60/NG40: (φ_z/φ_z0) = 0.995 (α₀/α)^−1.08

(20)

D40/NG60: (φ_z/φ_z0) = 0.993 (α₀/α)^−1.97

(21)

D20/NG80: (φ_z/φ_z0) = 1.013 (α₀/α)^−3.63

(22)

The denominators to the value of (α₀/α) in Equations (19)–(22) are comprised of numbers that do not differ from each other by more than ±1%. Thus, by assuming that it is equal to 1, this enables us to obtain the generalized equations for all investigated D/CNG ratios. To graph data on a line plot, we put engine’s operational characteristics at various advanced degrees of SOI angle relative to TDC (from −1 to −13 degCA) obtained for three different values of the indicated mean effective pressure. The results revealed that the excess air coefficient has a primary influence on the duration of the heat release since the points form clearly defined lines where the variables can be an example of a deterministic relationship, defined by the expression (φ_z/(φ_z0) = f(α₀/α). This coincides very well with the conclusion derived before that η_i = f(α).

The duration of heat release for D/CNG blends increases exponentially with the reduction in the mean effective pressure (at low and average loads) and the excess air coefficient, as demonstrated in Figure 4. This is in line with the findings published in Refs. [45,46]. This is mainly due to the increase in the flame propagation speed of the natural gas–air mixture [46]. The study explains that the combustion duration decreases significantly with increasing %CNG due mainly to the increase in the laminar flame propagation speed of the natural gas–air mixture. Concerning the duration of combustion, we observe that it is highest under dual fuel operation [47] at P_mi = 4.4 bar: as the amount of primary fuel increases from 40%CNG to 60%CNG and later to 80%CNG, the relative duration of heat release increases by a factor of ~2.5, 3.0, and 3.5 compared to diesel fuel, respectively. The further increase in the indicated mean effective pressure to 6.2 bar reduces this difference to ~1.5, 1.9, and 2.1 times, accordingly.

In Equation (17), the suitability of the original Woschni and Anisits function [48] to describe the shape of the instantaneous air pressure rise rate after the compression versus induction period during dual fuel engine operation is assessed, taking into account the changes in fuel composition (see Figure 5). The relative shape factor m is varied using the function of Woschni and Anisits and depends on variables on the right side of the equation (pressure and temperature when the inlet valve closes, as well as the rotation speed of the crankshaft [49]. To avoid the scattering between positive and negative values for the shape factor m, the left-hand side of Equation (17) is modified to (1 + m)/(1 + m₀).

A mathematical representation of the results, the actual heat release diagrams are replaced by simplified ‘Wiebe’ heat release diagrams, which have the same beginning and duration of combustion; the shape, however, is simplified and chosen so that if they are used for cycle simulations, the calculated values of peak pressure, power output, and fuel consumption are in agreement with the measured data. Such a simplified Wiebe heat release diagram is characterized by four parameters: the beginning and duration of combustion, the Wiebe parameter m, and the equivalence ratio.

The relative ratio of the heat release duration φ_z/φ₀ as a function of the relative ratio of the excess air coefficient increased to a certain degree and were obtained as a result of the Woschni correlation Equation (23) and derived from Figure 4 (Equation (24)), presented below:

$$\frac{{\phi }_z}{{\phi }_{z0}}={\left(\frac{{\alpha }_0}{\alpha }\right)}^{0.6}$$

(23)

$$\frac{{\phi }_z}{{\phi }_{z0}}={\left(\frac{{\alpha }_0}{\alpha }\right)}^{0.5}$$

(24)

Within the investigated range of α₀/α, the difference in calculation outcomes obtained by using the Equations (23) and (24) comprises only 1.5%. The right-hand of Equation (23) has been transformed from the original variables by adding the newly derived degree (β) to which the excess air coefficient will be increased in case of engine evaluation when D/CNG is used in dual fuel mode:

$$\frac{{\phi }_z}{{\phi }_{z0}}=1·{(\frac{{\alpha }_0}{\alpha })}^{\left(-4.83\beta +0.603\right)}$$

(25)

where: β is a constituting part of the CNG available in a blend (β = 0 stands for pure diesel fuel, β = 1 stands for pure CNG).

Henceforth, for the wide-scale adoption of the single-zone model to establish the heat balance characteristic X = f(φ) of an ICE working in a dual fuel mode, we highly recommend applying a newly derived Woschni correlation:

$$\frac{{\phi }_z}{{\phi }_{z0}}=1·{\left(\frac{{\alpha }_0}{\alpha }\right)}^{\left(-4.8\beta +0.6\right)}·{\left(\frac{n}{n_0}\right)}^{0.5}$$

(26)

It is worth mentioning that for the model-based assessment of dual-fuel engines, the second original equation by Woschni (see Equation (17)) will remain unmodified, since the modelling results exhibited a high-level correspondence with the experimentally derived values (see Table 8 and Figure 5).

4. Validation of the Newly Derrived Equation

The last step of this study was to examine the precision of the empirical formula presented in the form of an Equation (26), together with the often-used formula (Lukachev formula [50]) for computing the dependences of the laminar flame speed propagation on the pressure and initial temperature during the combustion of methane. A brief review of the Lukachev formula is described below.

In Ref. [50], the laminar flame speed (S_l) values were obtained experimentally [51,52,53,54,55,56] and calculated using the original formulas and using the detailed kinetic mechanism for methane/air combustion (GRI-Mech 3.0 detailed mechanism) [57]. Lukachev et al. (2016) [50] concluded that the results of the calculation of S_l according to the optimized mechanism designed to model natural gas combustion (GRI-Mech 3.0) agree well with the experimental data and can be used to refine the calculations of the dependence of the components S_l = f (φ, T_k, p_k) (see

Table 9. Refined calculations of the dependence of the components S_l = f (φ, T_K, P_K) [50].

Lukachev’s Formula [50]
$S_l=S_{l0}{\left(\frac{T_K}{T_0}\right)}^{\alpha }{\left(\frac{P_K}{P_0}\right)}^{\beta }$(27)
$S_{l0}=\left(145{\phi }_f^3-850{\phi }_f^2+1265{\phi }_f-325\right)$
$\alpha =a_1+a_2{\phi }_f+a_3{\phi }_f^2$	$\beta =b_1+b_2{\phi }_f+b_3{\phi }_f^2$
$T_0=800$	$P_0=\text{101,325}$
$a_1=7.6$	$b_1=-0.966$
$a_2=-10.8$	$b_2=1.1$
$a_3=5.1$	$b_3=-0.4715$

).

An accuracy of the newly derived Woschni correlation (Equation (26)) was compared with the results obtained using the original formula proposed by Lukachev et al. (2016) [50] (see

Table 10. Comparison of the calculated results by using an originally derived Equation (26) and reference Equation (27).

	Parameters	D20/NG80
Equation (26)	Pmi, bar	8.2	6.2	4.4
	Pk, Pa	6,520,000	5,860,000	5,120,000
	Tk, K	1100	1040	1030
	α	1.84	2.01	2.28
	α0/α	1.00	0.90	0.80
	φz/ φz0	1.00	1.55	2.20
Equation (27)	SL	45.640	31.808	21.815
	SL−1 (ϕz)	0.022	0.031	0.046
	φ = 1/α	0.543	0.498	0.439
	φz/φz0	1.00	1.43	2.09

).

Very good agreement was obtained between the relative duration of heat release values for both scenarios. Equation (26) demonstrated 5–7% lower φ_z/φ_z0 values for P_mi = 6.2 bar and P_mi = 4.4 bar compared to the results for the same loads obtained by Equation (27).

5. Conclusions

The research led to the successful derivation of the necessary correlations to model the heat release characteristics of an unmodified ICE operated in dual fuel mode when a single-zone thermodynamic model is used. The proposed correlation is simpler than models with the comparable uncertainty described in the literature (for example, in CONVERGE CFD software, the laminar flame speed can be calculated by one of three different approaches: (i) Metghalchi and Keck correlation (1982); (ii) Gulder correlation (1984); and (iii) user-supplied data tables [38]. The main findings can be summarized as follows:

(1)

It was found that the use of the single-zone model for the detailed assessment of an unmodified ICE in dual fuel mode gives a calculation error of 29% for the establishment of the total energy losses of heat energy generated during the combustion, compared to Q_{total_corr} (see Table 5). Furthermore, the difference between the calculated and experimentally established φ_z (duration of heat release) may exceed 307% (see Table 8), while the indicators of the diesel engine (p_me, η_e, etc.) may differ by more than two times (see Table 7).

(2)

The findings revealed that the linear increase in %CNG causes the corresponding reduction in the indicated thermal efficiency, where the variation of the excess air coefficient has the shortest range. From the process modelling point of view, these differences are attributed to the tangible change in combustion duration and shape factor m of the Vibe function, which are coupled to each other by scaling factors. The reduction in η_i/α ratio with an increase of %CNG is predetermined by the disproportionate increase in the mass fraction of natural gas entering the cylinder.

(3)

The results revealed that the excess air coefficient has a primary influence on the duration of the heat release since the points form clearly defined lines where the variables can be an example of a deterministic relationship defined by expression (φ_z/(φ_z0) = f(α₀/α).

(4)

Very good agreement was obtained between the relative duration of heat release values calculated using the newly derived Woschni correlation and the original formula proposed by Lukachev. The difference varied from 5.3% at P_mi = 6.2 bar to 8.3% at P_mi = 4.4 bar. Through the model-based assessment of dual-fuel engine, it is recommended to leave the second original equation by Woschni unmodified.

The proposed method allows one to predict the combustion and performance characteristics of a dual-fuel turbocharged engine with high accuracy. The minor disadvantage of the proposed correlation is the limitation of experiments to only one particular type of diesel engine.