1.

Engine Model

Operation of 4-stroke cycle engine:

Intake Compression Combustion Exhaust

25
Cylinder Pressure (bar) 20
XC125 4000rpm WOT 15
50 10
45 5
40 0
35 0 90 180 270 360 450 540 630 720 Crank
30
Angle (deg)

1
Using MATLAB/SIMULINK to simulink real engine, refer to Simulink Demonstrationas shown

in below Fig. The engine model includes：(1)set up the timing, (2)calculation of mass air

flow rate from throttle and manifold, (3)identify the air mass for one cycle, (4)calculation of torque,
and (5)vehicle dynamics. Input throttle positioncan estimate engine speed.

A practical engine model includes the subs-systems: (1) charging, (2) combustion,(3) heat
transfer, (4) friction, and (5) work done. Please refer to IMECE 2003 paper, by YY Wu et al.,
“Engine Modeling with Inlet and Exhaust Wave action for Real Time Control”

2
Practice:
5. To run the SIMULINK demonstration model “[Link]” for exercise.
(1) Change the throttle angle to: 10-20, 30-50, 60-90
(2) Change the “Drag Torque”.
 (3) Change the throttle with a “constant” of 10, 20, 30, 40, 50, 60, 70, 80, 90.
(4) Change the throttle with a “ramp” as slope 10.

Simplify the engine model “[Link]” by deleting the sub models: „valve timing‟,
„compression‟, „intake‟, „vehicle dynamics‟ and „drag torque‟. The new model is
“[Link]”. The engine model ([Link]) is used for torque calculation only. The
inputs are throttle angle and engine speed (rpm), the output is torque (Nm). It is like an engine
tested on a dynamometer.

Practice:
6. To run the simplified model “[Link]” for exercise.
(1) Change the throttle angle to: 10-20, 30-50, 60-90
(2) Change the “rpm” to: 3000, 4000, 5000, …, 10000.
(3) Change the throttle with a “constant” of 10, 20, 30, 40, 50, 60, 70, 80, 90.
(4) Change the throttle with a “ramp” as slope 10.

7. Using several „step‟ blocks to set engine speed. (1000, 3000, …,

10000rpm) parameters of „Step‟: Step time, Initial value, Finial value.
parameters of „Sum‟: Icon shape, List of signs.

3
 4

2. Subs-Systems of Engine Model:

˙Throttle

˙Intake manifold

˙Mass flow rate

˙Torque generation

 − mao
 ( ) mai mao

mai
 Intake
ManifoldCylinder

Fi g. M od e l fo r ca l cu l ati n g o f in t ak e ma s s fl o w r at e

：inlet air flow rate (g/s)

mai

：outlet air flow rate (g/s)

mao

Throttle:

()
23
mg
= ( )⋅ 2.821− 0.05231θ + 0.10299θ − 0.00063θ ai P m

P
= 1 for
g(Pm) 2amb
P≤
m

2
2 P
for
g=− m 2amb

P>
p PPP
() m m
P m amb m

Pm=manifold pressure (bar); Pamb=ambient pressure (bar)

Intake manifold:
d
=−
RT dt
P
mm m V

( ) ai ao m

where:
Pm：manifold pressureVm：manifold volume R：Air gas constant (287
J/kg/K) T：manifold temperature

Mass flow rate:

5
 2
2

= −0.366+ 0.08979 − 0.0337 + 0.0001

mao NPm NPm N Pm
where:
N = engine speed (rad/s)

Torque generation:
torque = -181.3 + 379.36*ma + 21.91*AF - 0.85*AF*AF + 0.26*S - 0.0028*S*S + 0.027*N -
0.000107*N*N + 0.00048*N*S + 2.55*S*ma - 0.05*S*S*ma
ma=mass of air (g/cycle),
mf=mass of fuel (g/cycle),
S=spark advance (deg CA),
N=engine speed (rad/s)

Practice:
8. To build up engine intake flow model

9. To build up engine torque model

6
7

3. Modeling a Real Engine:

3.1 Charging Model (Moskwa’s Ph.D. thesis, Heywood’s, chapter 14.3.3,
Appendix C)
Four types of models for calculating details of intake and exhaust flows have been developed
and used:
(1) Quasi-steady model: used for control, such as EGR valve control. The disadvantage is too many
empirical equation s.
(2) Filling and empting model: used for performance prediction. Using mass continuity and energy
equation or entropy only. Do not need momentum equation.
(3) Wave action model: One dimensional flow with momentum equation. The calculating time is
pretty long. When the ratio of length/diameter is greater than 10, the effect of wave action should
not be neglected.
(4) CFD model: 3-dimensional viscous flow. The calculating time is very long. Used for a new
engine design.

Take the engine mass flow rate as a function of A (area), Po (ambient pressure), and Pm/Po
(pressure ratio of manifold pressure to ambient pressure).

m = f A P
P

(,,)
m
o
P
o

If the manifold pressure Pm is steady state, as only function of engine speed (rpm), then it is a
“quasi-steady model”.
If the manifold pressure Pm is unsteady, then it is a “filling and empting
model”. If the pressure wave in the pipe is considered, then it is a “wave action
model”.
Filling and Empting Model

 − mao
 ( ) mai mao

mai

Intake
ManifoldCylinder

T h e a ir m a ss f lo w in to en gin e m od el

(g/s)；mao
mai Throttle:
(g/s)

C A( ) k PRI
P0
1

mai d = ⋅ ⋅ ⋅ ⋅2

θ
RT
0
1
1
k
−
1

P: 2
k k k
⎜ ⎟ ⎪ ⎢ ⎡ ⎟ ⎥ ⎪
⎜ ⎝⎛ ⎟ ⎠⎞ ⎨⎧ ⎢⎢ ⎣ ⎟ ⎠⎞ ⎥⎥ ⎦⎤ ⎬⎫
− 2 P
2 ⎜
P ⎜ ⎝⎛
t
⎜ ⎛ ⎟ t

⎝ + ⎠⎞
If 1
>k t = ⋅ ⋅−1
 PRI
−
0
Pk 1
⎪
Pk P ⎩ ⎪
01
2 + ⎭
− 0

⎜ ⎛ ⎟
⎝ + ⎠⎞
k 1

P 2( 1)
:
2− PRI
=k
⎟
⎜ ⎛ ⎠⎞
⎝ +
If 1
≤k t k 1
Pk
01

Where
Pt：Throttlepressure at the minimum cross-sectional area
P0：Stagnation Pressure
Cd：discharge coefficient

A(θ )：the flow area of the orifice

PRI：Affect the nature of the pressure ratio of the flow rate
2 2 2
⎥ ⎢ ⎡ ⎟ θ ⎪
dDd⋅ ⎪ ⎧ ⎥
⎥ ⎦⎤ ⎢ ⎣ ⎟ ⎠⎞ dDd ⎬⎫
⎥
1 ⎨ ⎥ ⎦⎤
⎥ ⎦⎤
1
2
2 12
2
cos( )
⋅ ⎜
⎝⎛
⎢ ⎡⎟
⎢ ⎣ ⎠⎞ ⎢ ⎡⎟
D d ⎢ ⎣ ⎠⎞

=− ⋅−1 + ⋅−⋅
−
()θ 1
⎜
0 ⎝⎛
0
2 A D
⎜ ⎛
⎜ ⎝ +
D D ⎪
1 cos( ) 2 ⎪ ⎭
⎩
+ ⋅ − sin 1
2
θθ
2
 1
cos( ) θ θ ⎪ ⎧ ⎥ ⎢ ⎡ ⎟
⎪
⎨ ⎥ ⎦⎤ ⎢ ⎣ ⎟ ⎠⎞ cos( )
⎬⎫
+ 2
2

Dd
θ

−⋅
0 10
−

⋅−⋅
2 sin 1 ⎜ ⎛
⎜ ⎝ +
cos( )

When ⎢
d
⎥
θ ⎪ cos( ) θ θ ⎪
− ⎣⎡ 0 ⎩ ⎦⎤ ⎭
1
D 0

( ) 00
2
2
21
1

θ cos cos θ −θ ⎪
≥⋅ ⎬⎫
D
2
⎪ ⎧ ⎥
⎨ ⎥ ⎦⎤

2 ⎥
( )2 ⎥ ⎦⎤

⎢ ⎡⎟ dDd⋅
D d ⎢ ⎣ ⎠⎞

⎢ ⎡⎟
⎢ ⎣ ⎠⎞
⎜
⎝⎛
1
Aθ ⎜
− ⎝⎛

= ⋅ − sin 1 − ⋅−⋅
 D 2 D
⎪ ⎪
⎩ ⎭ 1

D: throttle bore diameter (D=0.024m)

d: throttle shaft diameter
(D=0.00694m)
RT
θ 0: zero flow angle ( P

d ( ) ai ao
Intake Manifold: =  −  dt
mm m V
θ 0=0.926o)
m

9
PV
⋅

⋅⋅
mωη= π

md
ev
ao
RT⋅⋅

ηv=(-0.366/ωe+ 0.08979Pm- 0.0337Pm2 + 0.000001ωe Pm)*15

Heywood‟s p.312 need to be corrected.

Where
Pm：Manifold pressure
Vm：Manifold volume (127 cc)
R：Air gas constant (287 J/kg/K)

T：Intake manifold temperature

ωe：engine speed
ηv：volumetric efficiency
Vd：displaced volume (m3)
：Inlet mass flow rate
mai

：outlet mass flow rate

mao
Wave Action Model

Continuity equation,
ρρρ

+
u u dF
∂ ∂

+=0
t .
∂ ∂ …………………
………
x F dx

Momentum equation,
22

ρρρ
∂
G()()

+ ∂+up
u u
dF
ρ

++=
t 0 ………..
∂ ∂

a = kRT = k
Energy equation p F
2
ρ dx
x
 10

3.2 Combustion Model

(1) Torque Function (Moskwa)
Indicated Torque:
Indicated Torque Tind (Nm):
Tind = (TF)(MAC)
TF: torque function; MAC: mass of air per cylinder (g)
−
0.258
⎥
⎥ ⎦⎤

⎢ ⎡
⎢ ⎣ ⎟−

⎜ ⎞
⎝⎛

= − P TF ω ( )
− m 14.39*1.09/ 6* 44.8 3.723 e
0.088

Tbr = Tind - Tfr 9.77

1.013 ⎠

Tbr: brake torque (Nm); Tfr: friction torque (Nm)

Pm：manifold pressure；ωe：engine speed

(2) Wieve Function (Heywood, p.390)

 11
m
⎢ ⎡⎟ +1 ⎥
⎢ ⎣ ⎠⎞ θ θ ⎥ ⎦⎤

xb a = − −
0
⎜ ⎛
⎝ Δ− θ
1 exp

x= mass fraction burned, θ 0= start of combustion, Δθ= total combustion duration, b

a, m = adjustable parameters, typically, a = 5, m = 2

θθ θθ
dx 0 0 1 b − +

+−1
m

= =m
y a
m
( ) exp( ( ) )

−
d a
θθd θd θd

heat release rate

dQ
=()
hr
yQm
HV f
d
θ

(3) Flame Propagation Model (Lumley)

 qb

Tb, P, hb, mb

(dmu/dt)hu 3.3 Heat Transfer Model

Tu, P, hu, mu qu

A constant wall temperature Tw (K) and a varying mean charge temperature T (K). The heat transfer
rate is:

dQ
=−

()w
ht
hA T T
θ
d
where A is the surface area of combustion chamber. Mean charge temperature T is computed from
the state equation of gas with known cylinder pressure.
12
A lot of empirical equations can be used to calculate the instantaneous convective heat transfer
coefficient, h (W/m2/K). We compared several empirical equations with our experimental data and
obtained the following correlation with curve fitting. It is a function of cylinder pressure P (bar),
gas temperature T (K) and mean gas velocity Cm (m/s):

()
6 0.635 1.450 0.052
−

h p T Cm
=1.92*10 +1.4

The constants in this equation were obtained from experimental regressions. The mean gas
velocity Cm was obtained from Woschni‟s correlation:
=2.28*
Cm Sp

where
S pis piston mean speed (m/s).

Lumley, p.106
h = S uρc
tp
dQ
=−

( ) pw
ht
hA T T
θ
d

3.4 Friction Model

Friction Torque (Heywood, chapter 13, p.722):

2

⎥
⎥ ⎦⎤

⎢ ⎡⎟
⎢ ⎣ ⎠⎞
NN
⎜
⎝⎛
T ⎜
⎞ ⎝⎛

fr Vd π = +

⎟−
1.97 0.15 1000 0.01 1000 *100000* /(4 )

Friction Mean Effective Pressure (Honda’s Model):

The total friction loss of a four-stroke S.I. motorcycle engine proposed by Honda

is: ( ) 23
 2

pmf = C1N +C C

mf pis the total friction mean effective pressure (MPa), N is the engine speed (rpm).

Coefficient C1 is used to estimate the pumping loss, constant C2 is proportional to the viscosity of
lubrication oil, and C3 is dimensional coefficient. For a 125 cc single cylinder engine with two
valves: C1 =2.95*10-9(MPa/rpm2), C2 =0.185 (MPa), C3 =0.887.
The friction torque Tfr is obtained by:

τ ( )/(4π)

fr = pmfVd

13
Practice:
10. To set a new torque model and friction model ([Link]) by using Moskwa‟s
torque function and modified Heywood‟s friction torque model..

3.5 Work Done

Instantaneous cylinder pressure P and corresponding incremental change of cylinder volume V
can further express the incremental work done by the system:

dW
=
dV
P
θ dθ
d
Cylinder pressure is obtained by solving the following ordinary differential equation, which is
derived from equation (19).
dQ dQ
hr ht
dp 1 1
p dV
 ⎟ ⎜
( ) ⎠⎞ ⎝⎛

=−+−−
k k
θ θ θ dθ
d V d V d
The instantaneous volume of engine cylinder and its rate change can be expressed [10, 11],
respectively, by:
⎜ ⎛
⎜ ⎝ +−+−
V r
⎟
2 ( )⎟ ⎠⎞

d θθ

V
= 1 cos(2 )

−
2 1 4
Cr 1 cos( ) L

Vr
d
dV
⎜ θθ ⎟
⎝⎛ ⎠⎞

= + sin(2 )
sin( )
θL
d 2 2
where Vd is engine displacement volume, Cr is engine compression ratio, r is the crank radius of the
engine, and L is the length of connecting rod.
The indicated mean effective pressure and torque can be expressed by:

pmi W Vd
=/

τ ( )/(4π )

ind = pmiVd

τ = τ −τ
br ind fr
Practice:
11. To build up an engine model for a 125cc motorcycle.
14
% * Engine Parameters
Vd = 125e-6; % engine displacement volume [m^3]
Vm = 127e-6; % intake manifold volume [m^3]
% * Engine Throttle Parameters.
Cd=0.85; % discharge coefficient of throttle
D=0.024; % throttle bore diameter (m)
d=0.0069; % throttle shaft diameter (m)
theta0=0.926*pi/180; % zero flow angle, 0.926 (deg)
% * Environment Parameters.
k = 1.3; % Specific heat ratio of mixture
R=287; % Gas Constant of air [J/(kg-K)]
% * Ambient Initialization Parameters
T0=298; P0=1e5;

15
12. To calculate bhp [Link] (brake specific fuel consumption)

m
bsfc f
(kW);
bhp = 2πτN / 60000 bhp

= (g/h/kW)

13. Put the output data to workspace and then save the data to a file. (format short g)

16
Save the data in Workspace to an excel file: in “Command Window”, key
in save -ascii [Link] filename1

14. Plot bsfc map (contour of bsfc on the coordinate of torque vs. speed).
contour(Z) is a contour plot of matrix Z treating the values in Z as heights above a plane. A
contour plot are the level curves of Z for some values V. The values V are chosen automatically.
contour(X,Y,Z), X and Y specify the (x,y) coordinates of the surface as for SURF. contour(Z,N)
and contour(X,Y,Z,N) draw N contour lines, overriding the automatic value. contour (Z,V) and
contour(X,Y,Z,V) draw LENGTH(V) contour lines at the values specified in vector V. EXAMPLE:
contour(x, y, z)
contour(x, y, z, 10)
contour(x, y, z, [200 250 300 350 400 450 500 550 600 650])
clabel(contour(x, y, z, [200 250 300 350 400 450 500 550 600 650]))

reference：「bsfc_map.m」or 「bsfc_map_s41.m」

x=[1000 1600 2200 2800 3400 4000 4600 5200 5800];

y=[6.8 13.6 20.4 27.2 33.8 40.6 47.4 54.2 61 67.8 74.6 81.4];
z=[635.0 678.4 463.4 699.1 592.9 667.9 630.6 698.4 751.1
635.7 500.1 463.4 567.9 592.9 524.8 630.6 500.5 637.8
541.4 443.8 407.6 500.3 494.6 381.6 522.5 428.6 521.1
447.2 387.4 350.1 432.7 393.4 351.9 411.1 392.7 407.8
17
352.9 331.1 294.3 301.4 295.1 322.2 303.0 356.8 393.1 332.2 301.8
280.8 283.9 279.4 304.9 304.4 337.9 378.4 311.4 297.0 267.3 266.3
263.6 287.5 305.8 328.4 363.3 322.4 283.4 253.9 248.7 247.9 270.8
304.2 319.0 348.2 333.5 269.8 269.8 258.8 255.2 290.8 314.5 328.8
318.8 333.5 358.0 303.2 268.8 262.5 310.9 324.8 338.6 340.2 333.5
358.0 336.7 271.9 295.0 330.9 327.7 333.7 340.2 333.5 358.0 336.7
317.9 322.6 330.9 327.7 333.7 340.2];
clabel(contour(x,y,z, [250 300 350 400 450 500 550 600 650]))

18
Why do all engine manufactures use computer simulations?
1. Engine testing is very expensive ($2000 - $4000 per day).
2. Developing prototypes of new engines is expensive and time consuming.
3. Modeling can reveal the root cause of behavior where testing provides only
behavior.
4. Simulation development often leads to „new concepts‟ through a
better understanding of engine processes.

Why do they bother to test engines at all?

1. Some engine phenomenon are not yet well enough understood to capture in a
simulation (turbulence, chemical kinetics, etc.).
2. Complete and accurate modeling of some processes requires massive amounts of
computational time (i.e. Computational Fluid Dynamics).
3. Models and simulations often need some „calibration‟ data from real world tests.
4. Testing provides not only performance information, but also reveals design faults
and durability issues.

source: Dr. Joel Hiltner

19
ADVISOR

Simulink model: Advisor3.2/model/ BD_CONV.mdl

Data files:
1. FC_*.m Fuel Converter 2. ESS_*.m Energy Storage System 3. TX_*.m Transmission 4.
GC_*.m Generator/Controller 5. MC_*.m Motor/Controller 6. PTC_*.m Power Train
Control 7. TC_*.m Torque Coupler 8. VEH_*.m Vehicle 9. WH_*.m Wheel 10. ACC_*.m
Accessory 11. CYC_*.m Driving Cycle 12. EX_*.m Exhaust Aster Treatment

20