*************************************************************
    *                 A first draft of                          *
    *                                                           *
    *    THE MULTI-FLUID MODEL OF TURBULENT FLOW, MIXING AND    *
    *                                                           *
    *                    COMBUSTION                             *
    *                                                           *
    *                        by                                 *
    *                                                           *
    *              PROFESSOR BRIAN SPALDING                     *
    *                                                           *
    *                        of                                 *
    *                                                           *
    *         Concentration, Heat and Momentum Limited          *
    *                                                           *
    *                      May 1996                             *
    *************************************************************

                         ABSTRACT

    An outline is given of the physical concepts and mathematical
    formulations of a multi-fluid model (MFM) of turbulence.

    Similarities to, and differences from, earlier turbulence models
    are described.

    References to recent works on MFM are provided, but with only brief
    summaries of their contents, the aim of the present paper being to
    collect and explain all the relevant ideas, without distraction by
    particular examples.

    Consideration is given to the advantages likely to be derived from
    applying the MFM to practical problems of engineering and
    environmental science.

    Although MFM is already capable of being applied to many practical
    problems, especially those involving combustion processes or other
    chemical reactions, research and development work can be foreseen as
    likely to improve its acceptability, its realism and its economy. A
    classified list of such tasks is provided.

    --------------------------------------------------------------------

    CONTENTS

    1. An introduction to the multi-fluid-model concept
    1.1 Some origins of the multi-fluid idea
    1.2 Quantification of the multi-fluid concept (pdf's and fpd's)
    1.3 An example of a multi-fluid gas mixture
    1.4 The concept of dependent-variable discretization
    1.5 Discretization possibilities
    1.6 An example of a two-dimensional population
    1.7 Interactions between fluids
    1.8 Coupling and splitting
    1.9 Influences of hydrodynamic factors
    1.10 Why combustion modellers need MFM

    2. Mathematical  formulation
    2.1 Differential equation for the conservation of a fluid
    2.2 The source-sink term
    2.3 The differential equation for a within-fluid property
    2.4 Population-average quantities
    2.5 Computational aspects

    3. Recent work on MFM
    3.1 Reports and publications
    3.2 Unpublished work

    4. Relations to other turbulence models
    4.1 The four main classes of turbulence model
    4.2 Advantages and disadvantages of the four classes
    4.3 Similarities and differences between MFM and previous models
    4.4 Generalising the eddy-break-up model

    5. Application prospects
    5.1 The possibility of immediate use
    5.2 Applications to chemically-reacting flows
    5.3 Applications to hydrodynamic and heat-transfer phenomena

    6. Research and development tasks
    6.1 Mathematical and conputational tasks
    6.2 Comparisons with experiment
    6.3 Conceptual developments


    7. References

    --------------------------------------------------------------------

    1. An introduction to the multi-fluid-model concept
    ---------------------------------------------------

    1.1 Some origins of the multi-fluid idea
    ----------------------------------------

    The fragmentary, intermittent and near-random nature of turbulent
    flows, as a qualitative concept, has been long been observed and
    described qualitatively; yet to express this observation adequately
    in useful quantitative terms is still very difficult.

    In striving to do so, many authors have found it helpful to
    regard a turbulently-flowing system as consisting of an assemblage
    of fragments of many distinct fluids, which intermingle, exchange
    heat, mass and momentum, and possibly merge.

    For example, Reynolds (1874) explained the transfer of heat to
    a surface in terms of the rate at which a hotter fluid moved towards
    that surface while a cooler fluid moved in the opposite direction;
    and Prandtl (1925) conceived of the transfer of momentum across (ie
    the shear stress within) a turbulent jet or boundary layer in much
    the same way.

    Generalising Reynolds' concept, one might postulate a population
    of fluids distinguished by their temperatures.

    To generalise Prandtl's concept, one would imagine a population of
    which the distinguishing attribute was the main-flow-direction
    component of velocity.

    [More information on the origin and history of the multi-fluid idea
    can be found in section 4.1c below.]


    1.2 Quantification of the multi-fluid concept (pdf's and fpd's)
    ---------------------------------------------------------------

    A quantitative way of expressing the variability of a fluid
    attribute is by way of its probability-density function
    (abbreviated to pdf), the shape of which can be deduced from
    measurements of instantaneous local values by suitable processing
    of their variations with time.

    The counterpart of the pdf, in the multi-fluid approach, is the
    "fluid-population distribution" (fpd)

    The multi-fluid concept implies that the pdfs of any fluid
    attribute should be plotted as a series of spikes (ie delta-
    functions), corresponding to the fact that the population-
    defining space is discretized.

    Fdp's are more commonly plotted as histograms, such as the
    following diagram, in which the amount of each fluid in the
    population element is represented by the area of a rectangle
    standing on a base centered on the value of the attribute
    which defines the element:
             __
            |  |                     __
     amount |  |                  __|  |
          ^ |  |         __      |  |  |__          __
          | |  |      __|  |__ __|  |  |  |        |  |
          | |  |     |  |  |  |  |  |  |  |__      |  |
          | |  |   __|  |  |  |  |  |  |  |  |     |  |
            |  |  |  |  |  |  |  |  |  |  |  |     |  |
            |  |__|  |  |  |  |  |  |  |  |  |__   |  |
            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
            |  |  |  |  |  |  |  |  |  |  |  |  |__|  |
            |__|__|__|__|__|__|__|__|__|__|__|__|__|__|
                           ----------> attribute


    The fpd of Reynolds' and Prandtl's two-fluid models would have
    to be represented as:               __
                                       |  |
                                       |  |__
                                       |  |  |
                                       |  |  |
                                       |  |  |
                                       |  |  |
                                       |  |  |
                                       |__|__|
    or rather as:
                                        _____
                                       |  |  |
                                       |  |  |
                                       |  |  |
                                       |  |  |
                                       |  |  |
                                       |  |  |
                                       |  |  |
                                       |__|__|

    because neither author indicated that the two fluids might be
    present in unequal proportions.


    1.3 An example of a multi-fluid gas mixture
    -------------------------------------------

    If a turbulent flame results from the injection of a combustible gas
    into an air atmosphere, and the characterising attribute of the
    imaginary population is the fuel mass fraction, a multi-fluid
    model of the gas in the flame could be a population in which there
    intermingle fluid fragments consisting of, for example:

    - pure air;
    - pure fuel;
    - gas having (say) 10% fuel and 90% air;
    - gas having (say) 20% fuel and 80% air;
    - gas having (say) 30% fuel and 70% air;
      and so on.

    With the above numbers the model could be called an eleven-fluid
    one, with uniform sub-division of the fuel-fraction dimension.

    In reality, concentrations can take any value between zero and
    unity; but the discretization principle entails that attention is
    focussed on a finite number of discrete concentration values.
    Mixtures which in fact have concentrations lying between the two
    nearest discrete values must then be represented by interpolation.

    The material which carries one of the selected concentrations is
    regarded as a distinct fluid, and will be referred to as such below.


    1.4 The concept of dependent-variable discretization
    ----------------------------------------------------

    The example brings immediately into notice the arbitrariness of the
    multi-fluid concept; for three un-forced decisions have already been
    made, namely in respect of:-

    (1) choice of distinguishing attributes, (here fuel fraction,
        rather than temperature or velocity);

    (2) choice of number of fluids (here 11, whereas 110, 1100 or
        11000 might have been chosen);

    (3) choice of attribute difference between the fluids (here a
        uniform 10%, whereas one might have preferred to cluster all 9
        of the "impure" fluids closely around the stoichiometric fuel-
        air value.

    The arbitrariness of the discretization may be an obstacle to
    acceptance of the multi-fluid concept by purists, who may not at
    first recognise that the same objection can be raised to numerical
    simulations of flames in which GEOMETRICAL SPACE is discretized and
    TIME is also sub-divided into a finite number of intervals.

    This analogy can be instructive as well as persuasive; for it
    reminds doubters that grid-independence tests can be conducted for
    fluid-attribute "grids" as well as for geometric ones.


    1.5 Discretization possibilities
    --------------------------------

    The analogy  may also stimulate researchers to employ their skills
    in "adaptive gridding" so as to make the numerical computation of a
    fluid-population distribution both accurate and economical.

    Further, just as grids which sub-divide geometrical space can be
    one, two- or three-dimensional, so can those which are employed to
    distinguish the fluids in an imaginary multi-fluid population.

    Temperature, velocity and fuel-fraction have already been mentioned;
    and of course the number of usable distinguishing attributes is
    unlimited.

    Turbulence energy, vorticity, fragment size, degree of reactedness
    and smoke content are a few of the other attributes which spring
    immediately to mind.

    However, possibility and necessity are not the same: just as there
    is no need to use a multi-fluid model at all, so is there no
    necessity to discretize any more of the fluid-attribute variables
    than the modeller believes will improve flow-simulation capability.

    This is an important observation; for it distinguishes multi-fluid
    models from probabilistic models (see, for example, Pope, 1982),
    which use Monte-Carlo methods for computing the pdf in a multi-
    dimensional space.

    [More is written about the possible varieties of population grid in
    section 6.1 below.]


    1.6 An example of a two-dimensional population
    ----------------------------------------------

    Many turbulent combustion problems involve the separate entry of
    fuel- and oxidant-bearing streams. Gas-turbine combustors, for
    example, exhibit this feature.

    Because fluctuations of both fuel-air-ratio and reactedness may be
    expected to be significant, consideration of a two-dimensional fluid
    population is appropriate.

    The following diagram represents such a population schematically,
    as used for the numerical simulation of a turbulent Bunsen-burner
    flame (Spalding 1995b).

            _________________________________________
            |///////|///////|///////|       |       |
    f - mfu |///////|///////|///////|  16   |  20   |  * mfu stands for
       ^    | inaccessible  |_______|_______|_______|    mass fraction
       |    | fluid states  |       |       |       |    of unburned
       |    |///////|///////|  11   |  15   |  19   |    fuel;
       |    |_______|_______|_______|_______|_______|  * fluid 1 is
       |    |///////|       |       |       |       |    fuel-free air;
            |///////|   6   |  10   |  14   |  18   |  * fluids 13-16
            |_______|_______|_______|_______|_______|    are stoichio-
            |       |       |       |       |       |    metric;
            |  1    |  5    |  9    |  13   |  17   |  * fluid 17 is the
            |_______|_______|_______|_______|_______|    fuel-rich entry
                    --------> mixture fraction, f        stream.

    The vertical dimension (f - mfu, ie the mass proportion of fuel-
    stream-derived material in the local mixture, less the mass
    proportion of unburned fuel) measures the extent of the reactedness
    of the gas, being at its maximum (namely f) when no unburned fuel
    remains, and at its minimum (namely 0) when all the fuel is
    unburned.


    Other examples of two-dimensional populations might involve the
    attributes:-

    * temperature and salinity (ie salt content), for use in simulations
      of sea-water flows in estuaries;

    * along-wall velocity and normal-to-wall velocity in a boundary
      layer; and

    * temperature and upward-direction velocity in a buoyancy-driven-
      flow simulation.

    For example, a 2-dimensional concentration space, concerned with two
    mixture components A and B, might focus attention on 10 distinct
    values of the concentrations of each component. The whole set of
    relevant mixture-composition possibilities would then be represented
    by the 100 pairs of A and B concentration values.

    Once again, real mixtures possessing intermediate compositions would
    be represented by interpolation between the nearest 4 of the 100
    allowed possibilities.

    Three-, four- and more-dimensional models can of course be
    envisaged. There is no limit.


    1.7 Processes affecting the fluid-population distribution
    ---------------------------------------------------------

    (a) The atmosphere as a suggestive example
    ------------------------------------------

    Clouds in the sky may be regarded as easily recognisable
    manifestations of the multi-fluid character of the Earth's
    atmosphere, with water-vapour content as the distinguishing
    attribute.

    Observation of such clouds reveals that:-

    * they are convected across the sky by the wind, more or less at a
      uniform velocity;

    * they move upwards and downwards at what may sometimes be differing
      velocities; and

    * they may come into existence, grow in size, then diminish, and
      ultimately disappear.


    Similar renarks can be made about smoke plumes from chimneys, and
    about the much larger plumes which are created by forest fires. The
    processes so revealed must be expected to take place in all multi-
    fluid populations.

    (b) Convection and diffusion
    ----------------------------

    About the first-mentioned (uniform-convection) process, there is
    little to be said, for it takes place whether the atmosphere is
    turbulent or not.

    As to the second, two modes of relative motion between members of
    the fluid population can be distinguished, namely:-

    - turbulent diffusion, as a consequence of the gradients of
      concentration; and

    - "sifting convection", which is the relative motion according to
      which, by reason of gravity, lighter fluids move upwards more
      rapidly (or downwards more slowly) than heavier ones.

    The speed of relative motion can be reasonably supposed to be the
    outcome of the balance between the buoyancy-force difference (which
    increases the speed) and friction between the different-velocity
    fluids (which decreases it).

    Thus friction between fluids in relative motion is one of the
    interactions to be accounted for.

    Heat transfer and mass transfer between the fluids, by reason of
    their differences in temperature and composition, are also certain
    to take place to some extent.

    (c) Sources and sinks of individual fluids
    ------------------------------------------

    Finally, the creation, growth, diminution and disappearance of
    clouds is a reminder that a multi-fluid model of turbulence must
    make provision, since mass, energy and concentration must obey their
    relevant conservation laws, for the transfer of mass from one
    element of the population into another.

    Just as a human population which is decribed by its age-distribution
    (eg 50% under 30, and 25% over 60) or political allegiances
    (eg 40% left, 35% right, 20% centre and the rest don't-know) will
    change as its members grow older or change their opinions, so will a
    fluid-population-distribution change.

    Heating will effect this if temperature is one of the distinguishing
    attributes; momentum sources will do so if the attributes are
    velocities; and chemical reaction will bring it about if one
    attribute is the mass concentration of one of the reactants.

    Sources and sinks of the above kind can take place in a single
    fluid, on its own. However, multi-fluid turbulence modelling
    involves consideration also of sources and sinks which result from
    the interactions BETWEEN fluids. These are of such central
    importance as to need a section to themselves, which now follows.


    1.8 Coupling and splitting
    --------------------------
    (a) The general idea
    --------------------

    The present author, (Spalding, 1995a) has used the somewhat
    anthropomorphic term "coupling and splitting" to describe the
    interactions between fluids.

    This term suggests that two fluid fragments making contact act as
    parents; and that, after the contact, there exist (what is left of)
    the parents and a collection of offspring.

    The coupling process is conceived as the approach of two fragments
    of unlike fluids, followed by their partial coalescence and mixing.

    The splitting process is their subsequent break-up into new
    fragments which possess concentrations (or temperatures, velocities,
    etc), which are intermediate in value between those of the parent
    fluids.

    Although more sophisticated hypotheses may prove to be more in
    accord with reality, it is the simplest which will be described
    here, namely that which has been called "promiscuous" and
    "Mendelian".

    "Promiscuity" implies that any fluid will couple with any other, ie
    indiscriminately, at a rate depending only on its availability.

    "Mendelian" implies that the offspring may possess the character-
    istics of either parent, in any proportion.

    The following diagram may serve to convey this idea.


    ^  frequency in
    |  population
    |                    father                                  mother
    |           ______ /                                        /
    |          |      |          promiscuous              ******
    |        __|______|__ <------- coupling --------->   ********
    |           | .. |                                   *| .. |*
    |           |*  *|                                   *| -- |*
    |           /-**--/      ____  Mendelian  ______      /-----/
    |         / |/////| /   |      splitting        |   / |_____| /
    |           |/////|     v                       v     /     /
    |           |  |  |    _     _    _    _    _    _  /_________/
    |           |  |  |   | |   | |  | |  | |  | |  | |   |  |  |
    |           |  |  |   | |   | |  | |  | |  | |  | |   |  |  |
    -------------------------------------------------------------------
                        fluid attribute -------------------------->

    (b) The underlying mechanisms envisaged
    ---------------------------------------

    The way in which this redistribution of material between fluids is
    effected is imagined to be as follows:

    (1) Two fragments of fluid are brought into temporary contact by the
        random turbulent motion, thus:

                     ______________________________
                     xxxxxxxxxxxxxxx...............
                     xxxxxxxxxxxxxxx...............
                     xxxxxxxxxxxxxxx...............
                     xxxxxxxxxxxxxxx...............
                     xxxxxxxxxxxxxxx...............
                     xxxxxxxxxxxxxxx...............
                     ------------------------------


    (2) Molecular and smaller-scale turbulent mixing proceses cause
        intermingling to occur, with the result that, after some time,
        the distribution of x (father) and . (mother) material within
        the coupling fragments appears as:

                     ______________________________
                     xxxxx.xxxxx.x.x.x..x.....x....
                     xxxxxxxx.x.x.x.x.x..x.........
                     xxxxxx.xx.x.x.x.x.x..x..x.....
                     xxxx.xx.xxxx.x.x.x...x.x......
                     xxxxxxxxx.xxx.x.x.x.x.........
                     xxxxxxxx.x.x.x.x.x.x..x.......
                     ------------------------------

    (3) Before the intermingling is complete, however, the larger-scale
        random motions cause the fragments to be plucked away again, with
        the result that the amounts of the material having the
        compositions of the parent fluids (pure x, and pure .) are
        diminished, while some fluid  material of intermediate
        composition has been created, as shown below.

             pure   <------------ intermediate --------->   pure
             ____   ____   ____   ___   ___   ____   ____   ____
             xxxx   x.xx   xxx.   x.x   .x.   .x..   ...x   ....
             xxxx   xxxx   .x.x   .x.   x.x   ..x.   ....   ....
             xxxx   xx.x   x.x.   x.x   .x.   x..x   ..x.   ....
             xxxx   .xx.   xxxx   .x.   x.x   ...x   .x..   ....
             xxxx   xxxx   x.xx   x.x   .x.   x.x.   ....   ....
             xxxx   xxxx   .x.x   .x.   x.x   .x..   x...   ....
             ----   ----   ----   ---   ---   ----   ----   ----


    (c) Coupling and splitting for 2D populations
    ---------------------------------------------

    If the population is a two-dimensional one, as illustrated in the
    diagram of section 1.6 above, the promiscuous-Mendelian hypothesis
    would imply that the offspring of the coupling between parents such
    as F and M below would be shared among the cells enclosing some part
    of the straight line joining F to M, ie those containing asterisks
    in the diagram below.

            |_______|_______|_______|_______|_______|_______|_______|
            |       |       |       |       |       |       |  M    |
            |       |       |       |       |       |       |   *   |
            |_______|_______|_______|_______|_______|_______*_______|
            |       |       |       |       |       |   *   |       |
            |       |       |       |       |       *       |       |
            |_______|_______|_______|_______|___*___|_______|_______|
            |       |       |       |       *       |       |       |
            |       |       |       |   *   |       |       |       |
            |_______|_______|_______*_______|_______|_______|_______|
            |       |       |   *   |       |       |       |       |
            |       |       *       |       |       |       |       |
            |_______|___* __|_______|_______|_______|_______|_______|
            |  F    *       |       |       |       |       |       |
            |   *   |       |       |       |       |       |       |
            |_______|_______|_______|_______|_______|_______|_______|


    (d) Alternative nomenclature and hypotheses
    -------------------------------------------

    The fluid-interaction process can also be usefully called the
    "micro-mixing process", so as to distinguish it from the "macro-
    mixing" ones of convection and turbulent diffusion, which re-
    distribute fluids in space without however increasing the intimacy
    of their contact.

    The corresponding term employed by the users of probabilistic models
    (eg Pope 1982) is simply "the mixing model".

    Alternative fluid-interaction processes, which remain to be explored
    are:-

    * parental bias; and

    * differentially-diffusive.

    These will be discussed below, in section 6.3, as will their
    probabilistic-model counterparts.



    1.9 Influences of hydrodynamic factors
    --------------------------------------

    Before the coupling-and-splitting idea can be used in a quantitative
    method of turbulent-flow prediction, some means has to be provided
    for determining its rate. This must be connected in some way with
    quantitative properties of the hydrodynamic flow, for example the
    sizes of the fluid fragments and their relative velocities.

    One of the simplest hypotheses, which has its origin in the so-
    called "eddy-break-up" formula (Spalding,1971a) [see section 4.4
    below], is that the rate is proportional to:

        the mean fluid density,
          times
        the product of the mass fractions of the two colliding fluids,
          times one or other of:
        the mean velocity-gradient
          or
        the energy-dissipation-rate / turbulence-energy

    The first of the two alternatives may be appropriate when the
    hydrodynamic calculations are based on the Prandtl (1925) mixing-
    length theory; and the second may be more suitable when some variant
    of the "k-epsilon" turbulence model (Harlow and Nakayama, 1968) is
    employed.

    Such presumptions are however merely starting points, proposed
    because they provide links with conventional turbulence models.

    The MFM allows insights which suggest more sophisticated formulae,
    especially those which recognise the influences on relative motion
    of the pressure gradient, which may be taken as being the same for
    all fluids in the mixture, and the density, which often differs
    greatly from fluid to fluid.

    The Reynolds number is of course a further probable source of
    influence, as are the "Densitometric Froude Number" and other
    quantities expressing body-force action.

    These matters will be further discussed in section 6 below.


    1.10 Why combustion modellers need MFM
    --------------------------------------

    The following section is inserted at the present point in order to
    provide motivation, at least for combustion specialists, to proceed
    further. It has been adaptted from an earlier publication (Spalding,
    1995a).

    (a) The practical background
    ----------------------------

    Laminar flames are, in principle, easy to simulate by numerical
    calculations; for the equations governing them are known, and
    chemical-kineticists have provided data on the rates of the
    (admittedly very numerous) chemical reactions which take part. To
    simulate a laminar flame therefore, all that is needed is a very
    powerful computer, equipped with the appropriate software.

    However, most flames of industrial and environmental importance,
    whether in furnaces, engines, oil-platform explosions, or forest
    conflagrations, are turbulent.

    (b) The theoretical difficulty
    ------------------------------

    Turbulent flames, by contrast with laminar ones, present
    difficulties of principle which have not so far been adequately
    resolved.  They result from the facts that:

    (1)  turbulent fluctuations of temperature and concentration are of
    too small a scale (in both space and time) to be resolved by any
    currently available computer, so that local-average values of
    temperature and concentration must be used;

    (2)  the relation between the required local-average reaction
    rates and the local-average temperatures and concentrations
    is therefore not what it would be in the absence of fluctuations,
    ie in the circumstances for which chemical kineticists have reported
    their findings.

    In order to understand this mathematically, suppose that a reaction
    between species A and B proceeds at a rate proportional to their
    instantaneous concentrations [A] and [B].

    In the absence of fluctuations, the rate is thus:

                    constant * [A] * [B]

    However, if (to take an extreme view) fluctuations cause [A] to be
    finite only when [B] is zero, and vice versa, the time-average
    reaction rate, which is proportional to the time-average product of
    [A] and [B], is zero.

    The time-average reaction rate is even harder to predict when the
    effects of temperature and its fluctuations are brought into
    consideration; for reaction rates are non-linearly dependent on
    temperature, rising steeply with it.  Thus, if the rate is
    proportional to:

            T ** n,

    where T stands for absolute temperature and n is an exponent between
    6 and 10, it is easy to work out that

            0.5 * (T1 ** n + T2 ** n)

    has a very different magnitude from (0.5 * (T1 + T2) ) ** n.

    The relation between reaction rate and temperature is in fact more
    more complex than a simple power law. The following diagram shows
    the typical shape. Its rise and fall result from the fact that the
    increase of temperature is accompanied by decreases in the amounts
    of available reactants; and, without them, reaction cannot proceed.

------------------------------------------------------------------------
              Variation of reaction rate, oxygen content
              and unburned-fuel content for a sub-stoichiometric
              fuel-air mixture


          ^
          |*      oxygen---> x               # #
          |    *                x           #   #
          |        *               x       #     #
          |           *               x   #       #
     rate |               *              x         #<--reaction
          | unburned fuel---> *        #    x       #  rate
          |                       *   #        x     #
          |                         # *           x
          |                       #       *          x#
          |                     #             *         x
          |                  #                    *     #
          |           #                               *  #
          ------------------------------------------------*
          0    reactedness ( = (T - Tu)/(Tb - Tu)  )      1
------------------------------------------------------------------------

    (c) A homely illustration
    -------------------------

    The point can also be explained and understood as follows:

    when the man is a night-worker and his wife a day-worker, their
    contribution to the population explosion is hard to estimate.

    Being present at the same place AND at the same time is as
    important to interactions between chemical species as it is to
    those between humans.

    (d) Conclusion
    --------------

    The multi-fluid model has been created, in part, to enable
    turbulent-flow phenomena involving chemical ractions to be
    understood, and then simulated numerically.

    However, that is far from being its only purpose; for it can,
    potentially, provide simulations of a pure hydrodynamic or heat-
    transfer character which take better account of the physics than
    do conventional models.

    --------------------------------------------------------------------
    2. MATHEMATICAL FORMULATION
    ---------------------------
    2.1 Differential equation for the conservation of a fluid
    ---------------------------------------------------------

    Let the mass fraction of a specific fluid in the local population
    of fluids be denoted by ma, where a represents the attribute, or set
    of attributes, which characterise this fluid.

    Then the usual considerations of time-dependence, convection,
    diffusion and source-sink action, lead to the following differential
    equations for ma:

           ma_t + V.grad ma = v.div (G.grad ma) + madot    (1)

    where: ma_t  stands for the differential of ma with respect to time,
           V     is the velocity vector of the mixture,
           v     is its specified volume,
           G     is the turbulent exchange coefficient and
           madot is the net rate of creation of ma per unit volume.

    In this equation, V and G are taken as having the same values for
    all fluids of the population.

    If however different values are appropriate for different fluids
    (as, for example, where one is heavier than the other and therefore
    moves at a different velocity), such fluid-specific features will
    be expressed by modification of the source-sink terms, madot.



    2.2 The source-sink term
    ------------------------
    (a) The various kinds of source and sink
    ----------------------------------------

    Apart from the just-mentioned effects associated with velocity and
    diffusion-coefficient disparities, the phenomena to be expressed by
    madot are those which cause the fluid with attribute a to be created
    or destroyed. These may include the following:

        (1) inter-fluid mixing;
        (2) chemical reaction;
        (3) within-fluid sources and sinks.

    These will now be discussed in turn.

    (b) (1) Exchange of mass between fluids
    ---------------------------------------

    The individual fluid fragments which constitute the postulated
    population can be imagined to come into contact with each other in
    a random manner; and, while in contact, to interact by transfer of
    mass, momentum and energy.

    Suppose, for example that the distinguishing attribute chosen for a
    population is temperature T, and that the discrete values of
    temperature which are chosen to describe the population are: Ta, Tb,
    Tc,....... Tz, where: Ta < Tb < Tc <.....< Tz.

    Let the fluids be referred to by name as: A, B, C,.....Z.

    Consider now the coupling of fragments of A with those of, say, E.
    It is probable on physical grounds that, after the contact is over,
    some fluid at temperature Ta and some at temperature Te will still
    exist; but fluids at the intermediate temperatures Tb, Tc and Td are
    also likely to have come into existence.


    Thus, the contact will have reduced the values of ma and me,
    because the parents have given up some of their mass; but the
    values of mb , mc and md will have been increased, and in such
    amounts that the sum ma + mb + mc + md + me is unchanged.

    In order that the source-sink terms ma, mb, mc etc can be computed,
    mathematical expressions must be devised to dictate:

    - at what rate fragment-types a and e "couple" with each other; and

    - how the offspring which they produce are distributed along the
      fluid attribute (ie temperature) scale.


    (c) (2) Chemical-reaction sources
    ---------------------------------
    Let it now be supposed that the fluids in question can sustain
    an exothermic chemical reaction, so that fluid C, for example, would
    tend (if it were allowed), to increase its temperature at the rate
    Tc degrees per second.

    It is of course NOT allowed, because the discretised nature of the
    multi-fluid model does not permit any fluid to exist with
    temperatures other than Ta, Tb, Tc, Td, etc.

    How can this difficulty, which appears to be rather serious for a
    model which purports to be able to represent combustion processes,
    be resolved?

    The answer is: by regarding the heating process as reducing the mass
    fraction of C and increasing that of D. This may be expressed
    quantitatively as follows. If the temperature rise at rate Tcdot
    continued for long enough, the temperature of the fluid would reach
    the next-allowed temperature Td.

    How long is "long enough"? The answer is: (Td-Tc)/Tcdot.

    The heating process can therefore be represented as causing a sink of
    C, at the rate:

         mcdot = - mc . Tcdot/(Td-Tc)       (2)

    simultaneous with a source of D given by

         mddot = mc . Tcdot/(Td-Tc)       (3)


    (d) (3) within-fluid sources and sinks
    --------------------------------------
    The same technique can be used for any other source which tends
    to change the state of fluid in the discretized-attribute direction.

    Such sources include:

    - if the attribute is temperature, the emission or absorption of
      radiation;

    - if the attribute is velocity, the resolutes of pressure gradient
      and body forces.

    In all cases the mdot expressions are of the form:

      source / attribute difference.


    2.3 The differential equation for a within-fluid property
    ---------------------------------------------------------

    The temperature-distinguished fluids just discussed can, so far as
    the multi-fluid concept is concerned, each have their own
    composition, and other fluid properties such as velocities, etc.

    These too need to be computed; and a differential equation is
    available for each.

    If Fb stands for one of the properties of fluid b, the equation is:

      (mb.Fb)t + V.grad (mb.Fb) = v.div (G.grad (mb.Fb) + Sfb (4)

    where b is used in place of a, so as to ease later discussion.

    The terms on the left-hand side of the equation require little
    comment; for they amount to saying that the amount of property F
    per unit mass of fluid b, namely Fb, remains unchanged along a
    particle track in space and time.

    The first term on the right-hand side can be best interpreted by
    decomposing the term inside the bracket as follows:

        grad (mb.Fb) = mb . grad Fb + Fb . grad mb             (5)

    The first term on the right stands for the within-fluid-b
    diffusion of the property F, while the second term implies that
    any of fluid b which changes location as a result of random
    turbulent motion carries its F property with it.

    The term SFb requires more discussion. It will ordinarily consist of
    two parts; the first corresponding to the within-fluid source (for
    example the rate of creation of oxides of nitrogen appropriate to
    the temperature of the fluid and to its chemical composition), and
    the second to the consequences of the coupling-and-splitting
    process.

    In the latter regard, the simplest hypothesis is that the
    "offspring" possess the arithmetic - mean F values of their parents.

    Thus, if the source of B (say) resulting from the couplings between
    A and D is mb(AD), the coupling-splitting contribution to Fb
    consists of terms such as:

      mb(AD)(Fa+Fd)/2

    Of course, since B can be formed from the coupling of A with any
    item in the alphabet from C to Z, the number of possible terms may
    be large.

    There is something else to be said: whereas couplings between A and
    B can be left out of account when the production of offspring is
    concerned, because there are no allowable values of temperature
    between TA and TB, such couplings can produce sources and sinks of
    Fb.

    Their magnitudes will be proportional to (TA-TB), and to the
    coupling rate.


    2.4 Population-average quantities
    ---------------------------------
    The population-average value of a discretized variable, such as
    temperature in the above example, is given by:

     Taz = Sum az , ma.Ta      (6)

    That of a non-discretized variable, such as F in the above example,
    is given by:

      F az = Sum az , ma.Fa      (7)

    Of course, whereas it is possible to compute the population-average
    quantities from knowledge of the properties of all the individual
    fluids, the reverse is not true.

    Many other population-characterizing quantities can be computed once
    the full fluid-population distribution (ie the set of ma's) and the
    properties of each fluid (ie the set  of Fa's) are known. Of
    especial interest are the "moments", ie such quantities as:

      T  az(n)  =  Sum  az  (maTan) (8)

    These are useful when it is desired to transfer information from one
    fluid-attribute grid to another, for example as part of a computer-
    time-economising adaptive-grid procedure.



    2.5  Computational aspects
    --------------------------
    Many computer programs are available which solve, by finite-volume
    or finite-element techniques, differential-equation sets such as (1)
    or (4).

    In order to use these, it is necessary only to provide sufficient
    storage space for (what may become) an unusually large number of
    dependent variables (the ma's and  Fb's), and also to provide coding
    for the source terms.

    The amount of computation increases more than linearly with the
    number of variables, because the number of coupling-splitting
    events to be considered may increase as N*(N-1), where N is the
    number of fluids.

    The word "may" is used, because it will often be possible to find
    physically- or numerically-based reasons for not considering all the
    theoretically possible events.

    It is therefore probable that use of multi-fluid models for
    practical turbulent--flow computations will often cause computer
    time to become a matter of concern.

    Lest it be thought, however, that it will be as burdensome as is
    that associated with the Monte-Carlo approach referred to above, the
    following considerations should be weighed:

    (1) Experience has shown that even a very coarse discretization of
    fluid-attribute space can produce considerable increases in realism.

    Thus the 25-year-old "eddy-break-up model" (Spalding,1972) which is
    a two-fluid model in the present nomenclature, is still in use
    because it is much more realistic than would be the use of a one-
    fluid model (ie complete neglect of the fluctuations), which is what
    it replaced.

    Further, as has recently been demonstrated (Spalding, 1995a), the
    modest increase from two to four fluids brings significant increases
    in realism.

    There appears to be no such coarse-grid possibility in the Monte-
    Carlo method, which requires a large number of particles in order
    to produce meaningful results at all.

    (2) Serious attention to computations based on larger numbers of
    fluids (eg 100) has only just started (Spalding, 1995b).

    Therefore no significant use of the many available economising
    techniques has been made. Doing so is only a matter of time,
    funding, and a sufficient drive from the practical users.


    --------------------------------------------------------------------

    3. Recent work on MFM
    ---------------------

    3.1 Reports and publications
    ----------------------------
    (a) The first MFM report
    ------------------------

    The paper presented at the 2nd Helsinki Colloquium on Process
    (Spalding,1995a) was the start of a series of publications on MFM,
    each of which has carried the idea a little further. It showed how
    a four-fluid model of premixed-fuel combustion allowed phenomena to
    be predicted, and specifically those concerning the spread of a
    steady flame confined in a duct (Williams, Hottel and Scurlock,
    1953) which were beyond the capability of the "two-fluid" eddy-
    break-up model. Much of this demonstration is reproduced in section
    4.4 below.

    A ten-fluid model of an axi-symmetrical turbulent jet was also
    presented.

    This paper already introduced the notions of "promiscuous coupling"
    and "Mendelian splitting", which have proved to be central
    ingredients of later multi-fluid models; and it was there pointed out,
    but not illustrated, that two-, three- and more-dimensional
    population grids could be envisaged and used.

    The essential idea was recognised as that of "discretization" of
    fluid-attribute space, which can be carried out in numerous ways,
    all of which should be equivalent if the "population grid" is made
    sufficiently fine.

    The "fluid-population distributions" of multi-fluid models and the
    "probability-density functions" of probabilistic methods (Dopazo and
    O'Brien, 1974) were regarded as being equivalent, again in the
    fine-grid limit.


    (b) The first two-dimensional population
    ----------------------------------------

    At the CTAC-95, Conference in Melbourne, a paper was presented
    (Spalding,1995b) which applied MFM to the turbulent Bunsen burner.
    This required the use of a two-dimensional fluid-population
    distribution, indeed the one which appears in section 1.6, above.

    The results were plausible, but no attempt was made to ensure grid-
    independence, or to make comparisons with experiment.


    (c) Detailed fpd studies for the "well-stirred reactor"
    -------------------------------------------------------

    A 100-fluid model was used (Spalding, 1995c) for the study of the
    shapes of the fluid-population distributions, both one- and two-
    dimensional, which appear (according to MFM) when combustible gases
    flow steadily into a reactor which is so "well-stirred" that
    macro-mixing is perfect.

    It was shown there that the one-dimensional fdps could take a great
    variety of different shapes, far exceeding those which have been
    imagined by the practitioners of the "presumed-pdf" approach
    [see section 4.1f below].

    (d) Steady one-dimensional flame propagation
    --------------------------------------------

    The same conclusion about pdf shapes emerged from a later study
    (Spalding, 1996a), in which MFM modelling was applied to a one-
    dimensional steadily-propagating turbulent pre-mixed flame.

    This study further demonstrated that the concept of grid refinement,
    which is familiar to all flow-simulation experts in connexion with
    the sub-division of time and of geometrical space, is equally
    applicable to the population grid. It was further shown that, for
    the case in question, results of adequate accuracy could be obtained
    with 10- or 20-fluid models; litlle further increase of accuracy
    resulted from using the 100-fluid model.

    (e) Transient flame propagation
    -------------------------------

    The four-fluid model has been employed in two transient explosion
    computations. One of these concerns a laboratory experiment in which
    flame propagates through a pre-mixed gas in a duct on the walls of
    which are a series of turbulence-promoting baffles. The results,
    reported by Freeman and Spalding (1995) are in good qualitative and
    semi-quantitative agreement with experiment.

    The second application has been to the explosions recently carried
    out at Spadeadam, in Northern England, by the UK Steel Construction
    Institute. The results have not yet been published. It can however
    be said that the computations based upon the 4-fluid model were far
    from providing the worst agreement with experimental data.

    Later research has shown that it would be wiser to employ a 10-fluid
    model, in order that the population-grid effects should not be
    troublesome.

    However, the ways in which the presence of solid obstacles influence
    the intensity and length-scale of turbulence, and so of coupling-
    splitting rates, is at the present a matter of guess-work.
    No-one living in the world today, in the present author's belief,
    can quantify this influence reliably.


    3.2 Unpublished work
    --------------------

    More recently (Spalding,1996d) the following studies have been made,
    albeit not yet completed:-

    (1) grid-refinement studies for the well-stirred reactor, operating
        in the mode in which a single cold-gas stream enters and its flow
        rate is increased until sudden extinction occurs;

    (2) simulations of transient macro- and micro-mixing in a paddle-
        stirred reactor containing initially-separated acidic and
        alkaline liquids.

    (3) simulations of combustion of fuel and air in a gas-turbine-type
        chamber; and

    (4) a study of the plane turbulent non-reacting mixing layer.

    The first study has confirmed the validity, in different
    circumstances, of the above-mentioned findings about grid
    refinement.

    The second, which probably breaks new ground for chemical engineers
    concerned with reactor-vessel design, reveals than it is much more
    important to simulate the micro-mixing than it is the macro-mixing.

    The third, by introducing an idealised smoke-production reaction,
    shows how different are the amounts of this production predicted by
    single- and multi-fluid models.

    The fourth study may be regarded as being of especial interest
    because, whereas all the others relied on a hydrodynamic model of
    turbulence, either mixing-length or k-epsilon, to provide the
    multiplier of the coupling-and-splitting-rate term, this one did
    not.

    Instead, the rate was deduced from the fluid-population distribution
    of longitudinal velocity, and from a prescribed length scale. The
    first step towards a stand-alone MFM for turbulent flows has thus
    been taken. The next, namely the removal of the need for a length-
    scale prediction, is not far away.

    --------------------------------------------------------------------

    4.   Relations of MFM to other turbulence models
    ------------------------------------------------
    4.1 The four main classes of turbulence model
    ---------------------------------------------
    (a) Introduction
    ----------------
    Turbulence models, in the sense used here, are collections of
    concepts about the nature of turbulent fluids which may be expressed
    in mathematical form in such a way as to constitute a method of
    prediction.

    Questions to be answered by such a method include:

      for a given set of flow-defining conditions (inlet geometry,
      outlet location and shape, internal baffles and sources of heat
      and momentum, fluid density, viscosity and heat capacity, number
      of thermodynamic phases, etc),

      - will the flow be turbulent at all?

      - if so what will be the intensity of the turbulence?

      - what will be the consequent exchange rates of mass, momentum and
        energy between the fluid and the solids with which it is in
        contact?

      - what will be the rates of chemical reaction within the fluid?

    The answers to all these questions are of course different for each
    location in space and instant in time.

    Methods of answering such questions are numerous and varied in
    nature; but it appears useful to separate them into four main
    classes, namely as:

      I   single-fluid models,
      II  multi-fluid models,
      III probabilistic models,
      IV  direct-numerical-simulation models.

    These terms will now be explained.

    (b) Single-fluid models; Class I
    --------------------------------

    Single-fluid models are those in which the fluid condition
    is characterised by average values of velocity, temperature, etc, at
    each location and time instant, and by statistical representations
    of the fluctuations about those averages.

    This class contains most of the models which are described in the
    text-books and are in use today. They include, to use the
    nomenclature of Launder and Spalding (1972):

    - zero-(differential)-equation models, such as those of Boussinesq
      (1877) and Prandtl (1925);

    - one-equation models, such as those of Prandtl (1945) and Bradshaw
      (1967);

    - two-equation models, such as those of Kolmogorov (1942), Harlow
      and Nakayama (1968), Saffmann (1970), Spalding (1969), Wilcox
      (1993), Yakhot and Orszag (1986) and many others;

    - multi-equation models, such of those concerned with multiple
      levels of energy or scale (Elhadidy, 1980; Kim and Chen, 1989) or
      with further statistical quantities (Daly and Harlow, 1970;
      Spalding, 1971b; Naot, Shavit and Wolfshtein, 1974).

    The methods of the Reynolds-stress kind fall into the last-named
    grouping.


    (c) Multi-fluid models; Class II
    --------------------------------

    Multi-fluid models express the fluctuations by representing them as
    though there were many different fluids mingling within the same
    space, each with its own local and instantaneous velocities and
    temperatures.

    The mingling of steam bubbles with water in a kettle is an extreme
    example; but the fluids are more usually of a single thermodynamic
    phase, as when tongues of flame rise above a garden bonfire.

    The basic concept was already present in the writings of Reynolds
    (1874) and Prandtl (1925); for they both conceived of the relative
    motion of fluid fragments, which could be of significant size.

    Scientists concerned with combustion (Shchelkin,1943; Wohlenberg,
    1953; Howe and Shipman, 1965; Kuznetzov,1979) have found the concept
    much to their liking; for it is only by taking into account the
    fragmentariness of the burning gases that observed phnomena can be
    explained, even qualitatively.

    The author would place in the two-fluid category:

    - his own "eddy-break-up" model (Spalding, 1971a; Mason and
      Spalding,1973),

    - the "eddy-dissipation concept" of Magnussen and Hjertager (1976),
      which was derived from it, and

    - the Bray-Libby-Moss "flamelet" model (Bray and Libby,1981; Moss,
      1980).

    A more extensive exploration of the capabilities of two-fluid models
    was made later by the author and his co-workers (Spalding, 1982,
    1983, 1984, 1987, 1993; Andrews, 1986; Fueyo, 1992; Malin, 1986; Xi,
    1986; Wu, 1987; Ilegbusi and Spalding, 1987a,b; Pericleous and
    Markatos, 1991).

    In these works, the two fluids were distinguished in a variety of
    ways, eg turbulent/non-turbulent, hotter/colder, upward/downward-
    moving; and two sets of Navier-Stokes equations were solved
    simultaneously.

    Still more recently, as is discussed elsewhere in this paper,
    the multi-fluid idea has been extended to 4-, 14-, and many-fluid
    models (Spalding, 1995a,b,c). However the first publications on the
    population-of-fluids concept may have been those arising from study
    of the long-forgotten "ESCIMO" model of turbulent combustion
    (Spalding, 1979; Noseir, 1980; Tam, 1981; Sun, 1982).


    (d) Probabilistic models; Class III
    -----------------------------------

    Probabilistic models represent the near-randomness of turbulence by
    introducing some randomness of their own. Specifically, they employ
    Monte-Carlo methods to establish the probable distribution of fluid
    attributes within a multi-dimensional space of which the coordinates
    include the components of velocity, the temperature, etc.

    This class appears to have originated in the chemical-engineering-
    science field with the publications by Curl (1963) and Dopazo and
    O'Brien (1974).

    More recently Pope (1980, 1982, 1985), Chen and Kollmann (1988,
    1990) and others have conducted a vigorous research campaign, which,
    to judge from the proceedings of a recent conference (ICOMP, 1994),
    is gathering momentum.

    Once again, it is the desire to be able to simulate combustion
    phenomena with quantitatuve realism which provides much of the
    motivation.


    (e) DNS models; Class IV
    ------------------------

    Finally, direct-numerical-simulation models revert to the single-
    fluid approach, but without any built-in averaging. Their users
    solve the Navier-Stokes equations with extremely fine sub-divisions
    of space and time, deriving averages and statistical measures of
    fluctuations only after the solutions have been obtained.

    Direct-numerical simulation, has been subjected to increasing
    attention as computers have become more and more powerful. Pioneers
    include Schumann (1973) and Reynolds (1975).


    (f) "Presumed-pdf" models
    -------------------------

    There is a class of model featuring "probability-density functions"
    which the present author would place in the single-fluid category,
    because the shapes of the pdfs are presumed rather than computed.
    Models of this kind may perhaps have started with the author's own
    paper (Spalding,1971b), in which a "two-spike" double-delta-function
    presumption was made.

    Subsequently, Lockwood and Naguib (1975) made the "clipped-Gaussian"
    presumption; they were followed by others, including Kent and Bilger
    (1976), Kolbe and Kollmann (1980), Rhodes et al (1974) and Gonzalez
    and Borghi (1991).

    Employment of these presumptions gave rise to additional
    computational expense, because all reaction-rate terms (for example)
    required the evaluation of integrals for each cell at each
    iteration; but no clear advantage in respect of generality of
    agreement with experiment (in the present author's opinion) ever
    emerged.

    As is shown elsewhere in this paper, the multi-fluid model makes it
    possible to calculate the shape of the pdf (or FPD, to accord with
    MFM nomenclature) if the number of fluids is very large; and this
    shape depends greatly upon the dimensionless parameters which relate
    the coupling rates of the fluids on the one hand, and their chemical
    reaction rates on the other, to the mass flow rate into the cell.

    A detailed parametric study of these effects by way of the multi-
    fluid model might perhaps lead to formulae with the aid of which one
    might know what shape to presume in particular circumstances.
    Presumption of the shape without such guidance is however likely to
    do as much harm as good.


    4.2 Advantages and disadvantages of the four classes
    ----------------------------------------------------
    (a) Ease of use
    ---------------

    CLASS I must be regarded as the easiest to use, because many members
    of the class are embodied in widely-available software packages. For
    example, the Shareware (and therefore freely available) version of
    PHOENICS contains zero-, one-, and two-equation models; and, because
    of its open-source character, it allows other models of Class I to
    be built into it.

    The latest versions of PHOENICS have more than a dozen identifiable
    Class-I turbulence models, and indeed many more if the possibility
    of sub-model interchanges betwen them are considered.

    Most other general-purpose computer codes (eg FLUENT, FIDAP,
    FLOW-3D, CFX, TASC) now have a similarly wide range of models.

    Ease of use is greatest for the 0-, 1-, and 2-equation models, and
    least for the Reynolds-stress models, for which convergence is not
    always guaranteed.

    CLASS II models are less widely available in packages, although
    Shareware PHOENICS does possess a 2-fluid model which solves two
    sets of Navier-Stokes equations.

    Multi-fluid models of this class are not yet widely available. Those
    which will be described below have been implemented by use of the
    open-source facilities of the latest PHOENICS; and it is probable
    that the implementation will be attached to the next or next-but-one
    release of the code.

    There are no publicly available codes which embody the probabilistic
    models of CLASS III, but the US Government makes some available to
    its contractors. Although the methods may be simple for those who
    have become used to them, or to Monte-Carlo methods in other
    circumstances, getting started with Class-III methods is difficult.

    Recently, however, a paper has been published which shows how the
    general-purpose code PHOENICS can be used as the basis for Monte-
    Carlo calculations of turbulent-mixing and -combustion processes
    (Fueyo, Larroya, Valino and Dopazo,1995)

    Class IV methods would probably be regarded as the easiest to use,
    were it not for the fact that not even the largest computers in
    the world are large enough. On the surface, they appear to require
    nothing more than the ability to make a time-dependent non-
    turbulent flow simulation. But large-grid-size problems bring
    their own special difficulties, and means for resolving them; so
    this is also not an easy field of research to get into.


    (b) Ease of understanding
    -------------------------

    For the reasons just alluded to, Class-IV methods can be regarded
    as the easiest to understand, because the physical assumptions are
    those of laminar flow; and indeed there is no "modelling" (in the
    sense of substituting guesswork for ignorance) at all.

    The present author would rank the other classes in the order
    I II III, in respect of understandability, with the qualification
    that some of the more complex single-fluid multi-equation models
    may be more difficult to grasp than the simpler two-fluid models.

    It has to be admitted, however, that the two-fluid concept has not
    so far achieved much popularity; and this may in part result from
    reluctance to make the imaginative leap from steam/water mixtures
    on the one hand to hot-air/cold-air mixtures on the other.

    Another misgiving has perhaps deterred more reflective persons,
    namely the thought: Why only TWO fluids? Perhaps the demonstration
    in the present paper that one can indeed handle any number of fluids
    will bring some reassurance.

    The probabilistic models of Class III, are not, to this author's
    mind, easy to grasp at all. It is not so much the idea of multi-
    dimensional space that is difficult, but rather the esoteric
    mathematical language and symbols which, perhaps necessarily, the
    practitioners of these models employ.

    To those who make the effort, Classes III and II can at first
    appear to be almost the same; so the next difficulty is that of
    understanding the differences. Of these, a crucial one is:-

    * Class-II methods discretize some (but not all) fluid-attribute
      dimensions, in the manner familiar to all users of finite-volume
      or finite-element codes; and they can obtain useful results from
      very coarse discretizations, as use of the two-fluid model has
      shown.

    * Class-III methods work in a non-discretized, multi-dimensional,
      all-fluid-attribute space; "coarsening" for economy appears to
      require reduction of the number of the attributes (eg species
      concentratons) which are considered.

    It might be said that the Monte-Carlo approach of probabilistic
    methods is akin to the use of a shot-gun: many pellets are needed if
    the target is to be hit; whereas Class II methods work with rifles
    which, if carefully aimed, need discharge fewer bullets.


    (c) Extent of validation
    ------------------------

    Class I is by far the most extensively validated of the four; and
    although no widely-agreed answer can be given to the newcomer's
    reasonable question, "Which model is best?", advice can be given,
    based upon experimental knowledge, as to when it is permissible to
    use Prandtl's rather-simple mixing-length model, and when only a
    Reynolds-stress model will suffice.

    However, it is just because of this extensive validation campaign
    that it is possible to assert that NONE of the Class-I models will
    suffice for the explanation of commonly-observed and practically-
    important phenomena. This has been suspected for many years; and
    the passage of time justifies the substitution of conviction for
    suspicion.

    None of the other three models have been subjected to quantitative
    tests on the same scale. Probably the work on the two-fluid model
    by the author's former students at Imperial College represents the
    most systematic campaign; but this came to an end in 1988.

    Comparison of Class-III predictions with experiments is being made
    in the USA for both simple and complex (eg gas-turbine-combustor)
    circumstances. The former comparisons appear to have drawn
    attention to a deficiency of all of the "mixing-model" assumptions
    which have been tried so far: they do not produce the right
    (Gaussian) probability-density-function shape in the late-time limit.

    It is not clear whether this difficulty can be resolved while
    retaining the Monte-Carlo framework of calculation.



    4.3 Similarities and differences between MFM and previous models
    ----------------------------------------------------------------

    (a) Introduction
    ----------------

    The brief review, in sections 4.1 and 4.2, of the whole turbulence-
    modelling scene, has been provided in order that the multi-fluid
    turbulence model can be understood in context.

    MFM is not entirely new; but it has some significant points of
    novelty and advantage. The purpose of the present section is to
    bring these to the reader's attention, by pointing out similarities
    and differences.


    (b) Comparison with single-fluid models
    ---------------------------------------

    The basic mathematical structure of MFM is similar to that of
    single-fluid models of turbulence, in that differential equations
    are formulated and solved in order to represent the influences of:

    - time-dependence,
    - convection,
    - diffusion and
    - sources and sinks.

    MFM may also make use of single-fluid concepts such as "effective
    viscosity", "epsilon/k", mixing length, and root-mean-square
    fluctuations; however it also possesses other means of describing,
    calculating and responding to the influences of the fluctating
    nature of turbulent flows.

    Like single-fluid turbulence models, MFM has its empirical
    "constants", which may turn out to be functions of the local
    Reynolds number of turbulence and of other dimensionless quantities;
    but, so far, little effort has been expended on determining what
    they are.

    MFM has many more equations to solve than single-fluid models do; so
    computer times are bound to be greater, if the same computational
    grid is used. However, MFM uses the available computer power in a
    different way; for it devotes as much attention to differences in
    fluid properties prevailing at a single location as it does to
    differences from place to place. Therefore MFM may use coarser
    geometrical grids than have become usual for single-fluid models.


    (c) Comparison with the 1982 two-fluid model
    --------------------------------------------

    MFM differs from the two-fluid model worked on by the author and
    his colleagues in the early 1980s in two main ways, neither of
    which are essential.
                                  ________________________________
    The first is that the         |                              |
    states of the two fluids  ^   |                  Fluid 2     |
    could roam freely over    |   |                    *         |
    the "phase space"         |   |                              |
    defined by the fluid     vert-|                              |
    attributes, whereas      ical |                              |
    MFM restricts the        attr-|                              |
    states of its fluids     ibute|         *                    |
    to discrete locations.        |     Fluid 1                  |
                                  |                              |
    As soon as adaptive           |______________________________|
    gridding is introduced           horizontal attribute ----->
    however, MFM can enjoy
    the same freedom.

    The second difference is that two complete sets of Navier-Skokes
    equations were solved for the two-fluid model, an achievement which
    may never be attainable for MFM, and has certainly not been
    attempted.

    Nevertheless, the discrete fluids of MFM can indeed move relative to
    one another under the influences of body forces and fluid-to-fluid
    friction; and it appears that to allow each fluid its own Navier-
    Stokes equations would be needless indulgence. A "drift-flux" or
    "algebraic-slip" approximaton is likely to suffice.

    It is worth noting that MFM has adopted one valuable feature of two-
    fluid thinking, namely the distinction between the fluid attributes
    which DISTINGUISH one fluid from the other and those which are
    merely POSSESSED by the fluid.

    Thus, in the two-fluid model, the vertical velocity may be that
    which distinguishes the fluids; but each fluid still has its own
    distinct values of all three velocity components.


    The difference between the models that IS essential, of course, is
    that between two and many; inevitably, MFM is the more capable of
    describing real flows.


    (d) Comparison with the presumed-pdf models
    -------------------------------------------

    MFM focusses attention on the probability-density functions of the
    fluid attributes for some of the same reasons as do users of the
    presumed-pdf approach, who nowadays form the majority of those who
    attempt to simulate turbulent combustion processes; for it is only
    knowledge of those functions which enables chemical-kinetic
    information to be rationally introduced.

    However, MFM calculates the pdfs (ie its fluid-population
    distributions) and shows that the functions may take very varied
    shapes; whereas the basis of selection used by the shape-presumers
    is far from reliable.

    Moreover, MFM finds it as easy to calculate two-dimensional pdfs as
    one-dimensional ones; presumed pdfs, by contrast, are perhaps always
    one-dimensional.

    Of course, there can be no certainty, in advance of systematic
    comparison with experimental data, that the pdfs calculated by MFM
    in its present form are correct; and there is good reason to believe
    that the promiscuous-Mendelian coupling-splitting hypothesis is not
    always the most realistic (see section 6.3).

    However MFM is at the beginning of its research life; whereas the
    presumed-pdf approaches may be near the end of theirs. It appears
    certain that MFM has the greater potential for development and
    improvement.


    (e) Comparison with the probabilistic models
    --------------------------------------------

    Here the abbreviation CPDFM (for Computed PDF model) is used in
    place of the longer "probabilistic model".

    There is certainly a more than superficial relation between the PDF
    of the CPDFMs and the FPD of the MFM. In the limit of an infinite
    number of fluids (for the latter) and of particles (for the former),
    they may become quantitatively identical.

    However, they ARE somewhat different concepts; they are computed by
    different mathematical techniques; they embody different physical
    hypotheses; and (possibly) they appeal differently to persons of
    differing educational backgrounds.

    In respect of mathematics, CPDFMs employ Monte-Carlo methods,
    whereas MFMs use direct numerical solution. This means that CPDFMs
    focus on particles, whereas MFMs focus on cells. A consequence may
    prove to be that, for problems for which they can be reasonably
    compared, CPDFMs will prove to use more computer time, and MFMs to
    use more storage.

    In respect of physics, CPDFMs employ something akin to "coupling",
    but not, it appears to "splitting".  Nothing like the Mendelian
    concept appears so far to have been introduced; and perhaps indeed
    to do so would be impossible within the Monte-Carlo framework.

    The 'mixing models' of the PDFMs are recognised, even by their
    strongest proponents, as being their major source of weakness; so it
    may be the constraints imposed by the mathematical technique which
    cause "relaxation to the local mean" to be used as the mixing model
    when a two-dimensional bluff-body flow is simulated (Correa and
    Pope, 1992), despite its known inadequacies.

    There are also practical differences: MFMs can provide meaningful,
    and perhaps adequate, results with rather coarse population grids,
    eg no more than 20 fluids; but it is not clear that a corresponding
    saving is possible with CPDFMs.

    Moreover, MFMs have always possessd, and need to retain, the notion
    that it is only few of the fluid properties which need to be
    discretised; whereas CPDFMs appear to adopt an all-or-nothing
    principle.

    Finally, to bring to an end a comparison which deserves much more
    study and space, papers describing CPDFMs are written in language
    of a highly mathematical character, whereas descriptions of MFMs use
    more words and fewer symbols. Some users of turbulence models may
    base their choice between MFM and CPDFM on which is described in the
    language which is more congenial to them.

    In summary, there appear to be sufficient differences between CPDFMs
    and MFMs, and (in the author's view) sufficient advantages on the
    side of the MFMs, to render the latter well worthy of the attention
    of those researchers who aim to provide engineers with the best
    predictive tools for turbulent combustion.


    (f) Comparison with the eddy-break-up model
    -------------------------------------------

    The eddy-break-up model (Spalding, 1971a) has had a surprisingly
    (in view of its crudeness) long-lasting influence: it inspired the
    invention of  the eddy-dissipation concept (Magnussen,1976); it
    influenced the development of Pope's probabilistic model (1982); and
    it can usually be found to underly many models of the presumed-pdf
    kind.

    It is now possible also to discern that the multi-fluid model which
    is described in the present paper is merely the long-delayed
    extension of the ideas of the original paper.

    All persons with experience in computational fluid dynamics know
    that, in order to obtain improved accuracy, it is necessary to
    "refine the grid" sufficiently. MFM, it can be said, is merely a
    grid-refined eddy-break-up model.

    In order to make clear the steps which led from EBU to MFM, the next
    section summarises part of the relevant publication (Spalding,
    1995a).


    4.4 Generalising the eddy-break-up model
    ----------------------------------------

    (a) The two-fluid nature of the eddy-break-up model
    ---------------------------------------------------

    The central idea of the EBU was that, where the time-average
    temperature of the gas was above that of the fully-unburned mixture,
    TU, and below that of the fully-burned mixture, TB, this was the
    consequence of its being made up of colder fragments inter-mingled
    with hotter fragments.

    Moreover the cold fragments were as cold as they could be, namely of
    temperature TU; and the hot fragments were as hot as they could be,
    namely at temperature TB.

    The mass fractions of cold and hot gas, MU and MB, were given by:

          MU = 1 - MB

             = (TB - TU) / (TB - TU).

    Of course, combustion could not take place in the cold fragments;
    because they were too cold: nor in the hot fragments, because they
    were too hot.

    Since combustion undoubtedly does take place, it was supposed that
    it did so at the inter-faces between the two types of fragments; but
    these were supposed to occupy only a small proportion of the volume.

    (b) The rate-determining process
    --------------------------------

    The suppositions were then made that:-

    (1) the rate of combustion per unit volume
            was proportional to

        the rate of intermingling of the two types of fragments; and

    (2) this rate was proportional to:

        MB * MU * MIXRATE
            and

    (3) the magnitude of the variable MIXRATE
            was proportional to either:

        VELOCITY_GRADIENT * MIXING_LENGTH
           in the first model (Spalding, 1971b)

           or to:

        TURBULENCE_ENERGY_DISSIPATION_RATE / TURBULENCE_ENERGY
            in a later version (Mason and Spalding,1973).

    (c) Successes and failures of EBU
    ---------------------------------

    The model was successful in explaining certain otherwise
    inexplicable experimental findings, for example the fact that the
    angle subtended by the flame anchored in a plane-walled channel was
    nearly independent of approach-gas velocity.

    However, it had no means for expressing the influence of chemical
    kinetics. Yet such an influence does exist, as witness the fact
    that, when the approach-gas velocity becomes very large, flame
    propagation abruptly ceases.

    This latter shortcoming was addressed by the author already in the
    original reference, and later by Magnusson and Hjertager (1976), in
    rather similar ways; but, at least in the present author's opinions,
    nothing truly satisfactory emerged.

    (d) The four-fluid model
    ------------------------

    Nearly twenty-five years elapsed before (what now seems) an obvious
    next step was taken, namely to increase the realism of the model by
    increasing the number of fluids, ie, in the terms of the present
    paper, to refine the population grid.

    The central idea was to suppose that the turbulent reacting mixture
    consisted of fragments of more than two kinds. Four being the
    minimum number which, it then appeared, would explain all the
    qualitatively observed facts, it was the number first used.

    Specifically, the mixture was supposed to consist of four fragment
    classes, namely:

    A,  comprising fully-unburned gas;
    B,  comprising a mixture of unburned gas and combustion products of
        too low a temperature for the chemical reaction rate to proceed
        with significant speed;
    C,  comprising a mixture of unburned gas and combustion products of
        higher temperature, at which the chemical reaction rate is
        significant; and
    D,  comprising fully-burned gas.

    The eddy-break-up model had allowed for only A and D; the new
    4-fluid model allowed for B and C in addition.


    (e) The rates of production of A, B, C and D
    --------------------------------------------

    By analogy with the formulation successfully used for the EBU, it is
    postulated for the 4-fluid model that, per unit mass of mixture:

    *  fluid B is produced from fluids A and D at the rate:
            0.5 * MA * MD * MIXRATE

    *  fluid B is produced from fluids A and C at the rate:
            MA * MC * MIXRATE

    *  fluid C is produced from fluids A and D at the rate:
            0.5 * MA * MD * MIXRATE

    *  fluid C is produced from fluids B and D at the rate:
            MB * MD * MIXRATE

    In addition to these rates of production of one fluid (and
    corresponding diminution of the others):

    *   fluid D is produced from fluid C (only) at the rate:
            MC * CHEMRATE

    CHEMRATE is, of course, the rate of chemical conversion of the
    unburned fuel in fluid C to combustion products. It is the element
    of the model which allows (as EBU did not) chemical-kinetic
    effects to be introduced rationally.


    (f) The first application of the 4-fluid model
    ----------------------------------------------

    The first application of the 4-fluid model was to the same problem
    as that for which the EBU was invented, namely that of steady
    turbulent flame spread in a plane-walled duct, through the upstream
    end of which flows both a stream of unburned combustible mixture,
    and a separate but thinner stream of fully-burned products to serve
    as an igniter, shown below.
------------------------------------------------------------------------

    Steady turbulent flame spread in a plane-walled duct

    Entrance ##################### Channel wall ############## Exit
  unburned gas->                         ::::::::::::          ->
      "     --->               ::::::::::                      -->
      "     --->      :::::::::  flame spreading to wall       --->
      "     --->::::::                                         ---->
  burned gas -->::  - --- - --- - --Symmetry plane - --- - --- ----->


------------------------------------------------------------------------

    The first remarkable feature of the experimental results is that
    the rate of spread of the flame, as measured by its angle for a
    fixed inlet-stream velocity, is very little dependent on the fuel-
    air ratio of the incoming mixture; until, that is to say, this
    ratio becomes too rich or weak to sustain combustion at all.

    The second remarkable feature is that the angle of the flame depends
    very little on the inlet-stream velocity. This implies that the
    rate of combustion somehow keeps pace with an increased supply of
    reactants.


    (g) The numerical simulation
    ------------------------------

    This process was simulated by activating the 4-fluid model within
    the PHOENICS computer program operating in "parabolic" (ie marching-
    integration) mode.

    The magnitude of the MIXRATE quantity was computed in the same
    manner as for the eddy-break-up model, ie as proportional to the
    square root of the turbulence energy divided by the length scale;
    and

    the quantity CHEMRATE was taken as equal respectively to

        10000, 1000, 100 and 10 times

    that which would suffice, if it were Fluid C which filled the duct,
    to consume all the fuel by the time that the duct exit was reached.

    The following figures show, for the four different CHEMRATE values,
    the profiles of concentration of fluids A, B, C and D across the
    duct, at a distance half-way between the entrance and the exit.


------------------------------------------------------------------------
    Profiles of A, B, C, and D for CHEMRATE = 10000
              The symmetry plane is on the left, the wall on the right

       VARIABLE        D          C          B          A
         MINVAL=  0.000E+00  0.000E+00  0.000E+00  0.000E+00
         MAXVAL=  1.000E+00  1.000E+00  1.000E+00  1.000E+00
         CELLAV=  3.992E-01  7.687E-04  1.871E-01  4.126E-01

     1.00 +....+....+....+....+....+....+....+....+A.A.A.AA.A
          .                                       A         .
     0.90 D DD                                              +
          .    D D                              A           .
     0.80 +        DD                                       +
          .           D                       A             .
     0.70 +             D                   A               +
          .               D                                 .
     0.60 +                D               A                +
          .                  D                              .
     0.50 +                    D         A                  +
          .                     D                           .
     0.40 +                       D B  A                    +
          .                    BB B D AB B                  .
     0.30 +                B B      A D    B                +
          .             B B       A         B               .
     0.20 +         B B         A      D      B             +
          .      B B           A         D                  .
     0.10 B BB B           A A             DD   B           +
          .         A A A A                   D D B         .
     0.00 A.AA.A.A.AC.C.C+CC.C+CC.C+C.CC+C.CC+C.C.CB.B.B.BB.B
          0   .1   .2   .3   .4   .5   .6   .7   .8   .9  1.0
     the abscissa is       Y  .  min= 8.33E-04 max= 4.92E-02

------------------------------------------------------------------------

    Profiles of A, B, C, and D for CHEMRATE = 1000
              The symmetry plane is on the left, the wall on the right
       VARIABLE        D          C          B          A
         MINVAL=  0.000E+00  0.000E+00  0.000E+00  0.000E+00
         MAXVAL=  1.000E+00  1.000E+00  1.000E+00  1.000E+00
         CELLAV=  3.722E-01  7.029E-03  1.883E-01  4.322E-01
     1.00 +....+....+....+....+....+....+....+....+A.A.A.AA.A
          .                                       A         .
     0.90 D D                                   A           +
          .  D D                                            .
     0.80 +      D D                          A             +
          .         D                       A               .
     0.70 +           D                                     +
          .             D                  A                .
     0.60 +               D                                 +
          .                D             A                  .
     0.50 +                  D                              +
          .                    D       A                    .
     0.40 +                     D B B A                     +
          .                  B BB D A BB                    .
     0.30 +               BB        D    B                  +
          .           B B         A   D    B                .
     0.20 +         B           A      D    B               +
          .    B B B         A A         D    B             .
     0.10 B BB            AA               D    B           +
          .        AA A A                   D D   B         .
     0.00 A.AA.A.A.CC.C.C+CC.C+CC.C+C.CC+C.CC+C.C.CB.B.B.BB.B
          0   .1   .2   .3   .4   .5   .6   .7   .8   .9  1.0
     the abscissa is       Y  .  min= 8.33E-04 max= 4.92E-02

------------------------------------------------------------------------

    Profiles of A, B, C, and D for CHEMRATE = 100
              The symmetry plane is on the left, the wall on the right
       VARIABLE        D          C          B          A
         MINVAL=  0.000E+00  0.000E+00  0.000E+00  0.000E+00
         MAXVAL=  1.000E+00  1.000E+00  1.000E+00  1.000E+00
         CELLAV=  2.234E-01  3.924E-02  1.919E-01  5.451E-01
     1.00 +....+....+....+....+....+....+....+A.A.AA.A.A.AA.A
          .                                 A               .
     0.90 +                                A                +
          .                              A                  .
     0.80 +                                                 +
          D D                          A                    .
     0.70 +  D                        A                     +
          .    D                                            .
     0.60 +      D                  A                       +
          .        D                                        .
     0.50 +         D             A                         +
          .           D         A                           .
     0.40 +             B BB B BB                           +
          .           B        A  B                         .
     0.30 +        BB     D  A      B                       +
          .      B         A          B                     .
     0.20 +  B B          A  D         B                    +
          B B           A      DD        B                  .
     0.10 +  C C C CA A C C       D        B                +
          C CA A A A       C C CC C C DD    B               .
     0.00 A.A..+....+....+....+....+..CC+C.CC+B.B.BB.B.B.BB.B
          0   .1   .2   .3   .4   .5   .6   .7   .8   .9  1.0
     the abscissa is       Y  .  min= 8.33E-04 max= 4.92E-02

------------------------------------------------------------------------

    Profiles of A, B, C, and D for CHEMRATE = 10
              The symmetry plane is on the left, the wall on the right

       VARIABLE        D          C          B          A
         MINVAL=  0.000E+00  0.000E+00  0.000E+00  0.000E+00
         MAXVAL=  1.000E+00  1.000E+00  1.000E+00  1.000E+00
         CELLAV=  8.441E-02  4.799E-02  8.105E-02  7.864E-01
     1.00 +....+....+....+...A+AA.A+A.AA+A.AA+A.A.AA.A.A.AA.A
          .                A                                .
     0.90 +               A                                 +
          .                                                 .
     0.80 +                                                 +
          .             A                                   .
     0.70 D                                                 +
          .                                                 .
     0.60 + D         A                                     +
          .                                                 .
     0.50 +         A                                       +
          .  D                                              .
     0.40 +                                                 +
          .    D B AB                                       .
     0.30 +    B                                            +
          .  B C A    B                                     .
     0.20 C C    C      B                                   +
          . B  A   C                                        .
     0.10 B  A      C     B                                 +
          . A         C C  B                                .
     0.00 A....+....+...D+CC.B+BB.B+B.BB+B.BB+B.B.BB.B.B.BB.B
          0   .1   .2   .3   .4   .5   .6   .7   .8   .9  1.0
     the abscissa is       Y  .  min= 8.33E-04 max= 4.92E-02

------------------------------------------------------------------------

    (h) Discussion of the results
    -----------------------------

    Despite the crudeness of the line-printer plots, the following
    observations can be made:

    (1) In the first three figures, the differences between the
        profiles of the unburned-gas (A) and the burned gas (D)
        concentrations are rather small. Evidently the change of
        CHEMRATE from 10000 down to 100 has had little effect of the
        angle of the flame, in qualitative agreement with experiment.

    (2) However, the average values of the concentration of the reactive
        gas (C) (see CELLAV at the tops of the figures) are very small
        only when CHEMRATE is high (7.87E-04 for 10000, 7.03E-03 for
        1000). For CHEMRATE = 100, the CELLAV of C has risen to 3.92E-02,
        signifying that this gas is not now so easily able to convert
        into products (D) all that is created by micro-mixing

    (3) When, for the fourth figure, CHEMRATE drops to 10, the profiles
        become radically different: the burned-gas (D) concentration
        becomes small everywhere, signifying that the incoming hot-gas
        stream has simply been mixed with the colder unburned-gas stream
        without promoting any significant further reaction; and the
        concentration of C rises as a result of this mixing, without any
        of the diminution which would result from such reaction.

    (i) Concluding remarks
    ----------------------

    The above example has been discussed at some length because it shows
    how "refinement of the composition grid" enables new insights to be
    gained and improved realism to be added to a numerical simulation.
    It also perhaps has some historic importance.

    Later studies have shown that, as any CFD expert would expect,
    grid refinement from 2 to 4 is not enough for grid-independent
    results to be attained; but (depending on the presumption made for
    the chemical kinetics) as few as 10 fluids will give acceptable
    accuracy.

    --------------------------------------------------------------------
    5. Application prospects
    -----------------------------
    5.1 The possibility of immediate use
    ------------------------------------

    (a) The implications of already-reported results
    ------------------------------------------------

    In the author's opinion, the multi-fluid model is ready for use,
    right now, for all simulations in which chemical reactions in
    turbulent fluids are of prime importance. This opinion holds both
    for chemical-engineering equipment such as the paddle-stirred
    reactor, and for combustion chambers of engines and furnaces.

    In so far as a turbulence model such as k-epsilon can be relied on
    for simulating the hydrodynamics, there is only one constant needed
    by the model, namely that which relates the "fluid-coupling" rate to
    epsilon/k; and this is known to be of the order of unity.

    The results for the rate of reaction in the stirred vessel in
    (Spalding, 1996b) show very clearly the danger of using a single-
    fluid model. In particular, the common belief  that single-fluid
    models are usable, when reaction is slow compared with turbulent
    mixing, is no longer tenable.

    It has also been shown that the computer-time penalty associated
    with use of MFM is not a serious one.

    Therefore, it can be plausibly argued, some form of MFM is the
    ONLY turbulence model which can be justified for use with
    homogeneous-reaction simulation from now on.

    (b) The available alternatives
    ------------------------------

    A counter-argument may be raised by users (who are numerous) of
    the "presumed-pdf" approach, to the effect that they are already
    taking sufficient account of fluctuations.

    To rebut this counter-argument, one has only to inspect the computed
    fluid-population distributions which have appeared in the recent
    publicationss cited in the present paper.

    These FDPs have the same significance as the PDFs; and they are of
    such varied shapes that to suppose that anyone could safely
    "presume" them would be extremely rash.

    How about the use of a probabilistic method instead? If the computer
    times are low enough, and the computer programs easy enough to use,
    the case for using such methods is as strong as for MFM. Certainly,
    although the underlying micro-mixing models are different, there is
    no quantitative evidence which yet favours one rather than the
    other.


    5.2 Applications to chemically-reacting flows
    ---------------------------------------------

    (a) A list of the main application areas
    ----------------------------------------

    Certainly the reacting-flow applications of MFM are of high
    importance. They include:-

    (1) the calculation of the efficiency of combustion of:
       gaseous, liquid and/or solid fuels in furnaces, whether for:
       - production of electricity via steam-raising;
       - the production of materials such as steel, glass, bricks,
         pottery or petroleum products;
       - the destruction of waste; or
       - the heating of buildings or cities;

    (2) the calculation of the efficiency of combustion of:
       gaseous and/or liquid fuels in engines, whether:
       - gas turbines for aircraft propulsion;
       - gas turbines for power production at ground level;
       - gasoline-fired spark-ignition engines; or
       - oil-fired diesel engines;

    (3) the calculation of the rates of production of (usually)
        undesirable side-products such as:
       - oxides of nitrogen;
       - smoke;

         in any of the above types of equipment;

    (4) the simulation of reactions other than combustion in reactors
        of the kinds employed in the chemical, pharmaceutical and
        food-manufacturing industries;

    (5) the simulation of flame propagation in hazardous circumstances
        such as those arising in off-shore oil-platform explosions,
        where the information sought is the pressure-versus-time profile
        on bounding walls.


    (b) The relevance of MFM to the computation of combustion efficiency
    --------------------------------------------------------------------

    Let item (1) be considered first; and let it be supposed that the
    temperature is so high that, at the prevailing pressure (usually
    atmospheric), all energy-producing reactions proceed very rapidly
    towards completion.

    If a single-fluid model is employed, the reaction will be supposed
    to be completed at the downstream section at which the time-average
    fuel/air mixture fraction is stoichiometric or less.

    A multi-fluid model will reveal that some unburned fuel is present
    at locations downstream of that section; and the equipment designer
    seeking 100% combustion efficiency will want to know how much larger
    he must make his combustion chamber.

    The answer will of course depend on the details of the geometry of
    the system; which is why a detailed three-dimensional analysis is
    needed.

    MFM is already able to provide this, with some margin of error, of
    course.

    MFM, in order to provide reliable predictions for most
    circumstances, will need to take account of the fact that light-gas
    fragments are accelerated (or decelerated) by pressure gradients
    more easily than dense-gas fragments, so giving rise to what is
    sometimes called "counter-gradient diffusion". The phenomenon is
    perhaps most important (and least taken into account) in simulations
    of reciprocating-engine (ie gasoline or diesel) combustion.

    In gas turbines, perhaps, a one-dimensional fluid population will
    suffice, the distinguishing attribute being the fuel-air ratio.

    Gasoline engines may also be simulated adequately by a one-
    dimensional population, with however reactedness (or mixture
    fraction minus unburned fuel) as the distinguishing attribute.

    Diesel engines, on the other hand, because the mixing and the
    ignition processes are not separated in time, will probably require
    two-dimensional populations of the kind discussed in section 1.6.


    (c) The computation of pollutant production
    -------------------------------------------

    There are innumerable papers which describe how, once the pdf is
    known, the rate of production of oxides of nitrogen, smoke, and
    other undesirable side-products, can be computed. To this author's
    surprise, many of these claim to have demonstrated sufficient
    agreement between predictions and experimental data to justify the
    assumptions.

    The methods of calculation are ingenious, and frequently very
    time-consuming; yet they almost invariably rest on precarious
    foundations. One of these is usually some variant of the two-fluid
    eddy-break-up model; and the other will be some presumption about
    the probability-density of fuel/air ratio or reactedness, but not
    usually of both.

    Sometimes the pdf is taken to be that which corresponds to an
    assemblage of flamelets, ie one-dimensional steadily-propagating
    pre-mixed flames.

    The author's current view is that the PDFs computed by the MFM as it
    stands today, ie without any refinements whatsoever, are likely to be
    as least as realistic as any that are currently in use.

    His recommendation is therefore to combine the best chemical-kinetic
    models of the current literature with the pdfs, one- or
    two-dimensional according to circumstance, computed by the MFM.


    (d) Flame propagation in explosion circumstances
    ------------------------------------------------

    As already mentioned, the four-fluid model has been used with some
    success; and it can be be expected that further insight will be
    gained when:

    (1) 10 or more fluids are employed, so that the chemical-kinetic
        data can be more accurately embodied, and

    (2) the fluids are allowed to move relative to one another under the
        influence of pressure gradients.

    It is of course the second of these features which will allow the
    acceleration of the flame propagation resulting from its
    semi-confinement to be properly represented.

    [Further applications to reacting-flow situations remain to be
    discussed]

    5.3 Applications to hydrodynamic and heat-transfer phenomena
    ------------------------------------------------------------

    The "MFM revolution" will not be complete until the reliance on an
    underlying hydrodynamic model, for example the k-epsilon model which
    has been used in several of the examples referred to above, has been
    totally dispensed with.

    Here it needs to be noted the k is nothing but a root-mean-square
    kinetic-energy fluctuation; and MFM is well able to compute RMS
    values, as has been seen above.

    All that is needed therefore is for velocity or energy to be made
    the discretized variable, and for something similar to be done about
    the length scale; then the underlying hydrodynamic model can be
    dropped.

    Work along these lines has already enabled simple shear flows to be
    simulated by MFM in a stand-alone manner (Spalding, 1996b). It will
    be reported elsewhere.

    It is flows with strong body-force effects which are likely to
    provide the greatest rewards; for single-fluid models simulate them
    least well. Therefore atmospheric and oceanographic flows, where
    gravitational effects are significant, and turbo-machinery flows,
    where centrifugal forces are large, may be the first to be subjected
    to MFM analysis.

    In order that MFM should contribute significantly to heat-transfer
    engineering, it will have to throw quantitative light on behaviour
    close to walls, where viscous effects become dominant. Nothing has
    been done to explore these effects, in the MFM context, so far.


    --------------------------------------------------------------------
    6. Research and development tasks
    ---------------------------------

    The multi-fluid model of turbulence, even though its essential ideas
    have been known for many years, has not yet been subjected to
    intense and prolonged research and development efforts such as have
    been lavished on the single-fluid models.

    Now is therefore the time to consider what form such efforts might
    take, and to initiate the appropriate actions.

    The author's current thoughts on these now follow.

    6.1 Mathematical and conputational tasks
    ----------------------------------------

    (a) Grid-refinement
    -------------------

    It is necessary to establish, for a much wider range of problems
    than has been investigated so far, that refining the population grid
    (ie increasing the number of fluids) does indeed always lead to a
    unique solution, and then to establish how the accuracy of solution
    depends upon the grid fineness.

    This needs to be done for two- and (at least a few) three-
    dimensional population grids as well as one-dimensional ones, before
    the soundness of the basic ideas and solution algorithms can be
    regarded as secure.

    (b) Consistency tests
    ---------------------

    When the micro-mixing constant becomes very large, the range of
    fluids present at a particular location reduces to 2 for a 1D
    population, to 4 for a 2D one and to 8 for a 3D one. It needs to be
    checked that the system of equations does indeed exhibit this
    behaviour, and that the right values are calculated for the fluid
    concentrations.

    A prudent and insightful investigator will wish to devise other
    tests of a self-consistency or conservativeness character before he
    wholly trusts his mathematical apparatus.

    (c) Non-uniform, unstructured, and self-adaptive grids
    ------------------------------------------------------

    So far, the author has used only uniform and structured grids, which
    have remained the same throughout the computation; yet there are
    obvious advantages to be gained by relaxing these restrictions.

    Thus, it may be convenient to define the attribute to be discretised
    as the ratio of the current temperature to the maximum temperature
    in the field, which temperature may change throughout the period of
    the calculation. Such a grid would have to be self-adaptive.

    Another kind of self-adaptive grid would change its uniformity as
    the calculation proceeded, so as to capture with maximum accuracy
    the shape of the fluid-population distribution.

    Yet another would insert new subdivisions of the discretised
    attribute in regions of steep variation, and perhaps remove them
    from regions of less interest.

    The best way of summarising the possibilities is to say that
    probably all the ingenious devices which specialists have invented
    for the better representation of variations in geometric space are
    likely to have their population-grid counterparts.

    (d) Computational improvements
    ------------------------------

    Although the computer-time burden is not yet great, it may become
    considerable when fine-grid three-dimensional transient simulations
    have to be undertaken.

    Fortunately, there is much which can be done to reduce the load. For
    example:-

    * population grids may be made non-uniform, self-adaptive and
      unstructured, as just explained;

    * the multi-fluid model can be confined to only those parts of the
      flow domain in which unmixedness is of physical significance;

    * cell-wise simultaneous-variable-adjustment procedures can be used,
      with multi-grid (in both senses) acceleration devices;

    * the employment of multi-block solution techniques.

    It will also be necessary to work out what is the optimal discretiz-
    ation principle; for certainly it is not optimal to employ a grid
    which is extremely fine in geometric space but has too small a
    population grid (ie too few fluids) to represent adequately the
    fragmentariness of the flow.


    6.2 Comparisons with experiment
    -------------------------------

    In order to establish the predictive capabilities of MFM now and in
    the future, it is of course necessary to make comparisons with
    experimental data.

    Three kinds of comparisons should be distinguished, namely:

    (1) those most desired by the potential end-users, for example
        between measured and predicted total smoke production for a
        particular gas-turbine combustor, or yields of main- and
        side-products for a paddle-stirred reactor;

    (2) those with well-known data on the velocity, temperature
        and concentration profiles in jets, wakes, boundary layers
        and simple turbulent flames; and

    (3) those which provide detailed measurements of probability-density
        functions of important fluid variables.


    Were there no urgency about solving the practical problems of
    engineering and the environment, a rational research program would
    probably concern itself with kind (2) first of all, in order to
    establish a convincing prima facie case for MFM; then experiments
    would be conducted systematically so as to permit kind-(2)
    comparisons, which might lead to improvements to some components of
    the model.

    Finally comparisons of kind (1) would be made, whereafter, if the
    comparisons were successful, industrial use of MFM would begin.

    However, perhaps fortunately, there IS much urgency; for the models
    which are currently used for predicting turbulent chemical reaction
    in particular are far from satisfactory. It therefore seems
    reasonable that some comparisons of kind (1) should begin in the
    near future.


    For example, it may be found that the Mendelian principle does not
    always apply; and it is certainly to be expected that, when the
    FPD's of two velocity components are in question, the coupling of
    two fragments with differing x-direction velocities will lead to
    offspring with differing y-direction velocities also. (See section
    6.3 below.)



    6.3 Conceptual developments
    ---------------------------

    (a) "Parental bias" in the "splitting" process
    ----------------------------------------------

    As has been mentioned above, the promiscuous-Mendelian hypothesis
    is unlikely to prove to be the best; and it is not hard to improve
    upon it, if the underlying mechanism of section 1.8b above is
    believed, as follows.

    The profile of concentration (say) in the conjoined fragments, after
    some time will have some such "error-function-like" shape as
    indicated below on the left, to which coresponds a pdf of the
    shape shown on the right.

    Clearly the pdf exhibits the largest frequencies near the
    extreme, ie the "parental" concentrations.

      |******                                ^        |******
      |      *****                           |        |*****
      |           ****                       |        |****
      |               ***                 conc-       |***
      |                  **               entr-       |**
      |                    *              ation       |*
      |                     **               |        |**
      |                       ***            |        |***
      |                          ****                 |****
      |                              *****            |*****
      |                                   ******      |******
      |-------------------------------------------    |---------------
            ----------- distance ------->               - frequency--->

    Only if the profile of concentration were linear with distance would
    the Mendelian assumption be correct.

    (b) The effects of laminar Prandtl and Schmidt numbers
    ------------------------------------------------------

    Let it now be supposed that the two fragments differ both in
    temperature and salinity. Then, while they are in contact, the
    profile of temperature will broaden much more rapidly than the
    profile of salinity.  As a consequence, the offspring do NOT lie on
    the diagonal joining the two parental locations as indicated in
    section 1.8c above.


            |_______|_______|_______|_______|_______|_______|_______|
            |       |       |       |       |       |       |  M    |
      ^     |       |       |       |       |       |       |   *   |
      |     |_______|_______|_______|_______|_______|_______|_______|
      |     |       |       |       |       |       |       |  *    |
    temp-   |       |       |       |       |       |       |       |
    erat-   |_______|_______|_______|_______|_______|_______|_*_____|
    ure     |       |       |       |       |       |       |       |
      |     |      *|       |       |       |       |       |       |
      |     |_______|_______|_______|_______|_______|_______|_______|
      |     |     * |       |       |       |       |       |       |
      |     |       |       |       |       |       |       |       |
      |     |____*__|____ __|_______|_______|_______|______ |_______|
            |  F    |       |       |       |       |       |       |
            |   *   |       |       |       |       |       |       |
            |_______|_______|_______|_______|_______|_______|_______|

                          ----- salinity------------>

    Instead, the offspring are more likely to lie along the lines of
    asterisks shown above, with much smaller salinity changes than
    temperature changes.

    (c) The coupling-splitting hypothesis for a two-velocity population
    -------------------------------------------------------------------

    As a final example of how the coupling-splitting hypothesis can be
    modified in the direction of greater realism, let it be imagined
    that the attributes of a 2D population grid are the horizontal
    velocity U and the vertical velocity V.

    Then let the collision be imagined of two fluid fragments which
    have the same values of V but differing values of U. So the father
    and mother lie on the same horizontal. Where do the offspring lie?


            |_______|_______|_______|_______|_______|_______|_______|
            |       |       |       |       |       |       |       |
      ^     |       |       |       |   *   |       |       |       |
      |     |_______|_______|_______|_______|_______|_______|_______|
      |     |       |       |       |       |       |       |       |
            |       |   *   |   *   |   *   |   *   |   *   |       |
      V     |_______|_______|_______|_______|_______|_______|_______|
            | F     |       |       |       |       |       | M     |
      |     |   *   |   *   |   *   |   *   |   *   |   *   |   *   |
      |     |_______|_______|_______|_______|_______|_______|_______|
      |     |       |       |       |       |       |       |       |
      |     |       |   *   |   *   |   *   |   *   |   *   |       |
      |     |_______|____ __|_______|_______|_______|______ |_______|
            |       |       |       |       |       |       |       |
            |       |       |       |   *   |       |       |       |
            |_______|_______|_______|_______|_______|_______|_______|

                          -----     U   ------------>

    The answer, it appears to the present author, must be "not only on
    the horizontal line FM"; for colliding fluid fragments are likely
    to generate motion in lateral directions as well as being checked or
    accelerated in the direction of their velocity difference.

    Therefore some of the offspring must be deposited into boxes above
    and below the horizontal line, of course in such a way as to
    preserve momentum and to ensure that there is at least no gain of
    energy.

    The above sketch illustrates this by its band of asterisks; but,
    before the idea can be expressed in a computer program, a precise
    offspring-distribution formula must be settled.

    (d) Other matters
    -----------------

    Among the questions still to be addressed are:-

    * How is the inter-fluid friction which opposes motion caused by
      differential body forces to be formulated?

    * As the Reynolds number becomes lower, fluctuations die out. This
      can be ensured by increasing the MIXRATE quantity. But what
      dependence of MIXRATE on Reynolds number will fit the experimental
      data?

    * In the mixing-layer model alluded to in section 3.2, the root-
      mean-square of the velocity fluctuations has been employed in a
      formula for the effective viscosity, in place of the usual k**0.5.
      However, MFM provides so much more information about the velocity
      pdf than is needed for the RMS value. In what way can that
      information be exploited so as to give more realistic simulations?

    No answers are proposed at the present time.









    --------------------------------------------------------------------
    7.   References
    ---------------

    MJ Andrews (1986) "Turbulent mixing by Rayleigh-Tay lor
         instability"; PhD Thesis, London University
    WH Avery & RW Hart (1953) Ind Engg Chem vol 45 p 16 34
    J  Boussinesq (1877) "Theorie de l'ecoulement tourbillant";
        Mem, Pre. par. div. Sav. Paris, vol 23
    P  Bradshaw, DR Ferriss & NP Atwell (1967) "Calculation of
        boundary-layer development using the turbulent energy equation";
        J Fluid Mech, vol 28, p 593
    KNC Bray and PA Libby "Counter-gradient diffusion in pre-mixed
         turbulent flames"; AIAA J vol 19, p205, 1981
    JY Chen & W Kollmann (1988) "PDF modelling of non-equilibrium
         effects in turbulent non-premixed hydrocarbon flames";
         22nd Int. Symp. on Combustion, Combustion Inst pp 645-653
    JY Chen & W Kollmann (1990) "Chemical models for PDF modelling
        of hydrogen-air non-premixed turbulent flames";
        Combustion and Flame, vol 79, pp 75-99
    SM Correa and SB Pope (1992) "Comparison of a Monte Carlo PDF/
        finite-volume model with bluff-body Raman data"
        Twenty-Fourth International Combustion Symposium
        The Combustion Institute, pp279-285
    RL Curl (1963) AIChE J vol 9, p 175
    BJ Daly & FH Harlow (1970) "Transport equations in turbulence";
         Phys Fluids, vol 13, p 2634
    C  Dopazo and EE O'Brien (1974) Acta Astronautica vol 1, p1239
    MA Elhadidy (1980),'Applications of a low-Reynolds-number
        turbulence model and wall functions for steady and unsteady
        heat-transfer computations', PhD Thesis, University of London
    DJ Freeman and DB Spalding (1995) "The multi-fluid turbulent
         combustion model and its application to the simulation of gas
         explosions"; The PHOENICS Journal (to be published)
    N  Fueyo (1992) "Two-fluid models of turbulence for
         axi-symmetrical jets and sprays"; PhD Thesis, London University
    N Fueyo, JC Larroya, L Valino, C Dopazo (1995) "A combined
        CFD-Montecarlo method for the solution of the scalar PDF
        equation in turbulent reaction"
        European PHOENICS User Conference, Trento, Italy
    M  Gonzalez & R Borghi (1991) "A Lagrangian intermittent model
         for turbulent combustion; theoretical basis and comparison with
         experiments", Turbulent Shear Flows 7, Springer Verlag
    FH Harlow & PI Nakayama (1968) "Transport of turbulence-energy
        decay rate"; Los Alamos Sci Lab U Calif report LA 3854
    NM Howe and CW Shipman (1965) "A tentative model for rates of
         combustion in confined turbulent flames"
         10th International Symposium on Combustion, p 1139
         The Combustion Institute, 1965
    ICOMP (1994) ICOMP-94-30; CMOTT-94-9; "Industry-wide workshop on
         computational modelling turbulence"; NASA Conference
         Publication 10165
    JO Ilegbusi and DB Spalding (1987) "A two-fluid model of
         turbulence and its application to near-wall flows"
         IJ PhysicoChemical Hydrodynamics , vol 9, pp 127-160
    JO Ilegbusi and DB Spalding (1987) "Application of a two-fluid
         model of turbulence to turbulent flows in conduits and shear
         layers" I J PhysicoChemical Hydrodynamics, vol 9, pp 161-181
    JH Kent & RW Bilger (1976) "The prediction of turbulent
         diffusion flame fields and nitric oxide formation"
         16th International Symposium on Combustion, The Combustion
         Institute p 1643
    SW Kim and CP Chen (1989), 'A multi-time-scale turbulence
         model based on variable partitioning of the turbulent kinetic
         energy spectrum', Numerical Heat Transfer, Part B Vol 16 pp193
    W  Kolbe & W Kollmann (1980) "Prediction of turbulent diffusion
         flames with a four-equation turbulence model"
         Acta Astronautica, vol 71 p 91
    AN Kolmogorov (1942) "Equations of motion of an incompressible
        turbulent fluid"; Izv Akad Nauk SSSR Ser Phys VI No 1-2, p56
    VR Kuznetsov (1979); USSR Fluid Dynamics vol 14, p328
    BE Launder and DB Spalding (1972)
        "Mathematical Models of Turbulence",  Academic Press
    FC Lockwood & AS Naguib (1975) "The prediction of the
         fluctuations in the properties of free, round-jet, turbulent,
         diffusion flames", Comb & Flame, vol 24 p 109
    JP Longwell (1954) "Selected combustion problems"
         Butterworths, p 508
    BF Magnussen and BH Hjertager (1976) "On mathematical modelling of
         turbulent combustion with special emphasis on soot formation
         and combustion"
         16th International Symposium on Combustion, pp 719-729
         The Combustion Institute
    MR Malin (1986) "Turbulence modelling for flow and heat
         transfer in jets, wakes and plumes"; PhD Thesis, London
         University
    HB Mason and DB Spalding (1973) "Prediction of reaction rates in
         turbulent pre-mixed boundary-layer flows"
         Combustion Inst European Symposium, pp 601-606
    JB Moss (1980) "Simultaneous measurements of concentration and
         velocity in an open pre-mixed flame", Combustion Science and
         Technology, vol 22, pp 115-129
    D  Naot, A Shavit M Wolfshtein (1974);"Numerical calculation of
         Reynolds stresses in a square duct with secondary flow"; Waerme
         u Stoffuebertragung, vol 7, p151
    M Noseir (1980) "Application of the ESCIMO theory of turbulent
         combustion"; PhD Thesis, London University
    KA Pericleous and NC Markatos (1991) "A two-fluid approach to
         the modelling of three-dimensional turbulent flames", in
         Proc. Eurotherm Seminar 17, Springer Verlag
    SB Pope (1982) Combustion Science and Technology vol 28, p131
    SB Pope (1985) Progr Energy Combust Sci  vol 11, pp119-192
    SB Pope (1990) "Computations of turbulent combustion; progress
         and challenges"
         Twenty-Third International Symposium on Combustion,
         The Combustion Institute, pp 591-612
    L  Prandtl (1925) "Bericht ueber Untersuchungen zur ausgebileten
        Turbulenz"; ZAMM vol 3, pp 136-139, 1925
    L  Prandtl (1945) "Ueber ein neues Formelsystem fuer die
        ausgebildete Turbulenz", Nachr. Akad. Wiss. Goettingen
    O  Reynolds (1874) "On the extent and action of the heating
        surface of steam boilers";
        Proc. Manchester Lit Phil Soc, vol 8, 1874
    RP Rhodes, PT Harsha & CE Peters (1974) "Turbulent-kinetic-
        energy analyses of hydrogen-air diffusion flame"
        Acta Astronautica vol 1 p 443
    R  Sala and DB Spalding (1973) "A mathematical model for an
        axi-symmetrical diffusion flame in a furnace". La Rivista di
        Combustibili, vol 27, pp 180-186
    PG Saffmann (1970) "A model for inhomogeneous turbulent flow";
        Proc Roy Soc London vol A317 pp 417-433
    S Shaanen, JH Ferziger, WC Reynolds (1975) "Numerical
         Simulation of turbulence in the presence of shear"; Dept Mech
         Engg Stanford U, Rep TF-6
    KI Shchelkin USSR J Tech Phys vol 13, p 520, 1943
    U Schumann (1973) "Ein Verfahren zur direkten numerischen
         Simulation ..." Thesis 1854 U Karlsruhe
    DB Spalding (1969) "The prediction of two-dimensional steady
        turbulent elliptic flows" ICHMT Seminar, Herceg Novi, Yugoslavia
    DB Spalding (1971a) "Mixing and chemical reaction in confined
         turbulent flames"
         13th International Symposium on Combustion, pp 649-657
         The Combustion Institute
    DB Spalding (1971b) "Concentration fluctuations in a round turbulent
         free jet"; J Chem Eng Sci, vol 26, p 95
    DB Spalding (1979) "The influences of laminar transport and
         chemical kinetics on the time-mean reaction rate in a
         turbulent flame";
         17th International Symposium on Combustion, pp431-440,
         The Combustion Institute.
    DB Spalding (1982) "Chemical reaction in turbulent fluids"
         IC CFDU Report 82/8, June 1982
    DB Spalding (1983) "Towards a two-fluid model of turbulent
         combustion in gases, with special reference to the spark-
         ignition engine" Inst Mech Eng, London, paper c53/83, 1983
    DB Spalding (1984)"The two-fluid model of turbulence applied to
         combustion phenomena"
         22nd AIAA Meeting, Reno, Nevada
    DB Spalding (1987) "A turbulence model for buoyant and combusting
         flows"; International J. for Numerical Methods in Engineering
         vol 24, pp 1-23
    DB Spalding (1993) Lecture delivered at "First International
         Conference on Air Pollution", Monterrey, Mexico. The text is
         included in "Lecture on the two-fluid turbulence model"
         supplied with the PHOENICS software package
    DB Spalding (1995a) "Models of turbulent combustion"
         Proc. 2nd Colloquium on Process Simulation, pp 1-15
         Helsinki University of Technology, Espoo, Finland
    DB Spalding (1995b) "Multi-fluid models of turbulent
         combustion"; CTAC95 Conference, Melbourne, Australia
    DB Spalding (1995c) "Multi-fluid models of Turbulence",
       European PHOENICS User Conference, Trento, Italy
    DB Spalding (1996a) "Older and newer approaches to the numerical
         modelling of turbulent combustion". Keynote address at 3rd
         International Conference on COMPUTERS IN RECIPROCATING
         ENGINES AND GAS TURBINES, 9-10 January, 1996, IMechE, London
    DB Spalding (1996b) "Multi-fluid models of Turbulence; Progress and
         Prospects; lecture to be presented at CFD 96, the Fourth Annual
         Conference of the CFD Society of Canada, June 2 - 6, 1996,
         Ottawa, Ontario, Canada
    DB Spalding (1996c) "Progress report on the development of a multi-
         fluid model of turbulence and its application to the paddle-
         stirred mixer/reactor", invited lecture at 3rd Colloquium on
         Process Simulation, Espoo, Finland, June 12-14
    DB Spalding (1966d) Unpublished work at CHAM Ltd
    RLT Sun (1982) "Application of the ESCIMO theory to turbulent
         diffusion flames"; PhD Thesis, London University
    LT Tam (1981) "The theory of turbulent flow with complex
         chemical kinetics"; PhD Thesis, London University
    ST Xi (1986) "Transient turbulent jets of miscible and
         immiscible fluids"; PhD Thesis, London University
    V  Yakhot & S Orszag (1986) "Renormalization group analysis of
       Turbulence. I. Basic theory"
       J. Sci Comp, vol 1, no 1, pp 1-51
    DC Wilcox (1993) "Turbulence modelling for CFD", DCW Industries,
        La Canada, California
    GC Williams, HC Hottel and AC Scurlock,(1949), 3rd Symposium on
        Combustion, p21, Williams and Wilkins
    WJ Wohlenberg "Minimum depth for flame front for stable
         combustion ... in a gaseous system at constant pressure"
         4th International Symposium on Combustion p 796
         Williams and Wilkins, Baltimore, 1953
    JZ Wu (1987) "The application of the two-fluid model of
         turbulence to ducted flames" IC CFDU Report, June 1987