In-Cylinder Combustion Analysis of a Yamaha R6 FSAE Engine: 
Achieving Improved Engine Performance and Efficiency  
through Burn Rate Combustion Diagnostics 
 
by 
 
Jacob Timothy Fredrick 
 
 
 
 
A thesis submitted to the Graduate Faculty of 
Auburn University 
in partial fulfillment of the 
requirements for the Degree of 
Master of Science 
 
Auburn, Alabama 
May 9, 2011 
 
 
 
 
Keywords: internal, combustion, reciprocating, engine, burn rate analysis, FSAE 
 
 
Copyright 2011 by Jacob Timothy Fredrick 
 
 
Approved by 
 
Peter D. Jones, Chair, Woltosz WEMS Professor of Mechanical Engineering 
Daniel W. Mackowski, Associate Professor of Mechanical Engineering 
Jay M. Khodadadi, Alumni Professor of Mechanical Engineering 
Roy W. Knight, Assistant Professor of Mechanical Engineering 
 
 
 
   
 ii
Abstract 
 
 
The actual three-dimensional, turbulent combustion process that occurs in a modern, 
spark-ignited internal combustion engine is very complex and difficult to model. However, 
since cylinder pressure is directly related to the combustion event, in-cylinder pressure traces 
sampled and averaged over multiple cycles can be utilized to provide the burn rate of the air 
and fuel mixture in the cylinder with respect to crankshaft position. A design of experiment 
was created to examine the influence of ignition timing, fuel rail pressure, injection timing, 
and air to fuel ratio on the start of combustion, burn duration, and combustion centroid at 
various engine operating points. Once the factors controlling the three burn rate profile 
parameters were understood, they were used to manipulate the combustion event and 
improve engine performance and efficiency.  Through the in-cylinder pressure analysis 
performed, an indicated specific fuel consumption of 0.469 lbm/hp-hr and an indicated mean 
effective pressure of 43 psi were realized for the 15% load, 3,000 rpm operating point.  An 
indicated specific fuel consumption of 0.356 lbm/hp-hr and an indicated mean effective 
pressure of 88 psi were achieved for the 60% load, 3,000 rpm operating point.  The study was 
performed on a 599 cc Yamaha R6, four stroke, inline four cylinder, naturally aspirated, 
spark ignition engine with a restricted intake for use in Formula SAE competition.  
   
 iii
Acknowledgements 
 
 
I am very grateful to my parents, Tim and Wanda Fredrick, who have pushed me to 
give nothing but my best in everything.  It is by their example that I acquired the discipline 
and determination necessary to succeed in graduate school.  I would also like to show my 
appreciation to my teammates on Auburn FSAE.  They have been a great group to work with, 
and I?m fortunate to have been one of their teammates.  In particular, I would like to thank 
Dr. Peter Jones, Brian Keyser, and Sean Baker for their dedication, support, and expertise 
given throughout this experiment as well as Kistler and BEI Encoders for their generous 
sponsorship of this research.  In addition, the service of Dr. Daniel Mackowski, Dr. Jay 
Khodadadi, and Dr. Roy Knight on my advisory committee was much appreciated.  I would 
also like to acknowledge Paul Wiczynski, Chris Byus, Ryan Siefring, Arvind Suresh, Pete 
Woon, Joe Reynolds, and Chris Pollitt for the exposure to internal combustion engine testing 
they granted me in my time at Cummins, Inc.  Last, but certainly not least, I would like to 
thank Greg Nelson for all of the guidance he gave me during my young engineering career at 
Harley-Davidson Motor Company and for encouraging and greatly influencing my decision 
to pursue my Master?s Degree in Mechanical Engineering here at Auburn University. 
 
   
 iv
Table of Contents 
 
 
Abstract ..................................................................................................................................... ii?
Acknowledgements .................................................................................................................. iii?
List of Tables ........................................................................................................................... vi?
List of Figures ......................................................................................................................... vii?
List of Abbreviations .................................................................................................................x?
Chapter 1: Background ..............................................................................................................1?
Chapter 2: Working Principles of the System............................................................................5?
2.1 Spark Ignition Engine Pre-Mixed Flame Development ...................................................5?
2.2 Calculating Performance and Efficiency Indicators from the P-V Diagram ...................6?
2.3 Quantifying Combustion Variability ...............................................................................8?
Chapter 3: Modeling Techniques and Equations .....................................................................10?
3.1 Combustion Modeling Approaches ...............................................................................10?
3.2 Derivation of the Mass Fraction Burn Rate ...................................................................13?
Chapter 4: Instrumentation & Programming ...........................................................................20?
4.1 Piezoelectric Pressure Sensors .......................................................................................20?
4.2 Incremental Angular Encoders ......................................................................................27?
4.3 Data Acquisition ............................................................................................................28?
4.4 Data Integrity and Analysis ...........................................................................................30?
Chapter 5: Component Machining & Dynamometer Setup .....................................................34?
   
 v
5.1 Cylinder Head Drilling Procedure .................................................................................34?
5.2 Angular Encoder Component Machining and Dynamometer Setup .............................39?
Chapter 6: Design of Experiment ............................................................................................42?
6.1 15% Load, 3,000 RPM DOE .........................................................................................42?
6.2 60% Load, 3,000 RPM DOE .........................................................................................47?
Chapter 7: Experimental Results .............................................................................................54?
7.1 15% Load, 3,000 RPM Experimental Results ...............................................................54?
7.2 60% Load, 3,000 RPM Experimental Results ...............................................................60?
Chapter 8: Conclusions ............................................................................................................68?
References ................................................................................................................................70?
Appendix ..................................................................................................................................72?
LabView 2009 Data Acquisition Block Diagram ................................................................72?
LabView 2009 Data Acquisition Front Panel ......................................................................73?
MATLAB R2009a Data Analysis Code ..............................................................................74?
MATLAB R2009a JANAF Table Thermodynamic Data Code ..........................................87?
   
 vi
 List of Tables 
 
 
Table 1-1: Yamaha R6 engine specifications ........................................................................... 1?
Table 4-1: Measurement task and recommended crank degree resolution ............................. 27?
Table 4-2: Test readiness and verification instructions .......................................................... 31?
Table 6-1: 15% load, 3,000 RPM DOE matrix ....................................................................... 42?
Table 6-2: 15% load, 3,000 RPM SOC ANOVA ................................................................... 44?
Table 6-3: 15% load, 3,000 RPM combustion centriod ANOVA .......................................... 45?
Table 6-4: 15% load, 3,000 RPM combustion duration ANOVA .......................................... 46?
Table 6-5: 60% load, 3,000 RPM DOE matrix ....................................................................... 47?
Table 6-6: 60% load, 3,000 RPM SOC ANOVA ................................................................... 48?
Table 6-7: 60% load, 3,000 RPM combustion centriod ANOVA .......................................... 50?
Table 6-8: 60% load, 3,000 RPM combustion duration ANOVA .......................................... 51?
Table 7-1: 15% load, 3,000 RPM performance and efficiency comparison ........................... 58?
Table 7-2: 60% load, 3,000 RPM performance and efficiency comparison ........................... 64?
   
 vii
List of Figures 
 
 
Figure 1-1: Location of combustion on in-cylinder pressure trace ........................................... 2?
Figure 1-2: Combustion energy release aspects ........................................................................ 3?
Figure 2-1: Four-Stroke P-V Diagram ...................................................................................... 6?
Figure 3-1: Cylinder control volume ...................................................................................... 14?
Figure 4-1: Kistler 5001 charge amplifier and 6052A high temperature pressure sensor ...... 21?
Figure 4-2: Standard deviation of pressure as a function of crankshaft position .................... 23?
Figure 4-3: Variation in pressure readings from 40 degrees ATDC to EVO ......................... 24?
Figure 4-4: Thermal shock distortion on the log-log P-V diagram ........................................ 24?
Figure 4-5: Thermal shock effect on the mass fraction burned curve .................................... 25?
Figure 4-6: National Instruments USB-6210 data acquisition unit ........................................ 29?
Figure 4-7: MoTeC M400 ECU .............................................................................................. 29?
Figure 4-8: IMEP control chart ............................................................................................... 31?
Figure 4-9: Data analysis program flowchart ......................................................................... 33?
Figure 5-1: Tilt table and cylinder head jig adjustment on mill ............................................. 35?
Figure 5-2: Mill setup for second cylinder head machining operation ................................... 36?
Figure 5-3: Pressure sensor communication passage .............................................................. 37?
Figure 5-4: Mounting sleeve location in cylinder head .......................................................... 37?
Figure 5-5: Disassembled upper crankcase, transmission, & power cylinder components .... 38?
Figure 5-6: Disassembled lower crankcase, cylinder head, & clutch components ................. 38?
   
 viii
Figure 5-7: Crankshaft encoder assembly ............................................................................... 39?
Figure 5-8: Crankshaft encoder bracket construction ............................................................. 40?
Figure 5-9: Crankshaft encoder mating shaft .......................................................................... 41?
Figure 5-10: Engine dynamometer test stand ......................................................................... 41?
Figure 6-1: 15% load, 3,000 RPM SOC half-normal plot ...................................................... 43?
Figure 6-2: 15% load, 3,000 RPM combustion centroid half-normal plot ............................. 44?
Figure 6-3: 15% load, 3,000 RPM combustion duration half-normal plot ............................. 45?
Figure 6-4: 15% load, 3,000 RPM DOE mass fraction burned curves ................................... 46?
Figure 6-5: 60% load, 3,000 RPM SOC half-normal plot ...................................................... 48?
Figure 6-6: 60% load, 3,000 RPM combustion centroid half-normal plot ............................. 49?
Figure 6-7: 60% load, 3,000 RPM combustion duration half-normal plot ............................. 50?
Figure 6-8: 60% load, 3,000 RPM DOE mass fraction burned curves ................................... 51?
Figure 6-9: Mass fraction burned dependency on ignition timing .......................................... 52?
Figure 6-10: Mass fraction burned dependency on lambda .................................................... 53?
Figure 7-1: 15% load, 3,000 RPM ISFC as a function of SOC placement ............................. 54?
Figure 7-2: 15% load, 3,000 RPM IMEP as a function of SOC placement ........................... 55?
Figure 7-3: 15% load, 3,000 RPM ignition delay ................................................................... 55?
Figure 7-4: 15% load, 3,000 RPM ISFC as a function of combustion centroid placement .... 56?
Figure 7-5: 15% load, 3,000 RPM IMEP as a function of combustion centriod placement .. 57?
Figure 7-6: 15% load, 3,000 RPM combustion centroid placement characteristics ............... 57?
Figure 7-7: 15% load, 3,000 RPM fired and motored cycle ................................................... 58?
Figure 7-8: 15% load, 3,000 RPM P-V diagram .................................................................... 59?
Figure 7-9: 15% load, 3,000 RPM log-log P-V diagram ........................................................ 59?
   
 ix
Figure 7-10: 15% load, 3,000 RPM mass fraction burned curve ............................................ 60?
Figure 7-11: 60% load, 3,000 RPM ISFC as a function of SOC placement ........................... 61?
Figure 7-12: 60% load, 3,000 RPM IMEP as a function of SOC placement ......................... 61?
Figure 7-13: 60% load, 3,000 RPM ignition delay ................................................................. 62?
Figure 7-14: 60% load, 3,000 RPM ISFC as a function of combustion centroid placement .. 62?
Figure 7-15: 60% load, 3,000 RPM IMEP as a function of combustion centriod placement 63?
Figure 7-16: 60% load, 3,000 RPM combustion centriod placement characteristics ............. 63?
Figure 7-17: 60% load, 3,000 RPM fired and motored cycle ................................................. 64?
Figure 7-18: 60% load, 3,000 RPM P-V diagram .................................................................. 65?
Figure 7-19: 60% load, 3,000 RPM log-log P-V diagram ...................................................... 66?
Figure 7-20: 60% load, 3,000 RPM mass fraction burned curve ............................................ 66?
   
 x
List of Abbreviations 
 
Acronyms 
 
AF  air to fuel ratio 
ANOVA analysis of variance 
BDC  bottom dead center 
COV  coefficient of variation 
CR   compression ratio 
DOE  design of experiment 
ECM  engine control module 
FSAE  Formula Society of Automotive Engineers 
IMEP  indicated mean effective pressure  
ISFC  indicated specific fuel consumption 
LNV  lowest normalized value 
LPP  location of peak pressure 
MAP  manifold absolute pressure 
MBT  max brake torque 
NA  naturally aspirated 
PCP  peak cylinder pressure 
P-V  pressure vs. volume diagram 
SI  spark ignition 
   
 xi
SOC  start of combustion 
STD  standard deviation 
TDC  top dead center 
 
Symbols 
a   Wiebe function efficiency parameter 
A   instantaneous combustion chamber surface area 
ch
A   cylinder head combustion chamber surface area 
f
A   area of flame front 
p
A   piston crown surface area 
B   cylinder bore 
v
c   specific heat at constant volume 
chv
c
 
  charge specific heat at constant volume  
rv
c
 
  residual specific heat at constant volume 
C   Woschni heat transfer calibration constant 
1
C   average cylinder gas velocity empirical calibration parameter 
2
C   average cylinder gas velocity empirical calibration parameter 
e   Reynolds number exponent 
ff   turbulent flame factor 
h   heat transfer coefficient averaged over combustion chamber surface area 
i   index of summation 
t
i    ignition timing in terms of crank angle BTDC 
   
 xii
l   connecting rod length 
L   piston stroke 
m   mass 
Wiebe function form factor 
b
m   burned mass 
ch
m   charge mass 
f
m   mass of fuel inducted per cycle 
r
m   residual mass 
n   upper bound of summation 
N   crankshaft rotational speed 
P   instantaneous cylinder pressure 
m
P   instantaneous motored cylinder pressure 
r
P   working fluid reference pressure 
Q   heat transfer 
ht
Q   spatially averaged combustion chamber heat transfer  
HV
Q   heating value of fuel 
r   crankshaft throw 
s   distance between crankshaft axis and piston pin axis   
p
S   mean piston speed 
t   time 
T   mean gas temperature 
r
T    working fluid reference temperature 
   
 xiii
wall
T   mean combustion chamber wall temperature 
u   specific internal energy 
bg
u   specific internal energy of burned exhaust gas 
ch
u   specific internal energy of unburned charge 
U   internal energy 
1
U   laminar flame speed 
t
U   turbulent flame speed 
V   instantaneous cylinder volume 
c
V   cylinder clearance volume 
d
V   displaced cylinder volume 
r
V   working fluid reference volume 
s
V   swirl velocity 
W   work transfer 
ic
W
,
  indicated work per cycle 
b
x   mass fraction burned 
f
?   fuel conversion efficiency 
?   crank angle 
o
?   start of combustion crank angle 
??   crank angle duration of combustion 
?   relative air/fuel ratio 
u
?   density of unburned gas 
   
 xiv
?   characteristic charge velocity 
 
Subscripts 
i   index of summation 
min   minimum 
 
   
1 
 
Chapter 1: Background 
 
 
Measuring pressure inside the cylinder of a reciprocating engine was first utilized 
around the late 1700?s when James Watt and others were developing steam engines.  The 
early engine indicators used during that time were primitive mechanical devices.  Although 
the measurement technology has advanced greatly, the goal is still the same: understand the 
energy release process in the cylinder to optimize power and efficiency.  Utilizing in-cylinder 
pressure measurements is a widely used, established method of developing internal 
combustion engines in industry, but has yet to be used by Formula SAE (Society of 
Automotive Engineers) teams due to the cost and complications associated with performing 
combustion diagnostics.  Therefore, being able to comprehend how design choices affect the 
in-cylinder energy release is a substantial advantage over the competition.  The specifications 
for the Yamaha R6 engine under study are listed in Table 1-1.  To comply with FSAE rules, 
the engine was run with a 20 millimeter restriction between the throttle and plenum. 
Table 1-1: Yamaha R6 engine specifications 
Engine Configuration Inline four cylinder, 4-stroke 
Engine Displacement 599 cc 
Compression Ratio 13.1:1 
Valves per Cylinder 4 
Bore x Stroke 67 x 42.5 mm 
Intake Valve Opening 39 degrees BTDC 
Intake Valve Closing 65 degrees ABDC 
Exhaust Valve Opening 64 degrees BBDC 
Exhaust Valve Closing 24 degrees ATDC 
   
2 
 
There are two main approaches to analyzing cylinder pressure data.  One is a burn 
rate analysis and the other is a heat release analysis.  The burn rate analysis performed in this 
experiment is the standard choice for use with gasoline engines and produces a normalized 
quantity on a scale of 0 to 100 percent mass fraction burned.  Heat release analysis is 
commonly used with diesel engine combustion studies and produces an absolute energy 
quantity [1].  The heat release analysis tends to be more sensitive to measurement errors than 
the burn rate analysis since the magnitudes of the errors are reduced by normalizing.  
The pressure rise in an internal combustion reciprocating engine is due to three main 
factors: combustion, change in cylinder volume, and heat transfer to or from the containing 
surfaces [2].  The pressure change due to cylinder volume can be taken out of the equation by 
comparing the engine motoring curve (engine is externally powered) to the firing curve.  The 
difference in pressure between these two diagrams is the effect of combustion.  The 
crankshaft position where the two curves first start to diverge is the start of combustion 
(SOC). 
-400 -300 -200 -100 0 100 200 300 400
0
100
200
300
400
500
600
700
Crankshaft Angle (Degrees)
C
y
l
i
n
d
e
r
 P
r
e
ssu
r
e
 (
p
si
)
Pressure vs. Crankshaft Position
 
 
Fired Cycle
Motored Cycle
 
Figure 1-1: Location of combustion on in-cylinder pressure trace 
The effect of heat transfer in all of this is something that can be estimated using 
various methods.  Quite often the heat transfer is modeled as convective heat transfer while 
employing a heat transfer coefficient that includes the effects of radiation.  The amount of 
Combustion 
Intake 
  Compression Expansion 
Exhaust 
   
3 
 
heat transfer happening in the cylinder is very complex to model since it is dependent on the 
in-cylinder combustion species, swirl and squish, and flame front geometry.  The empirically 
based and proven Woschni correlation [3], which is still a preferred way to model the in-
cylinder heat transfer, was utilized in this analysis.  
A First Law of Thermodynamics energy balance performed on the cylinder provides a 
differential equation from which the in-cylinder energy release profile can be calculated.  The 
energy-release aspects of combustion can be characterized by the parameters shown in  
Figure 1-2.  The SOC, combustion duration, and centroid of combustion (crankshaft angle 
where 50% burn occurs, often abbreviated CA50) are the critical parameters.  A shorter burn 
duration will provide a higher torque output and reduced combustion variability.  The 
placement of the centroid and SOC has a great effect on ISFC (indicated specific fuel 
consumption) and IMEP (indicated mean effective pressure).  The mass burn rate at the 
centroid location is also important to note. 
-30 -20 -10 0 10 20 30
0
0.2
0.4
0.6
0.8
1
Crankshaft Angle (Degrees)
Ma
s
s
 F
r
a
c
ti
o
n
 B
u
r
n
e
d
, X
b
Mass Fraction Burned vs. Crankshaft Position
 
Figure 1-2: Combustion energy release aspects 
Combustion 
Centroid 
Combustion 
Duration 
SOC
   
4 
 
Combustion in internal combustion engines occurs in a three-dimensional, time 
dependent, turbulent flow with a fuel containing a blend of hydrocarbons which have poorly 
understood combustion chemistry [4].  Even while running at constant speed, load, and 
temperature, there will be considerable variations in performance from cycle-to-cycle.  
Cycle-to-cycle and cylinder-to-cylinder variation are due to variations in the amount of air 
and fuel delivered, spark plug, injector signal wire resistance (from temperature variation), 
valve timing (spring rate variation, natural frequency, etc.), fuel density (temperature and/or 
air bubbles), and buildup on the injector tip.  This variability can result in compromised spark 
timing, torque and power, lower fuel economy, and engine roughness.  Taking the average of 
the pressure traces over a span of 100 cycles is the common practice used to minimize the 
effects of the variations on test measurements.   
Although the combustion process in an internal combustion engine is quite 
complicated and riddled with variation, utilizing in-cylinder pressure traces to perform a 
quasi-steady, closed loop analysis offers insight on what is happening inside the cylinder.  By 
examining the key combustion performance parameters (IMEP, combustion phasing, cyclic 
variation, and energy release), cylinder pressure combustion analysis allows for the 
evaluation of various air handling, fueling strategies, combustion chamber shapes, 
compression ratios, spark plug parameters, valve timing, and ignition timing in order to make 
sound design decisions that maximize the engine?s potential for improved performance, 
economy, drivability, reliability, and cost. 
   
5 
 
Chapter 2: Working Principles of the System 
 
2.1 Spark Ignition Engine Pre-Mixed Flame Development 
 
 The Yamaha R6 engine under study is naturally aspirated (NA) meaning that it draws 
air into the intake system at atmospheric temperature and pressure.  Once the air-fuel mixture 
has been inducted into the cylinder and compressed, combustion initiates at the spark plug 
and a plasma channel is created which is roughly at 60,000 K, 300 bar and expands at 
supersonic velocity [5].  The chemical reactions only become possible on the kernel surface 
at temperatures below 3,000 K.  Within roughly 0.1 ms, the kernel has grown to a diameter 
of 1 mm, and grows at a laminar velocity.  The kernel eventually becomes wrinkled and 
turbulent after its diameter reaches the local integral scale of turbulence.  The kernel must 
reach a critical diameter (the local integral scale of turbulence) in order for self-sustained 
combustion to begin.  It is because of these events early on in the combustion process that 
ignition delay exists.  Turbulent combustion is a wrinkled flame that behaves as a laminar 
flame locally.  The laminar flame speed is dependent on the type of fuel, AF ratio, and 
pressure and temperature of the mixture [6].   
 The ignition delay and laminar flame speed come into play when determining spark 
timing.  For a given mixture of fuel and air inside of a cylinder, there exists an optimum 
spark timing for which the engine will produce maximum torque.  This timing is referred to 
as maximum brake torque (MBT) timing.  It is a compromise between starting combustion 
   
6 
 
too early in the compression stroke when the work transfer is actually to the cylinder gases 
(as well as potentially subjecting the cylinder to knock) and completing combustion too late 
in the expansion stroke, thus lowering peak expansion stroke pressures.  Knock is the pre-
ignition of the end gas (mixture on the opposite side of the flame front from the spark plug) 
that produces a very violent and undesired detonation and supersonic pressure oscillations 
that are very destructive to an engine?s power cylinder components.  The end gas auto-ignites 
from abnormally high temperatures and pressures in the cylinder caused by aggressive 
timing.  Ignition timing is a balance between the work transferred to the piston before TDC 
and the work transferred to the piston before BDC.  The specific timing is dependent on 
engine speed and load. 
2.2 Calculating Performance and Efficiency Indicators from the P-V Diagram 
 In-cylinder pressure data acquired over the operating cycle of an engine can be used 
to calculate the work transfer from the gas to the piston.  A pressure-volume (P-V) diagram 
of the engine cycle is very useful in calculating the indicated work per cycle.  It is simply the 
integral of the P-V diagram given by Eq. (2.1) [6]. 
0 2 4 6 8 10 12 14 16 18 20
0
100
200
300
400
500
600
700
Cylinder Volume (cubic inches)
C
y
l
i
nd
er
 P
r
e
s
s
u
re 
(ps
i
)
P-V Diagram
 
Figure 2-1: Four-Stroke P-V Diagram 
Expansion 
Compression 
 
Intake 
 
Exhaust 
Combustion 
   
7 
 
 
?
? PdVW
ic,
     (2.1) 
 
 Gross indicated work per cycle is only the work done on the piston over the positive 
work loop bound by the compression and expansion strokes.  Net indicated work per cycle 
includes the losses in the pumping loop (charge induction and exhaust) as well.  The 
pumping loop is negative in NA engines since some amount of work is required to intake 
charge and exhaust burnt combustion gases.  Indicated quantities only include the impact of 
the compression, combustion, and expansion processes on engine performance.  It is a 
measure of the available power and torque at the crankshaft and does not include frictional 
losses through the drivetrain.  A useful way to compare closed-cycle engine performance 
between different engines is to calculate the IMEP by Eq. (2.2) [6].  IMEP takes the engine?s 
displacement into account in order to provide a fair comparison between engines of different 
sizes. 
     
d
ic
V
W
IMEP
,
?      (2.) 
 The fuel conversion efficiency is the ratio of work output to the total available energy 
of the inducted fuel and is defined by Eq. (2.3) [6].  The indicated specific fuel consumption 
(ISFC) is the fuel flow rate per unit power output.  Eq. (2.4) expresses this quantity on a per 
cycle basis. 
     
HVf
ic
f
Qm
W
,
??      (2.3) 
     
ic
f
W
m
ISFC
,
?     (2.4) 
   
8 
 
 Tuning for MBT can be accomplished by examining the average IMEP versus the 
crankshaft angle corresponding to the centroid of combustion and SOC.  To achieve the 
maximum efficiency, each cylinder should be optimized individually if the ECM (engine 
control module) is capable of providing different ignition timings to individual cylinders.  
However, there will always be some phasing efficiency loss that exists.   
2.3 Quantifying Combustion Variability 
Combustion variability can be quantified by examining the IMEP of each cycle.  The 
most basic approach is to calculate the standard deviation (STD) of the IMEP which gives a 
representation of cyclic variation and engine stability.  The variability can also be expressed 
as a coefficient of variation (COV) by relating the STD as a percentage of the mean IMEP in 
Eq. (2.6) [7]. 
   
?
?
?
?
?
n
i
i
n
IMEPIMEP
nIMEPSTD
1
2
)1(
)(
    (2.5) 
   %100
 
 ??
IMEP
IMEPSTD
IMEPCOV     (2.6) 
 Typically COV?s of IMEP above 3-5% will result in drivability issues that are 
detectable by the driver.  Misfires and partial burns can be detected by calculating the lowest 
normalized value (LNV) of IMEP in Eq. (2.7) [7].  The lowest IMEP is normalized by the 
mean IMEP.  An LNV less than 80% indicates a partial burn, and a LNV less than 0% 
indicates a misfire. 
   %100 
min
??
IMEP
IMEP
IMEPLNV     (2.7) 
 Misfires occur when the flame propagation is never initiated properly.  This may be 
due to insufficient ignition energy, conditions at the spark plug at the time of ignition, 
   
9 
 
excessive residuals, AF ratios outside of combustion limits, high compression pressures, low 
temperatures, high mean flow around the spark plug, and excessive plug fouling.  Partial 
burns occur when the flame never fully propagates across the combustion chamber.  Partial 
burns can be a result of burn durations that are too long, low compression pressures, retarded 
spark timing, and poorly mixed charge.  Partial burns can also cause irregular misfires [7]. 
   
10 
 
Chapter 3: Modeling Techniques and Equations 
 
3.1 Combustion Modeling Approaches 
 
Many different types of combustion models have been conceived.  They are all 
subject to a tradeoff between the simplicity of the model and the accuracy of the results that 
it produces.  There are three categories of approaches to model the combustion process [3]: 
? Zero-Dimensional Models 
? Use empirical heat-release model 
? Time is the only independent variable 
 
? Quasi-Dimensional Models 
? Use a separate submodel for turbulent combustion 
 
? Multi-Dimensional Models 
? Numerically solve equations for mass, momentum, energy, and species 
conversion to predict flame propagation 
? 1, 2, or 3-dimenssional      
The zero and quasi-dimensional models are simplistic, easy to employ into an engine 
model, and are quite useful for parametric engine modeling.  However, they are not linked to 
combustion chamber geometry explicitly.  When combustion chamber geometry is important 
or being developed, a multi-dimensional model should be used. 
   
11 
 
A rather well known zero-dimensional combustion model was nicely presented by 
Heywood [6].  The model uses the three zones including an unburned gas zone, an adiabatic 
burned gas core, and a burned gas thermal boundary layer.  Conservation of mass and energy 
are utilized to determine the conditions in the burned and unburned gas.  In order to do this, 
however, the mass fraction burned must be calculated.  This can be done with the Wiebe 
function [6].  The function can be fitted to experimental data to greatly reduce the number of 
combustion model inputs into engine simulation software. 
    
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
?1
exp1
m
o
b
ax
?
??
   (3.1) 
Where: 
  ? = crank angle 
  ?o = start of combustion 
  ?? = total combustion duration 
  a = efficiency parameter 
  m = form factor 
Ferguson also utilizes the Wiebe function in his engine model [8].  Ferguson?s model 
is much more complex and includes the effects of heat loss, blow by, and equivalence ratio.  
Ultimately, six ordinary differential equations are required to be simultaneously, numerically 
integrated.  The six ordinary differential equations are the change of pressure, burned and 
unburned temperature, work, heat loss, and enthalpy loss with respect to crank angle.  
The requirement for arbitrarily determined burn rate information for these models is 
not very convenient.  Quasi-dimensional models try to predict the burn rate information by 
assuming a spherical flame front geometry, and by using information about the turbulence as 
an input.  The equation at the forefront of this technique proposed by Stone [4] is shown in 
   
12 
 
Eq. (3.2).  The product of the density of the unburned gas, area of flame front, and turbulent 
flame speed equals the product of the density of the unburned gas, area of flame front, and 
laminar flame speed multiplied by a turbulent flame factor (ff). 
ffUAUA
dt
dm
futfu
b
1
?? ??     (3.2) 
This simple approach can be made much better with additional turbulence 
consideration.  The effects of ignition delay and the pressure rise due to compression on the 
length scale of turbulence as well as squish and swirl need to be analyzed as well. 
Since even very complex multi-dimensional models in commercial codes have 
trouble producing accurate simulations, a very powerful tool needed for the improvement of 
combustion analysis is the incorporation of experimental cylinder pressure data.  The 
pressure inside of the cylinder is directly related to the combustion process which makes it 
very useful in determining the mass fraction burn profile.  One of the earliest models was 
created by Rassweiler and Withrow in 1936 [9].  Despite its simplicity, it is still widely used.  
The approach only requires cylinder pressure and volume geometric details.  However, this 
leads to several shortcomings.  The model assumes that the polytropic constant is constant, 
while in reality, it varies from compression to expansion and changes during the combustion 
process.  Since the model is ultimately a scaling relationship, it doesn?t take into account 
incomplete combustion or the accuracy of cylinder pressure and mass flow rate data.  
Therefore, the effect of heat transfer is not included in the model.  The flow of cylinder 
charge into crevices such as between the piston, rings, liner, head gasket gap, and spark plug 
threads is also ignored.   
Because of the shortcomings of the Rassweiler and Withrow model, other two region 
models with a burned and unburned gas region have since been proposed.  However, they are 
   
13 
 
computationally complicated and time-consuming.  A simple, more accurate one-zone model 
was created by Gatowski in 1984 [10].  Cheung and Heywood then tested the validity, 
accuracy, and robustness of the model in 1993 [11].  This experiment uses a similar approach 
that tracks the charge mass burned and residual mass fraction on the basis of an individual 
cycle.  It permits the investigation of partial burns, misfires, and heat losses through 
empirical correlations.  The degree of burn rate variability from cycle-to-cycle can also be 
more accurately quantified. 
3.2 Derivation of the Mass Fraction Burn Rate 
 Since an internal combustion engine is not a closed system, it is not a heat engine in 
the thermodynamic sense of the term.  An engine is best analyzed as an open system that 
exchanges heat and work with the surrounding environment.  The reactants (fuel and air 
mixture) flow into the cylinder, and the products (exhaust gases) flow out [6].  An equation 
for the energy-release of a reciprocating internal combustion SI engine can be derived from 
the First Law of Thermodynamics [12]. 
    
dt
dV
P
dt
dQ
dt
dW
dt
dQ
dt
dU
????    (3.) 
 Where: 
  
dt
dU
 = time rate of change of internal energy of the system 
  
dt
dQ
 = time rate of change of heat transfer (heat-release from combustion  
   minus the heat transfer to chamber walls) 
  
dt
dW
 = time rate change of work done by the system 
  
dt
dV
P  = time rate change of work done under quasi-static expansion 
   
14 
 
 
Figure 3-1: Cylinder control volume 
 
The temperature changes which occur around minimum and maximum cylinder 
volumes are not primarily a result of heat transfer interactions [6].  During combustion, the 
temperature rise occurs to stabilize the internal energy as the species change from reactants to 
products.  The heat transfer model used by Heywood [6] that is based off of the Woschni 
correlation [3] is a convective heat-transfer rate to the combustion chamber walls shown in 
Eq. (3.4).  It includes the effects of piston movement, swirl, and combustion.  This 
relationship is correlated with experimental data and gives an estimation of the heat flux. 
    )(
wall
ht
TTAh
dt
dQ
??      (3.4) 
Where: 
  A = instantaneous combustion chamber surface area  
  h = heat transfer coefficient averaged over the chamber surface area  
  T = mean cylinder gas temperature  
  Twall = mean chamber wall temperature  
Heat 
Loss 
CV 
Work 
Heat 
Addition 
   
15 
 
 The heat transfer coefficient (in W/m?K) is calculated by taking the cylinder bore as 
the characteristic length, calculating the average in-cylinder gas velocity, utilizing the ideal 
gas law, and including curve fit parameters in Eq. (3.5) [6].  The mean cylinder charge 
temperature can be calculated through the ideal gas law since the molecular weight of the 
reactants and products are nearly identical [13].  This gives a temperature very close to the 
actual mass averaged cylinder temperature. 
    
eeee
TPCBh
62.175.01
26.3
??
? ?    (3.5) 
Where: 
  e = Reynolds number exponent (unit less data fit parameter) 
  B = cylinder bore [m] 
  P =  cylinder pressure [Pa] 
  T = mean cylinder gas temperature [K] 
  C = calibration constant (unit less correction factor) 
  ? = characteristic charge velocity (piston movement, swirl, and combustion) 
 Eq. (3.6) [6] calculates the average cylinder gas velocity (in meters per second), 
which includes the effects of swirl in higher-speed engines.  It has been shown to provide 
acceptable velocity predictions in both compression and spark ignition four-stroke, water-
cooled engines.  Three separate equations emerge for the combustion and expansion period, 
gas exchange period, and the rest of the cycle.   
    )(
21 m
rr
rd
p
PP
VP
TV
CSC ?
?
?
?
?
?
?
???    (3.6) 
Where: 
  C1 = 2.28, C2 = 3.24 x 10
3?
  - combustion and expansion period   
   
16 
 
  C1 = 
p
s
S
V
417.028.6 ? ,  C2 = 0  - gas exchange period  (3.7) 
  C1 = 
p
s
S
V
308.028.2 ? ,  C2 = 0  - rest of the cycle   (3.8) 
  Sp = 2LN - mean piston speed [m/s]     (3.9) 
   (L = piston stroke, N = crankshaft rotational speed) 
  
s
V = swirl velocity [m/s] 
  Vd = displaced volume [m?] 
  Tr = working fluid temperature at start of combustion [K] 
  Pr = working fluid pressure at start of combustion [Pa] 
  Vr = working fluid volume at start of combustion [m?] 
  Pm = motored cylinder pressure [Pa] 
  P = cylinder pressure [Pa] 
 The time rate of change of work transfer done by the system can be calculated by 
converting the rate of change of cylinder volume with respect to time to the rate of change of 
cylinder volume with respect to crankshaft angle. 
    
?
?
d
dV
NP
dt
dV
P 2?      (3.10) 
 To calculate the cylinder volume, piston position, and combustion chamber area at 
any crankshaft position, the following three equations can be utilized respectively [6]. 
    )(
4
2
srl
B
VV
c
????
?
    (3.1) 
    
2/1222
)sin(cos ?? rlrs ???    (3.12) 
    )( srlBAAA
pch
????? ?     (3.13) 
   
17 
 
 
Where: 
  Vc = cylinder clearance volume 
  B = cylinder bore 
  l = connecting rod length  
  r = crankshaft throw 
  s = distance between crankshaft axis and piston pin axis 
  Ach = cylinder head combustion chamber surface area  
  Ap = piston crown surface area 
 The internal energy of the cylinder gases can be split by the burned and fresh charge 
mass fractions.  Then the time rate of change of internal energy can be related to the change 
in crankshaft angle. 
    
chbbgb
uxmumxmuU )1( ????    (3.14) 
    ? ?
chbbgb
uxmumx
d
d
N
dt
dU
)1(2 ???
?
?   (3.15) 
Where: 
  m = mch + mr  - mass of fresh charge plus residual mass 
  xb = mass fraction burned 
  ubg = specific internal energy of burned exhaust gas 
  uch = specific internal energy of unburned charge 
 The product rule can be used to expand Eq. (3.15) by relating the rate of change of 
mass and internal energy to the rate of change of crankshaft angle. 
  ??
?
?
?
?
?
?
??????
chbbgbchbbgb
uxux
d
d
m
d
dm
uxuxN
dt
dU
)1())1((2
??
?  (3.16) 
   
18 
 
   
??
????
?
d
dT
dT
du
x
d
dx
u
d
dT
dT
du
x
d
dx
u
uxux
d
d
ch
b
b
ch
bg
b
b
bg
chbbgb
)1(
)1(
????
???
  (3.17) 
 If blowby is ignored, the equation is simplified greatly.  Taking dm/d? = 0 reduces the 
rate of change of internal energy equation further. The definition of specific heat at constant 
volume given by Eq. (3.19) can be used to simplify the equation even more.  This reveals the 
heat-release in Eq. (3.20).   
  
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?????
??
?
d
dT
dT
du
x
dT
du
xm
d
dx
uumN
dt
dU
ch
b
bg
b
b
chbg
)1()(2  (3.18) 
     
dT
du
Tc
v
?)(      (3.19) 
   
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
?
d
dT
mc
releaseheat
d
dx
uumN
dt
dU
v
b
chbg
???????
)(2    (3.20) 
 Now equations (3.4), (3.10), and (3.20) can be combined to form a relation for the 
mass fraction burn rate with respect to the change in crankshaft angle. 
  
?
?
??
?
d
dV
NPTTAh
d
dT
mc
d
dx
uumN
wallv
b
chbg
2)()(2 ???
?
?
?
?
?
?
??  (3.21) 
 Eq. (3.21) can be solved for the mass fraction burn rate with respect to change in 
crankshaft position relation shown below. 
   
)(
)(
2
chbg
rvrchvch
ht
b
uum
d
dT
cmcm
N
dt
dV
P
dt
dQ
d
dx
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
  (3.22) 
   
19 
 
 A Taylor Series was utilized to numerically integrate the differential equation and 
obtain the mass fraction burned curve with respect to crankshaft position.  The 
thermodynamic properties of the unburned and burned gas mixtures necessary to calculate 
the burn rate were acquired from a separate program written by Dr. Peter Jones [14] which 
uses polynomial curve fits from the NASA Equilibrium Code for JANAF table 
thermodynamic data. 
 
   
20 
 
Chapter 4: Instrumentation & Programming 
4.1 Piezoelectric Pressure Sensors 
By the late 1800?s the need for engine indicators capable of accurately measuring 
cylinder pressure in internal combustion engines with higher peak cylinder pressures gave 
way to optical indicators.  Devices like the planimeter were later invented to find the area of 
the P-V diagrams created by the indicators.  The electronic age gave way to new engine 
indicator devices becoming commercially available.  Most of them consisted of a cathode ray 
oscilloscope with additional hardware for measuring cylinder pressure and crankshaft 
position.  Improvements in digital electronics allowed for high-speed analogue to digital 
conversions with electronic storage.  Once data processing via a computer was possible, 
better instrumentation and signal conditioning systems gave way to more accurate data and 
analysis.  These advanced electronic systems are now standard equipment for manufacturers 
and researchers of internal combustion engines, and the analysis methods used are widely 
used and proven. 
Piezoelectric pressure sensors are the preferred equipment in acquiring reciprocating 
engine cylinder pressure.  They are constructed of quartz which makes them very durable 
even in the harshest of environments, and they are capable of measuring highly dynamic 
pressures in the combustion chamber with great accuracy and repeatability.  When the quartz 
crystal (SiO2) deforms due to a pressure change in the cylinder, it sends a charge to the 
   
21 
 
charge amplifier through low noise co-axial cables.  The charge is a result of the change in 
electric polarization which is proportional to the mechanical deformation.  This is called the 
piezoelectric effect which was discovered by Jacques Curie in 1880 [15].  The charge 
amplifier which was patented by Kistler in the 1950?s, converts the charge into a voltage that 
the data acquisition unit can read.  However, the pressure signal is AC coupled.  This means 
that the shape of the signal and relative amplitude are known, but the DC offset is not.  
Because of this, a procedure known as pegging is needed.  Cylinder pressure pegging takes 
the pressure signal and corrects it to some know pressure at a given crank angle in order to 
convert it into an absolute pressure.  Most often the MAP (manifold absolute pressure) at 
BDC will be used to correct the pressure signal as was done in this experiment, but the 
important thing is to be consistent in whatever approach is used.  The magnitude of PCP 
(peak cylinder pressure) is so much greater than that of the air handling events that the data 
can contain large errors in the pumping loss loop in the P-V diagram without a pegging 
method.  Absolute pressure is necessary to accurately obtain PCP, polytropic coefficients, 
temperature estimates, burn rate calculations, and pumping losses.  However, it is not 
necessary to determine the location of peak pressure (LPP), IMEP, or the COV of IMEP [7]. 
   
Figure 4-1: Kistler 5001 charge amplifier and 6052A high temperature pressure sensor 
 
   
22 
 
There are two types of piezoelectric pressure sensors, dry and wet transducers.  The 
wet transducers utilize water cooling to minimize the effect of thermal shock on the pressure 
transducer.  Thermal shock is the short-term drift of the sensor from exposure to hot 
combustion gases.  As the temperature of the transducer body changes, it deforms, causing a 
fake signal to be superimposed on the real signal.  This will show up as distortion in the 
cylinder pressure trace and will be evident in the mass fraction burned rate (particularly in 
the tail) and in the expansion and exhaust events on the log-log P-V diagram.  The larger 
thread size of the water cooled design makes the sensor more difficult to fit into a small bore 
cylinder head (like passenger car engines).  Modern dry sensors are designed to minimize the 
effects of thermal shock through various crystal cuts, and offer a temperature range for which 
the transducer sensitivity published on the calibration certificate is valid.  However, it can 
never be completely compensated for and was experienced during this experiment. 
The smaller size of dry transducers, allows for minimal intrusion on the actual 
combustion chamber geometry.  The smaller crystal size also results in a lower natural 
frequency.  This enhances the transducer?s ability to measure knock and combustion noise as 
well as its ability to acquire measurements at higher engine speeds [15].  However, they are 
subject to greater changes in sensitivity at higher operating temperatures due to the sensor?s 
construction-related factors like diaphragm deformation. 
The degree of error due to thermal shock is also associated with the operating point in 
the engine?s map is being run (ignition timing, engine speed and load).  Changes in engine 
load will also result in zero-shift of the pressure signal in dry transducers.  This is an issue 
when switching from fired to motored curves and vice versa [15].  Thermal shock can be 
detected by a simple test.  The variation of pressure at -180? and +180? crank angles (in 
   
23 
 
relation to TDC) can be compared.  If the variation at +180? is much greater than at -180?, 
thermal shock is occurring.  Thermal shock will also cause unreasonably low pressure 
readings in the P-V diagram pumping loop.  Locating the transducer?s measuring face near 
the intake valves substantially reduces the heat flow load on the transducer when compared 
to mounting it near the exhaust valves.  In this experiment, the transducer was mounted in 
between the intake and exhaust valve.  Figure 4-2 shows a typical plot for the standard 
deviation of the pressure trace for the 15% load, 3,000 RPM case.  Little difference was seen 
in the variation of pressure between the -180? and +180? crankshaft positions.  The major 
variation occurred shortly after TDC. 
-400 -300 -200 -100 0 100 200 300 400
0
5
10
15
20
25
Crankshaft Angle (Degrees)
C
y
l
i
nder P
r
es
s
u
r
e
 S
t
a
ndard D
e
v
i
at
i
o
n
 (
p
s
i
)
Standard Deviation of Pressure vs. Crankshaft Position
 
Figure 4-2: Standard deviation of pressure as a function of crankshaft position 
 It seems as though thermal shock effects from the burned gases have the most 
influence on random sensor variation during the expansion stroke in between 40 degrees 
ATDC and exhaust valve opening (EVO, 64 degrees BBDC).  Figure 4-3 displays the mean, 
median, and mode of 100 consecutive pressure traces in which this phenomenon can clearly 
be seen. 
   
24 
 
-400 -300 -200 -100 0 100 200 300 400
0
50
100
150
200
250
300
Crankshaft Angle (Degrees)
M
ean,
 M
edi
an
,
 M
ode of
 
C
y
l
i
nder
 P
r
es
s
u
r
e
 (
p
s
i
)
Statistical Check for Cylinder Pressure Data
 
 
Mean
Median
Mode
 
Figure 4-3: Variation in pressure readings from 40 degrees ATDC to EVO 
The variations seen in the expansion stroke pressure values were evident in other 
plots as well.  The log-log P-V diagrams all had some distortion in the later portions of the 
expansion stroke and into the beginning of the exhaust event as seen in Figure 4-4.  The 
expansion stroke should have a linear slope slightly higher than the compression stroke slope.  
However, with the thermal shock present, the ideal values for the expansion polytropic 
exponents were not realized due to unwanted curvature at the end of the expansion stroke. 
-0.2 0 0.2 0.4 0.6 0.8 1 1.2
1
1.5
2
2.5
Log of Cylinder Volume (cubic inches)
Log of
 
C
y
l
i
nder P
r
es
s
u
r
e
 (ps
i
)
Log-Log P-V Diagram
 
Figure 4-4: Thermal shock distortion on the log-log P-V diagram 
   
25 
 
 The thermal shock also affected the mass fraction burned curve.  The thermal shock 
to the pressure sensor resulted in the mass fraction burned curves lacking a tail at the end of 
combustion (EOC).  This can be seen in Figure 4-5 where experimental data has been fitted 
to a Wiebe function for reference.  However, since this experiment was mainly concerned 
with the location of SOC and the combustion centroid, the problematic readings near EOC 
were tolerable. 
-15 -10 -5 0 5 10 15 20 25 30 35
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Crankshaft Angle (Degrees)
M
a
s
s
 F
r
a
c
ti
o
n
 Bu
r
n
e
d
, X
b
Mass Fraction Burned vs. Crankshaft Position
 
 
Experimental Data
Wiebe Function
 
Figure 4-5: Thermal shock effect on the mass fraction burned curve 
Errors in the cylinder pressure readings can also arise from transducer drift which is 
the undesirable change in transducer output over time that is not a function of the measured 
variable.  The transducer will drift to ?10 volts over some period of time on the order of 
hours, minutes, or even seconds depending on the condition of the cabling, connections, and 
charge amplifier used.  Any condensation that is allowed to collect in the cables and 
   
26 
 
connections will cause excessive drift.  The cables can be baked for two hours at 125?C to 
eliminate any condensation.  Dirty or corroded connections can also add paths for current 
leakage from the measurement system.  These must be kept very clean with contact cleaner 
to avoid issues.  A well maintained measurement system will take hours to experience this 
drift and not effect measurements.  A well implemented pegging procedure will eliminate the 
effects of long-term drift as well.  Typically, the resistance between the center conductor and 
the outer shield should be kept greater than 
13
10  ohms [7].  The choice of the charge 
amplifier?s time constant is also important.  Both drift and time constant simultaneously 
affect the voltage output of the charge amplifier [17].  One or the other will dominate.  
Charge amplifiers will have short, medium, and long time constant settings to choose from.  
Most often the short time constant is utilized for monitoring the cylinder pressure of a 
modern engine.  There was very little sensor drift experienced throughout this experiment. 
Electrical noise will also cause issues in acquiring good data.  A filter is sometimes 
required to be implemented into the data acquisition to smooth the pressure trace and 
eliminate noise spikes in the data.  The noise will also be visible in the P-V diagram in the 
form of oscillations, especially in the low pressure regions of the diagram like the pumping 
loop.  However, care must be taken not to filter the pressure signal so much that knock is no 
longer detectable.  Mechanical noise from the engine?s reciprocating parts can also cause 
problems with transducer readings.  Accelerations as high as 1000 g?s are commonly 
encountered in modern engines [15].  Even a transducer that is over-torqued the slightest bit 
will exhibit a changed sensitivity which throws off the calibration of the transducer. 
 
   
27 
 
4.2 Incremental Angular Encoders 
The angular encoder that monitors crankshaft position mounts to the side of the 
engine with a specially machined bracket.  The optical encoder used detects radial position 
using a light gate and optical pickup.  This type of encoder has become the industry standard 
when accurate crankshaft displacement measurements are required.  The optical signal 
received by the optical pickup is converted into a digital voltage signal which the data 
acquisition unit can read.  The encoder utilized for this testing had a resolution of 360 pulses 
per revolution which is preferred for general combustion work.  Encoders with up to 3600 
pulses per revolution may be required when analyzing engine knock or combustion noise.  
Table 4-1 below gives Rogers? [15] recommended resolution for various combustion work.   
Table 4-1: Measurement task and recommended crank degree resolution 
Task Resolution (crank degrees) 
Knock <0.2 
IMEP, Friction Measurements 1 
Combustion Noise <0.2 
Heat Release Calculations 1 
Calculation of Polytropic 1 
Injection Timing Curve Edges 0.1 
Direct Analysis of Pressure Curve 1 
Rotational Vibration Analysis <1 
 
The encoder?s shaft is coupled to the crankshaft through a bellows assembly to 
accommodate any angular misalignment due to the mating shaft run out.  However, any 
mechanical resonance in the shaft encoder can cause oscillations to show up in the mass 
   
28 
 
fraction burned curve.  Large diameter bellows can help to reduce these effects, but have a 
shorter life than smaller diameter bellows.  The data processing program must know the 
angular offset of the encoder relative to true TDC.  The encoder sends two pulses to count 
every degree increment and to determine the direction of rotation.  Another pulse occurs once 
every cycle on the third channel to reset the degree position counter in the data acquisition 
program.  The zero point that the encoder sees must be phased to true TDC.  The encoder 
also only records 360 degree positions before resetting while a 4-stroke engine cycle operates 
over a duration of 720 degrees.  Therefore, the data processing program had logic embedded 
in it to determine which revolution is which in the engine?s cycle.  TDC is most often found 
by locating the PCP on the motored pressure trace.  PCP actually occurs slightly before TDC 
due to heat transfer.  This phasing difference is called the ?loss angle?.  This can be 
anywhere from 0 to 0.5 degrees at lower engine speeds and up to 3 degrees at higher speeds.  
The loss angle was found by adjusting the pressure trace so that the compression stroke 
polytropic exponent on the log-log P-V diagram matched the calculated ratio of specific 
heats.  
4.3 Data Acquisition 
It is necessary for the in-cylinder analysis that the in-cylinder pressure signal be 
coupled with crankshaft angular position since most of the information derived from the 
pressure trace is related to cylinder volume and crankshaft position.  Therefore, the data 
acquisition program is setup to use each pulse from the crankshaft encoder (per single degree 
increment) as a trigger to acquire a cylinder pressure data point, and then write them both to a 
ASCII file format.  This allows for cylinder pressure data to be captured in the angular 
domain.  The data can then be scanned at a constant rate fast enough to avoid aliasing and 
   
29 
 
sampled at a variable frequency that automatically matches engine speed.  The user simply 
puts in the maximum sample rate necessary in order to not miss any events along with the 
number of cycles to sample.  According to common sampling theorems, the data must be 
sampled at a frequency at least twice that of the highest frequency component of the signal of 
interest.  A sample frequency of up to ten times the highest frequency component may be 
necessary to avoid aliasing the signal [15].  The National Instruments USB-6210 data 
acquisition unit utilized to acquire the cylinder pressure and crankshaft position data was 
fully capable of providing such a sample frequency with a maximum scan rate of 250kS/s.  A 
MoTeC M400 ECU acquired all the other data pertinent to monitoring engine operation. 
 
Figure 4-6: National Instruments USB-6210 data acquisition unit 
 
 
Figure 4-7: MoTeC M400 ECU 
   
30 
 
4.4 Data Integrity and Analysis 
Daily instrumentation checks are necessary to achieve high data integrity from day-
to-day.  It is a good idea to record combustion data at the same test condition each day to 
ensure that there are no instrumentation issues [7].  Ideally, the chosen test condition should 
be fairly representative of the test points being run that day, and the combustion data at 
motored and fired conditions should always be taken in the same order.  An engine warm-up 
procedure must also be created and strictly followed.  Coolant and oil temperatures should 
also be kept constant.  The intake and exhaust restrictions should also be kept the same when 
verifying the chosen test condition each day. 
Data from the daily instrumentation checks can be used to construct a control chart.  
Typically, control charts are created for monitoring data such as IMEP, PMEP, polytropic 
coefficients, PCP, and LPP during both the firing and motoring condition check.  PCP 
measurements provide a good cylinder pressure transducer check while the LPP confirms 
encoder phasing in relation to TDC.  Data such as the MAP, AF ratios, and the heat-release 
centroid should be recorded for the firing condition check [7].  The daily test readiness and 
verification instructions that were performed before every extensive test are listed in Table 
4-2.  An example of a control chart for IMEP is shown in Figure 4-8. 
   
31 
 
 
Figure 4-8: IMEP control chart 
Table 4-2: Test readiness and verification instructions 
1 Check that pressure transducer does not drift over a ten minute period 
2 Warm up engine at idle until engine coolant temperature is at 75 ?C for five minutes 
3 Shut down engine, disconnect the cylinder #1 spark plug and injector connectors 
4 Motor engine at 3,000 rpm, 75 ?C coolant temperature  
5 Take data point 
6 Shut down engine, reconnect spark plug and injector 
7 Examine LPP (location of peak pressure) 
8 
Take engine to 3,000 rpm, 15% load, 75 ?C coolant temperature, ignition timing: 15 
deg. BTDC, rail pressure: 60 psi. injection timing: 150 degrees BTDC, lambda: 0.96 
9 Take data point 
10 Shut down engine 
11 Compare variation in pressure STD between -180? and +180? ATDC data points 
12 Examine IMEP 
 
During this experiment, there was no way to control the atmospheric pressure, 
temperature, and humidity.  The variations in atmospheric conditions are primarily 
responsible for the variations in performance.  To counter act this, any data used for 
   
32 
 
comparison on performance and efficiency plots was all recorded on the same day to avoid 
outside factors such as the weather having an influence on the results. 
The data analysis program shown in Figure 4-9 processes the data.  Since the cycle-
to-cycle variation in the pressure readings can be high, around 100 cycles were analyzed and 
averaged into a single data set to be used in the analysis.  Because the sample sizes are quite 
large, it is important to monitor the standard deviation of the pressure data.  Samples that 
display high standard deviations in the pressure data are a red flag that something is wrong.  
This allowed for the detection of the thermal shock experienced by the sensor throughout the 
experiment.  It can also indicate when the engine may not be in a steady state condition.  
Both fired and motored cylinder pressure traces were used in a series of calculations to 
provide data integrity plots, P-V and log-log P-V diagrams, and the burn rate profile.  The   
P-V diagram is integrated using the trapezoid method to provide indicated work, MEP, 
power, torque, and SFC.  The log-log P-V diagram is utilized to determine the polytropic 
exponent of compression and expansion as an error check.  If a differing value from what is 
expected is obtained, this may indicate an incorrect compression ratio or TDC location input. 
If the compression line has a curve to it, that indicates the engine has not been allowed to 
obtain a stable operating temperature and heat transfer effects are playing a role.  The log-log 
P-V also displays the intake, compression, expansion, and exhaust portions of a given cycle 
more clearly than a standard P-V diagram.  The program returns other values of interest as 
well, such as the PCP and LPP of the motored and fired pressure curves, thermal shock test, 
STD, COV and LNV of IMEP, and the SOC, combustion centroid and duration. 
   
33 
 
 
Start 
Import ASCII Data 
(motored & fired)  
Specify compression & 
expansion start & end (for 
? calculation later) 
Specify engine parameters 
(geometry & valve timing) 
Specify fluid properties 
(intake manifold, cylinder 
wall, EGT temps, Pref) 
Specify Woschni model 
parameters (swirl ratio, 
swirl velocity, Reynolds 
number, calibration values 
Allocate space for matrices 
Shift pressure traces to true 
TDC, sort into 720 degree 
cycles, eliminate double 
counted positions 
Calculate mean, median, 
mode, & standard dev. of 
pressure values 
Determine PCP, locate 
position (LPP) 
Conduct thermal shock test 
(display STD of pressure at 
--180? & +180?) 
Calculate cylinder volume 
matrix (from angular 
position) 
Integrate P-V diagram 
(averaged over 100 cycles) 
Calculate work, IMEP, 
indicated power & torque 
Construct log-log P-V 
diagram (average P?s) 
Calculate polytropic 
exponents, ?, for 
compression & expansion 
Convert pressure & volume 
data to SI units 
Integrate P-V diagram (for 
each sampled cycle) 
Calculate STD, COV, & 
LNV of IMEP 
Calculate charge and 
residual mass 
Calculate temp profile 
using ideal gas law 
Calculate mass fraction 
burned profile (Xb vs. ?) 
Calculate combustion 
centroid location
Plot Data
End 
Program Outputs: 
? Fired & motored PCP & 
LPP 
? Thermal shock check 
? Work, IMEP, indicated 
torque & power 
? STD, COV, & LNV of 
IMEP 
? Polytropic exponent 
(compression & 
expansion) 
? SOC, combustion 
centroid and duration 
? Fuel conversion 
efficiency 
? Indicated specific fuel 
consumption 
 
Program Plots: 
? Averaged pressure trace 
? Mean, mode, median of 
pressure trace 
? Standard deviation of 
pressure trace 
? Averaged P-V diagram 
? Averaged log-log P-V 
diagram 
? Mass fraction burned  
Specify rapid burn period 
(for combustion centroid 
calculation later) 
Specify operating conditions 
(rpm, MAP, ?, ect.) 
Peg pressure traces to MAP 
Establish SOC & EOC, 
calculate mass fraction 
burned with the Rassweiler 
& Withrow method (no 
partial burn detection or 
heat transfer consideration) 
Calculate ISFC
Calculate rate of work 
transfer done by the system 
Calculate Woshni heat 
transfer coefficient 
Figure 4-9: Data analysis program flowchart
   
34 
 
Chapter 5: Component Machining & Dynamometer Setup 
 
5.1 Cylinder Head Drilling Procedure 
 In most cases, the combustion chamber of an internal combustion engine can be 
accessed rather easily with either a special spark plug or glow plug which has a piezoelectric 
sensor already integrated into the plug.  However, these units can be quite expensive so the 
Kistler 6052A sensor was utilized instead.  This sensor requires that the cylinder head be 
drilled and fitted with a special sensor mounting sleeve. The sleeve allows the transducer to 
access the combustion chamber by sealing the water jacket passages that the transducer 
mounting hole must be drilled through.  It also separates the transducer from the cylinder 
head material.  This isolates the transducer from deformation stresses in the head that can 
cause a shift in sensitivity during operation.   
When mounting the pressure transducer in the head, a flush mount with minimal 
communication passage length is preferable.  Although a large communication passage 
reduces the effects of thermal shock on the transducer, it is also susceptible to errors from 
pressure waves developing in the passage.  The longer the passage, the less responsive the 
pressure sensor is to the actual cylinder pressure due to the ringing effect.  The shortest 
passageway and smallest possible cavity volume along with the largest passage diameter is 
desired to provide the highest natural frequency.  Care must be taken when drilling the 
transducer passage that no portion of the casting becomes so thin that a premature structural 
failure of the head results. 
   
35 
 
Since dimensional drawings of the Yamaha R6 cylinder head were not available, 
much planning went into the machining process.  First a pilot hole was drilled from the 
combustion chamber out to ensure that the sensor mounting sleeve would not interfere with 
the engine valves.  The hole was drilled in between the intake and exhaust valves on the 
outside wall of cylinder number one.  The pilot hole was drilled with the mill head at a 30 
degree angle and the turn table at a 5.5 degree angle.  For the second operation, a centering 
insert was spun on the lathe that once inserted into the pilot hole, could be used to find the 
center of the hole and put it in the same plane as the mill head using the tilt table setup shown 
in Figure 5-1.  With the centering insert placed inside the pilot hole, a dial indicator was used 
to adjust the top surface of the insert so that it was perpendicular with the vertical axis of the 
mill.  The counter bored hole in the center of the insert was used to center-find the pilot hole.  
The cylinder head was secured to the tilt table by a jig plate which contained three tapped 
holes for fastening the head to the jig plate through the stock head bolt holes.   
 
Figure 5-1: Tilt table and cylinder head jig adjustment on mill 
 
   
36 
 
 Figure 5-2 shows the tilt table dialed in on the mill, ready to drill.  A through hole 
was drilled first for tapping the M7 x 0.75 threads for the sensor mounting sleeve.  Then an 
end mill and 10 mm reamer were utilized to machine the counter bore for the sensor 
mounting sleeve.   
 
Figure 5-2: Mill setup for second cylinder head machining operation 
 
The mounting sleeve sealed the inner combustion chamber wall from the water jacket 
by a 0.11 inch thick custom made copper washer.  The o-ring on the mounting sleeve along 
with the aid of semi-liquid Honda Bond was used to seal the outer wall of the water jacket.  
The fitment of the mounting sleeve in the R6 head is shown in Figure 5-3 and Figure 5-4. 
   
37 
 
 
Figure 5-3: Pressure sensor communication passage 
 
 
Figure 5-4: Mounting sleeve location in cylinder head 
 
 The engine run on the dynamometer for this testing was completely torn down prior 
to testing both to install the instrumented cylinder head and to refurbish the worn and abused 
power cylinder and clutch components.  The teardown also offered an opportunity to 
Mounting Sleeve 
   
38 
 
accurately measure critical geometric details of the engine?s cylinder head and power 
cylinder components.  
 
Figure 5-5: Disassembled upper crankcase, transmission, & power cylinder components 
 
Figure 5-6: Disassembled lower crankcase, cylinder head, & clutch components 
 
 
   
39 
 
5.2 Angular Encoder Component Machining and Dynamometer Setup 
Since the stock trigger wheel did not provide enough resolution for combustion 
analysis, an industrial encoder was mounted to the engine and directly coupled with the 
crankshaft.  The BEI H25 incremental encoder provided 360 degrees of resolution for 
monitoring the crankshaft position of the engine.  The picture below shows the bracket 
constructed for mounting the encoder to the engine as well as the mating shaft (coming out 
from the crankshaft), custom oil seal housing, and the flexible bellows used to couple the 
mating and encoder shafts. 
 
Figure 5-7: Crankshaft encoder assembly 
 
 The encoder bracket had to be constructed without the CMM coordinates of the 
center of the crankshaft with respect to the crankcase cover mounting holes since the engine 
Oil Seal 
Housing 
Mating 
Shaft 
Flexible 
Bellows 
Angular 
Encoder 
Encoder 
Bracket 
   
40 
 
was not able to be measured in time by an outside source.  Therefore, the feet of the bracket 
and a rectangular plate were machined separately.  Pins for the mounting holes were spun on 
the lathe and used to hold the assembly together while being welded.  Once the bracket was 
welded together, the center of the inspection hole on the crankcase cover was used to locate 
the mounting holes for the angular encoder on the plate so that the encoder shaft would be 
concentric with the crankshaft.  After the drilling was complete, the corners of the plate were 
trimmed off to lighten the bracket assembly. 
 
Figure 5-8: Crankshaft encoder bracket construction 
  
The stock flange bolt that secures the trigger wheel to the crankshaft was replaced 
with a custom made mating shaft to provide a direct coupling to the crankshaft.  The mating 
shaft shown in Figure 5-9 was spun on the lathe out of ?? stainless steel hex stock.  All of 
the special components used for this testing were made in-house at the SAE student machine 
shop.  The complete engine dynamometer setup is shown in Figure 5-10.  The engine was 
managed by a MoTeC 400 ECU and was coupled to a DYNO-mite water brake 
dynamometer.  A heat exchanger mounted underneath the engine provided the necessary heat 
rejection from the engine coolant system. 
   
41 
 
 
Figure 5-9: Crankshaft encoder mating shaft 
 
Figure 5-10: Engine dynamometer test stand 
Trigger Wheel 
Stock Flange 
Bolt 
Mating Shaft 
?? Hex for installation 
& removal 
Flat for coupling to 
flexible bellows 
M8 x 1.25 threads 
for insertion into 
crankshaft 
   
42 
 
Chapter 6: Design of Experiment 
 
6.1 15% Load, 3,000 RPM DOE 
 A two-level factorial design of experiment (DOE) was created to investigate the burn 
rate curve?s response to certain engine operating parameters.  Ignition timing, fuel rail 
pressure, injection timing, and the relative air/fuel ratio, lambda, were the four factors chosen 
to be investigated.  The start of combustion (SOC), combustion centriod, and combustion 
duration were the three mass fraction burned curve responses that were examined.  The 
experiment was designed as an eight-run, half factorial analysis to examine both the main 
and two factor interaction effects.  Table 6-1 below lists the results of the DOE. 
Table 6-1: 15% load, 3,000 RPM DOE matrix 
15% Load, 3,000 RPM Main Effects Response 
Standard Run 
A: Ignition 
Timing 
B: Rail 
Pressure 
C: Injection 
Timing 
D: 
Lambda SOC Centroid Duration 
5 1 10 50 200 0.98 6 8.1 24 
4 2 17 70 100 0.93 12 11.3 34 
6 3 17 50 200 0.93 13 9.9 33 
1 4 10 50 100 0.93 5 10.2 25 
7 5 10 70 200 0.93 8 9.1 27 
3 6 10 70 100 0.98 8 7.9 27 
8 7 17 70 200 0.98 14 8.3 33 
2 8 17 50 100 0.98 11 8.9 30 
SOC Effect 5.75 1.75 1.25 0.25 9.6   
Centroid Effect 0.78 0.13 0.72 1.83 9.2  
Duration Effect 6.75 2.25 0.25 1.25     29.1 
  
Interaction Effects 
  
      
A*B A*C A*D       
SOC Effect 0.75 0.75 0.25   9.6     
Centroid Effect 0.53 0.28 0.17     9.2   
Duration Effect 0.25 0.75 0.75       29.1 
   
43 
 
 
 Half-normal plots were constructed to examine the significance of each factor on the 
three mass fraction burned curve parameters.  Figure 6-1 displays the effects on the location 
of SOC.  At first glance, ignition timing, fuel rail pressure, injection timing, and the 
interaction of ignition timing with injection timing and fuel rail pressure seem significant.  
However, the analysis of variance (ANOVA) shown in Table 6-2 revealed that only ignition 
timing, fuel rail pressure, and injection timing were important.  The probability value of 
0.0955 for the two interaction effects was greater than the 0.05 cutoff [18], and thus was 
considered to have possibly been a result of experimental noise.  The model F-value of 124 
implies that the model is significant.  There is only a 0.80% chance that a model F-value this 
large could occur due to noise. 
Half Normal plot
H
a
l
f
 N
o
rm
al
 %
 
p
r
oba
bi
l
i
t
y
|Effect|
0.00 1.44 2.88 4.31 5.75
0
20
40
60
70
80
85
90
95
97
99
A
B
C
AB
AC
 
Figure 6-1: 15% load, 3,000 RPM SOC half-normal plot 
A: Ignition Timing
B: Rail Pressure
C: Injection Timing
D: Lambda
   
44 
 
 
Table 6-2: 15% load, 3,000 RPM SOC ANOVA 
Source Sum of Squares DF Mean  Square F  Value Prob > F
Model 77.63 5 15.53 124.20 0.0080
A 66.13 1 66.13 529.00 0.0019
B 6.13 1 6.13 49.00 0.0198
C 3.13 1 3.13 25.00 0.0377
AB 1.13 1 1.13 9.00 0.0955
AC 1.13 1 1.13 9.00 0.0955
Residual 0.25 2 0.13  
Corrected Total 77.88 7  
 
 In the same fashion, a half-normal plot and ANOVA analysis were constructed for the 
effects on the location of the combustion centriod.  Only ignition timing, injection timing, 
and lambda were found to be significant factors in the ANOVA analysis.  The model F-value 
of 29 indicates that there is only a 0.98% chance that it could occur due to noise. 
Half Normal plot
H
a
l
f
 N
o
r
m
al
 %
 
pr
ob
ab
i
l
i
t
y
|Effect|
0.00 0.46 0.91 1.37 1.83
0
20
40
60
70
80
85
90
95
97
99
A
C
D
AB
 
Figure 6-2: 15% load, 3,000 RPM combustion centroid half-normal plot 
A: Ignition Timing
B: Rail Pressure
C: Injection Timing
D: Lambda
   
45 
 
 
Table 6-3: 15% load, 3,000 RPM combustion centriod ANOVA 
Source Sum of Squares DF Mean  Square F  Value Prob > F
Model 9.46 4 2.37 29.12 0.0098 
A 1.20 1 1.20 14.78 0.0310 
C 1.05 1.05 12.94 0.0368 
D 6.66 1 6.66 81.98 0.0028 
AB 0.55 1 0.55 6.78 0.0800 
Residual 0.24 3 0.081   
Corrected Total 9.71 7    
  
The half-normal plot and ANOVA analysis for the combustion duration yielded 
ignition timing, fuel rail pressure, and lambda as the key factors.  The model F-value of 171 
indicates that there is only a 0.58% chance that it could occur due to noise. 
Half Normal plot
H
a
lf
 N
o
rm
a
l
 %
 p
r
o
b
a
b
ilit
y
|Effect|
0.00 1.69 3.37 5.06 6.75
0
20
40
60
70
80
85
90
95
97
99
A
B
D
AC
AD
 
Figure 6-3: 15% load, 3,000 RPM combustion duration half-normal plot 
A: Ignition Timing
B: Rail Pressure
C: Injection Timing
D: Lambda
   
46 
 
Table 6-4: 15% load, 3,000 RPM combustion duration ANOVA 
Source Sum of Squares DF Mean  Square F  Value Prob > F
Model 106.63 5 21.32 170.60 0.0058 
A 91.12 1 91.12 729.00 0.0014 
B 10.12 10.12 81.00 0.0121 
D 3.13 1 3.13 25.00 0.0377 
AC 1.13 1 1.13 9.00 0.0955 
AD 1.12 1 1.12 9.00 0.0955 
Residual 0.25 2 0.12   
Corrected Total 106.88 7    
 
The average location of SOC for the eight DOE runs was 9.6 degrees BTDC.  The 
average combustion centroid was 9.2 degrees ATDC, and the average combustion duration 
was 29.1 degrees.  Figure 6-4 displays the burn rate profiles for each of the DOE runs.  
 
Figure 6-4: 15% load, 3,000 RPM DOE mass fraction burned curves 
 
At the 15% load, 3,000 RPM operating point, a more homogeneous charge from 
higher fuel rail pressures resulted in less ignition delay.  It also allowed for a more complete 
burn and higher fuel conversion efficiency.  Commanding the fuel to be injected earlier in the 
   
47 
 
cycle also provided a decreased ignition delay as well as a combustion centroid closer to 
TDC.  The mixtures burning close to stoichometric at a 0.98 relative air/fuel ratio provided a 
centroid further from TDC.   
6.2 60% Load, 3,000 RPM DOE 
Another experiment was conducted to determine the effect of the four engine 
operating parameters on the three mass fraction burned profile quantities under study at a 
higher engine load.  Table 6-5 lists the results of the DOE run at the 60% load, 3,000 RPM 
operating point. 
Table 6-5: 60% load, 3,000 RPM DOE matrix 
 
 The half-normal plot shown in Figure 6-5 displays the effects on the location of SOC.   
The ANOVA revealed that ignition timing was the main factor of importance.  The 
probability value of 0.1161 for fuel rail pressure and the interaction between fuel rail 
pressure and ignition timing indicated that those effects could have possibly been a result of 
60% Load, 3,000 RPM Main Effects Response 
Standard Run 
A: Ignition 
Timing 
B: Rail 
Pressure 
C: Injection 
Timing 
D: 
Lambda SOC Centroid Duration 
5 1 15 50 200 0.96 8 9.8 26 
4 2 22 70 100 0.87 14 9.2 30 
6 3 22 50 200 0.87 14 8.3 29 
1 4 15 50 100 0.87 7 10.9 26 
7 5 15 70 200 0.87 10 10.6 29 
3 6 15 70 100 0.96 9 12.2 31 
8 7 22 70 200 0.96 15 7.2 30 
2 8 22 50 100 0.96 15 9.8 29 
SOC Effect 6.00 1.00 0.50 0.50 11.5     
Centroid Effect 2.25 0.10 1.55 0.00 9.8  
Duration Effect 1.50 2.50 0.50 0.50     28.8 
  
Interaction Effects 
  
      
A*B A*C A*D       
SOC Effect 1.00 0.50 0.50   11.5     
Centroid Effect 0.95 0.20 0.25     9.8   
Duration Effect 1.50 0.50 0.50       28.8 
   
48 
 
experimental noise.  This varies slightly from the light load case in that the effects of fuel rail 
pressure and injection timing have been drowned out.  The model F-value of 51 gives a 
0.12% chance that the model could have occurred due to noise. 
Half Normal plot
H
a
lf
 N
o
r
m
a
l
 %
 p
r
o
b
a
b
ilit
y
|Effect|
0.00 1.50 3.00 4.50 6.00
0
20
40
60
70
80
85
90
95
97
99
A
B
AB
 
Figure 6-5: 60% load, 3,000 RPM SOC half-normal plot 
 
Table 6-6: 60% load, 3,000 RPM SOC ANOVA 
Source Sum of Squares DF Mean  Square F  Value Prob > F
Model 76.00 3 25.33 50.67 0.0012
A 72.00 1 72.00 144.00 0.0003
B 2.00 1 2.00 4.00 0.1161
AB 2.00 1 2.00 4.00 0.1161
Residual 2.0 4 0.50  
Corrected Total 78.00 7  
A: Ignition Timing
B: Rail Pressure
C: Injection Timing
D: Lambda
   
49 
 
 
 
 Ignition timing, injection timing, and the interaction between ignition timing and fuel 
rail pressure were found to be significant factors in the placement of the combustion centroid.  
At the richer lambda limits, lambda was not found to be a significant factor as in the light 
load case (where lambda ranged from 0.93 to near stoichometric at 0.98).  This difference 
will be discussed later in the chapter.  The model F-value of 99 indicates that there is only a 
0.03% chance that it could occur due to noise. 
Half Normal plot
H
a
l
f
 N
o
r
m
al
 %
 pr
obabi
l
i
t
y
|Effect|
0.00 0.56 1.13 1.69 2.25
0
20
40
60
70
80
85
90
95
97
99
A
C
AB
 
Figure 6-6: 60% load, 3,000 RPM combustion centroid half-normal plot 
A: Ignition Timing
B: Rail Pressure
C: Injection Timing
D: Lambda
   
50 
 
Table 6-7: 60% load, 3,000 RPM combustion centriod ANOVA 
Source Sum of Squares DF Mean  Square F  Value Prob > F
Model 16.74 3 5.58 99.17 0.0003 
A 10.13 1 10.13 180.00 0.0002 
C 4.81 1 4.81 85.42 0.0008 
AB 1.80 1 1.80 32.09 0.0048 
Residual 0.22 4 0.056   
Corrected Total 16.96 7    
  
The half-normal plot and ANOVA analysis for the combustion duration yielded 
ignition timing, fuel rail pressure, and the interaction between ignition timing and fuel rail 
pressure as the key factors.  The model F-value of 14 indicates that there is only a 1.32% 
chance that it could occur due to noise. 
Half Normal plot
H
a
l
f
 N
o
r
m
al
 %
 pr
obabi
l
i
t
y
|Effect|
0.00 0.62 1.25 1.87 2.50
0
20
40
60
70
80
85
90
95
97
99
A
B
AB
 
Figure 6-7: 60% load, 3,000 RPM combustion duration half-normal plot 
A: Ignition Timing
B: Rail Pressure
C: Injection Timing
D: Lambda
   
51 
 
Table 6-8: 60% load, 3,000 RPM combustion duration ANOVA 
Source Sum of Squares DF Mean  Square F  Value Prob > F
Model 21.50 3 7.17 14.33 0.0132 
A 4.50 1 4.50 9.00 0.0399 
B 12.50 12.50 25.00 0.0075 
AB 4.50 1 4.50 9.00 0.0399 
Residual 2.00 4 0.50   
Corrected Total 23.50 7    
 
The average location of SOC for the eight DOE runs was 11.5 degrees BTDC.  The 
average combustion centroid was 9.8 degrees ATDC, and the average combustion duration 
was 28.8 degrees.  Figure 6-8 displays the burn rate profiles for each of the DOE runs for the 
60% load, 3.000 RPM operating point.  
 
Figure 6-8: 60% load, 3,000 RPM DOE mass fraction burned curves 
 
At the 60% load, 3,000 RPM operating point, a more homogeneous charge from 
higher fuel rail pressures allowed for a more complete burn and higher fuel conversion 
efficiency.  Commanding the fuel to be injected earlier in the cycle provided a faster burn 
   
52 
 
and thus a centroid closer to TDC.  Mixtures burning closer to stoichometric at a 0.96 relative 
air/fuel ratio exhibited a combustion centroid further from TDC.   
The influence of ignition timing on the mass fraction burned curve is illustrated in 
Figure 6-9.  In the plot, the ignition timing of 20 degrees BTDC corresponds to MBT timing.  
This is where the steepest rapid burn angle is seen.  If the timing is advanced or retarded 
slightly from MBT, the rapid burn angle is reduced and a larger amount of the mixture goes 
unburned.  If the timing is retarded substantially, a substantial incomplete burn is seen.  If the 
timing is advanced even further from what is shown, the engine may experience knock. 
 
Figure 6-9: Mass fraction burned dependency on ignition timing 
 
 The influence of the excess air ratio on the mass fraction burned profile was clearly 
identified in the 15% load, 3,000 RPM DOE.  However, in the 60% load case, the effects of 
lambda were not as clear.  Further testing at this operating point was performed to construct 
Figure 6-10 which confirms that lambda does indeed influence the curve at the higher load 
spark advance 
   
53 
 
operating point.  A change in lambda did not affect the rapid burn angle, but did shift the 
combustion centriod further from TDC as well as extend the duration slightly. 
 
Figure 6-10: Mass fraction burned dependency on lambda 
 
 Similar mass fraction burned curve deformations due to independent changes in 
ignition timing and excess air ratio have been found by V?vra and Tak?ts [19] in their study 
of a SI turbocharged ML636NG six cylinder engine.  These general trends discovered by the 
DOE and subsequent testing were utilized to uncover optimum combustion processes at both 
operating points under examination. 
leaner mixture 
   
54 
 
Chapter 7: Experimental Results 
 
7.1 15% Load, 3,000 RPM Experimental Results 
 The DOE was constructed with limits such that the optimum engine operating points 
for best ISFC and IMEP would be bracketed.  The statistical results were used to understand 
how to manipulate the mass fraction burned curve to improve engine performance and 
efficiency.  Figure 7-1 and Figure 7-2 show the ISFC and IMEP of the engine versus the start 
of combustion (SOC) location.  Best ISFC occurs with a SOC of 9 degrees BTDC while 
maximum IMEP occurs when SOC is retarded one more degree to 8 degrees BTDC.   
 
Figure 7-1: 15% load, 3,000 RPM ISFC as a function of SOC placement 
   
55 
 
 
 
Figure 7-2: 15% load, 3,000 RPM IMEP as a function of SOC placement 
 The relationship between the location of SOC and ignition timing, 
t
i  (in terms of 
crank angle BTDC), is shown in the ignition delay plot below.  Deviation from the curve fit 
is caused by changes in other factors like fuel rail pressure, injection timing, and excess air 
ratio since these factors were not held constant for each data point.  Thus, a general equation 
to estimate the ignition delay at this engine speed and load can be obtained.  
 
Figure 7-3: 15% load, 3,000 RPM ignition delay 
25.1 234.0 ??
t
iDelayIgnition  
   
56 
 
 The effect of the combustion centroid placement on ISFC and IMEP can be seen in 
Figure 7-4 and Figure 7-5.  Best ISFC occurs with a combustion centroid of 8.7 degrees 
ATDC while maximum IMEP occurs when the combustion centroid is located at 9.1 degrees 
ATDC.  Since ISFC is in terms of the fuel flow rate per unit power output, it is also useful to 
compare the shear amount of fuel being injected per cycle.  The optimum ISFC point 
consumes 8.5% less fuel per cycle than the maximum IMEP point.  There is always some 
kind of tradeoff to make between torque and fuel economy throughout the engine map.  
However, under this light load, low speed engine operating point, it is highly advantageous to 
choose to run at the minimum ISFC point and enjoy the maximum fuel conversion 
efficiency. 
 
Figure 7-4: 15% load, 3,000 RPM ISFC as a function of combustion centroid placement 
   
57 
 
 
Figure 7-5: 15% load, 3,000 RPM IMEP as a function of combustion centriod placement 
 
 Figure 7-6 displays the deformation of the mass fraction burned curve due to changes 
in the combustion centroid location.  Both the closest and furthest centroids from TDC as 
well as the optimum centroid from Figure 7-4 and Figure 7-5 are shown. 
 
Figure 7-6: 15% load, 3,000 RPM combustion centroid placement characteristics 
   
58 
 
 
 The mass fraction burned curve exhibiting the optimum centriod location boasts a 
steeper rapid burn rate and more complete burn over the other two curves.  The tabulation 
below displays the Wiebe function curve fit parameters and performance and efficiency 
comparisons between the three curves. 
Table 7-1: 15% load, 3,000 RPM performance and efficiency comparison 
Combustion 
Centriod 
Max Mass 
Burn Rate 
Efficiency 
Parameter, a 
Form 
Factor, m 
ISFC IMEP 
7.9 deg. ATDC 0.033 3.30 1.10 0.631 34.3 
8.7 deg. ATDC 0.039 3.60 1.05 0.469 43.2 
11.2 deg. ATDC 0.031 2.75 0.70 0.720 31.9 
 
 The fired and motored averaged pressure traces for the optimum engine tune are 
shown in Figure 7-7.  For the fired cycle, peak cylinder pressure was 265 psi.  SOC occurred 
at 9 degrees BTDC.   The largest standard deviation in cylinder pressure readings for the 
fired cycle occurred around 60 degrees ATDC where it peaked at 33 psi due to thermal 
shock.  The standard deviations seen during the exhaust and intake events were minimal. 
-400 -300 -200 -100 0 100 200 300 400
0
50
100
150
200
250
300
Crankshaft Angle (Degrees)
C
y
l
i
nder P
r
es
s
u
r
e
 (ps
i
)
Pressure vs. Crankshaft Position
 
 
Fired Cycle
Motored Cycle
 
Figure 7-7: 15% load, 3,000 RPM fired and motored cycle 
 
 Figure 7-8 displays the P-V diagram for the pressure trace shown in Figure 7-7.  
Integration of the diagram provided 32.9 ft-lbf of indicated work per cycle per cylinder and 
   
59 
 
an IMEP of 43.2 psi.  The indicated torque was 10.5 ft-lbf, and the indicated power at 3,000 
RPM was 6.0 hp for the four cylinders combined. 
0 1 2 3 4 5 6 7 8 9 10
0
50
100
150
200
250
300
Cylinder Volume (cubic inches)
C
y
l
i
nder
 P
r
es
s
u
r
e
 (
p
s
i
)
P-V Diagram
 
Figure 7-8: 15% load, 3,000 RPM P-V diagram 
 The log-log P-V diagram revealed a polytropic exponent of 1.30 on the compression 
stroke which matches the calculated ratio of specific heats.  The line does not exhibit any 
curvature which implies that the engine had been allowed to reach a stable operating 
temperature at the time of the data acquisition. 
-0.2 0 0.2 0.4 0.6 0.8 1 1.2
0.8
1
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
Log of Cylinder Volume (cubic inches)
Log
 of
 C
y
l
i
n
der P
r
es
s
u
r
e
 (
p
s
i
)
Log-Log P-V Diagram
 
Figure 7-9: 15% load, 3,000 RPM log-log P-V diagram 
 
   
60 
 
 The mass fraction burned curve derived from the in-cylinder pressure trace offered 
the desired SOC of 9 degrees BTDC and combustion centroid of 8.7 degrees ATDC.  This 
resulted in a fast burn rate in the rapid burn portion of the curve as well as a fairly complete 
burn which translated into an improved fuel conversion efficiency. 
-15 -10 -5 0 5 10 15 20 25 30 35
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Crankshaft Angle (Degrees)
M
a
s
s
 F
r
ac
t
i
on
 B
u
rn
ed
,
 X
b
Mass Fraction Burned vs. Crankshaft Position
 
 
Experimental Data
Wiebe Function
 
Figure 7-10: 15% load, 3,000 RPM mass fraction burned curve 
 
7.2 60% Load, 3,000 RPM Experimental Results 
In the same fashion, ISFC and IMEP as a function of SOC and combustion centroid 
location were examined for the 60% load, 3,000 RPM case.  Figure 7-11 and Figure 7-12 
show the ISFC and IMEP of the engine versus the SOC location.   The relationship between 
the location of SOC and ignition timing, 
t
i  (in terms of crank angle BTDC), is shown in 
Figure 7-13.   
   
61 
 
 
 
Figure 7-11: 60% load, 3,000 RPM ISFC as a function of SOC placement 
 
 
Figure 7-12: 60% load, 3,000 RPM IMEP as a function of SOC placement 
   
62 
 
 
Figure 7-13: 60% load, 3,000 RPM ignition delay 
 The effect of the combustion centroid placement on ISFC and IMEP can be seen in 
Figure 7-14 and Figure 7-15.  Best ISFC and IMEP occur with a combustion centroid of 8.8 
degrees ATDC.   
 
Figure 7-14: 60% load, 3,000 RPM ISFC as a function of combustion centroid placement 
  
08.1 299.0 ??
t
iDelayIgnition  
   
63 
 
 
Figure 7-15: 60% load, 3,000 RPM IMEP as a function of combustion centriod placement 
 
Figure 7-16 displays the deformation of the mass fraction burned curve due to 
changes in the combustion centroid location.  Both the closest and furthest centroids from 
TDC as well as the optimum centroid from Figure 7-14 and Figure 7-15 are shown. 
 
Figure 7-16: 60% load, 3,000 RPM combustion centriod placement characteristics 
   
64 
 
 
 The mass fraction burned curve exhibiting the optimum centriod location again has a 
steeper rapid burn angle and more complete burn over the other two curves.  The tabulation 
below displays the Wiebe function curve fit parameters and performance and efficiency 
comparisons between the three curves. 
Table 7-2: 60% load, 3,000 RPM performance and efficiency comparison 
Combustion 
Centriod 
Max Mass 
Burn Rate 
Efficiency 
Parameter, a 
Form 
Factor, m 
ISFC IMEP 
8.3 deg. ATDC 0.034 3.10 1.20 0.423 77.5 
8.8 deg. ATDC 0.038 3.30 1.40 0.356 87.5 
10.9 deg. ATDC 0.032 2.60 0.90 0.493 68.6 
  
The fired and motored averaged pressure traces for the optimum engine tune are 
shown in Figure 7-17.  For the fired cycle, peak cylinder pressure was 337 psi.  SOC 
occurred at 13 degrees BTDC.   The largest standard deviation in cylinder pressure readings 
for the fired cycle occurred around 40 degrees ATDC where it peaked at 13 psi due to 
thermal shock.  The standard deviations seen during the exhaust and intake events were 
minimal. 
-400 -300 -200 -100 0 100 200 300 400
0
50
100
150
200
250
300
350
Crankshaft Angle (Degrees)
C
y
l
i
n
de
r
 P
r
es
s
u
r
e
 
(
p
s
i
)
Pressure vs. Crankshaft Position
 
 
Fired Cycle
Motored Cycle
 
Figure 7-17: 60% load, 3,000 RPM fired and motored cycle 
   
65 
 
 
 Figure 7-18 displays the P-V diagram for the pressure trace shown in Figure 7-17.  
Integration of the diagram provided 66.7 ft-lbf of indicated work per cycle per cylinder and 
an IMEP of 87.5 psi.  The indicated torque was 21.2 ft-lbf, and the indicated power at 3,000 
RPM was 11.5 hp for the four cylinders combined. 
0 1 2 3 4 5 6 7 8 9 10
0
50
100
150
200
250
300
350
Cylinder Volume (cubic inches)
C
y
l
i
nd
er P
r
e
s
s
u
re
 (ps
i
)
P-V Diagram
 
Figure 7-18: 60% load, 3,000 RPM P-V diagram 
 The log-log P-V diagram showed that thermal shock played more of a role during the 
expansion event than in any of the light load cases (more curvature).  However, the 
compression line does not exhibit any curvature which implies that the engine had been 
allowed to reach a stable operating temperature at the time of the data acquisition. The 
polytropic exponent of 1.28 on the compression stroke matches the calculated ratio of 
specific heats. 
   
66 
 
-0.2 0 0.2 0.4 0.6 0.8 1 1.2
1
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
Log of Cylinder Volume (cubic inches)
Lo
g of
 C
y
l
i
n
der
 P
r
es
s
u
r
e
 (
p
s
i
)
Log-Log P-V Diagram
 
Figure 7-19: 60% load, 3,000 RPM log-log P-V diagram 
 
 The mass fraction burned curve derived from the in-cylinder pressure trace offered 
the desired SOC of 13 degrees BTDC and combustion centroid of 8.8 degrees ATDC.  This 
resulted in a fast burn rate in the rapid burn portion of the curve as well as a more complete 
burn which resulted in an improved fuel conversion efficiency and mean effective pressure. 
-15 -10 -5 0 5 10 15 20 25 30 35
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Crankshaft Angle (Degrees)
M
a
s
s
 F
r
ac
t
i
on B
u
r
ned,
 X
b
Mass Fraction Burned vs. Crankshaft Position
 
 
Experimental Data
Wiebe Function
 
Figure 7-20: 60% load, 3,000 RPM mass fraction burned curve 
 
   
67 
 
 The 15% load, 3,000 RPM study provided an ISFC of 0.469 lbm/hp-hr and an IMEP 
of 43.2 psi, while the 60% load, 3,000 RPM study provided an ISFC of 0.356 lbm/hp-hr and 
an IMEP of 87.5 psi.  The optimum 15% and 60% load operating points provided mass burn 
rates of 0.039 and 0.038 respectively.  These values are comparable to the maximum mass 
burn rate of 0.041 found by Attard, Konidaris, Hamori, Toulson, and Watson [20] for their 
430 cc parallel twin PFI WATTARD engine operating at 10,000 RPM, WOT, MBT spark 
timing. 
   
68 
 
Chapter 8: Conclusions 
 
 The in-cylinder pressure analysis performed in this experiment proved to be more 
capable of improving engine performance and efficiency over tuning techniques previously 
used by the Auburn FSAE team.  However, this type of analysis is quite computationally 
expensive as the data sets are large, and the analysis of the data sets must be done with care 
to ensure quality in the results.  Few problems were experienced with the instrumentation 
with the exception of the thermal shock which continually occurred late in the expansion 
stroke.  With a dry transducer, the only thing that can be done to reduce thermal shock errors 
is to increase the length of the communication passage between the pressure sensor and 
combustion chamber.  This may help reduce thermal shock, but may also introduce other 
problems like ringing effects from resonance within the longer communication passage. 
 For the 15% load, 3,000 RPM case, an ISFC of 0.469 lbm/hp-hr and an IMEP of 43 
psi were realized.  This was done by adjusting the ignition timing to 14 degrees BTDC, fuel 
rail pressure to 70 psi, injection timing to 125 degrees BTDC, and relative air/fuel ratio to 
0.97.  For the 60% load, 3,000 RPM case, an ISFC of 0.356 lbm/hp-hr and an IMEP of 88 psi 
were realized.  This was achieved by altering the ignition timing to 20 degrees BTDC, fuel 
rail pressure to 70 psi, injection timing to 175 degrees BTDC, and relative air/fuel ratio to 
0.95.  The adjusted settings for both cases resulted in a faster burn rates in the rapid burn 
portion of the curve as well as a more complete burn which translated into an improved fuel 
conversion efficiency and mean effective pressure. Wiebe functions fitted to the experimental 
   
69 
 
data utilized efficiency parameters of 3.60 and 3.30 and form factors of 1.05 and 1.40 
respectively.  Both studies shared a similar optimum combustion centroid at 8.7 and 8.8 
degrees ATDC and a maximum mass burn rate of 0.039 and 0.038 correspondingly, showing 
good agreement between the two tests.   
   
70 
 
References 
 
 
1. Brunt, M., Rai, H., and Emtage, A., ?The Calculation of Heat Release Energy from 
Engine Cylinder Pressure Data,? SAE paper 981052, 1998 
 
2. Plint, M., and Martyr, A., Engine Testing Theory and Practice, 1995, Woburn, MA: 
Butterworth-Heinemann  
 
3. Woschni, G., ?A Universally Applicable Equation for the Instantaneous Heat Transfer 
Coefficient in the Internal Combustion Engine,? SAE paper 670931, 1967 
 
4. Stone, R., Introduction to Internal Combustion Engines, 1992, Warrendale, PA: SAE, 
Inc. 
 
5. Iorio, B., Giglio V., Police G., and Rispoli, N., ?Methods of Pressure Cycle Processing 
for Engine Control,? SAE paper 2003-01-0352, 2003 
 
6. Heywood, J.B., Internal Combustion Engine Fundamentals, 1988, New York, NY: 
McGraw-Hill, Inc. 
 
7. Atkins, R., An Introduction to Engine Testing and Development, 2009, Warrendale, PA, 
SAE International 
 
8. Ferguson, C.R., Internal Combustion Engines ? Applied Thermosciences, 1986, Toronto, 
Canada: John Wiley & Sons, Inc. 
 
9. Rassweiler, G.M., and Withrow, L., ?Motion Pictures of Engine Flames Correlated with 
Pressure Cards,? SAE paper 670931, 1967 
 
10. Gatowski, J.A., Balles, E.N., Chun, K.M., Nelson, F.E., Ekchian, J.A., and Heywood, 
J.B., ?Heat Release Analysis of Engine Pressure Data,? SAE paper 841359, 1984 
 
11. Cheung, H.M., Heywood, J.B., ?Evaluation of a One-Zone Burn-Rate Analysis Procedure 
Using Production SI Engine Pressure Data,? SAE paper 932749, 1993 
 
12. Cengel, Y. and Boles, M., Thermodynamics An Engineering Approach, 2008, New York, 
NY: McGraw-Hill, Inc. 
   
71 
 
 
13. Eriksson, L., ?Requirements for and a Systematic Method for Identifying Heat-Release 
Model Parameters,? SAE paper 980626, 1998 
 
14. Jones, P., MECH 4410 Engines, Summer 2010 course notes, Auburn, AL, Auburn 
University, Department of Mechanical Engineering 
 
15. Rogers, D., Engine Combustion: Pressure Measurement and Analysis, 2010, Warrendale, 
PA: SAE, International 
 
16. Kistler Instruction Manual, 2003, Winterthur, Switzerland, Kistler Instrument 
Corporation 
 
17. ?The Piezoelectric Effect, Theory, Design, and Usage?, 
http://www.designinfo.com/kistler/ref/tech_theory_text.htm, Winterthur, Switzerland, 
Kistler Instrument Corporation 
 
18. Anderson, M. and Whitcomb, P., DOE Simplified: Practical Tools for Effective 
Experimentation, 2000, Portland, OR: Productivity, Inc. 
 
19. V?vra, J. and Tak?ts, M., ?Heat Release Regression Model for Gas Fuelled SI Engines,? 
SAE paper 2004-01-1462, 2004 
 
20. Attard, W., Konidaris, S., Hamori, F., Toulson, E., Watson, H., ?Compression Ratio 
Effects on Performance, Efficiency, Emissions and Combustion in a Carbureted and PFI 
Small Engine,? SAE paper 2007-01-3623, 2007 
 
 
 
 
 
 
 
   
72 
 
Appendix 
 
LabView 2009 Data Acquisition Block Diagram 
 
   
73 
 
LabView 2009 Data Acquisition Front Panel 
 
 
 
 
 
 
 
   
74 
 
MATLAB R2009a Data Analysis Code 
 
%% R6_Peg_Shift_Awesome.m 
% Performs post processing of text file output from LabView 
% NOTE: first import fired pressure data file into MATLAB 
% 1.  Go to File\Import Data... 
% 2.  Set delimiting type (tab) 
% 3.  Set "Number of text header lines" to 23 for typical file 
% 4.  Ensure that the first values of the file are included in the preview  
%     (i.e. crank=1 deg.), if not, adjust the "Number of text header 
lines" to the appropriate value 
% 5.  Click "Next"; check "create variables matching preview" 
% 6.  Remove checks from everything BUT the data matrix; click "Finish" 
% 7.  Repeat steps 1-6 for the motored pressure data file, then right-
click & rename the data matrix as "data2"   
% 8.  Edit pressure and position data column callouts below to obtain data 
in the sorting loop (if necessary) 
% 9.  Enter engine speed, number of cycles sampled, & lambda for the test 
% 10. Enter any other appropriate values (i.e. Tim, Twall, EGT, mfdot, 
%     Woschini calibration parameters, etc.) 
% 11. Adjust any other parameters as necessary 
  
clc; 
  
% global statements 
global Ru fuel phi EGT Tim 
  
% Data column callouts (for both fired & motored) 
Pcol=2;  % Pressure data column 
Ccol=3;  % Position data column 
  
% Enter combustion analysis start and end values [degrees from TDC] 
polystrt1=-91;  % start for analyzing polytropic index 
polyend1=-28;  % end for analyzing polytropic index 
polystrt2=(-1*polyend1)+3;  % start for analyzing polytropic index 
polyend2=(-1*polystrt1)-39;  % end for analyzing polytropic index 
  
% tolerance for SOC & EOC detection  
tolerance=0.005; 
% duration for combustion centroid calculation 
rapdbrnstrt=6;   
rapdbrnend=13; 
combstrt720=349; 
combend720=378; 
  
% Enter engine speed 
rpm=3150; 
  
% Enter MAP, plenim temperature, & volumetric efficiency 
MAP=89;  % [kPa] 
Trefc=20;  % intake reference temperature [C] 
etav=0.33;  % volumetric efficiency 
  
 
   
75 
 
% Enter # of cycles sampled 
cycles=99; 
mcycles=99; 
  
% Enter lambda value 
lambda=0.96; 
  
% Enter angular encoder shift 
shift=166;  % [degrees] 
mshift=166;  % [degrees] 
  
% Engine Geometry 
bore=6.7;  % [cm] 
stroke=4.25;  % [cm] 
Rc=13.1;  % compression ratio 
numcyl=4;  % number of cylinders 
Vd=(pi*((bore/2)^2)*stroke);  % displaced volume [cc] 
Vdm=Vd/(100^3);  % displaced volume [m^3] 
Vc=Vd/(Rc-1);  % clearance volume [cc] 
Vcm=Vc/(100^3);  % clearance volume [m^3] 
  
% Valve Timing (adjusted to Theta matrix "1-720 degrees") 
IVO=721-39; 
IVC=181+65; 
EVO=541-64; 
EVC=721-24; 
  
% Fluid Properties 
Timf=164;  % intake manifold temperature [F] 
Twallf=500;  % cylinder wall temperature [F] 
EGTf=1200;  % exhaust gas temperature [F] 
Tim=((Timf-32)/1.8)+273;  % intake manifold temperature [K] 
Twall=((Twallf-32)/1.8)+273;  % cylinder wall temperature [K] 
EGT=((EGTf-32)/1.8)+273;  % exhaust gas temperature [K] 
fuel=2;  % fuel type (1-CH4; 2-C7H17; 3-C14.4H24.9; 4-CH3OH; 5-CH3NO2) 
Ru=8314.34; % universal gas constant [J/kmolK] 
Pref=MAP*1000;  % intake reference pressure [Pa] 
Tref=Trefc+273;  % intake reference temperature [K] 
phi=1/lambda;  % fuel/air equivalency ratio 
  
% Woschini's correlation for engine heat transfer 
Rs=2.5;  % swirl ratio 
N=rpm/60;  % angular velocity of crank [rot/s] 
strokem=stroke/100;  % converting to meters 
borem=bore/100;  % converting to meters 
ws=Rs*2*pi*N;  % angular velocity of fluid [rad/s] 
vs=(borem*ws)/2;  % swirl velocity term [m/s] 
Sp=2*strokem*(rpm/60);  % mean pistion speed [m/s] 
c11=6.18+(0.417*(vs/Sp));  % calibration parameter (gas exchange period) 
c12=2.28+(0.308*(vs/Sp));  % calibration parameter (rest of cycle) 
c2=3.24e-03;  % calibration parameter (only valid for combustion & 
expansion period) 
Rem=0.5;  % Reynold's number exponent 
cal1=1;  % Woschini convection coefficient calibration parameter 
cal2=1;  % Woschini convection coefficient calibration parameter 
   
76 
 
  
% Allocating 23for space 
Pm=zeros(720,cycles); 
Pdata=zeros(720,cycles); 
Cdata=zeros(720,cycles); 
Mdata=zeros(720,mcycles); 
MCdata=zeros(720,mcycles); 
Pmeas=zeros(720,cycles); 
RunPm=zeros(720,cycles); 
RunTm=zeros(720,cycles); 
Pmot=zeros(720,mcycles); 
Pmeasmot=zeros(720,mcycles); 
RunPmot=zeros(720,mcycles); 
RunTmot=zeros(720,mcycles); 
Diff=zeros(5,1); 
Diffmot=zeros(5,1); 
Diffavg=zeros(5,1); 
Diffmotavg=zeros(5,1); 
allval1=zeros(1,cycles); 
allval2=zeros(1,mcycles); 
Pavg=zeros(720,1); 
Pavgmot=zeros(720,1); 
Pmode=zeros(720,1); 
Pmedian=zeros(720,1); 
Pstd=zeros(720,1); 
Theta=zeros(720,1); 
Vol=zeros(720,1); 
fpos=zeros(180,1); 
fneg=zeros(179,1); 
wpos=zeros(179,1); 
wneg=zeros(178,1); 
fposR=zeros(180,cycles); 
fnegR=zeros(179,cycles); 
wposR=zeros(179,cycles); 
wnegR=zeros(178,cycles); 
MEPR=zeros(1,cycles); 
Pslp1=zeros(abs(polystrt1)-abs(polyend1)+1,1); 
Vslp1=zeros(abs(polystrt1)-abs(polyend1)+1,1); 
Pslp2=zeros(abs(polystrt2)-abs(polyend2)+1,1); 
Vslp2=zeros(abs(polystrt2)-abs(polyend2)+1,1); 
T=zeros(720,1); 
w=zeros(720,1); 
h=zeros(720,1); 
Acyl=zeros(720,1); 
dQdT=zeros(720,1); 
PdVdt=zeros(720,1); 
VdPdt=zeros(720,1); 
  
% Shift pressure values to TDC 
j=1; 
for i=1:1:shift 
    Pdata(i,1)=data(i,Pcol); 
    Cdata(i,1)=data(i,Ccol)-shift; 
    j=j+1; 
end 
j=1; 
   
77 
 
for i=1:1:mshift 
    Mdata(i,1)=data2(i,Pcol); 
    MCdata(i,1)=data2(i,Ccol)-mshift; 
    j=j+1; 
end 
for i=shift+1:1:720 
    Pdata(i,1)=data(i,Pcol); 
    if j<=360 
    Cdata(i,1)=data(i,Ccol); 
    else 
    Cdata(i,1)=data(i,Ccol)+360;  
    end 
end 
for i=mshift+1:1:720 
    Mdata(i,1)=data2(i,Pcol); 
    if j<=360 
    MCdata(i,1)=data2(i,Ccol); 
    else 
    MCdata(i,1)=data2(i,Ccol)+360; 
    end 
end 
% Calculate # of double crank positions for fired cycle 
fend=(720*cycles); 
    for i=721:1:fend 
        if data(i,Ccol)==data(i-1,Ccol) 
            fend=fend+1; 
        end 
    end 
% Sort reamaining data values for fired cycle 
j=1; 
k=2; 
    for i=721:1:fend 
        if data(i,Ccol)>data(i-1,Ccol) 
            Pdata(j,k)=data(i,Pcol); 
            if j<=360 
            Cdata(j,k)=data(i,Ccol); 
            else 
            Cdata(j,k)=data(i,Ccol)+360;  
            end 
                if j<720 
                j=j+1; 
                else 
                j=1; 
                k=k+1; 
                end 
        elseif data(i,Ccol)==1 
            Pdata(j,k)=data(i,Pcol); 
            if j<=360 
            Cdata(j,k)=data(i,Ccol); 
            end 
                if j<720 
                j=j+1; 
                else 
                j=1; 
                k=k+1; 
                end 
        end 
   
78 
 
    end 
% Calculate # of double crank positions for motored cycle 
mend=(720*mcycles); 
    for i=721:1:mend 
        if data2(i,Ccol)==data2(i-1,Ccol) 
            mend=mend+1; 
        end 
    end 
% Sort reamaining data values for motored cycle 
j=1; 
k=2; 
    for i=721:1:mend 
        if data2(i,Ccol)>data2(i-1,Ccol) 
            Mdata(j,k)=data2(i,Pcol); 
            if j<=360 
            MCdata(j,k)=data2(i,Ccol); 
            else 
            MCdata(j,k)=data2(i,Ccol)+360;  
            end 
                if j<720 
                j=j+1; 
                else 
                j=1; 
                k=k+1; 
                end 
        elseif data2(i,Ccol)==1 
            Mdata(j,k)=data2(i,Pcol); 
            if j<=360 
            MCdata(j,k)=data2(i,Ccol); 
            end 
                if j<720 
                j=j+1; 
                else 
                j=1; 
                k=k+1; 
                end 
        end 
    end 
  
%Pegging Pressure Traces 
for i=1:1:cycles 
    for j=178:1:182 
        Diff(j-177,i)=Pdata(j,i)-(MAP*(14.696/101.325)); 
    end 
end 
  
for i=1:1:cycles 
    for j=1:1:5 
        Diffavg(j,i)=mean(Diff(j,i)); 
    end 
end 
  
for i=1:1:mcycles 
    for j=178:1:182 
        Diffmot(j-177,i)=Mdata(j,i)-(MAP*(14.696/101.325)); 
    end 
   
79 
 
end 
  
for i=1:1:mcycles 
    for j=1:1:5 
        Diffmotavg(j,i)=mean(Diffmot(j,i)); 
    end 
end 
  
Peg=mean(Diffavg); 
Peg2=mean(Peg); 
Pegmot=mean(Diffmotavg); 
Pegmot2=mean(Pegmot); 
  
for i=1:1:cycles 
    for j=1:1:720 
        Pm(j,i)=Pdata(j,i)-Peg2; 
    end 
end  
  
for i=1:1:mcycles 
    for j=1:1:720 
        Pmot(j,i)=Mdata(j,i)-Pegmot2; 
    end 
end  
  
% Calculating mean, median, mode, & std 
for k=1:1:720 
    for m=1:1:cycles 
    allval1(1,m)=Pm(k,m); 
    end 
    for j=1:1:mcycles 
    allval2(1,j)=Pmot(k,j); 
    end 
    Pavg(k,1)=mean(allval1); 
    Pavgmot(k,1)=mean(allval2); 
    Pmode(k,1)=mode(allval1); 
    Pmedian(k,1)=median(allval1); 
    Pstd(k,1)=std(allval1); 
end 
  
% Calculating PCP and LPP 
Peak=0; 
Peakmot=0; 
for k=1:1:720 
    if Pavg(k,1)>Peak 
        Peak=Pavg(k,1); 
        Locatefire=k-1; 
    end 
    if Pavgmot(k,1)>Peakmot 
        Peakmot=Pavgmot(k,1); 
        Locatemot=k-1; 
    end 
end 
  
Locatecorr=Locatefire-360; 
   
80 
 
Locatemotcorr=Locatemot-360; 
  
fprintf('\n') 
fprintf('Fired PCP                     %6.1f psi \n',Peak)  % fired PCP   
fprintf('Fired LPP                     %6.0f degrees \n',Locatecorr)  % 
fired LPP  
fprintf('Motored PCP                   %6.1f psi \n',Peakmot)  % motored 
PCP    
fprintf('Motored LPP                   %6.0f degrees \n',Locatemotcorr)  % 
motored LPP  
  
Pstd180p=Pstd(180,1);  % standard deviation of pressure @ -180 degrees 
Pstd180m=Pstd(540,1);  % standard deviation of pressure @ +180 degrees 
  
fprintf('\n') 
fprintf('STD of Pressure -180 degrees  %6.1f psi \n',Pstd180p)  % standard 
deviation of pressure @ -180 degrees   
fprintf('STD of Pressure +180 degrees  %6.1f psi \n',Pstd180m)  % standard 
deviation of pressure @ +180 degrees 
  
% Creating theta & volume matrices 
r=stroke/2;  % crankshaft throw [cm] 
l=9.1;  % connecting rod length [cm] 
for i=-360:1:359 
    q=i+361; 
    Theta(q,1)=i; 
    s(q,1)=(r*cos(Theta(q,1)*(pi/180)))+(((l^2)-
((r^2)*((sin(Theta(q,1)*(pi/180)))^2)))^0.5); 
    Ach=(pi*((stroke/2)^2))+(2*pi*(stroke/2)*0.087785); 
    Ap=(pi*((stroke/2)^2)); 
    A(q,1)=Ach+Ap+(pi*bore*(l+r-s(q,1))); 
    Vol(q,1)=(Vc+(((pi*(bore^2))/4)*(l+r-s(q,1))))*(1/(2.54^3)); 
end  
  
% Integrating P-V diagram (averaged data) 
j=1; 
for i=1:1:180 
    fpos(j,1)=abs(Pavg(i+361,1)-Pavg(361-i,1)); 
    j=j+1; 
end 
j=1; 
for i=1:1:179 
    wpos(j,1)=abs((Vol(i+361,1)-
Vol(i+362,1))*((fpos(i,1)+fpos(i+1,1))/2)); 
    j=j+1; 
end 
j=1; 
for i=1:1:179 
    fneg(j,1)=-1*(abs(Pavg(181-i,1)-Pavg(541+i,1))); 
    j=j+1; 
end 
j=1; 
for i=1:1:178 
    wneg(j,1)=-1*abs((Vol(i+179,1)-
Vol(i+178,1))*((fneg(i,1)+fneg(i+1,1))/2)); 
    j=j+1; 
   
81 
 
end 
W=0; 
for i=1:1:178 
    W=W+wpos(i,1)+wneg(i,1); 
end 
Work=W/12;  % per cylinder 
MEP=W/(Vol(181,1)-Vol(1,1));  % total for 4 cylinders 
Power=(rpm*numcyl*Work)/(60*2*550);  % total for 4 cylinders 
Torque=(Power*550*60)/(rpm*2*pi);  % total for 4 cylinders 
fprintf('\n') 
fprintf('Work/Cycle/Cylinder           %6.1f ft-lbf \n',Work)  % per 
cylinder 
fprintf('IMEP                          %6.1f psi \n',MEP)  % total for 4 
cylinders 
fprintf('Indicated Power               %6.1f hp \n',Power)  % total for 4 
cylinders 
fprintf('Indicated Torque              %6.1f ft-lbf \n',Torque)  % total 
for 4 cylinders  
  
% Integrating P-V diagram (raw data) 
for c=1:1:cycles 
    j=1; 
    for i=1:1:180 
        fposR(j,c)=abs(Pm(i+361,c)-Pm(361-i,c)); 
        j=j+1; 
    end 
    j=1; 
    for i=1:1:179 
        wposR(j,c)=abs((Vol(i+361,1)-
Vol(i+362,1))*((fposR(i,c)+fposR(i+1,c))/2)); 
        j=j+1; 
    end 
    j=1; 
    for i=1:1:179 
        fnegR(j,c)=-1*(abs(Pm(181-i,c)-Pm(541+i,c))); 
        j=j+1; 
    end 
    j=1; 
    for i=1:1:178 
        wnegR(j,c)=-1*abs((Vol(i+179,1)-
Vol(i+178,1))*((fnegR(i,c)+fnegR(i+1,c))/2)); 
        j=j+1; 
    end 
    WR=0; 
    for i=1:1:178 
        WR=WR+wposR(i,c)+wnegR(i,c); 
    end 
    MEPR(1,c)=WR/(Vol(181,1)-Vol(1,1));  % total for 4 cylinders 
end 
  
IMEPstd=std(MEPR);  % standard deviation of IMEP 
COVIMEP=(IMEPstd/MEP)*100;  % coefficient of variation of IMEP (COV) 
IMEPmin=min(MEPR);  % minimum IMEP value 
LNVIMEP=(IMEPmin/MEP)*100;  % lowest normalized value of IMEP (LNV) 
  
 
   
82 
 
fprintf('\n') 
fprintf('STD of IMEP                   %6.1f psi \n',IMEPstd)  % per cycle 
fprintf('COV of IMEP                   %6.1f percent \n',COVIMEP)  % per 
cycle 
fprintf('LNV of IMEP                   %6.1f percent \n',LNVIMEP)  % per 
cycle 
  
% Creating log-log P-V data 
Plog=log10(Pavg); 
Vlog=log10(Vol); 
  
% Calculating polytropic index using least squares regression 
n=polyend1-polystrt1+1; 
sumx=0; 
sumy=0; 
sumxy=0; 
sumx2=0; 
for i=polystrt1+361:1:polyend1+361 
    sumx=sumx+Vlog(i,1); 
    sumy=sumy+Plog(i,1); 
    sumxy=sumxy+(Vlog(i,1)*Plog(i,1)); 
    sumx2=sumx2+(Vlog(i,1)^2); 
end 
gamma1=-1*((n*sumxy)-(sumx*sumy))/((n*sumx2)-(sumx^2)); 
  
sumx=0; 
sumy=0; 
sumxy=0; 
sumx2=0; 
for i=polystrt2+361:1:polyend2+361 
    sumx=sumx+Vlog(i,1); 
    sumy=sumy+Plog(i,1); 
    sumxy=sumxy+(Vlog(i,1)*Plog(i,1)); 
    sumx2=sumx2+(Vlog(i,1)^2); 
end 
gamma2=((n*sumxy)-(sumx*sumy))/((n*sumx2)-(sumx^2)); 
  
fprintf('\n') 
fprintf('Polytropic Index 1            %6.2f  \n',gamma1) 
fprintf('Polytropic Index 2            %6.2f  \n',gamma2) 
  
% Calculating SOC  
slope=-0.5; 
start=polyend1+361; 
while slope>-gamma1-tolerance 
    slope=(Plog(start+1,1)-Plog(start,1))/(Vlog(start+1,1)-Vlog(start,1)); 
    start=start+1; 
end 
  
% combstrt720=start-1; 
SOC=combstrt720-360; 
  
% Calculating EOC and combustion duration  
slope2=1; 
endburn=polystrt2+361; 
   
83 
 
while slope2>-gamma2-tolerance 
    slope2=(Plog(endburn-1,1)-Plog(endburn,1))/(Vlog(endburn-1,1)-
Vlog(endburn,1)); 
    endburn=endburn-1; 
end 
  
% combend720=endburn-1; 
EOC=combend720-360; 
combdur=combend720-combstrt720; 
fprintf('\n') 
fprintf('SOC                           %6.0f degrees \n',SOC) 
fprintf('EOC                           %6.0f degrees \n',EOC) 
fprintf('Combustion Duration           %6.0f degrees \n',combdur) 
  
% Allocating for space 
Xb=zeros(combdur,1); 
deg=zeros(combdur,1); 
Xb2=zeros(combdur,1); 
HR2=zeros(combdur,1); 
deg2=zeros(combdur,1); 
deg3=zeros(combdur,1); 
Wiebe=zeros(combdur,1); 
  
% Calculating mass fraction burned (Rassweiler & Withrow Method) 
j=1; 
Theta1=SOC; 
for i=combstrt720:1:combend720 
    Xb(j,1)=(((Pavg(i,1)^(1/gamma1))*Vol(i,1))-... 
        ((Pavg(combstrt720,1)^(1/gamma1))*Vol(combstrt720,1)))/... 
        (((Pavg(combend720,1)^(1/gamma1))*Vol(combend720,1))-... 
        ((Pavg(combstrt720,1)^(1/gamma1))*Vol(combstrt720,1))); 
    deg(j,1)=Theta1; 
    j=j+1; 
    Theta1=Theta1+1; 
end 
  
% Beginning 1st Law of Thermodynamics analysis on system 
Psi=Pavg.*(101325/14.696);  % convert pressure to SI units [Pa] 
Pmsi=Pavgmot.*(101325/14.696);  % convert pressure to SI units [Pa] 
Volsi=Vol.*((2.54^3)/(100^3));  % convert volume to SI units [m^3] 
  
% Calculating charge mass 
[FA,Mch,cvch,cPch,uch,hch,soch,dsodTch]=properties(fuel,phi,0,Tim); 
Rch=Ru/Mch;  % specific gas constant for charge [J/kgK] 
mch=etav*((Pref*Vdm)/(Rch*Tref));  % charge mass [kg] 
  
% Calculating residual mass 
[FA,Mr,cvr,cPr,ur,hr,sor,dsodTr]=properties(fuel,phi,1,EGT); 
Rr=Ru/Mr;  % specific gas constant for residual [J/kgK] 
mr=(Psi(EVC,1)*Vcm)/(Rr*EGT);  % residual mass [kg] 
xr=mr/(mr+mch);  % residual mass fraction 
xch=1-xr;  % charge mass fraction 
     
% Calculating temperature profile using ideal gas law 
mt=mch+mr;  % total mass in cylinder [kg] 
   
84 
 
for i=1:1:720 
    T(i,1)=((Psi(i,1)*Volsi(i,1))/((Rch*mch)+(Rr*mr)));  % [K] 
end 
  
% Woschni's correlation for engine heat transfer 
Ahead=pi*(borem/2)^2;  % combustion chamber area of head [m^2] 
Apiston=pi*(borem/2)^2;  % area of piston crown [m^2] 
for i=1:1:720  
    % Calculating average cylinder gas velocity [m/s] 
    if i>combstrt720 && i<EVO  % (combustion & expansion period) 
        w(i,1)=(c12*Sp)+(cal2*((Vdm*T(combstrt720,1))/... 
            (Psi(combstrt720,1)*Volsi(combstrt720,1)))*(Psi(i,1)-
Pmsi(i,1))); 
    elseif i>IVC && i<combstrt720  % (compression period) 
        w(i,1)=(c12*Sp); 
    else  % (gas exchange period) 
        w(i,1)=c11*Sp; 
    end 
    % Calculating Woshini's convection coefficient [KJ/kg] 
    h(i,1)=(3.26*cal1*(borem^Rem-1)*(Psi(i,1)^Rem)*(T(i,1)^(0.75-
(1.62*Rem)))*(abs(w(i,1))^Rem))/1000; 
    % Calculating area of cylinder wall [m^2] 
    Acyl(i,1)=2*pi*(borem/2)*((strokem/2)-
((strokem/2)*(cosd(Theta(i,1))))); 
    % Calculating the heat transfer rate across the system boundary [KW/s] 
    dQdT(i,1)=(h(i,1)*(Acyl(i,1)+Ahead+Apiston)*(T(i,1)-Twall)); 
end 
  
% Calculating the rate of work transfer done by the system [Nm/s] 
for i=1:1:719 
    PdVdt(i,1)=Psi(i,1)*((Volsi(i+1,1)-Volsi(i,1))/(Theta(i+1,1)-
Theta(i,1)))*(2*pi*N); 
end 
  
% Calculating mass fraction burned (1st Law of Thermodynamics Approach) 
Xbo=0; 
j=1; 
Theta2=SOC; 
for i=combstrt720:1:combend720 
    [FA,Mch,cvch,cPch,uch,hch,soch,dsodTch]=properties(fuel,phi,0,Tim); 
    [FA,Mbg,cvbg,cPbg,ubg,hbg,sobg,dsodTbg]=properties(fuel,phi,1,EGT); 
    [FA,Mr,cvr,cPr,ur,hr,sor,dsodTr]=properties(fuel,phi,1,EGT); 
    dXbdo=(((dQdT(i,1)-PdVdt(i,1))/(2*pi*N))-
(((mch*cvch)+(mr*cvr))*((T(i+1,1)-T(i,1))/... 
        (Theta(i+1,1)-Theta(i,1)))))/((mt)*(ubg-uch)); 
    if j==1 
    Xb2(j,1)=0; 
    else 
    Xb2(j,1)=Xb2(j-1,1)+dXbdo; 
    end 
    deg2(j,1)=Theta2; 
    j=j+1; 
    Theta2=Theta2+1; 
end 
  
 
   
85 
 
% 50% Burn Location Calculation 
n=rapdbrnend-rapdbrnstrt+1; 
sumx=0; 
sumy=0; 
sumxy=0; 
sumx2=0; 
for i=rapdbrnstrt-SOC:1:rapdbrnend-SOC 
    sumx=sumx+deg2(i,1); 
    sumy=sumy+Xb2(i,1); 
    sumxy=sumxy+(deg2(i,1)*Xb2(i,1)); 
    sumx2=sumx2+(deg2(i,1)^2); 
end 
rapdbrnslp=((n*sumxy)-(sumx*sumy))/((n*sumx2)-(sumx^2)); 
rapdbrnoffst=(1/n)*(sumy-(rapdbrnslp*sumx)); 
Centroid=((0.50-rapdbrnoffst)/rapdbrnslp); 
  
fprintf('Combustion Centroid           %6.1f degrees \n',Centroid) 
 
% ISFC calculation 
AF=lambda*14.7;  % air/fuel ratio 
mf=mch/AF;  % mass of fuel inducted per cycle 
TotE=mf*44E6;  % total energy of fuel (J) 
Wc=(Work*12*2.54*4.4482)/100;  %  work per cycle (J) 
nf=Wc/TotE;  % fuel conversion efficiency 
sfc=(2545/(nf*19000));  % specific fuel consumption (lbm/hp-hr) 
nfperc=nf*100;  % fuel conversion efficiency to percentage 
  
% Wiebe Function 
a=2.5;  % curve fit efficiency parameter 
m=0.7;  % curve fit form factor 
durfac=14; 
Theta3=SOC; 
j=1; 
for i=combstrt720:1:combend720+durfac 
    Wiebe(j,1)=1-exp(-a*(((Theta3-SOC)/(combdur+durfac))^(m+1))); 
    deg3(j,1)=Theta3; 
    Theta3=Theta3+1; 
    j=j+1; 
end 
Wiebedur=j-1; 
  
fprintf('\n') 
fprintf('Fuel Conversion Efficiency    %6.1f percent \n',nfperc)  % per 
cylinder 
fprintf('ISFC                          %6.3f lbm/hp-hr \n',sfc)  % per 
cylinder 
fprintf('Mass of Fuel per Cycle        %e kg \n',mf)  % 1 cylinder 
fprintf('\n') 
fprintf('Wiebe Efficiency Parameter    %6.2f  \n',a)  % a parameter 
fprintf('Wiebe Form Factor             %6.2f  \n',m)  % b parameter 
fprintf('Wiebe Duration                %6.0f degrees \n',Wiebedur)  % 
Wiebe combustion duration 
  
 
 
 
   
86 
 
% Plotting Data 
figure (1) 
plot(Theta,Pavg,'.b'); 
hold on 
plot(Theta,Pavgmot,'.r'); 
xlabel('Crankshaft Angle (Degrees)'); 
ylabel('Cylinder Pressure (psi)'); 
title('Pressure vs. Crankshaft Position'); 
legend('Fired Cycle','Motored Cycle'); 
hold off 
figure (2) 
plot(Theta,Pavg,'.k'); 
hold on 
plot(Theta,Pmedian,'.b'); 
plot(Theta,Pmode,'.r'); 
xlabel('Crankshaft Angle (Degrees)'); 
ylabel('Mean, Median, Mode of Cylinder Pressure (psi)'); 
title('Statistical Check for Cylinder Pressure Data'); 
legend('Mean','Median','Mode'); 
hold off 
figure (3) 
plot(Theta,Pstd,'.b'); 
xlabel('Crankshaft Angle (Degrees)'); 
ylabel('Cylinder Pressure Standard Deviation (psi)'); 
title('Standard Deviation of Pressure vs. Crankshaft Position'); 
figure (4) 
plot(Vol,Pavg,'.b'); 
xlabel('Cylinder Volume (cubic inches)'); 
ylabel('Cylinder Pressure (psi)'); 
title('P-V Diagram'); 
figure (5) 
plot(Vlog,Plog,'.b'); 
xlabel('Log of Cylinder Volume (cubic inches)'); 
ylabel('Log of Cylinder Pressure (psi)'); 
title('Log-Log P-V Diagram'); 
figure (6) 
% plot(deg,Xb,'.k'); 
% hold on 
plot(deg2,Xb2,'.r'); 
hold on 
plot(deg3,Wiebe,'b'); 
hold off 
ylim([0 1]) 
xlim([-15 35]) 
xlabel('Crankshaft Angle (Degrees)'); 
ylabel('Mass Fraction Burned, Xb'); 
title('Mass Fraction Burned vs. Crankshaft Position'); 
legend('Experimental Data','Wiebe Function',2); 
 
 
 
 
 
 
 
 
 
   
87 
 
MATLAB R2009a JANAF Table Thermodynamic Data Code 
 
%% properties.m [14] 
 
function [far,molwt,cv,cP,u,h,so,dsodT]=properties(fuel,phi,xb,T) 
global Ru 
Ru      = 8314.34;  % universal gas constant [J/kmolK] 
global atomCO2  aCO2 
global atomH2O  aH2O 
global atomN2   aN2 
global atomO2   aO2 
global atomCO   aCO 
global atomH2   aH2 
global atomH    aH 
global atomO    aO 
global atomOH   aOH 
global atomNO   aNO 
global atomfuel afuel 
atomCO2 =[ 0.10000e+01  0.00000e+01  0.20000e+01  0.00000e+01]; 
aCO2    =[ 0.24008e+01  0.87351e-02 -0.66071e-05  0.20022e-08  0.63274e-15 
-0.48378e+05  0.96951e+01; 
           0.44608e+01  0.30982e-02 -0.12393e-05  0.22741e-09 -0.15526e-13 
-0.48961e+05 -0.98636e+00]; 
atomH2O =[ 0.00000e+01  0.20000e+01  0.10000e+01  0.00000e+01]; 
aH2O    =[ 0.40701e+01 -0.11084e-02  0.41521e-05 -0.29637e-08  0.80702e-12 
-0.30280e+05 -0.32270e+00; 
           0.27168e+01  0.29451e-02 -0.80224e-06  0.10227e-09 -0.48472e-14 
-0.29906e+05  0.66306e+01]; 
atomN2  =[ 0.00000e+01  0.00000e+01  0.00000e+01  0.20000e+01]; 
aN2     =[ 0.36748e+01 -0.12082e-02  0.23240e-05 -0.63218e-09 -0.22577e-12 
-0.10612e+04  0.23580e+01; 
           0.28963e+01  0.15155e-02 -0.57235e-06  0.99807e-10 -0.65224e-14 
-0.90586e+03  0.61615e+01]; 
atomO2  =[ 0.00000e+01  0.00000e+01  0.20000e+01  0.00000e+01]; 
aO2     =[ 0.36256e+01 -0.18782e-02  0.70555e-05 -0.67635e-08  0.21556e-11 
-0.10475e+04  0.43053e+01; 
           0.36220e+01  0.73618e-03 -0.19652e-06  0.36202e-10 -0.28946e-14 
-0.12020e+04  0.36151e+01]; 
atomCO  =[ 0.10000e+01  0.00000e+01  0.10000e+01  0.00000e+01]; 
aCO     =[ 0.37101e+01 -0.16191e-02  0.36924e-05 -0.20320e-08  0.23953e-12 
-0.14356e+05  0.29556e+01; 
           0.29841e+01  0.14891e-02 -0.57900e-06  0.10364e-09 -0.69354e-14 
-0.14245e+05  0.63479e+01]; 
atomH2  =[ 0.00000e+01  0.20000e+01  0.00000e+01  0.00000e+01]; 
aH2     =[ 0.30575e+01  0.26765e-02 -0.58099e-05  0.55210e-08 -0.18123e-11 
-0.98890e+03 -0.22997e+01; 
           0.31002e+01  0.51119e-03  0.52644e-07 -0.34910e-10  0.36945e-14 
-0.87738e+03 -0.19629e+01]; 
atomH   =[ 0.00000e+01  0.10000e+01  0.00000e+01  0.00000e+01]; 
aH      =[ 0.25000e+01  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.25472e+05 -0.46012e+00; 
           0.25000e+01  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.25472e+05 -0.46012e+00]; 
atomO   =[ 0.00000e+01  0.00000e+01  0.10000e+01  0.00000e+01]; 
aO      =[ 0.29464e+01 -0.16382e-02  0.24210e-05 -0.16028e-08  0.38907e-12  
0.29148e+05  0.29634e+01; 
   
88 
 
           0.55420e+01 -0.27551e-04 -0.31028e-08  0.45511e-11 -0.43680e-15  
0.29231e+05  0.49203e+01]; 
atomOH  =[ 0.00000e+01  0.10000e+01  0.10000e+01  0.00000e+01]; 
aOH     =[ 0.38376e+01 -0.10779e-02  0.96830e-06  0.18714e-09 -0.22571e-12  
0.36413e+04  0.49370e+00; 
           0.29106e+01  0.95932e-03 -0.19442e-06  0.13757e-10  0.14224e-15  
0.39354e+04  0.54423e+01]; 
atomNO  =[ 0.00000e+01  0.00000e+01  0.10000e+01  0.10000e+01]; 
aNO     =[ 0.40459e+01 -0.34182e-02  0.79819e-05 -0.61139e-08  0.15919e-11  
0.97454e+04  0.29975e+01; 
           0.31890e+01  0.13382e-02 -0.52899e-06  0.95919e-10 -0.64848e-14  
0.98283e+04  0.67458e+01]; 
atomfuel=[ 0.10000e+01  0.40000e+01  0.00000e+01  0.00000e+01; 
           0.70000e+01  0.17000e+02  0.00000e+01  0.00000e+01; 
           0.14400e+02  0.24900e+02  0.00000e+01  0.00000e+01; 
           0.10000e+01  0.40000e+01  0.10000e+01  0.00000e+01; 
           0.10000e+01  0.30000e+01  0.20000e+01  0.10000e+01]; 
afuel   =[ 0.19713e+01  0.78716e-02 -0.10486e-05  0.00000e+00  0.00000e+00 
-0.99304e+04  0.88737e+01; 
           0.40652e+01  0.60977e-01 -0.18801e-04  0.00000e+00  0.00000e+00 
-0.35880e+05  0.15450e+02; 
           0.79710e+01  0.11954e+00 -0.36858e-04  0.00000e+00  0.00000e+00 
-0.19385e+05 -0.17879e+01; 
           0.17798e+01  0.12625e-01 -0.36249e-05  0.00000e+00  0.00000e+00 
-0.25254e+05  0.15088e+02; 
           0.14126e+01  0.20871e-01 -0.81421e-05  0.00000e+00  0.00000e+00 
-0.10264e+05  0.19171e+02]; 
if T<1000  
    arange=1; 
else  
    arange=2; 
end 
t        = T/1000; 
lnK      = 2.743-1.761/t-1.611/t^2+0.2803/t^3; 
K        = exp(lnK); 
as       = atomfuel(fuel,1)+atomfuel(fuel,2)/4-atomfuel(fuel,3)/2; 
if phi<1 
    nCO2 = atomfuel(fuel,1); 
    nH2O = atomfuel(fuel,2)/2; 
    nN2  = atomfuel(fuel,4)/2+3.773*as/phi; 
    nO2  = as*(1/phi-1); 
    nCO  = 0; 
    nH2  = 0; 
else 
    d1   = 2*as*(1-1/phi); 
    a1   = 1-K; 
    b1   = atomfuel(fuel,2)/2+K*atomfuel(fuel,1)-d1*(1-K); 
    c1   =-atomfuel(fuel,1)*d1*K; 
    nCO  = (-b1+sqrt(b1^2-4*a1*c1))/2/a1; 
    nCO2 = atomfuel(fuel,1)-nCO; 
    nH2O = atomfuel(fuel,2)/2-d1+nCO; 
    nN2  = atomfuel(fuel,4)+3.773*as/phi; 
    nO2  = 0; 
    nH2  = d1-nCO; 
end 
yfuelch  = 1/(1+4.773*as/phi); 
yO2ch    = as/phi/(1+4.773*as/phi); 
   
89 
 
yN2ch    = 3.773*as/phi/(1+4.773*as/phi); 
yCO2bg   = nCO2/(nCO2+nH2O+nN2+nO2+nCO+nH2); 
yH2Obg   = nH2O/(nCO2+nH2O+nN2+nO2+nCO+nH2); 
yN2bg    = nN2/(nCO2+nH2O+nN2+nO2+nCO+nH2); 
yO2bg    = nO2/(nCO2+nH2O+nN2+nO2+nCO+nH2); 
yCObg    = nCO/(nCO2+nH2O+nN2+nO2+nCO+nH2); 
yH2bg    = nH2/(nCO2+nH2O+nN2+nO2+nCO+nH2); 
molwtfuel= 
atomfuel(fuel,1)*12.01115+atomfuel(fuel,2)*1.00797+atomfuel(fuel,3)*15.999
4+atomfuel(fuel,4)*14.0067; 
molwtCO2 = atomCO2(1)*12.01115 +atomCO2(2)*1.00797 +atomCO2(3)*15.9994 
+atomCO2(4)*14.0067; 
molwtH2O = atomH2O(1)*12.01115 +atomH2O(2)*1.00797 +atomH2O(3)*15.9994 
+atomH2O(4)*14.0067; 
molwtN2  = atomN2(1)*12.01115  +atomN2(2)*1.00797  +atomN2(3)*15.9994  
+atomN2(4)*14.0067; 
molwtO2  = atomO2(1)*12.01115  +atomO2(2)*1.00797  +atomO2(3)*15.9994  
+atomO2(4)*14.0067; 
molwtCO  = atomCO(1)*12.01115  +atomCO(2)*1.00797  +atomCO(3)*15.9994  
+atomCO(4)*14.0067; 
molwtH2  = atomH2(1)*12.01115  +atomH2(2)*1.00797  +atomH2(3)*15.9994  
+atomH2(4)*14.0067; 
molwtch  = yfuelch*molwtfuel+yO2ch*molwtO2  +yN2ch*molwtN2; 
molwtbg  = yCO2bg*molwtCO2  +yH2Obg*molwtH2O+yN2bg*molwtN2... 
          +yO2bg*molwtO2    +yCObg*molwtCO  +yH2bg*molwtH2; 
xfuelch  = yfuelch*molwtfuel/molwtch; 
xO2ch    = yO2ch*molwtO2/molwtch; 
xN2ch    = yN2ch*molwtN2/molwtch; 
far      = xfuelch/(xO2ch+xN2ch); 
xCO2bg   = yCO2bg*molwtCO2/molwtbg; 
xH2Obg   = yH2Obg*molwtH2O/molwtbg; 
xN2bg    = yN2bg*molwtN2/molwtbg; 
xO2bg    = yO2bg*molwtO2/molwtbg; 
xCObg    = yCObg*molwtCO/molwtbg; 
xH2bg    = yH2bg*molwtH2/molwtbg; 
xfuel    = (1-xb)*xfuelch; 
xO2      = (1-xb)*xO2ch+xb*xO2bg; 
xN2      = (1-xb)*xN2ch+xb*xN2bg; 
xCO2     = xb*xCO2bg; 
xH2O     = xb*xH2Obg; 
xCO      = xb*xCObg; 
xH2      = xb*xH2bg; 
molwt    = 1/(((1-xb)/molwtch)+(xb/molwtbg)); 
cPfuel   = 
afuel(fuel,1)+afuel(fuel,2)*T+afuel(fuel,3)*T^2+afuel(fuel,4)*T^3+afuel(fu
el,5)*T^4; 
cPfuel   = cPfuel*Ru/molwtfuel; 
hfuel    = 
afuel(fuel,1)+afuel(fuel,2)/2*T+afuel(fuel,3)/3*T^2+afuel(fuel,4)/4*T^3+af
uel(fuel,5)/5*T^4+afuel(fuel,6)/T; 
hfuel    = hfuel*Ru/molwtfuel*T; 
sofuel   = 
afuel(fuel,1)*log(T)+afuel(fuel,2)*T+afuel(fuel,3)/2*T^2+afuel(fuel,4)/3*T
^3+afuel(fuel,5)/4*T^4+afuel(fuel,7); 
sofuel   = sofuel*Ru/molwtfuel; 
dsodTfuel= cPfuel/T; 
   
90 
 
cPCO2    = 
aCO2(arange,1)+aCO2(arange,2)*T+aCO2(arange,3)*T^2+aCO2(arange,4)*T^3+aCO2
(arange,5)*T^4; 
cPCO2    = cPCO2*Ru/molwtCO2; 
hCO2     = 
aCO2(arange,1)+aCO2(arange,2)/2*T+aCO2(arange,3)/3*T^2+aCO2(arange,4)/4*T^
3+aCO2(arange,5)/5*T^4+aCO2(arange,6)/T; 
hCO2     = hCO2*Ru/molwtCO2*T; 
soCO2    = 
aCO2(arange,1)*log(T)+aCO2(arange,2)*T+aCO2(arange,3)/2*T^2+aCO2(arange,4)
/3*T^3+aCO2(arange,5)/4*T^4+aCO2(arange,7); 
soCO2    = soCO2*Ru/molwtCO2; 
dsodTCO2 = cPCO2/T; 
cPH2O    = 
aH2O(arange,1)+aH2O(arange,2)*T+aH2O(arange,3)*T^2+aH2O(arange,4)*T^3+aH2O
(arange,5)*T^4; 
cPH2O    = cPH2O*Ru/molwtH2O; 
hH2O     = 
aH2O(arange,1)+aH2O(arange,2)/2*T+aH2O(arange,3)/3*T^2+aH2O(arange,4)/4*T^
3+aH2O(arange,5)/5*T^4+aH2O(arange,6)/T; 
hH2O     = hH2O*Ru/molwtH2O*T; 
soH2O    = 
aH2O(arange,1)*log(T)+aH2O(arange,2)*T+aH2O(arange,3)/2*T^2+aH2O(arange,4)
/3*T^3+aH2O(arange,5)/4*T^4+aH2O(arange,7); 
soH2O    = soH2O*Ru/molwtH2O; 
dsodTH2O = cPH2O/T; 
cPN2     = 
aN2(arange,1)+aN2(arange,2)*T+aN2(arange,3)*T^2+aN2(arange,4)*T^3+aN2(aran
ge,5)*T^4; 
cPN2     = cPN2*Ru/molwtN2; 
hN2      = 
aN2(arange,1)+aN2(arange,2)/2*T+aN2(arange,3)/3*T^2+aN2(arange,4)/4*T^3+aN
2(arange,5)/5*T^4+aN2(arange,6)/T; 
hN2      = hN2*Ru/molwtN2*T; 
soN2     = 
aN2(arange,1)*log(T)+aN2(arange,2)*T+aN2(arange,3)/2*T^2+aN2(arange,4)/3*T
^3+aN2(arange,5)/4*T^4+aN2(arange,7); 
soN2     = soN2*Ru/molwtN2; 
dsodTN2  = cPN2/T; 
cPO2     = 
aO2(arange,1)+aO2(arange,2)*T+aO2(arange,3)*T^2+aO2(arange,4)*T^3+aO2(aran
ge,5)*T^4; 
cPO2     = cPO2*Ru/molwtO2; 
hO2      = 
aO2(arange,1)+aO2(arange,2)/2*T+aO2(arange,3)/3*T^2+aO2(arange,4)/4*T^3+aO
2(arange,5)/5*T^4+aO2(arange,6)/T; 
hO2      = hO2*Ru/molwtO2*T; 
soO2     = 
aO2(arange,1)*log(T)+aO2(arange,2)*T+aO2(arange,3)/2*T^2+aO2(arange,4)/3*T
^3+aO2(arange,5)/4*T^4+aO2(arange,7); 
soO2     = soO2*Ru/molwtO2; 
dsodTO2  = cPO2/T; 
cPCO     = 
aCO(arange,1)+aCO(arange,2)*T+aCO(arange,3)*T^2+aCO(arange,4)*T^3+aCO(aran
ge,5)*T^4; 
cPCO     = cPCO*Ru/molwtCO; 
   
91 
 
hCO      = 
aCO(arange,1)+aCO(arange,2)/2*T+aCO(arange,3)/3*T^2+aCO(arange,4)/4*T^3+aC
O(arange,5)/5*T^4+aCO(arange,6)/T; 
hCO      = hCO*Ru/molwtCO*T; 
soCO     = 
aCO(arange,1)*log(T)+aCO(arange,2)*T+aCO(arange,3)/2*T^2+aCO(arange,4)/3*T
^3+aCO(arange,5)/4*T^4+aCO(arange,7); 
soCO     = soCO*Ru/molwtCO; 
dsodTCO  = cPCO/T; 
cPH2     = 
aH2(arange,1)+aH2(arange,2)*T+aH2(arange,3)*T^2+aH2(arange,4)*T^3+aH2(aran
ge,5)*T^4; 
cPH2     = cPH2*Ru/molwtH2; 
hH2      = 
aH2(arange,1)+aH2(arange,2)/2*T+aH2(arange,3)/3*T^2+aH2(arange,4)/4*T^3+aH
2(arange,5)/5*T^4+aH2(arange,6)/T; 
hH2      = hH2*Ru/molwtH2*T; 
soH2     = 
aH2(arange,1)*log(T)+aH2(arange,2)*T+aH2(arange,3)/2*T^2+aH2(arange,4)/3*T
^3+aH2(arange,5)/4*T^4+aH2(arange,7); 
soH2     = soH2*Ru/molwtH2; 
dsodTH2  = cPH2/T; 
cP       = 
xfuel*cPfuel+xCO2*cPCO2+xH2O*cPH2O+xN2*cPN2+xO2*cPO2+xCO*cPCO+xH2*cPH2; 
h        = 
xfuel*hfuel+xCO2*hCO2+xH2O*hH2O+xN2*hN2+xO2*hO2+xCO*hCO+xH2*hH2; 
so       = 
xfuel*sofuel+xCO2*soCO2+xH2O*soH2O+xN2*soN2+xO2*soO2+xCO*soCO+xH2*soH2; 
dsodT    = 
xfuel*dsodTfuel+xCO2*dsodTCO2+xH2O*dsodTH2O+xN2*dsodTN2+xO2*dsodTO2+xCO*ds
odTCO+xH2*dsodTH2; 
cv       = cP-Ru/molwt; 
u        = h-Ru/molwt*T;