Natural gas engines can be classified into two basic types based on the combustion method: single-fuel engines using spark plug ignition and dual-fuel engines that primarily use natural gas as the main fuel with high-reactivity fuels such as diesel as pilot fuel. The inclusion of pilot fuel in the latter addresses the challenge of ignition difficulty, especially under high fuel replacement rates, ensuring smoother and more reliable engine operation. In dual-fuel engines for vehicles, two main natural gas injection technologies are employed: intake manifold pre-mixed injection and in-cylinder direct injection. There are three popular types of dual-fuel engines for marine applications: two-stroke high-pressure injection, two-stroke low-pressure injection, and four-stroke intake manifold injection25. Figure1 illustrates the main structure of the marine dual fuel engine. The model comprises key components such as the intake and exhaust modules, fuel supply module, cylinder module, supercharger module, and air cooler module for dual-fuel engines. Table 1 presents the main technical specifications of the MAN 8L51/60DF dual-fuel engine26.
Structure and working principle of marine dual-fuel engine.
Basic equation for thermodynamics
The operation of dual fuel engine involves a series of complex chemical combustion, physical transmission, gas/liquid flow and heat transfer processes, which can be simplified by the following equations.
Energy conservation equation:
$$dU = dW + \sum d Q_{i} + \sum\limits_{j} {h_{j} } \cdot dm_{j}$$
(1)
where hj⋅dmj is the energy entering or exiting the system, kJ; Qi is the heat exchanged through the system, kJ; W is the work exerted on the piston by the expansion of the gas, kJ; hj is the enthalpy, kJ/kg; U is the internal energy of the system, kJ.
Conservation of mass equation:
$$\frac{dm}{{d\varphi }} = \frac{{dm_{s} }}{d\varphi } + \frac{{dm_{e} }}{d\varphi } + \frac{{dm_{B} }}{d\varphi }$$
(2)
where me is the outflow cylinder mass, kg; mB is the instantaneous fuel mass in the cylinder, kg; ms is the inflow cylinder mass, kg.
Ideal gas equation of state:
$$pV = mRT$$
(3)
where p is the pressure, MPa; V is the volume, T is the temperature, K; R is the gas constant, kJ/(kg.K).
Cylinder mathematical model
The cylinder module stands out as the most intricate and pivotal component within the entire engine modeling framework. The modeling of cylinder module comprises three aspects: the cylinder working volume, combustion, and cylinder perimeter wall heat transfer. This module plays a crucial role in capturing and simulating the complex dynamics occurring within the engine cylinder27. The cylinder working volume and volume change rate can be presented as follows:
$$Vs = \frac{{\pi D^{2} }}{4}\left\{ {\frac{S}{\varepsilon - 1} + \frac{S}{2}\left[ {\left( {1 + \frac{1}{\lambda }} \right) - \cos \left( {\frac{\pi }{180}\varphi } \right) - \frac{1}{\lambda }\sqrt {1 - \lambda^{2} \sin^{2} \left( {\frac{\pi }{180}\varphi } \right)} } \right]} \right\}$$
(4)
$$\frac{dVs}{{d\varphi }} = \frac{{\pi^{2} D^{2} S}}{8 \times 180}\left[ {\sin \left( {\frac{\pi }{180}\varphi } \right) + \frac{\lambda }{2} \cdot \frac{{\sin \left( {\frac{\pi }{180} \cdot 2\varphi } \right)}}{{\sqrt {1 - \lambda^{2} \sin^{2} \left( {\frac{\pi }{180}\varphi } \right)} }}} \right]$$
(5)
where Vs is cylinder working volume, L; λ is connecting rod to crank ratio; φ is crankshaft angle; S is fuel injection penetration; D is the diameter of the jet hole; ε is compression ratio.
In the context of non-predictive combustion models, the combustion rate is solely influenced by the crankshaft rotation angle, irrespective of whether the cylinder's operating conditions meet the combustion demand. Given the substantial impact of air–fuel ratio and fuel substitution rate on the combustion rate, coupled with the direct influence of injection parameter settings on the combustion process within the cylinder, the MAN 8L51/60DF dual-fuel engine employed in this study integrates the predictive combustion model DI-Pluse from Matlab/Simulink. This choice facilitates an exploration of the effects of these variables on power output, fuel consumption, and exhaust emissions in a dual-fuel engine. Distinct thermodynamic state parameters, compositions, and concentrations characterize the discrete thermodynamic partitions, and their independence from one another is a notable feature. The primary unburned zone, or first partition, comprises gases remaining in the cylinder after closure. The second and third partitions consist of the injected unburned zone-composed of inhaled fuel and gases-and the injected burned zone, representing combustion products. The mathematical model used to simulate this process is as follows.
Divided length:
$$L_{b} = 0.7D\left( {1 + 0.4\frac{r}{D}} \right)\left( {\frac{{\rho_{g} }}{{\rho_{l} v_{i}^{2} }}} \right)^{0.05} \left( \frac{L}{D} \right)^{0.13} \left( {\frac{{\rho_{l} }}{{\rho_{g} }}} \right)^{0.5}$$
(6)
where L is the length of the jet hole; vi is the jet velocity; ρl is the density of the injected fuel; and ρg is the density of the surrounding gas.
Angle of the spray cone:
$$\theta = 83.5\left( \frac{L}{D} \right)^{ - 0.22} \left( {\frac{D}{{D_{0} }}} \right)^{0.15} \left( {\frac{{\rho_{g} }}{{\rho_{1} }}} \right)^{0.26}$$
(7)
Fuel injection penetration:
$$S = 0.39\left( {\frac{2\Delta p}{{\rho_{l} }}} \right)^{0.5} t,\quad 0 < t < t_{b}$$
(8)
$$S = 2.95\left( {\frac{2\Delta p}{{\rho_{l} }}} \right)^{0.25} (Dt)^{0.5} ,\quad t > t_{b}$$
(9)
$$t_{b} = 28.65\frac{{\rho_{1} D}}{{\left( {\rho_{g} \Delta p} \right)^{0.5} }}$$
(10)
where tb is the oil beam splitting time; \(\Delta\)p is the jet pressure difference; D0 is initial diameter of spray hole.
Fire delay duration:
$$I = \int {\frac{1}{\tau }} dt$$
(11)
$$\tau = ap^{b} \varphi^{c} \exp (d/Z)$$
(12)
where a is the ignition delay factor; b is the pressure delay exponent; c is the empirical constant; d is the delayed activation temperature; Z is the ignition delay period; φ is the equivalence ratio.
Rates of combustion and volumetric absorption can be defined as:
$$\frac{{dm_{e} }}{dt} = \sum\limits_{j} {\left( {\rho_{f,j} ,\rho_{a} } \right)^{1/2} } 2\pi R_{j} \left( {x_{j + 1} - x_{j} } \right)C_{x} \left| {v_{j} - v_{a,j} } \right|$$
(13)
$$\frac{{dm_{B} }}{dt} = C_{1} \frac{{m_{c} }}{{\tau_{c} }}$$
(14)
$$\tau_{c} = \frac{{l_{T} }}{{S_{L} }}$$
(15)
where xj is the coordinate of the jth cell along the spray axis; Rj is the radius of the jth cell; vj is the velocity of the local jet; va, j is the velocity of the air; Cx is the empirical constant; SL is the laminar flame velocity.
Diffusion combustion models for computing volumetric absorption and burning rates:
$$\frac{{dm_{e} }}{dt} = C_{e} \frac{{m_{e} }}{{\tau_{e} }}$$
(16)
$$\tau_{e} = \left( {\frac{{l_{I} }}{\varepsilon }} \right)^{1/3}$$
(17)
$$l_{I} = C_{\mu }^{3/4} \frac{{k^{3/2} }}{\varepsilon }$$
(18)
$$\frac{{dm_{B} }}{dt} = \frac{{m_{e} }}{{\tau_{c} }}$$
(19)
where lI is the integrated length scale, Cu is 0.09, and Ce is the pre-combustion coefficient before the ignition point of each partition, k is insulation index.
Heat exchange occurs throughout the entire operation of the engine, such as between the external environment and the cylinder body, as well as between the cylinder body and the cylinder mixture. Heat losses are calculated through a heat transfer model, which transforms the complex heat transfer processes into a heat transfer coefficient model. The Woschni, Eichelberg, Annand, and Sitkei formulas can all be used to calculate the instantaneous heat transfer coefficient. In this paper, the Woschni heat transfer model is employed to calculate the heat transfer coefficient. This model is widely used in the field of internal combustion engine simulation modeling due to its simplicity and accuracy. The calculation formula is as follows28:
$$\alpha_{g} = 820 \cdot D^{ - 0.2} \cdot P_{z}^{0.8} \cdot T_{z}^{ - 0.53} \left[ {c_{1} \cdot c_{m} + c_{2} \cdot \frac{{V_{s} T_{1} }}{{P_{1} V_{1} }}\left( {P_{z} - P_{o} } \right)} \right]^{0.8}$$
(20)
where Pz is the instantaneous cylinder pressure, MPa; Tz is the instantaneous cylinder temperature, K; P1, T1, V1 is the pressure, temperature, and volume at the starting point of compression; cm is the average piston velocity, m/s; c1 is the gas velocity coefficient; c2 is combustion chamber shape coefficient.
Intake and exhaust mathematical model
Due to the influence of exhaust gas pressure fluctuations on the efficiency of the turbocharger and the residual gas quantity in the cylinder, and the direct impact of intake manifold pressure fluctuations on the intake volume within the cylinder, it is necessary to perform detailed calculations of the gas flow in the dual fuel engine’s intake and exhaust pipes.
In the simulation and modeling calculations, it is common to simplify the three-dimensional non-isotropic flow simulation of the dual fuel engine's intake and exhaust processes into an unsteady one-dimensional flow. During the engine's operation, the relevant parameters of the intake and exhaust systems are constantly changing. The pipe diameter is much smaller than the pipe length, resulting in radial flow dominating the entire flow process. Therefore, the intake and exhaust processes are treated as one-dimensional, and the simulation calculation method adopts the one-dimensional finite volume method.
Firstly, it is necessary to establish the differential equations relating the state parameters such as velocity, temperature, and pressure of the pipe segment to the corresponding crankshaft rotation angle and axial direction. The solution of the continuity equation, energy equation, and momentum equation are crucial for the mathematical model of the intake and exhaust systems.
The continuous equation is presented as:
$$\frac{dm}{{dt}} = \sum {m_{in} } = \sum \rho Au$$
(21)
where min is the inlet mass flow rate, kg/s; m is the mass of the work fluid in the boundary, kg; ρ is the density of the fluid in the pipe, kg/m3. A is the corresponding cross-sectional area of the pipe, m2; u is the velocity of the control body boundary, m/s.
The energy equation is expressed:
$$\frac{d(me)}{{dt}} = p\frac{dV}{{dt}} + \sum {\left( {m_{in} \cdot H} \right)} - hA_{s} \left( {T_{{\text{fluid }}} - T_{{\text{wall }}} } \right)$$
(22)
where As is the heat transfer surface area of the inner wall of the tube, m2; p is the pressure, MPa; H is the total enthalpy corresponding to the fluid, J/kg; h is the heat transfer coefficient, W/(m2.K); Tfluid is the fluid temperature, K; Twall is the wall temperature of the tube wall, K.
Turbocharger mathematical model
In a waste gas turbocharged engine, the waste gas turbine converts the energy from the engine's exhaust gases into mechanical work through a turbine. This mechanical work is then transferred to the compressor, which compresses the air to increase the volume of air entering the engine and raise the intake pressure. Throughout this process, the turbine and compressor must achieve power balance, speed balance, and flow balance29.
The power balance is expressed as:
$$W_{K} = W_{T} \cdot \eta_{TKm}$$
(23)
where WK is the work consumed by the compressor, kW; WT is the turbine work, kW; and ηTKm is the supercharger efficiency.
The rotational speeds of the compressor, turbine and supercharger are the same:
$$n_{K} = n_{T} = n_{TK}$$
(24)
where nK is the speed of the compressor; nT is the speed of the turbine; and nTK the shaft speed of the supercharger.
The exhaust gas flow rate is the sum of the air and fuel flow rate:
$$\dot{m}_{T} = \dot{m}_{K} + \dot{m}_{f}$$
(25)
where \(\dot{m}_{T}\) is the exhaust gas flow rate; \(\dot{m}_{K}\) is the air flow rate; and \(\dot{m}_{f}\) is the fuel mass flow rate.
The exhaust gas turbine plays a crucial role in powering the supercharger, and its model should encompass key features such as air flow, boost ratio, speed, and efficiency. The supercharger model, in turn, is characterized by two input parameters: air flow rate and supercharger speed, while the output parameters include temperature (To) and outlet pressure (Po).
The pressurization ratio and the turbine's rotor speed can be calculated using the principles of thermodynamics and Newton's laws. Specific flow rate and efficiency data for the compressor can be acquired by referring to relevant tables. Subsequently, the compressor's efficiency can be determined using the one-dimensional isentropic adiabatic flow theory. Finally, the outlet temperature (To) and outlet pressure (Po) can be calculated using the one-dimensional isentropic adiabatic compression process, with the arithmetic formula provided as follows.
$$P_{o} = P_{i} \pi_{c}$$
(26)
$$T_{o} = T_{e} \left\{ {1 + \frac{1}{{\eta_{c} }}\left[ {\pi_{c} \left( {\frac{k - 1}{k}} \right) - 1} \right]} \right\}$$
(27)
where Pi is the atmospheric pressure, Pa; Te is the environmental temperature, K.
The pressurization ratio and isentropic efficiency of the compressor are contingent upon both the air flow rate and the rotational speed of the compressor. These functional dependencies can be expressed through a mathematical relationship:
$$\pi_{c} = f\left( {n_{c} ,q_{c} } \right)$$
(28)
$$\eta_{c} = f\left( {n_{c} ,q_{c} } \right)$$
(29)
The above functional relationships are typically determined based on characteristic curves of a specific model of a turbocharged gas engine. Due to the absence of characteristic curves for the turbocharger, a modeling approach is adopted using the characteristic curves of a diesel engine. Specifically, the pressurization ratio and efficiency of the turbocharger are formulated as a one-dimensional table that varies with the load of the dual-fuel engine, represented using the one-dimensional lookup table module in Simulink. The modeling process of the turbocharger involves the utilization of additional structural or performance parameters, with specific numerical values as detailed in Table 2.
Intercooler model
The intercooler serves to cool the fresh air entering the cylinder, thereby improving the turbocharging efficiency of the engine. The input parameters for the intercooler model include the air temperature, air pressure, and air mass flow rate at the outlet of the turbocharger. As air flows through the intercooler, there is a pressure loss and temperature decrease. By introducing the intercooler effectiveness coefficient (i.e., the temperature at the outlet of the intercooler), the outlet temperature of the intercooler is calculated using a simplified method. By introducing the intercooler effectiveness coefficient (β) and combining it with experimental data, the outlet temperature of the intercooler is calculated using a simplified method, with β set to 0.88, resulting in the intercooler outlet temperature.
$$T_{j} = T_{o} (1 - \beta ) + \beta T_{w}$$
(30)
where Tw is cooling water inlet temperature, K; β is cooling coefficient is generally taken as 0.7–0.9.
Intercooler pressure loss Ps is defined as:
$$P_{s} = P_{s0} \left( {\frac{{q_{m} }}{{q_{m0} }}} \right)$$
(31)
where qm is actual flow rate of the intercooler, kg/s; Ps0 is pressure loss of the intercooler at design operating conditions, Pa; qm0 is intercooler design flow rate, kg/s.
Then the air pressure at the outlet is:
$$P_{r} = P_{o} - P_{s}$$
(32)
The values of the intercooler parameters are shown in Table 3.
Cylinder intake model
The cylinder charging coefficient is a parameter that quantifies the amount of fresh air remaining in the cylinder at the end of the intake stroke relative to the volume of the intake manifold. Research results indicate that there is a functional relationship between the cylinder charging coefficient and the engine speed. Therefore, this article uses the least squares method to fit the engine speed and boost coefficient data to obtain the boost coefficient curve:
$$\eta_{v} = a_{0} + a_{1} \times n + a_{2} \times n^{2}$$
(33)
where a0, a1, a2 are the fitted constants.
Since the model established is an average model for studying air–fuel ratio control methods, it overlooks the residual exhaust gases in the cylinder after each engine cycle. According to the ideal gas equation, the intake air flow can be calculated as:
$$q_{m1} = \frac{{\eta_{v} nVsP_{r} }}{{120RT_{j} }}$$
(34)
Fuel sub-model
For the MAN 8L51/60DF LNG/diesel dual-fuel engine, a fuel sub-model is established based on the propulsion characteristics. Under propulsion characteristics, the relationship between power and speed is given by:
$$P = F \cdot n^{3}$$
(35)
where P is engine power, kW; F is scale factor.
The corresponding values of propulsion characteristic speed and load are shown in Table 4.
To accurately establish the fuel supply sub-model for marine liquefied natural gas/diesel dual-fuel engines, it is necessary to explore the fuel injection quantity under propulsion characteristics30. In dual-fuel engines, both diesel and natural gas are introduced into the engine cylinder for combustion. The calculation necessitates the incorporation of the substitution ratio, which denotes the reduction in the amount of diesel consumed by the engine in dual-fuel operating mode compared to the consumption in pure diesel operating mode31. This ratio is expressed as a percentage of the diesel consumption in pure diesel operating mode32.
$$T_{m} = \frac{{B - B_{o} }}{B} \times 100\%$$
(36)
where Tm is diesel replacement rate, %; B is diesel consumption in diesel-only mode, kg/h; Bo is diesel consumption in dual-fuel mode, kg/h.
The natural gas entering the cylinder precisely substitutes the reduced diesel fuel in the combustion process, and the substitution rate equals the doping rate. The calculation of the substitution rate is as follows33. Under specific operating conditions and at a specific substitution rate, the cylinder injection volume per cycle and the natural gas injection volume are determined based on theoretical considerations.
$$\left\{ {\begin{array}{*{20}l} {H_{ul} \cdot B_{l} = \left( {B - B_{o} } \right) \cdot H_{uo} } \hfill \\ {\frac{{B - B_{o} }}{B} = T_{m} } \hfill \\ \end{array} } \right.$$
(37)
where Bl is natural gas consumption in dual fuel mode, g/cycle; Hul is low calorific value of natural gas, MJ/kg; Huo is low calorific value of diesel, MJ/kg.
Based on the propulsion characteristics, the circulating injection volume and natural gas injection volume are calculated under various working conditions and different fuel substitution rates, following the aforementioned calculation process. Some of the calculation results are presented in Table 5.
Emissions calculation
The calculation of NOX and HC emissions is as follows:
$${\text{MNO}}_{{\text{x}}} = {\text{NO}}_{{\text{x}}} {\text{W}} \times \left[ {\begin{array}{*{20}l} {0.001586 \times (1 - Tm) + \left( {\frac{{{\text{CH}}_{4} }}{100} \times 0.001621 + \frac{{{\text{C}}_{3} {\text{H}}_{8} }}{100} \times 0.001603} \right.} \hfill \\ {\left. { + \frac{{{\text{C}}_{4} {\text{H}}_{10} }}{100} \times 0.001600} \right) \times Tm} \hfill \\ \end{array} } \right] \times \frac{{\dot{m}_{T} }}{1000}$$
(38)
$${\text{MHC}} = {\text{HCW}} \times \left[ {\begin{array}{*{20}l} {0.000479 \times (1 - Tm) + \left( {\frac{{{\text{CH}}_{4} }}{100} \times 0.000558 + \frac{{{\text{C}}_{3} {\text{H}}_{8} }}{100} \times 0.000512} \right.} \hfill \\ {\left. { + \frac{{{\text{C}}_{4} {\text{H}}_{10} }}{100} \times 0.000505} \right) \times Tm} \hfill \\ \end{array} } \right] \times \frac{{\dot{m}_{T} }}{1000}$$
(39)
where NOXW represents the wet volume concentration of nitrogen oxides in parts per million (ppm), HCW denotes the wet volume concentration of hydrocarbons in parts per million (ppm), MNOX signifies the mass flow rate of nitrogen oxides in grams per kilowatt-hour (g/kWh), and MHC denotes the mass flow rate of hydrocarbons in grams per kilowatt-hour (g/kWh). The detailed calculation and derivation of these parameters can be found in literature34. Table 6 provides information on the main components of the fuel along with their respective proportions.
