Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Experimental and numerical simulation of charge motion in internal combustion engines Anetor, Lucky 1994

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1994-953047.pdf [ 5.64MB ]
JSON: 831-1.0088033.json
JSON-LD: 831-1.0088033-ld.json
RDF/XML (Pretty): 831-1.0088033-rdf.xml
RDF/JSON: 831-1.0088033-rdf.json
Turtle: 831-1.0088033-turtle.txt
N-Triples: 831-1.0088033-rdf-ntriples.txt
Original Record: 831-1.0088033-source.json
Full Text

Full Text

EXPERIMENTAL AND NUMERICAL SIMULATION OF CHARGEMOTION IN INTERNAL COMBUSTION ENGINESByLucky AnetorB. Sc. (First Class Hons.) (Mech. Eng.) University of Lagos, NIGERIA, 1982MEAC Royal Naval Engineering College, Manadon, UK, 1983M. Sc. (Marine Engineering) University of London, UK, 1986A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESMECHANICAL ENGINEERINGWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAAugust 1994© Lucky Anetor, 1994In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)Department of Mechanical EngineeringThe University of British ColumbiaVancouver, CanadaDate October 6, 1994DE-6 (2188)AbstractA combined study based on Laser Doppler velocimetry measurements and numericalsimulation was used to investigate the effects of piston-bowl geometry and intake swirllevel on the structure and evolution of the flow field near the top dead centre positionof an internal combustion engine. The University of British Columbia Rapid Intakeand Compression Machine (UBCRICM), was used for the experimental studies, whilenumerical calculations were performed with the KIVA-Il code. The KIVA-Il code is an in-cylinder fluid dynamics flow computer program, incorporating the following submodels;spray dynamics, species transport, mixing, chemical reactions and heat release rate.The effects of turbulence are represented by either the k-e or the subgrid scale (SGS)turbulence models. The k-e turbulence model option was used in the present study.The results of this study show that high-shear regions and the consequent turbulenceproduction occurred near the bowl entrance around top dead centre of compression.These regions were created before top dead centre in both chambers by the interactionof swirl and squish. The flow parameters (mean radial and tangential velocities andturbulence intensity) resulting from the re-entrant bowl configuration were found to beconsistently higher than those generated by the bowl-in-piston arrangement. The computational results were compared to the experimental data. A very good agreement wasobtained between predicted and measured mean radial and tangential velocity distributions. The numerical code was found t under estimate the measured turbulence intensityvalues in both chamber designs, however, a reasonable level of agreement was obtainedwith both chamber configurations.This study also highlighted the importance of accurately specifying the initial and11boundary conditions. It was found that the very impressive level of agreement betweenmeasured and computed results was largely due to the fact that the initial and boundaryconditions were accurately modelled to mimic actual experimental conditions.The investigation of the evolution of certain flow and turbulence quantities whichcould not be obtained experimentally was also undertaken with the KIVA-Il code. Acomprehensive discussion of these flow quantities was presented and the various waysthey are likely to affect engine performance in terms of combustion were highlighted.‘UTable of ContentsAbstractTable of Contents ivList of Tables xList of Figures xiList of Symbols xviAcknowledgement xxiv1 INTRODUCTION 11.1 Background 11.2 Objectives of Present Study 62 LITERATURE REVIEW 82.1 Introduction 82.2 Experimental Studies 102.3 Summary 232.4 Numerical Simulation Studies 273 EXPERIMENTAL AND DATA ANALYSIS METHOD 373.1 Experimental System 373.2 Modifications made to the UBCRICM 39iv3.2.1 Intake Arrangement .3.2.2 Optical Window3.3 Optical Setup3.4 Data Acquisition System3.5 Experimental Procedure3.5.1 Techniques Used to Ensure Data3.6 Data Analysis3.6.1 Issues3.7 Review of Data Analysis Methods3.7.1 Ensemble-Averaging Analysis3.7.2 Cycle-Resolved Technique3.7.3 Filter Implementation3.7.4 Cutoff-Frequency Selection3.7.5 Missing Data3.7.6 Conditional Sampling3.7.7 Autocorrelation Functions3.8 Data Analysis Method of Present Work . .4 NUMERICAL MODELLING4.1 General Background4.1.1 General Solution Procedure4.1.2 Turbulence Models4.1.3 Mesh Generation Scheme4.2 The Governing Equations4.2.1 The Fluid Phase4.3 The Computational DomainQuality39404041414344444546485051535455606161626364646570v4.4 Boundary Conditions4.4.1 Walls4.4.2 Velocities and Pressure4.4.3 Stagnation Enthalpy4.4.4 Turbulent Kinetic Energy and Dissipation4.5 The Numerical Scheme4.5.1 Temporal Differencing4.5.2 Spatial Differencing4.5.3 Lagrangian Phase Difference Equations4.5.4 Mass Density Equations4.5.5 Momentum Equation4.5.6 Internal Energy Equation4.5.7 Turbulence Kinetic Energy4.5.8 Dissipation Rate4.5.9 Equation of State4.6 Solution Procedure for Implicit Phase B Equations4.7 Accuracy Conditions and Automatic Timestep Control5 EXPERIMENTAL RESULTS AND5.1 Measurement Locations and Engine5.2 Results5.2.1 Mean Radial Velocity. .5.2.2 Mean Tangential Velocity5.2.3 Mean Velocity Close to Bowl Entrance Region5.2.4 Velocity Evolution on Cylinder Axis5.2.5 Radial Turbulence Intensity71717173Rate 74757576808182838484858586DISCUSSION 91Operating Conditions 91939495969798vi6 NUMERICAL RESULTS AND DISCUSSION6.1 Methods6.2 Modelling Approach6.2.1 Flat Top Piston6.2.2 Re-Entrant Bowl Piston6.2.3 Bowl-In-Piston6.2.4 Initial Conditions6.2.5 Numerical Accuracy6.2.6 Normalization of Results6.2.7 Parametric Studies .6.3 Results and Discussion6.3.1 Mean Radial Velocities6.3.2 Mean Tangential Velocities6.3.3 Intake Swirl Decay Rate6.3.4 Turbulence Intensity .6.3.5 Turbulence Length Scale6.3.6 Turbulence Diffusivity6.4 Discussion of the Computed Flow Fields6.4.1 Velocity Vector Fields6.4.2 Swirl Velocity Fields5.2.6 Tangential Turbulence Intensity5.2.7 Turbulence Intensity Close to Bowl Entrance5.2.8 Turbulence Intensity On The Cylinder Axis5.3 Measurement Uncertainty5.4 Discussion of ResultsRegion9899100100102106106107107108108109110111111112113113114115116116117117119vii6.4.3 Turbulent Kinetic Energy 1206.4.4 Turbulent Length Scale 1216.5 Discussion 1227 COMPARISON OF NUMERICAL AND EXPERIMENTAL RESULTS 1257.1 Modification Technique 1257.1.1 Mean Radial Velocities 1277.1.2 Mean Tangential Velocities 1277.1.3 Turbulence Intensities 1287.2 Discussion 1298 CONCLUSIONS AND RECOMMENDATIONS 1338.1 Conclusions 1338.1.1 Experimental 1338.1.2 Computational 1358.1.3 Data Comparison 1378.2 Recommendations 1388.2.1 Experimental 1388.2.2 Numerical 1398.2.3 Data Comparison 139References 140• Appendices 235A GOVERNING PARTIAL DIFFERENTIAL EQUATIONS 235A.1 General Framework 235A.2 Equations Governing Fluid Mechanics 236viiiA.3 Equations for Mean Quantities 237A.3.1 Turbulence Modelling 238B PROCEDURE FOLLOWED TO CONVERT KIVA-Il CODE 241B.1 CRAY Fortran Intrinsics 241B.2 Block Data Subroutine 242B.3 Real Variables are Double-Precision 242B.4 Comdeck Statement 242B.5 Debugging, Organizing and Testing of the Converted Code 243C DETERMINATION OF SWIRL RATIO 244C.1 Swirl Measurement 244ixList of Tables3.1 Rapid Intake and Compression Machine Operating Parameters 1504.1 Standard k - e Turbulence Model Constants 151xList of Figures2.1 Clerk’s Experiment To Demonstrate Role of Turbulence, Clerk (1921). . 1522.2 Ricardo’s Chamber For Producing Squish Flow (see Amann (1985)). . . 1522.3 Correlation of Flame Speed ratio with Turbulence Intensity, Lancaster(1976), Lancaster et al (1976) 1533.1 Schematic Diagram of UBCRICM, Dohring (1986) 1543.2 Schematic Diagram of Intake Arrangement 1553.3 Idealized Combustion Chamber Flow Configuration. Dimensions are inmillimetres 1563.4 Block Diagram of model 1400A Laser Velocimeter Interface 1573.5 Schematic diagram of Experimental Apparatus 1583.6 Schematic of Jet created by Flow through the Intake Valve indicating itsTurbulent Structure, Reynolds (1980) 1594.1 Typical Finite-difference Cell 1604.2 The portion of momentum cell (i, k) lying within regular cell (i, k).The three momentum cell faces lying within the regular cell are shaded. 1604.3 Cell Faces Associated with cell (i, k) 1614.4 The six points used to define the gradient of cell-centred quantity Q oncell face a 1614.5 Cell shapes before (dashed lines) and after (solid lines) the Lagrangia.nphase in a plane parallel flow 162xi5.1 Measurement Coordinate System Relative To The Piston at TDC. Dimensions in millimetres 1635.2 Bowl-in-Piston and Re-entrant bowl pistons. Dimensions are in millimetres. Crosses indicate LDV measurement locations 1645.3 Evolution of Measured Mean Radial Velocity with Crank Angle. CrossesIndicate Me.surement Location in millimetres 1655.4 Evolution of Measured Mean Radial Velocity with Crank Angle. CrossesIndicate Measurement Location in millimetres 1665.5 Variation of Measured Mean Tangential Velocity with Crank Angle. CrossesIndicate Measurement Location in millimetres 1675.6 Variation of Measured Mean Tangential Velocity with Crank Angle. CrossesIndicate Measurement Location in millimetres 1685.7 Mean Radial and Tangential Velocity Results Measured Close to BowlEntrance Region. Crosses Indicate Measurement Location in millimetres. 1695.8 Mean Radial and Tangential Velocity Results Measured on the CylinderAxis. Crosses Indicate Measurement Location in millimetres 1705.9 Evolution of Measured Radial Turbulence Intensity with Crank Angle.Crosses Indicate Measurement Location in millimetres 1715.10 Evolution of Measured Radial Turbulence Intensity with Crank Angle.Crosses Indicate Measurement Location in millimetres 1725.11 Evolution of Measured Tangential Turbulence Intensity with Crank Angle.Crosses Indicate Measurement Location in millimetres 1735.12 Evolution of Measured Tangential Turbulence Intensity with Crank Angle.Crosses Indicate Measurement Location in millimetres 174xii5.13 Radial and Tangential Turbulence Intensity Results Measured Close toBowl Entrance Region. Crosses Indicate Measurement Location in millimetres 1755.14 Radial and Tangential Turbulence Intensity Results Measured on the Cylinder Axis. Crosses Indicate Measurement Location in millimetres 1766.1 Finite Volume Mesh at 25 degree BTDC 1776.2 Finite Volume Mesh at 25 degree BTDC (3 - Dimensional) 1786.3 Computational Mesh near TDC 1796.4 Computational Mesh near TDC (3 - Dimensional) 1806.5 Effect of Bowl Shape on Radial Velocity 1816.6 Effect of Bowl Shape on Tangential Velocity 1826.7 Effect of Bowl Shape on Intake Swirl Decay Rate 1836.8 Effect of Bowl Shape on Turbulence Intensity 1846.9 Effect of Intake Swirl on Turbulence Intensity 1856.10 Effect of Intake Swirl on Turbulence Intensity 1866.11 Effect of Intake Swirl on Turbulence Intensity 1876.12 Effect of Bowl Shape on Turbulence Length Scale 1886.13 Effect of Bowl Shape on Turbulence Diffusivity 1896.14 Velocity Vector (cm/s) Plot across j = 1 Plane at 25 degree BTDC . 1906.15 Velocity Vector (cm/s) Plot across j = 1 Plane at 20 degree BTDC . 1916.16 Velocity Vector (cm/s) Plot across j = 1 Plane at 15 degree BTDC . 1926.17 Velocity Vector (crn/s) Plot across j = 1 Plane at 10 degree BTDC . 1936.18 Velocity Vector (cm/s) Plot across j = 1 Plane at 5 degree BTDC. . 1946.19 Velocity Vector (cm/s) Plot across j = 1 Plane at TDC 1956.20 Velocity Vector (cm/s) Plot across j = 1 Plane at 5 degree ATDC. . . . 1966.21 Velocity Vector (cm/s) Plot across j = 1 Plane at 10 degree ATDC. . . 197xlii6.22 Velocity Vector (cm/s) Plot across j = 1 Plane at 15 degree ATDC. 1986.23 Velocity Vector (cm/s) Plot across j = 1 Plane at 20 degree ATDC. 1996.24 Velocity Vector (cm/s) Plot across k = 6 Plane at 20 degree BTDC. 2006.25 Velocity Vector (cm/.s) Plot across k = 12 Plane at 20 degree BTDC. 2016.26 Velocity Vector (cm/s) Plot across k = 6 Plane at 15 degree BTDC. 2026.27 Velocity Vector (cm/s) Plot across k = 12 Plane at 15 degree BTDC. 2036.28 Velocity Vector (cm/s) Plot across k = 6 Plane at 5 degree BTDC. . 2046.29 Velocity Vector (cm/s) Plot across k = 12 Plane at 5 degree BTDC. 2056.30 Velocity Vector (cm/s) Plot across k = 6 Plane at TDC 2066.31 Velocity Vector (cm/s) Plot across k = 12 Plane at TDC 2076.32 Velocity Vector (cm/s) Plot across k = 6 Plane at 5 degree ATDC. . . 2086.33 Velocity Vector (cm/s) Plot across k 12 Plane at 5 degree ATDC. . 2096.34 Swirl Velocity (rad/s) Contour across j = 1 Plane at 20 degree BTDC. 2106.35 Swirl Velocity (racl/s) Contour across j = 1 Plane at 15 degree BTDC. 2116.36 Swirl Velocity (rad/s) Contour across j = 1 Plane at 10 degree BTDC. 2126.37 Swirl Velocity (rad/s) Contour across j = 1 Plane at 5 degree BTDC. 2136.38 Distribution of Turbulence Kinetic Energy across j = 1 Plane at 20 degreeBTDC 2146.39 Distribution of Turbulence Kinetic Energy across j = 1 Plane at 15 degreeBTDC 2156.40 Distribution of Turbulence Kinetic Energy across j = 1 Plane at 10 degreeBTDC 2166.41 Distribution of Turbulence Kinetic Energy across j = 1 Plane at 5 degreeBTDC 2176.42 Distribution of Turbulence Length Scale across j = 1 Plane at 20 degreeBTDC 218xiv7. of ComputedComparison of ComputedComparison of ComputedComparison of ComputedComparison of ComputedComparison of ComputedComparison of ComputedComparison of ComputedComparison of ComputedComparison of ComputedComparison of ComputedComparison of ComputedComparison of Computedand Measuredand Measuredand Measuredand Measuredand Measuredand Measuredand Measuredand Measuredand Measuredand Measuredand Measuredand Measured222223224225226227228229230231232233234CRAY and SUN SPARC Velocity Vector at 54 degree BTDCSchematic diagram of steady flow swirl apparatus6.43 Distribution of Turbulence Length Scale across j = 1 Plane at 15 degreeBTDC 2196.44 Distribution of Turbulence Length Scale across j = 1 Plane at 10 degreeBTDC 2206.45 Distribution of Turbulence Length Scale across j 1 Plane at 5 degree221and Measured Mean Radial VelocityMean Radial VelocityMean Radial VelocityMean Radial VelocityTangential VelocityTangential VelocityTangential VelocityTurbulence IntensityTurbulence IntensityTurbulence IntensityTurbulence IntensityTurbulence IntensityTurbulence IntensityB.1C.1247248xvList of SymbolsA Area of cell faceAa Area of subface, aACF Autocorrelation functionALE Arbitrary-Lagrangian-EulerianA0 Switch for selecting turbulent flow calculations(ANC),k Alternate node couplerATDC After top dead centreB Combustion chamber bore, Constant in law-of-the-wall formulationBTDC Before top dead centreCAD Crank angle degreeCd Courant numberC1, C2, c3 Turbulence modelling constantsCL Cc2, c, c, c3 Turbulence modelling constantsc Constant in law-of-the-wall formulationSpecific heat at constant pressure of species, mC3 Sound speed Courant numberxviSpecific heat at constant volumeD Diffusion coefficientDiameter of measurement volumeP Rate of momentum gain per unit volume due to sprayf0 Cut-off frequencyf3 Fringe spacingf Heating rate per unit area of wallSpecific body forceh Half distance between top of the cylinder head and piston head attop dead centreHWA Hot wire anemometerHz Hertz, cycles per secondi Particular engine cycle from a large number of cycles which makeup an ensembleI Specific internal energyI Unit dyadicSpecific internal energy of species, m, k Unit vectors in x, y and z directions respectivelyIMEP Indicated mean effective pressurexviiISFC Indicated specific fuel consumptionISNO Indicated specific nitrogen oxideI Heat flux vectorJ Wall Heat fluxk Turbulent kinetic energyThermal conductivity£ Turbulence length scaleIntegral length scaleLength of measurement volumeMicro length scaleLSGS Length scale in subgrid scale turbulence modelMa Mach numberMBF Mass burn fractionMBT Mean time for best torqueMEAC Marine Engineering Application CourseMHz Mega HertzMass of cell (i, j, k)Mass of momentum cell (i, j, k)xviiiMe” Mass of left cell-face control volume of cell (i, j k)n Cycle numberN Number of samplesN Number of cyclesNE Engine speedN Number of velocity measurements recorded in the cycleNm Number of measurementsTotal number of valid cross product pairsN Total number of measurementsp PressurePGS Pressure gradient scalingPr Prandtl number, ci/kPrE Prandtl number in turbulent kinetic energy dissipation equationPrk Prandtl number in turbulent kinetic energy equationQ Thermodynamic quantity, p, p, T, I, Pm, k andQA Computed value of Q at end of phase AQc Source term due to chemical heat releaseSource term due to sprayxixr Distance far less than integral length scalePosition vectorPosition vector of vertex (i, , k)Reference measurement coordinateR Constant in law-of-the-wall formulationRe Reynolds number, TTB/zirma Root mean squarerpm revolutions per minuteUniversal gas constantR(r) Spatial autocorrelation coefficientS Swirl numberSc Schmidt number, v/DSR Swirl ratio,w3/2IrNEt Time, Transpose of a matrixT Fluid temperatureTDC Top dead centreT Wall temperatureii Velocity vectorxxu’ Turbulence intensityu Shear speed, u/.,/c7U Measured radial or tangential component of velocityEnsemble averaged velocityU& In-cycle mean velocityUBCRICM University of British Columbia Rapid Intake aad Compression MachineUzf Low frequency velocity fluctuationEnsemble average of low frequency velocity fluctuationu1 Low frequency velocity fluctuation intensityuhf High frequency velocity fluctuationEnsemble average of high frequency velocity fluctuationsU’hf High frequency velocity fluctuation intensityv Cylinder volumeV Mean piston speed, 2LNEW Ensemble Averaged Mean Swirl Velocity at radius i’Wm Molecular weight of species, mWwajj Speed of wall in the z directionW8 Interaction with sprayxxix, y, z Cartesian coordinatesYm Mass fraction of species, mDistance from the wallGreek Symbolsa Dimensionless quantity used in pressure gradient scalingFaces of momentum control volume (i, , k)Dirac delta functionDissipation rate of turbulent kinetic energyReynolds number based on the gas velocity relative to the wallKolmogorov microscale8 Crank angle in degreesVon Kármán constantSecond coefficient of viscosity, eigen value in Bessel functions, Taylor MicroLength scale in Eqn. (3.24)p Coefficient of viscosity (Molecular)Effective viscosity (p + Pt)Pt Eddy viscosityKinematic viscosityxxiiV0 Diffusivityp Total mass densityPm Mass density of species, mWall shear stressT Lag or delayIntegral timescaleTk Kolmogorov microscaleTM Taylor micro length scaleqD Variable implicitness parameterVariable implicitness parameter for differencing the pressure gradient termEngine angular velocityLz Cell sizeLtC Convection timestepTime stepI6 Crank angle window in degreesWidth in delayV Del operatorxxiiiAcknowledgementFirst and foremost, I would like to express my special thanks to the Chief of NavalStaff, Nigerian Navy, for nominating me for this course. I wish to express my sinceregratitude to all my friends who have contributed in one way or another to the successfulcompletion of this research. I would also like to acknowledge with appreciation thesuggestions, counsel, moral support, equipment purchase and encouragement from myresearch supervisor, Professor R.L. Evans. The incisive comments of Professor P.G. Hill ofthe Department of Mechanical Engineering, University of British Columbia are gratefullyacknowledged. A special thanks to Mr. Tom Nicols of the University Computer Servicesfor his cheerful assistance during the KIVA-Il code conversion exercise. I would also liketo thank Dr. Zia Abdullah of the Department of Mechanical Engineering, University ofBritish Columbia, Dr. P.J. O’Rourke of Los Alamos National Laboratory, Los Alamos,NM. and Mr. T.L. Mckinley of Cummins Engine Co., Inc. Columbus, IN., for theirwillingness to always discuss various aspects of the KIVA-Il code with me. I also like tothank the staff of the Mechanical Engineering Workshop, especially Mr. Tony F. Besié,who fabricated all the experimental apparatus. To my wife Mrs. Ijeoma Kate Anetorand my children Omoaga, Ehidiamen, Yvonne and Olatokunbor and my nephew Marvis,I owe a debt of gratitude for all of the hours stolen from them. Their forbearance isgreatly appreciated. The typing of the manuscript was expertly done by my wife, Mrs.Ijeoma Kate Anetor, for this I say a big thank you.The financial support for this study was provided by the Canadian CommonwealthScholarship and Fellowship Plan, and the Ministry of Defence, Lagos, Nigeria.xxivChapter 1INTRODUCTION1.1 BackgroundOver the years there has been a constant decline in fossil fuel reserves and various Environmental Protection Agencies (EPAs) have been promulgating more stringent emissioncontrol standards for Internal Combustion Engines (ICEs). Another potentially harmfulfactor which is capable of aggravating the already deteriorating situation even further,is the political differences which often manifest themselves from time to time betweenmajor oil consumer nations of the West and leading petroleum producing countries of theMiddle East. The 1970 Arab-Israeli War, the just concluded Gulf War and the recentenvironmental concerns reflected in various national statutes, like the Clean Air Act lendfurther credence to the issues stated above. The question widely asked these days, is nolonger whether standards will be raised, but by how much. The available evidence (data)strongly suggest that the requirements concerning automobile emissions will definitelybe up-graded.It is estimated that in the United States of America alone the automobile accountsfor 50% of the petroleum consumption, 20% of the total energy consumption and 50% ofthe pollutants entering the atmosphere, Butler et al (1981). Considering these statisticsand the many more automobiles in use world-wide, it does not take much to realize thata relatively minor improvement in efficiency and emissions of ICEs will yield significantbenefits. Moreover, the concerns raised above combine to determine the market share,1Chapter 1. INTRODUCTION 2or simply put, the continued survival or otherwise of ICE manufacturers.The response of both ICE researchers and manufacturers in tackling these problemshas been unanimous. The mutual agreement from both parties has resulted in an unprecedented impetus for research into ways of alleviating these difficulties.It has long been realized that fluid motion (turbulence) within the combustion chamber of an internal combustion engine is one of the most important factors controlling thecombustion process. The degree and details of this influence on combustion are still notvery well understood. However, it can be safely asserted that the major roles of fluidmechanics in spark and compression ignition engines are first and foremost to preparethe mixture of fuel and air for combustion and secondly, to control flammability throughlarge and small scale mixing. The mixture preparation in conventional spark ignitionengine is achieved mainly by the induction system, which generates the bulk motionand turbulence required to mix both air and fuel prior to ignition. In diesel engines,the induction system (for example, swirl, tumble) and chamber geometry (for example,squish) control the air motion; whereas the details and mechanics of the fuel injectionsystem governs the spray dynamics (for example, penetration, atomization, break-up andevaporation), air utilization and pollutant formation. The flow pattern also affects therate of heat transfer to the cylinder walls, turbulent transport of mass, momentum andenergy.It is well known that turbulence intensity in an engine has a profound effect oncombustion rate. Turbulent flame speeds in a Spark Ignition (SI) engine are many timesthe laminar flame speed of the same fuel-air mixture at normal engine temperatures andpressures, Groff and Matekunas (1980). A better understanding of the characteristicsof engine turbulence and interactions of turbulence and bulk flow with the combustionprocess may help designers improve future engines. The nature of in-cylinder flows isvery complex — turbulent, unsteady, reactive and three dimensional. Due to the complexChapter 1. INTRODUCTION 3nature of combustion processes in an engine, insufficient knowledge has resulted in designswhich are not necessarily optimal. The design engineer has many options open to explorein order to maximize engine performance. Some of the alternatives are variation in intakegeometry, spark configuration, bore-stroke ratio, chamber geometry and ignition method.All of these can affect the flow-field and significantly influence engine performance.A Recent trend among engine designers is the development of lean-burn engines; thedirect-injection stratified-charge engine is one example. Lean-burning engines have several potential advantages; namely, they are thermodynamically more efficient because ofa higher ratio of specific heats due to composition and lower combustion temperatures,Matekunas (1983), they guarantee sufficient oxygen for complete combustion and theyresult in reduced pumping losses due to higher volumetric efficiency. A reduction in NOemissions is also another advantage, possibly due to the reduced operating temperature.Inspite of these attractive benefits, most engines operate with near stoichiometric mixtures because lean burn mixtures result in poor driveability due to slower combustionrate, increased cycle-to-cycle variation, increased hydrocarbon emissions, due to loweroperating temperatures, which are not sufficient to cause complete oxidation, and in theextreme case misfire can occur when the mixture is too lean.One possible method of alleviating these shortcomings, particularly slow burningrate, is by increasing both the swirl and turbulence levels. A fast burn engine hasseveral attractive attributes, namely; a more favorable tradeoff between NO emissionsand fuel economy, a tolerance to greater dilution and a reduction of cyclic variability.Furthermore, a fast burn engine allows delayed ignition which places the early flame kernelin an environment more conducive to strong propagation, that is, higher temperaturesand pressures. Although a fast burning engine is, in general, desirable, a too rapidpressure rise (greater than about 240 kPa per crank angle at 200-2000 revolutions perminute) has been found to result in undesirable roughness, Andon and Marks (1963). TheChapter 1. INTRODUCTION 4manifestation of this phenomenon could be noticed when automotive engine structuresbecome generally noisy.While fast burning of the air-fuel mixture is a fundamental approach for reducingfuel octane requirements, there are various design techniques available for achieving thisfast burn process. Higher turbulence in the engine cylinder at the time of ignition can inprinciple be achieved through;• Higher shear in the annular jet(s) entering the cylinder through the intake valves,• Squish-generated shear in certain combustion chamber geometries, and,• Compression of large scale structures (for example, Swirl and Tumble) persistingthrough induction and their subsequent decomposition just before ignition.Extensive research of the flow structure prevailing in the combustion chambers of internal combustion engines has shown that for practical valve lifts, the induction-generatedturbulence decays by the time of ignition to levels determined by the engine speed, (thatis, between 0.4- 0.5 times the mean piston speed) Arcoumanis and Whitelaw (1987).On the other hand, squish-generated turbulence peaks at about 340- 350 degrees crankangle after top dead centre of induction and is localized around the periphery of thebowl lip, which therefore reduces its effectiveness as a promoter of faster burning rates ininternal combustion engines. The third method, which offers real promise for turbulenceenhancement in internal combustion engines is the formation of swirl and tumble motionsduring induction. The swirl motion spins-up during compression in order to conserve itsangular momemtum. This results in higher shear within the combustion chamber, withthe consequent production of turbulence. On the other hand, the tumble motion breaksdown during the second half of the compression stroke and releases its stored kineticenergy. Since the time remaining up to ignition is inadequate to cause significant decay,Chapter 1. INTRODUCTION 5the turbulence levels that the flame encounters during propagation are higher than thoseprevailing in conventional engines.Several attempts, Arcoumanis (1990), Haworth et al (1990), Ahmadi-Befrui (1989),Kent et cii (1989), Uzman and Hazelton (1986), Gosman et cii (1985) have been made todemonstrate the advantages of both swirl- and barrel-generating induction systems, asmeans of promoting turbulence and increasing burning rates. However, there is incomplete knowledge of the detailed mechanisms of formation and destruction of large scalestructures such as tumble and its interaction with swirl. These two flow structures arecomplimentary to each other and their combined effect is likely to offer the best promisefor stable lean engine operation.During the last 20 years, a lot of effort has been devoted to developing multidimensional fluid dynamics codes capable of predicting flows in engines, with the hopeof reducing engine development time and cost. In order for these codes to become acceptable as a legitimate design tool, two conditions must be fulfilled. First and foremost,the accuracy of the numerical code must be demonstrated to meet some required standards. For example, in order to predict indicated fuel consumption, the model must becapable of predicting. the cylinder pressure history. This in turn requires accurate predictions of the fuel burning rate, and at least reasonable estimates of the heat transferfrom the cylinder gases. The burning rates and heat transfer to the cylinder wall are intimately linked to the cylinder flow fields. Since the cylinder pressures and temperaturesare coupled through the equations of state and the first law of thermodynamics, theseare also prerequisites for predicting emissions from engines. This perspective is largelybased on the accepted observation that pollutants, such as N0, are sensitive to the localtemperature-time history prevailing in the engine cylinder.The second condition which numerical codes must satisfy is the ability to make predictions with reasonable turn-around times on affordable computers. This criterion canChapter 1. INTRODUCTION 6be met through a combination of fast, inexpensive computers and efficient computationalalgorithms. The recent development of a new class of relatively cheap, high performanceSUN — SPARC work stations has helped in meeting computer speed and cost considerations. There have also been a number of improvements in numerical algorithms, Amsdenet al (1985b), Amsden el al (1989). Another technique for reducing the turn-around timeand computational cost associated with using these combustion codes is to solve simplified problems. For example, solving two-dimensional rather than three dimensionalproblems can easily reduce the run time by an order of magnitude. Furthermore, areduction in both computer time and cost can be made by solving non-reacting flows.From the foregoing, we see that detailed mechanism(s) of swirl/squish interaction andthe consequent turbulence generated is still lacking. Moreover, works featuring comprehensive code validation exercises where real boundary conditions have been utilized ininitializing the code have been very few. It is considerations such as these that dictatethe impetus for the present study.1.2 Objectives of Present StudyThe global objective of the present study is to provide further understanding into two ofthe areas generally recognized as important for further research in the internal combustionengine field. These areas are:• Swirl-squish interaction and the consequent turbulence generated,• To develop an experiment which has similar initial and boundary conditions (thatis, fixed swirl rate) to those encoded in KIVA-Il code. The results obtained fromthis specially designed experiment will provide a consistent set of data for validatingthe KIVA-Il code.Chapter 1. INTRODUCTION 7Another goal of the present study is to conduct some parametric studies with theKIVA-JI code. The KIVA-Il code is an in-cylinder fluid dynamics flow code incorporatingthe following submodels: spray dynamics, species transport, mixing chemical reactionsand heat release rate.Although this study is based on a spark ignition engine, some of the results may proveuseful in the diesel engine environment, where the characteristics of the flow field areimportant factors influencing such parameters as, spray penetration, break-up, mixing,evaporation rate and so on.Chapter 2LITERATURE REVIEW2.1 IntroductionThis literature survey is undertaken in order to illustrate the progress made so far, asregards meeting the dual objectives of fuel economy and emission reduction. The mainfocus is on. the various efforts aimed at understanding the interaction between swirl andsquish as a means of enhancing near top dead centre turbulence and consequently engineperformance.It has long been realised that the characteristics of turbulence around top dead centreposition is very important for internal combustion engine performance, since this is whatcontrols the air-fuel mixing and burning rates. In cylindrical chambers with flat pistons inthe absence of squish, the overall turbulence decays continuously during the compressionprocess despite the turbulence production by compressive strain, probably due to theeffects of various dissipative flow mechanisms, Gosman et al (1985). At top dead centreof compression and before ignition, turbulence levels of 0.3 to 0.5 times the mean pistonspeed have been measured under different conditions, Arcoumanis et al (1985, 1982),Hayder et al (1985), with tendencies towards homogeneity and isotropy. These resultsindicate that in the absence of strong squish and long-lived rotational motions, the turbulence intensity at top dead centre of compression has a maximum limit of approximately0.5 times the mean piston speed and the influence of the inlet boundary conditions onthe flow field around top dead centre is relatively weak, Arcoumanis et al (1982). From8Chapter 2. LITERATURE REVIEW 9the foregoing, it can be concluded that in non-squish configurations, the maximum swirland turbulence intensity that can be obtained before ignition are controlled by angularmomentum decay and turbulence dissipation which impose limits to the desired top deadcentre flow field. In order to increase induction generated swirl and turbulence intensity,bowls of various shapes have been incorporated into the piston or the cylinder head toenhance fluid motion at the point in the engine cycle where it is mostly needed— thepoint of ignition. The diameter of the piston bowl in cylindrical configurations and thebowl-entry diameter in re-entrant-type geometries play an important role in determiningthe in-bowl flow field (see Arcoumanis and Whitelaw (1987)). The bowl-entry diameterirrespective of the bowl eccentricity from the cylinder axis, determines the amplificationfactor of the induction swirl, since conservation of angular momentum dictates an increase in angular velocity proportional to the square of the ratio of the cylinder to thebowl diameters. This implies that the difficulties associated with generation of adequateswirl during induction can be resolved by a judicious choice of piston bowlIn practice, the swirl-generating characteristics of inlet ports are usually determinedeither by paddle wheel or swirl-torque meter methods and these are usually carried outunder steady flow conditions. These type of measurements can be used to derive a numbersuch as, swirl number or swirl ratio with which to rank the swirl-generating capabilityof the port. Some experience based type of recommendations can also be made fromsuch measurements, Monaghan and Pettifer (1981) or they could be used as a boundarycondition for estimating in-cylinder swirl throughout the engine cycle, Davis and Kent(1979). Although these methods can be applied quickly and provide useful criteria forcomparing different port designs, they yield no detailed structural information about theswirl. A theoretical study conducted by El-Tahry (1982) suggests that the swirl structureis of similar importance to the level. Experimental observations have also shown that theswirl structure at top dead centre contains features present at inlet valve closing, VafidisChapter 2. LITERATURE REVIEW 10(1984).Some experimental studies in the combustion chambers of real engines, Monaghanand Pettifer (1981), Matsuoka et at (1980), model engines, Morse et at (1980) and insteady flow rigs, Kajiyama et at (1984) have provided valuable spatial- and temporal..resolved velocity data. In addition they also pointed out the complexity and difficulty ofstudying such flows.One possible way of obtaining a more comprehensive flow information is to extendthe scope of these studies to full three-dimensional mappings. It is however, necessaryto point out that the enormity of such a task and the optical inaccessibility of largeportions of the cylinder make this an arduous undertaking. An alternative route that iscurrently being pursued is to calculate the details of the flow field by numerically solvingthe differential equations which govern its unsteady, three-dimensional and turbulentnature. Although this is also not an easy task, such an approach does, however, offer theprospect of obtaining reasonably complete flow-structural information that is sufficientfor most engineering purposes provided the results can be believed that is, if validatedto a reasonable degree.The reviews to follow will be treated under two broad categories, namely, experimentaland numerical studies. This order was chosen in view of the historical developments ofboth approaches as techniques for studying ICE combustion performance.2.2 Experimental StudiesThere is a vast number of reported studies pertaining to experimental investigations ofvarious aspects of the internal combustion engine. The following is a summary of someof the experimental studies relevant to the present study.Chapter 2. LITERATURE REVIEW 11Clerk (1921) conducted an engine experiment to demonstrate the importance of turbulence in hastening combustion flame speed. He employed a 229 x 432 mm (bore xstroke) engine running at 150 revolutions per minute and took a normal pressure-volume(p — v) diagram shown in Figure 2.1, with ignition occurring at a and combustion completed at b. Another set of experiments was conducted with a freshly inhaled charge,but with the valves deactivated. The trapped charge was then compressed and expandedthree times before igniting it at a’. Because of the decay of turbulence during those threecycles, combustion lasted until b’. The effect of lower turbulence level has resulted inlonger combustion time, (approximately doubled).Clerk’s conclusion regarding the data from his experiment was correct and to thepoint. He attributed the decrease in combustion time with increasing engine speed tothe turbulence or eddying caused by the rush of gases into the cylinder during the suctionstroke, which persisted during the compression stroke. It was also observed that theintake turbulence does not remain constant, but that it undergoes some degradationduring compression. As a result of this finding, Ricardo (see Amann (1985)) successfullycompensated for this decay in turbulence level with squish flow from the space betweenthe piston crown and cylinder head in his chamber design shown in Figure 2.2.The first attempt to look at the fluid motion inside the engine cylinder appears to havebeen by, Lee (1939), who photographed the movement of white feathers down throughthe glass cylinder of a motored research engine. His interest was the effect of shroudedinlet valves and directed inlet ports on the air movement within the cylinder. Differentconfigurations of these inlet assemblies were examined. The main observations from thisinvestigation can be summarized as follows:• For a plain intake valve the air motions created in the cylinder during intake persisted throughout the compression stroke with slowly decaying velocities. TheseChapter 2. LITERATURE REVIEW 12motions were, however, observed to die out soon after the start of the expansionstroke.• The tangential air flow created by the shrouded valve caused the air to rotaterapidly about the cylinder axis, this rotation persisted right up to the end of theexhaust stroke.• The induced velocities were approximately proportional to the engine speed andinversely proportional to the intake valve flow area.Semenov (1958) conducted an experimental investigation using a. hot-wireanemometer to measure instantaneous velocities in a motored engine. The mean velocitiesas a function of crank angle were obtained by averaging the signals from the anemometerover 24-degree crank angle intervals and then over 20 to 40 cycles. The fluctuation• velocity component was separated from the mean velocity signal with the aid of bandpass filters spanning five frequency ranges from 850 to 6200 Hz, thereby giving thespectral composition of the fluctuations. The intensities or root mean square (rms)values of the fluctuations were obtained over the same number of cycles as for the meanvelocities. The conclusions from this study were:• The gas velocity during intake varied considerably across the combustion chambers,reflecting the jet-like nature of the incoming flow. The velocities varied linearly withengine speed and decreased with distance from the cylinder axis. The large velocitygradients so formed were thought to be the major source of turbulence generation.• It was observed that the mean velocities and the fluctuating components decreasedrapidly during the first third of the compression stroke, after which both remainedfairly constant for the rest of the stroke.Chapter 2. LITERATURE REVIEW 13• When the compression ratio was varied between 4 to 9.5, the mean velocities attop dead centre were decreased by 20% and were found to vary with the square ofthe engine speed.• The turbulence spectrum indicated that most of the energy during the compressionstroke lies in the frequency range below 500 Hz.• The macro-scale of turbulent eddies formed during intake were estimated to be ofthe order of 1.8 - 2.2 mm.Lancaster (1976), employed tn-axial hot-wire anemometer to conduct a systematicstudy of the effects of turbulence on flame speed in an engine. Two types of intakevalves; namely, non-shrouded and swirl-producing shrouded intake valves were used inturn. Concurrently, a computer model was used to deduce flame speed from measurements of instantaneous cylinder pressure acquired with a piezoelectnic pressure transducer, Lancaster et at (1976). In order to separate the effects of turbulence from those ofthermodynamics and chemistry, Lancaster worked with flame speed ratio. Flame speedratio is defined as the ratio of turbulent burning velocity (relative to the unburned chargeahead of the flame) to the laminar flame speed. Laminar flame speed is a defined property of the mixture established by experimental measurement that depends on chemicalcomposition and thermodynamic state (see Heywood (1988)).From Lancaster’s calculated flame speeds, flame travel was divided into four zones,namely:• a flame-initiation period (during which the initial flame kernel is formed at thespark plug and grows into a recognizable flame front);• a flame-development period;Chapter 2. LITERATURE REVIEW 14• a fully developed period (during which flame speed ratio is reasonably constant);and,• a termination period (during which wall interference effects become significant).It was found that the largest fraction of the charge mass was consumed during the fullydeveloped period; this was perceived to be of major significance.When flame speed ratio was plotted against turbulence intensity for the fully developed period, a singe linear correlation resulted that fitted both non-shrouded andshrouded data acquired at different speeds, degrees of throttling, mixture ratios, sparkadvances and compression ratios, as shown in Figure 2.3. This correlation has been incorporated into a computer model that has proven useful for a priori evaluation of differentcombustion-chamber designs.Witze (1982) conducted an experimental study to investigate the effect of spark location on combustion in a variable swirl engine. A single cylinder engine of a Teledyne-Wisconsin type AENL with an air-cooled L-head block was used for the study.The important findings and conclusions from this study could be summarized asfollows:• It was observed that spark plug location has a significant effect on burning rates.In general, ignition at the centre of the chamber was found to yield optimal burnrate, except for condition of very high swirl (swirl number, S = 8.3), for whichside-wall ignition produced the fastest burn. The study also showed that, with noswirl present, the burn duration was simply a function of flame travel distance,exhibiting a nearly linear reduction in burn rate as the spark plug was moved awayfrom the cylinder axis.• It was noticed that with high swirl and side-wall ignition, the flame was driven byconvection around the perimeter of the combustion chamber. Side-wall ignition wasChapter2. LITERATURE REVIEW 15found to produce a faster burn than central ignition when the convection processis faster than the mechanism for a turbulent diffusion burn.• It was found that as the swirl rate was increased, the turbulence level in the cylinderwas reduced. The corresponding reduction in the diffusive flame propagation speedwas attributed to this observation. It was believed that the net effect of this wasthe two very distinct minima in the burn duration observed when the swirl wasvaried. The fastest burning rates occurred for the lowest swirl level used (S = 3.2)and the highest (S = 8.3). The slowest burning rates occurred for the conditions ofno swirl and S = 6.6. The explanation for the latter result was that the turbulencelevel was too low for a fast diffusive burn, and the swirl was not high enough for afast convection-driven burn.• A comparison between protruding and surface-gap spark plugs showed that coolwall-boundary-layers act to hold the flame at the wall in a partial slip condition.However, it was noticed that the wake produced by a protruding electrode was moreefficient at holding the flame. It was found that boundary-layer stabilization wasa viable mechanism for increasing flame speed in a swirling flow. The importanceof holding the flame in a swirling flow, regardless of the details of the mechanism,was shown to be particularly significant as the equivalence ratio is reduced.• The final conclusion from this investigation was that cyclic variation tends to increase with increased burn duration and decreased swirl level.This study shows that there is an optimal level of intake swirl, below or beyond whichengine performance deteriorates rapidly.Chapter 2. LITERATURE REVIEW 16Mikulec et al (1988) carried out an experimental investigation, with the aim of studying the separate effects of swirl and inducted kinetic energy on combustion in a homogeneous charge engine. Combustion measurements were made for two spark locations;namely, central and peripheral. The swirl ratio used in the study ranges from zero to2.8; this was achieved by gradual rotation of the intake port about the acis of the intakevalve.The main results from this study can be summarized as follows:• It was observed that changing the swirl ratio from zero to 2.8, while holding theinduced kinetic energy constant reduce the ignition delay period, (0-10%) massburn fraction (MBF) by about 25%, decrease the combustion duration, (10-90%)MBF, by approximately 10%, and achieved a significant improvement in combustion stability.• Combustion duration and ignition delay for peripheral ignition were found to belonger than for central ignition over the range of swirl ratios tested.Some of the reasons adduced for the observations mentioned above were:1. The turbulence intensity increased by about 10%, thereby lowering the combustionduration for either spark location;2. The combination of increased flame area ( 7%) for peripheral ignition and increased turbulence level were found to produce a slightly larger sensitivity of combustion duration to swirl than when ignition was centrally initiated.3. The decrease in length scale ( 10%) coupled with increased turbulence intensity resulted in a reduction of ignition delay period and a consequent increase incombustion stability.Chapter 2. LITERATURE REVIEW 17• It was noticed that for either spark location and an air/fuel ratio as large as 27:1,the combustion stability was greatest for the maximum swirl ratio (S = 2.8) tested.The combustion stability was however found to deteriorate rapidly as the mixturebecame leaner than 27:1 for the peripherally ignited case.• The combustion stability was found to correlate with the ignition delay period,(0-10%) MBF, with spark location and swirl ratio playing a predominant role.• It was observed that increased swirl and central ignition were capable of improvingboth the indicated specific fuel consumption (ISFC) and the indicated specificnitrogen oxide (ISNO) emissions.Liou et al(1983) conducted an experimental study to investigate the effect of swirl andno swirl on various flow quantities. The measurements were made using a back scatterlaser Doppler velocimetry arrangement. Measurements were taken in the mid-plane ofthe top dead centre clearance height from 34° BTDC to 30° ATDC in the no swirl casesand from 44° BTDC to 20° ATDC in the measurements involving swirl. The engine wasoperated at three speeds namely, 1200, 1800 and 2400 revolutions per minute. The dataacquired were reduced by the cycle-resolved analysis technique.The various conclusions arrived at are summarized below:• It was found that the ensemble averaged velocity scaled with the engine speed whenswirl was present, but there was no scaling in the absence of swirl.• The swirl number was found to decay during the compression stroke and to decreasewith increasing engine speed. The swirl number is an expression of the torqueexerted upon the flow-straightening honeycome of a swirl meter, whose axis isparallel to that of the engine cylinder, normalised by engine speed.Chapter 2. LITERATURE REVIEW 18• The turbulence intensity near top dead centre was found to be relatively homogeneous and to scale in an approximately linear fashion with engine speed with andwithout swirl. The turbulence intensity was also found to decay during compressionboth with and without swirl. However, the magnitude of the turbulence intensityfor a given engine speed was found to be from 25% to 50% greater when swirl waspresent.• It was also found that the turbulent energy shifted to higher frequencies with increasing engine speed irrespective of whether swirl was present or not, however, fora given speed more energy was found at higher frequencies with swirl than without.• The integral time scale was found to decrease with increasing engine speed bothwith and without swirl. When swirl was present the integral time scale attained aminimum value at top dead centre, whereas in the absence of swirl it was found toincrease up to and beyond the top dead centre position.• The cyclic fluctuations of the mean velocity at top dead centre were found to benearly twice as large without swirl than when swirl was present. It was also foundthat with swirl the cyclic fluctuation reached a minimum at top dead centre whereaswithout swirl they were found to decay up to and beyond top dead centre.From the various findings stated above, it was concluded that top dead centre turbulenceis relatively insensitive to the details the intake process.Hadded and Denbratt (1991) carried out an experimental study with the aim ofdetermining the influence of tumble motion on in-cylinder turbulence characteristics.• The investigation was carried out in a single cylinder research engine which was motoredat 1500 rpm. Four cylinder heads with varying tumble magnitudes were evaluated usingconventional and scanning laser Doppler anemometry measurements. The data acquiredChapter 2. LITERATURE REVIEW 19were analysed with a specially designed algorithm, which has the ability to account forthe effects of mean flow cyclic variations and system noise.The conclusions from this study are as follows:• It was found that steady flow derived tumble ratio serves as a good guide forincreasing the in-cylinder turbulence. The tumble ratio is an expression of thetorque exerted upon the flow-straightening honeycomb of a swirl meter, whose axisis normal to that of the engine cylinder, normalised by engine speed.• The experimental results show that increasing the tumble motion at inlet valveclosure by shrouding the inlet valves does not result in turbulence improvementsnear top dead centre. This was attributed to the early decay of some of the tumblevortex which was assumed to be due to strong valve-to-valve flow interaction resulting from reduced curtain area, which ultimately led to the formation of secondaryvortices.• The turbulence generated was found to be bi-modal for tumble motion in a pent-roof configuration. The initial vortex was found to breakdown at approximately1100 BTDC, this led to reduced turbulence generation around the top dead centreposition.• It was found that the influence of tumbling motion on turbulence integral lengthscales was negligible. The magnitudes of the turbulence length scales were observedto be generally smaller than those measured in a disc chamber of comparable clearance height.• When the effects of noise were corrected for, the turbulence intensities and integrallength scales show that the resulting turbulence structure was neither homogeneousnor isotropic.Chapter 2. LITERATURE REVIEW 2O• A strong correlation was observed between the corrected turbulence intensity at• the spark plug at ignition and the rate of early combustion. The ignition- 1% andignition- 10% mass burn fraction were found to have a positive correlation withturbulence intensity.• The experimental data showed that there was no apparent correlation between therms of the mean flow cyclic variations and the combustion cyclic variability.Arcoumanis et al (1990) conducted an experimental study to investigate the abilityof certain induction systems to enhance turbulence levels at the time of ignition, throughformation of long-lived tumbling vortices on the plane of the valve and cylinder axis.A horizontally-opposed twin-cylinder Citroen 2CV gasoline engine was employed forthe study. The inlet port was first inclined 90 degrees to the conventional directedport arrangement to generate pure tumble along the plane passing through the centresof the inlet and exhaust valves. In the second part of the investigation, the port wasinclined 45 degrees to the cylinder head in order to generate both tumble and swirl. Thelatter configuration represents an intermediate stage between the zero degree port which• produces a standard swirl and the 90 degree orientation which generates pure tumble(zero swirl). The engine has a compression ratio of 7.4 and was run at 1000 rpm. TwoLDV arrangements were used to measure the three velocity components.The conclusions from this study could be summarized thus:• The intake port generating pure tumble was found to produce 42% turbulenceimprovement over the directed port.• The probably ffiore practical configuration of combined tumble and swirl generateda 24% turbulence enhancement over the directed port intake assembly.• The 45 degree intake port arrangement resulted in a flow having a swirl ratio of 1.0Chapter 2. LITERATURE REVIEW 21which was found to be slightly less than the 1.2 swirl ratio produced by the directedport. Swirl ratio is an expression of the torque exerted upon the flow-straighteninghoneycomb of a swirl meter, whose axis is parallel to that of the engine cylinder,normalised by engine speed.The implications of pure tumble and tumble/swirl combination for combustion werealso discussed. It was stipulated that generation of a strong tumbling motion means thatthe port be oriented as close as possible to the plane passing through the centres of thevalves. This implies that there is less chance of any motion present in the tangentialplane surviving to the time of ignition. The stronger the tumbling motion, the moreturbulent kinetic energy is released during its break down and this release takes placelater during compression relative to a weaker motion. These two factors were thoughtto be capable of ensuring higher turbulence levels during the time of ignition but notnecessarily a proportional enhancement of the burning rates since during the combustionperiod turbulence will be in a state of continuous decay in the absence of any othermechanism of turbulence production.It was also stated that, although the combination of tumble and swirl as generatedby the 45 degree port produces inferior turbulence levels at the time of ignition whencompared to the pure tumble motion, they may allow for additional turbulence productionby swirl after top dead centre and or enhancement of the burning rate through strand soonhing of the flame.The authors indicated that, in spite of the fact that in disk-type chambers a solid bodytype of swirl may not generate more turbulence than that persisting through intake, pentroof chamber geometries with four valves might offer more promise for achieving fasterburning rates through;• Generation of tumbling vortices, and through them enhanced turbulence at theChapter 2. LITERATURE REVIEW 22time of ignition, and,• Distortion of swirl by chamber geometry at top dead centre and release of turbulentenergy at a time when the tumble-generated turbulence has started decaying tolevels insignificant to affect combustion.On the basis of the observations stated above, a hypothesis was propounded which statesthat, the combination of tumble and swirl may offer more promise as a mechanism forenhancing burning rates in spark-ignition engines.Finally, it was recommended that a frequency analysis of the flow be carried out toconfirm whether there are any quasi-periodic structures affecting the ensemble-averageinterpretation of the flow and to examine the cycle-resolved stability of the breakdownprocess of the stronger versus the weaker tumbling vortices, which may profoundly affectthe operation of lean burn engines.Kent et a! (1989) undertook an experimental investigation with the aim of studyingthe effects of intake-generated flow fields and the resulting combustion characteristics, ina reciprocating piston water analogue flow apparatus and in a firing engine respectively.Three different cylinder head/port geometries capable of generating six different levelsof tumble and swirl motions were employed in the study. For the flow visualization experiments, water containing suspended flow tracer particles (0.5 mm diameter polystyrene,specific gravity = 1.03) was substituted for the normal working medium — air.The major results from the investigation can be summarized as follows:• The correlation between increasing swirl and decreasing burn duration, which hasbeen observed in other studies, was seen to hold in this investigation as long asswirl was the major component of the large-scale motion. It was also observedthat the correlation between burn duration and swirl ratio was better when thepiston-driven water analog swirl data were used rather than those obtained fromChapter 2. LITERATURE REVIEW 23steady-flow swirl rig.• The burn duration decreased as the bottom dead centre tumbling strength increased.• The in-cylinder flow fields at bottom dead centre of intake were found to exhibitthree dimensional characteristics. For example, tipped swirl was observed. It wastherefore speculated that this observation was a more general case encompassing• both swirl and tumble motions. It was then inferred that, in addition to the strengthof vortical motion, the angle between its axis and the cylinder axis (zero for swirland 90 degrees for tumble) may be an important parameter of intake-generatedfluid motion for influencing the combustion rate in homogeneous charge engines.The first conclusion from this study reveals the inherent flaw in using steady-flowswirl results as a basis for deducing swirl ratios in engine studies.2.3 SummaryIt is generally very difficult to assess the findings of different investigators on a comparative basis, since the operating conditions considered are as vast as the engine types(configurations) and performance parameters of interest. However, there appears to bea common agreement on the following points:• The velocity and spatial velocity gradients during the intake process in valved engines vary considerably from one point to another within the combustion chamber.This suggests that the intake event is highly anisotropic and non-homogeneous andtherefore cannot be studied by using a single point measurement. Instead a largenumber of spatial data points are necessary.Chapter 2. LITERATURE REVIEW 24• The intake process is the major source of turbulence generation in simple chamberconfigurations, Clerk (1921). This suggests that it will be appropriate to studyamong other useful performance parameters, the effect of intake geometry (forexample, swirl generating ports, tumble inducing ports, swirl/tumble generatingintakes and so on) on TDC turbulence.• The turbulence generated during the intake process decays rapidly during compression for simple chambers, Clerk (1921). For complex chamber configurationsthe turbulence is somewhat enhanced during the compression stroke (see Amann(1985)). The inference that could be drawn from these observations is that it wouldbe a worthwhile exercise to investigate the effect of chamber configurations (for example, divided chamber, bowl-in-piston, bowl-in-head, wedged-chamber, and so on)on top dead centre turbulence intensity.• In disc chambers, as top dead centre position is approached the turbulence intensitytends to be uniform in space (homogeneous) and insensitive to direction (isotropic).This profound observation shows that it is sufficient to characterize TDC turbulence intensity by a single measurement at the plane of interest (that is within theclearance volume). The observation above refers only to the TDC turbulence intensity. In order to fully realize the notion of homogeneous turbulence, the uniformityof turbulence length scales must be obtained from experimental studies. This requires a direct measurement of the spatial correlations, rather thanattempting todeduce them from Taylor’s Hypothesis, Taylor (1938) that is, via time scales andmean convective velocity. The length scales calculated from Taylor’s Hypothesiswill most often be in error, since the ratio of turbulence intensity to mean convectivevelocity in an ICE is always greater than 0.2, (that is, u’/U 0.2). For Taylor’sHypothesis to be valid, the ratio should be less than 0.2, (that is, u’/tf < 0.2).Chapter 2. LITERATURE REVIEW 25• The instantaneous amplitude of turbulence intensity tends to be proportional tothe engine speed. This has been found to be true for both pancake and complexcombustion chambers, Liou (1983).• In pancake combustion chambers, it has been observed that near TDC of the compression both the mean velocity and turbulence intensity are nearly independent ofthe compression ratio.• The integral time scale of turbulence near TDC is approximately of the order ofone millisecond in the range 1000- 2000 revolutions per minute and the microtime scale is approximately one order of magnitude smaller. With these integraland micro time scales, some workers have computed length scales of turbulenceby using Taylor’s Hypothesis. They found that the computed integral length scalenear TDC is of the order of clearance height and that the ratio of the integral to themicro length scales was approximately 3 - 5, Lancaster (1976), Dent and Salama(1975), Tabaczynski (1976).• The number of experimental studies that has been conducted using the same enginefor both motoring and firing over a broad range of engine speeds and locations hasbeen very low. However, the paucity of data that are available from simultaneousstudies of motoring and firing operating conditions (in the same engine), suggestthat for engine speeds up to 2700 rpm, the following conclusions can be drawn:1. The gas velocities under firing conditions are greater than those measured for motored operations under the same conditions of intake arrangements.2. Close to the TDC position, but prior to ignition, the motored and fired gas velocitieswere noticeably very similar, but were found to differ considerably only when theflame was near the measurement point (see Brandl et at (1979)).Chapter 2. LITERATURE REVIEW 263. During the expansion and exhaust strokes, substantial differences were noticedbetween measurements taken in motored and firing engines.The second conclusion stated above is worthy of particular note, since it implies thatthe combustion during a cycle has little or no effect on the gas mixture velocity duringcompression in a subsequent cycle. It can then be argued that irrespective of the modeof measurement ((LDV) or Hot-Wire Anemometer (HWA)) employed to take data in amotored engine, such measurements will adequately reproduce the conditions prevalentin a firing engine during compression prior to ignition.All measurements that have been made to date in motored engines have thereforerelied on this profound conclusion to justify the extrapolation of motored turbulencemeasurements to the performance of firing engines. From these types of investigation,Lancaster et al (1976), Lancaster (1976) concluded, that there is a positive correlationbetween turbulence intensity and flame speed in valved engines.• Intake-generated swirl flow is very important in many types of diesel and gasolineengines. This conclusion is supported by some of the studies already reviewed andthe following:1. Brandl et al(1979) found that in small direct-injection diesel engines, swirl increasesthe rate of fuel-air mixing, thus leading to a reduction in combustion duration atretarded injection timings.2. In a study conducted by, Arcoumanis et al (1983), Gosman and Johns (1978), itwas concluded that swirl interacts with the compression-induced squish flow tomodify the flow structure, thereby increasing the levels of turbulence within thecombustion bowl.Chapter 2. LITERATURE REVIEW 273. Sung and Patterson (1982) found that in large uni-flow scavenged diesel engines,the swirl flow not only enhances the combustion rate, it aids the scavenging processas well.4. In a study conducted by Nagayama et al (1977), it was found that the effect ofswirl in gasoline engines tended to extend the air-fuel ratio at which misfire occursand Witze and Vilchis (1981) found that it reduced cyclic variability and improvedfuel economy (via enhanced operability with lean/dilute mixture).• Intake-generated tumbling motion (also called vertical or barrel swirl) has beensuggested as being an alternative to conventional swirl as a way of improving therate of cyclic stability of combustion in homogeneous charge engines under diluteoperating conditions. The impetus for this appears to be a study, Gosman etal (1985) in which enhanced turbulence at top dead centre of compression wasobserved for motoring operation when tumbling was produced by a central shroudedvalve. Similar observations were also recorded by, Arcoumanis et al (1990).From the general conclusions stated above, it becomes obvious that a thorough knowledge of the mechanics of top dead centre turbulence generation, swirl and tumble formation and their influence on engine performance is highly desirable.2.4 Numerical Simulation StudiesThe findings from various experimental studies mentioned above strongly suggest thatcombustion and flame propagation inside the engine cylinder are very complex and intimately dependent on turbulent flow fields. It therefore follows that parametric studiesof engine flow fields through experimental methods will not only be time consuming butvery expensive as well. It is this realism that has informed the use of mathematicalChapter 2. LITERATURE REVIEW 28models in guiding the development of internal combustion engine flow dynamics. Sincethe main objective of this section is to review some modelling works featuring the effectsof swirl and squish on top dead centre turbulence, other aspects, such as the analysis ofspray and combustion dynamics, will not be examined any further.Numerous investigative studies have been carried out with the k- e turbulence model.These studies involve the generation of cold- and combustion-flow predictions and comparison with experimental data for both two-dimensional, Li et al (1989), Fukumoto etal (1985), Abraham et al (1985) and three-dimensional, Naitoh et al (1990), Haworth etal (1990), Le Coz et al (1990), Henriot et al (1989), Errera et al (1988), Abraham andBracco (1988), Amsden et al (1985) cases. It is worth mentioning however, that theseclassification apply more rigorously to the calculations than to the experiments, since theassumption of two-dimensionality was not confirmed in most of the experiments.Some works drawn from recent ICE studies where particular attention was paid tothe choice of data and in minimizing as far as possible the prediction errors from sourcesother than those due to turbulence modelling, initial and boundary conditions are nowbriefly reviewed.Le Coz et al (1990) conducted an investigative study with the aim of determiningthe flow fields generated by different types of intake configurations. The study wassupplemented with a three-dimensional calculation, so as to aid in data interpretation.Three intake configurations (called A, B and C) differing in their valve layout were chosenfor the study. In configuration A, both inlet valves operate. In configuration B, one of thevalves was blocked and the other shrouded in order to generate a high level of swirl. Theintake valve was equipped with a 120 degrees valve mask, oriented towards the nearestcylinder wall. The last configuration, C, is similar to B except that the inlet valve wasunshrouded, and therefore gives a lower level swirl. The swirl ratios of the three engineheads, defined as the equivalent angular velocity of the cylinder flow field divided by theChapter 2. LITERATURE REVIEW 29engine speed, were determined on a steady-flow bench. They were found to be 0.0, 2.5and 1.0 for configurations A, B and C respectively. The data acquired were analyzed bycycle-resolved techniques.The predicted results were obtained by performing a three-dimensional numericalcalculation with the computer code KIVA. Details about the code are given in chapter4 of this thesis; nevertheless, some of the techniques employed to carry the calculationsthrough are briefly described. The authors successfully developed a new mesh generatorand an interface to exchange data between KIVA and their grid generator. The intakeprocess was computed by using a moving valve, which was modelled as a moving obstacleopening and closing the inlet pipe. The valve trajectories and velocity were set to theactual values.The initial gas temperature and pressure were assumed to be homogeneous and theseparameters were estimated by using a two-zone zero-dimensional code. The initial valuesof the turbulent kinetic energy, k, and its dissipation rate, 6, were assumed to be homogeneous. The turbulent kinetic energy was set to 10% of the kinetic energy based onmean piston speed while its dissipation rate was set k3/2h, where h is the half distancebetween top of the cylinder head and the piston head at TDC. The turbulent boundarylayer was modelled with a law-of-the-wall assuming a logarithmic velocity profile. At theintake port, the velocity profile was assumed to be uniform across the cross-section ofinterest, and the mass flow rate was calculated from the pressure difference between thecombustion chamber and the upstream conditions.With all these modifications/assumptions, the calculations were made from intakeTDC to compression TDC, for the same set of conditions as those of experimentalengine. The important findings resulting from this work are summarized below;• In the two configurations for which only one intake valve operates one inclinedChapter 2. LITERATURE REVIEW 30tumbling motion was generated throughout the combustion chamber. In the basicconfiguration, for which two intake valves operate, the experimental data werefound insufficient to adequately describe the resulting motion, due to the verycomplex nature of the flow field. The results obtained from the numerical simulationwere employed to facilitate the interpretation of the resulting flow field. From thecomputer simulation it was found that the flow field is symmetrical about themedian plane of the intake valves. In each half-chamber the flow structure wasobserved to be close to an incline tumble. These inclined tumble motions wereobserved to interfere with one another in the symmetry plane.• The turbulence structures generated within the combustion chamber were tracedto two sources, iamely; those due to the intake velocities and those resulting fromthe break down of tumbling motions at the end of compression.• The experimental data were cycle-resolved and from the results, it was deducedthat the ensemble fluctuations were for the most part situated in the low frequencydomain. It was further shown that the ratio between the low frequency and thehigh frequency contributions were practically constant for a given cut-off frequency.Based on these observations, it was concluded that when turbulence generationoccurs, all the frequencies are amplified including the large- and small-scale fluctuations. From this conclusion, it was inferred that an optimal turbulence level mayexist that is suitable for operating an engine on a lean mixture.• The spectral analysis of the instantaneous velocity showed that the turbulencestructure is roughly in an equilibrium state, the spectra being almost self-similar.This implies that the turbulent kinetic energy is always distributed in the samemanner in the frequency domain.Chapter 2. LITERATURE REVIEW 31• The horizontal velocity components from both the experimental and numericalcalculations were compared. The agreement was found to be very poor. This wasattributed to numerical truncation errors. However, the general trends of the fluidmotion were correctly predicted.A very important observation worth noting, is in the area of data comparison. The authors compared cycle-resolved velocity components with predicted velocity values. Thisis an inconsistent comparison, because the predicted values are based on ensemble averaging techniques, whereas, the experimental values are not — they are cycle-resolved values.The philosophy underlying both techniques are fundamentally different and unrelated.Bizhan et al (1989) conducted a detailed multidimensional study of the influence ofintake port design on the in-cylinder flow structure during the intake and compressionstrokes, the mixing of residual gas and a non-premixed intake charge and the extent andpattern of charge inhomogeneity near the time of ignition.The engine simulated in the investigation is typical of the current lean-burn designand the inlet ports studied are:• An idealized directed (zero-swirl) and,• A helical (swirl-producing) port.The FMCS numerical code was used for this study. The code solves the densityweighted, ensemble-averaged differential conservation equations of mass, momentum andstagnation enthalpy, in addition to transport equations for the turbulence kinetic energy,k and its dissipation rate, e. In this study two additional differential equations for the fuelmixture fraction and the residual gas mixture fraction are solved in order to investigatethe in-cylinder mixing process.In order to accurately represent curved boundaries, the transport equations mentionedabove were recast into their general curvilinear orthogonal form and to allow for theChapter 2. LITERATURE REVIEW 32movement of boundaries (for example, the piston and valves) an Eulerian-Lagrangianco-ordinate transformation was applied.The differential equations were discretized by the finite-volume method. The Euler-implicit temporal and hybrid of upwind/central spatial differencing schemes were employed in order to ensure numerical stability. The solution is advanced in time (crankangle) by solving the system of algebraic equations in an iterative fashion at successivetime steps.The thermodynamic, initial and the boundary conditions immediately upstream ofthe valve, in addition to the instantaneous mass flow rate during the induction stroke,were obtained from unsteady gas dynamic calculation.The characteristics of the intake port was imposed on the intake flow by prescribingthe spatial velocity distribution around the valve periphery and within the gap, throughout the intake. A circumferentially uniform plug flow velocity distribution, zero swirl anda flow entrance angle equal to the valve seat angle was assumed for the idealized directport. In the case of the helical port, LDV measurements of the velocity distribution wereemployed to prescribe the velocity boundary conditions.Calculations were commenced at 112° crank angle and terminated at 3200 crank angle.The initialization time of the calculation was chosen to correspond to inlet valve closurewhile the termination time was selected to coincide with the initiation of combustion. Acomputational time step equivalent to 1/16 crank angle degree was used.The important findings from this investigation were:• It was found that the structure and evolution of the in-cylinder flow field is dependent on the flow characteristics at the inlet valve and the cylinder geometry.• The flow in the cylinder at the time of ignition (3200 crank angle) was found to bepredominantly characterized by the evolution of the intake generated flow.Chapter 2. LITERATURE REVIEW 33• The flow structure in the piston-bowl was found to be dependent on the interactionbetween the in-cylinder flow structure and squish motions during the compressionstroke.• The numerical simulation showed that the flow in the piston-bowl at the time ofignition was still not being governed by the dynamic forces (that is, squish, swirlmomentum transfer, and so on) and that it exhibited some memory influences.The implications of the findings stated above with respect to the simulation of in-cylinder flow processes were summarized thus:1. Studies that attempt to simulate only the compression process, by starting thecomputations with assumed and grossly simplified initial conditions at inlet valveclosure, are very likely to lead to erroneous conclusions.2. With regards to the initialization practices, particularly those employed in thisstudy, it was recommended that the validity of the assumptions and the extent oftheir influence be carefully examined and if possible alternative practices sought.The findings concerning the in-cylinder mixing process and its dependence on thecharacteristics of the intake flow are:• The flow characteristics at the intake valve strongly influence the intake flow structure and the mixing process.• The intake swirl promotes mixing by enhancing the interaction of the in-cylindercharge with the intake jet and the primary recirculations.• In the absence of the above interaction, intake swirl tends strongly to preserve theaxial charge stratification (for instance, after the completion of the intake process).Chapter 2. LITERATURE REVIEW 34• The extent of charge inhomogeneity at the time of ignition was found to be similarfor the direct (zero swirl) and helical (swirl) intake ports, but the distributions weresignificantly different. The charge distribution was found to be markedly irregularfor the direct port and axially stratified in the case of the helical port.• It was found that with the present fuel injection timing, the direct port providesa better fuel distribution, however, the expected susceptibility of the charge distribution to small changes in the flow boundary condition may produce large cycle-to-cycle variations of combustion.• From the calculations it was inferred that the helical port with a late (duringintake) fuel injection timing may provide a stable and favourably-stratified chargedistribution for combustion.• The study also revealed that, in investigations of the influence of inlet port designon the combustion characteristics, it is important to separate the influence of theflow structure from its charge preparation characteristics.It is however, pertinent to mention that the various conclusions/findings reviewedabove be taken with caution, since the code has not been validated against any set ofreliable experimental data.Uzkan and Hazelton (1986) carried out a theoretical study with the aim of studyingthe effects of engine swirl level on the mixing of the fresh charge with the residual gases• during compression in a uniflow type two-cycle diesel engine geometry.The CONCHAS-Spray code, Butler et al (1979) was used in this study. The valveswere treated as annular slots on the head with valve opening simulated through the use ofcoefficients as a function of valve lifts. Subsonic flow was assumed at the port and valves.Inlet ports were simulated as circular openings on the liner with effective area controlledChapter 2. LITERATURE REVIEW 35by piston motion. The compressible flow equations with constant flow coefficients wereused to calculate the flow rates to and from the cylinder. Back-flow to the air box wasallowed at the inlet port. During inflow at the port, the flow vector direction was assumedto have a given effective angle with respect to radial direction and this was named theport angle. The effect of engine swirl on the flow field characteristics were studied bychanging the effective port angle within 10 to 30 degrees range at intervals of 5 degrees.The initial conditions of the radial, axial, swirl velocities and the turbulence lengthscale were prescribed. The initial thermodynamic states of the working fluids at the airbox, exhaust manifold and in the chamber were predicted through a two-zone thermodynamic engine simulation model.The computations were started at the exhaust valve opening and continued up to thetop dead centre position.The important conclusions arrived at in the study were:• Residual gases do not mix completely with the fresh air charge during the scavenging and compression period. It was inferred that a considerable stratification offresh charge/residual gas exists at top dead centre.• It was found that the inlet port angle was an important parameter affecting theprocess of mixing the residual gas with the fresh air. It was also found to havea considerable effect on both the magnitude and the stratification pattern of thefresh air charge distribution at top dead centre.• The fresh air charge stratification patterns at top dead centre were found to bedifferent for various inlet port angles.• The calculated fresh air and residual gas conditions at top dead centre for the threedifferent port angles indicated that the trapping and scavenging efficiencies as wellChapter 2. LITERATURE REVIEW 36as the 10%, 90% and global fresh air concentrations increase with an increase inthe port angle.• A new vortex pattern named Head Vortex which was generated by the swirl velocityfield was identified, and a mechanism for its formation was also suggested.The work under review attempted to demonstrate that it is possible to study thecomplex three dimensional, time dependent intake-generated flow structure prevailing inan engine configuration having a contemporary inlet port design, with a two-dimensionalnumerical code. A number of questions remain unanswered as to the validity of thek-c turbulence model for this type of flow processes and the spatial accuracy of thenumerical solution, considering the fact that a two-dimensional code was used to studya three-dimensional process.It is however, necessary to mention that, the final assessment of the accuracy of theresults reviewed above must be based upon detailed comparison with reliable sets ofexperimental measurements.Chapter 3EXPERIMENTAL AND DATA ANALYSIS METHODIn this chapter a detailed description of the experimental system and the procedurewhich was followed to make the reported measurements are described. A detailed description of the data acquisition system and laser optics can be found in Dostek and TSIUser’s Manuals respectively. Also included in this chapter are a review of various datareduction methods and the factors considered in selecting the data analysis technique ofthe present study.3.1 Experimental SystemThe main equipment selected for the realization of the objectives of the present studyis the UBCRICM. This machine was selected because of its ease of use, coupled withthe fact that desired operating conditions can be closely achieved — this is vital for thesimulation phase of the present study, since one of the goals of the experimental portionof this thesis was to setup boundary conditions which closely mimic those in KIVA-Ilcode. This has the advantage of providing experimental measurements which can bemuch more closely related to KIVA-Il generated results than in studies where arbitraryboundary conditions are invoked.The UBCRICM was designed and built by Dohring (1986). A detailed description ofthe machine can be found in Dohring (1986). The UBCRICM was built for the purposesof studying the comparative effects of intake and compression stroke generated turbulenceon the combustion duration in a spark ignition engine. The UBCRICM can be operated37Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 38in either of two modes, namely, with both intake and compression strokes active or withonly the compression stroke in operation. In either mode of operation, the UBCRICMoperates at a maximum simulated engine speed of 1000 revolutions per minute, withthe piston starting and stopping at rest, this ensures that combustion takes place atconstant volume. In the two stroke mode, the piston starts at TDC before intake andthe intake mixture is drawn from a small tank which is initially at ambient temperatureand pressure. Two stroke mode was employed in the present study.The UBCRICM is driven by compressed air acting on a pneumatic cylinder on oneside of the rack and controlled through a hydraulic ram on the other side of the rack.The rack is meshed with a pinion connected to the crank shaft which drives the piston.A schematic diagram of the UBCRICM is shown in Figure 3.1.The rapid acceleration of the rack resulting from triggering a fast response solenoidvalve is kept constant by a constant flow rate orifice on the hydraulic side of the machine.A braking pin, which blocks the flow through the orifice, was designed to bring the rack toa rapid stop at the end of the compression stroke without a destructive collision betweenthe ram and the orifice. A slot in the rack was milled so that a pin riding in the slotwould control the intake valve through a linkage mechanism.The piston, cylinder and head assembly were designed to be readily changeable, sothat various combustion chamber configurations could be easily realized. The pistonwas designed in such a way that the crown could be easily screwed-in or unscrewed.The bowl-in-piston and re-entrant combustion chamber configurations were used in thepresent study. The cylinder was made from a pipe available off-the-shelf and micro-honedon the inside surface. The head assembly was machined from mild steel and no coolingwas required, since the engine does not run continuously. The combustion chamberspecifications and the UBCRICM operating parameters are given in Table 3.1.A soft pot Si series optical shaft encoder was used to monitor crank shaft angularChapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 39position. The output from the encoder was sent to the Laser Velocimeter Interface1400A (LVI 1400A) data acquisition system. The encoder had one crank angle degreeresolution which can be improved to a half-degree resolution by setting the relevant• jumper accordingly, but one-degree crank angle resolution was found to be adequate forthe present study.3.2 Modifications made to the UBCRICMAs mentioned previously, the main equipment selected for the realization of the objectivesof the present study is the UBCRICM. For the UBCRICM to be suitable for the LDVexperimental studies,. some modifications had to be made. The relevant modificationsare as follows:3.2.1 Intake ArrangementThe intake system was redesigned in order to allow for a ported induction system. Amethod of generating the required level of swirl motion in the combustion chamber of theUBCRICM was designed and installed on the machine. The design entails drilling fourreinforcing tangential ports on the wall of the combustion chamber. This arrangementhas the advantage of closely duplicating the initial conditions assumed in KIVA-Il code,• that is, uniform intake flow quantities, and moreover, there is no volumetric efficiencypenalty incurred. Intake air was supplied from a compressed air bottle. A schematicdiagram of the intake arrangement is shown in Figure 3.2 and the subsequent interaction(idealized) of swirl and squish is shown in Figure 3.3.Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 403.2.2 Optical WindowIn order to facilitate data acquisition by Laser Doppler Velocimetry (LDV) method, thecombustion chamber head was redesigned to accommodate a quartz glass window. Thequartz glass has a thickness of 2 inches (50.8 mm) and a diameter of 4 inches (100.0mm). With this set-up a full field (100%) view of the combustion chamber was obtained.This represents one of the largest in-cylinder optical views that has been reported for aninternal combustion engine experiment. The generous optical view was realized, becausevalves were not located on the cylinder head.3.3 Optical SetupThe back-scatter mode of the LDV data acquisition system was used for the reportedmeasurements. One of the commonly cited disadvantages in back-scatter set-ups is thatthe beams which enter the engine through the head window are diffusely reflected fromthe top of the piston as it approaches TDC. This is believed to cause increasing noisein the LDV measurements. This shortcoming was eliminated by painting the top of thepiston flat black.When optical access is limited, (as it always is in the engine environment) and oroptical alignment becomes critical, the back-scatter method becomes the preferred optionsince optical access is required from only one side. Furthermore, the same lens is used forcrossing, focussing and collection of the scattered beam, this makes the back scatter modeessentially self-aligning. With the most serious drawback of the back-scatter methodremoved, the benefits were then explored to the fullest in acquiring data. A 500 mm lenswas used to cross, focus and collect the scattered beam.Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 413.4 Data Acquisition SystemA high-performance bOS-Tek IBM-PC compatible Laser Velocimeter Interface (LVI)model 1400A, specifically designed for LDV data acquisition was used to acquire datain the present study. The basic LVI model 1400A interface is a single full-size IBM-PCcompatible circuit board (card). It acquires data from a single LDV signal processorwhile simultaneously recording the angular position of each event (trigger). A detaileddescription of the LVI model 1400A can be found in Dostek user’s manual. A schematicdiagram of LVI model 1400A is shown in Figure Experimental ProcedureTwo components of velocity (radial and tangential) measurements were made in thenon-firing UBCRICM. Both velocity components in planes (called A, B and C) of thepiston crown were obtained by simply rotating the beam splitter through 90 degrees. Thereported measurements were made at a simulated speed of 1000 rpm, and at various radialpositions and three axial planes. Two intake swirl levels (0.5 and 2.0) were investigatedcoupled with the effects of squish. The back-scatter LDV mode of operation was usedfor all experimental conditions. A schematic diagram showing the complete engine-LDVsetup is shown in Figure 3.5. The optical components of the system, including the 2Watt Spectra-Physics laser were mounted on an optical table supplied by TSI Inc.The laser beam was incident into the modular optical units, which essentially consistedof prisms and lenses. The laser beam then traversed a number of optical modules whereit was split, frequency shifted, colour separated, expanded and aligned. The two parallelbeams emerging from the optical arrangements were steered to the engine test sectionby two mirrors mounted on a vertical frame attached to the optical table. A 500 mmfocusing lens, mounted on a graduated traversing unit was used to cross and focus theChapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 42beam to the desired measurement location. The dimensions of the measurement volumewere of, length, 1m = 2.68 mm, diameter, dm = 0.235 mm and the fringe spacing, f3, was0.295 sum.Since LDV does not measure the velocity of the fluid itself, but only determines thespeed of scattering particles, it therefore follows that such particles must be added to theintake air of the engine. For this reason a specially designed canister was used to storethe magnesium oxide (MgO) particles used for seeding. The particles had a nominaldiameter of 2 microns. In order to remove the water vapour in the intake air, so asto avoid the danger of condensation on the walls and window and agglomeration of theseeding particles, an air drier was installed in the compressed air supply line, upstreamof the canister holding the seeding particles. Dried compressed air was passed throughthe canister and it entrained particles on its way to the combustion chamber. A steadyflow method was used to estimate the average intake swirl level.During experimentation, the scattered light was reflected backward along the sameoptical path as those of the transmitted beam (backscatter). A pinhole located in theback scatter unit acts as a spatial filter by removing or reducing the background noise,thereby improving signal quality. The filtered back scattered light was then focused onto a one micron pin hole and finally passed through to the photomultiplier, which thentransformed the light intensity into an electrical signal. The signals were band-pass filtered to remove low frequency pedestal and high frequency electronic noise, amplified andtransmitted to a TSI frequency counter through a gating circuit after validation. Theywere also checked for modulation depth on a Nicolet 3091 digital storage oscilloscope.The 16-bit digital outputs of the counter and the crank angle information (resolution1.0 degree crank angle) were transferred to an IBM-PC compatible digital computer.A special LDV DOS-Tek interface commercially available and described extensively inthe Dostek user’s manual was used for the data acquisition, analysis and transfer to theChapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 43computer.A graphics monitor and a Nicolet 3091 digital storage oscilloscope were also connectedto the whole set-up. These peripherals were very important for rapid data display,monitoring and analysis.3.5.1 Techniques Used to Ensure Data QualityThe following tests were performed to check the quality of the data being acquired:• Just before data acquisition, one of beams was blocked to verify that the data-ratefell to zero. The blockage of one of the beams meant that no fringe patterns wereformed at the measurement volume position.• Operation of the system and counter set-up were check by confirming that no-particle-seeding corresponded to zero data rate.• Measurements were taken using different seeding densities to ensure that the seedingrate was not affecting the flow.• While taking measurements, the signal quality of the Doppler burst was monitoredon a storage oscilloscope to ensure that counter gain and photomultiplier voltageswere not set too high. Excessively high settings resulted in Doppler bursts withirregularly shaped envelopes and jagged oscillations.• In order to verify that the data acquisition card was taking data accurately, a knownconstant-frequency (50 kHz) sine-wave signal was applied to the photomultiplierinput and sampled by the card via the signal processor counter.Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 443.6 Data AnalysisThere are four categories of data analysis methods presently being used in characterizing turbulence intensity and mean velocity from measurements taken in ICEs. Thesetechniques are ensemble-averaged, cycle-resolved, conditional and autocorrelation methods. These techniques are reviewed below and considerations leading to the choice ofthe preferred option are stated. Before describing the procedures any further, a reviewof the underlying issues hampering a general characterization of in-cylinder turbulenceis discussed.3.6.1 IssuesIn general, it would be preferable to describe a turbulent velocity field in terms of itsspatial characteristics. However, available techniques to measure simultaneous two-pointinstantaneous spatial-velocity distributions in engines are still in their infancy, althoughsome degree of mean-velocity information has been obtained, Marko et aX (1986), Bates(1988). Quantitative turbulence information are presently unavailable. Direct length-scale measurement using multi-point LDV is also in a comparatively early stage of development, Ikegami et aX (1987), Glover et al (1988), Corcione and Valentino (1991).Single-point velocity-crank angle data are more readily acquired. Even in steady flows,these single-point measurements are very limited in their ability to adequately characterize the velocity fluctuations.The flow field within the cylinder is non-stationary due presumably to the dynamicallychanging boundary conditions resulting from piston and valve motions. Other factors thatcould be responsible for the unsteady in-cylinder flows can be attributed to fuel-injectionand combustion in fired engines. The frequency spectrum of the resulting unsteady meanor bulk flow is therefore expected to overlap the low-frequency end of the spectrum of theChapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 45turbulent velocity fluctuations. Furthermore, the in-cylinder flow field may also be actedupon by various instabilities, such as intake- or squish-jet, swirl precession, Arcoumaniset al (1987). These are often considered to produce a. mean or bulk flow field that isnot statistically repeatable from cycle to cycle or that involves deterministic motionsthat are not in phase with the piston motion. It is however, pertinent to mention thatthe experimental results of Abdel-Gayed et al (1984), do not support these suggestions.It is often argued that, for single-point velocity measurements, such processes lead tofluctuations (about the many-cycle ensemble average) which should not be regarded asturbulence in the classical sense of small-scale random motion superimposed on a welldefined, deterministic mean flow. Even if one does not fully subscribe to this point of view,it is obvious that large-scale, low-frequency processes will affect mixing and combustionvery differently from small-scale, highly chaotic fluctuations. The ensemble-rms velocity-fluctuation intensity, however, includes all deviations from the ensemble-mean velocity,no matter what their origin, spatial/temporal scale, or degree of deterministic character.Cycle-resolved velocity data are usually acquired with the hope that they can lead to amore comprehensive description of in-cylinder processes, particularly true turbulence.This expectation is far from being universal. An alternative point of view is thatconventional ensemble averaging provides the appropriate characterization of the flowfield, apart from quasi-periodic systematic processes such as vortex shedding or swirlprecession, Gosman (1985b).3.7 Review of Data Analysis MethodsThere are four basic approaches for extracting quantitative information from instantaneous velocity-crank angle measurements made in internal combustion engine. They areChapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 46the phase time-ensemble-averaging, cycle-resolved mean (bulk) velocity, conditional sampling and non-stationary velocity autocorrelation techniques. The cycle-resolved mean(bulk) velocity, conditional sampling, (this technique involves recording in permanentform the output from. one or more sensors thereby permitting fluid, thermodynamic andchemical information to be obtained. Subsequently the data are interpreted (processed),in a variety of useful ways subject to some conditions) and the non-stationary velocityautocorrelation techniques were founded on the philosophy of possessing the ability toaccurately distinguish between true turbulent velocity fluctuations and cycle-to-cycle orlow-frequency large-scale fluctuations in the mean or bulk velocity. In the discussions tofollow, these various methods of obtaining quantitative information from instantaneousvelocity-crank angle measurements are described. In section 3.8, considerations leadingto the choice of the preferred method of data analysis (ensemble averaged) of the presentstudy are presented.3.7.1 Ensemble-Averaging AnalysisConventional time-averaging of the measured instantaneous velocity is not applicable tothe unsteady turbulent flows encountered in internal combustion engines, where changesin the mean flow can occur on a time scale of milliseconds. The equivalent method wouldbe a phase average, whereby the instantaneous value at a specific crank-angle position isaveraged over many cycles of the engine. This type of phase-averaged quantification ofturbulent flow is appropriate for processing hot-wire anemometer measurements, wherethe response to the signal is continuous. However, this procedure is not suitable forprocessing LDV signals, because data are acquired at random intervals. Therefore, itbecomes necessary to use a phase-time average, where the time average is taken over afinite crank-angle width, 9.If the instantaneous velocity, U, at an engine crank angle, 6, is defined to be madeChapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 47up of a mean component Ti and a fluctuating component u, then;U(6, 6, i) = TT(6, zi6) + u(6, 6, i) (3.1)where, 6 is the crank angle of the measurement, LI6 is the width of the crank anglewindow in which the measurement was made, and, i is a particular engine cycle froma large number of cycles which make up the ensemble. The symbol U (unadorned byany sub- or super-scripts) denotes either of the two measured velocity vector components(radial or tangential) of the instantaneous velocity at crank angle, 6 in engine cycle, i.The overbar on any quantity represents phase-time average. The phase-time averagemean (ensemble average) velocity can be evaluated by the expression;1 N—1N1 / L6TJ(e, M) = U, (6±—i--) (3.2)t =o3i \ /where N1 is the number of velocity measurements recorded in the i’ cycle, N the number of cycles and N the total number measurements. The corresponding turbulencefluctuation intensity, u’(6, z6), is defined as;1 N—1N / Z61/2u’(6, t6) = u2j (6 ± (3.3)t i=Oj=1 \ /where,(e ± = U1, (s ± — TT(6, M) (3.4)Because the fluctuations are distributed randomly, their ensemble mean vanishes, thus;U1,2 (e ± = 0 (3.5)The essential drawback in estimating a phase average from a phase-time average is theresulting degradation in crank-angle resolution. For most fluid dynamic motions that occur within an engine, 10-degree resolution has been found to be adequate. The exceptionChapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 48occurs in a combusting engine, where crank-angle resolution on the order of 1-2 degreesis required. A 1.0 CAD width was found to be satisfactory for the LDV measurementsreported in this thesis.3.7.2 Cycle-Resolved TechniqueIn the cycle-resolved technique, Liou and Santavicca (1983), some form of low-pass filtering is used to estimated the mean velocity in each engine cycle. The rms deviations ofthe instantaneous velocities from the in-cycle mean velocities are then computed over allthe cycles and usually interpreted as a measure of the intensity of true turbulent velocityfluctuations. Conversely, the rms fluctuations about the ensemble-average velocity areusually considered as been made up of contributions from the trite turbulence and fromcycle-to-cycle variations in the mean or bulk velocity.As already stated above, the filtering approach to cycle-resolved-velocity analysisbegins by low-pass filtering the instantaneous velocity-crank angle traces in order toestimate the individual-cycle mean velocity, Ub(8 i). This is the same as stating that thefluctuations about the ensemble-mean velocity can be split into low- and high-frequencycomponents relative to a filter cut-off frequency, Fansler and French (1988). Under thisassumption, we have:u(8, 6, i) = Ujf(9, z8, i) + uhf(6,8, i) (3.6)This is equivalent to a three-term decomposition of the instantaneous velocity, as shownbelow:U(9, zSO, i)=TJ(8, z8) + ujf(O, z6, i) + uhf(S, z8, i) (3.7)The filtered or in-cycle mean velocity, Ub(O, i), is the sum of the ensemble-mean velocityand the low-frequency fluctuation term; that is,Ub(S, 8, i)=T(8, z6) + u,(8, tO, i) (3.8)Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 49The ensemble averages of the low- and high-frequency fluctuations should vanish individually since, by definition, the ensemble-mean velocity (8, M) already contains allvelocity contributions that are phase-locked with respect to crank-angle Therefore,lj1(8,9,i) =izhf(e,/.o,i) = 0 (3.9)Intuitively, the ensemble average of the in-cycle mean velocities will be the same as theconventional ensemble-mean velocity. This is readily evident from Equations 3.8 and 3.9,thus have:Ti(6, 9, i) = U(8, z8) (3.10)The intensities of the high- and low-frequency fluctuations are characterized by their rmsdeviations about an appropriate mean velocity, using relations analogous to Equations3.3 and 3.4. The high-frequency velocity-fluctuation intensity is the rms deviation of theinstantaneous velocity about the in-cycle mean velocity.N 21/2uf(8,M)=/hf(8,M,i)= [[-U(8±_i)_ub(8±Ti)]] (3.11)Similarly, the low-frequency velocity-fluctuation intensity is evaluated as the rms departure of the in-cycle mean velocity from the ensemble-mean velocity:1 N Ls.821/2u1(O,1X8) =.1,/jj(8,M,i) = [[>Ub (e±—u (o± .-_)] ] (3.12)The conventional ensemble-rms fluctuation intensity is related to the low- and highfrequency intensity through Equations 3.3 and 3.6, thus:u’(S) = J(uij(8, M, i) + uhf(O, M, i))2= /u(8, z.6) + uf(8, t9) + 2u,(8, tO,i)uhf(6t6,i) (3.13)The low-frequency and high-frequency fluctuations are usually assumed to be statisticallyindependent (see Fansler and French (1988)), so that the cross-term in Equation 3.13 isChapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 50zero, that is,Ujf(8,zS,i)Uh (S,X)= 0 (3.14)and,u’(G) = ,/u(6, 8) + u’,1(O 8) (3.15)This expression is not always true. An example of a situation where this expression isnot correct is when making measurements near the edge of a flapping intake jet. Theinnermost part of the jet will be substantially more turbulent than the outer sections,Arcoumanis and Whitelaw (1987). Therefore, the high-frequency intensity u1(,8)at a spatially fixed measurement point could depend on how low-frequency large-scalefluctuations affects the jet location. The high- and low-frequency fluctuations intensitiesuhf(8,i.x8) and u1(8, z8) are sometimes referred to as the turbulence intensity and thecyclic variation in the mean (bulk) velocity, respectively, Liou and Santavicca (1983).Although filtering schemes may possibly be useful, they can not rigorously decomposethe instantaneous velocity into mean velocity, cyclic variations and turbulent fluctuations.This is due to the following reasons;• The frequencies of the larger-scale turbulent velocity fluctuations must be expectedto overlap those of the unsteady mean motion and of any instability-related phenomena (flapping jets, precessing swirl, etc) that may be occurring in the combustion chamber;• Filtering cannot ascertain whether the velocity fluctuations are random or deterministic.3.7.3 Filter ImplementationThe low-pass digital filtering necessary to estimate the in-cycle mean velocity, (4(8, L6, i)for each engine cycle can be accomplished in a variety of ways. The simplest approachesChapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 51are a moving- or sliding-window average in the time (crank-angle) domain and a hardcutoff in the frequency domain, Liou and Santavicca (1983).The moving-window average is implemented thus: U&(, M, i) is estimated by taking alocal time average (uniformly weighted) over several adjacent values of the instantaneousvelocity, symmetrically disposed about 8. This averaging window is then stepped successively through the cycle. The length of the averaging window determines the effectivecutoff frequency.The frequency-domain filter is equally as complicated as the moving-window averagingtechnique. The first step in the implementation of the frequency-domain filter is to takethe Fourier transform of the instantaneous velocity-crank angle record, U(8, 8, i), thisresults in a series of Fourier coefficients representing its frequency content. The low-pass filter is then applied in the frequency domain by setting all the Fourier coefficientsabove the selected cutoff frequency, f, to zero and retaining the coefficients below, funchanged. Finally, the inverse Fourier transform of this truncated series back to the time(crank-angle) domain gives the desired estimate of the in-cycle mean velocity, Ub(8L8, i).The two filters just described represent opposite extremes in their time- and frequency-response characteristics. The moving-average filter’s hard cutoff in the time-domainproduces a relatively broad frequency-response curve with ringing or side lobes. On theother hand, the sharp-cutoff frequency-domain filter can lead to distortion or over-shootin its time domain output. However, some more elaborate filtering schemes can be usedto soften the hard cutoffs and reduce their undesirable features.3.7.4 Cutoff-Frequency SelectionThe success or otherwise of the frequency-domain filtering technique of data analysisdepends on the accuracy with which the cutoff frequency is specified. Presently, there isno algorithm to determine the cutoff frequency, f0, a priori. Ideally, the cutoff frequency,Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 52fe,, should be determined from the characteristic scales of the physical processes beingstudied, but in many cases (for example, turbulent flame propagation, jets, etc) these• processes are poorly understood. Due to insufficient knowledge concerning these practicalprocesses, various ad hoc methods are often used to determine the cutoff frequency. Someof the approaches that have been used to select cutoff frequency are;• Determination of cutoff frequency based upon the highest frequency in the ensemblemean velocity signal, Liou and Santavicca (1983), Fansler and French (1988). Theargument for this line of approach was based on the premise that if a frequencyexists in the ensemble mean velocity signal, it will also be present in each of theinstantaneous velocity signals, and moreover, that it should correspond to a partof the bulk velocity. Based on the foregoing, the highest frequency in the ensemblemean velocity signal is taken to be an estimate of the cutoff frequency. The majordrawback here is that some low frequency turbulence contributions from uhf(S, i)would be excluded. Too low a cutoff frequency will also lead to in-cycle mean-velocity estimates U(8, i) which are too smooth.• The cutoff frequency can be taken as the inverse of the decay time of the randomcomponent in the autocorrelation functions. In this approach, the autocorrelationwill have to be calculated, and the decay time (or integral scale) is defined as thefirst lag at which the autocorrelation function crosses the time or crank-angle axis.This approach is based on the assertion that the bulk velocity differs from thecycle-resolved turbulence because of its coherence characteristic. That is, the bulkvelocity is thought of as being deterministic, since it is a large-scale process whichis relatively well organized, whereas, the cycle-resolved turbulence is taken to be arandom, small-scale process generated by the shearing motion the larger scales.Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 53• To assign a cutoff frequency on the order of the engine cycle frequency; for example,the inverse stroke duration, that is, f0 = 2NE, where NE, is the engine speed. Inthis approach, the frequency content of the ensemble-mean velocity can extend tosignificantly higher frequencies than 2NE, as for example, during the passage offlame front or strong squish flow.3.7.5 Missing DataIn order to estimate the in-cycle mean (bulk) velocity, Ub(O, z6, i), each of the velocity-crank angle records must contain a single value of velocity in the window of interest. Thiscondition does not come naturally from the LDV data acquisition system, since velocitymeasurements are recorded only when a particle crosses the measurement volume. Consequently, the actual velocity-crank angle record features some windows which do nothave any data point in them, and some windows which contain more than one measurement. Simply substituting zeros for missing data points can produce large discontinuitiesin the velocity-crank angle measurements and can lead to chaotic response when such arecord is Fourier transformed.Several algorithms have been proposed for dealing with this situation. One approachis to interpolate (linearly) a missing data from its neighbours, Liou and Santavicca (1983).A sample-and-hold scheme is another method used to assign a data point to an emptywindow at the beginning or end of record. This method entails replacing a missingpoint by the immediately preceding valid measurement, which is retained until the nextvalid measurement is encountered. If more than one measurement is recorded in a givenwindow, the values are arithmetically averaged to obtain a single value for the window.Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 543.7.6 Conditional SamplingThe conditional sampling technique is based on the philosophy that a set of acquired datacan be conditioned, that is, can be subjected to suitable discrimination criteria so thata subset of data can be identified and subsequently processed separately. This techniquewas originally applied to distinguishing between turbulent and irrotational fluids at theouter edges of a turbulent shear flow. From these studies, there result for example:• The experimental determination of intermittency,— the fraction of time the flow isturbulent at a given spatial location;• Zone averages, the average values within either the turbulent fluid or irrotationalfluid; and• The crossing frequency, the number of crossings from turbulent to irrotational fluidper unit time. These findings have added a lot to our understanding of the fluidmechanics of turbulence.The failure so far to identify a definitive procedure for removing the cyclic-variation biasfrom time-phase-average turbulence measurements (if indeed any exists) has led someinvestigators in the field of internal combustion engine to embrace this technique.In order to process the acquired data using the conditional sampling technique, theinvestigator first decides on which parameter(s) to use as a conditioner. In the internalcombustion engine environment one of the most favoured conditioning parameters isthe flame arrival time at an ionization probe. Once the measurement location has beenchosen, the acquired data are then grouped into bins based on arrival time. Each group ofdata bins is interpreted as been subjected to identical conditions. The ensemble averagingtechnique is then applied to each bin, with the hope of obtaining true turbulence intensityand mean velocity.Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 55One of the drawbacks in the conditional sampling technique is the difficulty in selecting the most representative conditioning parameter, since those which are easily measurable, like peak pressure, flame speed and flame arrival time, are time or space-averagedquantities. The simultaneous sampling of a number of these parameters could minimizethe integration effect, Witze et al (1984). A limiting factor in conditional sampling analysis is that only the most probable events are statistically represented. For example, itis not possible to capture early flame development and late end-gas burn-up processes.3.7.7 Autocorrelation FunctionsThere are a number of length scales which characterize different aspects of turbulent flow.The largest eddies in the flow are constrained in dimension by the system boundaries,and they are highly anisotropic. These large-scale nearly two-dimensional eddies arefed by the continual merging of small-scale motions which are nearly two-dimensional incharacter (see Reynolds (1980)).In internal combustion engines, the eddies responsible for most of the turbulenceproduction during the intake process are the large eddies in the conical inlet jet flow.They are characterized by the integral length scale, £. The eddies are roughly equal insize to the local jet thickness and they are approximately of equal size as the valve lift.Furthermore, they become larger in size as they move away from the valve in a wideningjet flow. This situation is illustrated in Figure 3.6.The large-scale turbulence have superimposed on them a range of eddies of smaller andsmaller sizes. These small scale eddies derive their energy from the continual breakdownof larger three-dimensional eddies. Since the smaller eddies respond more quickly to localchanges, they tend to be more isotropic in structure and more homogeneous in nature.On this basis, the small-scale structure of turbulence in an internal combustion engineis assumed to be very similar to those present in other turbulent flows, and this has• Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 56a profound implication for turbulence modelling in ICE. The smallest scale structuresappear to be very thin shear layers distributed throughout the flow domain like a tangleof ribbons, Tennekes (1968). The dissipation of turbulent kinetic energy into heat by theaction of kinematic viscosity takes place in these thin layers.The classic statistical tool for assessing the relative contributions of deterministic andrandom processes is the autocorrelation function (ACF). The velocity autocorrelationfunction expresses the extent to which the velocity fluctuations at one spatial location ortime are coupled to or influenced by those at another position or time, or, loosely, howwell the velocity fluctuations retain their memory.Velocity measurements taken at two points in fluid flow separated by a distance r,significantly less than the integral length scale, L, will correlate with each other; whereasthose made for r far greater than the integral length scale, L, will not correlate. Thespatial autocorrelation coefficient, R(r), of the fluctuating velocity at two adjacent pointsin the flow with respect to the variable distance r between the points is given by;— 1 u(ro)u(ro+r) 316— Nm — 1 u’(ro)u’(ro + r)where Nm is the number of measurements, r0 is the measuring coordinate and r is thespatial separation.The integral length scale, L, is defined in terms of the area under the normalizedcorrelation curve; thus:Li=JR(r)dr (3.17)The integral length scale is determined primarily by the extent of coherent motions andhence is a measure of the size of large eddies. This technique for determining the integralscale requires simultaneous measurements at two points. This is not an easy task toaccomplish in the internal combustion engine, hence most attempts to determine lengthscales have first used correlations to determine the integral time scale, TI.Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 57The integral time scale, TI, of turbulence is defined by the autocorrelation betweentwo velocities at a fixed point in space, but separated in time; thus:ri__f R(t)dtwhere,R(t)= N,n 1 : u(to)(uU( + t) (3.18)Under conditions where the turbulence pattern is convected past the observation pointwithout significant distortion (frozen field hypothesis), and turbulence is homogeneousand relatively weak, (u’/U << 0.2), Taylor (1938), the integral length and time scales arerelated by:Lr = T7r1 (3.19)In engine flows these conditions are not satisfied, so the application of Equation 3.19 inthe engine environment is doubtful.A third scale which is useful for characterizing turbulence is the Taylor microscale,£M. The micro length scale, tM is defined by relating the fluctuating strain rate of theturbulent flow field to the turbulent intensity. It can be determined from the curvature ofthe autocorrelation curve at the origin, Tennekes and Lumley (1972). Usually, the microtime-scale, TM, is determined from the autocorrelation function of Equation 3.18; thus:= [8t] (3.20)at2 t=oFor turbulence which is homogeneous (has no spatial gradients) and is isotropic (has nopreferred orientation), the microscales tM and TM are related by;£M=UTM (3.21)A measure of the size of the smallest turbulence structures is determined by the rateat which turbulent kinetic energy must be dissipated per unit mass, e, and the fluidChapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 58kinematic viscosity, v. This is called the Kolmogorov microscale, 7k. The Kolmogorovrnicroscale is usually obtained from dimensional analysis as:1/477k= (--) (3.22)Another useful timescale is the Kolmogorov time scale, 7k,,1/2Tk=(3.23)This characterizes the momentum-diffusion time of the smallest structures.The final quantity of interest is the turbulent Reynolds number, Re,,, defined as:Re = (3.24)where the kinematic viscosity, v, is usually computed using the Sutherland’s model andA is the Taylor micro length scale. This particular Reynolds number is an importantparameter to the description of the structure of turbulence, and may be interpreted asthe ratio of the turbulence time scale to the molecular time scale. This might also beconsidered as the ratio of the apparent turbulent viscosity to the molecular viscosity. Theturbulent Reynolds number characterizes the degree of separation between the largestand smallest eddies in a turbulent flow. The separation between scales increases withReynolds number, such that at very large Reynolds number, for example Re,, > 4000for stationary turbulent flows, an inertial subrange comes into being between the largerenergy containing eddies and the smaller dissipation range. The inertial subrange isimportant in that the turbulence in this region is statistically independent of both thelarger energy-containing eddies and the smaller viscous-dissipation eddies.For the non-stationary in-cylinder velocity, these ideas are generalized by replacingtime with crank angle and defining the non-stationary velocity autocorrelation function,R(8,T) as, Bendat and Piersol (1986);R8 u(8)u(O + r) 3 25“ “h’— u’(8)u’(O+r)Chapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 59where,u(S) = U(6,i)— (3.26)is the fluctuation about the (time-varying) ensemble-mean velocity, u’(S) is the corresponding ensemble-rms fluctuation intensity, and S is referred to as the reference crankangle and r as the delay or lag.In contrast to Equation 3.18, notice the dependence on S and the normalization by thevalues of u’ at the reference and delay crank angles. These features explicitly reflect thenon-stationary situation. Several methods exist for evaluating Equation 3.25 from theraw LDV data. One such procedure is the discretized lag product or slotted autocorrelationscheme, Roberts and Ajmani (1986), Glover (1986). This technique entails using a seriesof discrete windows or slots of width LT in the delay T to estimate the autocorrelationfunction. When this algorithm is applied to Equation 3.25, we have;R6 — 1N u(S, i)u(O + r, 3 27‘TI— N(S, r) i=1 u’(S)u’(S + T)Contributions to the summation in Equation 3.27 exist only when both velocity estimatesu(9, i) and u(6+r, i) in cycle i are valid. The quantity N(6, T) is the total number of validcross-product pairs for each (6, T) slot over the entire ensemble of N engine cycles. Missingvelocity estimates are usually not interpolated from neighbouring ones, neither are theysupplied artificially. However, missing velocity values affect the statistical uncertainty ofR(S,r), which depends on [N(6,r)]112.Once Equation 3.27 has been successfully evaluated, such parameters as the integrallength scale, and so on can then be obtained, via the Taylor hypothesis, Taylor (1938).From the foregoing, we notice that the relevant conditions (Taylor Hypothesis) for thevalidity of various concepts postulated here are violated in the internal combustion engineenvironment. Consequently, the use of the resulting turbulence parameters for characterizing the in-cylinder turbulent flow should only be applied for qualitative purposes andChapter 3. EXPERIMENTAL AND DATA ANALYSIS METHOD 60results interpreted with extreme caution.3.8 Data Analysis Method of Present WorkAlthough the cycle-resolved and conditional sampling methods claim to have the abilityto reveal the difference, if any, between cyclic variations of the mean flow and low-frequency turbulence, they involve arbitrariness in the choice of cut-off frequency andin the selection of the conditioning parameter(s). We have also seen that the variousturbulence parameters emanating from the autocorrelation function can not be rigorously used to characterize in-cylinder turbulence, since the assumptions on which theirestimates are based are not strictly valid in the ICE environment. In particular, theconditions stipulated in the Taylor hypothesis are not strictly met in ICE combustionchambers. Conversely, ensemble averaging is unambiguous but suppresses phase information about the flow in a particular engine cycle.From a mode1lin point of view, it is very important to choose a data reductionmethod that has a modelling counterpart. Current multidimensional models, however,solve the ensemble-averaged equations of motion and they do not possess the capability toextract any information out of the other methods, (cycle-resolved, conditional samplingand autocorrelation function).In view of the foregoing, vis--vis the objectives of the present work, the ensemble averaged technique was selected for processing the data acquired,— since it has a modellingcounterpart (KIVA-Il).Chapter 4NUMERICAL MODELLINGIn this chapter, a brief description of the KIVA-Il numerical code is presented. Alsoincluded in this chapter are brief summaries of the numerical procedures that are used toaccomplish the solution. Since the objectives of this thesis are in essence the study of fluidmotions in internal combustion engine chambers, spray dynamics and combustion are notexplicitly considered in the expositions given here. A more comprehensive treatment ofall the transactions stated here can be found in Amsden et al (1989).4.1 General BackgroundThe in-cylinder flow fields of advanced internal combustion engines, such as the direct-injection stratified-charge (DISC) engine is characterized by a number of complex, intimately linked physical and chemical processes. These processes consist of the transientthree-dimensional dynamics of evaporating fuel sprays interacting with flowing multi-component gases undergoing mixing, ignition, chemical reactions and heat transfer. TheKIVA-Il numerical code, Amsden et al (1989) is capable of calculating these types ofcomplicated flow dynamics occurring in engine cylinders with arbitrarily shaped pistongeometries, including the effects of turbulence and wall heat transfer.KIVA-Il can be used to solve the unsteady equations of motion of a turbulent, chemically reactive mixture of ideal gases, coupled to the equations for a single-component vaporizing fuel spray. The gas-phase solution procedure is based on a finite volume methodcalled the arbitrary Lagrangian-Eulerian (ALE) method. Spatial differences are formed61Chapter 4. NUMERICAL MODELLING 62on a finite-difference mesh that subdivides the computational domain into a number ofsmall cells that are hexahedrons. The corners of the cells are called vertices, and thepositions of the vertices may be arbitrarily specified functions of time, thereby allowinga Lagrangian, Eulerian or mixed mode description. The arbitrary mesh can conform tocurved boundaries coupled with the capability to follow variations in combustion chamber geometry. An attribute of this method is that the mesh need not be orthogonal. Thedifferencing is made conservative wherever possible. The method employed in KIVA-IT isto difference the fundamental equations in integral form, with the volume of a typical cellused as the control volume, and with divergence terms transformed to surface integralsusing the divergence theorem, Hirt et al (1974).The cartesian components of velocity vectors are stored at cell vertices, and the momentum equations are differenced in a strictly conservative manner. Cell-faced velocitiesare used during a portion of the computational cycle, this greatly reduces the susceptibility of the ALE method to parasitic velocity modes, thereby considerably eliminatingthe need for a node coupler.4.1.1 General Solution ProcedureThe unsteady solution is marched out in a sequence of finite time increments called cyclesor timesteps. During each cycle the values of the dependent variables are calculated fromthose of the previous cycle. Each cycle is split into two phases, namely; a Lagrangianphase and a rezone phase. In the Lagrangian phase the vertices move with the fluidvelocity, and there is no convection across cell boundaries. In the rezone phase, theflow is frozen, the vertices are moved to new user-specified positions, and the flow-fieldis remapped or rezoned into the new computational domain. The rezoning is accomplished by convecting material across the boundaries of the computational cells, whichare regarded as moving relative to the flow field.Chapter 4. NUMERICAL MODELLING 63The temporal difference scheme in KIVA-Il is largely implicit. As a result of this, thetimesteps used by KIVA-Il are calculated based on accuracy and not stability criteria.Since larger timesteps can be employed without the risk of computations going unstable,a considerable savings in computational time is therefore possible. In the Lagrangianphase, implicit differencing is used for all the diffusion terms and the terms associatedwith pressure wave propagation. The coupled implicit equations are solved by a methodsimilar to the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE), Patankar(1980).Explicit methods are used to calculate convection in the rezone phase; but the convection calculation can be subcycled an arbitrary number of times, and thus the maincomputational timestep is not restricted by the Courant stability condition of explicitmethods. The convection timestep is a submultiple of the main computational timestepand does satisfy the Courant condition. The KIVA-Il code has various schemes for differencing the convection term. These include partial donor cell, interpolated donor cell andQuasi-Second-Order Upwind (QSOU) differencing schemes. The QSOU is monotoneand approaches second-order accuracy when convecting smooth profiles. Quasi-second-order-upwind scheme is more accurate than the partial donor cell differencing method,but is time-consuming. Therefore, in situations where quick and rough estimates arerequired, it is preferable to use the partial donor cell differencing scheme. Interpolated-donor cell differencing was used in the present study.4.1.2 Turbulence ModelsThe effects of turbulence are represented by two types of models in KIVA-Il. They are thestandard k- turbulence model, Launder and Spalding (1972), modified to include volumetric expansion effects and spray/turbulence interactions and the subgrid scale (SGS)Chapter 4. NUMERICAL MODELLING 64turbulence model. In the subgrid scale model, a transport equation for the turbulent kinetic energy associated with subgrid scale motions are solved. The SGS model however,reduces to the k-c model near walls where all turbulence length scales are too small tobe resolved by the computational mesh. Boundary layer drag and wall heat transfer arecalculated by matching to a modified turbulent law of the wall.4.1.3 Mesh Generation SchemeBecause of improvements to the code’s ease-of-use and versatility for many applications,all required geometrical specifications, initial conditions and boundary conditions can bespecified by using the standard input alone. This is especially true for internal combustionengine applications.The mesh generation logic allows the computational domain to include bowl-in-pistonsand domed cylinder heads and to offset these relative to the axis of the cylinder. In addition to two- and three-dimensional cartesian and cylindrical meshes, the code allows thecalculation of the flow in a single sector of certain three-dimensional cylindrical configurations in which there is an n-fold symmetry about the axis of the cylinder. This symmetryis often found in engine cylinders with multihole injectors. For initial conditions, one canspecify an axisymmetric swirl-velocity field with a Bessel function profile and a specifiedswirl ratio. Standard boundary conditions and rezone logic allow the mesh to follow themotion of a piston.4.2 The Governing Equations• The equations of motion for the fluid phase and the boundary conditions are given inthis section. For compactness these equations are written in vector notations. The unitvectors in the , y and z directions are denoted by 3 and k respectively. In order toChapter 4. NUMERICAL MODELLING 65preserve all the vital features of KIVA-Il code, the equations stated here follow those inAmsden et al (1989) very closely. However, a detailed derivation of the k-e turbulencemodel is presented in appendix A.The position vector i is defined as:= xi+y+zkthe vector operator V is given by;.a a -aV =z— +j— + k—ox Oy Ozand the fluid velocity vector iiis given by;it = u(x,y, z, t)+ v(x,y, z, t)+ w(x,y, z,where t is time.4.2.1 The Fluid PhaseThe KIVA-Il equations can be used to solve for both laminar and turbulent flows. Themass, momentum and energy equations for the two cases differ primarily in the formand magnitude of the transport coefficients (that is, viscosity, thermal conductivity andspecies diffusivity), which are much larger in the turbulent case because of the additionaltransport caused by turbulent fluctuations. In the turbulent regime, the transport coefficients are derived from a turbulent diffusivity that depends on the turbulent kineticenergy, k, and its dissipation rate, .Continuity EquationThe continuity equation for species m is:8Pm+ V • (pmii) V [pDV(1)] +, + m1 (4.1)Chapter 4. NUMERICAL MODELLING 66where Pm 15 the mass density of species m, p, the total mass density and ii, the fluidvelocity. Fick’s Law is assumed for diffusion processes, with a single diffusion coefficient,D. The equation for D will be given later. Species 1 is that for which the spray dropletsare composed, and 5ij is the Dirac delta function. By summing Equation 4.1 over allspecies, we obtain the total fluid density equation, thus;(4.2)since mass is conserved in chemical reactions. Where ‘ is the spray density. For coldflow simulations, which is the subject addressed in this thesis; = 0.Momentum EquationThe momentum equation for the fluid mixture is;ô(pu)+V.(piii)=—-V_AoV(pk) +v.&+i+pj (4.3)where p is the fluid pressure. The dimensionless quantity a is used in conjunctionwith the Pressure Gradient Scaling (PGS) Method (see Amsden et al (1989)). This isa method for enhancing computational efficiency in low Mach number, Ma, flows wherethe pressure is nearly uniform.In Equation 4.3 the quantity A0 is zero in laminar calculations and unity when oneof the turbulence models is invoked.The viscous stress tensor, ö, is Newtonian in form, therefore;= LLff [v’i + (Vi)tj + v • (4.4)The first and second coefficients of viscosity, /1eff and ) are defined later. The superscript t denotes the transpose and I is the unit dyadic. F5 is the rate of momentum gainper unit volume due to the spray. The specific body force, , is assumed constant.Chapter 4. NUMERICAL MODELLING 67Energy EquationThe energy equation is;ôCoI)+ v • (pill) = —pV • ii + (1 — A0) : Vil — V • 1+ Aop€ + + C (4.5)where I is the specific internal energy, exclusive of chemical energy. The heat fluxvector I is the sum of contributions due to heat conduction and enthalpy diffusion;1= _KVT_PDhmV(1) (4.6)where T is the fluid temperature, hm, the specific enthalpy of species, m, Qc and Q8are the source terms due to chemical heat release and spray interactions respectively. Inthe present study, Qc = = 0.Turbulent Kinetic Energy and Dissipation ModelsWhen one of the turbulence models is in use (A0 = 1), two additional transportequations are solved for the turbulent kinetic energy, k and its dissipation rate, e, thus;+ V • (pilk) = —pkV . ii + o• : il + . [(tt.) Vkj — p6 + W (4.7)and;+ v.(p = — (c — p6V • il + . [(v’) e] +.[co. : Vil — c2pe + cWi (4.8)These are standard k-6 equations, Launder and Spalding (1972) with some addedterms. The source term (CE3 — )V • ii in the 6-equation accounts for length scalechanges when there is velocity dilation. Source terms involving the quantity W8 arisedue to interaction with the spray. In the present study (W’ = 0).Chapter 4. NUMERICAL MODELLING 68The quantities c61, c, C63, Prk and Pr6 are constants whose values are determined• from experiments and some theoretical considerations. Standard values of these constantswere used in the present study, and are given in Table 4.1.The physical meaning of the various terms in Equation 4.7 is as follows. The termV • (piZk) simply represents convection of the turbulence by the resolved velocity field.The term —(2/3)pkV • il is a compressibilty term that represents the turbulent analogof pdV work. The term c3- Vii represents the production of turbulence by shear inthe resolved velocity field, and V •[(6ff/Prk)Vkj represents the self-diffusion of theturbulence with a diffusivity of efi/P. The term —pe represents the decay of turbulenceinto thermal energy.When the SGS turbulence model is used, the value of e is constrained to satisfy theinequality:,1C61)] LsGs(4.9)LSGS is an input SGS length scale whose value is typically taken to be 4Sx, where ax isa representative computational cell dimension. Inequality 4.9 is enforced by integratingEquations 4.7 and 4.8 in time at all points and then setting € equal to the right-hand sideof Equation 4.9 at points where the inequality is violated. Since k3/2€ is proportional tothe ic-c length scale, Equation 4.9 is a constraint that the turbulent length scale be lessthan or equal to LSGS. Near rigid walls this constraint is always satisfied and thus thestandard k-c equations are solved near walls. In regions where the length scale is LSGS,the model reduces to a one-equation model.State RelationsThe state relations are assumed to be those of an ideal gas mixture. Therefore,(4.10)Chapter 4. NUMERICAL MODELLING 691(T) = () Im(T) (4.11)c(T) = > () cpm(T) (4.12)and,hm(T) = Im(T) + (4.13)where .l? is the universal gas constant; Wm, the molecular weight of species m; Im(T),the specific internal energy of species m; and c, the specific heat at constant pressureof species m. The values of hm(T) and cpm(T) are taken from the JANAF tables, Stulland Prophet (1974).The transport coefficients in KIVA-TI are taken to be,k2ILefI = (1.0 — Ao)pvo + ,Uair +A0c,—=K — ILeffCppScand,D = (4.14)pScThe diffusivity v0 is an input constant, and c is an empirical constant with a standard• value of 0.09. A Sutherland formula is used for iajr:A1T312= T + A2(4.15)where A1 and A2 are constants. The constant A3 is taken to be— in calculationsof turbulent flow, but can be arbitrarily specified in laminar flows. The Prandtl andSchmidt numbers, Pr and Sc, are input constants.Chapter 4. NUMEPJCAL MODELLING 704.3 The Computational DomainThe computational domain over which solutions of the differential equations describedin section 4.2 are required is that lying between the cylinder surfaces, cylinder head andthe piston face.In computational fluid dynamics parlance, these surfaces are referred to as walls. Inreciprocating engines two of these walls move— the piston face moves continuously intime and the valves move intermittently. If a finite-difference computational mesh is usedwhich is fixed in space, much of the grid will not be utilized during parts of the pistoncycle, and in order to obtain accurate results when the piston is near top dead centrea very fine mesh is required near the cylinder head. A mesh which does not translatewith the piston motion is therefore very wasteful of both computer time and storage.This problem may be overcomed by defining a mesh which moves with the motion of thepiston such that the whole mesh is always confined between the cylinder and piston heads.Various ways of accomplishing this task have been proposed. For example, the governingequations can be transformed into a Lagrangia.n form thereby allowing the grid nodes totranslate with the fluid (see Hirt et al (1970)); however, for multi-dimensional flows thismethod suffers from inaccuracies due to mesh distortion. This drawback can be removedby means of a hybrid Eulerian-Lagrangian scheme, such as the Implicit Continuous-FluidEulerian (ICE)-ALE method of Hirt et al (1974). This method uses a finite differencemesh with vertices that may move with the fluid (Lagrangian), be held fixed (Eulerian),or be moved in any other prescribed way. Because of this flexibility the method is referredto as an Arbitrary-Lagrangian-Eulerian (ALE) technique. This attribute is particularlyuseful for those flows in which the movement of the boundary cannot be predicted, as forexample, in the calculation of the position of the interface between two different liquidswhich are subjected to external forces.Chapter 4. NUMERICAL MODELLING 71The advantages of the ICED-ALE method include its ability to resolve arbitrary confining boundaries, to have variable zoning for purposes of obtaining optimum resolution,to be almost Lagrangian for improved accuracy in problems where fully Lagrangian calculations are not possible and to operate with time steps many times larger than possiblewith explicit methods. This is the approach adopted in the numerical scheme (KIVA-Ilcode) used in the present study. Another method of achieving mesh translation withthe moving wall is by co-ordinate transformation, Gosman (1985b). The instantaneousposition of the moving boundaries (piston and valves) are often used for scaling.4.4 Boundary ConditionsIn order to completely specify a mathematical problem it is necessary to supply theconditions at the boundaries of the solution domain for the variables. These conditionsare usually obtained by either specifying the value of the dependent variable at theboundary, or the value of the associated flux, or a relation between the two.The specifications of the relevant boundary conditions in KIVA-Il code also followthese basic rules.4.4.1 Walls4.4.2 Velocities and PressureVelocity boundary conditions on rigid walls are introduced either by imposing the valueof the velocity on walls or the value of the wall stress ö, = • , where is the unitnormal to the wall. On no slip walls, the gas velocity is set to the wall velocity;UWwallk (4.16)where the wall is assumed to be moving with speed in the z-direction. The wallChapter 4. NUMERICAL MODELLING 72stress is then determined implicitly through Equation 4.4. On free-slip and turbulent law-of-the-wall boundaries the normal gas velocity is set equal to the normal wall velocity,U•fl=Wwallk•fl (4.17)and the two tangential components of ow are explicitly specified. For free-slip walls thetangential components of are zero.Turbulent Law-of-the-Wall ConditionsFor turbulent law-of-the-wall conditions the tangential components are determinedby matching to a logarithmic profile:-- = + B for > R= .1/2 for < R (4.18)where= IpV (4.19)fLatr( )is the Reynolds number based on the gas velocity relative to the wall,=—Wwajlk (4.20)which is evaluated a distance, y, from the wall, and u’ is the shear speed, which is relatedto the tangential components of the wall stress by:(4.21)where,V = U— WwajllCIn Equations 4.19 and 4.21 it is assumed that y is small enough to be in the logarithmicregion or the laminar sublayer region of the turbulent boundary layer. The ReynoldsChapter 4. NUMERICAL MODELLING 73number, R defines the boundary between these two regions. The constants IG, c, R,‘and B in Equation 4.18 are related to the k-€ model constants by:— /1/2,Ic— v c — c61 4. 2and,B = R’12 — (4.23)For commonly accepted values of the k-e constants, B = 5.5 and cj = —0.15, we obtain,Ic = 0.4327 and R =114.When the velocities are specified at a boundary, the pressure need not be specifiedthere.4.4.3 Stagnation EnthalpyIn principle the boundary conditions for the enthalpy equation are simple to supply ifeither the inner wall temperature or the heat flux through the inner wall surface is fixed.In practical engines neither of these conditions apply; the heat flux and temperatureboth vary in space and time due to compression, combustion and cooling. Consequently,• recourse to experimental data may be necessary, or advantage may be taken of the factthat the outer wall temperatures are more precisely known than the inner wall valuesbecause they are in contact with the coolants; water for the cylinder head and wall, andoil for the piston head. The latter fact was made use of by Watkins (1977) in specifyingtemperature boundary conditions, via the solution to the one-dimensional transient heatconduction equation. However, this approach was not used in the KIVA-Il numericalcode. The procedure employed in KIVA-Il is as follows:Temperature boundary conditions on rigid walls are introduced by specifying eitherthe wall temperatures or the wall heat flux J = —KVT • . For adiabatic walls, isset to zero. For fixed temperature walls that are also either free slip or no slip, the wallChapter 4. NUMERICAL MODELLING 74temperature is prescribed, and L is determined implicitly from Equation 4.6. For fixedtemperature walls using the turbulent law-of-the-wall condition, J,, is determined fromthe modified Reynolds analogy formula:=1forCRPrç1= for c> .R (4.24)Pr[+Q?j_i)bY2jwhere T is the wall temperature and Prt is the Prandtl number of the laminar fluid.In addition to the wall heat loss, there is a source to the internal energy due tofrictional heating. Frictional heating occurs whenever turbulent law-of-the-wall velocityconditions are used and has the form:flU = •= p(u)2v (4.25)where, f is the heating rate per unit area of wall.4.4.4 Turbulent Kinetic Energy and Dissipation RateIn calculations of turbulent flow, boundary conditions are also needed for the turbulentkinetic energy, k and its dissipation rate, . The standard k-e model is not valid in thebuffer or the linear sublayer of the turbulent boundary layer because the local Reynoldsnumberk2/pe is low, in this regime. Consequently, boundary conditions for Equations4.7 and 4.8 as well as the ensemble-averaged conservation equations of k and e are appliedat a distance away from the wall in the fully turbulent region (for example, the log-lawlayer) where the local Reynolds number is sufficiently high. The boundary conditionsapplied at such a location away from the wall are referred to as wall functions.In the KIVA-Il coLle, these wall functions are given below.Since the shear stress near the wall is nearly constant, the wall function for turbulentChapter 4. NUMERICAL MODELLING 75kinetic energy, k is given by:Vk.=0 (4.26)and, in the log-law layer, the production of turbulence is approximately equal to its rateof dissipation (local isotropy). Therefore, the wall function for e is approximated by;k3126 = (4.27)lipwhere k and e are evaluated a distance y, from the wall, and,1/2= [T:_ CEJ]4.5 The Numerical SchemeKIVA-JI solves finite-difference approximations to the governing equations of section 4.2.The equations are discretized both in space and time. In order to fully understand thenumerical scheme, a brief discussion of some of its general features are presented. Fora more detailed description of the issues presented here, reference should be made toAmsden et cii (1989).4.5.1 Temporal DifferencingThe temporal differencing is performed with respect to a sequence of discrete times t’(n = 0,1,2,...). The time interval /.t” t’1 — t’ is the timestep, and the integeris the cycle number. The latter is displayed as a superscript, so that QTh denotes thedifference approximation to the quantity Q at time t’. When iXt appears without asuperscript, tt’ is understood. The difference approximation to the derivative ÔQ/Otis the first-order expression (Q’— Qj/Lt. In the KIVA-Il code, a cycle is performedin three stages or phase. Phases A and B together constitute a Lagrangian calculationin which computational cells move with the fluid. Phase A is a calculation of sprayChapter 4. NUMERICAL MODELLING 76droplet collision and oscillation/breakup terms and mass and energy source terms due tochemistry and spray. Phase B calculates in a coupled, implicit fashion the acoustic modeterms (namely the pressure gradient in the momentum equation and velocity dilatationterms in mass and energy equations), the spray momentum source term and the termsdue to diffusion of mass, momentum and energy. Phase B also calculates the remainingsource terms in the turbulence equations.In phase C, the flow field is frozen and rezoned or remapped onto a new computational mesh. Intermediate quantities that have been partially but not fully updated areidentified by superscripts A and B. For example, Q’ is the computed value of Q at theend of phase A. Notice however, that superscript C is not needed because it is equivalentto index n + Spatial DifferencingThe spatial differencing is based on the ALE method, which in three dimensions usesa mesh made up of arbitrary hexahedrons. Spatial difference approximations are constructed by the control-volume or integral-balance approach, which largely preserves thelocal conservation properties of the differential equations.The computational domain of interest is subdivided into a number of small cells orzones, the corners of which are the vertices. The total number of such cells make-upthe mesh with respect to which spatial differences are formed. The vertices need not bestationary, but may translate in any arbitrarily prescribed manner. The Lagrangian andthe Eulerian modes are special cases. In general, the cells are asymmetrical; a typical cellis shown in Figure 4.1, with the vertices conventionally numbered. The cells are indexedby integers (i,j, k) which may be regarded as coordinates in logical space. The indices(i, j, k) also label the vertices, with the understanding that vertex (i, j k) is vertex 4for cell (i,j, k). The cartesian coordinates of vertex (i,j, k) are (x13k, Ysjk, zI,k), which inChapter 4. NUMERICAL MODELLING 77general depend on time t. Thus the position vector of vertex (i, j, k) is;ik = Xi,kZ + YijkJ + Zkk (4.28)It is convenient to define auxiliary cells centered about the vertices. These cells are calledmomentum ceiTh since their primary function is in differencing the momentum equations.Momentum cell (i, j, k) is centered about vertex (i, , k). In contrast to regular cells, whichhave six faces, momentum cells have twenty-four faces, each of which is comparable in sizeto one-fourth of a regular cell face. Three of these twenty-four faces lie within each of theeight regular cells which share common volume with the momentum cells. The portion ofthe momentum cell (i, j, k) lying within regular cell (i, j, k) is shown in Figure 4.2. Thepoints of intersection of the momentum cell faces with the regular cell edges are definedas the midpoints of the regular cell edges. The points of intersection of the momentumcell edges with the regular cell faces are then defined implicitly by the requirement thatthe regular cell face be partitioned into four subfaces of equal area by the momentum cellfaces. The corners of the momentum cells are then implicitly defined by the requirementthat the overlap volume between a regular cell and a momentum cell centered at one of itscorners be one-eighth of the regular cell volume. Generally, the momentum cell corners donot coincide with the cell centers. The momentum cell corners and the intersection pointsof momentum cell edges with regular cell faces are not actually solved since they are notneeded. The location of velocities at cell vertices in the ALE method is convenient becauseno interpolation is required when determining vertex motion in the Lagra.ngian phase ofthe calculation, but it has a major drawback. The disadvantage is that the ALE methodsolutions are notoriously susceptible to parasitic modes in the velocity field. The mainreason for this is because pressure waves tend to propagate along cell diagonals insteadof via adjacent cells (see Amsden et a! (1989)). A checker boarding effect is therebycreated in the pressure fields, with associated irregularities in the velocity field that areChapter 4. NUMERICAL MODELLING 78usually suppressed by the introduction a numerical damping called node coupling, Hirtet al (1974). In a major improvement to the ALE method, Amsden et al (1985b), thesusceptibility to parasitic modes were alleviated by the introduction of velocities centeredon all faces. Vertex velocities are retained, and momentum associated with the vertices isconserved, but normal velocity components on cell faces are used to compute cell volumechanges in Phase B and fluxing volume in Phase C. The resulting scheme greatly reducesthe need for node coupling.Accelerations of cell-face velocities due to pressure gradients are calculated by constructing momentum control volumes centered about the cell faces. Like the momentumcells, the cell-face control volumes have twenty-four faces. Referring to Figures 4.1 and4.2, the cell-face control volume for the left face of cell (i, , k) is composed of those portions of the momentum cells of vertices, 3, 4, 7 and 8 that lie in regular cells (i, j k) and(i— 1,j, k). Faces associated with cell (i, j k) are shown in Figure 4.3. Control volumesassociated with the other faces are defined analogously.hi the finite-difference approximations of KIVA-Il, velocities are fundamentally located at the vertices, so that:U23k = u(xk,ytjk,zjjk) (4.29)Thermodynamic quantities are located at cell centers, such that:(4.30)where Q = p, p, T, I or Pm, as well as k and . Quantities needed at points where theyare not fundamentally located are obtained by averaging neighbouring values.Spatial differences are usually performed by integrating the differential term of interestover the volume of a typical cell (or momentum cell). Volume integrals of gradient ordivergence terms are usually converted into surface area integrals using the divergenceChapter 4. NUMERICAL MODELLING 79theorem. The volume integral of a time derivative may be related to the derivative ofthe integral by means of Reynolds’ transport theorem, Karamcheti (1980). Volume andsurface area integrals, are usually performed under the assumption that the integrandsare uniform within cells or on cell faces. Therefore area integrals over surfaces of cellsbecome sums over cell faces (or subfaces):a(4.31)When differencing diffusion terms for cell-centered quantity, Q, it is necessary to evaluate(VQ) • La. Referring to Figure 4.4, this quantity is approximated as follows. The pointsand i are the centers of the cells on either side of face a and i, i,, ir and i1 arethe centers of the four edges bounding face a. We first solve for coefficients a, at6, andaid such that:aj,.( — i.) + at6(i— 4) + ajd(i, — ij) = La (4.32)Note that since the mesh may be non-orthogonal, the vector (aij—) need not be parallelto La and therefore atb and aid may be nonzero. The finite-difference approximation to(VQ) • La is obtained by dotting both sides of Equation 4.32 with (VQ)0 and ignoringterms of second and higher order in the cell dimensions:aj,.(Qj’— Qr) +a6(Q — Qb) + af(Qf — Qd) = (VQ)a • La (4.33)In Equation 4.33, Qt is the simple average of the values of Q in the four cells surroundingcell edge t and Qb, Qj and Qd are defined analogously.Area integrals over momentum cell faces are ordinarily converted into integrals overregular cell faces by the following procedure. Let Q be a quantity that is uniform withinregular cells, and consider the volume overlap between regular cell (i,j, k) and the momentum cell associated with one of its vertices. Three faces of this overlap volume (callthem a, b, c) are faces of the momentum cell in question, with outward area vectorsChapter 4. NUMERICAL MODELLING 80• while the other three (call them d, e, f) are surfaces of regular cell (i, , k), with outwardarea vectors Aa. But the divergence theorem shows that the integral f dA over theentire surface of this overlap volume is zero, so that;—-+ 1 — -. —A + A, + A = —(Ad + Ae + A1) (4.34)Thus the integral f QdA over the three momentum cell faces in question may be represented by;J QdA = Qk( + + = Qjk(Xd + Ae + (4.35)so that the area vectors A never need to be explicitly evaluated. A similar procedure isused to express the outward normal areas A’ of faces of the cell-face control volumes interms of the regular cell face areas A’a.The mass of cell (i,j, k) is denoted by M13k and is given by;M13k = (4.36)The mass of momentum cell (i, k) is given by;M = + M_1,,k + M_1,_1k + M,_1,k -I- + M_1,J,k_1 ++ .M1,1,k—1) (4.37)The mass of left cell-face control volume of cell (i, j Ic) is given by;M7 = (M,k + M1_l,,,k) (4.38)and the mass of other cell-face control volumes are defined analogously.4.5.3 Lagrangian Phase Difference EquationsWith the background given above, the KIVA-Il difference equations can now be specified. For the purpose of brevity, the Lagrangian phase difference equations are specified,Chapter 4. NUMERICAL MODELLING 81comprehensive details concerning the Phase C or the rezone phase are contained in Amsden et a! (1989). The implicit methods are used to difference the terms associated withacoustic pressure wave propagation and diffusion of mass, momentum and energy.4.5.4 Mass Density EquationsThe Lagrangian phase difference approximation to Equation 4.1 is:(Pm)P,kV— (Pm)ikk= (pD’V [6DYmB + (1— cD)] • +[C1rn)tik + PijicSmi] Vik (4.39)The mass fractions are related to the densities by:(4.40)where x = n, A or B. The phase A mass densities are defined below. An importantfeature of Equation 4.39 is the use of variable implicitness parameter 4’D in differencingthe diffusion term. Parameter çb. varies in space and time and is defined at cell centers.Its value which lies between zero and one, depends on the local Courant number, Cd;1LLtCd p tx2where t2x is a measure of the cell size. When Cd is small compared to unity, i’D is zeroand a fully explicit difference approximation is used. When Cd is large compare to unity,i’D is close to unity and an implicit formulation is used.The Phase A density of species rn includes contributions due to chemistry and sprayevaporation.( ‘A ( ‘nIPmhjjc — .Pm.:,ic •clxt = (Pm)ijk + Ptjk8ml (4.41)where ()ijk and 13jk are the chemistry and spray mass source terms respectively.Chapter 4. NUMERICAL MODELLING 82Summing Equation 4.39 over all species m gives the following Lagrangian phase difference approximation to the mass density Equation 4.2:BiyB nPk ‘:jk —• T/n—Pi,kv,kThis shows that the total gas mass in a cell changes only due to the spray source.Similarly, summing Equation 4.41 over all species and comparing with Equation 4.42shows that Phase B cell masses are known after Phase A:B 17B — A un — AA —Pji ijk — Pi3k ‘ijk ijk— ijkEquation 4.43 can be combined with Equation 4.39 to 4.41 to obtain;Mk[(Ym)k_(Ym)tk]= E(pD):V [DY + (1 — 4’D)Y] . A (4.44)Equation 4.44 is solved in Phase B.4.5.5 Momentum EquationThe Lagrangian phase difference approximation to the momentum Equation 4.3 is thefollowing:= (n)2 > + (1 — bP)Pj (A’) + F+G + H (4.45)where,,rv-’ 2 ALa! Ainr — —ii“i3V1 )G = + [Dii + (1 — D)0(i)] . (A’)H = —1(Jc,k + S’ijkU,k) + g(M’), — (M’)jk(ANC):k/1tThe index 3 refers to the faces of momentum control volume (i, k), whose normalarea vectors at time tT’ are (A’). The pressure F, turbulent kinetic energy, k, and viscousChapter 4. NUMERICAL MODELLING 83stress tensor, ö, are regarded as being uniform within the regular cell in which face j3lies. Thus the A’ can be eliminated in terms of regular cell face areas A as describedabove (see Equation 4.34). This makes it convenient to evaluate momentum changes dueto surface stresses by sweeping over cells rather than vertices in the computer program.In differencing the pressure gradient term, variable implicitness parameter is used.Parameter 4, plays a role for the acoustic mode terms analogous to that of I’D for thediffusion terms. The coefficient I’, depends on the sound speed Courant number, C3,thus;C3 = (4.46)where c is the isentropic speed of sound.The Pressure Gradient Scaling (FGS) parameter a” which is used to artificially raisethe Mach number in far subsonic flows, depends on time but not on space.The viscous stress tensor, 6 is a weighted average of the stresses based on timelevel n velocities and those base on the Phase B velocities. This weighting employs thesame variable implicitness parameter, I’D that is used for the mass diffusion terms. Thequantities R’13k and S’13k are associated with the implicit coupling of the computed gasand drop velocities. The effects of the alternate node coupler are represented by the term(ANC)13k.4.5.6 Internal Energy EquationThe Lagrangian phase finite-difference approximation to the internal energy equation isthe following:M8 18— M” I” F” -- P8 V8 —i3k :jk ijk ijk :3k :3k :,k :3k- 2 It(1 — A0) [I’D&(jtB): ViiB + (1 — I’D)c7(i”) : ViZI•k JkChapter 4. NUMERICAL MODELLING 84+KV [DTB— 1DTJ .A+ (pD) [ hm(T:)v [DY,I + (1 — c6D)Yi’] ] .+AoMkEk + Vk(Qk + Qk) (4.47)4.5.7 Turbulence Kinetic EnergyThe Lagrangian phase difference approximation to the turbulent kinetic energy, Equation4.7 is:M8 kB— M k TJB — un3k ijk ijk t3k —“ B ‘tjk ijk Fri r \in iB 1 /TT r\—— Ji3k)r1.k Jiikl’.’ijkj .. V L’)i3kyzF [ç +(1 — D)kAj .A.+ (W8)jkVk (4.48)tjk4.5.8 Dissipation RateThe difference approximation to the dissipation rate; Equation 4.8 is;MB€—M’ ‘.‘13k :3k ijkCi,k— ( ‘ \ B uk ijk Fr i n : B——— Cea) Pu,kzt 1 — Jijk)Sjk Jtjk6:j+ > v [€ + (1— çbD)e4] • +—cE2Mk---ck + cs(W3)ujkVj7kj (4.49)tjk ijkThe quantity ft,k is zero or unity depending on the sign of the cell volume change;VB_Vijk tjkfu,k = 1, if jk—T’k > 0= 0 otherwise (4.50)The prescription for f1jk is chosen to avoid negative computed values of k and e whenthere are large volume changes during the Lagrangian phase.Chapter 4. NUMERICAL MODELLING 85Viscous dissipation of mean flow kinetic energy is represented by the term (VD)13k;which is given by:(VD)13k Vjk [cbD(u1) : Vu8 + (1 — kD)J(U) : ViiI..k (4.51)4.5.9 Equation of StateThe equations of state 4.10 and 4.11 are approximated by:‘ijk= [ (Prn)jk] R07’iik (4.52)and‘ijlc ‘uk + (cv)Jk(Tk— jk) (4.53)where,Tt — ‘(,r B r (‘pt1ijk —mand=m)tB,kvm(Tijk) (4.54)4.6 Solution Procedure for Implicit Phase B EquationsThe Phase B values of the flow field variables are found by solving the implicit differenceequations stated above. The method of solution is patterned after the SIMPLE procedurewith individual equations solved using the conjugate residual method, (see Amsden et al(1989)). The SIMPLE method is basically a two-step iteration scheme. After selecting apredicted value of the Phase B pressure pB, in step 1, the predicted pressure field is frozenand other flow quantities are re-solved using finite difference equations that differencethe diffusion terms implicitly. In the original SIMPLE method, the convection terms areChapter 4. NUMERICAL MODELLING 86also difFerenced implicitly, and their effects are included in step 1. In KIVA-Il, convection is calculated in Phase C in a subcycled implicit fashion that offers some significantadvantages over implicit methods. In step 2, the values of the diffusion terms obtainedin step 1 are frozen and the corrected pressure field are solved for using equations thatdifference pressure terms implicitly. Sometimes a Poisson equation for the pressure isderived and solved in step 2. In KIVA-Il, for step 2, the cell-face velocity equations, thevolume change equations and a linearized form of the equation of state are solved simultaneously. By algebraically eliminating the volumes and cell-face velocities from theseequations in favour of the pressures, it can be shown that the Poisson equation for thepressure in step 2 is being solved. Following step 2, the predicted and corrected pressuresare compared. If they agree to within a specified convergence tolerance, the equationshave been solved and computation proceeds to Phase C. If the difference between thepressure fields exceeds the convergence tolerance, the corrected pressure field becomesthe new predicted pressure field, and the calculations reverts back to step 1 and the wholeprocess is repeated. Each pass through the two steps is called an outer iteration. Theonly equations in the outer iteration are the momentum, internal energy and pressureequations.4.7 Accuracy Conditions and Automatic Timestep ControlThe timesteps Lt and zt are selected automatically at the beginning of each cycle.Because diffusion terms are differenced implicitly and convective terms are subcycled,there are no stability restrictions on t, but there are several accuracy conditions uponwhich the automatic selection of zt is based. The details can be found in Amsden et al(1989). Brief description of this is presented below. The convection timestep mustsatisfy the Courant condition for stability, and a description of how this is generalized toChapter 4. NUMERICAL MODELLING 87an arbitrary mesh can also be found in Amsden et al (1989), however, a brief descriptionof this is given in this section.The accuracy conditions used to determine L± cannot give a universally reliable selection because there are many accuracy conditions which were not taken into account inKIVA-Il code. It was found from this study, that the criteria we use for determining ztgive temporally accurate solutions in most of the calculations performed. It is however,necessary to mention that other users of the KIVA-Il code may find some other accuracyconditions relevant to their own applications. It is therefore very important to vary thetimestep in order to test for temporary accuracy.The first accuracy condition on Lt is that:(4.55)where fa is some positive real number of order unity and I.S.x is an average cell dimension.This condition arises because terms of order higher than Lit are ignored in the phase Bvertex positions, which are given as:B — n -BXjjk— Xi3k -r U1jkThe accuracy constraint equation 4.55 is the only one we used in which the cell size Lixappears. We note that Lit Lix”2 for condition 4.55, in contrast to explicit convectivestability criteria, which give Lit - Lix, and explicit diffusional stability criteria, whichgive Lit Liz2. Constraint (equation 4.55) is implemented by calculating a timestepLitacc, thus:A n+i — faLiXi,j,k——B .-1,3, j— U2jkwhere= [(Iii — x12 + jz — al2 + jx3 — (4.58)+ 1x2— £i12 + Ixs — iI2 + Is —x41)/6]”Chapter 4. NUMERICAL MODELLING 88and the subscripts in equation 4.59 refer to the vertices of cell (i,j, k) as numbered inFigure 4.1. We then constrained t”’ to be less than as described below. Thedefault value of f4 is 0.5.The second accuracy condition on Lt is that:Pit<fr (4.59)where f7 is of order unity and ). is an eigenvalue of the rate of strain tensor. Thiscriterion limits the amount of cell distortion that can occur due to mesh movementin the Lagrangian phase. When cells become very distorted, the spatial accuracy ofthe difference approximations suffers. One example of how equation 4.59 works is thefollowing. In a plane parallel shear flow, there is one non-zero eigenvalue (1/2 öu/ôy),where u is the streamwise velocity component and y is the cross-stream direction. Asdepicted in Figure 4.5 for a rectangular cell alligned with the flow, constraint (equation4.59) limits the distance upper and lower cell vertices can move relative to one another,divided by the cell height. Constraint (equation 4.59) is enforced by calculating a timestep= mm______(4.60)‘ (i,j,k) 2.Ja1k/3where= (Pk 3q,k) (4.61)and Pijk and qjjk are given in terms of the rate of strain tensor 3em in cell (i, i,Pt;k = (4.62)q,k =61tmn82t3 +6L2m8113 + €Lm331tS2m (4.63)In equation 4.63 ELmn is the alternating tensor, and 8m is given by:1 (Out OUm’\.9tm—+ J (4.64)2 \t9Xm OXLJ..kChapter 4. NUMERICAL MODELLING 89where the velocity derivatives are evaluated using time level n velocities. It can beverified from the formula for the roots of a cubic polynomial that the denominator ofequation 4.60 is greater than the magnitudes of all eigenvalues of 8Lm, and therefore if weselect tt the constraint 4.59 will be satisfied. Using f,. = 1// gave sufficientaccuracy in the calculations reported in this thesis.Two other accuracy criteria for t are obtained from the need to couple accurately theflow field and source terms due to chemical heat release and mass and energy exchangewith spray. For the chemical heat release, the requirement is that:v.nc________A fIf T LX<JchIVLjk1ijkwhere fth is an input constant typically taken to be 0.1. Equation 4.65 is the requirementthat the total heat release from all chemical reactions in a cell should not exceed a smallfraction fth of the total internal energy in the cell. To enforce equation 4.65 we calculatetimestep Ltth as:=[i] } (4.66)and choose /t’’ < This formula assumes that the fractional rate of energyrelease varies slowly from cycle to cycle. The main timestep used for cycle n + 1 is thengiven by:min(zt’, LtmX, IS.tmax) (4.67)The timestep tt1 limits the amount by which the timestep can grow:= 1.02& (4.68)Timesteps 1it and are, respectively an input maximum timestep and a maximum timestep based on an input maximum crank angle in engine calculations. TheChapter 4. NUMERICAL MODELLING 90initial guess for Lt on cycle 0 is given by input quantity DTI. This is then comparedwith Ltmx and Ltmxca to determine the actual initial timestep t°. If the DTIsupplied on a subsequent restart differs from the DTI at cycle 0, the current Lt will bereset to this new DTI.The convection timestep Lt is based on the Courant stability condition. In a rectangular mesh, this condition is:( Lx Iky______It mm ,,(4.69)\IU—bI Iv—bI w—b21where b, b and b are the components of the grid velocity Constraint (equation 4.69)limits the magnitude of the flux volume in any direction to a value less than the cellvolume. To generalize this to an arbitrary mesh, it is natural to replace equation 4.69 by• the similar criterion:zt—’ t’ mm Vk (4.70)where the SVa are the flux volumes calculated for cell (i, j, k) using timestep Inpractice, for accuracy we also reduce the timestep determined from equation 4.70 by afactor f, typically taken to be 0.2.Chapter 5EXPERIMENTAL RESULTS AND DISCUSSIONIn this chapter the results of a detailed experimental study of mean and turbulenceflow field near the top dead centre position are presented. Laser Doppler velocimetrymeasurements were performed in the single-cylinder of the University of British ColumbiaRapid Intake and Compression Machine (UBCRIM), shown schematically in Figure 3.1and described further in Table 3.1. Optical access for backscatter LDV is provided by alarge quartz window in the cylinder head. Four tangential ports drilled on the chamberwall serve for intake, while a fifth port also drilled on the chamber close to the top deadcentre position serves as the exhaust.The specific measurement locations and engine operating conditions are described insection 5.1. The engine was motored at 1000 revolutions per minute and the intake swirllevels were varied and their effect on near top dead centre mean and turbulence flows arediscussed in section Measurement Locations and Engine Operating ConditionsThe tangential (swirl) and radial (squish) components of velocity were separately measured at three depths below the cylinder head along a radial direction. The tangentialports produce counter clockwise intake swirl, (when viewed from the cylinder head). Thisis regarded as positive swirl, that is, V> 0, positive radial velocity, U> 0 is defined asinward towards the cylinder axis. Two swirl levels (0.5 and 2.0) characteristic of directinjection stratified charged (DISC) engines were used in the experiments. These intake91Chapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 92swirl levels were obtained by varying the inlet velocities through the ports. Their valueswere estimated by the steady flow method and the details of the calculations involvedare described in appendix C. The swirl level (ratio) is defined as the equivalent angularvelocity of the cylinder flow field divided by the engine angular velocity.The cylinder axis is defined as the z-axis with the cylinder head as z = 0. Axialdistances are non-dimensionalized by the distance H = 14.0 mm from the head to thebowl bottom when the piston is at top dead centre. Similarly, radial distances are scaledby the bowl radius, R = 29.0 mm for bowl-in-piston and R = 22.0 mm for re-entrantbowl. All measurement locations lie within the top dead centre piston-bowl volume.Figure 5.1 shows the locations of the measurement planes relative to the piston at topdead centre. The planes are referred to as A, B and C.The two piston geometries shown in Figure 5.2 were investigated. The simple standardbowl-in-piston Figure 5.2a is similar to the relatively deep combustion chambers of late-injection direct stratified charge engines. The re-entrant bowl Figure 5.2b correspondsto the geometry of some early-injection direct stratified charge engines. Also shown inFigure 5.2 is the grid of Laser Doppler velocimetry measurement points used to mapthe flow field in the three axial planes considered. It is pertinent to mention that themeasurement locations are fixed relative to the cylinder head and do not move withthe piston. For the standard bowl-in-piston, the measurement depths were z = 3.0, 6.0and 9.0 mm below the head, the same axial planes were used with the re-entrant bowl.These planes were carefully chosen so as to facilitate data (experimental) comparisonwith numerical results.Four combinations of bowl geometry and intake swirl level examined are:• The re-entrant bowl with low swirl (0.5),• The re-entrant bowl with high swirl (2.0),Chapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 93• The bowl-in-piston with low swirl (0.5), and,• The bowl-in-piston with high swirl (2.0)Laser Doppler Velocimetry data (10 cycles) were collected with a 1 .0-degree resolutionfrom 60 degrees before top dead centre to the top dead centre position, for all the casesexamined. The effects of sample size was investigated by Mawle (1989). Mawle (1989)• varied the sample size N between 2 and 15 cycles and found that the ensemble-averagedvalues (mean and turbulence intensity) remained essentially constant when N was 10or greater. The result of a single experimental run is thus a history of instantaneousvelocity versus crank angle throughout the engine cycle. The data collected were ensembleaveraged off-line using a dedicated computer and software supplied by Dostek Inc. Theuse of ensemble averaging over 10 cycles implies that the data do not characterize thevelocity behaviour within individual cycles.5.2 ResultsThe results of the temporal evolution of the local (radial and tangential) velocity components and the respective components of turbulence intensity, measured from inlet portclosure (60 degrees BTDC) to the top dead centre position (0 degree) for some of therepresentative cases investigated are presented in this section. Adjacent measurementtimes are separated by 5 degrees crank angle. The results presented for the measured radial and tangential components of velocity and their respective turbulence intensities are• expressed in terms of the mean piston speed, V, which at 1000 revolutions per minuteis equal to 2.33 metres per second.The discussion of the results are focused on three flow regions, namely;• Measurement position which is neither dominated by squish nor located on thecylinder axis;Chapter 5. EXPERIMENTAL RESULTS AND DISCUSSION• 94a• Measurement location close to the bowl entrance region, and;• Location on the cylinder axis.5.2.1 Mean Radial VelocityThe variation of measured mean radial (squish) velocity with crank angle at locationswhich are neither in the squish-dominated region nor on the cylinder axis are shown inFigures 5.3 and 5.4, for the two chamber configurations investigated. Squish is the namegiven to the radially inward or transverse gas motion that occurs toward the end of thecompression stroke when a portion of the piston face and cylinder head approach each• other closely. The amount of squish is often defined by the percentage squish area; thatis the percentage of the piston area, 7rB2/4, which closely approaches the cylinder head,see Figure 3.3. The measured squish velocities in both combustion chambers at variouscrank angles before 30 degree BTDC are just a few percentage points different from eachother in magnitude. This is to be expected since the effect of squish is not really much atthese points in the compression process. As the piston moves past the 30 degree BTDCposition, a steep rise in squish velocity was observed in both chamber configurations, withthe re-entrant bowl having superior squish generating characteristics. It is also interestingto note that swirl has little or no effect on the peak values of squish; moreover, higherswirl levels only smooth out the peak values of squish as shown in Figures 5.3b and 5.4b.These observations are consistent with the fact that in the absence of swirl (low swirl0.5), the strong radially inward flow (squish) dominates, while with high swirl (swirl level2.0), the squish and swirl flows interact forcing a redistribution of the flow field. Thecentrifugal forces associated with swirl oppose the radially inward squish flow, and resultis a very complicated flow pattern. More details of the resulting flow pattern betweenthe effects of squish and high swirl in particular is presented in chapter 6.Chapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 95Figures 5.3 and 5.4, also show that the maximum values of squish velocity occur atabout 10 degrees before top dead centre. The magnitude of the measured squish velocitydepends on the bumping clearance and squish area and is also affected by air leakagethrough the piston rings and heat transfer losses to the wall and cylinder head. It canalso be seen from Figures 5.3 and 5.4 that between 10 degree before top dead centre andtop dead centre, the radial velocities all fall-off as the piston reaches the end of its travel.5.2.2 Mean Tangential VelocityThe evolution of measured mean tangential velocity with crank angle at locations whichare neither in the squish-dominated nor on the cylinder axis are plotted in Figures 5.5and 5.6. The mean tangential velocities are similar in trend, except for the expecteddifference in magnitudes due to the two swirl levels and piston bowl configurations. Itcan be seen from these figures, that at the top dead centre position, the re-entrant pistonbowl configuration shows a value of tangential velocity which is about 36% more thanthat of the simple bowl-in-piston arrangement. This could be due to the fact that thebowl-in-piston has more surface area than the re-entrant bowl, hence higher frictionaldissipation rate. The somewhat slow response of the mean tangential velocity at between60 to about 30 degrees before top dead centre of compression could be attributed to thefact that the flow field was still undergoing a process of development. As the top deadcentre position is approached, angular momentum conservation causes the air to rotatefaster as it is been forced into the piston bowl and its moment of inertia reduced. Theeffects of dissipation (friction at the walls, turbulence dissipation within the fluid andvelocity gradient) and modification of the swirl profile by squish, however, contribute tothe decay of the mean tangential velocities after they reach their maxima at between 10and 8 degrees before top dead centre.Chapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 96Mean tangential velocity peaks occurring slightly before top dead centre have been observed in several other laser Doppler velocimetry studies, Monaghan and Pettifer (1981),Fansler and French (1988), Beacon et at (1981). Other investigations, Fansler and French(1987) including some hot-wire measurements, Brandi at at (1979), have shown that swirlvelocity peak occurs at or after top dead centre. Intuitively, however, one might expectdissipation to cause the tangential velocity peak to precede top dead centre somewhat.As the results in the present study show, the in-bowl flow field is sufficiently complexthat one should not be surprised if timing of the peak tangential velocity varies betweendifferent engines or even between different locations in the same piston bowl. Variationsin peak timing of mean tangential velocity were observed in Brandi et at (1979).5.2.3 Mean Velocity Close to Bowl Entrance RegionFigure 5.7 presents data on both the swirl and radial velocities for a location where thefluid motions are dominated by the effect of squish. The swirl behaviour shown in Figure5.7 is similar to those seen in Figures 5.5 and 5.6 in that between 60 and 30 degree beforetop dead centre of compression, the response to compression was somewhat sluggish,probably due to the fact that the flow was still undergoing a process of developmentcoupled with the fact that the effect of squish was not appreciable at this point in theengine cycle. From about 30 degrees before top dead centre of compression onwards theswirl velocity shows a sharp increase in magnitude peaking at about 10 degrees beforetop dead centre and then decaying to about 25 metres per second at top dead centre.The radial velocity shown in Figure 5.7, behaves precisely as one would expect of asquish flow. Beginning at about 32 degree BTDC for the re-entrant bowl configurationand later at about 26 degree BTDC for the bowl-in-piston chamber, a radially inward(U > 0) flow develops as the rising piston forces the air in the squish region into the bowl.The largest inward velocity of about 20 meters per second for the re-entrant bowl occursChapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 97at 10 degree BTDC, while that for the bowl-in-piston was about 16 metres per secondand it occured at 15 degree BTDC. In the re-entrant bowl the maximum squish and swirlvelocities occurred at about 10 degree BTDC, while the two velocity components wereout of phase in the bowl-in-piston design. It is generally desirable for the swirl and squishvelocity peaks to be out of phase, so that advantage can be taken of their turbulencegenerating abilities at different points within the engine cycle, particularly around the topdead centre position. The bowl-in-piston configuration clearly demonstrates this desiredfeature. From Figure 5.7, we see a decline in the squish velocity after the maximumvelocity point. From this point onwards the radial velocity decreases in magnitude toabout 14 and 9 metres per second at top dead centre for the re-entrant and bowl-in-pistonchamber respectively.5.2.4 Velocity Evolution on Cylinder AxisAt a location remote from the region dominated by the piston-bowl wall and the squishjet, the variation of the tangential and radial velocities with crank angle becomes morecomplicated as depicted in Figure 5.8. The magnitude of the velocity components measured on the cylinder axis were seen to be consistently lower than those acquired elsewhere(see Figures 5.3, 5.4, 5.5, 5.6 and 5.7). The velocity components were characterized bynoticeable oscillations which were somewhat in phase with one another throughout theengine cycle. These oscillations could be attributed to the combined effects of swirl centreprecession, Arcoumanis et al (1987) and the increasing angular velocity of the gas as thesquish flow forces it into the bowl. Another possible reason for this observation could bethe effect of recirculation, resulting from the interaction of the squish flow and the mainflow (see Figures 6.17 to 6.19).Chapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 985.2.5 Radial Turbulence IntensityThe variation of radial turbulence intensity with crank angle at locations which areneither in the squish-dominated region nor on the cylinder axis for the two chamberconfigurations are shown in Figures 5.9 and 5.10. The figures show that the turbulence levels decrease after inlet valve closure (60 degrees BTDC), until about 50 degreesBTDC, thereafter the turbulence levels in both chambers remain relatively constant until at about 25 degrees BTDC when the effect of squish takes over. It is also worthmentioning that the levels of turbulence were generally higher in the re-entrant bowlconfiguration. Figures 5.9a and 5.lOa also show that at low swirl levels, the levels of turbulence between 50 and 25 degree BTDC in both chambers were more or less the sameafter the initial decay experienced after inlet port closure. This observation reinforcesthe need of having a substantial intake swirl during engine operation.The results shown in Figures 5.9 and 5.10, show that the turbulence intensity levelsare higher in the lower planes (z = 9.0 mm) than at the upper planes (z 3.0 mm)when the intake swirl levels are low and higher in the upper planes when intake swirllevel is high. This observation could be due to the combined effect of high intake swirland squish. The practical importance of this finding could be in pointing the way onwhere to locate the spark plug or in identifying the optimum time within an operatingcycle at which to initiate ignition.5.2.6 Tangential Turbulence IntensityShown in Figures 5.11 and 5.12 are the evolution of measured tangential turbulenceintensity with crank angle at locations which are neither in the squish-dominated regionnor on the cylinder axis for the two chamber configurations investigated. Figure 5.11shows a steady drop in turbulence intensity after inlet port closure (60 degree BTDC)Chapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 99until about 45 and 35 degrees before top dead centre for the re-entrant and bowl-in-pistonchamber configurations respectively. At intake swirl level of 0.5, Figures 5.lla and 5.12ashow that between inlet port closure and 30 degree before top dead centre the turbulenceintensity was very much stable after the initial decay. This observation shows that lowlevels of intake swirl really do not contribute much to turbulence generation. The sharpincrease in turbulence intensity noticed from about 25 degree before the top dead centreonwards could be largely attributed to the effect of squish. On the other hand, Figures5.1 lb and 5. 12b show that high intake swirl levels are necessary for enhancing the levelsof turbulence throughout the engine operating cycle. Figures 5.11 and 5.12, show that apositive correlation exists between the levels of intake swirl and turbulence generation;moreover, higher intake swirl levels result in higher top dead centre turbulence intensity.The data presented in Figures 5.11 and 5.12 which correspond to measurements taken attwo different axial planes show that higher levels of turbulence were prevalent at a lowerplane (z = 9.0 mm) for high intake swirl levels, whereas turbulence intensity levels weremore or less the same at an upper axial plane (z = 3.0 mm) when the intake swirl levelis low. These observations could be attributed to the resulting vigorous activities whichensue as a result of the interaction between high intake swirl and the squishing action ofthe piston, whereas in the low intake swirl situation the only action of any significanceis squish and this is largely effective in the higher axial planes, hence the near identicalturbulence intensity measured.5.2.7 Turbulence Intensity Close to Bowl Entrance RegionThe evolution of radial and tangential turbulence components with crank angle measuredat a point where the effect of squish velocity is predominant in the two combustionchamber configurations studied are shown in figure 5.13. In Figure 5.13a, it can seenthat after the initial decrease in both components of turbulence intensity the squish andChapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 100intake swirl effects become more dominant and this results in continuous increase ofturbulence until top dead centre. Figure 5.13a also show that the tangential componentof turbulence fluctuates for quite some time (up to 20 degree BTDC) before it finallyshows a sharp increase which lasted until top dead centre. The fluctuation shown by thetangential component could be due to the effect of vortex precession, Arcoumanis et al(1987). In Figure 5.13b, the radial component of turbulence intensity was observed tooscillate over a considerable period of time during the engine operating cycle, but thetangential component was more or less steady with some occasional spikes.The data plotted in Figure 5.13 show that the radial component of turbulence wasconsistently higher than its tangential counterpart. This observation suggests that theflow field in the squish-dominated region is highly anisotropic.5.2.8 Turbulence Intensity On The Cylinder AxisThe variation of radial and tangential components of turbulence with crank angle measured on the cylinder axis for both combustion chambers studied are presented in Figure5.14. The trends shown by the data of Figure 5.14 are very different from those displayedin Figure 5.13. The turbulence intensity data presented in Figure 5.13 were generallyhigher than those of Figure 5.14, except at the top dead centre position where the results were found to be almost identical in the bowl-in-piston configuration (see Figures5.13b and 5.14b). The data shown in Figure 5.14 also support the observation that theturbulence intensity measured on the cylinder axis is non-isotropic.5.3 Measurement UncertaintyIf it is assumed that the system is measuring on a real Doppler burst, the error in thecounter measurement of time is 0.25% at 20 MHz, (TSI user’s Manual). The digitalChapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 101data system uses 16 bits, for approximately 0.25% accuracy. The beam angle was measured by Thermo-System Inc., (TSI) and specified to 0.1%. The Bragg-cell frequencyshift introduces some error due to the frequency uncertainty of the shift and this erroris approximately 0.1%. The total error arising from the Doppler burst measurement wastherefore approximately 0.7%, a very negligible value. The absolute crank angle is accurate to approximately 1.0 degree with 0.14 degree resolution. The uncertainty in themeasurement-location co-ordinates is about 0.2 mm.There were also some influence of both fringe and velocity bias on the measurements.• Fringe Bias — The probability of a particle generating a measurable signal dependson the direction of the particle relative to the laser beams (fringes). For fringe bias,an extreme example would be a particle moving parallel to the fringes. Under thissituation, no measurement would be recorded. Ideally, it is required that particlestraverse the fringes at right angles.• Velocity Bias — The probability (per unit time) of a particle entering the measurement volume depends on the fluid velocity. With a spatially uniform distributionof particles, it should be evident that faster moving particles will most likely crossthe measurement volume more frequently than slower moving particles.Extreme care was taken to ensure that the measurements reported in this thesis werecorrected for both forms of bias. This was achieved by ensuring that the fringe patternswere properly aligne4 in the desired orientation (vertical and horizontal). The sourceof error with the greatest potential for causing large uncertainty is the measurement ofnoise resulting from stray scattered light. This error source was minimized to the barestlevel by painting the piston top black.The ensemble averaged mean velocity and root-mean-square (rms) velocity fluctuations were calculated from a limited number of velocity measurements (10 cycles) andChapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 102were hence, subject to statistical sampling uncertainty. In general, the error in the meanis proportional to the ratio of the relative rms fluctuation to the square root of the numberof samples thus:meanu,,s.ertainty ()=) (5.1)Similarly, the error in the rms fluctuation is inversely proportional to the square root ofthe number of samples, that is,1 (5.2)As discussed above, the random uncertainties in the mean and root mean square(rms) velocities were estimated to be “-i 6% and 30% respectively.5.4 Discussion of ResultsThe experimental results presented so far provide evidence that turbulence is producedaround top dead centre in high-shear-regions formed near the piston-bowl rim. Thesehigh-shear regions are established by the interaction of squish and swirl. The resultsalso show that intake-generated turbulence has significantly decayed so that its influenceon top dead centre turbulence is negligible (see Figures 5.9 to 5.14). It was consistentlyevident that higher levels of shear and turbulence were obtained with the re-entrant bowlconfiguration. Generally, turbulence around top dead centre is more intense in the reentrant bowl under all swirl conditions than with either swirl level in the bowl-in-pistoncombustion chamber.These results are consistent with other Laser Doppler velocimetry studies, Vafidis(1984), Fansler and French (1988), which indicate the crucial role of bowl geometry indetermining the turbulence intensity near top dead centre position. Turbulence intensities produced around top dead centre of compression have important implications forChapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 103combustion, since they are essential for augmenting local transport and mixing processes;moreover, the right intensity and scale of top dead centre turbulence can reduce combustion duration and cycle-by-cycle variation.The measurements reported here and most other studies, Vafidis (1984), Fansler(1985) were all obtained with Laser Doppler velocimetry, with the data averaged overmany engine cycles (10 cycles in the present study). Neither length- nor time-scale datawere available in these studies. Furthermore, the data presented in these studies were notcycle-resolved. An interesting question that arises then, is to what extent do the measuredroot mean square velocities near the bowl entrance and elsewhere within the combustion chamber represent small-scale, high-frequency random fluctuation (turbulence) thatmight wrinkle a flame and enhance combustion during each cycle, rather than large-scale, low-frequency fluctuations (on the order of the chamber size and the engine-cyclefrequency). Furthermore, if indeed the measured quantities are true turbulence, thentwo other issues of interest are: whether they are homogeneous (that is, uniform) andwhether they are isotropic (that is, independent of direction). Before answering thesequestions, it is necessary to mention that the influence low-frequency or cycle-to-cycle velocity fluctuations are also important to engine operation, but the manner by which theyaffect engine combustion are different from those of small-scale, random high-frequencyfluctuations.In responding to the first question, it is necessary to note that, the geometry of thepiston-bowl must force the flow field near the wall to become a.xi-symmetric. The datapresented in this thesis do not contain measurements after top dead centre, but we can seefrom the results shown that some amount of symmetry particularly around the top deadcentre position would have resulted when the expansion stroke is factored in. It followsthen that ensemble root mean square velocity fluctuation near the bowl wall shouldhave minimal contributions from swirl-centre precession that is out of phase with theChapter 5. EXPERIMENTAL RESULTS AND DISCUSSION 104engine cycle or from low-frequency (cycle-to-cycle) random fluctuations in the swirl-centrelocation. Furthermore, if the low-frequency or cyclic-variation contributions dominate themeasured root mean square fluctuations, then the measured turbulence intensity shouldfollow the crank angle variations of the local spatial and temporal velocity gradients. Themeasurements reported at various depths and radial distances, however, indicate that theroot mean square fluctuation peaks show a tendency of being consistently narrower thanthe intervals over which significant gradients exist.Additional proof will now be presented to show that, at regions remote from the swirlcentre, the ensemble root mean square fluctuations presented in this thesis are indeedturbulence and are not dominated by large-scale, low-frequency fluctuations. The proofis based on the work of Mawle (1989).Mawle (1989) used cycle-resolved hot-wire anemometer measurements to investigatethe same engine employed in this study. Cycle-resolved velocity-fluctuation data formeasurement points identical to those considered in the present study were made (thatis, at locations near the bowl edge and at the centre of the bowl). The cycle-resolvedturbulence intensity was evaluated by averaging the individual-cycle velocity data with acut-off frequency of 500 Hz or a 12-degree wide window at 1000 revolutions per minute.The window averaged mean velocities for each cycle were taken to represent the meanvelocity at the centre of each window. A cubic-spline curve fitting routine with zerotension was then used to interpolate between the mean velocity points. The resultingcurve was interpreted as the mean velocity profile. The root mean square deviationof the instantaneous velocity within each window from this in-cycle mean velocity wascomputed over 10 cycles. The cycle-resolved velocity-fluctuation histories close to topdead centre presented in Mawle (1989) were found to be quite similar in magnitudes tothose obtained in this study (see Figures 4.18 and 4.19 in Mawle (1989) and Figures 5.13and 5.14 of present study).tdoiostiousnoauomotaat{tusiauajuqin1puflJD-uososje(eILoW2sini)pusadpij(unqan)suop,npn£pOpAUOpU1£2unbai;-tT‘ros-jrmspapuia&£pnssputpiodaisuui-ainsam&suu!aujnqintpj9flUOOUA&OsSUO!SSflSpt{tUOiJOTNOISSflOSIGUNYS17fl63717VLLN3Y1I?1JXIDChapter 6NUMERICAL RESULTS AND DISCUSSIONThe work reported here describes an exploratory step in combining boundary conditions obtained from laser Doppler velocimetry measurements with a multidimensionalflow calculation procedure to compute the three-dimensional air motion within the cylinder of UBCRICM. Also contained in this chapter are outlines of the methods and resultsobtained from the theoretical study.6.1 MethodsA detailed description of the numerical code used in the present study was presented inchapter 4 of this thesis, as a result only a brief summary will be given here. The KIVAII computer program solves the three-dimensional, unsteady conservation equations ofmotion of a chemically reactive mixture of ideal gases. It can also solve for the dynamics• of a liquid fuel spray and the coupling between the spray and the gas.The gas phase equations are solved by using finite difference approximations on acomputational mesh composed of arbitrary hexahedrons. The arbitrary mesh can conform to curved boundaries and can translate to follow changes in combustion chambergeometry. This method has two advantages, namely, the mesh need not be orthogonaland the same computer program can be used to compute efficiently flows of arbitraryMach number. In low Mach number problems, an acoustic subcycling method is used toimprove computational efficiency.The finite difference method consists of a process whereby each term in the partial106Chapter 6. NUMERICAL RESULTS AND DISCUSSION 107differential equations is replaced with an algebraic approximation obtained by integratingover a control volume of finite size. With the provision of suitable initial and boundaryconditions, the solution is advanced in time by solving the resulting system of algebraicequations at successive time steps. Further details of the current practices employed inthe KIVA-Il code can be found in Amsden et al (1989).In order to account for the effect of turbulence, two turbulence models, namely, k -and subgrid scale (SOS) models are embodied in the KIVA-Il code. When either of thetwo turbulence models is in use, two additional transport equations are solved for theturbulent kinetic energy, k, and its dissipation rate, e. The k-e turbulence model optionwas used in the present study.6.2 Modelling ApproachThe validated modified KIVA-Il code was used to compute the flow fields. It is pertinentto mention that the converted code retains all the salient features in Amsden et al (1989).The various modifications that were made to the code are described in appendix B.The first stage in the theoretical study was to obtain the various computationalmeshes for the relevant combustion chamber geometries. These were generated using themesh generating facilities of KIVA-Il code. Some of the mesh resolutions at 25 degreeBTDC in the circumferential, radial and axial planes for the re-entrant and bowl-inpiston combustion chambers are shown in Figures 6.1 and 6.2. The piston bowls weremodelled as follows:6.2.1 Flat Top PistonThe resolution in the r — z plane for the flat topped piston was 20 x 21 (radial x axial)mesh. The planes in the axial direction were varied from 21 at inlet valve closure (IVC)Chapter 6. NUMERICAL RESULTS AND DISCUSSION 108to 4 at top dead centre (TDC). The grid was designed to be uniformly spaced sincethere are no sharp discontinuities in the chamber shape.6.2.2 Re-Entrant Bowl PistonThe bowl was modelled using a 9 x 11 x 10 (radial x axial x circumferential) mesh.The region above the piston was modelled with a mesh that had 21 cells in the radialdirection, while the number of planes in the axial direction varies from 20 at inlet valveclosure to 3 at top dead centre.6.2.3 Bowl-In-PistonThe bowl was modelled using an 18 x 13 x 10 (radial x axial x circumferential) mesh.The region above the piston was modelled with a mesh having 27 cells in the radialdirection, while the number of planes in the axial direction varies from 29 at inlet valveclosure to 3 top dead centre.For ease of description of the results, planes of mesh cells in the axial, radial andcircumferential directions are referred to according to their I, J and K values respectively.Figures 6.3 and 6.4 show that near top-dead-centre only three planes of cells remainabove the piston face. As the piston moves upward, the axial dimension of cells abovethe piston face decreases. Each time this axial dimension becomes less than a smallpredetermined value (for example cell height at top dead centre), one plane of cells abovethe piston is deleted and the flow field is re-mapped onto the new mesh. This procedureimproves computational efficiency by maximizing the cell axial dimension and at thesame time reducing the number of redundant cells. Conversely, the deleted planes arerestored after top dead centre when the piston is on its downward stroke. Cells lying inthe piston bowl are simply translated with the piston motion.Due to severe memory limitations, the various resolutions described above were the• Chapter 6. NUMERICAL RESULTS AND DISCUSSION 109most practicable that could be realised on the SUN SPARC computer system.6.2.4 Initial ConditionsThe code was initialized at 60 degree BTDC and computations were terminated at 20degree ATDC. The initial temperature and pressure were set at 300°K and 1.0 barrespectively. The initial axial velocity distribution was assumed to vary linearly betweenthe piston and cylinder head, being zero at the cylinder head and equal to the pistonspeed at the piston head.At the boundaries of the computational domain, heat and momentum losses are calculated by assuming boundary layers are turbulent and matching to a law-of-the-wallvelocity profile. The details of this procedure are given in chapter 4 of this thesis.The initial swirl velocities, w(r), are related to the distance r from the cylinder axisby:w(r) = UtJi ( * (6.1)where .11 is the first order Bessel function, R, is the cylinder radius, and A is aconstant with a value of 3.11. The form of this profile and the value of A are taken fromexperimental measurements of swirling flows in cylinders (see Amsden et al (1985a)).The value of U is taken to be:624J2(A)where SR is the swirl ratio, is the engine angular velocity in radians per second(RPM *27r/60), and J2 is the Bessel function of order 2. This value of U gives the sameangular momentum as a solid body swirl with angular velocity (SR) * ftChapter 6. NUMERICAL RESULTS AND DISCUSSION 110The initial turbulent kinetic energy was assumed to be uniform and equal to 10% ofthe square of the mean piston speed, while the initial turbulent length scale was assumedto be 40% of the distance from each cell to the closest solid boundary. The k-c modelconstants were the same as those reported in, Amsden et al (1989).6.2.5 Numerical AccuracyIt is well known that the quality of flow code solutions can be strongly degraded bynumerical or false diffusion. In order for this to occur, three conditions must be met.First, the numerical diffusion coefficient must be of the same order or larger than thephysical diffusion coefficient. Secondly, gradients must exist in the true solution in orderto generate numerical diffusive fluxes. Finally, the diffusive fluxes must be of the sameorder or larger than other mechanisms in the governing equations (for example, convectivefluxes, source terms). As indicated by Roache (1982), the ratio of the numerical diffusioncoefficient to the physical diffusion coefficient for donor cell (upwind) differencing isproportional to the cell Reynolds number.Rejj= U XZ2 (6.3)where u’ is the turbulence intensity, tX the cell size and 11 the kinematic viscosity.For donor cell differencing, the numerical and physical diffusive fluxes are equal whenthe cell Reynolds number has a value of two. The cell Reynolds number is often used asa measure of numerical diffusion since if it is small (much less than two), false diffusionmust be negligible. However., if this condition is not met, false diffusion effects may ormay not be significant, depending on the other two criteria.The KIVA-Il code offers various types of convection schemes. These include quasisecond-order upwind (QSOU) differencing, interpolated donor cell and partial donor cellChapter 6. NUMERICAL RESULTS AND DISCUSSION 111(PDC) differencing schemes. The quasi-second-order scheme uses a flux limiting algorithm to vary the convective flux between the values given by donor cell and interpolateddonor cell approximations to make the scheme monotone. The calculations presentedin this thesis were performed using the interpolated donor cell technique, which has theadvantages of greater accuracy (second order versus first order) and smaller numericaldiffusion over the upwind scheme. As such, the Reynolds number at which the numericaland physical diffusive fluxes are equal is somewhat greater than two.In order to investigate possible degradation of the predictions due to numerical diffusion, computations were performed with grids of different sizes and tirnesteps. Themesh designs listed above produced predictions which were grid-resolution independent.A computational timestep of 1.04167 x 10 and smaller was found to produce stableresults.6.2.6 Normalization of ResultsIn order to facilitate comparison with data from present experimental studies and othermodelling studies, the computed results were normalized. The mean velocities and turbulence intensities were non-dimensionalized by the mean piston speed. The turbulencelength scale was normalized by the clearance height, while the turbulence diffusivity (viscosity) was scaled with the product of the clearance height and the mean piston speed.6.2.7 Parametric StudiesOn successful completion of the validation exercise, (see Figure B.1) the converted codewas then used to investigate the performance of three different combustion chamberdesigns. The indices of performance forming the basis of the present study were the meanradial and tangential velocities and various turbulence parameters (turbulence intensity,length scale and diffusivity) emanating from the resulting fluid flow as a consequenceChapter 6. NUMERICAL RESULTS AND DISCUSSION 112of the different combustion chamber configurations and intake swirl levels (0.5 and 2.0).These turbulence parameters were selected because various turbulent flame speed modelsand most qualitative descriptions of flame structure are based on them, Abdel-Gayed etal (1984), Andrews et al (1975). The converted code was also used to investigate theeffect of swirl ratio on the flow and turbulent quantities listed above.6.3 Results and DiscussionThe results of some of the numerical simulations obtained from the converted KIVA-Ilcode will be presented in this section. The sample calculations represent typical simulations of portions of the compression and power strokes of Spark Ignition (S.I) engines.On the intake stroke, pure air is drawn into the cylinder through four tangentially placedinlet ports, causing the air to swirl uniformly about the axis of the cylinder. The combustion model embodied in the KIVA-Il was disabled, hence giving rise to cold flowsimulations. Apart from switching-off the combustion model, all other engine operatingcycles, namely; compression, power and exhaust strokes were simulated. These calculations are for a spark ignition engine with two swirl levels (0.5 and 2.0) and percentagesquish of 80, for both the re-entrant and bowl-in-piston chamber configurations. Threecombustion chamber configurations, were investigated. They are:• Disc-shape chamber - that is, a chamber having both the cylinder head and thepiston crown fiat, this is also referred to as a simple pancake combustion chamber,• A complex combustion chamber - having a flat cylinder head and a re-entrant bowlin the piston crown, and,• A complex combustion chamber - having a flat cylinder head and bowl in the pistoncrown, this is the standard bowl-in-piston configuration.• Chapter 6. NUMERICAL RESULTS AND DISCUSSION 113The computational mesh of these pistons are shown in Figures 6.1 to 6.4. The radial andtangential velocity components and various turbulence parameters emanating from thesecombustion chamber configurations were used as the basis for evaluating the combustionchamber effectiveness. The velocities shown in the present work were area averaged inplanes between the piston crown and cylinder head.6.3.1 Mean Radial VelocitiesComputed mean radial (squish) velocities resulting from various combustion chamberconfigurations are plotted in Figure 6.5a. and 6.5b. From the figures we see that thedata from all the chambers have a similar trend. The calculations also show that themore complex a chamber is the higher the squish flow it can generate. This feature isvery important for mixing and can result in good fuel economy and lead to acceptableemission levels.The peak values of the computed squish velocities occur at about 10 degree BTDCand they are slightly higher than their measured counterparts (see Figures 5.3 and 5.4).A possible cause of this discrepancy could be attributed to the fact that, the effects ofleakage past piston rings, heat transfer, friction and gas inertia are not adequately takeninto consideration in modelling squish. From Figures 6.Sa and 6.5b we see that higherpeak values occur for low swirl flows. This is due to the fact that, in low swirl flows, theradial velocities encounter a lower centrifugal pressure barrier created by swirl.6.3.2 Mean Tangential VelocitiesMean tangential velocities generated by the re-entrant and bowl-in-piston combustionchambers are shown in Figure 6.6. It can be seen from Figures 6.6a and 6.6b that bothcombustion chambers have some amplifying effects on the tangential velocity. We alsonote from Figure 6.6, that the tangential velocity is lower in the bowl-in-piston than inChapter 6. NUMERICAL RESULTS AND DISCUSSION 114the re-entrant bowl piston configuration. This could be attributed to the fact that the• bowl-in-piston chamber presents more surface area to the flow field, hence the tangentialvelocity within it experiences a higher level of frictional force. The trends exhibited bythe computed mean tangential velocities are in very good qualitative agreement with theexperimental results (see Figures 5.5 and 5.6). As expected, higher swirl velocities wereobtained throughout the compression stroke when the intake levels of swirl are high.6.3.3 Intake Swirl Decay RateFigure 6.7 shows the rate of decay of intake swirl for the combustion chamber designs.The intake swirl decay rate is the rate of change of the ratio of the angular momentum ata given crank angle about the cylinder axis to the angular momentum at exhaust valveclosure. This parameter represents the fractional rate of decay of the initial angularmomentum remaining at each crank angle. Figure 6.7 shows that the decay rate decreases rapidly from inlet port closure (60 degree BTDC) until about 54 degree BTDC,thereafter the decay rate remained constant up to the top dead centre position for thebowl-in-piston and disc chamber designs. The sharp increase in the rate of decay noticed• between 20 and 5 degree BTDC for the re-entrant bowl could be attributed to the pronounced increase in the frictional forces acting on the swirl velocity around these pointsduring the compression stroke (see Figures 6.6a and 6.6b). The figures also show thatbeyon4 10 degree BTDC the rate of decay of intake swirl starts to decrease once againand this continues till the top dead centre position. This is presumably due to the relative decrease in the effect of wall friction in the squish area, since the angular momentumin the squish region is transferred into the bowl by the squish flow (see Figures 6.14 to6.23).Chapter 6. NUMERICAL RESULTS AND DISCUSSION 1156.3.4 Thrbulence IntensityComputed turbulence intensities for the various combustion chamber shapes are shownin Figure 6.8. In Figure 6.8, we see that the re-entrant and the bowl-in-piston chamberdesigns have an amplifying effect on turbulence and in particular, maintain the turbulence intensity levels at high values around the top dead centre position. Specifically, themaximum values attained by the turbulence intensities at the TDC position could beattributed to squish and spin-up mechanism, a process whereby the turbulent eddies tryto preserve their angular momentum upon being compressed. This fact is corroboratedby the corresponding decrease displayed by the turbulence length scale around the topdead centre position, (see Figure 6.12). The observed global increase in turbulence intensity could be attributed to the squishing effects induced on the charge by the complexcombustion chamber geometry (see Amann (1985)). The data in Figure 6.8 also showthat the combination of squish and high levels of intake swirl lead to higher values ofturbulence, particularly around the top dead centre position. Moreover, the data alsoindicate that the re-entrant bowl chamber consistently out-performed the bowl-in-pistondesign in terms of overall turbulence generation.Figures 6.9 to 6.11 show the effect of intake swirl level on domain averaged turbulenceintensity at various planes (A, B and C) in both the re-entrant and bowl-in-pistoncombustion chambers. The figures show that the levels of turbulence are higher in thelower planes. The data presented in Figures 6.9 and 6.11 also show that the peak valuesof turbulence intensities occur near the bowl lip. This observation could be attributed tothe fact that high levels of squish are generated near the bowl lip, and moreover, intakeswirl levels are higher in the lower planes, since they must spin faster in order to conservetheir angular momentum. It can also be seen from the figures that the distribution ofturbulence throughout the planes investigated are different.Chapter 6. NUMERICAL RESULTS AND DISCUSSION 1166.3.5 Turbulence Length ScaleThe values of turbulence length scale obtained for the three combustion chamber investigated are plotted against crank angle in Figure 6.12. Figure 6.12 shows that contrary toexpectations, the turbulence length scales resulting from the various chamber configurations do not scale with top dead centre clearance height, instead they were independentof it at the top dead centre position. The turbulence length scales were found to decrease monotonically over the crank angle interval investigated, with all of them showinga broad minimum between 5 degree before top dead centre and top dead centre. Thistrend is consistent with the observed maximum values of the turbulence intensities atthe top dead centre position (see Figure 6.8).6.3.6 Turbulence DiffusivityThe evolution of the turbulence diffusivity (locally proportional to the product of intensityand length scale) with crank angle is shown in Figure 6.13. Figure 6.13, shows that thediffusivity is highest for the re-entrant bowl configuration and least for the pancakecombustion chamber. This is not very surprising because the simple pancake chamberhas no additional means of supplementing the intake generated turbulence (see Clerk(1921)). The distribution of turbulence difFusivities as depicted in Figure 6.13, stronglysupport the idea of using complex combustion chambers, especially in diesel engineswhere turbulent diffusivities are very important in mixing during the ignition delay period(approximately 10 to 20 degree BTDC). The prevalence of larger diffusivities aroundtop dead centre could also have a beneficial effect in determining NO and hydrocarbon(HC) emissions. Furthermore, the large values of diffusivities from 10 degree BTDC toTDC could enhance a faster diffusion burn, and consequently improvement in specificfuel consumption arid smoke can be obtained, Abraham et al (1988). The trend shownChapter 6. NUMERICAL RESULTS AND DISCUSSION 117by the profiles in Figure 6.13 are consistent with the simplified analysis of diffusivity,D’—’u’t.6.4 Discussion of the Computed Flow FieldsA challenging aspect of three-dimensional calculations lies in visualizing the complexflow fields and in relating the different flow field properties to determine the causes ofthe observed flow features. In this section two-dimensional velocity vector plots and contour maps of scalar variables in selected planes that slice through the three-dimensionalflow field are presented. The velocity vector plots show the projections of the velocityvectors onto selected planes. By juxtaposing several such velocity plots, a reasonableunderstanding of the complex three-dimensional flow-field can be acquired. Figures 6.14to 6.45 show a sample of the various flow fields generated.6.4.1 Velocity Vector FieldsThe velocity fields generated by the re-entrant and bowl-in-piston combustion chambersfrom 25 degree BTDC to 20 degree ATDC are shown in Figures 6.14 to 6.33. Fromthese figures, we could see the strong squish motion generated by the re-entrant bowlconfiguration during compression stroke. These figures also show that, in the absence ofswirl (low swirl 0.5), the strong radially inward flow out of the squish region at thebowl circumference sets up a torroidal flow within the bowl as expected. This can beseen from Figures 6.14a to 6.21c for the low swirl cases. When swirl is present, however,the squish and swirl flows interact forcing a redistribution of flow’s angular momentum.The centrifugal forces resulting from the intake swirl oppose the radially inward squishflow, causing it to turn into the bowl earlier. The resulting vortex flow in the bowl thenrotates in opposite direction to the vortex in the low swirl case. These observations areChapter 6. NUMERICAL RESULTS AND DISCUSSION 118accurately captured in the data presented in Figures 6.14b to 6.21d.It is also evident from Figures 6.14a to 6.21b that, apart from the short squish penetration observed in the re-entrant bowl for high intake swirl levels, the resulting squishattaches itself to the bowl wall just below the lip. This re-attachment phenomena isnot particularly suitable for global turbulence generation, it only promotes localized highshear velocities around the lip of the bowl. The figures also show that because of theshort penetration range, there is hardly any interaction between the main flow and squishmotion, whereas a high level of interaction is evident between the main flow and squishfor the case of low or no intake swirl. A more serious drawback of this type of flow isthat, the fuel particles could become separated from the air, due to the fuel particlesin-ability (inertia) to follow such an abrupt change in air flow pattern. The evidenceof this could be seen as ifim of fuel deposition both on the wall and at the bottom ofthe cup. However, because of the somewhat narrow channel created by the front of thesquish motion, the main stream fluid, flows past in a jet-like fashion and hence generatessome turbulence.The observation raised above indicates that the re-entrant bowl configuration requiressome form of optimization as regards global turbulence generation. One possible way ofachieving this optimization is to design a re-entrant combustion chamber which generatessquish motions with zero re-attachment and longer penetration into the main stream flow.Where bowl modification is not feasible, then an optimal level of intake swirl must besought and used.From Figures 6.18a. and 6.18b, we notice that the re-attached squish motion hasdisintegrated into a form of clockwise motion, and at about the top dead centre position,Figure 6.19, the motion has reduced in both magnitude and intensity and is now uniformlydistributed throughout the combustion chamber.The squish motion generated by the bowl-in-piston configuration is not as intense asChapter 6. NUMERICAL RESULTS AND DISCUSSION 119that generated by the re-entrant bowl, but it does not re-attach itself to the chamberwall, moreover it has a longer penetration range than that prevailing in the re-entrantbowl. Because of its longer range it interacts more with the main stream flow and theresult is the very noticeable vortex generated in the fourth quadrant. The vortex motionis very suitable for enhancing the mixing of air-fuel mixture and it also promotes uniformturbulence generation. The mixing is very important as it improves fuel economy andreduces emissions.Figures 6.18c and 6.18d show the motion within the bowl-in-piston chamber levelling-off and at the top dead centre position, Figures 6.19c and 6.19d, show that the entirecombustion chamber is now filled with uniformly distributed motion.Figures 6.20 to 6.23 show the reverse squish motion during the expansion stroke. Thereverse squish motion is useful for scavenging operation as it helps clean the chamberprior to induction of fresh charge, thereby reducing charge dilution and consequent improvement in cycle-by-cycle variation. Figures 6.24 to 6.33 show the velocity vector plotsat two K planes, namely, 6 and 12. Figures 6.27b, 6.29b, 6.30b and 6.32b show thepronounced effect of squish in the re-entrant bowl configuration. This is evident fromthe slight inclination of the velocity vectors.6.4.2 Swirl Velocity FieldsFigures 6.34 to 6.37, show the swirl velocity contours for the various combustion chambersfrom 20 degree BTDC to 5 degree before top dead centre position. In these plots thecontour lines marked L show the location of the highest swirl velocities. As can be seenfrom the figures, the contour plots of the swirl velocity displays a characteristic rigid bodymotion (notice the almost parallel and equally spaced contour lines). The radially inwardtransport of fluid from the squish motion is responsible for the spin-up of fluid in thelower sections of the L contour line. The figures show that the swirl motion generated byChapter 6. NUMERICAL RESULTS AND DISCUSSION 120the re-entrant and bowl-in-piston chamber configurations are very similar in all aspects.6.4.3 Turbulent Kinetic EnergyThe temporal distributions of total turbulence kinetic energy integrated over the relevantchamber volume from 20 degree BTDC to 5 degree before top dead centre position areshown in Figures 6.38 to 6.41. In each of the figures, the contour line marked H representsthe maximum while that marked L corresponds to the minimum. The figures show thatthe maximum values of the turbulent kinetic energy extend throughout the central regionof the re-entrant bowl and are somewhat localized to interior portions of bowl-in-pistondesign throughout the compression stroke. Because of the shear flow induced by thesquish flow, one might guess that the highest turbulence intensities would be observednear the rim of the piston bowl. The plots show however, that the highest intensities arelocated in the interior of the piston bowl. This is because the dissipation rate of turbulentkinetic energy is inversely proportional to the length scale, £, through the relationship,k2/3e, and £ is very small near walls (see Figures 6.42 to 6.45). The observed highvalues of the turbulent kinetic energy around the lip of the re-entrant bowl design aredue to strong squish actions. Figures 6.38b to 6.41b, show that the area covered by themaximum value of turbulence kinetic energy in the bowl-in-piston configuration is verysmall. This observation could be attributed the relatively weaker squish generated by thebowl-in-piston configuration. At 10 degree BTDC, we see that the level of turbulencekinetic energy has undergone a drastic reduction in both combustion chambers, butstrong turbulence still persist around the lip of the re-entrant bowl configuration. Thisobservation could be attributed to the fact that, the squish motion generated by there-entrant combustion chamber is localized to the bowl lip region. The situation justdescribed is depicted in Figures 6.40 and 6.41.Figure 6.41 shows the distribution of turbulent kinetic energy in both combustionChapter 6. NUMERICAL RESULTS AND DISCUSSION 121chambers at 5 degree before top dead centre position. We notice that the turbulencekinetic energy has virtually dwindled to very low values everywhere within the combustionchambers except the nugget of high turbulence kinetic energy located around the bowl lipin the re-entrant chamber design. The figure also shows that though the intake generatedturbulence kinetic energy is almost completely dissipated by the action of viscous forces(diffusion), bowl geometries (particularly the re-entrant chamber configuration) was stillable to enhance turbulence somewhat. A possible explanation for this could be due tothe interaction of the stress field (squish) with the mean normal strain imposed by the• piston motion and by compression of the induction generated vortex.6.4.4 Turbulent Length ScaleThe contours depicting the distribution of turbulent length scale from 20 degree BTDCto 5 degree before top dead centre position are illustrated in Figures 6.42 to 6.45. Theprofiles of turbulent length scale show that the maximum values corresponding to the Hcontour lines occur at the centre of all the combustion chambers. The lowest values of theturbulence length scale are distributed around the periphery of the chambers close to thewalls. This is to be expected since turbulence length scales are determined by chambergeometry near top dead centre with the integral scale of the order of the clearance gapin disc chambers and of bumping clearance in squish chambers. The magnitude of theturbulence length scale is a very important parameter as regards early flame developmentin the sense that it can be consumed very rapidly thereby resulting in the formation of afavourable flame kernel. The concentration of low values of the turbulence length scalearound the lips of both the re-entrant- and the bowl-in-piston combustion chambers are• consistent with the high levels of turbulence intensities generated by the action of squishin these areas.Chapter 6. NUMERICAL RESULTS AND DISCUSSION 122Figures 6.18, 6.41 and 6.45 show the mean motion and turbulence kinetic energy (intensity) as decaying and the turbulence length scale increasing to the order of magnitudeof the bumping clearance, due to the action of turbulent dissipation in converting energyfrom the mean flow into heat and the preferential destruction of the smaller eddies, Tennekes and Lumley (1972). At this time, diffusion is causing the intensity field to becomemore homogeneous and the turbulence structure, as depicted by the length scale, is beingconstrained to assume the shape imposed by the geometry of the combustion chamber.6.5 DiscussionThe results presented so far show that multidimensional calculations could be very informative and useful as a tool for investigating engine flow fields. The results of the variouscomparisons (chapter 7) also indicate that it will be a very formidable task to adequatelyassess the accuracy of multidimensional models; on account of limited data and uncertainties about them. These constraints are due to difficult experimental conditions. Thepresent study and other investigations show that different experiments are of markedlydifferent value in assessing the validity of specific submodels.Although some difficulties exist concerning the ability of multidimensional modelsto accurately predict flow quantities, there utility in recent times, have reinforced theirversatility as a diagnostic method of the future. Some of the considerations supportingthis assertion are:• Reasonably-accurate predictions of the ensemble-average flow properties has beenshown to be possible by using existing statistical flux methods of turbulence modelling (k-e), provided that the numerical and boundary-condition uncertainties havebeen reduced to acceptable levels.• The speed, flexibility and accuracy of the numerical methods have steadily improvedChapter 6. NUMERICAL RESULTS AND DISCUSSION 123to the stage where three-dimensional calculations for practical chambers are nowfeasible; although costly and subject to some limitations.• The applications of multidimensional methods have increased and their ability toprovide information about trends and mechanisms which would be difficult andor costly to obtain by any other means, for example by experimentation has beendemonstrated.Some of the weaknesses of multidimensional methods are: inability to predict cycle-by-cycle variations, discretization errors, which require the use of fine grids even withhigher-order schemes, the problem of maintaining the required mesh connectivity properties, especially during the induction and exhaust phases and the burdensome computingoverheads of three-dimensional computations. A lot of these problems are computationalin nature and are not unique to engines, their counterparts exist in many other fields.For these reasons, research efforts are been intensified in these areas, so as to alleviatethe computational difficulties. An example of such research effort is in the area of multi-or adaptive-gridding.Other causes of error and confusion in multidimensional modelling are due to the factthat the complex set of coupled, nonlinear, partial differential equations do not exactlyrepresent the physical system under consideration. Consider for example, the use ofinput parameters such as chemical rates or diffusion coefficients. These quantities areoften used as submodels in multidimensional computations, and for them to adequatelyserve as an accurate input, they need to be derived from more fundamental theories,models or experiments. Experiences have however, shown that this is not the case— theymay not be known to any appreciable accuracy and often their values are simply guesses.Further challenges are posed by the general research areas of radiation transport andheterogeneous processes such as phase changes, droplet dynamics and soot formation.Chapter 6. NUMERICAL RESULTS AND DISCUSSION 124The most formidable of all the problems which multidimensional methods still have tocontend with is that of turbulence modelling, for which a general satisfactory frame-workfor its inclusion in hydrodynamics models is yet to be developed.Although multidimensional methods are presently limited in their ability to modelor predict engine flow dynamics due to various inadequacies already mentioned, theydo however, possess the potential for examining the trends of fluid-flow engine-geometryinteraction. The impact of these inadequacies is expected to decrease as these models andsolution techniques become more developed. It is encouraging to note that, in the longrun multidimensional methods are likely to be most useful in developing our conceptualunderstanding of the details of processes occurring inside the combustion chamber whichare important in determining engine operating characteristics. This optimism is rootedin the fact that, with time more knowledge and understanding regarding the constituentsubmodels will become available.In summary, the available evidence indicates that multidimensional techniques canbe regarded as being quite successful in providing useful insights into engine flow fielddiagnostics.Chapter 7COMPARISON OF NUMERICAL AND EXPERIMENTAL RESULTSIn this chapter computed values of mean radial and tangential velocities and turbulence intensities are compared with their measured counterparts. The reasons(s) forsome of the discrepancies between the experimental and computational results are given.The results obtained from this study were also compared (in terms of accuracy) to someinvestigations of similar scope and objective, with a view to determining whether thereare inherent flaws in the KIVA-Il code formulation or that various discrepancies thathave been reported during code validation are due largely to the inability to accuratelyspecify the initial and boundary conditions. The findings from the present comparison isexpected to shed more light on the accuracy or lack thereof of the fluid dynamic modelembodied in KIVA-Il code.7.1 Modification TechniqueBecause the k - turbulence model is based on the assumption that turbulence is isotropicand in view of the fact that only two velocity components were measured in the presentstudy, the measured turbulence intensity values had to be modified so as to providea consistent set of results for comparison with the numerical data. The modificationtechnique described by McKinley et al (1988) was adopted here. The details of themethod is as follows:Measurements were made of turbulence intensities in the radial and tangential directions. The turbulent kinetic energies corresponding to these intensities, k,.ad and ktan125Chapter 7. COMPARISON OF NUMERICAL AND EXPERIMENTAL RESULTS 126were then computed, thuskrad = 0.5 X (U’rad)2 (7.4)ICtan = 0.5 x (uan)2 (7.5)It was assumed that the axial turbulent kinetic energy, kaxjaj, was equal to the averageof the turbulent kinetic energies in the radial and tangential directions, hence;kaxiai = 0.5 X (k.rad + ktan) (7.6)With this artifice the total kinetic energy, k, was then calculated for each of themeasurement locations, thus:k = k7d + ktan + kaxwi (7.7)The experimental data were assumed to be symmetric; and anisotropy was accounted forby calculating the isotropic intensity, u50, corresponding to the total kinetic energy, k,thus:2 1/2= ( x k) (7.8)This modification shows that the comparison between the measured and theoretical data, is strictly speaking based on turbulent kinetic energy rather than turbulenceintensity. There was no attempt made to modify the mean velocities.Chapter 7. COMPARISON OF NUMERICAL AND EXPERIMENTAL RESULTS 1277.1.1 Mean Radial VelocitiesFigures 7.1 to 7.4 show the computed and measured mean radial velocities for the reentrant and bowl-in-piston chamber designs. From these figures one can conclude that thecomputed and measured mean radial velocities were generally within the experimentalaccuracy of each other, however, the measured values were 5 to 20% higher in magnitudethan the computed results. The figures also show that apart from the difference inquantitative values between the computed and measured results, the experimental trendswere well reproduced. The experimental investigations conducted by lijima and Bracco(1987) and Fansler and French (1988) showed measured squish velocities that were 20 to40% smaller than computed values and were asymmetric about top dead centre. Thesediscrepancies were attributed to inadequate heat transfer and blow by models embodiedin current flow codes,Iijima and Bracco (1987); viscous dissipation, the centrifugal effectof swirl, uncertainty in the actual top dead centre clearance height, the effects of valve-recess, and ring-land volumes were added to the list as possible culprits, Fansler andFrench (1987). The marked deviation of the measured velocities from simple two-zonemass transfer and two-zone squish models may not be due solely to the causes listedabove. Instead they may be the result of rapid decrease in the radial velocity as theradius decreases due to the abrupt turning of the squish jet into the bowl as top deadcentre is approached (see Figures 6.14a, 6.14b, 6.16a and 6.16b).• 7.1.2 Mean Tangential VelocitiesThe predicted and measured mean tangential velocities at two intake swirl levels for boththe re-entrant and bowl-in-piston chamber configurations are shown in Figures 7.5 to 7.7.The overall data presented show that the experimental trends are well reproduced by thenumerical simulation. The experimental values were however, between 10 to 25% higherChapter 7. COMPARISON OF NUMERICAL AND EXPERIMENTAL RESULTS 128than the predicted results. Figures 7.6 and 7.7 tend to suggest that, the agreementbetween the measured and computed values are better for deep planes. If this turns outto be the case, then it would be advisable to use measurements deeper in the bowls forcode validation.7.1.3 Turbulence IntensitiesThe computed and (modified) measured turbulence intensities were compared in planesA through to C for both chamber designs and intake swirl levels. Figures 7.8 to 7.13show the computed and measured turbulence intensities for both the re-entrant andbowl-in-piston combustion chamber designs.Figures 7.8a to 7.13a show that the experimental results are well predicted around theaxis of the re-entrant bowl configuration, but the accuracy gradually deteriorates towardsthe bowl lip. This could be attributed to the vigorous flow activities (shear and countershear) prevailing in this region. The figures also show that the location of maximumturbulence intensity were accurately predicted. The experimental trends were also wellreproduced by the numerical results. It is also worth mentioning that the distribution ofturbulence within the three planes investigated differ by between 10 to 15%.The data presented in figures 7.8b to 7.13b show that the values of turbulence intensities are poorly predicted on the axis of the bowl-in-piston chamber design. This couldbe due to complex flow pattern prevailing on the axis of the bowl-in-piston design (seeFigure 5.6b). As we approach the squish region, the experimental and computed turbulence intensity improves somewhat. We also see from these figures that the measuredtrends are reasonably well reproduced by the code. When the experimental uncertainties(‘s 30%) are taken into consideration, the glaring discrepancies between measured andcomputed results in Figures 7.8b to 7.13b are likely to reduce.Chapter 7. COMPARISON OF NUMERICAL AND EXPERIMENTAL RESULTS 1297.2 DiscussionThe studies conducted by Henriot et al (1989), showed measured turbulence intensitythat were smaller than predicted values. This discrepancy was attributed to such factorsas, truncation errors, uncertainty of inlet conditions and inadequate grid resolution. ZurLoye et al (1989) conducted a KIVA-TI code validation exercise similar to that of thepresent study, but instead of using measured mean velocities and turbulence parametersat inlet valve closure, in initializing the code, these quantities were estimated with the aidof zero-dimensional model. The results show a discrepancy of between 20 to 30% higherthan the measured values. The computed intensities were within 2% of the experimentalvalues at top dead centre and about 10% of the measured values during the expansionstroke.The differences between the computations and the measurements early in the compression stroke were thought to be due to the breakdown of the tumbling motion thatis present in the engine but not modelled in the computational code. It is however,pertinent to point out that, the effect of tumble breakdown would act to maintain theintensity during the compression stroke due to conversion of the mean flow kinetic energy into turbulence. Contrary to expectations, the model over-predicted the turbulenceintensity prior to top dead centre. This discrepancy was blamed on uncertainties in specifying the correct inflow length scale and turbulence kinetic energy at inlet valve closure.There was a close agreement between measured and computed results at top dead centreand beyond. This is to be expected, since the tumbling motion cannot contribute todeviations between the measurements and the computations, as it is expected to havebeen substantially dissipated at this point.Recent KIVA code validation studies undertaken by Amato et al(1988) showed resultswhich were very similar to those reported in this thesis. In the investigation of AmatoChapter 7. COMPARISON OF NUMERICAL AND EXPERIMENTAL RESULTS 130et al (1988) a method similar to that employed in this thesis (that is Laser Dopplermeasurements) was used for prescribing the numerical boundary conditions. This couldbe responsible for the very good accuracy between the experimental and numerical resultsobtained in the present study.Several other studies similar to the present investigation have attributed the discrepancies between the predicted and measured flow quantities to one or more of thefollowing;• Differences between the real and assumed inlet conditions;• Numerical truncation errors, and• Inadequacies in turbulence modelling, including the use of the law-of-the-wall asa bridging mechanism for near wall regions which are too small to be resolved inpracticable grid designs.Most workers using one form of numerical code or the other to predict engine flowand performance characteristics, have blamed the poor correlation between predictedand measured values on numerical truncation errors. The impetus for this has largelybeen based on the conventional method of assessing numerical accuracy, namely, gridrefinement tests. Systematic grid refinement tests for axisymmetric engine intake flowsusing numerical methods Harworth and El Tahry (1991) similar to those employed hereindicate that on the order of iO - i0 grid points (in three dimensions) may be required for mean velocity profiles to be everywhere within 10% of their converged values.Higher-order spatial differencing improves the situation somewhat. This example andsimilar studies Henriot et al (1989), Gosman et al (1984) lead to the conclusion that gridrefinement is an insensitive measure of numerical errors over the range of mesh densitythat can be afforded for three-dimensional calculations. In the work of Harworth andChapter?. COMPARISON OF NUMERICAL AND EXPERIMENTAL RESULTS 131El Tahry (1991) for example, it was found that the rate of convergence increases as themesh is refined; this implies that grid refinement tests on coarse grids could lead one tothe erroneous conclusion that numerical errors are small.Attempts to quantify numerical errors using Taylor-series-based methods are likelyto yield misleading results on coarse meshes (typically 3O or fewer grid points) thatare currently the norm in multidimensional in-cylinder flow modelling. However, Taylorseries analyses are appropriate only in the limit as the mesh spacing (or other parameterdetermining convergence) approaches zero. In complex geometries, such as in bowl-in-piston configurations, the effect of mesh topology (for example, spatial distribution ofgrid points, departure from orthogonality, etc) further complicates the interpretation ofgrid refinement tests.Another important source of error which is rarely mentioned is the various atypicalfeatures which are introduced into experimental engines. Some of these features include,mounting valves on the side of the cylinder, in the present study a view through thecylinder head was provided by replacing poppet valves with ports in the cylinder wall.This design was delibrate since the main objective of the present study was to create flowregimes which accurately model that embodied in the KIVA-Il numerical code. Squarepistons have also been used in some studies with hope of obviating optical distortions.It is worth mentioning however, that these various artifices create conditions which aredifferent from those prevailing in conventional engines and are therefore crucial sourcesof error.As already implied, reduction in numerical truncation errors is clearly a strong priority, but the results shown in Figures 7.1 to 7.13 of the present study and those of Amato etal (1988) strongly support the possibility that, if the initial and boundary conditions areaccurately specified, then numerical codes are capable of predicting engine flow dynamicsuioduiisomstl;suotpuo£ipunoqpuaujoUHPPOW1nDORtflSo1ompuajpojtds£jadoidmtaaa&&suottpuo3£iepunoqaiato&‘(6860ivo’ianz‘(686!)ivi1OUUH&qppnpuostpns-jatmuoujpM.£Jq?uosarsturgiunt1u3Tnidxapptoad‘psraid!SL7flS3&7VLLNI4II?IJXGNV7VOflTWflNJONOSflIVcTI4IOQ21OChapter 8CONCLUSIONS AND RECOMMENDATIONSA combined experimental and computational investigation of the in-cylinder flowfield near top dead centre in the UBCRICM was conducted. Detailed Laser Dopplervelocimetry measurements of the ensemble-mean velocity and root mean square (rms)velocity fluctuations in the motored engine was used to examined how two importantdesign parameters; namely: piston-bowl geometry and intake-swirl level affect the in-cylinder flow field near top dead centre of compression.In order to relate this study to real engines, for example, the direct injection stratifiedcharged engine (DISC) where the effects of intake swirl and squish are very important, abowl-in-piston and re-entrant bowl configurations, along with two intake swirl levels (0.5and 2.0) were used.The modified KIVA-II numerical code was also used to investigate various flow andturbulence quantities resulting from the different chamber configurations and intake swirllevels mentioned above.8.1 Conclusions• Based on these studies, the following conclusions were drawn:8.1.1 ExperimentalMean Flow133Chapter 8. CONCLUSIONS AND RECOMMENDATIONS 134• The squish flow at the bowl entrance is strongly affected by geometry and intakeswirl level.• The flow structure within the bowl is complex and depends strongly on both chamber geometry and intake swirl level.• The squish penetration is about 15% shorter for high intake swirl than for the lowintake swirl level.• The re-entrant bowl configuration was found to have superior flow generating capability, for example, squish and swirl were found to be 15 and 50% respectively moreintense in the re-entrant bowl piston geometry than in the bowl-in-piston assemblyat top dead centre.• Velocity measurements taken on the axis of both combustion chambers confirmsthe existence of swirl precession (see Figure 5.6).Turbulence Intensity• It was observed that the high-shear regions near the piston-bowl entrance producesignificant turbulence around top dead centre.• These shear regions were established before top dead centre by the interaction ofswirl and squish and these interactions were. further increased by the re-entrantgeometry (see Figures 5.3 and 5.4).• Under all intake swirl level conditions, the re-entrant chamber configuration wasfound to generate more turbulence than the bowl-in-piston assembly; at top deadcentre the swirl levels were about 30% higher.Chapter 8. CONCLUSIONS AND RECOMMENDATIONS 135• The velocity-fluctuations measured in this study are expected to contain a significant contribution from smaller-scale, higher-frequency fluctuations (turbulence)that can enhance combustion through increased transport, mixing and flame wrinkling. Evidence for this assertion comes mainly from the cycle-resolved measurements made in the same engine for a similar piston bowl shape, Mawle (1989).8.1.2 ComputationalThe modified KIVA-IT numerical code was also used to investigate the effects of differentcombustion chamber designs on various flow and turbulence quantities. Furthermore,the converted code was used to study the effects of swirl variation on these quantities.It is worth mentioning however, that the turbulence constants were not tuned. Based onthese studies, the following conclusions were arrived at:• The computer simulations show that the re-entrant bowl assembly generated moreturbulence than the standard bowl-in-piston configuration. In particular, Figure6.5 shows that the re-entrant bowl assembly generates about 35% more squish thanthe bowl-in-piston configuration at top dead centre. However, these squish flowsmay not result in the generation of uniform turbulence intensity within the wholedomain of the combustion chamber as they may become re-attached to the bowlwalls (see Figures 6.14a, 6.l4b, 6.15a and 6.15b).• The simulation results also indicate that squish penetration is about 20% more forlow intake swirl conditions (see Figures 6.14 to 6.18).• The response of the evolution of turbulence intensity to both combustion chamberconfigurations show that the re-entrant piston bowl assembly generates more turbulence. For both low and high intake swirl levels the re-entrant bowl produces aboutChapter 8. CONCLUSIONS AND RECOMMENDATIONS 13633% more turbulence at top dead centre than the bowl-in-piston arrangement (seeFigure 6.8).• The computations show that contrary to some experimental observations the turbulence length scaies do not necessarily have to scale with the top dead centreclearance height if piston geometries have bowls in them. The values of the lengthscales predicted at top dead centre position were consistent with the increasedvalues of turbulence intensity at the top dead centre position.• The results of the numerical simulation show that the re-entrant chamber assemblygenerated the highest level of turbulence diffusivity while the simple pancake combustion chamber configuration produced the least amount of turbulence diffusivity.This observation is consistent with the higher levels of squish motions and the resuiting vortex motion. These combinations help to diffuse turbulence much morefavourably. In particular, the diffusivity levels produced by the re-entrant bowlgeometry is about 10% higher than that produced by the bowl-in-piston chamberat low swirl levels and 50% more at high intake swirl levels (see Figure 6.10).• The computer simulations show that in general, higher intake swirl values lead tohigher flow and turbulence quantities.• The flow fields .generated by the KIVA-lI numerical code were very informativeand they were instrumental in visualizing trends which would have been impossibleotherwise. An example is the profound re-attachment phenomenon depicted by there-entrant-bowl assembly (see Figures 6.14a, 6.14b, 6.15a and 6.15b). An immediatebenefit that could be realized from this observation is to institute a program ofchamber optimization.Chapter 8. CONCLUSIONS AND RECOMMENDATIONS 137• After top dead centre, high shear and reverse-squish occurs in both chamber configurations, but these were more intense in the re-entrant bowl configuration. Thiscould be attributed directly to the re-entrancy features of the bowl. This flowcondition could enhance burnup of unburned fuel with a consequent reduction inhyrocarbon emissions.• In summary, we can conclude that the KIVA-Il numerical code is a versatile toolwhich can be cautiously employed in practical design problems to understand thevarious complex phenomena occurring in the internal combustion engine chamber.Some of the results also indicate that the KIVA-Il numerical code can be used tovisualize trends and if needed be used for chamber optimizations (see Figures 6.14a,6.14b, 6.15a and 6.15b).8.1.3 Data ComparisonSince one of the main objectives of the present study is to validate the KIVA-Il code,the predicted values were compared to the experimental data. Based on this comparison,the following conclusions were arrived at:• The computed and measured mean radial velocities were found to be generallywithin the experimental accuracy of each other, however, the measured values were5 to 20% higher in magnitude than the computed results.• The measured mean tangential velocities were found to be 10 to 25% higher thanthe predicted values.• The values of the turbulence intensity close to the bowl axis measured in the reentrant bowl chamber configuration were accurately predicted, but the quality ofprediction falls-off as the bowl entrance region is approached. The values of theChapter 8. CONCLUSIONS AND RECOMMENDATIONS 138turbulence intensity were however, poorly reproduced near the axis of the bowl-in-piston chamber assembly.• The values of turbulence in the three planes investigated were found to differ bybetween 10 to 15%.From the foregoing, it can be stated that, the KIVA-Il numerical code is capable ofaccurately reproducing the experimental data, provided the actual initial and boundaryconditions are precisely modelled.The results contained in chapter 7, where care was exercised in specifying the realinitial and boundary conditions, clearly produced data comparisons which were superiorto those obtained in actual engine studies Henriot et al (1989), zur Loye et al (1989). Inthese other studies the boundary conditions were not accurately specified, because theauthors found it very difficult to take measurements in the valve entrance area.The results of the present study also reveal that the fluid dynamics models embodied in the KIVA-Il code are sound, and any discrepancies between measurement andmodelling is most likely to be due to inaccuracies in specifying the correct boundaryconditions. This conclusion is clearly supported by the studies of Bizhan et al (1989).8.2 RecommendationsThe results of this study show that there are certain areas of in-cylinder flow field whichrequire further investigation. These areas are both experimental and numerical in nature.8.2.1 ExperimentalIn this study the effects of intake swirl and squish were considered, it would however,be necessary to investigate the effects of tumble on top dead centre turbulence intensity.Furthermore, a study of both types of in-cylinder bulk motion (swirl and tumble) isChapter 8. CONCLUSIONS AND RECOMMENDATIONS 139highly desirable with the aim of delineating their individual contributions to top deadcentre turbulence.It is also important to carry out length scale measurements in combustion chambersfeaturing swirl, tumble, squish or a combination of these. This will shed some lighton what length scales are most important in determining the structure of the flame.Moreover, the measurement of length scales will provide answers to the questions of theextent to which turbulence is homogeneous and isotropic during the period in whichcombustion occurs.Lastly, it is recommended that combustion studies be carried out in these complexchambers with swirl and tumble charge intake. This will further enhance our knowledgeof the various attributes of such combustion chambers cum intake flow configurations.8.2.2 NumericalThe various numerical results presented in this study has elucidated the amount of information that can be obtained from numerical modelling. However, some aspects ofnumerical modelling require further attention.The KIVA-Il code used in the present study does not have a model for tumble, andas such it is recommended that a model for tumble bulk motion be incorporated.It is also recommended that attempts be made to validate other submodels in theKIVA-Il code.8.2.3 Data ComparisonThe results in chapter 7 clearly highlight the importance of accurately specifying theinitial and boundary conditions in code initialization. In view of this, it is highly recommended that a lot of effort be expended in the accurate specification of the initial andboundary conditions.References[1] Abdel-Gayed, R.G., Al-Khishali, K.J. and Bradley, D., (1984) Turbulent BurningVelocities and Flame Straining in Explosions, Proc. R. Soc. London, vol. A391, pp.393-414.[2] Abdel-Gayed, R.G. and Bradley, D., (1981) A Two-Eddy Theory of PremixedTurbulent Flame Propagation, Phil. Trans. R. Soc. London, vol. A301, pp. 1-25.[3] Abraham, J., Wey, M.J. and Bracco, F.V., (1988) Pressure Non-Uniformity andMixing Characteristics in Stratified Charge Rotary Engine Combustion, SAE paper880624.[4] Abraham, J. and Bracco, F.V., (1988) Comparisons of Computed and MeasuredMean Velocity and Turbulence Intensity in a Motored Rotary Engine, SAE PaperNo. 881602.[5] Abraham, J., Bracco, F.V. and Reitz, R.D., (1985) Comparisons of Computed andMeasured Premixed Charge Engine Combustion, Combustion and Flame, vol. 60,pp. 309-322.[6] Ahmadi-Befrui Bizhan, Brandstätter Wilhelm, Kratochwill Heinrich and TrogerChristian, (1989) The Influence of Inlet Port Design on the In-Cylinder ChargeMixing, SAE 890842.[7] Amann, C.A., (1985) Classical Combustion Diagnostics for Engine Research, SAEpaper 850395.[8] Amato, U., Bertoli, C., Corcione, F.E., Valentino, G. and Vecchio, A., (1988)Experimental and Numerical Investigation of Air Flow Field in an Open ChamberDiesel Engine, SAE 880382.[9] Amsden, A.A., O’Rourke, P.J. and Butler, T.D., (1989) KIVA-lI: A Computer Program for Chemically Reactive Flows with Spray, Los Alamos National LaboratoryReport LA- 11560 - MS.[10] Amsden, A.A., Ramshaw, J.D., O’Rourke, P.J. and Dukowicz, J.K., (1985a) KIVA:A Computer Program for Two- and Three-Dimensional Fluid Flows with ChemicalReactions and Fuel Sprays, Los Alamos National Laboratory, Report LA-10245-MS.140References 141[11] Amsden, A.A., Ramshaw, .J.D., Cloutman, L.D. and O’Rourke, P.J., (1985b) Improvements and Extension to the KIVA Computer Program, Los Alamos NationalLaboratory report LA-16534-MS.[12] Andon, J. and Marks, C., (1963) Engine Roughness— The Key to Lower OctaneRequirement, SAE Transactions, vol. 72, pp. 636-658.[131 Andrews, G.E., Bradley, D. and Lwakabamba, S.B., (1975) Turbulence and Turbulent Flame Propagation — A Critical Appraisal, Combustion and Flame, 24, pp.285-304.[14] Arcoumanis, C., Ru, Z., Vafidis, C. and Whitelaw, J.H., (1990) Tumbling Motion:A Mechanism for Turbulence Enhancement in Spark-Ignition Engines, SAE 900060.[15] Arcoumanis, C., Hadjiapostolou, A. and Whitelaw, J.H., (1987) Swirl Center Precession in Engine Flows, SAE 870370.[16] Arcoumanis, C. and Whitelaw, J.H., (1987) Fluid Mechanics of Internal Combustion Engines — A Review, Proc. Inst. Mech. Engrs., vol. 201, No Cl.[17] Arcoumanis, C., Vafidis, C. and Whitelaw, J.H., (1985) Valve and In-CylinderFlow Generated by a Helical Port in a Production Diesel Engine, Proceedings ofSymposium on Flows in Reciprocating IC Engines, ASME-FED-28.[18] Arcoumanis, C., Bicen, A.F. and Whitelaw, J.H., (1983) Squish and Swirl-SquishInteraction in Motored Model Engines, J. Fluid Engineering, 105, pp. 105-112.[19] Arcoumanis, C., Bicen, A.F. and Whitelaw, J.R., (1982) Effects of Inlet Parameterson the Flow Characteristics in a Four-Stroke Model Engine, SAE 820750.[20] Bates, S.C., (1988) A Transparent Engine for Flow and Combustion VisualizationStudies, SAE Paper 880520.[21] Beacon, D.M., Remshaw, J. and Walker, K.L., (1981) the Use of In-Cylinder Modelling and LDA Techniques in Diesel Engine Design, in New Energy Conser-’uationTechnologies and Their Commercialization (eds. J.P. Milihone and E.H. Willis),Springer-Verlag, Berlin.[22] Bendat, J.S. and Piersol, A.G., (1986) Randon Data: Analysis and MeaswrementProcedure, Wiley-Interscience, New York.[23] Bizhan Ahmadi-Befrui, Wilhelm Brandstätter, Heinrich Kratochwill and ChristianTroger, (1989) The Influence of Inlet Port Design on the In-Cylinder Charge Mixing,SAE 890842.References 142[24] Blizzard, N.C. and Keck, J.C., (1974) Experimental and Theoretical Investigationof Turbulent Burning Model for Internal Combustion Engines, SAE Paper No.740191.[25] Bracco, F.V. and Grasso, F., (1983) Sensitivity of Chamber Turbulence to IntakeFlows in Axisymmetric Reciprocating Engines, AIAA J., vol. 21, pp. 607-640.[26] Bracco, F.V., (1974) Introducing a New Generation of More Detailed and Informative Combustion Models, SAE Paper No. 741174.[27] Brandi, F., Reverencic, I., Cartellieri, W. and Dent, J.C., (1979) Turbulent Air Flowin the Combustion Bowl of D.I. Diesel Engine and its Effect on Engine Performance,SAE 790040.[28] Bray, K.N.C., (1980) Turbulent Flows with Pre-Mixed Reactants, in TurbulentReacting Flows, P.A. Libby and F.A. Williams Eds., Springer.[291 Butler, T.D., Arnsden, A.A., O’Rourke, P.J. and Ramshaw, J.D., (1985) KIVA: AComprehensive Model for 2D and 3D Engine Simulations, SAE Paper No. 850554.[30] Butler, T.D., Cloutman, L.D., Dukowicz, J.K. and Ramshaw, J.D., (1981) Multidimensional Numerical Simulation of Reactive Flow In Internal Combustion Engines,Prog. Energy Combustion Science, vol. 7, pp. 293-315.[31] Butler, T.D., Cloutman, L.D., Dukowicz, J.K. and Ramshaw, J.D., (1979) CONCHAS: an Arbitrary Lagrangian-Eulerian Computer Code for Multi-ComponentChemically Reactive Fluid Flow at All Speeds, Los Alamos Scientific LaboratoriesReport No. LA-8129-MS.[32] Clerk, D., (1921) Cylinder Actions in Gas and Gasoline Engines, SAE Journal, vol.8, pp. 523-539.[33] Corcione, F.E and Valentino, G., (1991) Characterization of Fluid-Dynamic Behavior of Diesel Engines by LDA Technique, Laser Anemometry, vol. 1, ASME.[34] Davis, G.C. and Borgnakke, C., (1982) The Effects of In-Cylinder Flow Processes(Swirl, Squish and Turbulence Intensity) on Engine Efficiency - Model Predictions,SAE 820045.[35] Davis, G.C. and Kent, J.C., (1979) Comparison of Model Calculations and Experimental Measurements of the Bulk Cylinder Flow Processes in a Motored PROCOEngine, SAE 790290.[36] Dent, J.C. and Salama, N.S., (1975) The Measurement of Turbulence Characteristics in an Internal Combustion Engine Cylinder, SAE Paper No. 750886.References 143[37] Dohring, K., (1986) The Relative Effects of Intake and Compression Stroke Generated Turbulence on I.C. Engine Combustion Duration, M.A.Sc. Thesis, Universityof British Columbia, Report AFL-86-01.[38] El-Tahry, S.H., (1983) k - e Equation for Compressible Reciprocating Engine Flows,Journal of Energy, 7, pp. 345-353.[39] El-Tahry, S.H., (1982) A Numerical Study on the Effects of Fluids Motion at Inlet-Valve Closure on Subsequent Fluid Motion in a Motored Engine, SAE 820035.[40] Errera, M.P., Labbe, J. and Jerot, A., (1988) Three-Dimensional Numerical andExperimental Analysis of In-Cylinder Flow in an Internal Combustion Engine, SAEPaper No. 880106.[41] Fansler, T.D. and French, D.T., (1988) Cycle-Resolved Laser-Velocimetry Measurements in a Re-entrant- Bowl-in-Piston Engine, Trans. SAE, Paper No. 880377.[42] Fukumoto, A., Ohsawa, K. and Ohkubo, Y., (1985) An Assessment of a Multidimensional Numerical Method to Predict the Flow in Internal Combustion Engines,SAE Paper No. 850500.[43] Glover, A.R., Hundleby, G.E. and Hadded, 0., (1988) The Developement of Scanning LDA for the Measurement of Turbulence in Engines, SAE Paper 880378.[44] Glover, A.R., (1986) Towards Bias-Free Estimates of Turbulence in Engines, InThird Intl Symp. on Applications of Laser Anemometry to Fluid Mechanics, Lisbon,July 7-9.[45] Gosman, A.D., (1985a) Multidimensional Modelling of Cold Flows and Turbulencein Reciprocating Engines, SAE Paper No. 850344.[46] Gosman, A.D., (1985b) Computer Modeling of Flow and Heat Transfer in Engines, Progress and Prospects, In Proc. Intl. Symp. on Diagnotics and Modeling ofCombustion in Reciprocationg Engines, Tokyo, (September).[47] Gosman, A.D., Tsui, Y.Y. and Vafidis, C., (1985) Flow in a Model Engine witha Shrouded Valve- A Combined Experimental and Computational Study, SAE850498.• [48] Gosman, A.D., Tsui, Y.Y. and Watkins, A.P., (1984) Calculation of ThreeDimensional Air Motion in Model Engines, SAE Paper No. 840229.[49] Gosman, A.D. and Harvey, P.S., (1982) Computer Analysis of Fuel-Air Mixing andCombustion in an Axisymmetric Direct-Injection Diesel, SAE Paper No. 820036,vol. 89.References 144[50] Gosman, A.D., Johns, R.J.R. and Watkins, A.P., (1980) Development of Prediction Methods for In-Cylinder Processes in Reciprocating Engines, in CombustionModelling in Reciprocating Engines, J.N. Mattavi and C.A. Amann, Eds., PlenumPress, pp. 69-124.[51] Gosman, A.D. and Johns, R.J.R., (1978) Development of a Predictive Tool forIn-Cylinder Gas Motion in Engines, SAE 780315.[52] Groff, E.G. and Matekunas, F.A., (1980) The Nature of Turbulent Flame Propagation in a Homogeneous Spark-Ignited Engine, SAE Trans., vol. 89, pp. 740-764,Paper 800133.[53] Hadded, 0. and Denbratt, I., (1991) Turbulence Characteristics of Tumbling AirMotion in Four-Valve S.I. Engines and their Correlation with Combustion Parameters, SAE 910478.[54] Harworth, D.C. and El Tahry, S.H., (1991) Probability Density Function Approachfor Multidimensional Turbulent Flow Calculations with Application to In-CylinderFlows in Reciprocating Engines, AIAA Journal, vol. 29, No. 2 pp. 208-218.[55] Harworth, D.C., E1-Tahry, S.H., Huebler, M.S. and Chang, S., (1990) Multidimensional Port-Port-and-Cylinder Flow Calculations for Two- and Four-Valve-Per-Cylinder Engines: Influence of Intake Configuration on Flow Structure, SAE No.900257.[56] Hayder, M.E., Varna, A.K. and Bracco, F.V., (1985) A Limit to TDC TurbulenceIntensity in Internal Combustion Engines, AIAA J., Propul. Powr., 1(4), pp. 300-308.[57] Henriot, S., Le Coz, J.F. and Pinchon, P., (1989) Three Dimensional Modelling ofFlow and Turbulence in a Four-Valve Spark Ignition Engine — Comparison withLDV Measurements, SAE Paper No. 890843.[58] Heywood, John, B., (1988) Internal Combustion Engine Fundamentals, McGraw-Hill, Inc., p. 403.[59] Heywood, J.B., (1987) Flow Motion Within the Cylinder of Internal CombustionEngines — The 1986 Freeman Scholar Lecture. Journal of Fluid Engineering, vol.109.[60] Heywood, J.B., (1980) Engine Combustion Modelling— An Overview, in Cornbustion Modelling in Reciprocating Engines, Eds. J.N. Mattavi and C.A. Amann,Plenum Press, New York, 1980, pp. 1-38.References. 145[61] Hirt, C.W., Amsden, A.A. and Cook, J.L., (1974) An Arbitrary LagrangianEulerian Computing Method for All Flow Speeds, J. Comp. Phys., vol. 14, pp.227-253.[62] Hirt, C.W., Cook, J.L. and Butler, T.D., (1970) A Lagrangian Method For Calculating the Dynamics of an Incompressible Fluid with Free Surface, J. Comp. Phys.5, p. 103.[63] lijima, T. and Bracco, F.V., (1987) LDV Measurements in an Engine with Squareand Circular Piston Cups, SAE paper 872073.[64] Ikegami, M., Shioji, M. and Mishimoto, K., (1987) Turbulence Intensity and Spatial Integral Scale During Compression and Expansion Strokes in a Four-CycleReciprocating Engine, SAE Paper 870372.[65] Johns, R.J.R., (1984) A Unified Method for Calculation Engine Flows, ASMEPaper No. 84-DGP-18.[66] Jones, W.P., (1979) Models for Turbulent Flows with Variable Density and Combustion, Prediction Methods for Turbulent Flows, Von-Krmán Institute for FluidDynamics, Lecture Series 1979-2.[67] Kajiyarna, K., Nishida, K., Murakami, A., Arai, M. and Hiroyasu, H., (1984) AnAnalysis of Swirling Flow in Cylinder for Predicting D.I. Diesel Engine Performance, SAE 84b518.[68] Karamcheti, K., (1980) Principles of Ideal-Fluid Aero-dymamics, Robert E. KriegerPublishing Company, Malabar, Florida.[69] Kent, J.C., Mikulec, A., Rimai, L., Adamczyk, A.A., Mueller, S.R., Stein, R.A.and Warren, C.C., (1989) Observations on the Effects of Intake-Generated Swirland Tumble on Combustion Duration, SAE 892096.[70] Kent, J.C., Haghgooie, M., Mikulec, A., Davis, G.C. and Tabaczynski, R.J., (1987)Effects of Intake Port Design and Valve Lift on In-Cylinder Flow and Burnrate,SAE Paper No. 872153.[71] Lancaster, D.R., (1976) Effects of Engine Variables on Turbulence in a SparkIgnition Engine, SAE paper 760159.[72] Lancaster, D.R., Krieger, S., Sorenson C. and Hull, W.L., (1976) Effects of Turbulence on Spark-Ignition Engine Combustion, SAE Trans., vol. 85, pp. 689-764,Paper 760160.[73] Laser Anemometer Systems, Thermo-Systems Inc., St. Paul, Minnesota, 1976.References 146[74] Laser Velocimeter Interface - model 1400A, Owner’s manual, DosTek Inc., Kitchener, Ontario, Canada, 1989-1990.[75] Launder, B.E. and Spalding, D.B., (1972) Mathematical Models of Turbulence,Academic Press, New York.[76] Le Coz, J.F., Henriot, S. and Pinchon, P., (1990) An Experimental and Computational Analysis of the Flow Field in a Four-Valve-Spark Ignition Engine — Focuson Cycle-Resolved Turbulence, SAE Paper No. 90056.[77] Lee, D.W., (1939) A Study of Air Flow in an Engine Cylinder, NACA Report No.653.[78] Li, Z., Steinthorsson, E., Shih, T. I-P. and Nguyen, H.L., (1990) Modelling andSimulation of Wankel Engine Flow Fields, SAE Paper No. 900029.[79] Liou, T.M. and Santavicca, D.A., (1983) Cycle-Resolved Turbulence Measurementsin a. Ported Engine With and Without Swirl, SAE Paper 830419, (also SAE Transactions vol. 92, sec. 2, pp. 131-145).[80] Magnussen, B.F. and Hjertager, B., (1976) On the Mathematical Modelling ofTurbulent Combustion with Special Emphasis on Soot Formation and Combustion,Proc. 16th International Symposium on Combustion, pp. 719-727, The CombustionInstitute, Pittsburg.[81] Marko, K.A., Li, P., Rimai, L., Ma, T. and Davies, M., (1986) Flow Field Imagingfor Quantitative Cycle-Resolved Velocity Measurements in a Model Engine, SAEPaper 860022.[82] Matekunas, F.A., (1983) Modes and Measures of Cyclic Combustion Variability,SAE Paper No.. 830337.[83] Matsuoka, S., Nagakura, K., Kawai, T., Kumimoto, T. and Aoyagi, Y., (1980) Application of Laser Doppler Anemometry to a Motored Diesel Engine, SAE 800965.[84] Mattavi, J.N., Groff, E.G., Lienesch, J.H., Matekunas, F.A. and Noyes, R.N.,(1980) Engine Improvements Through Combustion Modelling, Combustion Modelling in Reciprocating Engines, Eds. J.N. Mattavi and C.A. Amann, Plenum Press,New York, pp. 537-588.[85] Mawle Craig, D., (1989) The Effects of Turbulence and Combustion Chamber Geometry on Combustion in a Spark Ignition Engine, M.A.Sc. Thesis, University ofBritish Columbia, Report AFL-89-02.References 147[86] McCuiston, F.D., Lavoie, G.A. and Kauffman, C.W., (1977) Validation of a Turbulent Flame Propagation Model for a Spark Ignition Engine, SAE Paper 770045.[87] McKinley, T.L., Primus, R.J., O’Rourke, P.J. and Butler, T.D., (1988) Comparisonof Computed and Measured Air Motion in Circular and Square Piston Cups, SAE881612.[88] Mikulec, A., Kent, J.C. and Tabaczynski, R. J., (1988) The Effect of Swirl onCombustion in a Pancake Chamber Spark Ignition Engine: The Case of ConstantInducted Kinetic Energy, Trans. SAE, Paper No. 880200.[89] Monaghan, M.L. and Pettifer, H.F., (1981) Air Motion and Its Effect on DieselPerformance and Emissions, SAE 810255.[90] Morse, A.P., Whitelaw, J.W. and Yianneskis, M., (1980) The Flow Characteristicsof a Piston-Cylinder Assembly with an Off-Centre, Open Port, Proc. Inst. Mech.Engrs. vol. 194, No. 31.[91] Moss, J.B., (1980) Simultaneous Measurements of Concentration and Velocity inan Open Premixed Turbulent Flame, Combustion Science and Technology, vol. 22,pp. 115-129.[92] Nagayama, I., Araki, Y. and lioka, Y., (1977) Effects of Swirl and Squish on SIEngine Combustion and Emission, SAE 770217.[93] Naitoh, K., Fujii, H., Urushihara, T., Takagi, Y. and Kuwahara, K., (1990) Numerical Simulation of the Detailed Flow in Engine Ports and Cylinders, SAE PaperNo. 900256.[94] Patankar, S.V., (1980) Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, Washington, D.C.[95] Rarnos, G.I. and Sirignano, W.A., (1980) Axisymmetric Flow Model With andWithout Swirl in a Piston-Cylinder Arrangment with Idealized Valve Operation,SAE Paper No. 800284.[96] Reynolds, W.C., (1980) Modelling of Fluid Motions in Engines — An IntroductoryOverview, in Combustion Modelling in Reciprocating Engines, J.N. Mattavi andC.A. Amann, Eds., Plenum Press, pp. 69 - 124.[97] Reynolds, W.C., (1976) Computation of Turbulent Flows, Ann. Rev. Fluid Mech.,vol. 8, pp. 1183-208.[98] Roache, P.J., (1982) Computational Fluid Dynamics, Revised Edition, HermosaPublishers, Albuquerque, New Mexico, pp. 64-67.References 148[99] Roberts, J.B. and Ajmani D.B.S., (1986) Spectral Analysis of Randomly SampledSignals Using a Correlation-Based Slotting Technique, TEE Proceeding, 133F: pp.153-162.[100] Semenov, E.S., (1958) Studies of Turbulent Gas Flow in Piston Engines, OtdelenieTekhnicheskikh Nauk, No. 8, 1958 (NASA Technical Translation F-97, 1963).[101] Stull, D.R. and Prophet, H., (1974) JANAF Thermochemical Tables, 2nd ed. (U.S.Department of Commerce/National Bureau of Standards, NSRDS-NBS 37, June1971). N.W. Chase et al., J. Phys. Chem. Ref. Data 3, p. 311.[102] Sung, N.W. and Patterson, D.J., (1982) Air motion in a Two Stroke Engine Cylinder - The Effects of Exhaust Geometry, SAE 820751.[103] Tabaczynski, R.J., (1976) Turbulence and Turbulent Combustion in Spark-IgnitionEngines, Frog. Energy Combustion Science, vol. 2, pp. 143-165.[104] Taylor, G.I., (1938) The Spectrum of Turbulence, Proceeding Royal Society, London, 164A.[105] Tennekes, H. and Luniley, J.L., (1972) A First Course in Turbulence, M.I.T. Press.[106] Tennekes, H., (1968) Simple Model for the Small-Scale Structure of Turbulence,Plays. Fluids, vol. 11, pp. 669-671.[107] Tindall, M.J., Williams, T.J. and El Khafaji, A.H.A., (1974) Gas Flow Measurements in Engine Cylinders, SAE paper 740719.[108] Traci, R.M., (1981) MAGIC CODE in Comparisons Between Measurement andAnalysis of Fluid Motion in Internal Combustion Engines, P.O. Witze, Ed., SandiaLaboratories Report SAND81-8242.[109] Uzkan Teoman and Hazelton, J.R., (1986) The Influence of Swirl on the FreshCharge Stratification in an IC Engine Combustion Chamber, SAE 860466.[110] Vafidis, C., (1984) Influence of Induction Swirl and Piston Configuration on AirFlow in a Four-Stroke Model Engine, Proc. I. Mech. E. 198C, pp. 71-79.[111] Watkins, A.P., (1977) Flow and Heat Transfer in Piston/Cylinder Assemblies,Ph.D. Thesis, University of London.[112] Witze, P.O., Martin, J.K. and Borgnakke, C., (1984) Conditionally-Sampled Velocity and Turbulence Measurements in a Spark Ignition Engine, Combustion Scienceand Technology, 36, pp. 301-317.References 149[113] Witze, P.O., (1982) The Effects of Spark Location on Combustion in a Variable-Swirl Engine, SAE Paper 820044.[114] Witze, P.O. and Vilchis, F.R., (1981) Stroboscopic Laser Shadowgraph Study ofthe Effect of Swirl on Homogeneous Combustion in a Spark-Ignition Engine, SAE810226.[115] zur Loye, A.O., Siebers, D.L., McKinley, T.L., Ng, H.K. and Primus, R.J., (1989)Cycle-Resolved LDV Measurements in a Motored Diesel Engine and Comparisonwith k-e Model.Predictions, SAE 890618.Tables 150Table 3.1: Rapid Intake and Compression Machine Operating ParametersBore 101.6 mmStroke 100.0 mmCompression Ratio 8:1Initial Pressure (1-Stroke) 50.8 kPaInitial Pessure (2-Stroke) AmbientBowl-in-Piston Squish Area 80%Re-entrant Piston Squish Area 80%Intake Port Opens 24 deg BTDCPort Closes 60 deg BTDCDriving Cylinder Pressure 900 kPaSimulated Engine Revolution 1000 rpmTables 151Table 4.1: Standard k -€ Turbulence Model Constantsc Cf2 C PTk P?’f1.44 1.92 -1.0 1.0 1.3Figures 152Figure 2.1: Clerk’s Experiment To Demonstrate Role of Turbulence, Clerk (1921).27- Z’////////////½Z //////////////////// /;Figure 2.2: Ricardo’s Chamber For Producing Squish Flow (see Amann (1985)).Figures 1536.05.04.0Iww3.0w2.01.01.0 2.0 3.0 4.0 5.0 6.0TURBULENCE INTENSITY (mis)Figure 2.3: Correlation of Flame Speed ratio with Turbulence Intensity, Lancaster (1976),Lancaster et a! (1976).I I I I ——B- B SHROUDEDI VALVEB---NON-SHROUDEDVALVEz A SPEEDQ • VOLUMETRIC EFFICIENCYE I COMPRESSION RATIO+ EQUIVALENCE RATIO220cmFigure3.1:SchematicDiagramofUBCRICM,Dohring(1986).Oq•HYDRAUUCBRAKINGCYUNDERHYDRAULICOILTANKANDASSOCIATEDGEARSTOPVIEWDRIVINGRACKPNEUMATICDRIVINGCYLINDERE 0 0BEDPLATERAPIDCOMPRESSIONMACHINECOMPRESSEDAIRC,’FiguresAIR INTAKEAIR INTAKE155EXHAUSTFigure 3.2: Schematic Diagram of Intake Arrangement.Figuresr///BORERe-Entrant PistonFigure 3.3: Idealized Combustion Chamber Flow Configuration. Dimensions are in millimetres.156SWIRLSQUISHt0/ TzFigure3.4:BlockDiagramofmodel 1400ALaserVelocimeterInterface.Oq’1OPTICALSHAFTENCODER 8-BITINPUTAI.’a__________REGISTERREGISTERREGISTER[ISTERIPROGRAMMABLEINTERVALTIMER486PERSONALCOMPUTERBUS01ImøQoc,sarn.rqFigures‘C‘C159Figure 3.6: Schematic of Jet created by Flow through the Intake Valve indicating itsTurbulent Structure, Reynolds (1980).Lr——————Figures 16083 .‘._2(xUk, Yuk zak) 1aFigure 4.1: Typical Finite-difference Cell./“E/(x&k, Yijk’ Zijk)Figure 4.2: The portion of momentum cell (i, k) lying within regular cell (, j, k). Thethree momentum cell faces lying within the regular cell are shaded.Figures 1618‘I5( jFRONT/4 %%% J/ 73/ BOOM /2Figure 4.3: Cell Faces Associated with cell (i, j, k).XtXd.. ‘EL.XeXrXbFigure 4.4: The six points used to define the gradient of cell-centred quantity Q on cellface a.Figures 162fyFigure 4.5: Cell shapes before (dashed lines) and after (solid lines) the Lagrangian phasein a plane parallel flow.Re-EntrantPistonBowl-In-PistonPlaneAPlaneBPlaneC‘1Figure5.1:Measurement CoordinateSystemRelativeToThePistonatTDC.Dimensionsinmillimetres.I—.Figures 164100.058.0.4— Quartz Glass Window.— Combustion Chamber Wall— Bowl-In-Piston100.044.0Quartz Glass Window. Combustion ChamberWall_____Re-entrant Bowl PistonFigure 5.2: Bowl-in-Piston and Re-entrant bowl pistons. Dimensions are in millimetres.Crosses indicate LDV measurement locations.Figures 1657.5r — 1.0cm Swirl Ratio —0.5 (a)30 mm0.CoC5.0 Re-entrant Bowl Chamber- - -- Bowl-In-Piston Chamber / NI+ - 10.0>+ Tz—3.02.5 - .._—_—a•a,0.0 I,.,,’60 50 40 30 20 10 0Crank Angle (degree BTDC)7.5r—1.Ocm SwirlRatio—2.0 (b)z— 3.0 mm0Ci) —-...-C..._• .—.—o_______________________I50 I_._.._ Re-entrant Bowl Chamber /I - - - - Bowl-In-Piston Chamber1 / \(a\ Q++ r—io.o2.5 +1Tz—3.o—____________a, -.---a--0.060 50 40 30 20 10 0Crank Angle (degree BTDC)Figure 5.3: Evolutioii of Measured Mean Radial Velocity with Crank Angle. CrossesIndicate Measurement Location in millimetres.Figures 1667.5.r—2.Ocm SwirlRatio—0.5 (a)6.0 mm0.C’)5.0 Re-entrant Bowl Chamber S,- - - - Bowl-In-Piston Chamber 7I.-+ r—20.0—— ‘‘a \2.5 __a- + T = 6.0— _0._... (0.0 .1 I,,,, I....60 50 40 30 20 10Crank Angle (degree BTDC)7.5r —2.0 ii Swirl Ratio —2.0 (b)z6.0mmCl)5.0__••__ Re-entrant Bowl Chamber 7 -- - -- Bowl-In-Piston Chamber /- -- \I __;‘‘2.5- I-.-+ -6.000.060 50 40 30 20 10 0Crank Angle (degree BTDC)Figure 5.4: Evolution of Measured Mean Radial Velocity with Crank Angle. CrossesIndicate Measurement Location in millimetres.Figures 167r 1.0 ii Swirl Ratio —0.5 (a)z— 3.0 mmU,C05.0-__••_ Re-entrant Bowl Chamber /- - - - Bowl-In-Piston Chamber /I_./-;::::- - --Os c... r - 10.0—F. — - F ÷ T z — 3.0E2.5-...-.—-... ..W- -.-I-0.060 50 40 30 20 10 0Crank Angle (degree BTDC)15.0r—1.Ocm SwirlRatio—2.0 (bz— 3.0 mm0.U,C12.5 --.-.o.._•____________________/_._.,._._ Re-entrant Bowl Chamber /--- Bowl-dn-Piston Chamber /ioo 0 r—100—-•— T + Tz—3.0C 7•5 -________5.0 .1...,60 50 40 30 20 10 0Crank Angle (degree BTDC)Figure 5.5: Variation of Measured Mean Tangential Velocity with Crank Angle. CrossesIndicate Measurement Location in millimetres.Figures 168r —2.0 ii Swul Ratio —0.5 (a-g z—6.OmmaCl)Co I— /_______________________I.S,O 5.0-_•_ Re-entrant Bowl Chamber.1__—j--_--- Bowl-In-Piston Chamber- + r - 20.0-S2.5- D - - a- - - + Z 6.0aCa,0.0 .1,,.. I.,.. I.... I....60 50 40 30 20 10 0Crank Angle (degree BTDC)15.0- Ir —2.0 Swirl Ratio —2.0z— 6.0 mma,C.C/)12.5 - /.......... Re-entrant Bowl Chamber If--- Bowl-In-Piston Chamber /- -,.- r-20.0—-+ T Z 6.0-________--5.060 50 40 30 20 10 0Crank Angle (degree BTDC)Figure 5.6: Variation of Measured Mean Tangential Velocity with Crank Angle. CrossesIndicate Measurement Location in millimetres.Figures 16914.0r — 2.2 ti Re-entrant Bowl Chamber (a)z —6.0 mm Swirl Ratio — -——- -- Tangential10.0• RadialCrank Angle (degree BTDC)14.0r — 2.9 fl Bowl-In-Piston Chamber (b)z = 6.0mm Swirl Ratio —2.0.,12.0Tangential10.0 -• Radial8.0 + r—29.O[ Tz-6OCrank Angle (degree BTDC)Figure 5.7: Mean Radial and Tangential Velocity Results Measured Close to Bowl Entrance Region. Crosses Indicate Measurement Location in millimetres.Figures 17012r — 0.0cm Re-entrant Bowl Chamber (a)10 z 6.0 mm Swirl Ratio 2.0-Tangential VelocityC.’) 8• Radial Velocity Fq (( r-O.OT—— 4. z—6.04 .D• ——___________—..—,.‘ _u II ()Crank Angle (degree BTDC)14 r—0.Ocm Bowl-In-PlstonChamber (b)z=6.Omm SwirIRatio2.012- - -- Tangential Velocityu, 10 • Radial Velocityq: — — S — — — r — 0.0— -,‘ ‘S..f I z —6.0Crank Angle (degree BTDC)Figure 5.8: Mean Radial and Tangential Velocity Results Measured on the Cylinder Axis.Crosses Indicate Measurement Location in millimetres.Figures 171r—1.Ocm SwirlRatlo-0.5 (a)0.75 z — 3.0 mm/1Re-entrant Bowl Chamber /—-- Bowl-in-Piston Chamber A - -- -- (( r -10.0—o.5o •‘.. _._._w\ -::::::::;1-- Tz-3.0I0.25O• 40 20 •iOCrank Angle (degree BTDC)1.25r—1.Ocm SwirlRatio—2.O (b)z—3.Omm0.CO /C /Re-entrant Bowl Chamber /1.00 --- Bowl-In-Piston Chamber / -\\\\\ ,v r-i0.o0.75 •—._. . —...—.—.—...—.—- + Tz-3.00.50.O 40 20 0Crank Angle (degree BTDC)Figure 5.9: Evolution of Measured Radial Turbulence Intensity with Crank Angle.Crosses Indicate Measurement Location in millimetres.Figures 1721.00r—2.Ocrn SwirlRatio—0.5 (a)z— 9.0 mm__._._ Re-entrant Bowl Chamber0.75- -- Bowl-In-Piston Chamber-W -+ r - 20.0J - U - - . .. .. L0+1TZ-90Crank Angle (degree BTDC)1.10—2.0cm Swirl Ratio —2.0 (b)z—9.Ommai.ooU)Re-entrant Bowi Chamber-- -- Bowl-In-Piston Chamber /0.90 .:-“ ;z:0.60 ‘I60 50 40 30 20 10 0Crank Angle (degree BTDC)Figure 5.10: Evolution of Measured Radial Turbulence Intensity with Crank Angle.Crosses Indicate Measurement Location in millimetres.Figures 1730.80r — 1.0cm Swirl Ratio 0.5 (a)z— 3.0 mmcO.70.. Re-entrant Bowl Chamber-- -- Bowl-In-Piston Chamber 70.60- - -r - 10.0d + Tz -3.0--a- - - - - - D1• 0.40tI-0.30 I’60 50 40 30 20 10 0Crank Angle (degree BTDC)1.00r — 1.0cm Swirl Ratio —2.0 (b)CO z—3.Omm /C I0.90a-_. Re-entrant Bowl Chamber A-- -- Bowl-In-Piston Chamber _10.80 ç4. r—1O.09-S.-— I—— -a. -c 0.70._________+ Tz-3.0I—0.60•C0.5040 o•• io••Crank Angle (degree BTDC)Figure 5.11: Evolution of Measured Tangential Turbulence Intensity with Crank Angle.Crosses Indicate Measurement Location in millimetres.Figures 1740.80r 2.0cm Swirl Ratio —0.5 (a)z—9.Omm0.70Re-entrant Bowl ChamberI - - - - Bowl-In-Piston Chamber II—° i J. r — -9.0!___0.4OCCI-0.30 I I I60 50 40 30 20 10 0Crank Angle (degree BTDC)1.00r — 2.0cm Swirl Ratio —2.0 (b)0 z—9.OmmF—/..._._ Re-entrant Bowl Chamber- - -- Bowl-In-Piston Chamber.-•080- - - - —- -: Q r -2000.70 - - -__ ______:5040 20 . -o+i1z_9.0Crank Angle (degree BTDC)Figure 5.12: Evolution of Measured Tangential Turbulence Intensity with Crank Angle.Crosses Indicate Measurement Location in millimetres.Figures 1751.20r —2.2 cm Re-entrant Bowl Chamber (a)z —6.0mm Swirl Ratio = 2.01.00I:::<::. ——+ T z—6.00.20o.oo I I60 50 40 30 20 10 0Crank Angle (degree BTDC)1.20r —2.9 cm BOwlinPiStOn Chamber (b)z —6.0 mm Swirl Ratio —2.0:.1.oo_______________C+0.40..______0200.00 I I I60 50 40 30 20 10 0Crank Angle (degree BTDC)Figure 5.13: Radial and Tangential Turbulence Intensity Results Measured Close to BowlEntrance Region. Crosses Indicate Measurement Location in millimetres.Figures 1761.40r — 0.0 ,i Me-entrant Bowl Chamber (a)z — 6.0mm Swirl Ratio ——- Tangential0• Radial01.00.j. r-O.O6o:3o20:oTz6OCrank Angie (degree BTDC)1.40r = 0.0 Bowl-In-Piston Chamber (b)z = 6.0mm Swirl Ratio — - - - Tangential1.00• Radial10.80j0.60 \.S S0.4060 50 40 30 20 10 0Crank Angle (degree BTDC)Figure 5.14: Radial and Tangential Turbulence Intensity Results Measured on the Cylinder Axis. Crosses Indicate Measurement Location in millimetres.oq I-. -S. I‘P CD ‘1 C,. 0 00 C,.0Oq 1/::::1:I:jI-’.-Joq I-.C.,I-.00oq Ct’ CA)0 0 CD CD 9 90qoq I-i CD 0 0 0 CD CD 9 CACD I—.00 QFigures 1817,5(a) Swirl Ratio = 0.5a.Cl)o i %5.0Re-entrant Bawl Chai1 / \-- Bowl-In-Piston Chamber f II!73I s— I s0—. .. ..> ,1’ I.‘ I•..-W’_,_D0 _.._._._._•_••••_..O’_.-o-_-- a - - - .-v.050 40 30 20 10 0Crank Angle (degree BTDC)7.5(b) Swirl Ratio -2.00. I’CoI \___________I sg 5.0,,.,. Re-entrant Bowl Chambé9 / \--- Bowl-In-Piston Chamber // \73 I’o I, IF’ I> / ‘I2.5—. ,..a—— .D. . — ,/— , ——C P——D...—C•‘I0.050 40 30 20 10 0Crank Angle (degree BTDC)Figure 6.5: Effect of Bowl Shape on Radial Velocity.Figuresa,aCo5.0a,0C2.5ICa,15.00a,a 12.5Cl)C0110.07.5a,>5.0I2.50.050 40 30 20 10 0Crank Angle (degree BTDC)Figure 6.6: Effect of Bowl Shape on Tangential Velocity.18250 40 30 20 10 0Crank Angle (degree BTDC)Re-entrant Bowl Chamber--- Bowl-In-Piston Chamber(b) swirl Ratio— 2.0—. •_..——Figures0.Cl)C010500.25Cl)0.75183Re-entrant Bowl Chamber-- Bowl-In-Piston Chamber(a) Swirl Ratio —0.5.1/ ••%F‘I‘I . S$1— -°- — / n_.-._._. p... .—— _-a..--50 40 30 20Crank Angle (degree BTDC)10 00aCl)II1. 40 30 20 10 0Crank Angle (degree BTDC)Figure 6.7: Effect of Bowl Shape on Intake Swirl Decay Rate.Figuresa,aC01.00.50.0-—- Bowl-In-Piston ChamberRe-entrant Bowl Chamber184(a) Swirl Ratio—0.5Re-entrant Bowl Chamber-- -- Bowl-In-Piston Chamber —-II/I50 40 20 10Crank Angle (degree BTDC)Swirl Ratio —2.0(b)•1-0- 0/ ___0—_ -..0-——- -—- --I,‘II,‘II2.5a,0.0C0 2.0aCa,15a,Ca) 40 30 20 10 0Crank Angle (degree BTDC)Figure 6.8: Effect of Bowl Shape on Turbulence Intensity.Figures 185I ._._. Swirl Ratio = 0.5I-—.-- SwirlRatio=2.0ISwirl Ratio= 3.010.25 0.50 0.75Non-dimensional Radial Distance2.5Plane A Re-entrant Bowl ChamberCrank Angle = 5 degree BTDC2.0 - (a)c 1.5 -_____________________— — —.- — —- -—C.-.1.0a(I)0(I)CU,CU,CU,C)CU,.01.1 00.00.00 0.25 0.50 0.75Non-dimensional Radial DistanceEffect of Intake Swirl on Turbulence Intensity.1.00Figure 6.9:Figures 1863.02.50.Cl)2.0_______________1.511. 0.25 0.50 0.75 1.00Non-dimensional Radial DistanceaU)IPlane B Re-enlrani Bowl ChamberCrank Angle = 5 degree BTDC(a)__. Swirl Ratio = 0.5- -. - Swirl Ratio = 2.0-._o._._ Swirl Ratio = 3.0-————- -0.00 0.25 0.50 0.75 1.00Non-dimensIonal Radial DistanceFigure 6.10: Effect of Intake Swirl on Turbulence Intensity.Figures 1873.5c 2.50cj0Cu 2.00U)C0C8c 1.00z.00.5fl-n1.5U)0o 1.0C00C00.52Plane C Re-entrant Bowl ChamberCrank Angle = 5 degree BTDC(a)-___._. Swirl Ratio = 0.5- - -- Swirl Ratio = 2.0.o—_._._._Swirl Ratio = 3.0•__.o.__-0.00 0.25 0.50 0.75 1.00Non-dimensional Radial DistancePlane C(b)Bowl-In-Piston ChamberCrank Angie =5 degree BTDC• Swirl Ratio = 0.51I__.__ SwirlRatio=2.0_.__._ Swirl Ratio= 3.01.1dø0.00.00 0.25 0.50 0.75Non-dimensional Radial Distance1.00Figure 6.11: Effect of Intake Swirl on Turbulence Intensity.Figures2.5a,2.01.5U)C1.0a,C)Ca,e 0.50.0a,Z 5.0U)2.53Ca,I—0.018850 40 30 20 10 0Crank Angle (degree BTDC)50 40 30 20 10 0Crank Angle (degree BTDC)Figure 6.12: Effect of Bowl Shape on Turbulence Length Scale.50 40 30 20 10 0Crank Angle (degree BTDC)s.... Disc Chamber-- Bowl-In-Piston ChamberRe-entrant Bowl Chamber50 40 30 20Crank Angle (degree BTDC)Figures 189(a)Swirl Ratio -0.5..Disc Chamber-- Bowl-In-Piston ChamberRe-entrant Bowl Chamber0.Oh. - --- -.-0I .1... .1 02.00a,x8a5° 1.5011.000.508Ca,=.D2o.0o7.5a,z8Ca,(aa,?.5.00.Cl)ga,0a,2.5- ,.. -8.*-- 4---4-- -0- - -.-. - -..fin- ---I-..(b)Swirl Ratio = 2.010 0Figure 6.13: Effect of Bowl Shape on Turbulence Diffusivity.(a)HighSwirl=2.0(b)LowSwirl=0.5___________________________(((((((‘(____________________•;•jIliiiIill,__,/fI!z:?lilpYi1j______,/1/‘.______-velocitycontourmaxspeed=1.991x—-:::‘%.•..I\\\\\\\\\%%__.__Figure6.14:VelocityVector(cm/3)Plotacrossj=1Planeat25degreeBTDC.Ii (c)Hi:ea=2.9xios]I.;;j!4f(d)(a)HighSwirl=2.0(b)LowSwirl=0.5I4III//_::::;1Iiiu)___Figure6.15:VelocityVector(cm/s)Plotacrossj=1Planeat20degreeBTDC.(c)H1;I!velocitycontourI imaxspeed231x103UIHflHI(d)‘%s’._,,/IJJ//!i,,,,/111velocitycontour___maxspeed=2.0834x103(a)HighSwirl=2.0(b)LowSwirl=0.5L.zI o.(C)I‘—//!!tIii,I,’’’velocitycontourmaxspeed.2.33x1O3(ci)Figure6.16:VelocityVector(cm/s)Plotacrossj=1Planeat15degreeBTDC.;;___.—*k:::-<4-L-._._z:•S.‘—.‘.I(//•.:.L(£—--...--‘•.I..’I,I,i’LIL1:I‘/1/f,._.•_,/f1••II\\....-—I,—I,,•,,maxspeed=2.216x103‘I-__Iy’,I-..-•_.,___,,_,,,//,,it‘H:::::,‘‘,‘,‘,‘I..,fIIt,•S•,/!Ilull.SIII.‘s‘.‘...•SI.•‘..,I/11111III,.siltIII’.‘.‘.\.‘“‘1velocitycontour‘.-“‘IItlIIi,,_______maxspeed=2.348xi0Figure6.17:VelocityVector (cm/s)Plotacroj=1Planeat10degreeBTDC.(a)HighSwirl=2.0(b)LowSwirl=0.5velocitycontour—I,-it,—maxspeed=2.37xi0(c)%*I’..IIII..--.III—4..•_;II(d)CD Cs)(a)HighSwirl=2.0o.•‘1/•I•p..,//1,.•,t,It,”...situ,,.Iii,•4I’•,iII!Ii#,..,,‘.,,I(IiI•.,i,IIIH4_,,I(iIiI ‘.‘‘_‘:iI,•I1Ill.•,Ii••l1Is.•1,ppII\.•,,..,,.,sii)1II’i.._l::‘.‘..‘_‘‘.‘_:‘11’”D’”•\\\%‘‘1IIIi..•IiII,,uiIt‘_•.,..,,,,iiIII“I—.5/1/1/‘—,fl’’..lIJ,IIII1IlIIlFigure6.18:VelocityVector(cm/a)Plotacrossj=1Planeat5degreeBTDC.(b)LowSwirl=0.5Cc)velocitycontourmaxspeed=2.38x iO(d)velocitycontourmaxspeed=2.348x103I-.CD(a)HighSwirl=2.0!HCCC‘%_._,—,-—•......•1l,..1’...,,,_.iIII,.I,.•.,iIlI,.11111,.I••,••L’..IIIi.1•velocitycontourmaxspeed=2.236xiOI-.•,III,..s%s..iII..11111..•IIIII.••.1111).•••JIliii,.•Ig,g1I,.•uI•:____Figure6.19:VelocityVector(cm/s)Plotacrossj=1PlaneatTDC.o.(b)LowSwirl=0.5(C)velocitycontourmax.speed=2.36x103(d)—.I•I——.—Il——I— 0’(a)HighSwirl=2.0(b)LowSwirl=0.511—/1I’\t,IImaxspeed—2.35x103::_‘‘_:(c)(d)JJWi:jiii:III•S_________1II....Fii..—_—...———_,_..I%%ItI•ss1I??....___s,.,,,,i,.,s_’’ss1Is....s.__Fl1Iit.•.sI5—.FiI,II,,),‘..—..——_,,_5....•_____:--—./velocitycontourmaxspeed=2.12x iOFigure6.20:VelocityVector(cm/a)Plotacrossj=1Planeat5degreeATDC.I-ICDDULLV011IJI=C9O.I?ojJ(/tL43)1Opj£pOIA:T9“agOTX961=pctsxw(.————..—•,__.....‘‘‘‘%I:•••II,,,,,,___s__t————s——F———S..—.V•.i#tlI,,._.s__s.••.,,Hl,•.•‘I‘Ii;;-.-.-;:;I,,,,,••.,.,*,II,,htS%115::,,F•_————.-.-I/—•-•-———-’i”•:’;•IF———_••I-‘II’S’’’’I0tX’=PdsXwII’’••‘••IlII•iI1’’’••I1’’’•’•Iji’’’•iV:;•II‘_ssi•II,,II•.11.1I••?I.II•‘‘.ii-.11_II,‘1/(q)O=T’!Sq!H()1flOUO4cpo1aA;III•I—I’••.‘‘—*5•••I(p)=I1!MSMO’JDULLVai2pi;uJdI=uomoJ(sj1u3)1OPA£;DO1A:91UL•.•..,,1!IL.._.‘\•.(p)I/frrrrmTiI1’I••••*‘‘‘‘‘ss*.sSIIIIIII,-11’•II7 YYiIW\II‘I111,,...,’IIII’..——I;;;:::;iI,,,_.•__II,,,ii11I11g..1II.Dll::;‘‘1I•I(q)=1’!S‘l!HsOP<L61=pdsxeui1no1uo(1po1aiiI‘-...—.—\\\1II,III1’1._sl,.IS,,.’’’IIS.:(D)•;•1I::::1OIX63=pdsxuiI.••1I’•II•I•I=T’!”SMO’J•.11I•‘.IIi••l..IIII.i.ItI.•••II.IIIII’.s••.ii‘‘111:.LL’maxspeed=2.25x1030q 1———‘‘gil‘‘Ill.•‘‘SIIu,’,.F,,,I’iIIIiijillIlIlIllivelocitycontourmaxspeed=1.97xLI‘‘‘‘‘‘\ I\xiIjIIjif•IJIjFigure6.23:VelocityVector(an/s)Plotacrossj=1Planeat20degreeATDC.1-’CD CD(a)HighSwirl=2.0(b)LowSwirl=0.5•••—444I—7’_,1.11._I,I—•gIII1—_I’’’(c)>>))-•..I’’’..’•,I!.‘Y.>___i(d)I0)I CD 0) z C, 0; 3 CD 10,0 -a I -v 0 rt 0 C, 0a 3 0 CD -St3 Q0CD a— 00 CD<-CD C,II 00 x2.II.<c3ô0ot%1x0’I CD S Di C, DI CD S1’a,0 -J -4 I -v -J. to 0 C, DI 3 0 CD S:i: 3 C.. o CD 50 w 0 -. I -o -J. U) 0 z C-., 0 CD 5U.c: cDCe,•cqo00)uCDU< 0 I-sxI CD r1 z C-, 0 0 CD -S03 0 -I -Q -I. 0 0 C-)CD -SCDCDa-..II< CDt.:o ‘1x0 I CD z rt C-) a,. 0 CDa, 0 I -Q 1. U) rf 0 z C-) a’ S CDQ bt.t•3CD 0)I— a C) .0) aiii0OCC.,. 2w 0 -I -J.U) 0 D a, a CD-zIC,’III. CD C) [CDI-.OqC.5 zCC.C.)C;-4B.s,- CDo C.) e-.C...0 -J ‘-I I -J. U) C. 0 C-) C) B C- CD 5co(D‘Jo aII<cA(.) C.,.0cJ,.0 I CD -S a, C) a’ 3 CD0 -I -4.1.0 0 D a, 3 0 CD -SCD.-.II< C-..x I.-’0•0Figures 208(a) Bowl—In—Piston Chamber velocity vector•maxspeed351x103• t1ttfS*• ••(b) Re—entrant ChamberFigure 6.32: Velocity Vector (cm/a) Plot across k = 6 Plane at 5 degree ATDC.0III<C3•1x0,I CD C. 0) C-s c, 0’ CD0 -J 1-4 I -v 0 C. 0 C, 0’ 0 CDQ CDFigures 210(a)(b)swirl velocitymm = -1.967x i0max = -1.612x 102L = -1.7864x 10H = -3.418x 102Delta = 1.806x 102Figure 6.34: Swirl Velocity (rad/3) Contour across j = 1 Plane at 20 degree BTDC.swirl velocitymm = -2.28x i0max -1.38x 102L = -2.06x103H = -3.52x 102Delta = 2.14x 102Figures 211(a)swirl velocity< K K (( )> > > 1mm -2.32x10mac=-1.54x102L = -2.lOx iO____________ _______________H = -3.70x10Delta = 2.16x 102(b)____> )Iswirl velocitymm=-2.19x103________________max= -l.33x 102L-l.987x 1OH=3.39x102Delta=2.059x10Figure 6.35: Swirl Velocity (rad/s) Contour across j = 1 Plane at 15 degree BTDC.Figures 212(a)swirl velocityI K < <K<ç)(f)/\7>>> > Imm = -2.35x iO\__________L = -2.13x103H = -3.28x10Delta = 2.25x10(b)1< < <swirl velocitymm = -2.31x iOmax = -0.838x10L = -2.084x103H = -3.06x 102Delta = 2.22x10Figure 6.36: Swirl Velocity (rad/s) Contour across j = 1 Plane at 10 degree BTDC.Figures 213swirl velocitymm = -2.313x10max=- -1.412x10L = -2.096x iOH = -3.585x10Delta = 2.17x 102Figure 6.37: Swirl Velocity (rad/s) Contour across j = 1 Plane at 5 degree BTDC.(a)I <-<<cçç 1(b)‘<1I/ swirl velocity/ miii = -2.36xI I max = -1.47x10f J L = -2.14x103/ / H=-3.68x102Delta = 2.22x10C’CD 0‘1 00 C.0 0 I-i 0 n CD-a.CD C.-a.txjCD n 0 0 0 IIC1llCDi,3eJ’I’C..,Zcx.\,-Oq0I-i0k-ICD CDUUccxx I.00CI,.CD C... U x 0cq.II II“CDXX00CI,tr3I•xI9oq I-.‘•lI-’. 0 0 0 C, CD CD c÷ -a.C, tn CD 3 C, ‘•l• 0CDIIiIll-I-’-_03_CDOOX)<xI-’CDoI CD CD0IIIIt3XXIIII‘ICDXXFigures 216(a)1tkemm = 1.58x iOmax = 5.07x 10L = 6.49x10_______________________H = 4.58x10Delta = 4.91x 10(b)I I___tke= 1.273x10max = 8.14x105L = 9.285x104H = 7.34x105Delta = 8.012x10Figure 6.40: Distributiàn of Turbulence Kinetic Energy across j = 1 Plane at 10 = 1.59x iomax = 3.Olx iOL = 4.44x10H = 2.72x105Delta = 2.85x 104Figure 6.41: Distribution of Turbulence Kinetic Energy across j = 1 Plane at 5 degreeBTDC.Figures(a)217-iL6tke(b)tkemm = 1.292x104max = 3.11iOL = 4.27x10H = 2.81x105Delta = 2.98x1049o.I-.C.,:i. 0 I-. 0 0 0 CD CD 1:-’ oq C,Ci)C, C, I-i 0 II.IIJI•CDàoXXI-•ocI-CDi1L-‘—.IIIIII-‘xo4 CAXxIt.3 I.p.,9I-. ICD I— p., Ii 00 I. x I-’IIII0oCxx’-’I— CDIIGO:•CDCD C.,I-’)—cD-‘xc I-QxII-A ‘-AIIIICO CD9I-.CD ‘—I’llc÷II><-4IIIlz0000-4eCAXi-a00xic•XI-I‘IIHIIP—.III’-‘-‘IIIIoqXt’.)I—’QxrICDII— CDIIoq00c,-‘OCDCD II I-’ Ci3Figures 221(a)[i)_______Ilength scalemm = 1.62x 10—2_max= 1.21L = 1.35x10’H=1.09Delta = 1.19x10(b)length scalemm = 5.77x 10—2____max = 7.59x10’L = 1.28x10H = 6.9x10Delta = 7.0x102Figure 6.45: Distribution of Turbulence Length Scale across j = 1 Plane at 5 degreeBTDC.Figures 22225 Plane A Re-entrant Bowl Chamber(a) Crank AngIe =5 degree BTDCSwirl Ratio = 2.0COMPUTED0.0- -.0.00 0.25 0.50 0.75 1.00 1.20Radial DistancelBowl RadiusPlane A Bowl-In-Piston ChamberCrank Angle = 5 degree BTDC2.5 (b)Swirl Ratio = 2.0U)C0_________________0.00 0.25 0.50 0.75 1.00 1.20Radial Distance/Bowl RadiusFigure 7.1: Comparison of Computed and Measured Mean Radial Velocity.Figures 2232.5 Plane B Re-entrant Bowl Chamber(a) Crank Angle =5 degree BTDCSwirl Ratio = 2.0_ COMPUTEDCO.OMEASUREMSURED-2.5 -_•__5._C0.00 0.25 0.50 0.75 1.00 1.20Radial DistancelBowl RadiusPlane B Bowl-In-Piston Chamber0 Crank Angle =5 degree BTDC2.5 (b)SwIrl Ratio = 2.00COMPUTED0.0- - - - - --2.5-5.00.00 0.25 0.50 0.75 1.00 1.20RadIal Distance/Bowl RadiusFigure 7.2: Comparison of Computed and Measured Mean Radial Velocity.Figures 224Plane C Re-entrant Bowl Chamber(a) Crank Angle =5 degree BTDCSwirl Ratio = 2.00.0 .-- -- COMPUTE]C-50‘U.75 I0.00 0.25 0.50 0.75 1.00 1.20Radial Distance/Bowl RadiusPlane C Bowl-In-Piston Chamberb Crank Angle = 5 degree BTDC2.5 (>Swirl Ratio =2.0Cl)0.00 0.25 0.50 0.75 1.00 1.20Radial Distance/Bowl RadiusFigure 7.3: Comparison of Computed and Measured Mean Radial Velocity.Figures 225Plane A Re-entrant Bowl Chamber2.5 (a) Crank Angle = TDCSwirl Ratio = 2.00C-7.50.00 0.25 0.50 0.75 1.00 1.20Radial Dlstance)Bowl RadiusPlane A Bowl-In-Piston Chamber0 Crank Angle = TDC(b)Swirl Ratio 2.00C0.00 0.25 0.50 0.75 1.00 1.20RadIal DistancelBowl RadiusFigure 7.4: Comparison of Computed and Measured Mean Radial Velocity.Bowl-lnPIston ChamberCrank Angle 5 degree BTDC0.50 0.75 1.00 1.20Radial Distance/Bowl RadiusFigures 226I10.0C000C0C00Radial DistancelBowi Radius20(b)Plane ASwirl Ratio 2.0_-_.COMPUTEDMEASUREDCl)C00C0080C0CC0015105- - — - -- aI0 0.25Figure 7.5: Comparison of Computed and Measured Tangential Velocity.12.510.0C8C00)CC0Figures 227Radial DistancelBowl Radius15Plane BSwirl Ratio = 2.0(b)Bowl-In-Piston ChamberCrank Angle = 5 degree BTDC10_-_- COMPUTEDMEASURED5— - - - -— — — —C 025 0.50 0.75Radial DistancelBowl Radius1.00 1.20Figure 7.6: Comparison of Computed and Measured Tangential Velocity.Figures 22815.0112.510.0C7.55.02.5Radial DistancelCylinder RadiusPlane CSwirl Ratio 2.0(a)Re-entrant Bowl ChamberCrank Angle = TDC- COMPUTED• MEASURED0.00.00 025 0.50 0.75Radial DistancelBowl Radius1.000.00 025 0.50 0.75 0.90Figure 7.7: Comparison of Computed and Measured Tangential Velocity.Figures 229Re-entrant Bowl ChamberCrank Angle = 5 degree BTDCPlane ASwirl Ratio = 0.5(a)I--ja-- COMPUTEDL. MEASURED1.501.2511.000.75jo.500.250.000.00- a--a---0.25 0.50Radial Distance/Bowl Radius0.75 1.Plane ASwirl Ratio = 0.5- (b)-1.50C00CCD1.00.CCDo.5oBowl-in-Piston ChamberCrank Angle =5 degree BTDC- .a - --0---9--___-COMPUTED• MEASURED- .9_ ---a--.0.00.00 0.25 0.50 0.75 1.00Radial Distance/Bowl RadiusFigure 7.8: Comparison of Computed and Measured Turbulence Intensity.Figures 230Plane B Re-entrant Bowl ChamberSwirl Ratio = 0.5 Crank Angle 5 degree BTDC(a)CoC___________________1.0 - -- -- COMPUTED• MEASURED0.00 0.25 0.50 0.75 1.00Radial Distance/Bowl RadiusPlane B Bowl-In-Piston ChamberSwirl RatIo =0.5 Crank Angie—S degree BTDC1.50 -COMPUTED• MEASURED.O.50-I-0.00 .1.. .1..,,0.00 0.25 0.50 0.75 1.00Radial Distance/Bowl RadiusFigure 7.9: Comparison of Computed and Measured Turbulence Intensity.Figures 231Plane CSwirl Ratio = 0.5(a)Re-entrant Bowl ChamberCrank Angle =5 degree BTDCaJ0.5V_D_ COMPUTED0.0 —0.00 0.25 0.50 0.75 1.1Radial DistanceiBowi RadiusPlane CSwIrl Ratio = 0.5.50 (b)Bowl-In-Piston ChamberCrank Angie = 5 degree BTDC2.00LI--- COMPUTED. MEASURED0.50-t0.00 —0.00— -—B—0.25 0.50Radial DlstancelBowl Radius0.75 1.00Figure 7.10: Comparison of Computed and Measured Turbulence Intensity.Figures 232Radial DislancelBowl Radius0.0 —0.00 0.25 0.50Radial DlstancelBowi Radius0.75 1.00Plane ASwirl Ratio = 2.0(a)Re-entrant Bowl ChamberCrank Angle 5 degree BTDCCOMPUtEbJ• MEASURED I2.0________e 0.5---a.-0.0 —0.00 0.25 0.50 0.75 1.00Plane ASwirl Ratio = 2.0(b)COMPUTED• MEASUREDBowl-In-Piston ChamberCrank Angle = 5 degree BTDC- .Q ---a--—-C--a-2.0:>0.5— a -Figure 7.11: Comparison of Computed and Measured Turbulence Intensity.0.50Radial Distance/Bowl RadiusFigures 233Plane B Re-entrant Bowl ChamberSwirl Ratio = 2.0 Crank Angle 5 degree BTDC2.0 - (a)-- -- COMPUTED1.5 - • MEASURED1.00.5 -0.0 —0.00--C.-- a- -— —a. —- a- -0.258.C,)C0(I,0CCD5;(0CCDC8CID1.0C8CCD.o 0.520.75 1.00Plane BSwirl Ratio = 2.0(b)Bowl-In-Piston ChamberCrank Angie = 5 degree BTDCCOMPUTED_. MEASURED0.0 —0.00 0.25 0.50 0.75 1.00Radial Distance/Bowl RadiusFigure 7.12: Comparison of Computed and Measured Turbulence Intensity.Figures 234Plane CSwIii Ratio = 2.0Re-entrant Bowl ChamberCrank Angle 5 degree BTDC2.01.51.0--a-- COMPUTED_____MEASURED- -0.5 -0.00.00 0.25-gU)C0Co0C0CoC0CCCo-e2.0[.51.0CCo0C0.0 0.5I-0.50Radial Distance/Bowl Radius0.75 1.00COMPUTED• MEASUREDPlane C Bowl-in-PIston ChamberSwirl Ratio = 2.0 Crank Angle = 5 degree BTDC0.00.00 0.25 0.50Radial Distance/Bowl Radius0.75 1.00Figure 7.13: Comparison of Computed and Measured Turbulence Intensity.Appendix AGOVERNING PARTIAL DIFFERENTIAL EQUATIONSAd General FrameworkIn common with most researchers, we work with average forms of the differential conservation equations of mass, momentum, energy and individual chemical species to escapethe burden of having to calculate the fine-scale turbulent motions. Here the equations aredensity-weighted and ensemble-averaged, the relations between the instantaneous value4 of a dependent variable and its density-weighted ensemble-average (DWEA) value,being:= n!o{2m1+ 1 n:h3(t +nTo)q5(,t +nTo)} (A.1)• where the summation is over an ensemble of samples taken at given position and timet in the engine cycle, whose period is T0. The symbol denotes the ensemble-averagedensity.This averaging process gives rise to additional unknowns in the equations, some ofwhich are conventionally interpreted as representing turbulent fluxes of the entity inquestion. These are presumed to obey constitutive relations similar to laminar flows, withthe molecular transport coefficients replaced by turbulent or eddy coefficients. The latterare linked to two DWEA parameters of the turbulence field (k and e), which representthe dependent variables of a turbulence model comprising two additional differentialtransport equations.235Appendix A. GOVERNING PARTIAL DIFFERENTIAL EQUATIONS 236A.2 Equations Governing Fluid MechanicsWithin the framework given above, the DWEA conservation equations of mass andmomentum can be derived as shown below.We start with the partial differential equations expressing mass conservation and thebalance of linear momentum; the compressible form of these equations is used. For compactness, the cartesian tensor notation is employed, with the usual summation conventionover repeated indices:Continuity Equation+ .9(pU1)= (A.2)atMomentum EquationOpU2 + 8(pU2U,) = + (A 3)at a a a,Here T32 is the deviatoric (zero trace) part of the stress tensor, p is the fluid density,p is the pressure, and U1 is the i’ component of the velocity vector ZT; body forces havebeen neglected. For a Newtonian fluid of viscosity , rjj is given by;1ou öU’ 2 c9Urjj = + — p-—S, (A.4)where Sj is the Kronecker delta. For compressible flows, these equations do notcomprise a closed set; an energy equation and equations of state (not discussed here)complete the set of fundamental equations.Appendix A. GOVERNING PARTIAL DIFFERENTIAL EQUATIONS 237A.3 Equations for Mean QuantitiesWe define two averaging operations to transform the above equations for instantaneousquantities into equations for ensemble-averaged quantities. Mean quantities are to beinterpreted as ensemble- averaged values at a given crank angle. The conventional averageis denoted by overbars while tilde is used for the density-weighted average (the latterbeing convenient in variable-density flows). Thus p, p and U can be decomposed intomean and fluctuating components as follows: for conventional averaging;p=+p’ (=O)p=+p’ (=O)U=U+u (=O) (A.5)where single primes (or lower case u1) denote the fluctuations; and, for density-weighted averaging,U1=U++u’ (=O) (A.6)The conventional mean of density-weighted fluctuations is not in general equal to zero.The mean continuity equation results from taking the mean of Equation (A2) (note thattaking the derivative with respect to t or X: commutes with either averaging procedure),(A.7)atwhile the mean momentum equation from Equation (A3) can be written as,+ O(U,U) =+ aIf (A 8at ôz, 8x ox,Appendix A. GOVERNING PARTIAL DIFFERENTIAL EQUATIONS 238Here the effective stress is the sum of a turbulent stress, or Reynolds stress,resulting from averaging the non-linear convective term in Equation (A3) and the meanviscous stress:(A.9)where,—uii’ (A.lO)A.3.1 Turbulence ModellingThe turbulent stress represents transport of mean momentum by turbulent velocity fluctuations. It is modelled by drawing an analogy between the turbulent transport ofmomentum and viscous transport. A turbulent viscosity is introduced (see EquationA4),2 10(13 8U1’\ 2 OUTt31 = —kS1+ p + — (A.11)where k u’u’/2 is the turbulent kinetic energy. The value of is related to localturbulence velocity and length scales. These turbulence scales are formed from k andits viscous dissipation rate e. A turbulence velocity scale is k1”2 while a length scaleis proportional to k2/3e. In this k - model, the turbulent viscosity and effectiveviscosity 1.Leff are,k2ptCpm PeffILt+IL (A.12)Appendix A. GOVERNING PARTIAL DIFFERENTIAL EQUATIONS 239where c is a model constant. To complete the model, equations governing the evolution of k and are needed. These equations are:03k 0U,k — a (‘ Ok 0(1, —j i +TtjP€ (A.13)at &, Ox, \ Ox,j Ox18€ &Ue 0 (Lff £96’ oU aU—— I —:--- I — C2P— c3pe— ci—Tt1Ox, Ox,‘a Ox, k Ox, Iv Ox1In Equations (A13) and (A14), c1, c2, c3, 0k and o are model constants: values usedin the present study are summarized in Table 4.1In Equation (A13), the penultimate term represents the production of turbulencekinetic energy by interaction of the turbulent velocity fluctuations with mean velocitygradients. It is useful to expand this term into three parts using Equations (All) and(A12), thus:Tt,1—Oxlou, OU OU 2 0(13 2 OUOU= t+ h—) -h—- — —= P1+ P2+ P3 (A.15)These three contributions, P1, P2 and P3 will be referred to as the first, secondand third production terms, respectively. The first (P1) represents conversion of meanto turbulent kinetic energy by mean shear, or shear production, it is usually positive,although may be negative in isolated regions of the flow. The second production termrepresents turbulence generation by dilatation of the velocity field; the sign of P2 isopposite to that of the local mean velocity divergence. Production term P3 (dilatation-squared term) is always negative.The equations above have their origins in a model described by Launder and Spalding86T)ia‘(i)sauo‘sjajias£qjqtss:duoooppux£junb-qnspuSMOJqsardmoDuTjoxuoUtpdopap£Tru10SMtpt{M‘(6I)SNOLLVflb7VLLN37J3TJIU7VLL2IVJDNIN?1AODVXrptTddVAppendix BPROCEDURE FOLLOWED TO CONVERT KIVA-Il CODEThe code conversion was necessary because the KWA-II code was written specificallyfor use on the CR1 cray family computers, operating under the Cray Time Sharing System(CTSS); and users who do not have access to a CRAY computer must have to undertakethe code conversion excercise in order to avail themselves of the opportunities of usingthe KIVA-Il code.Figure Bi shows a sample of the results obtained from both CRAY and SUN-SPARCcomputer systems. These results were obtained by using identical input data sets. FromFigure Bi, we see that the results obtained from the SUN computer system accuratelyreproduced those got from the CRAY computer system.B..1 CRAY Fortran IntrinsicsThe coding changes that took the greatest amount of time were due to the use of theCRAY intrinsic functions, namely; CVMGM, CVMGP, CVMGT and CVMGZ, aswell as their integer versions. These functions are non-standard fortran coding statementsand they therefore needed to be replaced. Versions of each were written in Fortran77(f77); however, the intrinsics in CRAY Fortran (cft) may have been written to allowvarying data types as arguments, since both reals and integers were used in differentinvocations of the same function. In any case, rather than provide code to allow differentvariable types as arguments, the function arguments were changed to double-precisionfor the real functions and integer for the integer functions. This entailed an examination241Appendix B. PROCEDURE FOLLOWED TO CONVERT KIVA-Il CODE 242of each CALL. The file containing these routines is called otherdp.f.Several other intrinsic functions required changes as well; specifically amaxl andaminl were changed to max and miri respectively. The KIVA-JI code has references toroutines timeh, dateh and second, these routines provide information about time anddate of a run, as well CPU accounting information. Fortran77 versions of these routineswere written, these are contained in file dummy.f.B.2 Block Data SubroutineCRAY Fortran allows variables to be initialized within a subroutine and then passedusing COMMON statements. Fortran77 does not allow this and so a number of arrayinitializations had to be moved to fuelib.f , the file containing the block data subroutine.They are marked with comments in that file.B.3 Real Variables are Double-PrecisionIn CRAY FORTRAN (cft), the default real data type is 64 bit, or what in f77 is doubleprecision. To make sure that all real variables were declared properly, an implicit double-precision statement was included in the comd.h include file. It was necessary to makesure that all subroutines contain this statement in the converted KIVA-Il code. KIVA-Iluses a number of routines from the SLATEC Mathematics Subroutine Library. The codedistributed with KIVA-Il contained single-precision versions of these routines; double-precision versions of these routines are contained in the subdirectory slatecdp.B.4 Comdeck StatementThe comdeck statement in cft is equivalent to include in f77; all references were thereforechanged from comdeck to include in the converted code.Appendix B. PROCEDURE FOLLOWED TO CONVERT KWA-Il CODE 243In subroutine begin, the cft intrinsic locf is used to zero all values in common blocks.These statements were commented out, since f77 initializes to zero by default. If thiswere not the case, either the scheme used in the original code must be set up or eacharray and variable must be set to zero individually.B.5 Debugging, Organizing and Testing of the Converted CodeThe original files were split into individual files for each subroutines, using the 3plit anduiility. Revision Control System (RCS) was used to keep track of code changes.Appendix CDETERMINATION OF SWIRL RATIOSwirl is usually defined as organized rotation of the charge about the cylinder axis.• Swirl is created by bringing the intake flow into the cylinder with an initial angular momemtum. While some decay in swirl due to friction occurs during the engine cycle, intakegenerated swirl usually persists through the compression, combustion and expansion processes. In engine designs with bowl-in-piston combustion chambers, the rotational motionset up during intake is substantially modified during compression. Swirl is used in dieselsand some stratified-charge engine concepts to promote more rapid mixing between theinducted air charge and the injected fuel. Swirl is also used to speed up the combustionprocess in spark-ignition engines. In two-stroke engines it is used to improve scavenging.In some designs of prechamber engines, organized rotation about the prechamber axis isalso called swirl. In prechamber engines where swirl within the precombustion chamberis important, the flow into the prechamber during the compression process creates therotating flow.C.1 Swirl MeasurementThe nature of the swirling flow in an actual operating engine is extremely difficult todetermine. Accordingly, steady flow tests are often used to characterize the swirl. Airis blown steadily through the inlet port or valve assembly in the cylinder head into anappropriately located equivalent of the cylinder. A common technique for characterizingthe swirl within the cylinder has been to use a light paddle wheel, pivoted on the cylinder244Appendix C. DETERMINATION OF SWIRL RATIO 245centreline (with low friction bearings), mounted between 1 and 1.5 bore diameters downthe cylinder. The paddle wheel diameter is close to the cylinder bore. The rotation rateof the paddle wheel is used as a measure of the air swirl. However, in the present work,the tangential velocity, v0, exiting the intake port was determined as shown below, andthis was used as a measure of inlet swirl. Schmatic diagrams depicting the determinationof swirl are shown in Figures 3.2 and C.1.The characteristic velocity, v0, can be derived from the pressure drop across the portusing an incompressible flow equation:1/2V0[2(o_c)] (C.1)or compression flow equation:=2-y — ()(Y_l)/Y1/2(C.2)((7—l)Po \Po 3where the subscripts 0 and c refer to upstream stagnation and cylinder values, respectively. The difference between equations C.1 and C.2 is usually small. Equation C.1 wasused to estimate the characteristic velocity in the present study.When swirl measurements are made in an operating engine, a swirl ratio, R, is nor-• ma.lly used to define the swirl. It is defined as the angular velocity of a solid-body rotatingflow w3, which has equal angular momentum to the actual flow, divided by the crankshaftangular rotational speed, WE, (2IrNE); that is,R= (C.3)2 7rNESteady-state flow rig data were used to estimate swirl, w3 in the following manner.Assuming that the port retains the same characteristics under the transient conditionsof the engine as on the steady-flow rig, the equivalent solid-body angular velocity w, atthe end of the intake process is given by:Appendix C. DETERMINATION OF SWIRL RATIO 246(C.4)where B is the cylinder bore. Equation C.3 was then used to estimate R9. The swirlratios used in the present study are (0.5 and 2.0).During the induction stroke in an engine the flow and the port open area, and consequently the angular momentum flux into the cylinder, vary with crank angle. Whereasin rig tests the flow and the port open area are fixed and the angular momentum passesdown the cylinder continuously, in the engine intake process the momentum producedunder corresponding conditions of flow and valve lift remains in the cylinder.The relationship between steady-flow rig tests (which are extensively used becauseof their simplicity) and actual egine swirl patterns is not fully understood. Steady-flowtests adequately describe the swirl generating characteristics of the intake port and valve(at fixed valve lift) and are used extensively for this purpose. However the swirling flowset up in the cylinder during intake can change significantly during compression.*1111,s%I,..IIIIiilli•....‘.‘iIiIl1??’....S%%%llIiliiI/lu•s%%‘Iiihi,S.’’’’IçIi,?..I.,,Iilsi,.....a.,..hut...,ai11,la..l%1lIIt..,.i,,,lliii%ffjI..a1*11111111taSiiii‘111111111111111111*SIIIIIIIlillIlltillIII11111111I‘CRA‘iiiUThISUNSPARCFigureB.1:CRAYandSUNSPARCVelocityVector at54degreeBTDC.fliiiillltThIIHIIIIIIhiiIIIIIIIITT1111111111117f[11‘.5.4‘‘‘*‘‘k’lII,’IIIu,UH1IZII1-4-4-1_IIhiiiiiiiWiUiJUII!T7CRAY11111111IIii’Oq ‘1 Cb 9:SUNSPARC/_-.—_____—.,,S............%%-..••..-—.1_a—I—I——-,,.{iI—._—————S.’.’/ I,,.5Ii,,...,I,alitIIIlI,III,iIIFigures (Appendix) 248COMBUSTIONCHAMBERCOMPRESSED AIRp0COMPRESSED AIRCYLINDERFigure C.1: Schematic diagram of steady flow swirl apparatus.


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items