Parametric Knocking Performance Investigation of Spark Ignition Natural Gas Engines and Dual Fuel Engines

: Both spark ignition (SI) natural gas engines and compression ignition (CI) dual fuel (DF) engines suffer from knocking when the unburnt mixture ignites spontaneously prior to the flame front arrival. In this study, a parametric investigation is performed on the knocking performance of these two engine types by using the GT-Power software. An SI natural gas engine and a DF engine are modelled by employing a two-zone zero-dimensional combustion model, which uses Wiebe function to determine the combustion rate and provides adequate prediction of the unburnt zone temperature, which is crucial for the knocking prediction. The developed models are validated against experimentally measured parameters and are subsequently used for performing parametric investigations. The derived results are analysed to quantify the effect of the compression ratio, air-fuel equivalence ratio and ignition timing on both engines as well as the effect of pilot fuel energy proportion on the DF engine. The results demonstrate that the compression ratio of the investigated SI and DF engines must be limited to 11 and 16.5, respectively, for avoiding knocking occurrence. The ignition timing for the SI and the DF engines must be controlled after -38°CA and 3°CA , respectively. A higher pilot fuel energy proportion between 5% and 15% results in increasing the knocking tendency and intensity for the DF Engine at high loads. This study results in better insights on the impacts of the investigated engine design and operating settings for natural gas (NG)-fuelled engines, thus it can provide useful support for obtaining the optimal settings targeting a desired combustion behaviour and engine performance while attenuating the knocking tendency.


Introduction
Natural gas (NG) has gradually been employed to internal combustion engines (ICE) as a substitution of the traditional liquid fuels because of its positive effect on reducing gas emissions [1]. Due to its higher octane number (120 in comparison with 87 for gasoline), natural gas exhibits a better antiknock performance, which allows the NG-fuelled engines to employ a greater compression ratio for improving their thermal efficiency [2,3]. The NG-fuelled engines are classified into the following three basic types: (a) spark ignition (SI) NG engines; (b) SI bi-fuel engines running on either gasoline or NG fuel (usually converted SI gasoline engines), and; (c) compression ignition (CI) dual fuel (DF) engines running on diesel and NG fuels (usually converted diesel engines) [3]. The latter can be further classified to premixed combustion DF engines and high-pressure direct gas injection DF engines (or gas diesel engines) according to the NG admission methods [4]. As the gasoline and NG fuels cannot be mixed and burnt simultaneously in the gasoline/NG bi-fuel engines [5], the thermal efficiency of these engines is constrained by the antiknock requirements of the gasoline fuel. This means that the NG advantages cannot be fully exploited in the bi-fuel engines. On the contrary, the SI natural gas engines [6] and diesel/NG DF engines [7] are designed with higher compression ratios to take advantage of the NG antiknock characteristics.
For the majority of the SI natural gas engines, the NG is injected in the intake manifold, where it mixes with the inflowing air forming a mixture, which enters into the engine cylinders. At the end of the compression process, the NG-air mixture is ignited by a spark plug, thus generating a flame front, which travels within the combustion chamber with an almost constant velocity [8]. The compression ratio of the SI natural gas engines is typically kept in the range from 8:1 to 13:1 due to the higher NG octane number [9]. Although a higher power output is expected for the SI NG engines (in comparison to the gasoline engines with similar dimensions), there are two main factors limiting these engines power. The first is the reduced engine volumetric efficiency due to the lower density of the injected NG at the inlet conditions (the natural gas occupies a part of the air−gas fuel mixture trapped within the engine cylinders) [10]. The second is the advanced spark timing, which can be up to 10°CA earlier compared to gasoline engines, due to the slower flame front propagation speed. This results in increased pressure at the latest part of the compression phase, thus reducing the engine power output [11]. In order to compensate the power loss caused by the reduced volumetric efficiency, forced induction by employing supercharging or turbocharging is usually applied to increase the trapped air-NG mixture amount [12,13]. Alternatively, NG direct injection into the engine combustion chamber can be also employed [14,15], however this is not a widely applied method due to the characteristics of the high-pressure gaseous fuel injection system. A measure that assists in increasing the flame propagation speed is the usage of a small amount of another fuel with higher flame speed (for example hydrogen) [16]. In this case, a retarded spark timing can be employed resulting in the engine power output increase.
For CI dual fuel engines, a small amount of diesel fuel is usually injected into the combustion chamber to initiate the NG ignition as the NG has a lower cetane number and can hardly be ignited by compression [15,17]. This pilot injection provides a high energy source to ensure reliable and powerful ignition of the mixture, which is needed when running with a high specific cylinder output and lean NG-air mixture [18]. In most premixed combustion DF engines, the NG-air mixture formation takes place outside the cylinder by inducting the NG into the manifold prior to the inlet valves closure [19]. This reduces the engine volumetric efficiency by 1-4% compared to the normal CI engines [17]. Direct gas injection technologies, including low-pressure gas injection and high-pressure gas injection, are applied to compensate the reduced volumetric efficiency by injecting the gaseous fuel directly into the cylinders. The low-pressure gas injection takes place during the compression phase, whilst the high-pressure gas injection initiates close to the Top Dead Centre (TDC) [4]. For the premixed combustion DF engines, the NG-air mixture is exposed to the compression process and hence to increasing pressure and temperature levels. This may lead to severe knocking and thus constrains the compression ratios of the premixed combustion DF engines, resulting in lower thermal efficiency compared to their diesel alternatives [20]. High-pressure direct injection makes the gaseous fuel follow a diesel combustion process [4], which renders the engine practically insensitive to the fuel knocking properties and achieves similar thermal efficiency with that of diesel engines. However, due to the complexity and associated cost of the fuel feeding system, direct gas injection engines are typically used for high power output ranges. The investigated DF engine herein refers to a premixed combustion DF engine with inlet port/manifold NG injection.
Although the SI natural gas engines and DF engines design and operating characteristics including their ignition method and volumetric efficiency are different, both engines suffer from combustion instabilities, namely, the misfiring and knocking phenomena. Misfiring is an abnormal engine operating condition with incomplete or no combustion for one or more engine cylinders and/or cycles. As the NG-air mixture exhibits a low flame front velocity [21], the flame in NG-fuelled engines propagates slowly, which causes combustion instability with high cyclic variation [22]. When the flame front propagation ceases, misfiring and incomplete combustion may occur in the end-gas zone, which results in low efficiency and high unburnt hydrocarbon (HC) emissions [23]. Knocking is associated with the auto-ignition of the unburnt fuel-air mixture ahead of the flame front [24]. It is a common instability in the premixed combustion engines [25], which may cause severe damage to the engine structure [26]. For avoiding these combustion instabilities, the NG-fuelled engines must operate at a relatively narrower air-fuel ratio range compared to diesel engines with similar dimensions. It is important to achieve the desired combustion behaviour whilst limiting the knocking occurrence at a minimum level. Since the knocking probability increases with higher in-cylinder temperature and pressure, the engine compression ratio must be kept below a maximum threshold for ensuring the engine knock-free operation. This, however, limits the engine efficiency.
For detecting the knocking phenomenon in internal combustion engines, the following two broad categories of experimental methods are used: (a) direct measurements by using an Intensified Charge Coupled Detector (ICCD) camera and Laser-induced Fluorescence (LIF) imaging [27,28]; (b) indirect measurements of in-cylinder pressure, cylinder block vibrations/noise and/or exhaust gas temperature [29,30]. Experimental studies require considerable resources and are costly, therefore they are not extensively used for the evaluation of the engine settings effect on knocking early in the engine design phase. In this respect, a more cost-effective method to predict the knocking onset is by using simulation methods [2], which can be generally classified into three categories: (1) detailed chemical kinetic mechanisms for the pre-flame (or low-temperature) reactions [31,32]; (2) simplified chemical kinetic mechanisms for the pre-flame reactions [33][34][35]; (3) empirical formulas based on Arrhenius expression [36].
Detailed chemical kinetic mechanisms [31,32] usually consider hundreds of species and thousands of elemental reactions, which provides the most detailed prediction of the local thermodynamic parameters including temperature, pressure and species concentration as well as the knocking onset. Nevertheless, this comes at the expense of huge computational cost, thus rendering this method not applicable in multi-variables and multi-objective optimisation studies based on engine cycle simulations [2]. Simplified chemical kinetic mechanisms were proposed to reduce the reactions number, which reduced the computational cost while maintaining acceptable prediction accuracy. Depending on the specific fuel type and application, the reactions number of these simplified mechanisms vary from several [33] to hundreds [34,35]. Computing power and computational cost are the major limitations of using detailed or simplified mechanisms for the prediction of autoignition of unburnt gases, especially when quasi-real time application is necessary. To overcome these limitations, a one-step reaction formula was proposed by Livengood and Wu [36] by integrating the reciprocals of the unburnt zone autoignition delay that was estimated by using an Arrhenius expression. Due to its simplicity and low computing cost, the Livengood-Wu integral has been widely used for the knocking onset prediction in engine cycle simulations. A number of studies investigated the refinement of the Arrhenius expression for more accurate estimation of the autoignition delay in the unburnt zone. Douaud-Eyzat [37], Franzke [38] and Worret [30] updated the Arrhenius expression by introducing the octane number (ON) and calculated the employed parameters values. Hann [39] modified the Arrhenius equation for predicting the knocking onset of a SI natural gas engine operating with fuel composition variations. Hoepke [40] developed an empirical formula by incorporating the EGR effects for SI engines operating at a stoichiometric air-fuel ratio. Chen [41] refined Hoepke's model by considering the effect of excess air ratio, which predicts the knock onset with high accuracy under the fuel enrichment conditions both with and without EGR. Bares [42] developed a new knock model, which combines the Douaud-Eyzat model and an exogenous noise disturbing the in-cylinder temperature, for estimating the knocking probability. Song [43] introduced an Arrhenius expression of a cool-flame eliminated ignition delay to improve the predictability of Livengood-Wu knock model [37] for spark-ignition engines operating with two-stage ignition fuels. Among all the above studies, the Worret knocking criterion [30] is one of the most extensively tested empirical formulae for both SI engines and dual-fuel engines.
As for the engine modelling and simulation, simpler zero-dimensional (0-D) up to more complicated three-dimensional models (3-D) can be employed for the knocking prediction [2]. The 0-D models employ the assumption of uniform variations of the working media state and concentration within the engine components, which are quite effective for engine performance prediction [44]. For predicting knocking [45,46], the in-cylinder needs to be divided into at least two zones (the burnt zone and the unburnt zone or end-gas zone), in order to characterize the flame propagation and end-gas auto-ignition. Thus, the 0-D combustion models for knocking prediction are commonly associated with the cylinder volume division similarly to the two-zone models [47,48] and multi-zone models [49][50][51]. Xiang et al. [47] proposed a two-zone model based on the Douaud-Eyzat correlation [37] to predict the knocking performance of the SI natural gas engines by assuming that the burnt and unburnt zones have a cylindrical shape. Javaheri et al. [49] studied the influence of the natural gas composition on the knocking combustion in SI gas engines with a three-zone approach, which characterised the reactions in the burnt zone by chemical equilibrium calculations. The 3-D Computational Fluid Dynamics (CFD) models are capable of providing the most detailed analysis of the engine in-cylinder combustion processes when knocking occurs, adequately predicting the pre-flame reaction [52], the heat release rate [53] and emissions formation [54]; therefore, they are appropriate for obtaining better insight of the involved thermo-physical-chemical processes. However, 3-D models have disadvantages of high complexity and massive simulation time [55], which may limit their application in knocking prediction. In order to improve the knocking prediction accuracy, 1-D models for modelling the manifolds are occasionally used coupled with knocking prediction models, as the knocking tendency is significantly affected by the air-fuel ratio obtained during the gas-exchange process [56]. Noda et al. [57] proposed a transient knocking prediction technique for SI engines by coupling a 0-D knocking model with chemical kinetics and a 1-D gas exchange model.
As concluded from the discussed literature [44][45][46][47][48], a two-zone 0-D model could be an effective tool for the knocking behaviour analysis in NG engines as it is capable of characterising the end-gas zone temperature with the simplest combustion zone division. In addition, a 1-D model [56,57] could be employed for simulating the gas exchange process in order to improve the simulation accuracy of the entire engine cycle. Thus, a two-zone 0-D combustion model coupled with a 1-D gas exchange model could effectively be employed for predicting both the engine performance as well as knocking compromising between the model complexity and the required computational time. Another conclusion drawn from the previous studies on knocking [47][48][49][50][51][52][53][54][55][56][57] is that few studies focused on comparing and comparatively analysing the SI natural gas engines and the diesel/NG DF engines, which may exhibit different knocking behaviour due to their distinct ignition methods. In particular, the flame propagation in SI natural gas engines initiates from a single ignition source [10], whilst that in CI dual fuel engines starts from multiple ignition points by employing a pilot fuel injection [18,19]. Hence, the knocking occurrence, which is highly dependent on the flame front propagation, could be distinct in SI natural gas engines and DF engines.
This study aims at the parametric investigation and comparative analysis of the knocking behaviour of these two types of NG-fuelled engines. Appropriate simulation models of the 1-D/0-D type are set up in the GT-Power software [58] for the investigated SI natural gas and diesel/NG DF engines. The in-cylinder processes were modelled by employing a two-zone 0-D model for the knocking prediction, whereas the gas exchange processes in the engine manifolds were modelled by employing a 1-D approach. The combustion is modelled by using a Wiebe functions model, which provides adequate accuracy for characterising the heat release rate (HRR) of SI and CI engines. The knocking onset is characterised by using the Worret knocking criterion [30], which is one of the most extensively tested correlations estimating the autoignition delay in the unburnt zone. The developed models are first validated against experimental measurements and subsequently are used for the parametric investigation of the effects of the engine settings on the engine performance and knocking behaviour. By considering the engine power output and the knocking performance, the derived results are discussed providing recommendations for the optimal settings of the investigated engines, which could form the basis for more systematic optimal engine design studies in the future.
The novelty of this study is summarised as follows: (a) an experiment-based interpolation model of the double Wiebe type is developed for characterising the combustion of a diesel/NG dual fuel engine; (b) a comprehensive parametric investigation for quantifying the effects of the engine design and settings on the engine power and knocking performance is conducted; (c) a comparative analysis of the knocking behaviour of a SI natural gas engine and a diesel/NG dual fuel engine is carried out.

Model Description
This section describes the mathematical principles of the developed models employed for the two investigated engines (SI natural gas engine and DF engine) performance predictions and knocking behaviour study. The GT-Power computational environment was used in the present study as it is widely acknowledged in both academia and industry [59]. A number of sub-models are used to simulate the engine processes including fuel injection, combustion, heat transfer and friction as well as the knocking behaviour. These are described in the following sections.
The following assumptions are considered for developing and setting up the two models: (1) The cylinder is modelled by employing two zones-the burnt zone and the unburnt zone.
(2) The pressure at any time is considered to be uniform within the cylinder.
(3) Each cylinder zone has its own uniform temperature. There is no heat transferred between the two zones.
(4) Gas in each zone is assumed to be homogenously distributed and the basic species are considered as ideal but non-perfect gases.
(5) Blowby and valves leakage in the engine cylinders are not considered.

Fuel Injector Model
The fuel injector models of SI natural gas engines and DF engines need to be selected according to the fuel type, the engine injection type and the fuel injection control method.

NG Injector Model for the Investigated SI NG Engine
For the investigated SI natural gas engine, the natural gas is injected into the manifold and sufficiently mixed with the inlet air prior to entering into the cylinders. The air-fuel equivalence ratio is primarily controlled by governing the throttle opening percentage and the gas injection mass, which indicates that the simulation accuracy does not benefit from the detailed modelling of the injection pulse. Thus, an imposed air-fuel equivalence ratio and a known injection rate are used to determine the injection pulse width in GT-Power.

NG and Pilot Diesel Injector Models for the Investigated DF Engine
The investigated DF engine is equipped with two fuel supply systems, one for the natural gas and another for the diesel pilot fuel. The operating conditions of the natural gas injector are characterised by the instantaneous injection rate. A constant injection rate was imposed. The pilot diesel injection is modelled by employing the injection timing and injection mass flow rate. This approach is commonly used for single pulse injection providing as input a periodic pressure or mass flow rate profile [58].

Combustion Model
Semi-empirical formulas are usually employed to simulate the combustion heat release rate (HRR). Typically, the HRR-determining methods [59] include the Triangular Exothermic function, the Polygon-hyperbola function and the Wiebe function, among which the Wiebe function is most widely used [60]. A single Wiebe function [61,62] is represented by the following equation: where, x is the burnt fuel fraction; φ is the crank angle (°CA); a is the coefficient related to the combustion efficiency, which is usually set at 6.9078 to maintain a combustion efficiency of 99.9%; τ is the normalised combustion time; m is the Wiebe exponent that determines the shape of the HRR curve.
The normalised combustion time (τ) is determined according to the following equation: where, SOC is the crank angle at the start of combustion (°CA); Δφ is the combustion duration (°CA). Prior to developing the combustion model, the number of the employed Wiebe functions need to be decided according to the fuel type and the engine combustion characteristics. The SI natural gas engines exhibit the main characteristics of premixed combustion, which can be adequately represented by employing a single Wiebe function. The Wiebe parameters (SOC, Δφ and m) of SI engines at varying operating points can be calculated according to following equations [60,63].
where, IGD is the ignition delay, which is used to calculate the SOC; subscript A represents the nominal operating point; fIT, h , hn and hwi are functions of ignition timing, residual gas portion, rotational speed and indicated work, which were reported in Witt [63]. DF engines, which employ a small amount of pilot fuel to ignite the gaseous fuel, exhibit characteristics of both the diffusion combustion and the premixed combustion, and therefore require more complex models rather than a single Wiebe function. As reported in pertinent literature [61,64], combustion models based on double-Wiebe functions exhibited sufficient accuracy for predicting the burned mass fraction, heat release rate and in-cylinder pressure of DF engines. The model can be developed into a widely applicable model if enough experimental data or CFD simulations are available to find the unique set (or sets) of parameters that minimise the error at most conditions. Hence, a double Wiebe function as calculated by Equation (6) is used in this study to represent the combustion of the investigated DF engine. Specifically, the first Wiebe function is used to represent the pilot diesel fuel combustion whereas the second is used to represent the NG combustion.
where, b1 represents the weight factor of the pilot diesel fuel combustion, 0 < b1 < 1; τ1 and τ2 are the normalised combustion time of pilot diesel and NG, respectively; m1 and m2 are the Wiebe exponents of pilot diesel and NG, respectively. Up to now, there is no published study which has presented equations calculating the double-Wiebe parameters for modelling the combustion of dual-fuel engines at varying operating points. The combustion model parameters (SOC, Δφ and m) of double-Wiebe function, which determine the combustion profile for both the diesel and the gas modes, were calibrated at different loads (at steady-state conditions), so that the predicted engine cylinder parameters (maximum pressure, compression pressure, IMEP) sufficiently matched their respective experimental values.
The calibrated values of the combustion model parameters are stored in a database as functions of the engine load and rotational speed. Quadratic interpolation is used for calculating the values of the respective combustion model parameters for each operating condition.
In the employed two-zone combustion model, the cylinder is divided into two zones at the start of combustion-in particular the unburnt zone and the burnt zone. The amount of fuel-gas mixture that is transferred from the unburnt zone to burnt zone during the combustion process is defined by the burnt rate, which is calculated by the employed Wiebe functions. Each zone temperature is calculated by employing the following energy conservation equations [47]. Unburnt zone: where, the subscript '1' denotes the unburnt zone; m1 is unburnt zone mass; u1 is unburnt zone energy;  ent g m is the mass flow from unburnt zone to burnt zone; ent g h is enthalpy of the gas that leaves unburnt zone;  exit sg m is the mass flow from burnt zone to unburnt zone; exit sg h is the enthalpy of the gas that leaves burnt zone; p is the cylinder pressure; V1 is the unburnt zone volume; Q1 is the unburnt zone heat transfer rate. Burnt zone: where, the subscript '2' denotes the burnt zone; ξ is the combustion rate; LHV is the lower heating value of the employed fuel.

Heat Transfer Model
The widely acknowledged Woschni heat transfer model is employed for the calculation of the instantaneous gas to wall heat transfer coefficient [58] according to the following equation: where, p is the in-cylinder pressure (bar); T is the in-cylinder temperature (K); D is the cylinder bore (m); pa, Ta and Va are the pressure (bar), temperature (K) and volume (m 3 ) at the inlet valve closure, respectively; Vs is the cylinder displacement (m 3 ); p0 is the in-cylinder pressure at motoring conditions (bar); cm is the average piston speed (m/s); C1 is the coefficient depending on the airflow velocity; C2 is the coefficient related to combustion chamber shape.

Knocking Prediction Model
There are three most commonly used knocking prediction means [2]: detailed chemical kinetic mechanisms, simplified chemical kinetic mechanisms and empirical formulae based on Arrhenius expression. The Worret knocking criterion [30] is one of the most extensively tested empirical formulae for knocking prediction in both SI engines and dual-fuel engines. In the Worret model, the knocking probability is estimated by means of an induction time integral (Ik) and a constant parameter K, which are calculated according to the following equations, as reported in [60]: where n is the engine speed (rpm); θCS signifies the crank angle of calculation start, whilst θK is the crank angle of knocking occurrence; a, b, c, aIK, aK, bK and cK are coefficients provided by [60]; θ50mfb and θ75mfb are the crank angle of 50% and 75% mass fraction burnt points, respectively; subscript ref represents values obtained at the reference operating condition. To obtain Ik and K, a set of reference parameters including the crank angles at 75% and 50% of burnt fuel mass fraction, the air-fuel ratio and a reference K-value that are calculated for a known cycle in the knocking boundary. Then the knocking sub-model calculates at each operating point, a characteristic crank angle (θE), based on the calculated K-value, as well as the crank angle (θk) at which Ik becomes equal to one. θE is calculated according to the following equation where θSOC is the crank angle where combustion starts; ΔθCD is the combustion duration.
Knocking is predicted to occur at θk with a probability when the following inequality holds: θk < θE. This probability derives from the model uncertainty that Ik and K have a dispersion range of 15% and 5%, respectively. In this respect, the knocking probability is calculated according to the following equation: In (13) where ΔθKDR is the crank angle resulting from the superimposition of Ik and K dispersion ranges; akp and bkp are coefficients provided by [60]. The Unburnt zone Mass Fraction (UMF), which refers to the mass fraction of the unburnt fuel at knocking onset, has been employed as a knock index to describe the knocking intensity in engine cycle simulations [65,66]. The UMF represents the energy released within a short time when the knocking occurs, thus it can be used to characterise the knocking intensity. Larger UMF values correspond to a stronger knocking intensity.

Model Layout
In this study, a 2135 SI natural gas engine and a YC6K dual fuel engine were used for model development representing the mentioned two types of NG-fuelled engines. To avoid mentioning them each time with their full name for a better reading, these two engines were listed as Engine A and Engine B, respectively. Figure 1 shows the GT-Power graphical representation of the developed simulation model for Engine A. The model includes sub-models for the air filter, inlet manifold, the natural gas injectors, the cylinders, the exhaust manifold and the engine crankshaft. As Engine A is naturally aspirated, the inlet manifold is connected to the ambient environment directly.

Engine B
The GT Power model for the Engine B is presented in Figure 2. Although the Engine B is turbocharged, the turbocharger unit model is not included in the presented engine modelling, as the detailed modelling of the turbocharger does not benefit the investigation in this study and increases the required computational time. The function of the absent turbocharger sub-model is replaced by using the measured thermodynamic parameters (temperature and pressure) downstream of the compressor and upstream of the turbine as boundary conditions.

Engine Setup and Model Validation
In order to validate the simulation models for Engine A and Engine B, engine tests were carried out to measure the engine operating parameters, including the in-cylinder pressure, fuel consumption rate, air mass flow, etc. The measured in-cylinder pressures were used for the validation of the two developed models, respectively. However, the knocking prediction models were not validated by using experimental data as the engine operation with knocking conditions was not tested to avoid jeopardising the engine integrity. To calibrate the knocking parameters, the knocking prediction models of Engine A and Engine B were run virtually at respective knocking boundaries. Table 1 shows the main characteristics of the investigated engines. Engine A is converted from a diesel engine by removing the diesel injection and installing a spark plug, whilst decreasing the compression ratio from 16 to 11. This engine is naturally aspirated and runs at a low air-fuel equivalence ratio of around 1.2. Engine B is turbocharged operating at a higher air-fuel equivalence ratio value (2.1 at nominal load). In this study, crank angles before top dead centre (BTDC) are defined as negative, whilst those after top dead centre (ATDC) are defined as positive. The ignition timing of Engine A and pilot injection timing of Engine B are 36 °CA BTDC and 5 °CA BTDC, respectively.

Experimental Set up
The experimentally measured parameters for the two investigated engines include in-cylinder pressure, fuel mass flow rate, inlet air mass flow rate, inlet/exhaust manifold pressure and emissions. Heat release rate (HRR) is calculated based on the measured in-cylinder pressure. The technical specifications of the employed sensors are shown in Table 2. Three flow meters are used to measure the mass flow rate of the pilot diesel fuel (AVL 735C), the natural gas (E+H 83F25-XRW2/0) and the inlet air (ABB FMT 700). A piezoelectric in-cylinder pressure sensor (AVL GU22CK) is installed on the cylinder head to obtain the in-cylinder pressure diagram, which is then transferred to the combustion analyser (KiBox To Go 2893) for heat release analysis. The emission analyser (AVL AMA i60 R1) receives exhaust gas samples from the gas sampling device and provides the concentrations of various exhaust gas species including NOx, CO and HC. A low-pressure indicating sensor (LP11DA) is used to measure the exhaust pressure at the exhaust gas manifold (upstream turbine for the Engine B). A natural gas chromatographic analyser [74] is used to measure the composition of the employed compressed natural gas (CNG), as shown in Table 3. The Lower Heating Value (LHV) is 47.18 MJ/kg. The diesel fuel used in the DF engine testing was the light fuel of grades 0# [75]. The diesel LHV is 42.652 MJ/kg whilst its Cetane Number is 50.

Model Validation
The in-cylinder pressure signals, as measured from engine facilities, are very erratic due to the inherent character of the combustion process. Thus, the method in [76] was employed to smooth the pressure signal, which is then averaged by 200 cycles for model validation. The residuals and the square of correlation coefficient (R-square) were used to evaluate the model accuracy. R-square can take on any value between 0 and 1, with a value closer to 1 indicting that a higher accuracy is achieved by the model. In addition, five in-cylinder parameters are chosen as criteria to verify the accuracy of the developed simulation model quantitatively; specifically, the cylinder peak pressure (pmax) and the corresponding crank angle (α1), the cylinder pressure at the spark timing (pcom), the pressure at the EVO (pEO) as well as the Indicated Mean Effective Pressure (IMEP). pcom, pmax and pEO are critical parameters that determine the diagram characteristic of in-cylinder pressure during the compression stage, combustion stage and expansion stage, respectively. IMEP is a valuable measure indicating an engine's power capacity independent of engine displacement, which should be considered for validating the combustion model.

Engine A
For Engine A, the NG is injected into the inlet manifold and mixed with the air forming a mixture that enters into the engine cylinders. The selected NG injection model keeps the air-fuel equivalence ratio at 1.2 for all the investigated cases. In order to develop the combustion model for the 2135 engine, the Wiebe model parameters (SOC, Δφ and m) for 100% load were calibrated by analysing the experimentally obtained HRR curve. The crank angle at the SOC and the end of combustion (EOC) at different operating conditions were obtained from the respective HRR curve. The Wiebe function exponent was obtained by analysing the HRR curve (which was estimated from the measured cylinder pressure) employing a curve fitting method [77]. Then, the Wiebe model parameters of other operating conditions (25% and 50% loads) were calculated by employing Equations (3)-(5).
The developed model was validated against the measured in-cylinder pressure for the Engine A at 25% and 50% loads. The parameters of the calculated single Wiebe combustion model at 25% and 50% loads are shown in Table 4. Figure 3 shows the simulated and measured in-cylinder pressure variations for these two operating points. The residuals for 25% load are between −1.3 bar and 0.8 bar with a R-square of 0.99, whilst the residuals for 50% load are within ±1.5 bar with a R-square of 0.98. As can be deduced from Figure 3, the simulated pressure diagrams sufficiently coincide with the respective measured data. The comparisons of pmax, α1, pcom, pEO and IMEP are shown in Table 5. The error of α1 is presented in the form of absolute difference (°CA), whilst those of pmax, pcom, pEO and IMEP are indicated by the relative percentage error (%). As can be deduced from Table 5, the relative errors of pmax, pcom, pEO and IMEP are below 3%, whilst the absolute difference of the peak pressure position is less than 1°CA. Thus, it can be inferred that the developed model is able to predict with adequate accuracy the performance of the Engine A.

Engine B
The NG fuel injection model for the Engine B is described by employing constant injection rates of 1.50 g/s and 2.33 g/s at 50% and 85% load, respectively. The pilot diesel injection timing is fixed at 5°CA before the TDC, whilst the injection mass rates are 0.21 g/s and 0.23 g/s at 50% and 85% load, respectively. In order to develop the combustion model for the YC6K engine, seven Wiebe model parameters (b1, SOC1, Δφ1, m1, SOC2, Δφ2 and m2) for four operating points (100%, 74%, 42% and 32% load) were calibrated by analysing the experimentally obtained HRR curve. The calibration procedure is described below.
The HRR obtained by the KiBox combustion analyser [68] is calculated from the measured pressure diagram by assuming that natural gas and diesel fuel burn proportionately according to their total injection mass ratio. Thus, the specific HRR of each fuel is not available, which makes it hard to obtain the respective Wiebe parameters. However, a number of parameters-particularly the Weight Factor and the SOC-can be obtained directly from the measured data. The Weight Factor b1, which is used to specify the energy percentages of pilot fuel, can be acquired by dividing the pilot fuel energy by the total energy of pilot diesel and NG. The energy provided by each fuel is calculated by multiplying the injection mass and the respective LHV. In order to reduce the Wiebe parameters, the pilot fuel and natural gas are assumed to start their combustion at the same crank angle. The SOC was obtained from the measured cylinder pressure diagrams (abrupt rising of pressure in comparison with the motoring cylinder pressure). As the natural gas is the primary fuel in the DF Engine and exhibits a much slower combustion speed than the pilot diesel, the total combustion duration obtained from the experimentally obtained HRR was used as the combustion duration for the natural gas. The remaining three Wiebe parameters, namely the combustion duration, the Wiebe exponent for diesel fuel and the Wiebe exponent of natural gas, can be obtained by a curve fitting method. Although the experimentally obtained heat release rate (HRR) diagram represents the

Residuals [bar]
R-square = 0.98 Residuals Zero Line overall HRR, the heat release contribution of pilot diesel and natural gas can be separated by the weight factor b1. After defining the exponential function shown in Equation (6), the remaining three Wiebe parameters were obtained by using the Fit function in MATLAB.
With the calibrated Wiebe combustion model parameters at the four operating points, an interpolation model was employed to predict the Wiebe combustion parameters at the other operating points. In this study, as the experimental data of Engine B were measured at propeller characteristics, there is a certain correspondence between the engine load and rotational speed. Thus, the calibrated values of the combustion model parameters are stored in a database as functions of the engine load. Then, quadratic interpolation is used for calculating the values of the respective combustion model parameters for each operating condition. The simulation model was validated against the measured in-cylinder pressures for the Engine B at 50% and 85% loads. Table 6 provides the parameters of the double Wiebe combustion model for the Engine B at four calibration points and two validation points. As can be deduced from Table 6, the SOC of the diesel fuel and natural gas is delayed from 2.9°CA ATDC to 4.0°CA ATDC when the operating condition varied from 100% load to 32% load. The combustion duration of the diesel fuel decreased from 14°CA to 7.5°CA, whilst that of natural gas reduced from 97°CA to 84.5°CA. The exponents for the diesel fuel and the natural gas Wiebe functions remained 0.5 in most operating points. The simulated and measured in-cylinder pressure variations are presented in Figure 4, which demonstrates that the simulated pressure diagrams exhibit a sufficient agreement with the measured ones. The comparisons of pmax, α1, pcom, pEO and IMEP are shown in Table 7. As the combustion in Engine B starts after the TDC, the pressure at TDC is selected as the pcom. As can be deduced from Table 7, the relative errors of pmax, pcom, pEO and IMEP are below 4%, whilst the absolute difference of the peak pressure position is less than 1°CA. Therefore, it is inferred that the developed model is capable of adequately predicting the performance of the Engine B, and can be used with fidelity for the knocking analysis presented in the next section.

Knocking Performance Parametric Investigation
The knocking tendency depends on the engine design and operating settings, which influence the end-gas zone temperature and pressure as well as the time spent at high levels of these two parameters prior to the flame arrival [24]. The compression ratio is one of the most key design parameters that determines the compression pressure and temperature, which exhibit apparent effects on the knocking tendency of gas engines [78]. The air-fuel equivalence ratio influences the flame propagation speed and energy density by changing the reactivity of the end-gas, thus affecting the knocking tendency [56]. The ignition timing is a critical parameter determining the combustion duration of NG-fuelled engines. If the ignition timing is advanced, the combustion duration would be extended thus providing the end-gas enough time for self-ignition, which renders the engine prone to knocking. Otherwise, the knocking phenomenon would be suppressed [79]. Based on the preceding discussion, the compression ratio, the air-fuel equivalence ratio and the ignition timing are selected for studying parametrically their influence on the knocking performance of the Engine A and the Engine B. For the DF engine, as the pilot fuel provides a multiple ignition source for the NG combustion, the pilot fuel mass variation affects the ignition timing and the flame propagation speed, which subsequently influence the knocking performance [80]. Thus, the effect of the pilot fuel mass variation on the knocking performance is also analysed for the Engine B.
As mentioned above, the knocking models were not validated by experimental data because pertinent testing might cause severe damage to the engine structure. To investigate the knocking performance parametrically, the knocking prediction models for the Engine A and Engine B were calibrated by virtually running the engine models at respective knocking boundaries. The spark timing of the 2135 engine is advanced by 5°CA at the 100% load operating condition, whilst the pilot injection timing of YC6K engine is advanced by 5°CA. A set of reference parameters including the crank angles at 75% and 50% of burnt fuel mass fraction, the air-fuel ratio and a reference K-value are obtained when the knocking probability reaches 100% by adjusting the model multiplier. Nevertheless, the Worret knocking criterion [30] is not a deterministic method that gives a certain indication whether knocking occurs inside the cylinder. It only indicates that the engine is more likely to operate at a hazardous condition of knocking occurrence when the knocking probability gets larger. According to reference [30], the knocking probability estimated by Worret criterion is approximately 20% when θk = θE. Thus, a knocking probability of 20% is set to be the knocking upper limit for the investigated parameters.

Engine A
The parametric investigation for the Engine A considers the compression ratio, the air-fuel equivalence ratio and the ignition timing variations within specific ranges around their baseline values, as shown in Table 8. The parametric settings are evenly distributed around respective original values in Table 1. In this section, the effects of compression ratio on the knocking phenomenon for the Engine A is studied. All the parameters except for the compression ratio remain the same with those at the engine nominal conditions. Figure 5a shows the effect of the compression ratio (CR) on the IMEP and the cylinder peak pressure for the Engine A. The IMEP varies from 4.9 to 5.72 bar achieving its maximum value (5.76 bar) when the CR equals 14. As it was expected, the cylinder peak pressure increases with the compression ratio increase. The peak pressure upper limit for this engine is 65 bar, which indicates an upper CR limit at around 14. As it is deduced from Figure 5b, the investigated SI natural gas engine runs into knock risk area (Pk > 20%) for CR values above 11. The knocking onset position is presented in Figure 5c (knocking position equal to 0 demonstrates the engine operating without knocking). As can be deduced from Figure 5c, the knocking onset is advanced from 1°CA BTDC (−1 °CA) to 7.5°CA BTDC (−7.5°CA) as the compression ratio increases from 10 to 16. Figure 5d shows the cylinder unburnt mass fraction (unburnt zone mass over the total mass) for the investigated CR range. The cylinder UMF varies from 0% to 75% as the CR increases from 7 to 16. UMF equal to 0 represents knocking free engine operation. Generally, a higher UMF value means that more energy will be released in a short period of time, thus indicating a stronger knocking phenomenon.

Air-Fuel Equivalence Ratio Effect on Knocking Performance
The knocking IMEP is used to evaluate the effect of the air-fuel equivalence ratio on the knocking onset. For a specific compression ratio, the knocking IMEP is obtained by increasing the injected natural gas mass while keeping a fixed air-fuel equivalence ratio until the knocking occurs. Figure 6 presents the knocking IMEP variation with the air-fuel equivalence ratio for different values of the compression ratio (ε). As can be deduced from Figure 6, for a given compression ratio value, the knocking IMEP varies with the air-fuel equivalence ratio without indicating a certain tendency. For a constant value of the air-fuel equivalence ratio, the knocking IMEP decreases with a larger compression ratio, indicating a higher tendency to knock. The nominal IMEP (at 100% load) of the 2135 engine is well below the knocking IMEP boundary of its original compression ratio (CR = 11), which indicates that the 2135 engine would operate knock free at any air-fuel equivalence ratios. In conclusion, the knocking boundary moves to lower IMEP values when the compression ratio increases and does not demonstrate a specific variation trend with the air-fuel equivalence ratio for the Engine A.

Ignition Timing Effect on Knocking Performance
The ignition timing for the Engine A is controlled by adjusting the spark timing. As the combustion model uses a single Wiebe function, the start of combustion (SOC) is assumed to be the same as the spark timing (the combustion ignition is assumed to occur at the spark timing). Figure  7a shows the effect of the SOC on the IMEP and the in-cylinder peak pressure. When the SOC is Nominal operation point retarded from −44°CA to −26°CA, the IMEP increases from 5.15 to 5.80 bar, whereas the cylinder peak pressure decreases from 48.2 to 35.5 bar. As presented in Figure 7b, the knocking probability decreases as the SOC is retarded, which means that the engine gets less prone to knocking. The more advanced ignition timing is associated with a higher knocking probability. In addition, the engine operates in the knocking safe area (Pk < 20%) when the ignition timing is retarded after −38°CA. Figure  7c indicates that the knocking occurs later when the ignition timing is retarded. In Figure 7d, the UMF decreases from 36% to 18% with the SOC retarded from −44°CA to −37°CA.

Engine B
In order to parametrically investigate the Engine B, the variation of compression ratio, air-fuel equivalence ratio, ignition timing and pilot fuel energy fraction were considered to be within the ranges shown in Table 9. Because the SOC can be apparently obtained by analysing the measured in-cylinder pressure, the combustion model of Engine B is developed by using the obtained SOC instead of employing the injection timing and the estimated ignition delay. Thus, the SOC variation is used to represent the ignition timing variation. As the LHV of the pilot fuel (42.65 MJ/kg) differs from that of the employed natural gas (47.18 MJ/kg), the energy proportion rather than the mass fraction is employed to investigate the effect caused by pilot fuel mass variation.   In order to study the effect of compression ratio on the knocking performance for the Engine B, all the parameters except for the compression ratio are kept the same with those at the nominal condition. As shown in Figure 8a, the IMEP of the Engine B increases from 18.2 to 21.5 bar when the compression ratio increases from 12 to 21. The peak pressure follows a similar trend with the IMEP, increasing from 92 to 185 bar. The peak pressure upper limit of the Engine B is 160 bar, which denotes a CR limitation at around 19. Considering the peak pressure upper limit and the knocking limit (Pk = 20%), the compression ratio of Engine B should be limited up to around 16.5 in order to operate safely. The knocking onset position is presented in Figure 8c, which indicates that the knocking occurrence is advanced from 25°CA ATDC to 3°CA ATDC with a higher compression ratio value. The unburnt mass fraction exhibits a monotonically increasing trend as indicated in Figure 8d, which indicates that the knocking gets stronger with a larger compression ratio.

Air-fuel Equivalence Ratio Effect on Knocking Performance
Similarly with the Engine A, the knocking IMEP versus the air-fuel equivalence ratio is used to evaluate the effect of the air-fuel equivalence ratio for the Engine B. Figure 9 presents the knocking IMEP variation with the air-fuel equivalence ratio for different compression ratio values. As can be deduced from Figure 9, the knocking boundaries under different compression ratios show the same descending trend with the increase in the air-fuel higher compression ratio. The nominal IMEP (at 100% load) for the Engine B is above the knocking IMEP boundary (when CR = 16.5), which denotes a knocking occurrence and suggests that the air-fuel equivalence ratio should be reduced from the original 2.1 to below 2. Nevertheless, the knocking prediction model for the Engine B was calibrated by virtually running the engine models at its knocking boundary, which is reached by advancing the injection timing by 5°CA. This was proposed by the laboratory manager and not validated by experimental data, resulting a manual estimation of the knocking prediction. Besides, the Worret knocking criterion is not a deterministic method that gives a certain indication of whether knocking occurs inside the cylinder. It only indicates that the Engine B is likely to operate at a hazardous condition of knocking occurrence at nominal operating conditions.

Ignition Timing Effect on Knocking Performance
As the employed combustion model is directly determined by the SOC rather than the injection timing and estimated ignition delay, the SOC is used to represent the ignition timing for the Engine B. For the current parametric investigation, the SOC is set varying from 6°CA BTDC to 8°CA ATDC. As can be deduced from the results presented in Figure 10a, both the IMEP and peak pressure of the Engine B decrease when the SOC is retarded from −6°CA to 8°CA ATDC. The knocking phenomenon is not likely to occur with a knocking probability smaller than 20% when the combustion starts after 3 °CA ATDC, as shown in Figure 10b,c. The knocking onset position and the unburnt mass fraction are presented in Figure 10c and Figure 10d, respectively. These results indicate that the knocking occurrence is delayed with an approximately constant knocking intensity when the SOC is retarded. The SOC is suggested to be limited after 3 °CA ATDC by jointly considering the peak pressure upper limit and the knocking limit.

Pilot Fuel Energy Proportion Effect on Knocking Performance
As the NG-air mixture is ignited by the multiple ignition sources provided by the pilot fuel, the mass variation of the pilot fuel may influence the combustion characteristics and the knocking performance of DF engines [80]. Thus, the effect of pilot fuel mass on knocking performance for the Engine B is investigated in this section. In order to maximize the influence of the combustion model parameters, the Wiebe functions parameters for the pilot diesel fuel and natural gas are calibrated or validated for six engine operating points, so that the engine cylinder pressure diagrams match the measured ones with adequate accuracy. The employed Wiebe parameters for the Engine B are shown in Table 6. Subsequently, the energy proportion of the pilot diesel fuel (b1 as in Equation (6)) varied, whilst the Wiebe functions parameters of both diesel and natural gas were kept the same with those provided in Table 6.
This parametric study considers the pilot diesel energy proportion varying from 5% to 15% under the six operating conditions shown in Table 6. The derived results are presented in Figure 11. It can be deduced from Figure 11a that the IMEP and the cylinder peak pressure slightly increase (by up to 5%) when the diesel energy proportion increases from 5% to 15%. In addition, the peak pressures at all the investigated cases are below the peak pressure upper limit (160 bar). Figure 11b presents the variations of knocking probability for the six investigated operating points, which exhibit a monotonically increasing trend at 100% load and remains zero for the other loads. The engine runs into knock risk area (Pk > 20%) 100% load when the pilot diesel energy proportion increases above 9%. In order to avoid knocking caused by the increase in pilot fuel energy proportion, the SOC must be retarded according to the analysis in Section 4.2.3. As can be deduced from Figure  11c, the knocking onset position at 100% load condition is advanced by around 5°CA when the diesel energy proportion varies from 8% to 15%. Figure 11d shows an increasing unburnt mass fraction from 48% to 57%, indicating a stronger knocking phenomenon with the increased diesel energy proportion.
(a) IMEP and peak pressure (b) Knocking probability (c) Knocking onset (d) Unburnt mass fraction Figure 11. Pilot fuel energy proportion effect on engine performance and knocking parameters (Engine B).

Discussion
As indicated by the results presented in the previous section, the increased compression ratio improves the IMEP and the peak pressure for both investigated engines. Meanwhile, higher compression ratios are found to aggravate the knocking occurrence, including advancing the knocking onset position and strengthening the knocking intensity. In order to avoid knocking operation, the compression ratio of Engine A should be limited up to 11, whilst that of Engine B is 16.5.
The air-fuel equivalence ratio effect on the IMEP knocking boundary for the Engine A is not explicitly indicated by the simulation results, whilst that for the Engine B exhibits a descending trend with the increase in the air-fuel equivalence ratio. The air-fuel equivalence ratio for the Engine B at a nominal operating point is suggested to be reduced from the original 2.1 to below 2. However, this requires additional experimental verification, which is out of the scope of the present study and will be considered in a future work.  The ignition timing delay has a suppressing effect on the knocking performance for both the investigated engines, i.e., retarding the knocking onset position and attenuating the knocking intensity. In order to operate safely, the ignition timings of Engine A and Engine B must be controlled after −38°CA and 3°CA, respectively. Within the knock-free operation range, the IMEP and the peak pressure of both the investigated engines decrease with the delay of ignition timing. Although the effect of pilot diesel fuel energy proportion is investigated under six operating points, the knocking phenomenon only occurs at the 100% load, which implies that the DF engine is prone to knock at a higher working load. The IMEP and the cylinder peak pressure for the Engine B increases within 5% range when the diesel energy proportion increases from 5% to 15%, which is attributed to the fact that more energy is released at the initial stages of the combustion with a larger weight factor b1. In addition, the knocking onset position is advanced with a stronger knocking intensity for the cases where the pilot fuel mass increases while keeping a constant total fuel energy input. A similar trend was experimentally observed in [80], where a greater pilot fuel quantity increases the knock tendency at high loads. It can be concluded that when DF engines operate at high load conditions, increasing the pilot fuel mass would not be recommended unless the SOC is retarded to suppress the knocking occurrence.
The obtained results provide suggestions for engine design settings (compression ratio) and operating settings (air-fuel equivalence ratio, ignition timing and pilot fuel mass) of the NG-fuelled engines by jointly considering the engine performance and the knocking behaviour. Nevertheless, the validity of the methodology presented is directly dependent on the existence of adequate experimental data in the examined operating conditions due to the absence of the physical and chemical process modelling that affects combustion. Although this type of model might not be adequate for determining the engine design settings (i.e., used for the design of a new engine/off-design investigation), it is of great value when used to explore the optimal settings of existing engine parameters.

Conclusions
This study aimed at the parametric investigation and comparative analysis of the knocking behaviour of an SI natural gas engine and a DF engine. Appropriate simulation models of the 1-D/0-D type were set up in the GT-Power software for both investigated engines. The combustion models were validated against experimental data, whilst the knocking models were calibrated by virtually running the engines at respective knocking boundaries. Subsequently, the engine models were used for the parametric investigation to study the effects of the compression ratio, the air-fuel equivalence ratio and the ignition timing for both engine performance as well as the effect of pilot fuel energy proportion for the DF engine performance.
The main findings of this study are summarised as follows.  Increasing the compression ratio results in higher cylinder peak pressure and IMEP for the SI natural Engine and the DF engine, and therefore increases the knocking occurrence tendency. The compression ratio upper limits of Engine A and Engine B are 11 and 16.5, respectively.  The knocking boundary for the SI natural gas engines does not exhibit an explicit variation trend with the air-fuel equivalence ratio, whilst that of the DF engine indicates a descending trend with the increase in the air-fuel equivalence ratio. The nominal air-fuel equivalence ratio for the Engine B is suggested to be reduced below 2, however experimental verification for this knocking limit is required.  The ignition timing delay exhibits a diminishing effect on the IMEP and peak pressure of both the investigated engines. The knocking occurrence for the Engine A and the Engine B could be suppressed by retarding the ignition timing after −38°CA and 3°CA, respectively.  A higher energy proportion (between 5% and 15%) of the pilot fuel increases the knocking tendency and intensity of DF engines at high loads. It is suggested to maintain the pilot fuel energy ratio less than 9% for Engine B.
The following recommendations were provided and can be used during the design and operation stages for SI natural gas engines and DF engines.
 Prior to attempting to improve the performance of an engine (SI NG or DF) by increasing the compression ratio or advancing ignition timing, knocking behaviour analysis must be performed by considering the knocking occurrence, onset position and intensity as well as the cylinder peak pressure upper limit.  The ignition timing of the SI natural gas engines and DF engines must be controlled precisely after the knocking limit in order to suppress knocking occurrence.  When operating at high load conditions, the pilot fuel mass of DF engine is suggested to remain at a low level to avoid deteriorating the knocking behaviour. This study presented a parametric investigation of the performance and knocking impacts of engine design and operation settings on two kinds of NG-fuelled engines. Therefore, this study provides the required numerical tools and framework for obtaining a better understanding and more thorough insights relating to the underlying parameters that influence the knocking behaviour of the NG-fuelled engine operation. The obtained results support the sustainability enhancement of NG-fuelled engines by identifying knock supressing settings, leading to safer operation and longer life cycle. The methodology presented is of great value to explore the optimal settings of the existing engine settings (inside the design space where the model parameters have been calibrated) and is ideal for digital twin applications and/or safety assessment or real-time condition monitoring. As a future study, it is expected that optimisation methods will be employed to obtain the optimal settings targeting a desired combustion behaviour and engine performance while holding the knocking tendency at a minimum level.