Accepted Manuscript

Implementation of EOQ and Lambert W function in 1-D Engine

simulation model for optimizing fuel injection in GDI engine

Pablo Santirso, Stephen Samuel

PII: S0307-904X(18)30407-4
DOI: https://doi.org/10.1016/j.apm.2018.08.018
Reference: APM 12432

To appear in: Applied Mathematical Modelling

Received date: 29 November 2017

Revised date: 13 August 2018
Accepted date: 20 August 2018

Please cite this article as: Pablo Santirso, Stephen Samuel, Implementation of EOQ and Lambert W
function in 1-D Engine simulation model for optimizing fuel injection in GDI engine, Applied Mathematical
Modelling (2018), doi: https://doi.org/10.1016/j.apm.2018.08.018

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service
to our customers we are providing this early version of the manuscript. The manuscript will undergo
copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please
note that during the production process errors may be discovered which could affect the content, and
all legal disclaimers that apply to the journal pertain.
 ACCEPTED MANUSCRIPT

Highlights
• A novel and unique approach based on EOQ for optimizing fuel injec-
tion is implemented

• A scheme for coupling crank angle resolved model with EOQ in MAT-

T
LAB SIMULINK

IP
• Predictive combustion model with EOQ is sufficient for optimizing fuel
quantity

CR
• Proposes EOQ with Lambert W function for model based engine cali-
bration approach

US
AN
M
ED
PT
CE
AC

1
 ACCEPTED MANUSCRIPT

Implementation of EOQ and Lambert W function in

1-D Engine simulation model for optimizing fuel
injection in GDI engine

T
Pablo Santirso, Stephen Samuel1,

IP
Oxford Brookes University, School of Engineering, Computing and Mathematics,
Wheatley Campus, Oxford, United Kingdom

CR
Abstract
US
The aim of this study was to implement Economic Order Quantity method
(EOQ) together with the Lambert W function in a 1-D engine simulation
AN
model in order to develop a fuel control strategy for a Gasoline direct injec-
tion (GDI) engine. Previous work of the co-author demonstrated the pos-
sibility of optimizing fuel injection quantity in GDI engine using the EOQ
that is commonly used in supply chain of perishable products. This work
M

extends the previous work and implements it in a 1-D, crank angle resolved,
engine simulation model for the application of model based calibration pro-
cess. The present work uses a validated engine simulation model, which is
ED

based on predictive combustion modelling approach, and couples the 1-D

engine simulation model with SIMULINK to add the evaporation, wall- wet-
ting and heat transfer models. It employs FORTRAN subroutines to modify
PT

the internal code of the 1-D simulation software in order to add crank angle
resolved evaporation model. Finally, EOQ with Lambert W function was
added to the model using MATLAB with special attention to the decimal
CE

control for the solution. This study demonstrated that EOQ and Lambert
W functions together are a suitable method to develop fuel control strategy
for a model based calibration procedure when implemented in crank angle
resolved 1-D simulation model.
AC

Keywords: Economic Order Quantity, Lambert W function, 1-D engine

simulation model, Predictive Combustion modelling, Fuel injection,

Email address: [email protected] (Stephen Samuel)

1
Corresponding Author

Preprint submitted to Applied Mathematical Modelling August 30, 2018

 ACCEPTED MANUSCRIPT

Gasoline Direct Injection Engine.

1. INTRODUCTION
One of the most important decisions taken by the engine design teams for

T
optimizing the performance of the engine and for meeting the emission tar-

IP
get is choosing appropriate injection strategy for gasoline injection engines.
This, in turn, decides the strategy for developing fuel injection mapping of

CR
the engine and powertrain calibration process. Improving the efficiency with
the Engine Control Unit (ECU) of the engine enables the designer to develop
fast and relatively less expensive modifications, which are easily adaptable to
different types of engines. Likewise, 1-D simulation modelling tools allow the

US
designer to develop complex algorithms for improving the fuel consumption,
numerically observe the behaviour of the engine and to find possible errors
before developing the testing schemes for calibration optimization.
AN
Historically, the fuel injector had been placed in the intake system of a
gasoline engine. The quantity of the fuel was metered based on the desired
air-to-fuel ratio in the cylinder considering wall-wetting and other losses.
M

Later, this system evolved into a multi-point port fuel injection (MPFI)
system where individual injectors were employed for each cylinder located
immediately before the intake valve. This system could meter the quantity
ED

of the fuel more accurately than the single point injectors for all cylinders.
Modern automotive engines have their injectors in the combustion chamber,
where fuel metering is done more precisely than the MPFI system and there
PT

is more freedom to vary the air to-fuel-ratio for different operating conditions
if required.
CE

However, this additional flexibility offered by the direct injection system

also brings additional challenges in relation to quantifying the amount of
fuel that is not available for combustion processes. For example, in a single
point throttle body injection system (TBI), the temperature of the air in
AC

the manifold upstream and downstream of the injection system will indicate
the amount of fuel evaporated in the manifold. Therefore, it was possible to
correlate the mixture temperature with the quantity of the fuel lost in the
system for that cycle [1, 2, 3].

3
 ACCEPTED MANUSCRIPT

Similarly, for MPFI engines the parameters required for estimating the
portion of the fuel available in the cylinder could be estimated using the tran-
sient wall wetting model [4, 5, 6, 7, 8, 9]. Additional measurements could be
carried out for validating the model. The accuracy required for predicting
the correction is very critical since the time available for fuel to evaporate

T
and mix with the incoming air stream is not long enough when compared to

IP
TBI system. However, the temperature and pressure, mass flow rate of air
in the manifold, parameters required for describing the fuel puddle could be
measured or observed using various methods in order to model the transient

CR
fuel correction [4, 10, 9, 11, 5]. 3-D computational fluid dynamic models
could be employed to visualise and also gain detailed insight into fuel pud-
dle and evaporation dynamics [12, 13] at different operating temperatures.

US
These studies could be validated using experimental observations [14, 15]
such as Shadowgraph [13] since the physical location of interest is accessible
for instrumentation.
AN
However, in GDI engines the time available for fuel injection, evaporation
and mixing is significantly lower than that of MPFI engine [5]. In addition,
most of the variables used for estimating fuel quantity available in gaseous
M

form for combustion and fuel lost during the injection and evaporation pro-
cess change constantly as a function of crank position. The variables such
as the charge temperature, fuel jet penetration distance, fuel velocity will
ED

change as a function of crank position. It is also difficult to use any direct

measurement or observation for verifying these instantaneously [16, 17, 18,
19] Therefore, developing physical models, which will inherently reflect the
main characteristics of the injection, wall wetting, evaporation and mixing, is
PT

paramount for predicting the amount of fuel available for combustion, which
is directly linked to the torque output of the engine.
CE

The concept of model based engine calibration is not a new field of study
[30, 69]; however, it has received more attention recently for the above stated
reasons. Some of the models rely heavily on experimental data and statis-
AC

tical inference, some rely on a semi-empirical model, [20, 21, 49] and some
are based on experimental data and statistical inference for training a neural
network based model [22]. The ultimate aim of these models is to predict the
injection quantity for achieving the required torque output from the engine
and emission levels, since this torque value determines the effort required to
drive the vehicle at a desired vehicle speed. If the predicted quantity of the

4
 ACCEPTED MANUSCRIPT

fuel is different from the quantity of the fuel actually available for combus-
tion at the start of the combustion process, it will have an adverse effect
on vehicle drivability, overall emission levels and vehicle performance. These
models should be able to predict the fuel quantity for steady-state as well as
transient operating conditions [23, 20].

T
IP
1-D simulation models to develop sufficient details for engine calibration
is one of the approaches used by engine manufacturers to reduce the calibra-
tion time and generate calibration sheets for achieving required performance

CR
and emission control process [24] using software such as GT-Power. These
types of simulations also use neural network for optimization using design
of experiments done numerically [25].These types of models carry inherent

US
limitations because they are not capable of capturing flow field characteris-
tics in the engine combustion chamber and have only limited capability for
capturing any physical process that is 3-D in nature. However, the use of 3-D
computational model for developing engine calibration sheets has not been
AN
fully realised so far mainly because of the computational cost. Therefore,
computational fluid dynamic (CFD) models are commonly used to study the
individual process such as wall wetting in the fuel injection engines [12]. In
M

addition, coupling a 1-D model with the 3-D model and using CFD for only
the critical components and control volume area where 3-D details are re-
quired for capturing physical processes is being tried by various researchers
ED

[26, 27]. However, validating a 3-D model for a combustion chamber and
flow field characteristics in the combustion chamber is a major challenge for
production type engines.
PT

Therefore, if detailed physics based models [28, 29] to capture the charac-
teristics of fuel injection and spray formation, wall wetting and evaporation
are chosen then it is possible either using 1-D or 3-D to estimate the fuel loss
CE

in the system and the amount of fuel available for combustion. Therefore,
it is possible to correlate the injected fuel mass with the torque output of
the engine provided a suitable crank angle resolved combustion model [30,
AC

31] is chosen. The choice of the combustion model is also crucial for relat-
ing injected fuel quantity with the torque output of the engine. If a forced
combustion model such as Wiebe function [32] based models are chosen then
the energy release characteristics are already imposed for a given quantity
of fuel and engine state. Any changes in engine temperature that will affect
the wall wetting and evaporation characteristics of the fuel which determine

5
 ACCEPTED MANUSCRIPT

the quantity of the fuel available for combustion cannot be linked to output
torque if Wiebe combustion model is used, unless a detailed map of Wiebe
constants [33, 32] which include the corrections for these types of changes are
used. However, if a suitable predictive combustion model is used which can
take the flow field characteristics and chemical kinetics of the fuel into ac-

T
count then correlating the injected mass with the torque output will be more

IP
accurate provided the physics model of injection, wall wetting and evapora-
tion are correctly coupled [14, 30, 34, 35, 34, 37, 38].

CR
Hence, a model based calibration has been considered a viable option by
various researchers for reducing the calibration time and to develop design
of experiments for capturing essential experimental data in order to optimize

US
the performance of the engines [39, 40, 41, 42, 43, 44]. Moreover, crank
angle resolved engine models [31] with varying degree of complexities [45]
in conjunction with statistical insight derived from experimental data and
Artificial Neural network based models are getting more attention for vehicle
AN
calibration and optimization recently [37]. These are for predicting cycle-by-
cycle indicated torque [46], detecting and estimating spark timing [47, 48] to
employ fuel injection strategy [50, 51, 45, 52, 53, 54], for cold start control
M

strategy [55] and for effective implementation of other controls such as knock
and for monitoring the performance of catalytic converters.
ED

Recently researchers are exploring the possibilities of combining various

simulation packages to capture essential characteristic of the flow field and
combustion but at the same time using the strength of different modelling
approaches to reduce the run time significantly. For example, Rask et al. [24]
PT

have demonstrated the use of GT-Power, a 1-D gas dynamic engine mod-
elling approach with experimental data for developing high fidelity model for
engine calibration. Sika et al. [56] have tried to couple GT-Power with MAT-
CE

LAB Simulink for developing control strategies. Similarly, Omekanda et al.

[28] have used 1-D thermal model of the engine using MATLAB Simulink and
Python in order to develop thermal management control. Recently Turkson
AC

et al. [25] have attempted to couple GT-Power with MATLAB for developing
optimization strategy.

However, the present study has taken completely different approach when
compared to the existing optimization strategies available in the public do-
main. The previous work of the author [57] demonstrated that fuel injection

6
 ACCEPTED MANUSCRIPT

and optimization could be treated as a supply chain problem and therefore,

the concepts such as Economic Order Quantity (EOQ) [73,59] commonly used
in determining the cost in the supply chain discipline could be used for op-
timizing the injected quantity of the fuel in Gasoline direct injection engine.
Recently, researchers have already demonstrated that borrowing concepts

T
from thermodynamics for modelling inventory systems and from economics

IP
and industrial engineering for modelling thermodynamic systems can enable
us to gain more insight for optimizing the performance of systems [70,71,72].
Therefore, in order to develop schemes for reducing the development time,

CR
especially for fuel injection and calibration process, the present study aims
to implement EOQ based optimization in a 1-D crank angled resolved engine
model.

2. SCHEME OF THE WORK US

This study is a continuation of the research work of Ventura et al. [57]
AN
which demonstrated that EOQ and Lambert W function could be used for
optimizing fuel injection process in GDI engine. It aims to implement EOQ
and Lambert W function by developing sub-model that can be linked to
M

crank angle resolved 1-D engine simulation software. The sub-models such
as fuel spray model, spray evaporation models and wall wetting and heat
transfer are developed in SIMULINK. A FORTRAN subroutine is used for
ED

interacting with 1-D simulation model in order to estimate crank angle re-
solved quantities. A script for the EOQ and a script to solve the Lambert W
function are developed in MATLAB based on Disney et al. [58]. Finally, the
results are correlated with the cylinder pressure data of the engine simulated
PT

and the EOQ function for different engine operating conditions.

3. ECONOMIC ORDER QUANTITY

CE

EOQ model is a well known model in inventory management [73, 74] and
it is used for estimating near-optimal order quantities. A brief summary of
AC

EOQ concept commonly used in supply chain management for optimizing

the cost [59] is given here for the sake of completeness. This work uses EOQ
principle for perishable products such as melons proposed by [58]. The model

7
 ACCEPTED MANUSCRIPT

minimises the cost of this supply chain using the following equation.
 
KD D αQ
TC = + DV − (pV exp(−βtj ) exp(αtr )) 1 − exp + cD + Cj
Q Qα P
(1)

T
Where T C is the total cost for harvest period, K is the batch transfer time
in currency, D is the total annual harvest, c is the scalar multiplier, Q the

IP
transfer batch size in cartons, V is the maximum value that a carton can
achieve at time t = 0, p is the picking rate in cartons/hour, α is the de-

CR
terioration rate of the product per hour outside the coolant storage, tr the
transfer time in hours from field to the cooling facility, β deterioration rate
of the value in the cooling facility, tj the time in the cooling facility and Cj
the cost of transportation to the retailer. If the cost is minimized and the
equation is adjusted:

Q=
p 
US
− k exp(
Qα
)−
p
(2)
AN
α p α
K
k = exp(αtr + βtj ). (3)
V
M

If the Lambert W function is applied in the equation, then:

   
∗ p kα 1
Q =− W−1 − +1 (4)
α pe e
ED

The solution exists for equation (4) if p, K and α are positives and V >
exp(αtr +βtj )
p
Kα.
In this study SIMULINK is used to estimate the solution of
PT

Lambert W function for the EOQ equation (4).

3.1. Analogy with Internal combustion engines

CE

Ventura and Samuel [57] established an analogy between the supply chain
of a perishable product and the fuel injection phenomenon in a GDI engine.
The analogy proposed by Ventura and Samuel with additional details linked
to the present work is shown in Figure 1.
AC

The final output when the equations are solved is the decision variable
Q*, it will be the mass of fuel that the injector has to reduce from the base
fuel mass quantity. The final equation in the GDI engine will be:
   
∗ p kα 1
Q =− W−1 − + 1 tj (5)
α pe e

8
 ACCEPTED MANUSCRIPT

T
IP
CR
US
AN
M
ED
PT
CE
AC

Figure 1: Introducing EOQ and Lambert W function in crank angle resolved 1-D engine
simulation software (part of the figure is redrawn from [57]

9
 ACCEPTED MANUSCRIPT

T
IP
CR
US
AN
M
ED
PT
CE

Figure 2: Cylinder pressure with the variables required for implementing EOQ parameters
in 1-D, crank angle resolved engine model
AC

10
 ACCEPTED MANUSCRIPT

Where tj is the duration of injection. Figure 2 shows the phases of the

injection, combustion and the main variables used in the EOQ function. For
estimating deterioration rate in evaporation, α, using EOQ it is essential to
include the spray evaporation model and wall-wetting model and heat trans-
fer model to the main engine model. It is also essential to define injection

T
event, start of combustion event in order to estimate the amount of fuel evap-

IP
orated before the start of combustion. The injection rate, p, was estimated
with the mass injected and the time of injection. The heating value, v, of the
fuel is assumed as constant and is 4.395 x 107 J/kg , tj and tr are the dura-

CR
tion of the injection process and wall wetting process and β the deterioration
rate of the heat transfer process. Therefore, the total cost of evaporation and
wall wetting, K, is the amount of energy lost because of the unburned fuel
quantity in kJ.
US
K = mU B ∗ v (6)
where mU B is the mass of fuel unburned fuel. Now, it is possible to calculate
AN
the Q∗ as shown in equation (5), this will be the amount of fuel that should be
subtracted from the original mass of fuel to obtain the most optimal amount
of fuel to be injected. Therefore, the iteration process will continue until
the system is stabilized and Q∗ become negligible as shown in Figure A1 in
M

Appendix A.

4. PHYSICS MODELS
ED

The EOQ function requires crank angle resolved quantities from spray
evaporation, wall wetting and heat transfer and combustion models. The
PT

evaporation of the spray is the rate of amount of fuel, which evaporates from
liquid to gas state in the spray jet. In addition, this spray can impinge on
the piston surface and create a puddle of fuel; this puddle has a different
CE

evaporation rate, which will be calculated in the wall-wetting model, finally,

the heat transfer model finds the amount of energy transfer by the walls of
the cylinder.
AC

4.1. Spray evaporation model

The direct injection system injects a jet of spray into the combustion
chamber. The evaporation rate of the spray is governed by the equations
proposed by Yoshizaki et al. [60]. These equations assume that the spray is
constituted by individual packages as shown in Figure 3.

11
 ACCEPTED MANUSCRIPT

T
IP
CR
Figure 3: Individual packages of spray as illustrated by model

US
The evaporation ratio of each droplet in a spray jet based on Lefebvre [61]
as a function of Sauter Mean diameter, D, (SMD) and the fluid properties
AN
are as follows :
∂mD Ka
= 2πD ln(1 + Bm ) (7)
∂t Cp,a
M

where Ka and Cp,a are the conductivity and specific heat at constant
pressure of the air and Bm the mass transfer number. Similarly, the number
of droplets in the spray could also be estimated by knowing the volume of
ED

fuel injected, SMD and the fuel density.

minj !
∂mD ρf uel
ṁspray = . 3 (8)
PT

∂t 4
π.(D/2)3
And the relationship between the droplet diameter and the temperature
of the droplet is:
CE

∂D 4Ka ln(1 + Bm )
= (9)
∂t ρf Cp,a D
∂T ∂mD Lf BT
AC

= ( − 1) (10)
∂t ∂t mD Cp,f Bm
Where ρf , Cp,f and Lf are the density, specific heat and latent heat of fuel
respectively and BT is the heat transfer number, the vapour concentration,
YF s , at the surface of the droplet. According to [62]:

12
 ACCEPTED MANUSCRIPT

YF s
Bm = (11)
1 − YF s
Cp,a (Tair − Tdroplet )
Bt = (12)

T
1 − YF s
 −1

IP
P Ma
YF s = 1+( − 1) (13)
PF s Mf

CR
Where Ma and Mf are the molecular weight of the fuel and air, P is the
pressure inside the cylinder and PF s is the fuel vapour pressure. Fuel vapour
pressure is estimated using the Antoine and Clausius-Claypeyron relation
[63]:

US B
PF s = 10(A− C+T ) (14)
AN
Where A, B and C are the Antoine coefficients. And the latent heat as
a function of temperature [62]is:
 −0.38
Tcr − Tdroplet
Lf = Lf,ref erence . (15)
M

Tcr − Tbn
Where Tcr and Tbn are the critical temperature and boil temperature. By
combining these principles, according to Hiroyasu et al. [64] the SMD is:
ED

SM DLS SM DHS
D = SM D = d.max{ , } (16)
d d
PT

SM DLS µf ρf
= 4.12Re0.12 W e−0.75 ( )0.54 ( )0.18 (17)
d µa ρa
CE

SM DHS µf ρf
= 0.38Re0.25 W e−0.32 ( )0.37 ( )−0.47 (18)
d µa ρa
Where, d is the injector diameter, Re the Reynolds number, W e the
AC

Webber number and µf and µa the static viscosity of fuel and air. The
superscript LS and HS mean low speed and high speed and represents the
incomplete and complete spray behaviour.

13
 ACCEPTED MANUSCRIPT

4.2. Wall wetting model

The wall wetting and evaporation happens when the fuel jet impinges
on the surface and creates a puddle of fuel on top of the piston surface.
The rate of evaporation of this fuel puddle is determined by the surface
temperature [4]. Curtis et al. [4] defines the evaporation rate of fuel puddle

T
for a port injection engine, and the same equation could be adapted for a

IP
direct injection engine:
Als ∆M F F

CR
ṁwall = Sh[ ρgm Df a Ln(1 + )] (19)
Dp 1 − MF F s
Where Sh is the non-dimensional Sherwood number, Als is the surface
area of the fuel puddle, ρgm is the density of the fuel, Df a the diffusion coef-

US
ficient between the fuel and the air, ∆M F F the difference in mass fraction
of fuel in the vapour phase at the liquid surface and free stream, M F F s is
the mass fraction of fuel in the vapour phase and Dp is the port diameter.
AN
In this work, the port diameter will be substituted by the diameter of the
cylinder.
The Sherwood number as function of Reynolds number, Re, and Schmidt
number, Sc [4] is:
M

Sh = 1 + 0.023Re0.83 Sc0.33 (20)

A MATLAB script was developed for estimating the surface of the puddle
ED

Als . This script assumes the spray jet as a cone and its penetration angle
was estimated with the equations of [65,61]as follows:
√ p
PT

4π 3/6 ρa /ρf
θspray = arctg( ) (21)
3 + 0.28.ln /d
CE

Where ln is the length of the nozzle and d is its diameter.

In the penetration model, the spray has two steps, the first is called pre-
breakup phase and the jet travels freely at constant velocity, when the fluid
achieves the break-up time (tb,k ) the penetration of the spray will change and
AC

it will be in the post-breakup phase:

s
2(Pinj − Pa
x(t)pre = Cd t (22)
ρ

14
 ACCEPTED MANUSCRIPT

s
q
Pinj − Pa 1/4 2(Pinj − Pa )
x(t)post = 2.95( ) d.(t − tb,k ) + Cd tb,k (23)
ρa ρf

T
ρ d 6−k
tb,k = 4.351 p f ( ) (24)

IP
Cd2 ρa (Pinj − Pa ) 5
Where Pinj is the pressure that the fluid is injected, Pa the pressure in the

CR
cylinder d the droplet diameter, k the radial index that it will be assumed as
1 and Cd is the discharge coefficient. This script calculates the surface from
the parametric equation of a cone and the general and implicit equations of
a plane. The equation for cone [66] :

x=
h−u
r. cos θ, y=
h−u
r. sin θ,
US z = u, 0 ≤ u ≤ h, 0 ≤ θ ≤ 2π
AN
h h
(25)
The implicit equation of a plane [66]:

x = x0 + u1 λ + v1 µ
M

y = y0 + u2 λ + v2 µ
ED

z = z0 + u3 λ + v3 µ
PT

−∞ ≤ λ ≤ ∞, −∞ ≤ µ ≤ ∞ (26)
The general equation of a plane [66]:
CE

Ax + By + Cz + D = 0 (27)

The script assumes that the fuel spray is in the shape of perfect cone
AC

and the top of piston is horizontal plane. The depth, angle of spray cone and
other relevant geometric characteristics of the spray cone are estimated using
spray model and simple trigonometric functions. Finally, MATLAB returns
the surface of the cut between the cone and the plane. Also, the script can
return a visual solution of the script to check that the results make sense:

15
 ACCEPTED MANUSCRIPT

T
IP
CR
US
AN
M
ED
PT
CE

Figure 4: Wall wetting area found with MATLAB

AC

16
 ACCEPTED MANUSCRIPT

4.3. Heat transfer model

The Woschni [67] heat transfer model was used to estimate the heat
transfer ratio:

Q = hc A(Tg − Tw ) (28)

T
Where A is the exposed area, Tg and Tw the temperatures of the gas in

IP
the combustion chamber and wall of the cylinder and hc is the convection
coefficient [32] which obey the following equation:

CR
hc = CB m−1 pm wm Tg−0.87 (29)
Where C and m are empirical coefficients which take the values of 0.0035

US
and 0.8 respectively (Heywood 1988). Also, P and Tg are the pressure and
temperature of the gas in the combustion chamber, finally w is the average
gas speed in the cylinder and it is estimated using [32]:
AN
Vd Tr
w = C1 Sp + C2 (P − Pm ) (30)
Pr Vr
Where C1 and C2 are coefficient and their value depends on the stroke ,
M

Vd is the displacement volume, Sp is the mean speed of the piston and Pm is

the motored pressure.
ED

5. SIMULATION
A 1-D simulation of a gasoline direct injection engine developed in GT-
Power was used in this study. The 1-D model is based on predictive com-
PT

bustion modelling approach and is a validated model [68] and the physics
models were coupled with this 1-D model. Detailed description of this model
and the strength of this model could be found in [68]. This coupling was car-
CE

ried out modifying the internal code of the 1-D simulation software through
FORTRAN.
AC

5.0.1. Coupling methodology

The 1-D simulation allows FORTRAN code to modify the internal code
for adding evaporation rate and the heat transfer coefficient for estimating
crank angle resolved quantities and couples with SIMULINK models as shown
in Figure 5. The FORTRAN code is used to read the variables produced by

17
 ACCEPTED MANUSCRIPT

SIMULINK. For the evaporation model, the FORTRAN subroutine will re-
ceive the evaporation rate from SIMULINK and the internal code needs to
know the amount of fuel evaporated in this time step for each component.
Fortunately, the valour of the time step is an internal variable in the FOR-
TRAN code and the operation of multiplying the evaporation rate by the

T
length of time step is trivial. On the other hand, for the heat transfer, the

IP
1-D simulation only allows modification of the convection coefficient. Then,
SIMULINK estimates the coefficient and the FORTRAN subroutine reads
it directly. Finally, the FORTRAN subroutines are compiled in a dynamic

CR
link library (dll) file, this file is called by the 1-D simulation software at the
beginning of the simulation and it is used by the subroutines in the 1-D sim-
ulation at each time step. MATLAB code developed for implementing EOQ
is given in Appendix A.
US
Two different SIMULINK models were developed in order to estimate
the evaporation rate, the convection coefficient and for applying EOQ and
AN
the Lambert W function. The first model calculates the physics phenomena
and the second the EOQ solution with the Lambert W function. These two
models use variables from the 1-D simulation to achieve the solution. These
M

models were compiled in C code to create a dynamic link library (dll) file.
These files are called from the 1-D simulation software for each time step.
The EOQ is referred to the whole cycle and not for each time step, therefore,
ED

the solution is calculated between 50 and 100 degrees before the firing TDC,
since the injection event should have been completely finished and therefore,
all the parameters required for estimating Q* will be available at this stage.
PT

The 1-D simulation was created from a gasoline direct engine with the
specifications shown in Table 1. The angles are referred to the top dead
centre (TDC) in the firing stroke in Table 1.
CE

After building the engine simulation, the 1-D variables are received by
coupling FORTRAN and SIMULINK subroutines through sensors and results
AC

variables (RLT) in the model. Also, this study assumes that the conditions
do not vary between the cylinders, therefore, the evaporation rate and heat
transfer in each cylinder is the same but the signal is delayed by the firing
order.

18
 ACCEPTED MANUSCRIPT

Bore 77mm
Stroke 85.8mm
Compression Ratio 10.5
displacement 1598 cc

T
Rated Power 171 bhp at 5500 rev/min

IP
Rated Torque 240 Nm at 1700-4500 rev/min
fuel injection type direct injection and wall-guided
maximum lift point for intake valve -450 CA

CR
Duration of intake 211 CA
Maximum lift point for exhaust valve 250 CA
Duration of Exhaust 205 CA
Rated Power
Rated Torque US 171 bhp at 5500 rev/min
240 Nm at 1700-4500 rev/min
Table 1: Specification of experimental engine
AN
M
ED
PT
CE
AC

Figure 5: 1-D engine simulation model coupled with Economic Order Quantity (EOQ)
and LambertW function using SIMULINK and FORTRAN

19
 ACCEPTED MANUSCRIPT

Speed(rev/min)/Torque(Nm) 1500 2000 2400 2800

20 91.6 93.1
40 97.8 96.7 92.1 94.8
60 94.0 98.3 93.5 94.1
80 95.9 98.9 95.8 95.2

T
100 94.0 94.5

IP
120 94.4 94.2
Table 2: Correlation coefficients in percentage for the model validated against the experi-

CR
mental in-cylinder pressure at different engine operating conditions.

6. RESULTS AND DISCUSSION

US
After coupling FORTRAN and SIMULINK models with 1-D crank angle
resolved engine simulation model, the results were obtained from the simula-
tion and compared with the data obtained in a test bench in the first stage.
AN
In the second stage, the EOQ results were evaluated.

6.1. Model Validation

To validate the model, the measured cylinder pressure data obtained from
M

the engine was used to compare with the results of the 1-D simulation model.
Cylinder pressure is the most important data (Heywood, 1988) in the engine
and, it enables us to derive significant amount of information in relation
ED

to combustion process. Experimental data from two different engine speeds

(1500 and 2000 rpm) and three loading conditions (40, 60 and 80 Nm) were
used for validating this approach. The model was also tested for two other
PT

engine speeds (2400 and 2800 rpm) and five loading conditions (20, 40, 60,
80 and 120 Nm).
Two different approaches to validate the data were used in the (Ventura
CE

and Samuel 2016) study: the peak pressure correlation and the area correla-
tion. Peak pressure and area under the Pressure-crank angle are chosen for
correlation since they are directly linked to energy release characteristics and
indicated mean effective pressure from combustion process. The correlations
AC

were compared between the data obtained in the laboratory and the data of
the 1-D simulation. In addition, the correlation between each sets of data
was estimated with MATLAB for each engine condition separately as shown
in Figure 6 and Table 2.

20
 ACCEPTED MANUSCRIPT

T
IP
CR
US
Figure 6: Model validation using Maximum Cylinder pressure and Area of pressure Vs
Crank angle trace
AN
The correlation coefficients for the two first approaches are 99.6% and
89.6%. Also, all the correlation coefficients evaluated in each engine condition
separately have values above 91%. This level of correlations is considered
M

sufficient for combustion analysis.

6.2. Results
ED

Finally, the fuel injection strategy is evaluated for each engine operating
condition. The final amount of fuel evaporated considering all the loses
was assessed as shown in Table 3. The Table 3 shows that the least fuel
evaporated is for high speed and low torque (2400 rpm and 40 and 60 Nm)
PT

operating condition. Although the amount of fuel evaporated is related to

the fuel saved in the EOQ function through the parameter, and it does not
hold a linear relationship, it is possible to see that the maximum fuel saved
CE

is for 2400 rpm and 20 Nm with a percentage of evaporation of 97.77%. The

mean fuel saving is 4.7% for all conditions. The Table 3 shows that, when
all the fuel is evaporated the solution of EOQ function is zero. Then, if the
AC

average is estimated only for the conditions that the fuel is not completely
evaporated, this average of fuel saved will be 6.91%. Table 4 shows inputs
and outputs of Lambert W function for each engine condition.

21
 ACCEPTED MANUSCRIPT

Speed Torque mass of fuel (mf) mf evaporated mass evaporated Q mf reduced

(rev/min) (Nm) injected (kg) (kg) % %
1500 40 1.31E-05 1.31E-05 100 0.00E+00 0.00
60 1.84E-05 1.50E-05 99.51 -1.08E-06 -5.83
80 2.40E-05 2.39E-05 81.42 0.00E+00 0.00
2000 40 1.43E-05 1.43E-05 100.00 0 0.00

T
60 1.65E-05 1.56E-05 94.45 -1.82E-06 -11.05
80 1.92E-05 1.54E-05 80.26 -1.04E-06 -5.44

IP
2400 20 1.82E-05 1.77E-05 97.77 -2.82E-06 -15.52
40 6.77E-06 4.18E-06 61.80 -3.28E-07 -4.85
60 7.79E-06 4.10E-06 52.61 -3.18E-07 -4.09
80 1.23E-05 1.03E-05 83.56 -6.92E-07 -5.62

CR
100 1.68E-05 1.67E-05 99.46 0 0.00
120 2.04E-05 2.02E-05 99.34 0 0.00
2800 20 1.89E-05 1.62E-05 85.89 -9.56E-07 -5.06
40 9.93E-06 7.46E-06 75.10 -5.01E-07 -5.04
60
80
100
120
1.14E-05
1.46E-05
2.20E-05
2.20E-05
US
9.02E-06
1.30E-05
2.11E-05
2.20E-05
79.13
88.80
95.82.80
100.00
-5.99E-07
-9.52E-07
-1.92E-06
0
-5.25
-6.52
-8.71
0.00
AN
Table 3: Percentage of fuel saved with EOQ and Lambert W function
M

7. Conclusion
This work has demonstrated that the EOQ and Lambert W function can
be implemented in a 1-D engine simulation model for developing fuel injection
ED

strategy for estimating optimum quantity of the fuel that should be injected
for achieving desired torque. The present work demonstrated that by cou-
pling 1-D engine simulation model with Simulink and Fortran sub-routines
PT

one can estimate the evaporation rate of fuel injected, and heat tranfer rate
in the engine as a function of crank angle in order to estimate the amount of
fuel available in the gas phase for combustion process.
CE

This novel approach based on the principles derived from EOQ for melon
picking has the potential for developing software in loop calibration process
AC

for optimizing fuel injection quantity in gasoline direct injection engines and
other liquid fuel and fuel neutral combustion engines.

The proposed approach could be extended across all type of combustion

devices provided the suitable physical models are identified.Therefore, EOQ
and Lambert W function based calibration process offer promising direction

22
 ACCEPTED MANUSCRIPT

T
IP
Speed(rev/min) Torque(Nm) Z W−1 W−1 +1

CR
1500 40 -0.36788 -1 0
60 -0.36788 -1.00022 -2.23E-04
80 -0.36788 -1 0
2000 40
60
80
US -0.36788
-0.36788
-0.36788
-1
-1.00008
-1.00017
0
-8.47E-05
-1.74E-04
AN
2400 20 -0.36788 -1.00005 -5.14E-05
40 -0.36788 -1.00016 -1.58E-04
60 -0.36788 -1.00019 -1.86E-04
80 -0.36788 -1.00014 -1.37E-04
M

100 -0.36788 -1 0
120 -0.36788 -1 0
2800 20 -0.36788 -1.00014 -1.38E-04
ED

40 -0.36788 -1.00013 -1.34E-04

60 -0.36788 -1.00013 -1.27E-04
80 -0.36788 -1.0001 -1.04E-04
PT

100 -0.36788 -1.00008 -8.14E-05

120 -0.36788 -1 0
Table 4: Inputs and outputs of Lambert W function for different engine operating condi-
CE

tions
AC

23
 ACCEPTED MANUSCRIPT

for improving the fuel economy, and to reduce combustion generated CO2
and harmful combustion generated pollutants.

The current work clearly demonstrated the potential of EOQ and Lambert
W function for improving the fuel economy, improving the combustion control

T
and therefore, reduce the pollutants, however, the implementation using real

IP
hardware and hardware in loop is yet to be done. Therefore, the future work
will implement it in real-hardware and thereby aim to reduce the development
time for meeting the emission regulations significantly.

CR
References

US
[1] Bowler, L.L., 1980. Throttle Body Fuel Injection (Tbi) an Integrated
Engine Control System. SAE.800164. http://dx.doi.org/10.4271/800164.

[2] Kirwan, J. E., Drallmeier, J.A., Coverdill, R. E., Crawford, R. R., Pe-
AN
ters, J. E., 1989. Spray Characteristics of Throttle Body Fuel Injection.
SAE.890318. http://dx.doi.org/10.4271/890318.

[3] Ali Khan,M., Watson, H., Baker, P., Liew, G.,Johnston, D.,2006. SI
M

Engine Lean-Limit Extension through LPG Throttle-Body Injection for

Low CO2 and NOx..SAE. 2006-01-0495. http://dx.doi.org/10.4271/2006-
01-0495.
ED

[4] Curtis, E.W., Aquino, C.F., Trumpy, D.K., Davis, G.C., 1996. A
New Port and Cylinder Wall Wetting Model to Predict Transient
Air/ Fuel Excursions in a Port Fuel Injected Engine. SAE.961186.
PT

http://dx.doi.org/10.4271/961186.

[5] Z.Han, Z. Xu, N. Trigui. Spray/Wall Interaction Models for Mul-

CE

tidimensional Engine Simulation. International Journal of Engine

Research 1, no. 1 (2000/02/01): 127-146. Accessed 2017/03/09.
http://dx.doi.org/10.1243/1468087001545308.
AC

[6] Senda, J., Ohnishi, M., Takahashi,T., Fujimoto, H., Utsunomiya, A.,
Wakatabe, M.,1999. Measurement and Modeling on Wall Wetted Fuel Film
Profile and Mixture Preparation in Intake Port of Si Engine. SAE.1999-
01-0798. http://dx.doi.org/10.4271/1999-01-0798.

24
 ACCEPTED MANUSCRIPT

[7] Dunne, W.J., Brace, C.J., Stodart, A., 2007. Automated Cal-
ibration of an Analytical Wall-Wetting Model SAE.2007-01-0018.
http://dx.doi.org/10.4271/2007-01-0018.

[8] Yu, J., Abe, M., Sukegawa, J., 2017. A Numerical Method to Simulate

T
Intake-Port Fuel Distribution in Pfi Engine and Its Application. SAE.2017-
01-0565. http://doi.org/10.4271/2017-01-0565.

IP
[9] Ye, Z., Mohamadian. H., 2013. Model Predictive Control on Wall Wetting

CR
Effect Using Markov Chain Monte Carlo. IEEE Latin-America Conference
on Communications.

[10] J. Zhang, T. Shen, G. Xu, J. Kako, Wall-Wetting Model Based Method

US
for Air-Fuel Ratio Transient Control in Gasoline Engines with Dual In-
jection System. International Journal of Automotive Technology 14, no. 6
(2013): 867-873. http://dx.doi.org/10.1007/s12239-013-0095-y.
AN
[11] B.Wang, S.Mosbach, S.Schmutzhard, S.Shuai, Y.Huang, M.Kraft, Mod-
elling Soot Formation from Wall Films in a Gasoline Direct Injection En-
gine Using a Detailed Population Balance Model. Applied Energy 163
(2/1/): 154-166 (2016) http://dx.doi.org/10.1016/j.apenergy.2015.11.011.
M

[12] Malaguti, S., Cantore, G., Fontanesi, S., Lupi, R., Rosetti,
A., 2009. CFD Investigation of Wall Wetting in a Gdi En-
ED

gine under Low Temperature Cranking Operations. SAE.2009-01-0704.

http://dx.doi.org/10.4271/2009-01-0704.

[13] Luijten, B., Adomeit, P., Brunn, A., Somers, B.,2013. Experimental
PT

Investigation of in-Cylinder Wall Wetting in Gdi Engines Using a Shad-

owgraphy Method. SAE. 2013-01-1604. http://dx.doi.org/10.4271/2013-
01-1604.
CE

[14] Ates, M., Matthews, R.D., Hall, M.J., 2017. A Full-Cycle Multi-Zone
Quasi-Dimensional Direct Injection Diesel Engine Model Based on a Con-
AC

ceptual Model Developed from Imaging Experiments. SAE.2017-01-0537.

http://doi.org/10.4271/2017-01-0537.

[15] Chen, L., et al., 2011. Spray Imaging, Mixture Preparation and Par-
ticulate Matter Emissions Using a GDI Engine Fuelled with Stoichiomet-
ric Gasoline/Ethanol Blends. In Internal Combustion Engines: Improving

25
 ACCEPTED MANUSCRIPT

Performance, Fuel Economy and Emission IMecHE, ISBN 978-0-85709-

20502. 43-52: Woodhead Publishing.

[16] Zhao, F.Q., Lai, M.C., Harrington, D.L., 1997. A Review

of Mixture Preparation and Combustion Control Strategies for

T
Spark-Ignited Direct-Injection Gasoline Engines. SAE.970627.
http://dx.doi.org/10.4271/970627.

IP
[17] M. Castagn, J.P. Dumas, S. Henriot, Ph. Pierre, Use of Advanced Tools

CR
for the Analysis of Gasoline Direct Injection Engines. Oil & Gas Science
and Technology - Rev. IFP 58, no. 1 (2003) 79-100.

[18] J. F. Le Coz, J. Cherel, S. Le Mirronet, Fuel/Air Mixing Process and

US
Combustion in an Optical Direct-Injection Engine. Oil & Gas Science and
Technology - Rev. IFP 58 (2003) no. 1 : 63-78.

[19] Pischinger, S., Schernus, C., Ltkemeyer, G., Theuerkauf, H.J.,

AN
Winsel, T., Ayeb, M.,2004. Investigation of Predictive Models
for Application in Engine Cold-Start Behavior. SAE.2004-01-0994.
http://dx.doi.org/10.4271/2004-01-0994.
M

[20] Indranil, B., 2011. Development of Dynamic Constraint Models

for a Model Based Transient Calibration Process. SAE.2011-01-0691.
http://dx.doi.org/10.4271/2011-01-0691.
ED

[21] Kitamura, Y., Sekikawa, A., Tokoro, M., Tomatsu, N., Uchimoto, T.,
Yuasa, H., Kato, A.,2011. Investigation About Predictive Accuracy of Em-
pirical Engine Models Using Design of Experiments. SAE.2011-01-4271.
PT

http://dx.doi.org/10.4271/2011-01-1848.

[22] Wilhelm, C., Winsel, T., Ayeb, M., Theuerkauf, H.J., Brandt, S.,
CE

Busche, E., Longo, C., Knoll, G.D., 2007. A Combined Physical / Neu-
ral Approach for Real-Time Models of Losses in Combustion Engines.
SAE.2007-01-1345. http://dx.doi.org/10.4271/2007-01-1345.
AC

[23] E.Corti, N.Cavina, A.Cerofolini, C.Forte, G.Mancini, D.Moro,

F.Ponti, V.Ravaglioli, Transient Spark Advance Calibration Ap-
proach. Energy Procedia 45 (2014/01/01/): 967-976 (2014)
https://doi.org/10.1016/j.egypro.2014.01.102.

26
 ACCEPTED MANUSCRIPT

[24] Rask, E., Mark, S.,2004. Simulation-Based Engine Calibra-

tion: Tools, Techniques, and Applications. SAE.2004-01-1264.
http://dx.doi.org/10.4271/2004-01-1264.

[25] R.F.Turkson, F.Yan, M.K.A.Ali, B.Liu, J.Hu., Modeling and

T
Multi-Objective Optimization of Engine Performance and Hydrocar-
bon Emissions Via the Use of a Computer Aided Engineering Code

IP
and the Nsga-Ii Genetic Algorithm. Sustainability 8, no. 1.(2016)
http://dx.doi.org/10.3390/su8010072.

CR
[26] T.Luigi, A.Roberto, N.Fabio, A 1D/3D Methodology for the
Prediction and Calibration of a High Performance Motorcycle
SI Engine. Energy Procedia 82 (2015/12/01/): 936-943 (2015)

US
https://doi.org/10.1016/j.egypro.2015.11.842.

[27] Torre, D., Montenegro, A.G., A. Onorati.A., 2017. Coupled 1d-Quasi3d

AN
Fluid Dynamic Models for the Simulation of Ic Engine Intake and Exhaust
Systems. In 17. Internationales Stuttgarter Symposium: Automobil- Und
Motorentechnik, edited by Michael Bargende, Hans-Christian Reuss, and
Jochen Wiedemann, 1461-1476. Wiesbaden: Springer Fachmedien Wies-
M

baden.

[28] Omekanda, S.O., Rahman, R., Lott, E.M., Rahman, S. S., Hornback,
D.E., 2017. Optimal Parameter Calibration for Physics Based Multi-Mass
ED

Engine Model. SAE. 2017-01-0214. http://doi.org/10.4271/2017-01-0214.

[29] Lumsden, G., Browett, C., Taylor, J., Kennedy, G., 2004. Mapping
PT

Complex Engines. SAE.2004-01-0038. http://dx.doi.org/10.4271/2004-01-

0038.

[30] R.Finesso, O.Marello, D.Misul, E.Spessa, M.Violante, Y.Yang, G.Hardy,

CE

C.Maier, Development and Assessment of Pressure-Based and Model-

Based Techniques for the MFB50 Control of a Euro VI 3.0L Diesel Engine,
SAE Int. J. Engines 10(4):1538-1555 (2017) https://doi.org/10.4271/2017-
AC

01-0794.

[31] Wurzenberger, J., Bartsch, P., and Katrasnik, T., 2010. Crank-
Angle resolved Real time Capable Engine and Vehicle Simulation-
Fuel consumption and Driving Performance. SAE.2010-01-0784.
http://dx.doi.org/10.4271/2010-01-0784.

27
 ACCEPTED MANUSCRIPT

[32] J.Heywood, Internal Combustion Engine Fundamentals.McGraw-Hill,

New York, 1988.

[33] J.I.Ghojel, Review of the development and applications of the

Wiebe function: A tribute to the contribution of Ivan Wiebe

T
to engine research,International Journal of Engine Research, Vol
11, Issue 4, pp. 297 - 312, First Published June 23 (2010),

IP
https://doi.org/10.1243/14680874JER06510

CR
[34] Kumar, P. and Makki, I., 2017. Model Based Control of a Three-Way
Catalytic Converter Based on the Oxygen Storage Level of the Catalyst.
SAE.2017-01-0960.

US
[35] Lee, H., Lee, B., Kim, S., Kim, N., Rousseau, A., 2017.
Model-Based Fuel Economy Technology Assessment. SAE.2017-01-0532.
http://doi.org/10.4271/2017-01-0532.
AN
[34] Watanabe, S., Miyata, Y., Ogata, Y., and Ivosic, V., 2017. Applica-
tion of Model-Based Development to Engine Restart Vibration after Idling
Stop. SAE .2017-01-1053. http://dx.doi.org/10.4271/2017-01-1053.
M

[37] C.Friedrich, M.Auer, G.Stiesch, Model Based Calibration Techniques for

Medium Speed Engine Optimization: Investigations on Common Modeling
Approaches for Modeling of Selected Steady State Engine Outputs, SAE
ED

Int. J. Engines 9(4):1989-1998 (2016), https://doi.org/10.4271/2016-01-

2156.

[38] Poetsch, C. and Katrasnik, T., 2016. Crank-Angle Resolved Model-

PT

ing of Fuel Injection, Combustion and Emission Formation for Engine

Optimization and Calibration on Real-Time Systems. SAE.2016-01-0558
http://dx.doi.org/10.4271/2016-01-0558.
CE

[39] Schlosser, A., Kinoo, B., Salber, W., Werner, S., Ademes, N., 2006.
Accelerated Powertrain Development through Model Based Calibration.
AC

SAE.2006-01-0858. http://dx.doi.org/10.4271/2006-01-0858.

[40] Stuhler, H., Kruse, T., Stuber, A., Gschweitl, K., Piock,
W., Pfluegl, H., Lick, P.,2002. Automated Model-Based Gdi En-
gine Calibration Adaptive Online Doe Approach. SAE .2002-01-0708.
http://dx.doi.org/10.4271/2002-01-0708.

28
 ACCEPTED MANUSCRIPT

[41] Jiang, S., Nutter, D., Gullitti, A., 2012. Implementation of

Model-Based Calibration for a Gasoline Engine. SAE.2012-01-0722.
http://dx.doi.org/10.4271/2012-01-0722.

[42] Suzuki, K., Nemoto, M., Machida, K., 2006. Model-Based Cali-

T
bration Process for Producing Optimal Spark Advance in a Gaso-
line Engine Equipped with a Variable Valve Train. SAE.2006-01-3235.

IP
http://dx.doi.org/10.4271/2006-01-3235.

CR
[43] Beattie, T., Osborne, R., Graupner, W., 2013. Engine Test
Data Quality Requirements for Model Based Calibration: A Test-
ing and Development Efficiency Opportunity. SAE .2013-01-0351.
http://dx.doi.org/10.4271/2013-01-0351.

US
[44] Onder, C., Geering, H., 1995. Model-Based Engine Calibration for Best
Fuel Efficiency. SAE.950983. http://dx.doi.org/10.4271/950983.
AN
[45] Wurzenberger, J. and Poetsch, C., 2015. Plant Modeling for Closed Loop
Combustion Control - a Thermodynamic Consistent and Real-Time Ca-
pable Approach. SAE. 2015-01-1247. http://dx.doi.org/10.4271/2015-01-
1247.
M

[46] A.S.Abbas, S.Saraswati, Cycle-by-Cycle Estimation of Cylinder Pres-

sure and Indicated Torque Waveform Using Crankshaft Speed Fluctua-
ED

tions. Transactions of the Institute of Measurement and Control 37,(2015)

no. 6: 813-825.

[47] Zhu, G., Daniels, C., and Winkelman, J., 2003. MBT Timing Detection
PT

and Its Closed-Loop Control Using in-Cylinder Pressure Signal. SAE.2003-

01-3266. http://dx.doi.org/10.4271/2003-01-3266.
CE

[48] Itshak,G., Powell, J.D., 1981. Optimal Closed-Loop Spark Control of an

Automotive Engine. SAE.810058. http://dx.doi.org/10.4271/810058.

[49] Kampelmhler, F. T., Paulitsch, R., Gschweitl, K., 1993. Automatic

AC

Ecu-Calibration - an Alternative to Conventional Methods. SAE .930395.

http://dx.doi.org/10.4271/930395.

[50] Rivard, J. G. Closed-Loop Electronic Fuel Injection Con-

trol of the Internal Combustion Engine. SAE.730005.(1973)
http://dx.doi.org/10.4271/730005.

29
 ACCEPTED MANUSCRIPT

[51] Schnorbus, T., Pischinger, S., Krfer, T., Lamping, M., Tomazic,
D., Tatur, M., 2008. Diesel Combustion Control with Closed-
Loop Control of the Injection Strategy. SAE.2008-01-0651.
http://dx.doi.org/10.4271/2008-01-0651.

T
[52] Franz, J., Schwarz, F., Guenthner, M., Reissing, J., Mueller, A.,
Donn, C., 2009. Closed Loop Control of an Hcci Multi-Cylinder

IP
Engine and Corresponding Adaptation Strategie. SAE.2009-24-0079.
http://dx.doi.org/10.4271/2009-24-0079.

CR
[53] Serizawa, K.,Ueda, D., Mikami, N., Tomida, Y., Weber, J., 2017.
Realizing Robust Combustion with High Response Diesel Injector with
Controlled Diffusive Spray Nozzle and Closed Loop Injection Control.

US
SAE.2017-01-0845. http://doi.org/10.4271/2017-01-0845.

[54] John, C., Rachel, T., 1975. Closed-Loop Electronic Fuel

AN
and Air Control of Internal Combustion Engines. SAE.750369.
http://dx.doi.org/10.4271/750369.

[55] Yinbo, C., Li, L., 2012. A Novel Closed Loop Control Based on Ion-
ization Current in Combustion Cycle at Cold Start in a GDI Engine.
M

SAE.2012-01-1339. http://dx.doi.org/10.4271/2012-01-1339.

[56] ika, Z., Valek, M., Florin, M., Macek, J., Polek, M.,2008. Multilevel Pre-
ED

dictive Models of Ic Engine for Model Predictive Control Implementation.

SAE.2008-01-0209. http://dx.doi.org/10.4271/2008-01-0209.

[57] V.Robert, S.Samuel, Optimization of Fuel Injection in GDI En-

PT

gine Using Economic Order Quantity and Lambert W Function.

Applied Thermal Engineering 101 (2016/05/25/): 112-120.(2016)
http://dx.doi.org/10.1016/j.applthermaleng.2016.02.024.
CE

[58] D.Stephen, R.Warburton, On the Lambert W Function: Economic Or-

der Quantity Applications and Pedagogical Considerations.International
AC

Journal of Production Economics Volume 140, Issue 2, (2012), 756-764.

[59] B.Joseph, G.Scudder, Supply Chain Strategies for Perishable Prod-

ucts: The Case of Fresh Produce. Production and Operations Man-
agement 18, (2009). no. 2: 129-137. http://dx.doi.org/10.1111/j.1937-
5956.2009.01016.x.

30
 ACCEPTED MANUSCRIPT

[60] Takuo, Y., Nishida, K., Hiroyuki Hiroyasu, H.,(1993). Approach to Low
Nox and Smoke Emission Engines by Using Phenomenological Simula-
tion.SAE.930612. http://dx.doi.org/10.4271/930612, 1993.

[61] A.Lefebvre, Atomization and Sprays. CRC Press, 1988.

T
[62] Robert, M., Heywood, J.B.,1999. Evaporation of in-Cylinder Liquid Fuel

IP
Droplets in an Si Engine: A Diagnostic-Based Modeling Study. SAE.1999-
01-0567. http://dx.doi.org/10.4271/1999-01-0567.

CR
[63] C.Antoine, Tensions Des Vapeurs; Nouvelle Relation Entre Les Ten-
sions Et Les Tempratures. Comptes Rendus des Sances de l’Acadmie des
Sciences (1888) 107: 681-684, 778-780, 836-837.

US
[64] Hiroyuki, H., Arai, M., Tabata, M., 1989. Empirical Equations
for the Sauter Mean Diameter of a Diesel Spray. SAE.890464.
http://dx.doi.org/10.4271/890464.
AN
[65] Dohoy, J., Assanis. D.N., 2001. Multi-Zone Di Diesel Spray Combustion
Model for Cycle Simulation Studies of Engine Performance and Emissions.
SAE.2001-01-1246. http://dx.doi.org/10.4271/2001-01-1246.
M

[66] I.N.Bronshtein , K.A.Semendyayev, G.Musiol , H.Mhlig, Handbook of

Mathematics. 5th ed.: Springer, 2007.
ED

[67] Woschni, G., 1967. A Universally Applicable Equation for the Instan-
taneous Heat Transfer Coefficient in the Internal Combustion Engine.
SAE.670931. http://dx.doi.org/10.4271/670931.
PT

[68] Divakera, Arjun, D., Samuel, S.,2011. Numerical Simulation of Adaptive

Combustion Control for Fuel-Neutral Smart Engines. SAE.2011-01-0848.
CE

http://dx.doi.org/10.4271/2011-01-0848.

[69] Spessa, E., D’Ambrosio, S., Iemmolo, D., Mancarella, A., Vi-
tolo, R., Hardy, G., 2017. Steady-State and Transient Operations
AC

of a Euro Vi 3.0l Hd Diesel Engine with Innovative Model-Based

and Pressure-Based Combustion Control Techniques. SAE.2017-01-0695.
http://dx.doi.org/10.4271/2017-01-0695.

31
 ACCEPTED MANUSCRIPT

[70] M. Y. Jaber, R. Y. Nuwayhid, M. A. Rosen, Price-driven eco-

nomic order systems from a thermodynamic point of view, Interna-
tional Journal of Production Research, 42:24, 5167-5184,(2007) DOI:
10.1080/00207540412331281971.

T
[71] H.Jawad, Y.M.Jaber, M.Bonney, The Economic Order Quantity model
revisited: an Extended Exergy Accounting approach. Journal of cleaner

IP
production, 105, 64-73, (2015). doi: 10.1016/j.jclepro.2014.06.079.

CR
[72] H.Jawad, M.Y.Jaber, M.Bonney, M.A.Rosen, Deriving an exergetic eco-
nomic production quantity model for better sustainability. Applied Math-
ematical Modelling, 40(11-12),(2016) 6026-6039.

US
[73] F.M.Harris, How Many Parts to Make at Once, Factory, The Magazine
of Management,10:2,(1913) 135136, 152. Reprinted in Operations Research
38:6 (1990), 947950.
AN
[74] L.B.Schwarz, The economic order-quantity (EOQ) model, chapter
8 from Building Intuition: Insights From Basic Operations Man-
agement Models and Principles,(2008) doi:10.1007/978-0-387-73699-0-8,
ISSN:0884-8289.
M
ED

8. Appendix A
PT
CE
AC

Figure A1 Overview of Flow of information and iterative process

32
 ACCEPTED MANUSCRIPT

Matlab Code

1 f u n c t i o n [ evap , a n g l e s p r a y , x s p r a y , m wall ]= c a l c u l a t e s (
t h e t a , P c y l i n d e r , T c y l i n d e r , mu f , mu a , k a , m traped ,
fuel rate injected , P rail , v injector ,

T
a v e r a g e g a s v e l o c i t y , dt , s u r f a c e )

IP
2 g l o b a l time
3 g l o b a l m spray acumulated
g l o b a l m wall acumulated

CR
4

5 g l o b a l Pa
6 global T drop previous
7 global flag
8

10
global d previous
global v previous
global m fuel traped
g l o b a l m spray
US
AN
11

12 global time injected

13 global t bk
14 % global m wall traped
M

15 global D droplet
16 global surface saved
17
ED

18 %d e f i n i t i o n o f n u m e r i c a l v a r i a b l e s ( r e q u i r e d by
MATLAB t o c o m p i l e t h e s c r i p t )
19 d v =0; d=0; v=0; m wall =0; d T drop =0; C drag =0;
Re=0; We=0; L f s =0; Bt=0; incr MFF =0; MFF=0;
PT

P f s =0; Yfs =0; Bm=0; C drag =0; Sc =0; Sh=0;

SMD ls =0; SMD hs=0; T drop =0; d d =0; x s p r a y =0;
a n g l e s p r a y =0; t i m e f r o m i n j e c t i o n =0;
CE

m s p r a y t o t a l =0; m suspended =0; N d r o p l e t s =0;

v o l s u s p e n d e d =0; P f s f r e e =0; Y f s f r e e =0; r =0;
Sc =0; A l s =0; m a i r t r a p e d =0;
AC

20

21 %c o n s t a n t s o f th e model
22 bore =77/1000; %b ore o f t h e c y l i n d e r
23 d i a m t e r i n j e c t o r= 0 . 3 / 1 0 0 0 ; %d i a m e t e r o f th e
injector
24 Cd= 0 . 7 ; %N o z z l e D i s c h a r g e

33
 ACCEPTED MANUSCRIPT

C o e f f i c i e n t of the i n j e c t o r
25 T i n j e c t o r =350; %I n j e c t e d F l u i d
Temperature
26 Ma= 2 9 . 9 7 ; %Molar we ith o f a i r
27 Mf= 1 1 4 . 2 2 8 5 ; %Molar weith o f f u e l

T
http : / / webbook . n i s t . gov / c g i / i n c h i ?ID=C111659&

IP
Mask=4#Thermo−Phase
28 Rg air =287.058; %Constant ga s o f a i r
Cp a =1024; %S p e c i f i c h e at o f a i r

CR
29

at constante p r e s s u r e
30 Cp f =2341; %S p e c i f i c h e at o f f u e l
a t c o n s t a n t e p r e s s u r e h tt p : / / webbook . n i s t . gov /

31

32 k=1;
o f break up
US
c g i / i n c h i ?ID=C111659&Mask=2#Thermo−Condensed

%R a d i a l i n d e x f o r time
AN
33 sigma = 0 . 0 2 1 6 2 ; %S u r f a c e t e n s i o n o f
ht tp : / /www. s u r f a c e −t e n s i o n . de /
34 Tcr = 5 6 8 . 9 ; %C r i t i a l t e m p e r a t u r e o f
M

o c t a n e ( f u e l ) ht tp : / / webbook . n i s t . gov / c g i / i n c h i
?ID=C111659&Mask=4#Thermo−Phase
35 Tbn= 3 9 8 . 7 ; %B o i l t e m e p e r a t u r e o f
ED

Octane ( f u e ) htt p : / / webbook . n i s t . gov / c g i /

i n c h i ?ID=C111659&Mask=4#Thermo−Phase
36 LTBN=350000; %Heat o f V a p o r i z a t i o n
a t 298K o f f u e l GT−POWER
PT

37 D fa =7.353E−7; %D i f f u s i o n c o e f i c i e n t
between o c t a n e and a i r htt p : / /www. j o c e t . o r g /
p a p e r s /012− J30011 . pdf
CE

38 r h o f =750; %D e n s i t y o f f u e l
39

40 r h o a=P c y l i n d e r / ( T c y l i n d e r ∗ R g a i r ) ; %d e n s i t y o f
AC

air
41 Ap=( p i /4 ) ∗ bore ˆ 2 ; %Area o f p i s t o n
42 % = 0.3∗ P r a i l ; %P r e s s u r e l o s t i n t h e
injection
43 l n=d i a m t e r i n j e c t o r / 1 0 ; %Lenght o f n o z z l e
44

34
 ACCEPTED MANUSCRIPT

45 a n g l e o f i n j e c t i o n =360;
46

47 %Antonie c o e f f i c i e n t s ht tp : / / webbook . n i s t . gov / c g i /

i n c h i ?ID=C111659&Mask=4#Thermo−Phase
48 i f T c y l i n d e r <297.1

T
49 A = 5.2012;

IP
50 B =1936.281 ;
51 C = −20.143 ;
else

CR
52

53 A = 4.04867;
54 B =1355.126 ;
55 C = −63.633 ;
56

57

58
end
US
%i n i c i a l i z a t i n g th e g l o b a l v a r i a b l e s i n each c y c l e
i f 20< t h e t a && t h e t a < 60 %we use a ra nge o f
AN
59

a n g l e s t o a v o i d t h a t a b i g t i m e s t e p jumps th e
iniciaization
60 f l a g =0; %f l a g used f o r s a v e th e
M

v a r i a b l e s at s t a r t of i n j e c t i o n
61 m f u e l t r a p e d =0; %f u e l acumulated i n t h e
cyinder
ED

62 m spray acumulated =0; %Acumulated f u e l

e v a p o r a t e d by th e s p r a y j e t
63 m wall acumulated =0; %Acumulated f u e l
e v a p o r a t e d by th e w a l l w e t t i n g
PT

64 %m w a l l t r a p e d =0; %f u e l acumulated i n th e
cyinder
D d r o p l e t=1E−3; %Diameter o f t he
CE

65

droplet
66 s u r f a c e s a v e d =0; %S u r f a c e o f w a l l
w e t t i n g i n th e p r e v i o u s t i m e s t e p ( i s used t o
AC

be s u r e t h a t t h e s u r f a c e o n l y w i l l i n c r e a s e
)
67 end
68

69 %A l s i s th e w a l l w e t t i n g s u r f a c e
70 %This p a r t a v o i d t h a t t h e w a l l w e t t i n g s u r f a c e

35
 ACCEPTED MANUSCRIPT

d e c r e a s i n g i n th e same c y c l e , we a r e assuming
t h a t once t h e s p r a y h i t s t h e p i s t o n , th e puddle
o n l y w i l l be a b l e t o i n c r e a s e
71 i f s u r f a c e s a v e d >s u r f a c e
72 A l s= s u r f a c e s a v e d ;

T
73 else

IP
74 A l s= s u r f a c e ;
75 s u r f a c e s a v e d=s u r f a c e ;
end

CR
76

77

78 %acumuated v a r i a b l e s
79 time=time+dt ; %time o f s i m u l a t i o n
80

81
US
m f u e l t r a p e d=f u e l r a t e i n j e c t e d ∗ dt + m f u e l t r a p e d
; %amount o f f u e l i n t he c y l i n d e r
m a i r t r a p e d=m traped−m f u e l t r a p e d ;
%amount o f a i r i n th e c y l i n d e r
AN
82 e v a p o r a t e d m a s s a c u m u l a t e d=m spray acumulated+
m wall acumulated ; %amount o f f u e l e v a p o r a t e d
83
M

84 i f m f u e l t r a p e d >0 %i f we a r e i n t h e e v a p o r a t i o n
phase , we do not want t o e s t i m a t e a n y t h i n g i s
t he exahust , t h i s s t a r t when t h e f u e l i s
ED

injected
85

86 %v a r i a b e s marked i n t h e i n j e c t i o n p o i n t
i f f l a g==0
PT

87

88 Pa=P c y l i n d e r ; %P r e s s u r e a t
injection
d= D d r o p l e t ; %Diameter o f t h e
CE

89

droplet
90 d p r e v i o u s=d ; %Diameter o f th e
d r o p l e t in the previous timestep
AC

91 v=v i n j e c t o r ; %Speed o f t h e
injection
92 v p r e v i o u s=v ; %Speed o f t h e
i n j e c t i o n i n th e p r e v i o u s t i m e s t e p
93 T drop=T i n j e c t o r ; %Temperature o f t h e
droplet

36
 ACCEPTED MANUSCRIPT

94 T d r o p p r e v i o u s=T drop ; %Temperature o f th e

droplet
95 t i m e i n j e c t e d=time ; %Time when t h e
s p r a y a r e i n j e c t e d , a f t e r we w i l l
c o r r e e c t i t with t i m e f r o m i n j e c t i o n t o

T
s t i m a t e t he t o t a l time from th e

IP
injection
96 t b k = 4 . 2 5 1 ∗ ( ( r h o f ∗ d i a m t e r i n j e c t o r ) / (Cdˆ2∗
s q r t ( r h o a ∗ abs ( P r a i l −P c y l i n d e r ) ) ) )

CR
∗((6 −k ) /5 ) ; %time f o r break−up
97 f l a g =1;
98 end
99

100
validated !
US
%Very h igh Re numbers i n combustion , but th e model i s

Re=r h o a ∗ a v e r a g e g a s v e l o c i t y ∗ bo re /mu a ;
%Re=( f l o w r e ∗ bo re ) / (Ap∗mu a ) ; %
AN
101

Reynols number
102 C drag =(24/Re ) ∗(1+((Re ˆ ( 2 / 3 ) /6 ) ) ) ;
%Drag c o e f f i c i e n t o f a
M

droplet
103 L f s=LTBN∗ ( ( ( Tcr−T d r o p p r e v i o u s ) / ( Tcr−Tbn) )
ˆ −0.38) ; %Latent h e at o f e v a p o r a t i o n
ED

104

105 %c o n s i d e r e r t o change a l l Pa by P c y l i n d e r
106 P f s f r e e =1E5 ∗10ˆ(A−(B/ ( T c y l i n d e r+C) ) ) ; %
P f s o f t he d r o p l e t f u e l e v a p o r a t e d
PT

107 Y f s f r e e =1/(1+((Pa/ P f s f r e e ) ∗(Ma/Mf) ) ) ; %

Yfs o f t he d r o p l e t f u e l e v a p o r a t e d
P f s =1E5 ∗10ˆ(A−(B/ ( T d r o p p r e v i o u s+C) ) ) ; %
CE

108

P f s o f t he d r o p l e t f u e l i n t he l i q u i d s t a t e
109 Yfs =1/(1+((Pa/ P f s ) ∗(Ma/Mf) ) ) ; %
Yfs o f t he d r o p l e t f u e l i n t h e l i q u i d s t a t e
AC

110

111 incr MFF=abs ( Yfs−Y f s f r e e ) ;

112 MFF=Yfs ;
113

114 %Bt and Bm we have a l o o p h e r e b e c a u s e th e

change o f t h e d r o p l e t d e p e r a t u r e depens on

37
 ACCEPTED MANUSCRIPT

Bt and Bm and , a t t h e same time Bt and Bm

depends on t h e t e m p e r a t u r e o f t h e d r o p l e t
115 i f Yfs==1
116 Bm=1;
117 else

T
118 Bm=Yfs /(1− Yfs ) ;

IP
119 end
120

Bt=Cp a ∗( T c y l i n d e r −T d r o p p r e v i o u s ) / L f s ;

CR
121

%Mass t r a n s f e r number o f
temperature
122

123

124
%D r o p l e t Temperature
US
%we have a l o o p h e r e b e c a u s e th e change o f t h e
d r o p l e t Temperature depens on Bt and Bm and ,
a t t he same time Bt and Bm depends on t h e
AN
t e m p e r a t u r e o f th e d r o p l e t
125 i f m spray acumulated >0
126 d T drop=m spray ∗( L f s / ( Cp f ∗
M

m spray acumulated ) ) ∗ ( ( Bt/Bm) −1) ; %

D i f f e r e n c i a l of temperature
127 %c o n t r o l t h a t t h e d i f f e r e n c i a l i s not ve ry
ED

b i g due t o i n c r r e c t t i m e s t e p s
128 i f d T drop >10000 | | d T drop <−10000
129 T drop =T d r o p p r e v i o u s ;
else
PT

130

131 T drop=T d r o p p r e v i o u s −(d T drop ∗ dt ) ;

132 end
%C o n t r o l t h a t t h e t e m p e r a t u r e o f t h e
CE

133

d r o p l e t i s not l e s s than 0
134 i f T drop <0
135 T drop=T i n j e c t o r ;
AC

136 end
137 T d r o p p r e v i o u s=T drop ;
138 end
139

140 %D r o p l e t Diameter
141 d d =(4∗ k a ∗ l o g (1+Bm) ) / ( r h o f ∗Cp a∗ d p r e v i o u s ) ;

38
 ACCEPTED MANUSCRIPT

%d i f f e r e n c i a l o f d i a m e t e r
142 d=d p r e v i o u s −d d ∗ dt ;
143 i f d>0 %no d r o p l e t n e g a t i v e s d i a m e t e r o r z e r o
144 d p r e v i o u s=d ;
145 end

T
146

IP
147 %D r o p l e t V e l o c i t y
148 i f d p r e v i o u s >0 %no d r o p l e t n e g a t i v e s d i a m e t e r
or zero

CR
149 d v =((3∗ r h o a ) /(4∗ r h o f ) ) ∗( C drag ) ∗ ( ( ( mu f−
mu a ) ˆ 2) / ( d p r e v i o u s ) ) ;
150 v=v p r e v i o u s −(d v ∗ dt ) ;
151

152

153
v p r e v i o u s=v ;
i f v<0
v=0;
v p r e v i o u s=v ;
US
AN
154

155 end
156 end
157
M

158 %D i m e n s i o n l e s s numbers
159 Sc = mu f / ( r h o f ∗ D fa ) ;
160 We =r h o a ∗( v ˆ 2 ) ∗ D d r o p l e t / sigma ;
ED

161 Sh =1+(0.023∗(Re ˆ ( 0 . 8 3 ) ) ∗( Sc ˆ ( 0 . 3 3 ) ) ) ;
162 %S a u t e r mean d i a m e t e r
163 SMD ls =4.12∗(Re ˆ 0 . 1 2 ) ∗(Weˆ −0.75) ∗ ( ( mu f /mu a )
ˆ0.54) ∗(( rho f / rho a ) ˆ0.18) ;
PT

164 SMD hs=0.38∗(Re ˆ 0 . 2 5 ) ∗(Weˆ −0.32) ∗ ( ( mu f /mu a )

ˆ 0 . 3 7 ) ∗ ( ( r h o f / r h o a ) ˆ −0.47) ;
D d r o p l e t= d i a m t e r i n j e c t o r ∗max( SMD ls , SMD hs )
CE

165

;
166

167 %Number o f D r o p l e t s
AC

168 i f m wall acumulated >0

169 m suspended=m f u e l t r a p e d ∗ 0 . 9 5 ; %C o r r e c t i o n
i f we have w a l l w e t t i n g
170 else
171 m suspended=m f u e l t r a p e d ;
172 end

39
 ACCEPTED MANUSCRIPT

173 v o l s u s p e n d e d=m suspended / r h o f ;

174 N d r o p l e t s=v o l s u s p e n d e d / ( ( 4 ∗ p i / 3 ) ∗( D d r o p l e t
/2) ˆ3) ;
175

176 i f Bm<0 %C o n t r o l t o a v o i d n e g a t i v e l o g

T
177 Bm=0;

IP
178 end
179

%e v a p o r a t i o n r a t e s

CR
180

181 m spray =2∗3.1415∗ D d r o p l e t ∗( k a / Cp a ) ∗ l o g (1+Bm)

; %i t i s n e c e s s a r y f o r s t i m a t e t he d r o p l e t
temperature , f o r t h i s r e a s o n i t i s a g l o b a l

182

183
variable
US
m s p r a y t o t a l=m spray ∗ N d r o p l e t s ;
m wall=Sh ∗( r h o f ∗ D fa ∗( A l s / bo re ) ∗ l o g (1+(
incr MFF/(1−MFF) ) ) ) ;
AN
184

185 evap=m s p r a y t o t a l+m wall ; %F i n a l e v a p o r a t i o n

rate
M

186

187 %Acumulated e v a p o r a t i o n s
188 m spray acumulated=m spray acumulated+
ED

m s p r a y t o t a l ∗ dt ;
189 m wall acumulated=m wall acumulated+m wall ∗ dt ;
190

%I f a l l f u e l were e v a o r a t e d , we can not

PT

191

continue evaorating
192 i f m f u e l t r a p e d <( m spray acumulated+
m wall acumulated )
CE

193 evap =0;

194 m spray acumulated=m spray acumulated−
m s p r a y t o t a l ∗ dt ;
AC

195 m wall acumulated=m wall acumulated−m wall ∗

dt ;
196 end
197

198 %P e n e t r a t i o n and a n g l e f o r t h e s p r a y t o s t i m a t e
t he A l s i n o t h e r s c r i p t ( s e e s i m u l i n k

40
 ACCEPTED MANUSCRIPT

blocks )
199 t i m e f r o m i n j e c t i o n=time−t i m e i n j e c t e d ;
200 i f ( t i m e f r o m i n j e c t i o n )<t b k
201 x s p r a y =(Cd∗ s q r t ( abs ( P r a i l −P c y l i n d e r ) /
r h o f ) ) ∗( t i m e f r o m i n j e c t i o n ) ;

T
202 else

IP
203 x s p r a y = 2 . 9 5 ∗ ( ( abs ( P r a i l −P c y l i n d e r ) / r h o a
) ˆ 0 . 2 5 ) ∗ s q r t ( d i a m t e r i n j e c t o r ∗ abs (
t i m e f r o m i n j e c t i o n −t b k ) ) +(Cd∗ s q r t ( abs (

CR
P r a i l −P c y l i n d e r ) / r h o f ) ) ∗( t b k ) ;
204 end
205 a n g l e s p r a y=atan ( ( 3 . 6 2 7 6 ∗ s q r t ( r h o a / r h o f ) )

206

207 else
evap =0;
US
/ ( 3 + 0 . 2 8 ∗ ( l n / d i a m t e r i n j e c t o r ) ) ) ∗180/ p i ;
AN
208

209 end
210 %{
211 subplot (4 ,4 ,1) ;
M

212 ho ld on
213 p l o t ( time , t h e t a , ’ ∗ ’ ) , t i t l e ( ’ t h e t a ’ )
214 subplot (4 ,4 ,2) ;
ED

215 hold on
216 p l o t ( time , m wall , ’ ∗ ’ ) , t i t l e ( ’ m wall ’ )
217 subplot (4 ,4 ,3) ;
hold on
PT

218

219 p l o t ( time , m s p r a y t o t a l , ’ ∗ ’ ) , t i t l e ( ’ m s p r a y t o t a l ’ )
220 subplot (4 ,4 ,4) ;
hold on
CE

221

222 p l o t ( time , m spray , ’ ∗ ’ ) , t i t l e ( ’ m spray ’ )

223 subplot (4 ,4 ,5) ;
224 hold on
AC

225 p l o t ( time , N d r o p l e t s , ’ ∗ ’ ) , t i t l e ( ’ N d r o p l e t s ’ )
226 subplot (4 ,4 ,6) ;
227 hold on
228 p l o t ( time , a n g l e s p r a y , ’ ∗ ’ ) , t i t l e ( ’Bm’ )
229 subplot (4 ,4 ,7) ;
230 hold on

41
 ACCEPTED MANUSCRIPT

231 p l o t ( time , x s p r a y , ’ ∗ ’ ) , t i t l e ( ’ Bt ’ )
232 subplot (4 ,4 ,8) ;
233 hold on
234 p l o t ( time , a v e r a g e g a s v e l o c i t y , ’ ∗ ’ ) , t i t l e ( ’
average gas velocity ’ )

T
235 subplot (4 ,4 ,9) ;

IP
236 hold on
237 p l o t ( time , incr MFF , ’ ∗ ’ ) , t i t l e ( ’ incr MFF ’ )
subplot (4 ,4 ,10) ;

CR
238

239 hold on
240 p l o t ( time ,MFF, ’ ∗ ’ ) , t i t l e ( ’MFF ’ )
241 subplot (4 ,4 ,11) ;
242

243

244
hold on

subplot (4 ,4 ,12) ;
hold on
US
p l o t ( time , A l s , ’ ∗ ’ ) , t i t l e ( ’ A l s ’ )
AN
245

246 p l o t ( time , Sh , ’ ∗ ’ ) , t i t l e ( ’ Sh ’ )
247 subplot (4 ,4 ,13) ;
248 ho ld on
M

249 p l o t ( time , Re , ’ ∗ ’ ) , t i t l e ( ’ Re ’ )
250 subplot (4 ,4 ,14) ;
251 ho ld on
ED

252 p l o t ( time , rho a , ’ ∗ ’ ) , t i t l e ( ’ r h o a ’ )

253 subplot (4 ,4 ,15) ;
254 hold on
p l o t ( time , m spray acumulated , ’ ∗ ’ ) , t i t l e ( ’
PT

255

m spray acumulated ’ )
256 subplot (4 ,4 ,16) ;
hold on
CE

257

258 p l o t ( time , m wall acumulated , ’ ∗ ’ ) , t i t l e ( ’

m wall acumulated ’ )
259 %}
AC

260 end
261

262 fu nct ion s u r f a c e = fcn ( angle spray , x spray , x p i s t o n )

263 %This s c r i p t s t i m a t e t h e w a l l w e t t i n g a r e a with t h e
p a r a m e t r i c s e q u a t i o n s o f a cone and a p l a n e
264

42
 ACCEPTED MANUSCRIPT

265 x p i s t o n=abs ( x p i s t o n ) ;
266 %c l e r a n c e = 0 . 0 0 1 ;
267 x0 =77/1000; y0 =0; z0 =0;%p o s i t i o n o f th e i n y e c t o r
268 i n c l i n a t i o n d e g r e e =70;%a n g l e o f i n c l i n a t i o n o f t h e
inyector

T
269

IP
270 t i t a = 90− a n g l e s p r a y ; %a n g l e o f s p r a y
271 h = x s p r a y ; %p e n e t r a t i o n o f t h e s p r a y
D=6.25/1000+ x p i s t o n ; %p o s i t i o n o f t h e p i s t o n ,

CR
272

be c a r e f u l l with th e c o o r d i n a t e s
273

274 r=h/ tan ( t i t a ∗ p i /18 0) ; %r a d i o o f t h e cone

275

276

277
m = h/ r ;
US
p hi = −i n c l i n a t i o n d e g r e e ∗ p i /180 ;
[ u ,A] = meshgrid ( l i n s p a c e ( 0 , r , 1 0 0 ) , l i n s p a c e ( 0 , 2 ∗ pi
,100) ) ;
AN
278 [ t t ]= meshgrid ( l i n s p a c e (−h , h , 1 0 0 ) ) ;
279

280 %p a r a m e t r i z a t i o n o f t h e cone
M

281 X cone = ( u . ∗ c o s (A) ∗ c o s ( p h i ) − ( h/ r ) ∗u∗ s i n ( p h i )+

x0 ) ;
282 Y cone = u . ∗ s i n (A)+y0 ;
ED

283 Z cone = u . ∗ c o s (A) ∗ s i n ( p h i ) + ( h/ r ) ∗u∗ c o s ( p h i )+

z0 ;
284

%p a r a m e z i t r a t i o n o f t h e cone ’ s c u t with t h e plane ,

PT

285

but , t h e cone has not end

286 k=(D−z0 ) . / ( c o s (A) ∗ s i n ( p h i ) + ( h/ r ) ∗ c o s ( p h i ) ) ;
x1=(k ) . ∗ c o s (A) ∗ c o s ( p h i ) −(h/ r ) ∗( k ) ∗ s i n ( p h i )+x0 ;
CE

287

288 y1=(k ) . ∗ s i n (A)+y0 ;

289 z1=D. ∗ ( ones ( 1 0 0 , 1 0 0 ) ) ;
290 %v e c t o r o f t he h i g h t o f t h e cone
AC

291 p=[−h∗ s i n ( ph i )+x0 ,0+ y0 , h∗ c o s ( p h i )+z0 ] ;

292 o r i g i n =[x0 , y0 , z0 ] ;
293 vd=p−o r i g i n ;%v e c t o r d i r e c t o r o f p l a n e o f t h e f i n a l
cone
294

295 D2=p ( 1 ) ∗vd ( 1 )+p ( 2 ) ∗vd ( 2 )+p ( 3 ) ∗vd ( 3 ) ;

43
 ACCEPTED MANUSCRIPT

296 %l i n e c r e a t e d by th e cu t with t h i s p l a n e and t h e

piston plane
297 vd2=c r o s s ( vd , [ 0 , 0 ,D] ) ;
298 p2 =[(D2−D) /vd ( 1 ) , y0 ,D ] ;
299

T
300 x l i n e 2=p2 ( 1 )+vd2 ( 1 ) . ∗ t t ;

IP
301 y l i n e 2=p2 ( 2 )+vd2 ( 2 ) . ∗ t t ;
302 z l i n e 2=p2 ( 3 )+vd2 ( 3 ) . ∗ t t ;

CR
303

304 %c r e a t e s u r f a c e d e l i m i t a t e d by t h i s l i n e and t h e
p r e v i o u s cu t o f t h e cone
305 c ut =[ x1 ( : , 1 ) y1 ( : , 1 ) z1 ( : , 1 ) x l i n e 2 ( 1 , : ) ’ y l i n e 2

306

307
(1 ,:) ’ z line2 (1 ,:) ’ ];
US
%i n i z i a l i z a t i o n o f f i n a l matrix
c u t f i n a l x 1=ones ( 1 0 0 ) ; c u t f i n a l y 1=ones ( 1 0 0 ) ;
c u t f i n a l z 1=ones ( 1 0 0 ) ;
AN
308

309 j =1; k=1;

310 %The l o o p f i n d th e s u r f a c e o f t he c u t u n t i l t he
M

line ,
311 %But , th e l i n e i s not p a i n t e d b e c a u s e we have th e
i n i c i a l i z a t i o n of the
ED

312 %matrx c u t f i n a l
313 f o r i =1:100
314 i f cu t ( i , 1 )<=cu t ( i , 4 )
c u t f i n a l x 1 ( i , : ) =c u t ( i , 1 ) ;
PT

315

316 c u t f i n a l y 1 ( i , : ) =c u t ( i , 2 ) ;
317 c u t f i n a l z 1 ( i , : ) =c u t ( i , 3 ) ;
j=i ;
CE

318

319 k=k+1;
320 end
321 end
AC

322 %we use t he p o i n t o f c u t t h e l i n e and extend t o th e

f i r s t and f i n a l matrix
323 %matlab w i l l j o i n t t h i s p o i n t s a u t o m a t i c a l l y
324 c u t f i n a l x 1 ( 1 : j−k + 1 , : ) = c u t ( j , 1 ) ;
325 c u t f i n a l y 1 ( 1 : j−k + 1 , : ) = c u t ( j , 2 ) ;
326 c u t f i n a l z 1 ( 1 : j−k + 1 , : ) = c u t ( j , 3 ) ;

44
 ACCEPTED MANUSCRIPT

327 c u t f i n a l x 1 ( j : 1 0 0 , : ) = cut ( j , 1 ) ;
328 c u t f i n a l y 1 ( j : 1 0 0 , : ) = cut ( j , 2 ) ;
329 c u t f i n a l z 1 ( j : 1 0 0 , : ) = cut ( j , 3 ) ;
330 %F i n a l l y , we w i l l c a l c u l a t e t h e matrix
331 s u r f a c e 1=p o l y a r e a ( c u t f i n a l x 1 , c u t f i n a l y 1 ) ;

T
332 s u r f a c e=s u r f a c e 1 ( 1 ) ;

IP
333 i f isnan ( surface )
334 s u r f a c e =0;
end

CR
335

336

337 end
338

339

340

341
US
f u n c t i o n [ hc ,w]= c o n v e c t i o n ( t h e t a , P c y l i n d e r , T c y l i n d e r ,
rpm , x p i s t o n )
AN
342 %Estimate t he h eat t r a n s f e r c o e f f i c i e n t and th e a v e r a g e
speed i n s i d e t he c y l i n d e r
343 g l o b a l Pi
M

344 g l o b a l Tr
345 g l o b a l Vr
346 g l o b a l time
ED

347 %d e f i n i t i o n o f n u m e r i c a l v a r i a b l e s ( r e q u i r e d by MATLAB
t o c o m p i l e t he s c r i p t )
348 w=0; Pm=0;
PT

349

350 %Constants
351 bore =77/1000;
s tr o ke =85.8/1000;
CE

352

353 C= 0 . 0 0 3 5 ; %E m p i r i c a l c o e f f i c i e n t
354 m= 0 . 8 ; %E m p i r i c a l c o e f f i c i e n t
355 gamma= 1 . 4 ; %A d i a b a t i c c o e f f i c i e n t
AC

356 r =10.5; %Compresion r a t i o

357 x p i s t o n=abs ( x p i s t o n ) ;
358

359 Sp=rpm∗ s t r o k e / 3 0 ; %Mean speed o f t h e p i s t o n

360 Vd=p i /4∗ s t r o k e ∗ bore ˆ 2 ; %D is pla ce men t volume
361

45
 ACCEPTED MANUSCRIPT

362 %C1 and C2 c o e f f i c i e n t s ( Heywood )

363 i f t h e t a >148 && t h e t a <554 %exchange between 148 and 554
364 C1= 6 . 1 8 ;
365 C2=0;
366 e l s e i f t h e t a >=554 && t h e t a <=698 %c o m p r e s s i o n between

T
554 and 698

IP
367 C1= 2 . 2 8 ;
368 C2=0;
e l s e %combustion and e x p a n s i o n

CR
369

370 C1= 2 . 2 8 ;
371 C2=3.24 e −3;
372 end
373

374

375
US
%P r e s s u r e t e m p e r a t u r e and volumen f o r a r e f e r e n c e s t a t e
i f t h e t a >300 && t h e t a <350 | | Vr==0
Pi=P c y l i n d e r ;
AN
376

377 Tr=T c y l i n d e r ;
378 Vr=p i /4∗( x p i s t o n + ( 6 . 2 5 / 1 0 0 0 ) ) ∗ bo re ˆ 2 ;
379 end
M

380

381 Pm=((Vr / ( p i /4∗( x p i s t o n + ( 6 . 2 5 / 1 0 0 0 ) ) ∗ bo re ˆ 2) ) ˆgamma) ∗ Pi

; %I s e n t r o p i c p r e s s u r e
ED

382 i f Pm==P c y l i n d e r
383 w=1;
384 else
w=C1∗Sp+C2 ∗ ( (Vd∗Tr ) / ( Pi ∗Vr ) ) ∗( P c y l i n d e r −Pm) ; %
PT

385

a v e r a g e speed i n s i d e th e c y l i n d e r
386 end
w=abs (w) ;
CE

387

388

389 hc=C∗( bor e ˆ (m−1) ) ∗( P c y l i n d e r ˆm) ∗(wˆm) ∗( T c y l i n d e r

ˆ ( 0 . 7 5 − 1 . 6 2 ∗m) ) ; % he a t t r a n s f e r c o e f f i c i e n t
AC

390 %{
391 subplot (4 ,4 ,1) ;
392 ho ld on
393 p l o t ( time , t h e t a , ’ ∗ ’ ) , t i t l e ( ’ t h e t a ’ )
394 subplot (4 ,4 ,2) ;
395 hold on

46
 ACCEPTED MANUSCRIPT

396 p l o t ( time , P c y l i n d e r , ’ ∗ ’ ) , t i t l e ( ’ P c y l i n d e r ’ )
397 subplot (4 ,4 ,3) ;
398 hold on
399 p l o t ( time , Vr , ’ ∗ ’ ) , t i t l e ( ’ Vr ’ )
400 subplot (4 ,4 ,4) ;

T
401 ho ld on

IP
402 p l o t ( time , w, ’ ∗ ’ ) , t i t l e ( ’w ’ )
403 subplot (4 ,4 ,5) ;
ho ld on

CR
404

405 p l o t ( time , hc , ’ ∗ ’ ) , t i t l e ( ’ hc ’ )
406 subplot (4 ,4 ,6) ;
407 ho ld on
408

409

410
subplot (4 ,4 ,7) ;
hold on
US
p l o t ( time , x p i s t o n , ’ ∗ ’ ) , t i t l e ( ’ x p i s t o n ’ )

p l o t ( time , P c y l i n d e r , ’ ∗ ’ ) , t i t l e ( ’ P c y l i n d e r ’ )
AN
411

412 subplot (4 ,4 ,8) ;

413 hold on
414 p l o t ( time , T c y l i n d e r , ’ ∗ ’ ) , t i t l e ( ’ T c y l i n d e r ’ )
M

415 subplot (4 ,4 ,9) ;

416 hold on
417 p l o t ( time ,Pm, ’ ∗ ’ ) , t i t l e ( ’Pm ’ )
ED

418

419 subplot (4 ,4 ,10) ;

420 ho ld on
p l o t ( time ,We, ’ ∗ ’ ) , t i t l e ( ’We ’ )
PT

421

422 subplot (4 ,4 ,11) ;

423 ho ld on
p l o t ( time , Re , ’ ∗ ’ ) , t i t l e ( ’ Re ’ )
CE

424

425 subplot (4 ,4 ,12) ;

426 ho ld on
427 p l o t ( time , A l s , ’ ∗ ’ ) , t i t l e ( ’ A l s ’ )
AC

428 subplot (4 ,4 ,12) ;

429 hold on
430 p l o t ( time , T d r o p p r e v i o u s , ’ ∗ ’ ) , t i t l e ( ’ T d r o p p r e v i o u s ’
)
431 subplot (4 ,4 ,13) ;
432 ho ld on

47
 ACCEPTED MANUSCRIPT

433 p l o t ( time , v p r e v i o u s , ’∗ ’ ) , t i t l e ( ’ v previous ’ )

434 subplot (4 ,4 ,14) ;
435 hold on
436 p l o t ( time , d p r e v i o u s , ’∗ ’ ) , t i t l e ( ’ d previous ’ )
437 subplot (4 ,4 ,15) ;

T
438 hold on

IP
439 p l o t ( time , Yfs , ’ ∗ ’ ) , t i t l e ( ’ Yfs ’ )
440 subplot (4 ,4 ,16) ;
hold on

CR
441

442 p l o t ( time , Pfs , ’ ∗ ’ ) , t i t l e ( ’ Pfs ’ )

443 %}
444 end
445

446
US
f u n c t i o n [ z , p a l f a , t i m e i n j e c t i o n 2 ] = EOQ( t h e t a ,
m a s s i n j e c t e d r a t e , Q lv , evap , m wall ,
b u r n e d f u e l f r a c t i o n , Q heat , dt , a n g l e s p a r k )
AN
447 g l o b a l time
448 g l o b a l time2
449 global p
M

450 global flag

451 global flag2
452 global flag3
ED

453 global t wall

454 global time injection
455 global tr
global tj
PT

456

457 global mass not evaporated

458 k=0; K=0;
%v a r i a b l e s
CE

459

460 v=Q lv ; a l f a =0;

461 %a l f a=evap ; %a l f a c o e f f i c i e n t
462 %beta=Q heat ; %be ta c o e f f i c i e n t
AC

463 time=time+dt ; %time from i n j e c t i o n t o s p a r k

464 time2=time2+dt ; %time o f s i m u l a t i o n f o r c o n t r o l
465 e =2.71828182;
466 Q r a t e =0; z =0; w lambert =0;
467 a n g l e s p a r k =720+ a n g l e s p a r k ;
468 %i n i z i a l i z a t i o n o f g l o b a l v a r i a b l e s

48
 ACCEPTED MANUSCRIPT

469 i f 20 < t h e t a && t h e t a < 60

470 p=0; %t o t a l f u e l i n j e c t e d
471 t w a l l =0; %time a t s t a r t o f w a l l w e t t i n g
472 t r =0; %d u r a t i o n o f w a l l w e t t i n g
473 t j =0; %d u r a t i o n between s t a r t o f

T
i n j e c t i o n and s p a r k

IP
474 t i m e i n j e c t i o n =0;
475 f l a g =0;
f l a g 2 =0;

CR
476

477 f l a g 3 =0;
478 end
479 %i n j e c t i o n
480

481 time =0;

US
i f m a s s i n j e c t e d r a t e >0 && f l a g==0
%d u r a t i o n between s t a r t o f
i n j e c t i o n and s p a r k
f l a g =1;
AN
482

483 end
484 %time from w a l l w e t t i n g
485 i f m wall >0 && f l a g 2==0
M

486 t w a l l=time ; %time a t s t a r t o f w a l l w e t t i n g

487 f l a g 2 =1;
488 end
ED

489

490 %s p a r k
491 i f t h e t a >a n g l e s p a r k && f l a g 3==0
f l a g 3 =1;
PT

492

493 t j=time ; %d u r a t i o n o f w a l l w e t t i n g
494 t r=time−t w a l l ; %d u r a t i o n o f i n j e c t i o n
end
CE

495

496

497 i f m a s s i n j e c t e d r a t e >0
498 p=p+m a s s i n j e c t e d r a t e ∗ dt ; %t o t a l
AC

mass i n j e c t e d
499 t i m e i n j e c t i o n=t i m e i n j e c t i o n+dt ; %d u r a t i o n
of the i n j e c t i o n
500 end
501 %between s p a r k and combustion
502 i f t h e t a <a n g l e s p a r k && p>0

49
 ACCEPTED MANUSCRIPT

503 m a s s n o t e v a p o r a t e d=p−evap ∗ dt ;
504 if m a s s n o t e v a p o r a t e d <0
505 m a s s n o t e v a p o r a t e d =0;
506 end
507 end

T
508

IP
509 i f p==0 | | t r==0 | | t j ==0
510 z=−1/e ;
else

CR
511

512 %Parameters f o r t h e i n p u t o f th e lambe rt W f u c t i o n ( z )

513 a l f a=m a s s n o t e v a p o r a t e d / t r ;
514 beta=Q heat / ( Q lv ∗ t j ) ;
515

516

517
US
K=b u r n e d f u e l f r a c t i o n ∗p ;
k=(K/v ) ∗ e ˆ ( a l f a ∗ t r+beta ∗ t j ) ;
z =(( a l f a ∗k ) / ( p∗ e ) ) −(1/ e ) ;
end
AN
518

519 i f a l f a==0
520 p a l f a =0;
521 else
M

522 p a l f a =(p/ a l f a ) ;
523 end
524 end
ED
PT
CE
AC

50