Modeling and Simulation
3.2 Main components in the modeling
The marine gas engine model studied in this project is developed by considering the following influential componens:
3.2.1 Combustion chamber
A mixture of air and fuel with an excess air ratio is captured in a non-fixed control volume, and it is assumed to be perfectly mixed. Several thermal zones are imposed for more accurate heat transfer calculation.
3.2.1.1 Heat transfer modeling
During the intake, compression, combustion, and exhaust, there is always an exchange of thermal energy between gases and boundaries.
Heat convection modeling in the main combustion chamber followed the method
proposed by Woschni [147] in Equation (3.1):
hc(Woschni)= K1P.8W0.8
B0.2TK2 (3.1)
whereK1 andK2 are constant, B bore and W is average cylinder gas velocity (m/s). The equation for the average in-cylinder gas velocity is:
W =C1Sp+C2VdTr
PrVr(P−Pm) (3.2)
C1 and C2 are constants, Sp mean piston speed (m/s), Tr unburned mixture temperature, P pressure, Pm motoring pressure, Prunburned mixture pressure, Vd
the total displacement and Vris the volume before combustion.
Woschni heat transfer model lumps the radiation portion, which is typically is estimated by the transmission of heat in the form of waves through space, into the convection portion inside the cylinder [148]. Moreover, radiative heat transfer from the hot burned gas in spark-ignition engines is tiny related to the convective heat transfer. Hence, Equation (3.1) is the only modeling theory for the heat transfer of the cylinder surfaces.
qconvection=hWoshni(Tgas−Twall) (3.3)
In Equation (3.3), Twallsplit into several zones.
The zones identified interior surfaces of the cylinder by: cylinder zone, piston zone, and head zone. Cylinder zone 1 is16 of the cylinder wall exposed when the piston is atBDC, cylinder zone 2 is13 of the cylinder wall exposed when the piston is atBDCand cylinder zone 3, covers the lower half of the cylinder wall exposed when the piston is atBDC.
The zones for piston surface and cylinder head are shown in Fig.3.2.
The temperature of the valve face surface is imposed and the total heat transfer is calculated for the valve surface area of:
areavalve =π
4Dvalve2+ 0.25Dvalvehvalve thickness (3.4) Where the average thickness of the valve head was chosen 16% of valve diameter.
3.2.1.2 Combustion
To predict the flame characteristics for the engine, the SI-turbulent flame model is employed. The model has the potential of predicting the combustion burning rate
Zone 1 Zone 2 Zone 3
(a)The three piston zones consisted of the piston zone 1 with 48%, piston zone 2 with 37% and piston zone 3 with 15% of the piston area.
Zone 1 Zone 2 Zone 3
valve face valve face
(b)The three head zones consisted of head zone 1 with 55%, head zone 2 with 35% and head zone 3 with 10% of the head area excluding valve faces.
Figure 3.2:Temperature zones of the cylinder head and piston top surface.
based on laminar and turbulent flame speed. It is assumed that there are two zones, burned and unburned zone, and the eddies of the unburned zone are entrained in the flame front at a turbulent velocity while the mixture of fuel and air is burning at laminar velocity. The rate of transformation of unburned to burned zone is calculated by Equation (3.5) and (3.6):
dMe
dt =ρuAe(ST+SL) (3.5)
dMb
dt = Me−Mb
τ (3.6)
where index b symbolizes burned classification, u unburned classification, and e entrained classification.
It is also assumed that the reaction is quite fast with a thin layer of the burning zone [9]. To calculate the laminar flame speed (SL), Heywood [109] recommended Equation (3.7) for several hydrocarbon fuels:
SL =
Bm+Bf(φ−φm)2Tu T0
αP P0
β
·fDilution (3.7)
where for natural gas fuel:
α= 0.68φ2−1.7φ+ 3.18 (3.8) β =−0.52φ2+ 1.18φ−1.08 (3.9)
and Bmis 0.490, Bfis -0.59, andφmis 1.390 [149].
fDilutionis the dilution effect and is employed as:
fDilution= 1−0.75A(1−(1−0.75A·Dilution)7) (3.10)
A is a multiplier, and Dilution is the mass fraction of the residuals in the unburned zone.
The time constant of combustion of fuel/air mixture entrained into the flame zone (τ)in Equation (3.6) was calculated by Equation (3.11):
τ= λ
SL (3.11)
with the Taylor length scale of:
λ=a Li
√Ret (3.12)
and turbulent Reynolds number:
Ret=ρuuL˙ i
µu (3.13)
Calculation of turbulent flame speed is done by Equation (3.14):
ST=b·u˙
1− 1
1 +c(RLf
i)2
(3.14)
where a, b, and c are constants. Moreover, u˙ and Li are turbulent intensity and turbulent length scale, respectively. They have been described in section 3.2.1.3. Calculation of surface area at the flame front (Ae) and flame radius (Rf) are presented in [150].
In order to calculate the concentrations of the species, chemical equilibrium is carried out for the entire lumped burned zone. This assumption is the state in which both reactants and products have no further inclination to change with time in the current pressure and temperature.
During the simulation, the main species and products of combustion are monitored in the main chamber by the relations between the composition of reactants and the composition of products by Equation (3.15):
x13
CnHmOlNk+n+m/4−l/2
F (O2+ 3.7274N2+ 0.044Ar)
−→
x1H+x2O+x3N+x4H2+x5OH+x6CO+x7N O+x8O2+ x9H2O+x10CO2+x11N2+x12Ar
(3.15)
Where x1 to x12 are mole fractions of the products and x13 is the mole of fuel, giving one mole of the product. The atom balances plus the products add up to unity provides six equations. To solve for the 13 unknowns, seven other equations are solved using the equilibrium provided by the following equations:
1
2H2↔H (3.16)
1
2O2↔O (3.17)
1
2N2↔N (3.18)
1 2H2+1
2O2↔OH (3.19)
1 2O2+1
2N2↔N O (3.20)
H2+1
2O2↔H2O (3.21)
CO+1
2O2↔CO2 (3.22)
Using the curve fitted from the JANAF thermodynamic table, the 13 equations of species are available and solvable in each pressure and temperature [151].
3.2.1.3 Flow model
The in-cylinder mean flow depends strongly on the chamber geometry, where the calculated values influence the burning rate. The sophisticated features of flow dynamics can be captured using detailed multidimensional Navier-Stokes models, but a simplified model of the combustion chamber is considered. The chamber is divided into three flow regions defined as the squish area above the piston crown, the cup volume, and the region above the cup lip. At each time step in each region, the mean radial velocity, axial velocity, and swirl velocity are calculated, taking into account the cylinder chamber geometry, the piston motion, and the incoming and exiting flow rates passes through the valves [152]. The turbulence model solves the turbulence kinetic energy equation and the turbulence dissipation rate equation. Two outputs of the model, turbulence intensity and turbulence length scale, are used directly in the turbulent flame velocity calculation in combustion modeling of Equation (3.5).
The turbulence length scale is calculated by:
Lm= Cµ
3 4.K32
(3.23)
Cµis constant. The detailed equations are presented in [152]:
3.2.1.4 Emission
Natural gas engines produce four main emission compounds throughout the combustion: NOX, UHC, CO2, and CO [153]. Generally, the concentration of these compounds differs in chemical equilibrium rather than the detailed chemical mechanism. Besides, implementing the detailed kinetics to find the emission level during combustion and post oxidation needs considering hundreds of reactions.
Therefore, a basic formation mechanism is provided for the pollutants.
3.2.1.4.1 NOX NOX is nitrogen oxide, NO, and nitrogen dioxide, NO2, but NOis the prominent oxide of nitrogen of the internal combustion engine [109].
Therefore, the mechanism of formation of NO proposed by Zeldovich have been extensively applied for description ofNOformation in engines [154]:
O+N2=N O+N (3.24)
N+O2=N O+O2 (3.25)
N+OH =N O+H (3.26)
The rate of formation ofNOvia reactions of (3.24), (3.25) and (3.26) is:
d[N O]
dt =k1+[O][N2] +k2+[N][O2] +k3+[N][OH]
−k1−[N O][N]−k2−[N O][N]−k3−[N O][H]
(3.27)
With the steady-state approximation for d[N]dt and equilibrium for other species, a simplified one-way equilibrium is achievable with the following reaction rates respectively:
k1 =F1·7.60·1010·e−38000
A1
Tb (3.28)
k2=F2·6.40·106·Tbe−3150
A2
Tb (3.29)
k3=F3·4.10·1010 (3.30)
F1, F2, F3, A1 and A2 are constant coefficients for tuning the modeling, and Tbis the burned zone temperature.
3.2.1.4.2 UHC The implemented equations forUHCmodeling to calculate the amount of methane slip is separately discussed inchapter 4.
3.2.1.4.3 CO2 CO2 is a greenhouse gas in the exhaust gases fromSIengines as a main product of the combustion. Natural gas has chemical properties with a high H/C ratio of around 3.7. Thus, by changing the fuel from diesel with a ratio of 1.8 to natural gas, an immediate reduction ofCO2 to half is achievable.
As a result, the gas engines produce less CO2 and even 20% less than similar gasoline engines [72]. Thus the relative amount of this emission depends on fuel properties than the combustion mechanisms. Therefore, the output of the chemical equilibrium is used for the validation of the combustion, and the formation of this emission is not discussed further.
3.2.1.4.4 CO At lean equivalence ratios, like for the lean burn spark-ignition engine, a low amount ofCOcan be achieved [155,156]. Moreover, the formation ofCOemission is high primarily when the engine operates close to stoichiometric, or during warm-up when the wall is cold. Since the primary output from theCOis not verified, any further presenting of the modeling and output is neglected.
3.2.1.5 Knock
In spark-ignition engines, knock represents one of the most critical issues to reach optimal thermal efficiency. Therefore, to deal with any probable knocking occurrence, an air/fuel ratio andEGRbased methodology aimed at evaluating the knock phenomenon. In the model presented here, the knock intensity is evaluated by the Kinetics-Fit-Natural-Gas model proposed by Gamma Technology and is based on the kinetic reaction mechanism of natural gas [157,158]:
KI =M.u(α)VTDC
V(α)exp(−6000 T(α) )·
max[0,1−(1−φ(α))2].Iindex
(3.31)
Where M is a multiplier and can be found during the validation step. Induction time integral, Iindex, is calculated by equation (3.32):
Iindex= Z 1
τdt (3.32)
and the induction time,τ, is defined by equation (3.33) : τ=M1−
"
10−9exp (18659 M2T )(M N
100 )0.978·
(F uel−0.578)(O2−0.28)(Diluent0.03)
# (3.33)
where MN is the fuel methane number, M1 and M2are multipliers, equalized to 1, and diluent concentration is the sum of concentrations ofN2, CO2, and H2O.
VTDC is cylinder volume at the top dead center andφis equivalence ratio of the unburned zone.u(α)is the percentage of cylinder mass unburned.
3.2.2 Power transmission
Sudden loading and unloading on the engine lead the engine to speed fluctuation and instability. Any idea to improve the engine stability is highly pertinent to the engine mass moment of inertia and the connected rotating components as power transmission [159,160]. The degree of change of the system frequency regarding the external forces is inversely proportional to the magnitude of the inertia:
Tb(t) =Ts(t)−I.ω˙ct(t) (3.34)
where,ω˙ct(t) is the instantaneous crank-train acceleration andI mass moment of inertia. Tb(t) represents the torque available at the flywheel, after accounting for all friction and attachment losses as well as the acceleration of the crank train inertia.
The entire propulsion shaft is modeled as a single rigid, and the shaft torque, Ts, is calculated by:
Ts(t) =Ti(t)−Tf(cyc) +Ta (3.35)
If there is any additional instantaneous torque of the other attachments, they can be added to the right of Equation(3.35) as Ta. Tf(cyc) is friction torque for the current cycle, and Ti(t) is the indicated torque and depends on the combustion chamber.
The pressure inside the main chamber, calculated by the flame modeling, provides forces on the x-y direction on connecting rod as:
Ti(t) =
9
X
1
(π
4)B2Rcrank(∆P(t) sinθi(t)−∆P(t) cosθi(t) tanαi(t)) (3.36)
j x ~ αi(t)
θi(t)
RCrank
Connecting road Bore (B)
Figure 3.3: The instantaneous indicated torque is a function of pressure in each crank angle and the in-cylinder surface area.
θi(t)is the instantaneous angle of crank i in degree andαi(t)is the instantaneous angle of connecting rod i in degree. The schematic shown in Fig. 3.3, indicates the pressure on the piston surface and the distribution of the forces.
The power transmission included the flywheel, crankshaft, main shaft, and gearbox to transmit the power to/from the engine. All are assumed to be rigid. Thus deformation of the crankshaft, corresponding to the immediate difference between the engine and load torque, is neglected. The inertia concerning the water resistance and the propeller is added to this model as well.
3.2.3 Controlling system
A proportional-integral-derivative controller (PID) consists of a proportional unit (P), an integral unit (I), and a derivative unit (D). The PID formula follows
Equation (3.37), and the schematic of the interconnections is shown in Fig.3.4.
y(t) =KPe(t) +Ki Z t
0
e(t)dt+Kdde(t)
dt (3.37)
where y(t) and e(t) are the process input and error, respectively. The PID controller is used when the system has a linear response with dynamic characteristics. In a complex system with more degrees of freedom and nonlinear dynamic processes in calibration, new calibration approaches are needed [161, 162].
In this study, for the engine modeling, three regular PI controllers (PID with a coefficient of D=0) were required:
• The fuel flow regulator with engine speed feedback
• TheVTGvane position with the feedback of air-fuel ratio
• The throttle opening angle with engine loading feedback
P kpe(t)
I kiR e(t)dt
D kde(t) dt
P Process
Setpoint +P error + +
+
Output
−
Figure 3.4: PID controller block diagram. Three parameters (KP, Ki, and Kd) can be manually or automatically tuned based on the setpoint and output.
3.2.4 Intake and exhaust manifold
Since the fluid flow is important in calculating the dynamic delay and consequently the turbocharger response, the equation of conservation are involved the conservation of momentum as well as Equation (3.38) and (3.39) [157]. The dimensions for all the pipes and junctions are implemented based on the designed
are chosen.
dm
dt = X
boundaries
˙
m (3.38)
dme
dt =−ρdv
dt + X
boundaries
( ˙mH)−hAs(Tfluid−Twall) (3.39)
dm˙
dt = dpA+P
boundaries( ˙mu)−4Cfρu|u|2 dxAD −KP(12pu|u|)A
dx (3.40)
Wheremis mass of volume,KPpressure loss coefficient,Dequivalent diameter, andm˙ is the boundary mass flux. Ais the cross-sectional flow area andAsis the heat transfer surface area.
In all discretized volumes, Equations (3.38) and (3.39) yield the mass and energy in each volume. With the available volume and mass, the density can be calculated.
Afterward, the equations of state define density and energy as a function of pressure and temperature, and the solution will be continued iteratively on pressure and temperature until they satisfy the density and energy already calculated for this time step.
3.2.4.1 Fuel system
To address the fuel system response, it is imperative for the modeling to pursue the same fluid dynamics approach as the air in the intake and exhaust manifold. This is attained by adopting the use of conservation of mass, energy, and momentum equations, providing a high fidelity of fuel flow model, cost-effective, and ensuring the dynamic system response. To simulate the fuel delivery system dynamics, all fundamental components of this system, including the fuel tank, orifices, fuel lines, and fuel valves, are assembled. For simplicity, the fuel tank and fuel pump are assumed to be constant pressure sources. The fuel, natural gas, is mixed with the air just upstream of the intake valves.
3.2.5 Intake and exhaust valves
The intake and exhaust valve capacities ultimately restrict the total flow of the engine. A feasible approach for compressible flow over a flow restriction is proposed by Heywood [109] with a pressure upstream and downstream of the valve:
˙
m=AeffρisUis=CDAR Pu
√RTu(Pd Pu)1/γ
r [ 2γ
γ−1(1−Pr)γ−1γ ] (3.41)
0 100 200 300 400 500 600 700 0
0.2 0.4 0.6 0.8 1
Crank angle
Lift(Normalized)
Inlet valve Exhaust valve
Figure 3.5: The intake and exhaust valves lift in an engine cycle. The values are normalized to be fitted in the maximum amount of one.
where u is upstream and d is downstream stagnation point. Here, CD is an experimentally determined discharge coefficient of the valve lift presented in Fig.
3.5.
For chocked flow, the equation modifies to:
ρis=ρo( 2
γ+ 1)γ−11 (3.42)
Uis=p
γRTo( 2
γ+ 1)12 (3.43)
3.2.6 Intercooler
Compressed air after turbocharger has higher internal energy with higher pressure and temperature. The drawback of the boosted fluid is now the low density due to the high temperature. For increasing the mass flow, it is essential to decrease the fluid density. An intercooler wastes the additional heat by multiple pipes. The infinite number of pipes is considered to provide an infinite heat sink for the modeling, and the outlet temperature leaving the intercooler equals the intercooler wall temperature and is equal to the measured value from the engine data. Calculation of total heat flux is performed by Equation3.44, where hgis the convection heat coefficient and is calculated using the Colburn analogy [163] by Equation3.45.
qconvection=hg(Tgas−Twall) (3.44) hg= (1
2)CfρUeffCPP r(−23 ) (3.45) where Cfis friction factor, Cpspecific heat, and Ueffis effective fluid velocity.
3.2.7 Turbocharger
A turbocharger consists of three components: The turbine, which converts the high enthalpy gas to rotating energy, the compressor, which utilizes the excess power from the turbine, and finally, the connecting shaft. In developing the models for the turbocharger, a performance map is an alternative solution. Utilizing this method, the compressor and turbine information as a function of speed, pressure ratio, mass flow rate, and efficiency are implemented in the format of look-up tables. Turbocharger speed and pressure ratio are predicted at each time-step, and two other unknowns are taken from the look-up table [164]. The calculation starts with the predicted pressure ratio, calculating total temperature by Equation (3.46) and Equation (3.47):
Ttotal,in =Tin+uin2
2cP (3.46)
∆hs=cPTtotal,in(P Rγ−1γ −1) (3.47)
where uin is inlet velocity. The outlet enthalpy will be calculated by Equation (3.48) and the compressor power by Equation (3.49) provides the torque of the compressor subsequently:
hout=hin+∆hs
ηs (3.48)
P = ˙m(hin−hout) (3.49)
Subtext script in, out and s stands for inlet, outlet and isotropic, respectively. By Equation (3.50), the new calculated rotational speed is provided:
∆ω=∆t(Tturbine−Tcompressor−Tfriction)
I (3.50)
A table of friction coefficients was considered for the turbocharger shaft for different rotational speeds, as Tfriction.
For the turbine, the following equations are employed instead of Equation (3.47) and (3.48):
∆hs=cPTtotal,in(1−P R1−γγ ) (3.51)
1
8 9
2 3
7
5 6
4
1 Compressor 2 Intercooler 3 Inlet manifold
4 Fuel 5 Block 6 Torque
7 Exhaust manifold 8 Turbine 9 Turbo Shaft
Figure 3.6: Engine modeling schematic. All essential elements influential on flow and dynamic are implemented in the engine model using two-zone zero-dimensional modeling for the combustion, one dimensional for pipes and connections, and a look-up table for the turbocharger.
hout=hin−∆hsηs (3.52)
The computation will be repeated until the predicted parameters reach convergence.
3.2.8 Boundary condition of inlet and outlet
Sea reference condition states the free water surface boundary condition on an ocean. A static atmospheric pressure, temperature, and zero amplitude are utilized for the compressor inlet and the turbine outlet.