TECHNICAL REPORT 1 SPARK IGNITION ENGINE SIMULATION MODELS C. Borgnakke P. Puzinauskas Y. Xiao Department of Mechanical Engineering & Applied Mechanics University of Michigan.For the Project: Development of a Spark Ignition Engine Simulation Model and Simulation of Intake Flow Heating Awarded by: The Chrysler Challenge Fund Project No.. 2000511 February 1986

1 <J 1 ) 3 i

Abstract The current project has been awarded to the Mechanical Engineering and Applied Mechanics Department of The University of Michigan by the Chrysler Corporation under the Challenge Fund Program. The purpose of the project is to develop a spark ignition engine simulation which predicts indicated performance over a wide range of operating conditions. The creation of such a simulation requires the development of models which describe all of the essential physical phenomena occuring in the combustion chamber, and structuring them to effectively interact with one another. The backbone of the simulation is a lumped thermodynamic model that gives an overall description of the changes in mass, composition, momentum and energy inside the combustion chamber'during the cycle. The analysis is mainly focused on energy related aspects such as volumetric efficiency, rate of heat release and heat losses. This information is utilized to calculate indicated work. Many of the different physical processes are described separately in the analysis, and in the model they are put together as interrelated modules. The present report gives the overall structure of the model together with a literature review of the current state-of-the-art in the modelling of the various processes. The basic fluid mechanics and thermodynamics is emphasized with little coverage of the detailed chemistry. The latter subject, because of its importance for the goal of this investigation, will be treated in a later report. i

SPARK IGNITION ENGINE SIMULATION MODELS CONTENT i ABSTRACT 1 GENERAL MODEL REVIEW 1 1.1 Various model formulations 1 1.2 Zero dimensional models 1 1.3 The mathematical formulation of multi-zone models 2 1.4 The simplified form for two zone models 5 2 SUBMODEL STRUCTURE AND GEOMETRY 12 3 THERMODYNAMIC PROPERTIES 18 3.1 The unburned gas properties 18 3.2 The burned gas properties 19 4 COMBUSTION MODELS 23 5 HEAT TRANSFER MODELS 30 6 VALVE FLOW MODELS 7 IN-CYLINDER MEAN FLOW MODELS 8 TURBULENCE MODELS 51.9 MANIFOLD FLOW MODELS 64 ii

NOMENCLATURE b any conserved property per unit mass BDC bottom dead center c speed of sound CR compression ratio e total specific energy (u + ~v2 + gz) h specific enthalpy hot specific stagnation enthalpy (h + 2v2) k conductivity K turbulent kinetic energy I turbulent length scale m mass P pressure q heat flux per unit area Q total heat flux r radial coordinate R gas constant s specific entropy Sb source of property "b" per unit mass Sb total source of property "b" SL laminar flame speed ST turbulent flame speed ST stroke t time T temperature TDC top dead center u specific internal energy u velocity component in x direction ut fluctuating component of u (rms) v velocity component in y direction V volume w velocity component in z direction x cartesian coordinate y cartesian coordinate z cartesian coordinate w angular velocity p density 0 crank angle in degrees r turbulent time scale rij turbulent shear stress iii

Subscripts b burned property ex exit condition FL flame location or condition i direction i or initial condition ig value at ignition in inlet condition j direction j L laminar property s value at spark timing T turbulent property u unburned property iv

Chapter 1 GENERAL MODEL REVIEW 1.1 Introduction The modelling of spark ignition engines has developed substantially from the thermodynamic analysis of the air standard Otto cycle. Charging of the cylinder with fresh mixtures, finite rate of combustion, heat losses and unburned charge fraction are a few of the essential aspects relating to the thermodynamic efficiency that have been added to the basic model. More recently, increased emphasis is also placed on the detailed chemistry to improve the understanding of the emissions or the undesirable combustion abnormalities. These abnormalities include misfire and knock which are caused by failure of combustion and excessively rapid combustion, respectively. The proper design of an engine and its control strategy often ends up as a trade off between its efficiency, emissions and many practical constraints. It is in this context that a reasonable model must try to predict the performance of a particular design to help in the development of engines. Many models have been developed to emphasize the influence of each specific process. Assembling different combinations has resulted in a substantial number of different engine simulation models. In the present report we will first classify the models according to the overall structure, which essentially becomes a mathematical classification. This has a profound impact on the method of solution and the phenomena that can be investigated with a given model. In a later section we then deal with each of the physical processes in more detail and the way they are being modeled. 1.2 Model Formulations Engine simulation models are build up around an analysis of what happens in the combustion chamber. This does not imply that the processes outside the cylinder are insignificant, but that these act as boundary or initial conditions for the combustion event. To perform such an analysis and formulate it as a mathematical model the 1

general laws are written down for one or more chosen control volume(s). These are the conservation equations for mass, momentum, energy and entropy, which are valid in all circumstances. These laws contain more unknown terms than equations, therefore to close the set of equations, we must introduce the auxiliary laws and the constitution equations. The auxiliary laws describe the terms which control some of the physical phenomena in the general laws: These are usually formulated as a model of a process and therefore are also referred to as submodels. Examples of these include the heat transfer model, valve flow, mean cylinder flow and turbulence. The constitutions are laws or correlations that expresses the behavior of a substance or system.Some examples are the equation of state, thermodynamic properties, chemical reaction rates and piston motion. The overall model structure is then built up around the general laws for a number of control volumes in the following manner; Control Volume 1 Control Volume 2 etc. General Laws General Laws Auxiliary Laws Auxiliary Laws Constitutions Constitutions - The form of the general laws will be shown in the following section and the auxiliary laws and constitutions will be shown in later chapters. 1.3 The General Laws. The general laws are express the conservation of several properties such as: 1. The continuity equation (m) 2. The momentum equations (u,v,w) 3. The energy equation (e) 4. The entropy equation (s) The second law of thermodynamics is not used in general as are the first three laws, but used only for a limited set of processes which we may regard as reversible. The equation can be used to calculate the increase in entropy to show that processes indeed are possible and the instantaneous entropy could express a possible availability. We know however, that several of the important processes (combustion, heat transfer) are irreversible and thus render such an analysis questionable. It will therefore be omitted in the following discussion and used later as appropriate. The major distinction between the model formulations comes from the choice of control volume. This can be taken as a differential volume in which case a full three 2

dimensional form appears. All the general laws are then partial differential equations of first order in time and up to second order in the spatial derivatives. The differential equations can also be integrated over part of the 3-D volume so the control volume(s) is finite in at least one direction. This corresponds to an averaging over part of the domain or direction and reduces the detail of the spatial variations that can be analyzed. The choice of model formulations are then classified according to the number of spatial directions in which the differential formulation is retained. We thus have 3,2,1 and O dimensional models ranked with a decreasing mathematical complexity. In the zero dimensional models there are no partial differential equations with all variables averaged over a finite volume, which is the integral method referred to as a thermodynamic model. This gives the following two classes of models: 1. Multi dimensional models 2. Zero dimensional models (Thermodynamic) The benefits of the thermodynamic models are simplicity and easy solutions - implementations, the advantages of the multi dimensional models are spatial resolution and determining the effect of possible complex geometries. The drawbacks of the thermodynamic models are oversimplifications (i.e. averaging of non-linear properties) and no consideration of geometric variance, whereas for the multi dimensional models the drawbacks are complex solutions and implementations and increased uncertainty in the modelling process, particularly for the turbulent reaction rate. The general laws have the same appearance whether they are formulated as differential or integral equations. They give the rate of change of a property in that control volume equal to changes due to flow in/out, diffusion of the property and a possible source (production) of that property. The form of these equations is illustrated below for the three dimensional case and the integral approach. The 1- and 2-dimensional cases result in a combination of these two forms. We show here the continuity equation and an equation for a property "b" that may stand for a velocity, energy or entropy per unit mass of the mixture. 3-D Differential Form Continuity Equation: -P + V(pvi) = (1.1) Conservation of b: ab Pt + V(pv'b) = V(DbVb) + psb (1.2) In these forms Db is a diffusion or transport coefficient and Sb is a source of b. Table 1 gives the diffusion coefficient and source terms, when the property b is concentration, 3

velocity and energy, respectively. In the case of a turbulent flow the diffusion coefficient becomes an effective transport coefficient, which using standard turbulent modelling is a gradient diffusion VT Db = Dbam + P- (1.3) Orb where VT is the turbulent kinematic viscosity and ab is the appropriate Prandtl or Schmidt number for the property b. Several examples of analysis using turbulence modelling can be found in the literature, see references where it is nearly exclusively done with a k - e model. For those models equations for the turbulent kinetic energy k and its rate of dissipation e are solved from similar equations. Integral Form Continuity Equation: dt im i- Ehez* (1.4) Conservation of b: dB T = E Zhinbin-E rzbez + Fb + Sb (1.5) Here, m and B are the total mass and amount of b in the control volume, the sums are taken over all boundaries with flow across and Fb is a diffusional flux of "b" at the boundary. The F and source term S are also shown in table 1 for various b". Table 1. Property b Db Fb Sb Sb Concentration c D 0 0 rc PAFLST/m Velocity u v Ar -dP/dx (Force). Energy e k Q - W rchRp pAFLSThRP All the models that are integrated in one or more directions may be formulated as a multi-zone model. This means that the averaging process is carried out in each zone separately with that zone as a control volume. In the limit of many zones with divisions in all three directions such a method is identical to a finite difference or finite element method in three dimensions. In this work we will use the zonal structure to handle volumes near the walls which are expected to behave differently from the bulk of the gases in the combustion chamber. 4

1.4 The Zero Dimensional Model When the control volume(s) is taken as the whole combustion chamber in the integral formulation, we have a zero dimensional model. This type of model is commonly referred to as a thermodynamic model or a phenomenological model. Most properties are averaged over the total cylinder volume and no spatial information is available. A survey of the thermodynamic models are presented by Blumberg, Lavoie and Tabaczynski [1], Mattavi et. al. [2] and Heywood [3] among others. These models rely on some understanding of the physics involved and try to capture the main features of the processes. By including the description of the most important aspects the models have performed surprisingly well and are idealy suited for parametric studies. Models in use of this type are implemented as multi-zone models in the recognition of the spatial distribution of properties. During combustion the charge in each control volume is divided into two zones one with unburned gases and one with burned gases. Before and after the combustion period the gases are assumed to be uniform over the control volume. The simplest model of this type has one control volume as the whole combustion chamber, with at most two zones present. Several different models have been proposed, see [1-10], and they show reasonable agreement with experimental data, when compared with overall measurements such as the pressure. Extending the model to contain two control volumes, the mean flow can be calculated as a non - solid body rotation [11-13] with the possibility of doing the same for the turbulence field. Other zones can be added for the more realistic computation of the wall heat transfer rates or a configuration with a prechamber. To cover those kinds of models the equations for the multi-zone models are shown in the next section. 1.5 Multi-Zone Model Formulation Assume the combustion chamber is divided into N control volumes, and let us write the continuity equation, energy equation and constitution for each of these. We will take the pressure to be uniform over the whole chamber, and assume a uniform condition over each of the zones. The integrated (lumped) formulation, with "i" denoting the zone i, then becomes: Continuity equations: dt mini- mhezi (1.6) Energy equations: dE iceVi dt =,rhinhi,i h mihez,i + Q P dt (1.7) Constitution: 5

PVi miRiTi (1.8) Constraint: Pvtot- = = Z(mi,Ti) (1.9) Integration in time of these equations will yield the mass from the continuity equation, the temperature from the energy equation and the chamber pressure from the constraint of the total volume. In order to solve these equations, we rewrite the energy equation in terms of temperature as follows: dEi dh V dmi dhi dVp dP d (mihi - PV ) =t- dt + mi — P - t (1.103 where dhi dTi dt = CP' dt Substituting this expression for the left hand side of the energy equation we may cancel the work term and substitute the rate of change of mass from the continuity equation. All this leads to an energy equation in the following form: dTj dP miCp,i d = E hin,i(hin-hi) + Qi - E -hez,i(h - hi) (1.11) Since the properties are not distributed over the control volume the value of the enthalpy that leaves the control volume is usually taken as the average, unless other information is available. Then the whole term with the mass flow out cancels, and the average enthalpy in the control volume is not affected by the mass leaving. To solve the system of equations, the rate of change of pressure is found from the constraint differentiated with respect to time. Summing over the control volumes to give the rate of change of the total volume given by the engine geometry, yields the equation d P~~i d+ ~ d i? T(1.12) andt dt -+= mPtot R'+ T+ (1.12) The variation in the gas constant R. is written from the sensitivity of R to the temperature and the pressure and finally the energy equations are substituted to elliminate the rate of Ti. The equation now has the form dVtot dP Pdat =A+B dt (1.13) where all the summations over "i" of the various terms are in the values A and B. This equation is solved for the time derivative of the pressure, which can then be substituted back into the energy equations and give the rate of change for the temperatures. 6

The formulation presented here can handle any specified number of zones, as long as all the terms expressing exchange of properties between the zones can be reasonably well modeled. These terms are the mass flow rates in and out of the zones, the diffusion of species and the conduction of energy. With many zones these quantities may be difficult to estimate if the velocity field is not found with any spatial details. The purpose of writing the general form is to include only a few zones for the bulk of the gases, where these terms can be modeled with some degree of confidence. The additional zones are planned to be for the wall boundary layers and the crevice volumes. The current knowledge indicates that these volumes are important for an understanding of the emissions and of course the heat losses in the engine. 1.6 The Two Zone Model Currently most thermodynamic models employ a two zone model analysis of the engine. The two zones are only present during combustion, where they are selected as the unburned and burned gas volumes. The separation between the zones is the propagating flame, which for the spark ignited gasoline engine have been shown to be reasonable well defined. Photographic recordings show a nearly spherical flame propagation from the spark location with a relatively thin flame structure. In engines with swirl or complex geometry the center of the spherical flame may not be the spark plug location, but a point that changes with time during the combustion period. If the overall flame propagation and the heat transfer is not strongly affected by the relative location of the spherical flame in the combustion chamber a model with a fixed origin will be quite sufficient. The formulation of the two zone model is a simplification of the more general multi zone model presented in the previous section. The equations are then obtained with the zone "i" indicating unburned gases, subscript "u", and burned gases, subscript'b". To present the model in a simple form let us assume that the gas constant, R, is a constant and not a function of temperature or pressure. There is then only one mass flow rate between the zones and that is the rate of combustion converting unburned gases to burned gases. Since the flame propagation is the rate at which the unburned mixture is converted into the composition of the burned gases there is no heat transfer between the zones, except if radiation is considered. For the gasoline engine radiation is not a major contribution to the change in the energy levels and will thus be neglected. The general laws can then be written for the two zones, with essentially three equations for each zone as the continuity, the energy and the constitution equations. Unburned gases: dmu dm= -rob (1.14) mrCpu dTuv (1.15) =Q +V dt

PVT, = m= RuTu (1.16) Burned gases: drb = rate of combustion = rhb (1.17) dTb dP mbCb dT= - mb(hu - hb) + Qb + Vb dP (1.18) ma~p~b dt dt PVb = mbRbTb (1.19) Constraint: dVtot u + (1.20) dt dt dt Differentiating the two constitutions and substituting them into the constraints now yields dVtot drb 1 1 P T1 Tb _t( - -U) -Vtot- + VuU + Vb (1.21) dt dt Pb i P TU Tb Substitution of the energy equations to elliminate the the rate of change of the temperatures provides the equation for the rate of change for the pressure as Ru Rb dP', rnu dVtot Ru (Vot-V -Vb Pb P — ) ) + QUC PIU Cp,b Pb P dt cp,u + [Qb + (hu - hb)b] b (1.22) The value of dP/dt can then be substituted into the two energy equations. As the equations are integrated in time the current values of mb, Tu and Tb are obtained from the integrator, whereas the pressure then is calculated as p mURuTu + mbRbTb (1.23) Vtot This is done to prevent any inconsistencies that might arise if the pressure was taken from the integration of the pressure rate equation listed above. It also ensures that the constraint is satisfied accurately. During all the other phases of the cycle there is only one zone present in the chamber. This means that there will be one continuity equation for the total mass, one energy equation for the homogeneous charge temperature and the volume constraint to determine the pressure. The mass in the chamber only changes during the scavenging period, 8

if the blowby is neglected, so the continuity equation does not have to be integrated when the valves are closed. This set of equations that determines the level of the properties must then be accompaniecd by the equations for all the submodels describing all the terms on the right hand side of these equations. These submodels will be covered separately in the following chapters.

References. Zero Dimensional Models. 1. Blumberg, P.N., Lavoie, G.A. and Tabaczynski, R.J, "Phenomenological Models for Reciprocating Internal Combustion Engines", Prog. in Energy and Comb. Sci., 23, 1979. 2. Mattavi, J.N., Groff, E.G., Lienesch, Matekunas, F.A. and Noyes, R.N., "Engine Improvements through Combustion Modelling", in Combustion Modeling in Reciprocating Engines, eds. Mattavi and Amann, Plenum Pres, 1980. 3. Heywood, J.B., "Engine Combustion Modeling-An overview', in Combustion Modeling in Reciprocating Engines, eds. Mattavi and Amann, Plenum Pres, 1980. 4. Lavoie, G.A., Heywood, J.B. and Keck, J.C., "Experimental and Theoretical Study of Nitric Oxide Formation in Internal Combustion Engines", Comb. Sci. and Tech. 13, 1970. 5. Annand, W.J.D., "A New Computational Model of Combustion in the Spark Ignition Engine", Proc. Inst. Mech. Engrs. 119, 1971. 6. Blizard, N.C. and Keck, J.C., "Experimental and Theoretical Investigation of Turbulent Burning Model for Internal Combustion Engines", SAE Paper No. 740191, 1974. 7. Zeleznik, F.J., "Combustion Modeling in Internal Combustion Engines", Comb. Sci. Tech., 59, 1976. 8. Tabaczynski, R.J., Ferguson, C.R. and Radhakrishnan, K., "A Turbulent Entrainment Model for Spark Ignition Engine Combustion", SAE Paper No. 770647, 1977. 9. Hires, S.D., Tabaczynski, R.J. and Novak, J.M., "The Prediction of Ignition Delay and Combustion Intervals for a Homogeneous Charge, Spark Ignition Engine", SAE Paper No. 780232, 1978. 10. Hiraki, H. and Rife, J.M., "Performance and NOx Model of a Direct Injection Stratified Charge Engine", SAE Paper No. 800050, 1980. 11. Borgnakke, C., Davis, G.C. and Tabaczynski, R.J.,- "Prediction of In-Cylinder Swirl Velocity and Turbulence Intensity for an Open Chamber Cup in Piston Engine", Trans. SAE 4, 1981. 12. Davis, G.C. and Borgnakke, C., "The Effect of In-Cylinder Flow Processes (Swirl, Squish and Turbulence Intensity) on Engine Efficiency - Model Predictions", Trans. SAE 91, 176, 1982. 13. Kono, S., Motooka, H. and Nagao, A., "Prediction of Combustion in Spark Ignition Engine by Simulation Model", Proceedings of Int. Symp. on Diagnostics and Modeling of Combustion in Reciprocating Engines (COMODIA-85), pp. 511, Tokyo, 1985. References. 1 Dimensional Models. Al. Bellan, J.R. and Sirignano, W.A., "A Theory of Turbulent Flame Development and Nitric Oxide Formation in Stratified Charge Internal Combustion Engines", Comb. Sci. and Tech., 8,51, 1973. 10

A2. Singh, T. and Surakomol, K., "Mathematical Modeling of Combustion Process in a Spark-Ignition Engine", SAE Paper No. 790354, 1979. A3. Sorenson, S.C. and Pan, S.S., "A One-Dimensional Combustion Model for a Dual Chamber Stratified Charge Spark Ignition Engine", SAE Paper No. 790355, 1979. References. 2 Dimensional Models. B1. Syed, S.A. and Bracco, F.V., "Further Comparisons of Computed and Measured Divided Chamber Engine Combustion", SAE Paper No. 790247, 1977. B2. Ahmadi-Befrui, B., Gosman, A.D., Lockwood, F.C. and Watkins, A.P., "Multidimensional Calculation of Combustion in an Idealized Homogeneous Charge Engine: A Progress Report", SAE Paper No. 810151, 1981. B3. Gosman, A.D. and Harvey, P.S., "Computer Analysis of Fuel-Air Mixing and Combustion in an Axisymmetric D.I. Diesel", SAE Paper No. 820036, 1982. B4. Diwakar, R., "Direct-Injection Stratified-Charge Engine Computations with Improved Submodels for Turbulence and Wall Heat Transfer", SAE Paper No. 820039, 1982. B5. Carpenter, M.H. and Ramos, J.I., "Modeling of a Gasoline - Injected Two - Stroke Cycle Engine", SAE Paper No. 860167, 1986. References. 3 Dimensional Models. C1. Markatos, N.C. and Mukerjee, T., "Three-Dimensional Computer Analysis of Flow and Combustion in Automotive Internal Combustion Engines", IMACS, 4, 1981. C2. Duggal, V.K., Kuo, T.W., Mukerjee, T., Przekwas, A.J. and Singhal, A.K., "Three-Dimensional Modeling of In-Cylinder Processes in DI Diesel Engines", SAE Paper No. 840227, 1984. 11

Chapter 2 SUBMODEL STRUCTURE AND GEOMETRY The models expressing or explaining the individual terms in the basic continuity and energy equations are considered as submodels. They are arranged, so each submodel contains the analysis and the description of a phenomena that can be regarded as a separate subject. In this way the models can easily be exchanged for better ones, as they become available, without changing the overall model structure. This kind of flexibility enables the model to stay current and be easily upgraded to serve new needs, as they appear. 2.1 Submodel structure The overall model subjects are shown in figure 2.1, where the various submodels correspond to the boxes as indicated. The connection between the subjects are also shown in the diagram to illustrate the dependence between them. The skeleton of the model is built around the continuity and energy equations as presented in chapter 1. The development of these models often leads to the necessity of still more models as new concepts are needed, for example the turbulent combustion model requires the knowledge of the turbulence, which in turn needs information about the mean flow. The presentation of all the submodels are arranged, so the most essential ones are given first and the more advanced given last. The simplest engine simulation model, besides the standard textbook air Otto cycle analysis, would entail the following models: 1.2

GENERAL SUBJECTS IN ENGINE ANALYSIS CHARGE PREPARATION INPUTS COMBUSTION' PRE-MIXED' DESIGN VARIABLES | PRE-MIXED (SI)' FUEL iNJECTION OPERATING VARIABLES j STRATIFIED (SI)' FUEL/AIR MIXING \ FUEL CHARACTERISTICS' / DIESEL * FUEL EVAPORATION | \ / KNOCK CYCLE ANALYSIS FLUID MOTION INDUCTON EMISSIONS 0 INDUCTION SWIRL * t COMPRESSION | HC, SQUISH r IEXPANSION I *CO TURBULENCE 0 EXHAUST * NOX L MANIFOLD DYNAMICS _ / _ __\ _,PARTI CULATES FRICTION PREDICTED OUTPUTS TN R HEAT TRANSFEI * WALL SHEAR / *FUEL ECONOMY CONDUCTION, MECHANICAL' EMISSIONS CONVECTION COMPRESSION ||~TORQUE ||~RADIA CONVECTION 1,COMPRESSION v TORQUE 0r

MINIMUM TWO ZONE MODEL Two zone geometry Thermodynamic properties Combustion model Heat transfer model Valve flow model A model with these building blocks could give a reasonable albeit simple description of engine performance. It would, however, lack a number of important aspects that are essential for both the calculations of engine performance and the emission analysis. The performance is generally acknowledged as being very influenced by the flow field set up in the combustion chamber. Fast burn engines rely for instance on the proper level of turbulence to ensure a high rate of combustion while the performance is then limited by the heat transfer rates to the walls. To predict some of these trade offs it is important to include the submodels for these effects. As an example of some of the models required in a more advanced analysis of the performance considerations the following models are considered: ADVANCED SUBMODELS Mean flow model Turbulence model Manifold flow model Boundary layer zone model A state-of-the-art thermodynamic model for the performance is achieved with these submodels included. From such a model several important design aspects can be analysed and optimized. Since many design constraints are imposed from the emnission standards, it is envisioned that some treatment of the emissions must be included for a realistic model to be usefull. The additional submodels for the analysis of the emissions will be given in a future report and they would include the following: EMISSION SUBMODELS NOx formation model Chemical equilibrium model Knock model Unburned HC model 14

For other engine configurations a number of models may be needed to simulate a direct injected gasoline engine or a diesel engine. These submodels are not planned for in the present program, but could fit in with the overall structure of the formulation. ENGINE CONFIGURATION MODELS Fuel injection model Fuel jet penetration and spread Droplet evaporation Prechamber model Mixing model Stochastic combustion model 2.2 Geometry of Two-Zone Model Two types of geometries must be considered for several of the submodels. The actual combustion chamber shape and the form of the two zones. The combustion chamber shape is important for the mean flow and turbulence in the chamber as well as the flame location and flame travel distances. During the combustion period two zones will fill the combustion chamber volume, one zone with unburned gases and another zone with the burned gases. These zones must be located in space to calculate the interface area and the area of wall contact for each zone as they are needed in the combustion and heat transfer models. The combustion chamber shape is selected as a flat head and a piston with a cylindrical cup or bowl. This is illustrated in Fig. 2.2 and approximates the actual shape, which often have curved surfaces and smaller recess volumes. Regions (1) and (2) are also indicated, as the volume over the piston top outside the radius of the bowl and the volume in and over the bowl, respectively. The following submodels are affected by the geometry with the influence indicated. Only the directly affected models are mentioned, with the other submodels only indirectly sensitive through the primary effects. Flame Front Geometry The geometry is important for several of the submodels, but its prime influence is on the flame travel distances. Assuming the flame propagation is spherical with the spark plug as the center of the sphere the spark plug location determines the longest distance the flame have to travel. In the given geometry the flame front area and the volumes of the unburned and burned gases can be calculated as a function of the spark plug location. In the present version of the model the spark location can be anywhere in the 15

clearence volume and there is a provision for up to two spark plugs. With two spark plugs the locations are restricted to a symmetric placement of the plugs along a diameter of the chamber. The resulting flame geometry will then be symmetric with two propagating flame that eventually meet. Mean Motion The swirling mean motion is considered as an axial symmetric motion and the geometry controls the distribution of the angular momentum. As the total angular momentum is integrated over the volume the clearence height as a function of radius is the important parameter. This also controls the squish motion induced by the piston motion, if there is a variation of the clearence height with radius as the case with a cup in the piston. The total wall surface area and its location with radius has an influence on the friction reducing the swirling motion. Heat Transfer Essentially only the amount of wall areas wetted by the burned gases are important during combustion, with the location less important. In the compression or expansion stroke the total wall area and not its location become a factor, since the charge is assumed uniform. The geometry could be'extended to a more general shape for both the head and the piston. Keeping a prescribed clearence height as a function of radius the calculation needed for the above submodels could be performed. It would also be possible to add volumes such as prechambers and crevice volumes around the spark plug and above the top piston ring. These can easily be handled by the multi zone formulation without introducing major changes to the geometry dependent calculations in the other zones. These extensions are possible future model developments that will be introduced as they become necessary. 16

EXHAUST I N TA K E me <gr i SPARK J TURBULENCE I TD( I ~ SQUISH I0PISTON I -/I Ei SWIRL,'/ 6~~~~~~~~~~~~~~~~~~~~~~~,,!

Chapter 3 THERMODYNAMIC PROPERTIES The thermodynamic properties for the unburned and the burned mixture are needed for the description of the gases. These are related to the mixture composition given by the mass fractions, Ci, or the mole fractions, Yi. From the composition and the thermodynamic state properties (T, P) we can find all the other properties. The methods employed are most accurate to the lean side of stoichiometry, since the rich mixture is chemically more complex. 3.1 The Unburned Mixture The unburned gases are treated with a model developed by Hires et al [1]. In this model the composition, the enthalpy and the gas constant is found from a fit to the properties listed in the JANAF tables. The fuel is considered to be a hydrocarbon with the following chemical composition Fuel: CxHyOzNw Any unburned fuel from residuals or exhaust gas recirculation (EGR) in rich combustion is considered to be decomposed into hydrogen, carbon oxides and nitrogen. The gas mixture is modeled with seven different species as C02, H20, CO, H2, 02, N2, Fuel with the rich and lean regimes treated separately. In the rich mixture the water gas reaction constant is fitted as a function of temperarture, so the composition is temperature sensitive, but not pressure sensitive. The lean mixture on the other hand has a constant composition. The mixture enthalpy, heat capacity and the gas constant are evaluated from the composition and the enthalpy data for each of the seven species. The enthalpy is curve fitted to be in the form 18

h(T) = ao + a1T + a2T2 + a3T3 + a4T4 + aT-1 (3.1) 3.2 The Burned Mixture The combustion products are modeled according to Martin and Heywood [2]. In this model the composition is not explicitly found, instead a correlation with temperature and pressure is found for two classes of dissociations. The split of thri-atomic molecules into diatomic molecules and the split of diatomic molecules into atoms are considered. This leads to two reaction equations with their equilibrium constants fitted as a function of both pressure and temperature. In this model the properties enthalpy, heat capacity, gas constant are found. Since the number of molecules may shift as a function of temperature and pressure, the gas constant, Rb, is a variable. This also implies that the enthalpy is a function of pressure. In the analysis in chapter 1 we needed the rate of change for the gas constant, so this model also finds the derivatives of the gas constant with respect to the temperature and pressure. CRT T (T) (3.2) P aR CRP = R(,;9P) (3.3) With these values, the rate of change of the gas constant can be evaluated from R T P -= CRT + CRP (3 4) which will lead to additional terms in the expressions for the A and B in Eq.(1.13). The derivative of enthalpy with respect to pressure at constant temperature and constant air-fuel ratio can also be calculated from this model. ah CT = ( (3.5) 19

3.3 The laminar flame speed The laminar flame speed is also a property of the gas mixture and it is essential for the turbulent combustions model. Finding the laminar flame speed from a theoretical analysis of the detailed chemestry is nearly impossible. It can formally be done by finding an eigenvalue of the differential equation for the propagating flame as shown by Williams [3]. However, the results are at best within the same order of magnitude as the experimental results. Such an analysis is better suited to show the general form for the dependence of the laminar flame speed on the composition, temperature and pressure. From theoretical considerations and experimental data the flame speed for hydrocarbon fuels have been found as St = C( ) [Y; YO exp(- TR )] (3.6) PO R m This was first done by Van Tiggelen and Deckers [4] and later refined from new measurements by Tabaczynski et. al. [5], Lavoie [6] and Ferguson and Keck [7]. In the above correlation the exponents ac, a and b are selected according to the fuel and the range of air-fuel ratios, a and b essentially indicates the order of the reaction with respect to the fuel and the oxygen; In reality several simultaneous reactions take place so the constants are fitted to data and the most important ones are the frequency factor C and the activation energy E. The temperature is the mean flame temperature, which can be found from the adiabatic flame temperature, Tbi, and the unburned gas temperature, Tu Tm = 0.74(Tb - Tu) + TU (3.7) 3.4 The quench distance The quench distance of the laminar flame is the distance from the wall that the flame will stop its propagation. The thin gas layer next to the wall is therefore not burned by a propagation mechanism, but rather in a subsequent diffusion process. For a long time this effect was thought to be responsible for the unburned hydrocarbon emissions and thus studied extensively. Today the consensus is that the crevice volumes are resposible, so the quench process may not be that important. In this model the quench distance is used for the unsteady heat transfer model and is thus computed. The correlations for the quench distance follows nearly the same procedure as for the laminar flame speed, but since it depends upon a heat loss to the wall the transport property is also important. From the previous work [6,7] a correlation was found as rdqunch = CT )( 2 (3.8) 20

In this form T1 is the lean limit flame temperature and k is the conductivity of the unburned mixture. Since Sl enters, the quench distance is a function of the same parameters as the laminar flame speed and when it is large the quench distance is small. Typically the distance is O.lmm or smaller growing to larger values as the lean flamability limit is reached. 21

REFERENCES 1. Hires, S.D., Ekchian, A., Heywood, J.B., Tabaczynski, R.J. and Wall, J.C., "Performance NOx Emissions Modeling of a Jet Ignition Prechamber Stratified Charge Engine", SAE Paper No.760161. 2. Martin, M.K. and Heywood, J.B., "Approximate Relationships for the Thermodynamic Properties of Hydrocarbon-Air Combustion Products", Comb. Sci. Tech., 15, 1, 1977. 3. Williams, F.A., "Combustion Theory", 2nd. Ed., Benjamin-Cummings, 1985. 4. Van Tiggelen, A. and Deckers, J., "Chain Branching and Flame Propagation", 6th. Int. Symp. on Combustion, 1957. 5. Tabaczynski, R.J., Ferguson, C.R. and Radhakrisnan, K., "A Turbulent Entrainment Model for Spark Ignition Engine Combustion", SAE Paper No. 770647, 1977. 6. Lavoie, G.A., "Correlations of Combustion Data for SI Engine Calculations - Laminar Flame Speed, Quench Distance and Global Reaction Ratesn, SAE Paper No. 780229, 1978. 7. Ferguson, C.R. and Keck, J.C., "On Laminar Flame Quenching and its Application to Spark-Ignition Engines", Comb. and Flame, 28, 197, 1977. 22

Chapter 4 COMBUSTION MODELING 4.1 Introduction Combustion is the single most important process in the engine. It will determine the timing of the heat release and therefore influences all the other phenomena that occur inside and outside the combustion chamber. The proper prediction of the combustion is thus essential for any attempt to forecast effects as the performance and emissions from the engine in detail. Trends such as the variation of efficiency with compression ratio, but not with speed (RPM), can be estimated by very simple models in a broad general sence. Assuming no heat transfer and a uniform gas as air, an instantaneous complete combustion at TDC will give the Otto cycle, and a constant pressure combustion over a finite time will give the Diesel cycle. However, these predictions are not satisfactory for todays sophisticated engines and the developments needed are far beyond such a simplistic description. The rate of combustion is strongly influenced by the flow field in the engine, which is always turbulent rather than laminar under normal running conditions. We are thus faced with the task of modelling the turbulent combustion process, which is vastly different from the purely chemically determined laminar combustion process. If this was not so, it can easily be shown that todays engines would not be able to operate with ordinary fuels since the laminar flames are an order of magnitude slower than the turbulent flames. Also, the engine would not be capable of operating at different RPM's due to the limited variation of the laminar flame speed with changes in the operating conditions such as air-fuel ratio, pressure and temperature. A completely analytical theory of turbulent combustion can only be developed via the extension of the statistical theory of turbulence to include all fluctuating quantities. This will lead to the formulation of a complex theory with conservation equations for turbulent quantities involving various double, triple and higher order correlations. However, the equations contain unknowns which have to be modeled in some way to close the formulation. Presently the theory has not been agreed upon for the unsteady 23

compressible conditions as found in the engine. The precise details of the flame propagation will remain uncertain in the absence of a better understanding and improved measurements of the microstructure of turbulence. In the meantime the practical need for the simulation of the engines is urgent, and to meet it a pragmatic approach is taken by inventing various'models', which will be presented next. 4.2 Combustion Rate Models One way to simulate the combustion in engines is an empirical specification of the combustion rate as a function of time or crank angle. Such functional relationships with the start of combustion, O., and the combustion duration, AOC, as parameters have been proposed for simplicity. Two examples of these functions are the cosine burning law [11] and the Wiebe function [2] as shown below. dXb 0 - 0, dab = 2/sin[lr(0 0d)] (4.1) dO 2AOc AO~ dXb (m + 1) ( - O mexp - (4.2) dO AO — AO - AO ] where Xb is the mass fraction burned, a is an efficiency perameter and m is a slope parameter. The adjustable parameters in these models can b'De selected from measurements by the correlation to pertinent operating conditions. The models and their parameters are typically not closely related to the physics, so it is difficult to express the correlation of these parameters in terms of physical quantities describing the engine conditions. Over small enough ranges of engine operating conditions one can always obtain a reasonable fit and then these models become very convenient tools to study all other engine related processes. In particular for the modeling efforts it is expedient to examine the behavior of the other submodels that influences the combustion, such as the turbulence, under well defined rates of combustion. An example is the study of NO formation done by Blumberg [1]. The first definite ideas about the mode of action of turbulence on a flame were formulated by Damk6hler in a well-known paper [3]. He suggested two possible ways in which turbulence can affect combustion: * The large scale effects, which lead to the development of the burning surface either by distorting it or by breaking the flame up into unconnected islands. Each element of the developed surface burns at the normal flame mechanism. * The small scale effects, which lead to an intensification of heat and mass transfer in the flame as a result of summation of the molecular and turbulent diffusion coefficients. The motion of the small scale eddies is emphasised as the major mechanism in flame propagation. 24

Following these ideas two categories of turbulent flame models have been proposed. One is based on the concept of a distorted wrinkled laminar flame front, this is a classical [4,5] mechanistic model type mainly focusing on the geometry of the flame surface. The other is developed from the progressing theory (actually modeling) of turbulence [6,7,8] assuming an infinite chemical reaction rate or zero chemical time compared to the finite turbulent mixing rate or time. This mixing con.troled combustion approach is more realistic for the engine environment as evident by several experimental investigations [9]. Small scale effects are local in nature, but they are highly influenced by the large scale effects. The turbulence enhances the local flame propagation, but the mean flow structures affect the total flame front area and the flame position within the combustion chamber. The total combustion duration is then determined by the combined effect of propagation speed and the necessary flame travel distances to burn the charge. Since the effect of the turbulence is to increase the combustion speed, the model given by Annand [10] was to simply multiply the laminar flame speed by a factor to obtain the turbulent flame speed, namely: St = CS, (4.3) The factor C, however, is to be determined from experiments and is a function of the engine type and running conditions. In an attempt to unify the model Lucas [II] specified the C as a function of the engine speed: C = 1 + 0.00197 rpm (4.4) However, they found that the model failed to predict the flame travel time variation over the complete equivalence ratio range. This indicates that the turbulent flame speed is not linearly related to the laminar flame speed. Another approach is to model the turbulent flame propagation process as an entrainment and subsquent laminar burnup of the discrete, coherent, turbulent eddies. Blizard and Keck [12] presented a model to calculate the mass burning rate in a spark-ignited engine, which can be discribed by the following equations: dm, d PUAfm (4.5) dmb me - mb (4.6 dt rb where me is the mass entrained by the flame front and mb the mass burned. The entrainment of the fresh charge with density pu takes place over the flame surface Af, with the entrainment velocity u,, followed by the burning of the mass in a characteristic time rb. In the model they assumed that entrained eddies having a characteristic radius I2 will subsequently burn in a constant characteristic time rb = I~/SR, where the SY is the 25

laminar flame speed for the fuel-air mixture. For the hydrocarbon fuel the laminar flame speed has been correlated as: s = C(Po) y E )] (4.7) by Van Tiggelen and Deckers [13]. In this expression the mean flame temperature Tm is: Tm = 0.74(Tb~ - Tu) + Tu (4.8) with the unburned gas temperature Tu and the adiabatic flame temperature Tbo. The power a and b are orders of the reaction with respect to fuel (concentration YF) and oxygen (concentration Yo), respectively. It is unlikely, however, that a laminar flame will be stable under the conditions inside the cylinder. They also assumed a constant entrainment velocity tul over the combustion period. This is not true since the flow will be accelerated due to the expansion of the burned gases. It is not surprising, therefore, that when McCuiston et al [14] used this model in an engine simulation they found that the predicted burning intervals increased with the engine speed contrary to the experimental data. An attempt has been made by Tabaczynski et al [15] to refine the Blizard and Keck model. Firstly, the eddy burning time which was assumed to be constant in the original model is modified to be dependent on the eddy size, laminar flame speed and turblence intensity: 2= 1 ~~Asls 1~~~ ~(4.9) where A is the Taylor microscale, I the eddy integral length scale. From the isotropic relationship the Taylor microscale is given by: A = I () ReT- (4.10) where the constant A has been set equal to the value 1, and the turbulent Reynolds number can be expressed as ReT = VT/V = utl/v. The second modification is that entrainment velocity us is assumed proportional to the turbulence intensity u': ue = St + u' (4.11) This allows for the acceleration of the burned gas if they had the proper turbulence intensity as input. They assumed that initially the integral length scale is proportional to the flow dimension and the turbulence intensity is proportional to the mean inlet velocity. And they used an isotropic rapid distortion assumption resulting in: 26bPO I= ~p- ) (4.12) 25

which gives an increase of tlu' as the unburned gas is compressed rapidly. Experiments, however, show that the turbulence level decreases early in the compression, see, for example, Lancaster [16] and Tabaczynski [17]. One way to describing the dynamic behavior of the turbulence is to use a set of transport equations for turbulent variables such as k E c, k, I, etc. ( For details see the chapter about turbulence modeling ), where k is the turbulent kinetic energy and E is the dissipation rate of k. From the solution of the equations and the dimensional analysis which gives the relation i/u' k/e, the necessary parameters for the modelling can be calculated. The models have been widely used in the past for a variety of shear flows, see Launder and Spalding [18]. The k e models have been most successful. They are fairly accurate and the most versatile models to date though some uncertainties still remain in the modelling of the terms, particularly in the e equation. These models have been used for the two-dimensional flow with swirl by Gosman [19,20], Saffman [21], Ramos [22]. The k e- models have been used for the combustion by Borgnakke et al [23-25] in the lumped form for the compressible unsteady flows without and with swirl as an extension of the incompressible models just mentioned. And they also have been used in the local differential form by Witze [26], Gosman [27], and Diwakar [28]. Although the models have shown excellent agreement with some experimental data, many problems still need to be studied. They are, for example, the question of homogeneity of the turbulence field for complex geometries involving squish and swirl, the use of the spherical flame front geometry, the heat transfer to the wall in a time and space dependent turbulent flow field, the role of the bulk flow in the determining the thickness and burnup of the HC quench layer, the aerodynamics of the intake and the compression process which determines the flow field at which combustion is initiated. From the survey we expect that some significant results can be obtained with considerable less sophistication than will be required to make a detailed multidimensional model. It is not apparent that the present phenomenological approach will be adequate. However, it is equally apparent that for the foreseeable future, the cost and the complexity of the detailed multidimensional computations will render them ineffective in dealing with the kind of parametric investigations typically required in the applied, industrial environment. From this perspective we can see the benefit of the phenomenological models. 27

REFERENCES 1. Blumberg, P.N. and Kummer, J., "Prediction of Nitric Oxide Formation in Spark Ignited Engines - An Analysis of Methods of Control", Comb. Sci. Tech., 4, 73-95, 1971. 2. Taylor, C. F., The Internal Combustion Engine in Theory and Practice, Vol. 2, MIT Press, Cambridge, MA, 1979. 3. Damkohler, G., Zeit. fur Electrochem., 46, 601, 1940. 4. Lewis, B. and Von Elbe, G., "Combustion, Flames and Explosions of Gases", Academic, N.Y., 1961. 5. Williams, F.A., Combustion Theory, Benjamin/Cummings Publ. Co. Inc., 1985. 6. Spalding, D. B., "Development of the Eddy Break-up Model of Turbulent Combustion", 16th Int. Symp. on Combustion, p. 1657, The Combustion Institute, 1976. 7. Libby, P.A., Bray, K.N.C. and Moss, J.B., "Effects of finite Reaction Rate and Molecular Transport in Premixed Turbulent Combustion", Comb. and Flame, 34, 285, 1979. 8. Pope, S.B., "A Monte Carlo Method for the PDF Equations of Turbulent Reactive Flow", Comb. Sci. Tech., 25, 159, 1981. 9. Cushing, B.S., Faucher, J.E., Gandbhir, S. and Shipman, C.W., "Turbulent Mass Transfer and Rates of Combustion in Confined, Turbulent Flames II", 11th. Int. Symp. on Combustion, p. 817, The Combustion Institute, 1967. 10. Annand, W.J.D., Research Notes, J.M.E.S., 12, no. 2, 1970. 11. Lucas, G. G. and James, E. H., "A Computer Simulation of a Spark Ignition Engine", SAE Paper No. 730056, 1973. 12. Blizard, N.C. and Keck, J.C., "Experimental and Theoretical Investigation of Turbulent Burning Model for Internal Combustion Engines", SAE Paper No. 740191, 1974. 13. Van Tiggelen, A. and Deckers, J., "Chain Branching and Flame Propagation", 6th Int. Symp. on Combustion, p. 61, The Combustion Institute, 1957. 14. McCuiston, F. D., Lavoie, G. A. and Kauffman, C. W., "Validation of a Turbulent Flame Propagation Model for a Spark Ignition Engine", SAE Paper No. 770045, 1977. 15. Tabaczynski, R. J., Ferguson, C. R. and Radhakrishnan, K., "A Turbulent Entrainment Model for Spark Ignition Engine Combustion", SAE Paper No. 770647, 1977. 16. Lancaster, D.R., Krieger, R.B., Sorenson, S.C. and Hull, W.L., "Effects of Turbulence on Spark-Ignition Engine Combustion", SAE Paper No. 760160, 1976. 17. Tabaczynski, R.J., "Turbulence and Turbulent Combustion in Spark-Ignition Engines", Prog. Energy Combust. Sci., 2, 143, 1976. 18. Launder, B.E. and Spalding, D.B., Mathematical Models of Turbulence, Academic Press, N.Y., 1972. 19. Gosman, A.D., R.J.R. Johns and A.P. Watkins, "Assessment of a Prediction Method for In-cylinder Processes in Reciprocating Engines", General Motors Research 28

Laboratorie Symposium on Combustion Modeling in Reciprocating Engines, Warren, Mi., Nov. 1978. 20. Gosman, A.D. and Harvey, P.S., "Computer Analysis of Fuel-Air and Combustion in an Axisymmetric D.I. Diesel", SAE Paper No. 820036, 1982. 21 Saffman, P.G., Studies in Applied Math., 53, no. 1, 17, 1974. 22 Ramos, J.I., Gany, A. and Sirignano, W.A., "Study of Turbulence in a Motored Four Stroke IC Engine", AIAA J., 19, 595, 1981. 23. Borgnakke, C., Arpaci, V.S. and Tabaczynski, R.J., "A Model for the Instantaneous Heat Transfer and Turbulent in a Spark Ignition Engine", SAhE paper No. 800287, 1980. 24. Borgnakke, C., Davis, G.C. and Tabaczynski, R.J., "Predictions of In-Cylinder Swirl Velocity and Turbulence Intensity for an Open Chamber Cup in Piston Engine", Trans. SAE, 90, 964, 1981. 25. Davis, G.C. and Borgnakke, C., "The Effect of In-Cylinder Flow Processes (Swirl, Squish and Turbulence Intensity) on Engine Efficiency - Model Predictions", Trans. SAE, 91, 176, 1982. 26 Witze, P.O., "Comparison Between Measurements and Analysis of Fluid Motion in Internal Combustion Engines", Sandia Rep. SAND81-8242, Combustion Facility, Livermore, 1981. 27 Gosman, A.D., Tsui, Y.Y. and Watkins, A.P., "Calculation of Three Dimensional Air Motion in Model Engines", SAE Paper No. 840229, 1984. 28. Diwakar, R., "Direct-Injection Stratified-Charge Engine Computations with Improved Sub-models for Turbulence and Wall Heat Transfer", SAE paper No. 820039, 1982. 29

Chapter 5 HEAT TRANSFER 5.1 Introduction In-cylinder heat transfer models have been developed for more than 20 years with only a modest degree of attention paid to the results. The reason for this lack of interest is that the early simple models provided adequate results for the equally crude engine simulations in which the models had been used. Recent advancements in the other submodels for internal comibustion engines require more accurate information about the heat transfer than previously available. In this chapter we will trace the development of the heat transfer models up to the latest, most accurate models, and try to give some insight into the in-cylinder heat transfer mechanism. Heat transfer takes place in thin layers between gases in the cylinder and the solid walls. These layers are the thermal boundary layers over which the temperature of the gases exibits a large change with the steepest gradient at the wall. Formally the heat transfer is given by this gradient and the thermal conductivity of the gases. However, since the gradient is dependent upon the transfer properties for energy in the outer part of the boundary layer it becomes influ enced by the velocities through the turbulent transport (random convection). The presence of the walls also has a profound influence on the momentum distribution in the gases since the mean velocities and the turbulence must vanish at the walls. The heat transfer at any given location is thus a function of both the thermal boundary layer and the momentum boundary layer. 5.2 Heat transfer models Early model [1-6] developments correlated instantaneous convection heat transfer coefficients in the same form as steady heat transfer correlations found in flows over flat plates or in pipes. In these convection dominated flows the boundary layers are fully developed and in equilibrium with each other, so an analogy between the transport of energy and momentum is valid. The correlations typically take the form: 30

Nu = CRe (5.1) where the Reynolds and Nusselt numbers are defined as: UL Re = (5.2) hL Nu =k (5.3) Here the characteristic length, L, is either taken as the current combustion chamber height or the bore, and the velocity is a function of the mean piston speed. Various assumption are usually made for the calculation of the transport coefficients k and v, since they are dependent upon the gas composition, pressure and temperature. In most analysis they are evaluated at the free stream conditions, but average film temperature or the wall temperature have also been used. The deficiency of the above type of heat transfer models is that they do not account for a number of effects that are present in the engine environment. Some of the most important of these are listed below * Unsteadiness. * Complex mean velocity field. * Compressibility. * Gas composition variations. * Wall deposits. The unsteadiness implies that the energy and the momentum boundary layers are not nescessarily in equilibrium or in a speudo steady state either with themselves or with each other. A strong compression in the engine affects all of the flow related properties and causes the boundary layers to change dimensions (thickness) due to the change in density; Since the velocity field is not as simple and organized as in the pipe flow or flow over a flat plate it cannot force the boundary layers to adjust quickly enough to achieve even near equilibrium conditions. Woschni [4] did improve the basic theory by altering the velocity, at which the Reynolds number is based, for different processes in an attempt to match the actual gas velocities. A corrected value based on the mean piston speed is used to account for scavenging, compression/expansion and combustion. The corrections are given in the following expressions.For the scavenging period the velocity is U = ClUmp, (5.4) 31

and during compression and the combustion phase it is U = C2 Umps + C3 Ucomb (5.5) with the combustion induced velocity Ucomb = VT~(P- Po) (5.6) pV~ The velocity during the scavenging period is then assumed to be different from the velocity in the compression stroke. It is expected to be larger, as the whole charge is driven by the flow through the valves, and this motion will decay with time after the valves are closed. During combustion the expansion of the burning gases induces a motion which he assumes is proportional to the increase in the pressure above the corresponding motored pressure. This additional velocity increases the modeled heat transfer and takes part of the compressibility, due to combustion, into account. The expansion after the completion of combustion will give the opposite effect than this model is predicting and the heat transfer will then be overpredicted. Woschni also included the effect of pressure and temperature on the transport property by writing the heat transfer coefficient as: h = 110 L-'2 p.8 U.8 T-'53 (5.7) The heat transfer data for a particular engine configuration was then, correlated and the calibrated constants were found to be C1 = 6.18; C2 = 2.28; C3 = 3.24 x 10-3 (5.8) These constants are only valid for his particular engine geometry and must be reevaluated for different configurations. After obtaining these constants he varied several operating parameters one by one and measured the heat transfer. His correlations matched the measurements reasonably well. Engine modelers decided that Woschni's correlations matched engine heat transfer data better than most other models and have used them almost exclusively ever since. Some authors [2,6,8] attempt to incorporate radiative effects in their heat transfer models. Woschni states the primary cause for variation in the radiation heat transfer in the internal combustion engine stems from variable emissivities caused by changes in soot concentration and flame volume, not from changing flame temperatures. The variation in soot concentration and flame volume is reasonably accounted for in the combustion velocity correction. Therefore, according to Woschni radiation effects are more accurately accounted for by being absorbed in the combustion velocity constant, instead of the fourth power temperature difference in the Stephan-Boltzmann equation. Annand [2,6] used the basic turbulent pipe flow correlation, but instead of using Woshni's variable velocities for Reynolds number calculations, he first [2] used mean piston speed

and later [6] backed out the velocity from the total kinetic energy calcultions described by Knight [3]. Also unlike Woschni, Annand [8] accounted for radiation using the black body.emmision given by the Stephan-Boltzmann equation. Because of the previously mentioned need for greater accuracy, improved models have recently been developed [11,12,15,16,21-24]. Borgnakke, Arpaci and Tabaczynski [16] have utilized the estimation -, - TW qwj = TgTw (5.9) f T where 8T = thermal boundary layer thickness keff = effective thermal conductivity Here keff and ET are determined from turbulent modeling theory. The thermal boundary layer thickness were determined from an energy equation written for the thermal boundary layer and thus changed with time. The results showed significant improvements in the calculations of instantaneous heat transfer. Others have also incorporated turbulence effects into their models, with the primary differences in their approaches lying in the method of determining the gas velocity used in the correlations. Poulos and Heywood [24] used the following expression; UeSf = [2K+2k1+ 24 ]2 (5.10) Here K and k are the kinetic energy of the mean flow and the turbulence, respectively. The turbulent kinetic energy is estimated by a k - c turbulence model, as shown in chapter 8. Morel et al [23], used the Coulburn anology, I 2 hg = Cf p Ueff Cp P 3 (5.11) where Cf and Cp are coefficients of skin friction and pressure respectively, Uff is the effective gas velocity and hg is the convective heat transfer coefficient used in the convection equation; q = hg(Tg - TW) (5.12) In this case the effective gas velocity is defined by; Uff = [U + U2 + 2k] (5.13) 33

Here U, and Uy are obtained from a mean flow model which accounts for piston speed, swirl and squish, k is again calculated using a k - E turbulence model. The gas temperature, Tg, is estimated from the amount and location of burned gasses, which are determined from a two zone combustion model. The development for the Morel model was specifically for the insulated diesel engine. The insulation drastically reduces convection heat transfer, which increases the significance of radiation. Therefore, radiation was treated completly seperatly in this case. Nikanjam and Greif [13] approached the problem by solving the energy and continuity equations for a rapid compression process. This approach matched experimental data well, but they have not yet expanded this to cover a firing engine or a turbulent flow field. Another numerical effort involving a motored engine [15] solved the axial and radial momentum equations along with the continuity and energy equations. This effort is along the lines of mutidimensional engine simulations which have not yet yielded practical applicable results. 5.3 Present Model Since only a few models have been proposed and they are very similar, three of them will be included in the overall engine simulation program. They can then be used for comparisons in the development of extended models. The models selected for study are the following: 1. Woschni [4] 2. Annand [6] 3. Borgnakke [16] The first two model types are the steady state correlations corrected for the unsteady velocities. The inclusion of the unsteady thermal boundary layer in the third model is expected to improve the prediction, particularly during the rapid compression in the combustion period. Further improvements on the mathematical models are closely dependant on obtaining an improved understanding of in - cylinder heat transfer through extensive experimental work. Recently, Alkidas [18,19] studied the effects of air - fuel ratio, load, engine speed, spark timing, and coolant temperature on engine heat transfer. His conclusions included that engine speed and load strongly influence the heat transfer, while air - fuel ratio and spark timing have less effect. Additionaly, variations of coolant temperature within practical operating limits have no measureable effect. Heywood and Lyford-Pike [20] have used schlieren photography to measure the thermal boundary layer thickness and were able to correlate the data with a simple unsteady boundary layer scaling law 34

since the wall heat transfer determines the wall temperature effects such as deposits will change the conditions. A study by Anderson and Prakash [21] showed the influence of a variable conductivity through the deposit layer. They developed an analysis and modeled the effects through correlations with experimental data. Efforts such as these are vital in improving and verifying future in-cylinder heat transfer models.

NOMENCLATURE h - heat transfer coefficient L - characteristic length k - thermal conductivity of gas k - turbulence kinetic energy K - kinetic energy in mean flow P - pressure Nu - Nusselt number Re - Reynolds number U - velocity v - kinematic viscosity C1 - constants Ump - mean piston speed V - volume Subscripts O- motoring i - m intake valve closure g - gas property mps - mean piston speed comb - combustion eff- effective property 36

REFERENCES 1. Overby, V.D., et al, "Unsteady Heat Transfer in Engines", SAE Trans., vol. 69, p. 461-494, 1961. 2. Annand, J.D., "Heat Transfer in the Cylinders of Reciprocating I.C. Engines", Proc. Inst. Mech. Eng., vol. 177, no. 36, 1963. 3. Knight, B.E., "The Problem of Predicting Heat Transfer in Diesel Engines", Proc. Inst. Mech. Eng., vol. 179, 1964-1965. 4. Woschni, G., "A Universally Applicable Equation for the Instantaneous Heat Transfer Coefficient in the I.C. Engine", SAE Paper no. 670931, 1967. 5. Lefeuvre, T., Meyers, P.S. and Uyehara, O.A., "Experimental Instantaneous Heat Fluxes in a Diesel Engine and Their Correlation", SAE Transactions, vol. 78, Paper no. 690464, 1969. 6. Annand, J.D. and Ma, T.H., "Instantaneous Heat Transfer Rates to the Cylinder Head Surface of a Small C.I. Engine", Proc. Inst. Mech. Eng., vol. 185, 1970-71. 7. Summers, I.G.S, "Convective Heat Transfer in a Rapid Compression Machine Simulating a Spark Ignition Engine", M.Sc. Thesis, University of Manchester, 1970. 8. Hassan, H., "Unsteady Heat Transfer in a Motored I.C. Engine Cylinder", Proc. Instn. Mech. Engr., vol. 185, p. 1139, 1971. 9. Dao, K., Uyehara,*O.A. and Meyers, P.S., "Heat Transfer Rates at Gas-Wall Interfaces in Motored Piston Engines", SAE Paper no. 730632, 1973. 10. Annand, J.D., "Heat Transfer from Flames in Internal Combustion Engines", in Heat Transfer from Flames, N.H. Afgan and J.M. Beer, eds., p. 377, Wiley, New York, 1974. 11. Chong, M.S., Milhers, E.E. and Watson, H.C., "The Prediction of Heat and Mass Transfer During Compression and Expansion in I.C. Engine", SAE Paper no. 760761, 1976. 12. Dent, J.C. and Suliaman, S.J., "Convective and Radiative Heat Transfer in a High Swirl Direct Injection Diesel Engine", SAE Paper no. 770407, 1977. 13. Nikanjam, M.- and Greif, R., "Heat Transfer During Piston Compression", J. Heat Transfer, vol. 100, 1978. 14. Siewart, R.M., "Engine Combustion at Large Bore-to-Stroke Ratios", SAE Paper no. 780968, 1978. 15. Greif, R., Namba, T. and Nikanjam, M., "Heat Transfer during Piston Compression Including Side Wall and Convection Effects", Int. Jour. Heat Mass. Transfer, vol. 22, p. 901, 1979. 16. Borgnakke, C., Arpaci, V.S. and Tabaczynski, R.J., "A Model for the Instantaneous Heat Transfer and Turbulence in a Spark Ignited Engine", SAE Paper no. 800287, 1980. 17. Annand, J.D. and Pinfold, D., "Heat Transfer in the Cylinder of a Motored Reciprocating Engine", SAE Paper no. 800457, 1980. 37

18. Alkidas, A.C. and Myers, J.P., "Transient Heat-Flux Measurements in the Combustion Chamber of a Spark Ignition Engine", Jour. of Heat Transfer, vol. 104, p.62, 1982. 19. Alkidas, A.C., "Thermal Loading of the Cylinder Head of a Spark-Ignition Engine", Heat Transfer Engineering, vol. 3, no. 3-4, p. 66, 1982. 20. Lyford-Pike, E.J. and Heywood, J.B., "Thermal boundary layer thickness in the cylinder of a spark-ignition engine", Int. Jour. Heat Mass Transfer, vol. 27, no. 10, p. 1873, 1984. 21. Anderson, C.L. and Prakash, C., "The Effects of Variable Conductivity on Unsteady Heat Transfer in Deposits", SAE Paper no. 850048, 1985. 22. Wang, C.S. and Berry, G.F., "Heat Transfer in Internal Combustion Engines", ASME Paper no. 85-WA/HT-23, 1985. 23. Morel, T., Blumberg, P.N., Fort, E.F. and Keribar, R., "Methods for Heat Analysis and Temperature Field Analysis of the Insulated Diesel". NASA CR-174783, 1984. 24. Poulos, S.G. and Heywood, J.B., "The Effect of Chamber Geometry on Spark Ignition Engine Combustion", SAE Paper no. 830334, 1983. 38

Chapter 6 VALVE FLOW MODELS 6.1 Introduction Knowledge of inlet and exhaust flow characteristics are essential for the prediction of volumetric efficiency and in-cylinder flow fields. Both of these are of importance for the combustion process. The parameters of greatest interest for the inlet flow are the mass flow rate and swirl generation, while for the exhaust process only the mass flow rate significantly affects subsequent processes. Present models (see any related reference) almost exclusively use isentropic compressible flow theory to predict ideal mass flows and correct these with a flow coefficient to predict actual mass flow rates. Swirl studies, on the other hand, differ on methods of characterization of the swirl generation [1]. The analysis usually lead to dimensionless number which correlates the experimental measurements of angular momentum in the flow to a group of characteristic parameters. 6.2 The Isentropic Mean Flow Model Flow calculations are based on the assumption that the valve assembly behaves like an isentropic nozzle. Valve flow predictions consist of theoretical calculation of ideal nozzle flows combined with an experimental correction obtained from a steady flow test of the valve and port in question. This procedure is outlined below with the continuity, energy and entropy equations written for a control volume around the throat section. The energy equation for a steady state nozzle flow with a constant mass flow rate, rh, is V2 V2 (h + )in + = (h +) (6.1) 2 26 where subscript "in" represents inlet (upstream) conditions and "ex" represents exiting (downstream) conditions. The heat transfer can be neglected due to the relatively small wall area in the throat, so the energy equation states that the stagnation enthalpy 39

V2 hst = h+ 2 (6.2) is constant through the throat section. For an ideal gas the change in enthalpy can be written as Cp times the change in temperature, so the velocity in the smallest area (throat) is given by the inlet velocity and the temperature difference as: Vc2- Vin + 2Cp(Tin - Tez) (6.3) To find the temperature of the gas in the exit we write the entropy equation for the throat section sez Sin + s", (6.4) where s"' is the specific entropy production in the flow. Assuming now that the flow has no friction the process can be regarded as reversible (also no q), so s" is zero. This gives a constant entropy for the flow, which for an ideal gas gives the desired relation to obtain the temperature T, (P.)Z (6.5) Tin Pin Here we have introduced the ratio of specific heats, k, which is calculated from the heat capacity, Cp, and the gas constant, R, in the following manner k = aCp — P- R) (6.6) When the temperature ratio is substituted into the energy equation the result is; k —1 z,e Vin + 2CpTin[l -( ) ] (6.7) Pin Assuming a uniform exit velocity over the area A,,, the ideal mass flow rate becomes r8, = pAezVez,s (6.8) with the exit density related to the inlet density and the pressure ratio we get,me = Pin( P ) zvAeV. (6.9) The ideal mass flow rate can be calculated from the stagnation conditions and the pressure ratio across the nozzle. The ratio between the actual mass flow rate and the ideal flow rate is known as the flow coefficient. Cf — trha (6.10) 40

or m~actual = Cf a (6.11) The above expressions must be modified if the flow through the throat becomes sonic which occurs for a pressure ratio of Pe. 2 kPin < ( 2 + 1 (6.12) where the mass flow rate is maximum and independent upon the pressure ratio. During normal engine operation the flow is only choked for very small valve lifts and for a very short time, see [12]. 6.3 Experimental Measurement It is necessary to carry out steady-flow measurements of the valve and port in question to obtain the flow coefficient. The value of the flow coefficient is obviously dependent on the choice of A,,. Two natural choices exist, the valve curtain area, which is the area from the valve seat to the rim of the valve, or the minimum cross sectional area of the port. The choice is unimportant as long as one is consistent during experimentation and model development and data input. Two main questions can be raised about this model. First is the question of an unsteady process being modeled by steady state experimentation and theory. This was examined by Fukutani and Watanabe [4] and they concluded that the effect is negligible if the gas velocity in the throat exceeds approximately 100 times the average velocity of the valve. This condition is satisfied most of the time and the effect is thus not important under standard engine operating conditions. The second problem is the affect of variable pressure ratio on flow coefficient. It has been observed by Woods and Khan [5] that at higher valve lifts the flow coefficient varied up to 20% for variations of pressure ratio from 8 to 5. This can either be ignored (or compensated for by correcting the constant flow coefficient) or corrected by using a flow coefficient which is both a function of valve lift and pressure ratio. A discussion of some of the operating effects are also presented by Sherman and Blumberg [12]. 6.4 Swirl Quantification Quantification of in cylinder swirl is necessary for prediction of; bulk fluid flows, fuel mixing in injected engines, and combustion rates. Current models utilizing swirl information [6] require knowledge of the bulk angular momentum of the fluid within the cylinder. Actual prediction of this information is still not possible without complicated, not yet practical three dimensional models. Consequently, models must rely on steady 41

state experimental results which determine a swirl coefficient in a similar manner as the flow coefficient is used to predict mass flow rates. This coefficient usually assumes a different form for every piece of experimental work, but the general idea is to compare measured angular momentum to some maximum possible value. Original measurements done by Alcock [7], Fitzgeorge and Allison [8], Jones [9], and Tanabe et al [10] used a paddle wheel inside the cylinder of a steady flow device. This configuration is inadequate because it's sensitivity to paddle locations and shapes cause inconsistent results from flow fields with equal angular momentum but different locations of swirl concentration. Thus, current measurements consist of a flow straightening honeycomb oriented on a steady flow rig such that the discharging flow is purely axial with respect to the cylinder. Consequently the entering flow which originally possesses some angular momentum is left with none. Measurement of the torque on the flow straightener required to achieve this yields the fluid's original angular momentum. This device was introduced by Tippleman [11]. The maximum possible influx of angular momentum has been expressed in several different ways. It is achieved when the flow enters perfectly tangent to the cylinder wall. Uzkan, Borgnakke and Morel [1] described it with the quantity MV;.B, here rh is the 2 mass flow rate, Vi. is the isentropic nozzle velocity, and B is the cylinder bore. This leads to the following definition of the swirl coefficient; 2T C B (6.13) Davis and Kent [6] exchanged B/2 with J,, which is the distance of the valve center from the cylinder axis. This makes their swirl coefficient equal to; T cB Nf= T, _= (6.14) rnVViR, - 2Ru To make use of the swirl and flow information, both coefficients must be available for the entire range of valve lifts used in the simulation. At each discrete point in time, the simulation calculates maximum flow and swirl and determines actual values by multiplying by the respective coefficient. 42

REFERENCES 1. Uzkan, T., Borgnakke, C., and Morel, T., "Characterization of Flow Produced by a High-Swirl Inlet Port", SAE Paper no. 830266, 1983. 2. Borgnakke, C., Davis, G.C. and Tabaczynski, R.J., "Predictions of In-Cylinder Swirl Velocity and Turbulence Intensity for an Open Chamber Cup in Piston Engine", SAE Paper no. 810224, 1981. 3. Davis, G.C. and Borgnakke, C., "The Effect of In-Cylinder Flow Processes (Swirl, Squish and Turbulence Processes (Swirl, Squish and Turbulence Intensity) on Engine Efficiency - Model Predictions", SAE Paper no. 820045, 1982. 4. Fukurni, I. and Watanabe, G., "Air Flow Through Poppet Inlet Valves-Analysis of Static and Dynamic Flow Coefficients", SAE Paper no. 820154, 1982. 5. Woods, W.A. and Khan, S.R., "An Experimental Study of Flow Through Poppet Valves", Proceedings IMechE, 180, Pt. 3N, pp. 32-41 and pp. 53-61 (Discussion), 1965-66. 6. Davis, G.C. and Kent, J.C., "Comparison of Model Calculations and Experimental Measurements of the Bulk Cylinder Flow Processes in a Motored PROCO Engine", SAE Paper 790290, 1979. 7. Alcock, J.F., "Air Swirl in Oil Engines", Proc. Inst. Mech. Engrs., 128, 123-193, 1934. 8. Fitzgeorge, D. and Allison, J.L., "Air Swirl in a Road-Vehicle Diesel Engine", Proc. Instn. Mech. Engrs., (A.D.) No. 4, 151-177, 1962-63. 9. Jones, P.E. "Induction System Development for High - Performance Direct - Injection Engines", Proc. Instn. Mech. Engrs., 180, Pt. 3N, 42 - 52, 1965 - 66. 10. Ohigashi, S., Hamamoto, Y. and Tanabe, S., "A New Method for Measuring Gas Flow Velocity by Electric Discharge", SAE Paper no. 690180, 1969. 11. Tippelmann, G., "A New Method of Investigation of Swirl Ports", SAE Paper no. 770404, 1977. 12. Sherman, R.H. and Blumberg, P.N., "The Influence of Induction and Exhaust Processes on Emissions and Fuel Consumption in the Spark Ignited Engine", SAE Paper no. 770880, 1977. 43

Chapter 7 MEAN FLOW MODELS 7.1 Introduction The mean flow in an actual engine cylinder takes place in an extreme 3 dimensional environment with unsteady inlet and exit flows and sharp changes in the gas velocity over space and time. These conditions make mathematical analysis of the mean flow very difficult. It is not clear from the theory how one can define the split of mean flow and turbulence, since they have time scales of the same order of magnitude. Several investigators have suggested a methodology to make this distinction between the turbulence and the mean flow [1,2,3,4] necessary for the data reduction in measurements. Though this problem is still unanswered, we will review the major methods for modeling the mean flow, where it is then assumed that the flow not included in the model is regarded as turbulence. If the mean flow is found from an analysis with reasonable simplifications that will show the main features of the flow field, it may be possible to examine the influence of geometry and other design variables. The influence of the mean motion is through the following processes * Flame front form * Flame front convection, i.e. burned gas location * Turbulence production, convection * Heat Transfer augmentation * Residual gas content All the above mentioned phenomena do have an important effect on the combustion and the efficiency. It is therefore desirable to predict the mean flow as accuratly as possible without undue complexity, as will be shown next. 44

7.2 Flow Models In the early engine simulation studies the mean flow was not considered at all, since the models for the above mentioned processes did not account for it explicitly. With a flat head and piston crown and symmetric intake configuration, no organized mean flow is set up during the intake stroke. The flow then consists of random vortices shed off the intake valve, decaying into turbulence [5,6,7] with very little mean flow before ignition. The combustion process will then force a slosh motion due to the expansion which depends upon the geometry and location of the flame front. With non flat head or piston surfaces a squish motion is introduced, which can be analyzed by the application of the continuity equation assuming a zonewise uniform density. The modern engine designs do have swirl generating port configurations and semispherical head geometry and in some cases a cup in the piston. These types of geometries do lead to a fairly well organized flow field with a rotation as the predominant motion which is overlaid with secondary vortices and squish and slosh motions [5,6,9,C1-C3,C8C10]. There is thus a need for the analysis of the mean flow to give a better prediction of the flame convection and the generation of turbulence. Several different methods have been used to develop such an analysis and get an understanding of the parameters controlling the flow. The greatest resolution of the mean flow has been obtained by the multi dimensional models that solve the unsteady momentum equations in up to three dimensions. Analysis of the flow have started with the cold flow in motored engines reviewed by Boni [7], Reynolds [8] and Gosman [10]. Very few cases have been reported for three dimensional flow calculations [A1-A5], which were only for the motored engine. Several studies have been made for the axi-symmetric two dimensional flow without combustion [B1-B7] and including combustion [B4]. These models have used various turbulence models ranging from a constant turbulent viscosity to the full Reynolds stress models, ( see some comparisons by El-Tahry [B7] and Ahmadi-Befrui et al [B6] ). The literature shows some of the difficulties involved in getting a proper resolution of the complex flow and the uncertainties in the models of turbulent flows. 7.3 Present Swirl Model The model that will be used in the present investigation belongs to the same class as the main engine simulation model. The combustion chamber is divided into a few control volumes and the flow is analyzed from the conservation of momentum (angular) for these control volumes. The rotating flow can then be calculated from a pseudo one dimensional formulation that includes variation with respect to radius [C5,C61. Based on several measurements in engine configurations that do have swirling ports, we may construct the model with the following assumptions: * axisymmetric flow. 45

* no spatial distribution of temperature and density over chamber. * no turbulent fluctuations in temperature or density. the integral equation for the angular momentum l = fier in the region'j' can be written as: dt= Tj + nifQi - Teif2e (7.1) where Ii= pQjdVj Tj = r.r2dSj (7.2) This equation has been used by Davis and Kent [C3] with V as the total combustion chamber volume. With only one control volume they assumed a velocity profile linear in radius, so it is as a solid body rotation. From the experimental data of Tanabe, Hamamoto and Ohigashi [C8], Johnston, Robinson, Rorke, Smith and Witze [C9] and Witze [C10], it appears that the velocity may be quite different from a solid body rotation. It is therefore desirable to characterize the velocity distribution by a profile with more flexibility than the linear variation with radius. From the work of Borgnakke et al [C5] and Davis and Borgnakke [C6] three simple profiles may be used, each with two parameters. These profiles have the flexibility to curve in the same way that the experimental data indicate and they are: u A ar r <rb (7.3) (ca - P)rb +pr r > rb = ar + fir2 (7.4) = _c +{rr < r+rb (7.5) The two parameters can be solved from two differential equations developed by using the angular momentum equation in two different regions. In these equations the inlet flow terms represents the process that generates the swirling motion inside the combustion chamber. In order to model those terms we use experimental investigations of the intake system, where the intake port and valve are characterized by a flow coefficient and a swirl coefficient. These are explained in more detail in the chapter about the valveflow, but will just for reference be shown here. The angular momentum flux can be correlated from data as done by Tippelman [C10], Novak [C11] and Uzkan et al [C12] in which case the inlet angular momentum can be expressed as: i = NfRV,-Ui (7.6) 46

where the coefficient Nf is a function of valve lift L,, and valve dimension DA. The friction from the wall surfaces are taken from standard steady state friction correlations for flows over flat plates and in channels. This may not be quite appropriate, but nearly no data is available for this confined rotating flow configuration. The velocity given by the profile at the wall location is used as the freestream velocity, U, so the shear stress, rw, becomes w = CfpU,2 Ren (7.7) The constant Cf and the exponent to the Reynolds number are used as 0.06 and -0.2, respectively. The Reynolds number is based on the bore and the velocity Uw, which works quite well, validated by the decaying swirl measured along the axial direction in a steady flow rig [C12]. Along the piston and head surface the wall stress is integrated over the area, since it will be zero at the centerline. The two regions used as the two control volumes are taken as the volume in and over the cup-in-piston and the volume outside as shown in chapter 2 with the geometry. For the case with no cup in the piston the two volumes are taken as equal, giving two concentric cylinders. 47

REFERENCES 1 Rask, R.B., "Comparison of Window, Smoothed-Ensemble and Cycle by Cycle Data Reduction Techniques for Laser Doppler Anemometer Measurements of In-Cylinder Velocity", in Fluid Mechanics of Combustion Systems, Morel, Lohmann and Rackley (Ed.), ASME, 1981. 2 Witze, P.O., Martin, J.K., Borgnakke, C., "Conditionally Sampled Velocity and Turbulence Measurements in a Spark Ignition Engine", Comb. Sci. and Tech., 36, 227, 1984. 3 Martin, J.K., Witze, P.O., Borgnakke, C., "Multiparameter Conditionally Sampled Laser Velocimetry Measurements During Flame Propagation in a Spark Ignition Engine", 20'th Int. Symp. on Comb., The Combustion Institute, 1984. 4 Liou, T.M., Santavicca, D.A., "Cycle Resolved LDV Measurements in a Motored IC Engine", Proc. of the 103rd. Winter Annual Meeting, ASME, 1982. 5 Dent, J.C. and Salama, N.S., "The Measurement of the Turbulence Characteristics in an Internal Combustion Engine", SAE Paper No. 750886, 1975. 6 Namazian, M., Hansen, S., Lyford-Pike, E., Sanchez-Barsse, J., Heywood, J. and Rife, J., "Schlieren Visualization of the Flow and Density Fields in the Cylinder of a Spark-Ignition Engine", SAE Paper No. 800044, 1980. 7 Boni, A.A., "Numerical Simulation of Flame Propagation in Internal Combustion Engines", SAE Paper No. 780316. 8 Reynolds, W.C., "Modeling of Fluid Motions in Engines - An Introductory Overview", in Combustion Modeling in Reciprocating Engines, Eds. Mattavi and Aman, Plenum Press, 1980. 9 Witze, P.O., "Measurements of the Spatial Distribution and Engine Speed Dependence of Turbulent Air Motion in an I.C. Engine", SAE Paper No. 770220, 1977. 10 Gosman, A.D., "Multidimensional Modeling of Cold Flows and Turbulence in Reciprocating Engines", SAE Paper No. 850344, 1985. Three Dimensional model references. Al. Griffin, M.D., Diwakar, R., Anderson, J.D. Jr. and Jones, E., "Computational Fluid Dynamics Applied to Flows in an Internal Combustion Engine", AIAA paper 78-57, 1978. A2. Markatos, N.C. and Mukerjee, T., "Three-Dimensional Computer Analysis of Flow and Combustion in Automotive Internal Combustion Engines", IMACS, Vol. XXIII, No. 4, p. 354, 1981. A3. Duggal, V.K., Kuo, T.W., Mukerjee, T., Przekwas, A.J. and Singhal, A.K., "Three - Dimensional Modeling of In - Cylinder Processes in DI Diesel Engines", SAE Paper No. 840227, 1984. A4. Gosman, A.D., Tsui, Y.Y. and Watkins, A.P., "Calculation of Three Dimensional Air Motion in Model Engines", SAE Paper No. 840229, 1984. 48

A5. Arcoumanis, C., Bicen, A.F., Vafidis, C. and Whitelaw, J.H., "Three - Dimensional Flow Field in Four - Stroke Model Engines", SAE Paper No. 841360, 1984. Two Dimensional model references. B1. Bracco, F.V., "Introducing a New Generation of More Detailed and Informative Combustion Models", Trans. SAE, Vol. 84, 1975. B2. Gosman, A.D. and Watkins, A.P., "A Computer Prediction Method for Turbulent Flow and Heat Transfer in Piston/Cylinder Assemblies", Proc. Symp. on Turbulent Shear Flows, Penn St. Univ., 1977. B3. Gosman, A.D. and Johns, R.J.R., "Development of a Predictive Tool for InCylinder Gas Motion in Engines", SAE Paper No. 780315, 1978. B4. Syed, S.A., Bracco, F.V., "Further Comparisons of Computed and Measured Divided Chamber Engine Combustion", SAE Paper No. 790247, 1979. B5. Ramos, J.I. and Sirignano, W.A., "Axisymmetric Flow Model With and Without Swirl in a Piston-Cylinder Arrangement With Idealized Valve Operation", SAE Paper No. 800284, 1980. B6. Ahmadi-Befrui, B., Gosman, A.D. and Watkins, A.P., "Prediction of In-Cylinder Flow and Turbulence with Three Versions of k - e Turbulence Model and Comparison with Data", Proc. ASME winter meeting, dec. 1984. B7. El-Tahry, S.H., "A Comparison of Three Turbulence Models in Engine - Like Geometries", Proc. of Int. Symp. on Diagnostics and Modeling of Combustion in Reciprocating Engines, p. 203, Tokyo, 1985. Zero Dimensional model references. C1. Fitzgeorge, D. and Allison, J.L., "Air Swirl in a Road Vehicle Diesel Engine", Proc. Instn. Mech. Eng., Vol. 25, 1962. C2. Dent, J.C. and Derham, J.A., "Air Motion in a Four - Stroke Direct Injection Diesel Engine", Proc. Instn. Mech. Engrs., 188, 269, 1974. C3. Davis, G.C. and Kent, J.C., "Comparison of Model Calculations and Experimental Measurements of the Bulk Cylinder Flow Processes in a Motored Proco Engine", SAE Paper No. 790290, 1979. C4.- Kido, H., Wakuri, Y., Ono, S. and Murase, E., "Prediction of In-Cylinder Gas Motion in Engines by an Energy Method", SAE Paper No. 800985, 1980. C5. Borgnakke, C., Davis, G.C. and Tabaczynski, R.J., "Predictions of In-Cylinder Swirl Velocity and Turbulence Intensity for an Open Chamber Cup in Piston Engine", SAE Paper No. 810224, 1981. C6. Davis, G.C. and Borgnakke, C., "The Effect of In - Cylinder Flow Processes (Swirl, Squish and Turbulence Intensity) on Engine Efficiency - Model Predictions", Trans. SAE, Vol. 91, p. 176, 1982. C7. Kono, S, Motooka, H. and Nagao, A., "Prediction of Combustion in Spark Ignition Engine by Simulation Model", Proc. of Int. Symp. on Diagnostics and Modeling of Combustion in Reciprocating Engines, p. 511, Tokyo 1985. 49

C8. Tanabe, S., Hamamoto, Y. and Ohigashi, S., "Swirl in a Four - Stroke Cycle Engine Cylinder", Bull. JSME, 24, No. 152, 287, 1978. C9. Johnston, S.C., Robinson, C.W., Rorke, W.S., Smith, J.R. and Witze, P.O., "Application of Laser Diagnostics to an Injected Engine", SAE Paper No. 790092, 1979. C10. Tippleman, G., "A New Method of Investigation of Swirl Ports", SAE Paper No. 770404, 1977. C11. Novak, J.M., "Simulation of the Breathing Processes and Air-Fuel Distribution Characteristics of Three-Valve, Stratified Charge Engines", SAE Paper No. 770881, 1977. C12. Uzkan, T., Borgnakke, C. and Morel, T., "Characterization of Flow Produced by a High Swirl Inlet Port", SAE Paper no. 830266, 1983. 50

Chapter 8 TURBULENCE MODELING 8.1 Introduction The fluid motion within the cylinder of a piston engine has a major influence on the performance. The general motion within the cylinder and the associated turbulence affect the charge stratification, combustion, heat transfer and scavenging processes. In order to predict the charge stratification, one must understand and be able to predict the turbulent mixing processes. Combustion analysis requires knowledge of the turbulence structure and the scales. The flow in wall boundary layers controls the heat transfer and quenching. Hence, a complete predictive capability requires a great deal of knowledge about the general circulation in the cylinder, the turbulence and the boundary layers. The purpose of this report is to survey the current turbulence models and provide the information of what is known about this phenomena. Assume the instantaneous velocity of the turbulence satisfies the usual continuity and momentum equations, ap apui dp + 0O (8.1) at ax6 Oui aui Oaij Ptu + _3 d ) = 0d (8.2) where oij is the stress tensor. Invoking Stokes hypothesis for an isotropic Newtonian fluid: aui auuj 2auk jk (8.3) ax6 Bai 3 azk where p is the hydrostatic pressure equal to the average trace of the stress tensor. Introducing the classical Reynolds decomposition as is usually done for incompressible flows: 51

ui = iai + ui' Xij = 2ij + a'ii (8.4) the averaged momentum equation is D.i a Dti =; a (i -puuP). (8.5) The appearence of the Reynolds stresses t = -puuj (8.6) gives extra unknowns, which have to be determined to get a closed set of equations. These terms must then be modeled in terms of other known quantities either by an algebraic equation or by a balance equation similar to Eq.(8.5) which is the so-called closure problem of turbulence. The turbulent kinetic energy is defined as the kinetic energy in the fluctuating motions. k - tuui. (8.7) From this definition it is observed that pk is equal to -- of the trace of the Reynolds stress tensor, which is twice the average of the normal stresses. This can be used to define a Reynolds shear stress tensor,rt, with zero trace, giving the deviatoric stresses in much the same way as Stokes, see Eq.(8.3). t 2 %. - -,pk8ij (8.8) The shear stresses defined this way express the departure from isotropy in the stress field including the normal stresses, so a completely isotropic field will have zero shear stresses. For compressible flows a mass weighted averaging as suggested by Le Favre 11] will give nearly the same final form of all equations as in the incompressible case. This may be the explanation for the appearent success of incompressible models applied to weakly compressible flows, where the divergency of the mean velocity is neglected. The internal flow as in the internal combustion engine is one example of a flow with a high compressibility effect, where the divergency of the mean velocity cannot be disregarded. Modeling of this phenomena is not firmly established. Presently there exist several different proposals for changes to the incompressible models to account for compressibility effects. These will be explored in a later section. 52

8.2 Algebraic Models The models of this class use the eddy-viscosity concept of Boussinesq [2], where all fluxes are modeled as a gradient diffusion of the driving property. This concept stems from the recognition that convection due to turbulent motion is radom in nature and behaves much like the diffusion process, such as conduction or viscous diffusion. The difference here is that the approbriate transport coefficient then becomes a flow property rather than a property of the fluid itself. The shear stress therefore becomes proportional to the turbulent viscosity and the mean velocity gradient as: T =-pu'v = (8.9) In this type of model only the turbulent viscosity, as a scalar quantity, is modeled. Relating the viscosity to a velocity scale, u', and a length scale, 1, it is most commonly written as vt = Cu'l (8.10) The two most familiar examples of this class are distinguished by the way in which Vt is calculated from length and velocity scales of the flow. * The constant eddy viscosity model. Trubchikov and Prandtl [3] have proposed, for free jet flows where the length scale I is usually taken as proportional to the width of the flow 6, the formula: Vt = C6S(maz - fimin) (8.11) All the quantities on the right hand side, except C, may be functions of longitudinal distance x. vt is supposed to be uniform over any cross section. * The celebrated mixing length hypothesis of Prandtl [4,5] is stated as: t = 121 a l (8.12) The length scale is then modeled as a fraction of the flow dimensions. This model is still the basis of many engineering calculations of the turbulent boundary layer which are carried out today. For wall-bounded flows, Prandtl introduced the linear relation I = icy, where ic - 0.41 denotes the experimentally determined Von Karman [5,6] length constant. He based this on the arguments that the mixing length must vanish at the wall and eddy sizes are controlled by their distance from the wall. In the outer region of the boundary layer (y/S > 0.2), or in free flows, I does not continue to increase linearly, but remains 53

approximately constant: I = AS. The values of the proportionality constant A are shown below: Flows Mixing Layer Plane Jet Radial Jet Round Jet A = 1/6 0.07 0.09 0.13 0.075 For many boundary layer flows, Prandtl's mixing length hypothesis works surprisingly well. However, the proportionality constants vary markedly from one case to another (see above table). This lack of universal applicability is an indication that the model lacks some important features of real flows. Actually, the mixing length hypothesis implies that generation and dissipation of turbulent energy are in balance everywhere, so the convection and diffusion of turbulent energy are ignored. Only the employment of differential equations for the turbulent quantities can overcome these restrictions. 8.3 One-Equation Models Most one-equation models make use of the equation for the turbulent kinetic energy per unit mass, k. An equation governing k can be developed by multiplying the instantaneous momentum equation for ui with ui, substracting from this the equation formed by multiplying the average momentum equation with uii, and averaging the result. This manipulation with the momentum equations will give the turbulent kinetic energy, i.e. k, equation in the following form 1( )=a(-plu + 2 4- p d- ) a i_ 2yss (8.13) where sij is the fluctuating rate of strain: Iau' au'. =i= 2(J i + a ) (8.14) Ng and Spalding [7] amongst others interpreted the equation as: Change = Diffusion + Production - Dissipation leading to a formulation of the model equation with three different types of terms on the right hand side Dk P- = Jk + Pk-PDk (8.15) The. total diffusion term is the combination of the laminar and turbulent diffusion fluxes, where the turbulent diffusion term can be modeled using the gradient flux in Eq.(8.9) with k as the driving property: 54

where ak is a turbulent Prandtl number for k. The production term Pk comes from the interaction between the Reynolds stresses and velocity gradients t &ai 2 Dp t'9 ai Pk-,t O =3 Dt +'r a'. (8.17) If the stress field is isotropic, so the shear stresses are zero, we see the influence of the dilatation alone. In an incompressible flow the first term on the right will vanish and the trace of the velocity gradient will be zero. With only small deviations in the normal stresses, so their influence can be neglected, the second term will only contain off diagonal contributions. This is usually a good approximation, but the anisotropy in the normal stresses may have to be accounted for under severe anisotropic or non-typical strain fields, as we will look at later. The dissipation term Dk = e is modeled from the assumption that the rate of the energy dissipation is controlled not by the viscosity, but by the rate at which the large eddies feed energy to the small dissipative scales, which in turn adjust in size to handle this energy. Thus, if I represents the integral length scale, e should scale on k and 1: e = CD I (8.18) Now with the various terms modeled, the k equation may be written as: Dk a Mt Ok r aui 2 _ (8.19) P~yjj~ =-Gu+-)-+r3s- + -k — PcD T8.19) P Dt dza Uk axj ax, 3 Dt I The coefficients CD = 0.09,crk = 1.0 are experimentally determined from the equilibrum turbulence, P = E, see Launder [3]. For the turbulent viscosity Vt, Kolmogorov [8], Prandtl [3,5] and Emmons [9] independently proposed the relation: Vt = C,k21 (8.20) where k is calculated from the differential equation so that vt becomes dependent upon the turbulent energy transport and 1 is modeled by an algebraic equation proportional to the flow dimension. If I is taken as the integral length scale the constant, Ca, can be chosen equal to unity without loss of generality. The mixing length which is typically used in these models differ from I by a constant as I = Clmiz (8.21) 55

Wieghardt [9], Glushko and Bushnell et al [11] applied this model to various boundary layer flows. Wolfshtein [12] used the model to calculate the impinging jet and Spalding [13] applied it to separated flows. The quality of the prediction depends, of course, on the adequacy of the expression for 1, but whenever the diffusion and convection of turbulence play a significant role, the superiority of this model over the algebraic models has been demonstrated. medskip There are other one-equation models. Nee and Kovasznay [14] postulated a transport equation for vt, which also needs an empirical algebraic expression for 1. Townsend [15] proposed the assumption that the ratio of uY.u to k might be a constant. Bradshaw et al [16] employed this assumption successfully in predictions for several uniform density flows, but the assumption is not realistic, because often the shear stress vanishes at locations where k remains finite. 8.4 Two-Equation Models The one equation models discussed in previous section have necessitated the use of an algebraic length scale specification. This specification must vary with the boundary conditions. There is little hope of achieving a universal model until I is itself determined from an independent differential equation. Kolmogorov [8] made the first suggestion of I this kind when he introduced a differential equation for w = k /1. This model was later redeveloped and extended by Saffman [17] in the analysis of several different effects not previously looked at. Rodi and Spalding [18] and Ng and Spalding [7] used the k kl model for a number of problems. Rotta [19] derived an equation for I from the NavierStokes equations. The ke E model, which was first proposed by Harlow and Nakayama [20] and also has appeared in the papers of Jones and Launder [21], has been explored and tested on engineering problems more than any other two equation models. It should be mentioned that all two-equation models of turbulence, where the variables are scalars, can be transformed mathematically from one to another. This can be done by the relation of the two scalar variables to the length scale and the velocity scale, where all models to date have used the k equation as one of the two equations. Such a transformation yields the same form except for minor variations in the diffusion term and the application of boundary conditions. It thus appears that the differences are more philosophical than substantial and we will here restrict the discussion to the k E c model, where the isotropic dissipation E is defined by: aul aui EV = 8; u (8.22) 8xj aXj With this model all turbulent properties are expressed in terms of k and c, in particular we have the length and time scales shown as: 56

-I =-C;- r-=- (8.23) A transport equation for e is obtained by the operation v azi on the momentum equation. The result for the incompressible flow is DE a k ax1 x p,+'u, ) u', auau'. au a De -- - (VU I Irt S + _ k ) _ 2Vt'- I — + _ Dt ark ka lj axj p ai asi ax1 oxj asi a k axk a', a au,.' a'I 2 -2v it k a 2(v' ) (8.24) axk axj axj axkax; which may be modeled into the standard form of a balance equation: p = J, + P, - pDe (8.25) The production term P, is often modeled from the expression similar to Pk by a simple scaling with E/k: P6 =' C l (t.a) + P (8.26) where the first term produces the proper behavior of homogeneous shear and the second term PS = E- t (see next section) modifies the behavior for homogeneous isotropic dilation. The diffusion term is modeled similarly to the term in the k equation using a Prandtl/Smith number for the dissipation: a /.It CaE Any <S +e) + j (8.27) The destruction term De is modeled as a decay term proportional to e divided by the turbulent time scale e2 De = CE2 (8.28) Summarizing the above discussion about each of the terms in the modeled equations they can be written: Dk a,t1 ak t. ail 2 Dp pDt = =-(~ + -)k + - i + k DP C (8.29) Dt axj k ax1 ax1 3 Dt(8.29) De a tl a r EaUi 4Dp E2 Pt as= a- ( + - ) + Cls 3 eDt k 57

where the turbulent viscosity is modeled as: k2 z = c — (8.31) The numerical values of constants are determined from experimental data. First C., is found by equating production and dissipation. The constants C 2, a, are obtained from experiments by considering the decay of homogeneous turbulence in the absence of production and diffusion. Then taking k = Ct-n and ignoring convection and viscous effects gives the relation: ]~2 C -= Ce2 - (8.32) The suggested model constants are: |k oa C,| Ce1 Ce2 1.0 1.3 0.09 1.45 1.9 The k - e models have been very successful and most popular in engine modeling. They are fairly accurate and the most versatile models to date though some uncertainties still remain in the modelling of the terms, particularly in the e equation. 8.5 Compressibility When the flow is compressible, so the divergence of the velocity does not vanish, the turbulence models must be modified accordingly. The turbulent kinetic energy equation needs the least modifications if the effect on diffusion is neglected. The productions term Pk is exact, given the mean flow gradients and the Reynolds stresses. An additional term is needed to accout for the changes in e produced by dilation of the flow. Reynolds [22] made an argument about this term. Consider the rapid spherical expansion of homogeneous isotropic turbulence, at the start of which puVu' p pkSHi (8.33) For a spherical expansion with - = rr V 1 a( r2) _ 3r (8.34).2 ar Then from the continuity equation, treating p = p(t), dp =_kV u = -3rp (8.35) 58

Consider only the production term linear in the mean strain, the k equation is dk 2 = — kV * u = -2kr (8.36) dt 3 In rapid distortion the angular momentum of the turbulence should be conserved. Hence, the product of the turbulence length and velocity scales should remain fixed during the expansion, which can be expressed from k and e as k2 — = const. (8.37) Thus the e equation for this case should be de e dk = 2 —-= -4Er (8.38) dt k dt This suggests that, under rapid distortion, P* should be: 4 1 4 c4Dp P =.-_V = Dp (8.39)' ~ 3 3p Dt Hence, in order to obtain the proper behavior of a compressible flow we should add a term Pe = 4DP to the P, term. Other corrections in the equations have been proposed by Morel and Mansour [23], but the effect was found to be small. Several investigations of the compressibility effect have been persued using a rapid distortion approximation. If the production term is much larger than the dissipation the latter may be neglected and the turbulence follows the production. These models have been examined by Townsend [24] for the general case and Hunt [25] presented the special case of an isotropic compression. 8.6 The Present Model 8.6.1 The General description Since the actual mean flow pattern in an engine is very complicated, as shown in chapter 8, it would also make the turbulence field complex. So the best spatial description of the turbulence is obtained when it is analyzed with the same resolution as the mean flow field. Notice, that because the turbulence is produced by the mean flow, recall Eqs.(8.17), it is not possible to describe it with better spatial resolution than the mean velocity field. The turbulence model is thus constructed by making some reasonable simplifications that at least will capture the most important processes in the flow field. We can then develop an understanding of the main features of the flow field and their influence on the combustion efficiency through design changes. The present model follows the formulation and outline from the model by Borgnakke [26]. This model describes- the turbulence by 59

an averaged turbulent kinetic energy model combined with a zonal model for the mean flow. The important aspects to be included in the present model are: * Swirling motion * Turbulence * Unsteadiness * Squish * Compressibility The chamber geometry is assumed to have a flat head with a symmetric cup in the piston. The whole chamber is divided into two regions. Region (1) is the volume over the piston outside the cup radius which is also called the squish region. Region (2) is the volume in and over the cup. This is an idealization of chamber shapes for most practical engine combustion systems. 8.6.2 Integrated Turbulence Model As previously indicated it is assumed that the turbulence field is uniform over the combustion chamber volume. This assumption is partially supported by a number of investigations as reviewed by Tabazynski [27], and Morse, Whitelaw and Yianneskis [28]. So that the spatial effects are not intended to be included in the present formulation. For the present mean flow model there is only one velocity gradient giving the shear stress: = pvt( - ) (8.40) After integration the total derivative will give the ordinary derivative and the convection associated with the mass flow rates in or out of the volume. The diffusion terms become boundary fluxes and are important for the boundary layers next to the walles, so the resulting equations can be written as: d-i 2',1e 27f,2k dk =vt le ) ne _ Vus + (ki - k) (8.41) dt V J o r r 3p dt m dE 4e dp e2 mi deig C eVt | dae _ dV + 6 dp_ C62 + m -) (8.42) dt L dr r 3pdt c +-(-) (8.42) with the turbulent viscosity vt given as: k2 vt - CD- (8.43) 60

The contribution to the turbulence from the inlet flow can be formulated as: 3 ki = (CkUi) =CD CL (8.44) based on the inlet velocity Ui and the inlet clearance L, with the adjustable constants Ck and C1. The above correlation is from the standard assumption, see Launder and Spalding [3], that the turbulent kinetic energy is proportional to the mean flow kinetic energy and the length scale of the turbulence is proportional to the flow dimension, which in this case is the valve lift. The solution of the integrated turbulence model combined with mean flow model will provide the necessary input such as turbulence intensity and turbulence length scales for the combustion, heat transfer and other submodels. Since the model does incorporate information about the inlet conditions and the geometry it will give an improved prediction of the engine performance and the influence of design parameters. 61

REFERENCES 1. Rubesin, M.W. and Rose, W.C., "The Turbulent Mean-flow, Reynolds-stress and Heat-flux Equations in Mass-averaged Dependent Variables", NASA TMX 62.248, March 1973. 2. Boussinesq, J., 1877, Mem. Pres. Acad. Sci. XXIII, 46, Paris. 3. Launder, B. E. & Spalding, D. B., 1972, Mathematical Models of Turbulence, Academic Press, N.Y. 4. Tennekes, H. and J. L. Lumley, A First Course in Turbulence, MIT Press, Cambridge, MA, 1972. 5. Schlichting, H., Boundary Layer Theory, McGraw Hill, 1984. 6. Van Driest, E. R., J. Aerosp. Sci, 23, 1007 (1956). 7. Ng, K.H. and Spalding, D.B., Phys. Fluids, 15, 20, 1972. 8. Kolmogorov, A.N., 1942, IV. Akad. Nauk. SSSR, Ser. Fiz. 56, 2 (see also Imperial College, Mech. Eng. Dept. Rep. ON/6 1968). 9. Hinze, J.O., Turbulence, McGraw Hill, 1975. 10. Bradshaw, P., Cebeci, T. Whitelaw, J. H., 1981, Engineering Calculation Methods for Turbulent Flow, Academic Press, New York. 11. Bushnell, D.M., Carey, A.M. and Holley, B.B., AIAA J., 13, 1119, 1975. 12. Wolfshtein, M., Imperial Coll., HTR section REP SF/R/2, 1967. 13. Spalding, D.B., J. Fluid Mech., 27, 97, 1967. 14. Nee, V.W. and Kovasznay, L.S.G., Phys. of Fluids, 12, 473, 1969. 15. Townsend, A.A., Jour. Fluid Mech., 26, 255, 1966. 16. Bradshaw, P., Ferriss, D.H. and Atwell, N.P., Jour. Fluid Mech., 28, 593, 1967. 17. Saffman, P.G., "Results of a Two Equation Model for Turbulent Flows and Development of a Relaxation Stress Model for Application to Straining and Rotating Flows", Proc. Heat transfer and Fluid Mech. Inst., pp. 191-231, 1978. 18. Rodi, W. and Spalding, D.B., Warme und Stoffubertragung, _, 85, 1970. 19. Rotta, J.C., Proc. 1st Symp. on Turbulent Shear Flows, Penn State Univ., 1977. 20. Harlow, F.M. and Nakayama, P.I., Phys. of Fluids, 10, 2323, 1967. 21. Jones, W.P. and Launder, B.E., Int. J. Heat Mass Transfer, 15, 301, 1972. 22.- Reynolds, W. C., "Modeling of Fluid. Motions in Engines - An Introductory Overview", in "Combustion Modeling in Reciprocating Engines", (Mattavi, J.N. and Amann, C.A., eds.), Plenum, N.Y., (1980). 23. Morel, T. and Mansour, N.N., "Modeling of Turbulence in Internal Combustion Engines", SAE Paper No. 820040, 1982. 24. Townsend, A.A., The Structure of Turbulent Shear Flow, Cambridge Press., 1976. 25. Hunt, J.C.R., "A Review of the Theory of Rapidly Distorted Turbulent Flows and Its Applications", XIII Biennial Fluid Dynamics Symp., Poland, 1977. 62

26. Borgnakke, C., Arpaci, V.S. and Tabaczynski, R.J., "A Model for the Instantaneous Heat Transfer and Turbulent in a Spark Ignition Engine", SAE paper 800287 (1980). 27. Tabaczynski, R.J., "Turbulence and Turbulent Combustion in Spark Ignition Engines", Prog. Energy Comb. Sci., 2, 143, 1976. 28. Morse, A.P., Whitelaw, J.H. and Yianneskis, M., J. Fluid Eng., 101, 208, 1979. 63

Chapter 9 MANIFOLD FLOW MODELS 9.1 Introduction The study of flow dynamics and gas exchange processes in an internal combustion engine lead to better design of the exhaust and inlet manifolds, turbo charger and catalytic converter. These design improvements serve to increase volumetric efficiency and create in-cylinder conditions favorable for complete efficient combustion. The development of the theoretical methods used to predict the gas exchange processes will be briefly reviewed here. This overview will be followed by a somewhat more detailed description of the Helmholtz resonator model. As early as 1964, Benson [1] used the method of characteristics to solve a variety of pipe flow problems which had direct application to internal combustion engine manifold systems. The method of charateristics is a first order finite difference technique and as an explicit method it allows direct calculation of dependent quantities. Subsequently, this investigation was extended to a real engine [2,3] and other investigations [4,51 provided a good understanding of the predictability of this method. Even though the method of characteristics was satisfactory in many cases, the improved accuracy of second order and higher explicit methods was desired. It was found in 1975 [4] and 1982 [6] that the two step second order Lax-Wendoff scheme was simpler to apply, required less CPU time and gave better agreement with experimental results. J.D. Ledger [7] considered several different finite difference schemes and found that the second order implicit centered difference method was most promising. An implicit solution scheme is characterized by the solution of simultaneous equations by an iterative method. This leads to good stability and the allowance of large grid and time step sizes. Currently third and fourth order finite difference methods are available and are very accurate, but high computation times restrict their use to applications requiring high precision. In 1979 Chapman [8] introduced the Dynamic Eularian, Lagrangian, Triangular Algorithm (Delta-4). The scheme approximates a solution to the equations of continuum 64

in two spacial dimensions and time. In the solution procedure the domain of the inte. gration is broken up into a number of triangular elements and finite difference equations are constructed to approximate the solutions to the governing equations. This method is very accurate, however it requires excessive CPU time to execute. In the same year Murthy [9] developed an implicit, finite element scheme using third and fourth order Runge Kutta numerical integration techniques to yield stable solutions. The major disadvantage of solving the equations of fluid motion in an implicit manner using any of the above mentioned techniques is the large amount of cpu time required to execute the solution code. 9.2 The Proposed Model With the need for a reasonably accurate, low cost gas exchange predictor model in mind, the recent developments by Chapman of the Helmholz resonator model warrent extended attention [10]. It is essentially an electrical analog of the gas flow within the manifold (exhaust or intake). A synopsis of the published article is listed below. A Helmholz resonator model can give an indication of where tuning peaks might occur in an engine speed range. It cannot predict volumetric efficiency or other engine parameters without being linked to a complete simulation code. The equations developed become non-linear through the incorporation of a wall friction term in the governing equation based on fully developed turbulent pipe flow. The simplicity of the model allows efficient machine computation for use with a full power train dynmaics simulation. The present model considers the inlet and exhaust system of a multicylinder internal combustion engine as a driven system of coupled non-linear Helmholz resonators and results in a coupled set of second order ordinary differential equations which may be solved efficiently on a digital computer. The resulting equations are a long wavelength approximation to the quasi-one dimensional gas dynamics equations. With the inclusion of wall friction effects, the equation of motion for the mass of air in the duct is given by: c2A2X piAL + iF(z) + Pa = (t) (9.1) where x = displacement of air in the duct L = length of duct Pa = ambient density c = accustic velocity V _ volume of the chamber F(x) = wall friction 65

A = duct cross sectional area The driving function g(t) is zero if the tube is open to the atmosphere and different from -zero if a driving pressure difference is present. The wall friction term F(i) is developed by assuming fully developed pipe flow with a standard correlation for the friction w = _CfP2 (9.2) Cf = CwRe-'25 (9.3) where r= wall shear C1 = skin fraction coefficiient Re = instantaneous duct Reynolds number based on diameter C,= cosntant proportional to the wall surface roughnessss Therefore the duct wall friction F(x) as a function of velocity is given by: F(z) = Cf Ap'2l (9.4) Aw = wetted wall area of the duct (9.5) The second order differential equation is integrated with the same integrator that is used for all the other rate equations, which are of first order. The equation is therefore written as two first order differential equations by using the velocity, U dx=U (9.6) and Eq.(9.1) is then determining the velocity dU c2A2x pAL dt = -UF(U) - Pa + g(t). (9.7) Though the model was presented as a Helmholtz resonator, it should be realized that Eq.(9.1) is the same as a momentum equation in integral form with the runner as a control volume. Now an arrangement of these simple resonators coupled with proper boundary conditions can be assembled to predict pressures in the intake and exhaust system. The boundary conditions are obtained implicitly from the valve flow models and explicitly from the atmospheric conditions. Knowledge of the pressure versus crank angle curve allows the location of a particular design's tuning peak and with the valve flow model can predict volumetric efficiency. 66

REFERENCES 1. Benson, R.S., Garg, R.D. and Woullatt, D.,'A numerical solution of unstead flow problems", International Journal of Mechanical Sciences, Vol. 6, No. 1, 1964. 2. Benson, R.S. and Galloway, "Gas exchange process in a two stroke engine", Proc. I. Mech., Vol. 183, 1964. 3. Blair, G.P. and Goulbourn, J.R.,'A pressure time history in the exhaust system of a reciprocating I.C. Engine", SAE Transactions 1968. 4. MacLaren, J.F.T., Tramschek, A.B., Sanjines, A. and Pastrana, O.F., "Comparisons of the numerical solutions of the unsteady flow equations applied to reciprocating compressor systems", Mechanical Engineering Science, 1975. 5. Low, S.C., Barush, P.C. and Winterbone, D.E., "Transportation of liquid fuel droplets in the pulsative air flow within the S.I. engine intake manifold", SAE paper 810497. 6. Takizawa, M., Uno, T., Oue, T. and Yura, T., "A study of gas exchange process simulation of an. automotive multi cylinder internal combustion engine", SAE paper 820410. 7. Ledger, J.D., "A finite difference approach for solving the gas dynamics in an engine exhaust". Mechanical Engineering Science, Vol. 17, #3, pp.125-132, 1975. 8. Chapman, M., "Two dimensional numerical simulation of inlet manifold flow in a four cylinder internal combustion engine", SAE paper 790244. 9. Murthy, B.S., Lakshminaranan, P.A., Janakiraman, P.A. and Gajendra Babu, M.K., "Prediction of gas exchange process in a single cylinder internal combustion engine", SAE paper No. 790359, 1979. 10. Chapman, M., Novak, J.A. and Stein, R.A., "A nonlinear accoustic model of inlet and exhaust flow in a multi-cylinder internal combustion engine", ASME paper. 67

UNIVERSITY OF MICHIGAN //1/Il///11//l/l//////1///l/// 11///I~//////ll//l /l11111111 111111 3 9015 02653 6006 THE UNIVERSITY OF MICHIGAN DATE DUE 1hi, g'q'1 K U,, 2D/ X< 7i) <.~# j-/ f.".'2t I ( C)