2. MATHEMATICAL FORMULATION

Contents

2.1 Differential equation for the conservation of a fluid
2.2 The source-sink term
2.3 The differential equation for a non-discretised attribute
2.4 Population-average quantities
2.5 Computational aspects
2.6 Auxiliary relationships
2.7 The length-scale equation

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.

Such attributes may be referred to as "population-distinguishing attributes", abbreviated to PDAs.

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

d_ma_d_t { time-dependence }
+ v_vec * grad ma { convection }
=
spvol * div (G * grad ma) { diffusion }
+ spvol * ma_dot { source of mass} _(Eq. 1)

where:

• d_ma_d_t stands for the partial derivative of ma with respect to time, t;
• v_vec is the velocity vector of the mixture,
• spvol is the specified volume of the mixture, i.e. the reciprocal of its density;
• G (the abbreviation of 'gamma') is the turbulent exchange coefficient of the mixture;
• ma_dot is the net rate of creation of the fluid associated with a;
• div and grad stand for 'divergence' and 'gradient' respectively;
• * signifies multiplication: and
• the dimensions of all terms are 1/time.

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

If however different values of v_vec and G are appropriate for different fluids (as, for example, where fluid is denser than another and therefore moves at a different velocity), such fluid-specific features will be expressed by modification of the source-sink term, ma_dot.

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 ma_dot are those which cause the fluid with attribute a to be created or destroyed. These are of two main kinds, namely:

1. exchange of mass between fluids;
2. translational movement of fluid in population space.

   These will be discussed in sections (b) and (c) below.
   Also to be considered are:

3. within-fluid sources and sinks.

   These will be discussed in the subsequent section 2.3.

(b) 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:

1. at what rate fragment-types A and E "couple" with each other; and

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

The first of these can be computed in various ways, of which two commonly-used ones are as follows:

• The rate per unit volume is proportional to the product of the mass fractions of the fluids in question, i.e. to:
  ma * me
  times
  epsilon / k,
  the reciprocal time constant of energy decay;

  or to:

• the same mass-fraction product times a quantity having the same dimensions as epsilon / k, derived from the population distribution itself.

The second, i.e. the mass-redistribution process, depends on the manner in which the individual fragments of fluid are supposed to encounter each other, and exchange mass and other attributes.

A simple example,the "promiscuous-Mendelian" hypothesis, has been mentioned in section 1.8 (a) above.

Another can be derived by analysing the time-dependent diffusion process which ensues when two semi-infinite bodies of material are placed suddenly in contact and remain so for a finite time.

Descriptions of this 'brief-encounter' process are to be found in a lecture on the relationship between MFM and the 'flamelet model' of turbulent combustion, and especially here and in associated illustrations concerned with momentum transfer and reactedness transfer.

The lecture in question contains an explanation of how epsilon / k is such a good indicator of the rate of micro-mixing.

(c) Movement in population space

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_t 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 conceptual difficulty 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 Tc_t 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) / Tc_t.

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

  mc_dot = - mc * Tc_t / (Td-Tc) _(Eq. 2)

  simultaneous with a source of D given by

  md_dot = + mc * Tc_t / (Td-Tc) _(Eq. 3)

• 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 concentration of chemical species, the rate of creation by chemical reaction;

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

• In all cases the m_dot expressions are of the form:

  source / attribute difference.

The following remarks concerning the choice of velocity as a PDA may be of interest at this point:

• Such a choice is appropriate when one is seeking a more profound insight into hydrodynamic turbulence than can be provided by a 'Kolmogorov-type' turbulence model, in which attention is focussed on statistical properties of the fluid mixture, for example k, epsilon or the 'Reynolds-stress' components.

• One may choose to discretise just one, or two or all three components of velocity. Even the first choice provides novel insigts, as is exemplified by the plane-mixing-layer study reported here and now made accessible as PHOENICS-Input-Library case L300.

• If only one velocity is discretised, turbulent-transport processes may still have to be modelled by some kind of 'effective-viscosity hypothesis'.
  However, if two components are discretised, for example u. the main-flow-direction (x) velocity, and v, the cross-stream-direction (y) velocity, a two-dimensional population results, from knowledge of which such important statistical properties as [u'*v'] (with [ ] indicating 'average for the whole population') can be derived, and the prime indicating 'difference from the local mean'.

• Of course, many variants are possible. One such is the discretization not of u and v, but of:
  • u', i.e. u - [u] , and
  • v', i.e. v - [v]
  i.e. the velocity differences from the local mean.

  By definition [u'] and [v'] are zero, which simplifies certain mathematical manipulations.

• At the time of writing (January 2003), almost no research into these modelling possibilities has been carried out.

2.3 The differential equation for a non-discretised fluid attribute

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 can be distinguished from PDAs by calling them "continuously-varying (i.e. not-discretised) attributes", abbreviated to CVAs.

These too need to be computed; and a differential equation must be solved for each.

Let cv1b stand for one of the CVAs of fluid B. Then the differential equation governing its variation with space and time is:

d_(mb*cv1b)_d_t { time-dependence }
+ v_vec * grad (mb*cv1b) { convection }
=
spvol * div (G * grad (mb*cv1b)) { diffusion }
+ spvol * (mb*cv1b)_dot { source of CVA1} _(Eq. 4)

The terms on the left-hand side of the equation require little comment; for they amount to saying that, in the absence of right-hand-side terms, the amount of property cv1 per unit mass of fluid B, namely cv1b, 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 contents of the bracket as follows:

grad (mb*cv1b) = mb * grad cv1b + cv1b * grad mb _(Eq. 5)

Discussion of Eq. 5

1. The first term on the right stands for the within-fluid-B diffusion of the property cv1, while

2. the second term implies that any of fluid B which changes location as a result of random turbulent motion carries its cv1 property with it.

3. The term (mb*cv1b)_dot 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

   • a second part representing the consequences of the coupling-and-splitting process.

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

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

   mb_dot,A,D * (cv1a + cv1d) / 2

5. 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.

6. 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 cv1b.

   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 = (ma*Ta + mb*Tb + mc*Tc + ....... + mz*Tz) _(Eq. 6)

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

cv1az = (ma*cv1a + mb*cv1b + mc*cv1c + .........+ mz*mcv1z _(Eq. 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:

Taz(n) = (ma Ta**n + mb*Tb**n + mc*Tc**n + ....... mz*Tz**n) _(Eq. 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 Eq. 1 or Eq. 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 (such as ma and cv1b), 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.6 Proportionality 'constants'

CONMIX

1. epsilon / k, or to
2. a rate quantity derived from the population distribution itself, for example the root-mean-square value of the velocity fluctuations, divided by some reference length.

In various implementations of MFM in PHOENICS, the name CONMIX has been used for this constant.

When epsilon / k is the rate term to be used, there are theoretical reasons for supposing that, in gases, the value of CONMIX is greater than 1.0 but probably not by much more than a factor of 10.

When some other rate term is used, of course, the value of CONMIX must not be expected to be the same. Experimental evidence must, in general, be sought for its determination.

In the mixing-layer example of the input-file library case L300, in which the rate is taken as the RMS value of the longitudinal-velocity fluctuations divided by a length scale, the value chosen for CONMIX, which fits the experimental data quite well, is 5.0.

VISCON

A typical example is MFM Input-file library case L300, which:

1. is concerned with the plane mixing layer, and
2. involves the postulate that the turbulent exchange coefficient (viz G in Eq. 1 above) is equal to:

   VISCON * RATE * MIXL ** 2

   where:

   1. MIXL represents the average length scale of the turbulent fragments, reminiscent of Prandtl's 'mixing length', and
   2. RATE stands for the local root-mean square value of the normalized.velocity fluctuations, divided by the length scale..

The value to be used for VISCON must depend upon how MIIXL and MNSQ are defined and derived.

Thus, the MIXL value used may actually be the Prandtl Mixing Length; or it may be derived from MFM's own length-scale equation, discussion of which now follows, and which involves its own dimensionless proportionality constant, LENCON.

2.7 The length-scale equation

In principle it can generate its own length scale from the postulate that the fluid fragments which constitute the turbulent mixture form a population of which size is the distinguishing attribute. However, this has not yet been done,

Of the many which can be imagined, only one length-scale equation has been investigated so far, namely that which is exemplified by the just-mentioned library case L300.

In this case, MIXL is supposed to be the dependent variable of a differential scalar-transport equation of which the only distinguishing features are:

1. a source term equal to:

   LENCON * MNSQ * POPMAX , where:

   • MNSQ stands for the normalised root-mean-square velocity fluctuations, and
   • POPMAX is the maximum value of the PDA of the population, i.e. the maximum longitudinal velocity;

   and

2. a rather low turbulent Prandtl number, namely 0.1

The source term implies that the length scale can only grow, except in so far as it is affected by communication with lower-length-scale boundary regions.

The low value of Prt is intended to ensure that non-uiformities of length scale across the mixing layer are small.

The invention and testing of alternative hypotheses remains as an interesting but unexplored research opportunity.