- General introduction
- IPSA in PHOENICS
- Interphase transport processes
- Interphase-Source Term
- Interface Value- PHINT
- Interphase Transfer Coefficient - FIP
- Reference Interphase Transport Coefficient - CFIPS
- Interphase Friction - Built-In Options
- Interphase Heat Transfer - CINT
- Interphase Heat Transfer - Built-In Options
- Interphase Mass Transfer - CMDOT
- Interphase Mass Transfer - Built-In Options
- "Virtual-mass" effects
- Introduction of Non-Standard Correlations

- Particle-Size Calculation
- Encyclopaedia and Help Entries
- Library Examples
- References

Multiphase flows are of great practical importance in many common engineering and environmental applications:

- Steam generators and condensers, steam turbines
- Coal fired furnaces
- Fluidized bed reactors
- Liquid sprays
- Separation of contaminants from a carrier fluid
- Free surface flows, where sharp interfaces exist

Multiphase flows are characterized by two or more fluids in motion relative to each other.

The fluids will also usually have different physical properties - temperature, density, conductivity.

A phase is an identifiable class of material in a mixture. Identifying characteristics may be:

- thermodynamic phase (solid, liquid or gas)
- size (e.g., all particles between 200 and 400 microns)
- geometric shape (e.g., all round particles)
- velocity (e.g., all upward moving particles)
- temperature (e.g., all particles above 500K)
- The most common identifier is the first, based on thermodynamic phase.

Flows where the phases are not mixed, but are separated by a sharp interface, can be also solved using multiphase techniques.

It is often more efficient to use modified single phase methods, such as the Scalar Equation or Height of Liquid methods.

The rest of this lecture concentrates on interspersed flows.

The prediction of multiphase phenomena involves computation of the values of:

- up to 3 velocity components for each phase,
**u**_{i} - 1 volume fraction for each phase,
**r**_{i}

and possibly:

- temperature
- chemical composition
- particle size
- turbulence quantities
- pressure

for each phase.

Specific features of the solution procedure are:

- Eulerian-Eulerian techniques using a fixed grid, and employing the concept of 'interpenetrating continua' to solve a complete set of equations for each phase present;
- The volume fraction, Ri, of phase is computed as the proportion of volumetric space occupied by a phase;
- It can also be interpreted as the probability of finding phase
at the point and instant in question;*i* - All volume fractions must sum to unity;

This method is embodied in PHOENICS as the Interphase Slip Algorithm, IPSA;

An alternative technique is to track individual particles in a Lagrangian framework. This technique, embodied in GENTRA, is beyond the scope of this lecture.

Each phase is regarded as having its own distinct velocity components.

Phase velocities are linked by interphase momentum transfer - droplet drag, film surface friction etc.

Each phase may have its own temperature, enthalpy, and mass fraction of chemical species.

Phase temperatures are linked by interphase heat transfer.

Phase concentrations are linked by interphase mass transfer.

Each phase can be characterized by a 'fragment size'. This could be a droplet or bubble diameter, film thickness or volume/surface area.

Phase 'fragment sizes' are influenced by mass transfer, coalescence, disruption, stretching etc.

Each phase may have its own pressure - surface tension raises the pressure inside bubbles, and interparticle forces prevent tight packing, by raising pressure.

The task is to provide equations, from the solution of which values of r_{i}, u_{i},
v_{i}, w_{i}, T_{i}, C_{i}, and so on can be deduced.

These equations include mathematical expressions for the rates of the interphase transport processes described above. These are sometimes called the 'constitutive relationships'. The expressions are often empirical, and mainly come from experiment.

The equations describing the state of a phase are basically the Navier-Stokes Equations, generalised to allow for the facts that:

- Each of the phases occupies only a part of the space, given by the volume fraction; and
- The phases are exchanging mass and all other properties.

The phase continuity is regarded as the equations governing the phase volume fractions. They are the counterparts of the single phase continuity equation.

d(r_{i}r_{i})/dt + div( r_{i}r_{i}V_{i} - G_{ri}grad(r_{i})
) = ñ_{ji}

__Transient__ __Convection__ __Phase Diffusion__ __Mass
Source__

Where:

- r
_{i}= phase volume fraction, m^{3}/m^{3} - r
_{i}= phase density, kg/m^{3} - V
_{i}= phase velocity vector, m/s - G
_{ri}= phase diffusion coefficient, Ns/m^{2} - ñ
_{ji}= net rate of mass entering phase i from phase j, kg/(m^{3}s)

Notes:

- The phase diffusion term models the turbulent dispersion of particles by random motion mechanism. It is not present in laminar flows.
- As the phases completely fill the available space, the volume fractions sum to unity:
r

_{1}+ r_{2}= 1

The conservation equations for any variable of phase i, f _{i},
read:

d(r_{i}r_{i}f_{i})/dt+div(
r_{i}r_{i}V_{i}f_{i}
-r_{i}G_{f}_{i} grad(f_{i})-f_{i}G_{ri} grad(r_{i}) ) = S_{i}+S_{ip}

__Transient__ __Convection__ __Within-phase__ and __Phase
Diffusion__ __Sources__

Where:

- r
_{i}= phase volume fraction, m^{3}/m^{3} - r
_{i}= phase density, kg/m^{3} - V
_{i}= phase velocity vector, m/s - G
_{f}i = within-phase diffusion coefficient, Ns/m^{2} - G
_{ri}= phase diffusion coefficient, Ns/m^{2} - S
_{i}= within-phase volumetric sources, kg×f/(m^{3}s) - S
_{ip}= interphase volumetric sources, kg×f/(m^{3}s)

Notes:

- The volume fraction multiplier allows for the 'diluting' effect of the other phase.
- The within-phase diffusion term represents the molecular and turbulent mixing present in the phase.
- The phase diffusion term represents the transport of f brought about by the turbulent dispersion of the phase itself.

IPSA has been implemented in PHOENICS for two phases only.

- The treatment in EARTH is impartial between the phases - if either phase volume fraction goes to 1.0, the equation for that phase reduces to the single phase form.
- The transient, convective and diffusive terms all contain the appropriate volume fraction multiplier, upwinded or averaged as required.
- The links between the phases - mass, momentum and heat transfer - are introduced via an interphase source.

IPSA is activated by placing in the Q1 file the statement:

ONEPHS = F (The default value is T)

Alternatively, in the VR-Editor, IPSA is activated from the Models panel of the Main Menu. The IPSA Settings panel allows all the control variables described below to be set. Solution of the required phase 1 and phase 2 variables is then automatic.

The Inlet- and Outlet-object attribute dialogs offer the choice of which phase is to enter or leave.

followed by:

- Solution of the required phase-1 and phase-2 variables is activated by the SOLVE
command:
SOLVE(P1, U1,U2, V1,V2, W1,W2, R1,R2, H1,H2, C1,C2, ...)

- Pressure (P1) is common to both phases. There is no provision for SOLVEing P2. However, if P2 is STOREd, it will be added to the pressure force acting on phase 2 only. It is left to the user to introduce GROUND coding which calculates values of P2. Some examples are present in GREX.
- Because R1 + R2 = 1, it is possible to SOLVE one of R1 or R2, and only STORE the other. This works well if both phases are present in finite quantities throughout the domain. If one of the phases can disappear, it is better to SOLVE both R1 and R2 lest divisions by zero occur.
- Second-phase density is set via RHO2, with the same built-in options as for RHO1.

A numerical technique called the Partial Elimination Algorithm (PEA), is used to procure convergence for tightly-linked variables. By default, the PEA is active between pairs of SOLVEd variables, as shown below:

Phase 1 Variable Phase 2 Variable

- U1 - - - - - - - - U2
- V1 - - - - - - - - V2
- W1 - - - - - - - W2
- H1 - - - - - - - - H2
- C1 - - - - - - - - C2
- C3 - - - - - - - - C4
- C5 - - - - - - - - C6
- etc.

Note that all odd numbered variables 'belong' to phase 1, and even numbered variables 'belong' to phase 2.

The TERMS command can be used to switch variables between phases.

Argument 6 is 'Variable belongs to first phase?'.

Thus:

TERMS(C2,P,P,P,P,Y,P)

would switch C2 from its default phase 2 to phase 1.

**Inlets**

Boundary conditions must be specified for whichever phase is to be affected by the PATCH in question. For mass flow conditions, the practice is to use P1 for phase 1 mass fluxes, and P2 for phase 2 mass fluxes.

A typical inlet specification where both phases are entering the domain is:

INLET(IN, LOW, 1, NX, 1, NY, 1, 1, 1, LSTEP) VALUE(IN, P1, R1in * RHO1in* W1in) VALUE(IN, P2, R2in * RHO2in * W2in) VALUE(IN, W1, W1in) VALUE(IN, W2, W2in)

wherein R1in, R2in are the volume fractions in the incoming stream.

Obviously, R1in + R2in = 1.0.

*Outlets*

If both phase densities are constant, and both phases can pass in or out, outflows can be simply specified using:

OUTLET(OUT, HIGH, 1, NX, 1, NY, NZ, NZ, 1, LSTEP) VALUE(OUT, P1, Pext) VALUE(OUT, P2, Pext)

- if Pext is not zero.

If only one phase is allowed to pass, then a PATCH must be used, with a COVAL for P1 or P2.

PATCH(OUT, HIGH, 1, NX, 1, NY, NZ, NZ, 1, LSTEP) COVAL(OUT, P1,1E3, Pext)

- if only Phase 1 is allowed to pass

or

COVAL(OUT, P2,1E3, Pext)

- if only Phase 2 is allowed to pass.

*Walls*

The WALL command assumes that phase 1 is the continuous phase, and will only generate wall functions for this phase. If wall effects are required for phase 2, these have to be explicitly added in.

Note that in the examples above, P2 is only a flag for phase 2 mass sources - it is NOT a 'real' pressure, and it does NOT have to be STOREd or SOLVEd.

For single phase flows, the mass flow associated with a PATCH and COVAL for P1 is:

m' = Type × C_{p1} × ( V_{p1} - P1 )

where Type is the geometrical multiplier specified through the PATCH type.

In two phase flows, the relationship for P1 (or P2) is:

m' = Type × C_{p1} × ( V_{p1} - P1 ) for P1 < V_{Pi}
{inflows}

m' = Type × R_{1} × C_{p1} × ( V_{p1} - P1 ) for P1 > V_{Pi}
{outflows}

i.e., it is the same for inflows, but has an additional volume fraction multiplier for outflows. This means that for equal phase COefficients and VALues, outflows are proportional to in-cell volume fraction.

Unfortunately, in many practical cases, outflow is more realistically represented as proportional to in-cell mass fraction. This can be achieved by setting the COefficient to the required COefficient multiplied by the in-cell phase density.

The OUTLET command thus arranges that the COVALs generated are:

COVAL(OUT,P1,1E3,Pext) or COVAL(OUT,P1,RHO1*1E3,Pext)

COVAL(OUT,P2,(RHO2/RHO1)*1E3,Pext) or COVAL(OUT,P2,RHO2*1E3,Pext)

If either RHO is set to a GRND number, both COefficients are set to 1E3, and it is left to the user to arrange for appropriate multiplication, either by average outlet densities in Q1, or by the actual in-cell densities in GROUND.

The normal within-phase diffusion treatment for two phase flow is very similar to the single phase. The diffusion coefficient for each phase is written as:

G_{f}i = RHO_{i} × ( ENUL / PRNDTL(f) + ENUT / PRT(f) )

The term includes both laminar and turbulent mixing. Note that there is only one laminar and one turbulent viscosity. These are the phase 1 viscosities. If different viscosities are required in the second phase, the Prandtl numbers can be used to introduce the appropriate ratios.

If one of the phases, normally phase 2, is a dispersed one, say droplets or particles, then within-phase turbulent effects are absent. PRT(f) should be set to a large number, say 1E10.

The phase diffusion term accounts for turbulent dispersion of particles or droplets due to turbulence in the continuous phase, and is present in both the phase continuity and phase phi equations. The diffusion coefficient is written as:

G_{ri} = RHO_{i} × ( ENUL / PRNDTL(R_{i}) + ENUT / PRT(R_{i})
)

Phase diffusion effects can be removed by setting both Prandtl numbers for both phases to large numbers, say 1E10.

The turbulence variables, KE and EP, are always phase-1 variables, and the generation rate is calculated from phase-1 velocity gradients. If required, it is possible to modify the GREX coding to switch the turbulence model to phase 2, though this should not normally be necessary.

The presence of particles will normally have a damping effect on the turbulence field. Strictly speaking, this effect is not completely modeled in the built-in k-E model. Only the generation term and dissipation terms are reduced by the volume fraction, but there is no explicit damping. This can be introduced via additional source/sink terms.

Two forms of these terms are provided in the GX routine GXDISP, which is activated by a PATCH whose name starts with 'KEDI'. COef and VALue for KE and EP are set to GRND1,0.0 or GRND2,0.0 to activate the models of Mostafa & Mongia or Chen & Wood respectively.

The interphase source contains the diffusive (e.g., friction, heat transfer) and convective (mass transfer related) links between the phases. It is formulated in the general form as follows.

*General Form *

S_{IP}=(f_{f,i}+ <m_{ji}>) (f_{i}^{int} -f_{i})

where:

f_{f,i} = interphase (diffusive) transfer
coefficient, kg/s;

m_{ji} = net mass transfer rate between the phases, kg/s;

<> = maximum of 0.0 and the quantity enclosed;

f_{i}^{int} = value of f
at the interface between the phases and

f_{i} is a bulk phase value of conserved variable f.

*Notes: *

- The unit of S
_{IP}is (kg/s) × (unit of f),i.e.- if f is velocity, m/s, the units of S
_{IP}are Newtons and - if f is enthalpy, J/kg, the units of S
_{IP}are Watts.

- if f is velocity, m/s, the units of S
- S
_{IP}, obviously, appears in the equations for each of paired variables, f_{1}and f_{2},i.eS

_{IP,1}=(f_{f,1}+ <m_{21}>) (f_{1}^{int}-f_{1})

S_{IP,2}=(f_{f,2}+ <m_{12}>) (f_{2}^{int}-f_{2})wherein m

_{21}= - m_{12}; - The last argument of the TERMS command determines whether the interphase source Sip is active between a pair of variables. Pairing is between odd and even numbered variables.

*The concept of the interface value*

f_{i}^{int} is the value of f_{i} at the i-j interface. It can be a function of space,
time, or local bulk phase values. It is a property, and not a variable obtained from a
conservation equation.

In terms of a steam - water system, PHINT(H1) and PHINT(H2) would represent the saturation enthalpies of steam and water at the local temperature and pressure. The near-interface situation is shown on the picture above.

In the case where CO2 diffuses between water and air, and C1 represents the concentration in the air, and C2 represents the concentration in the water, then PHINT(C1) and PHINT(C2) would represent the concentrations at the air-water interface.

In other cases, interphase momentum transfer for example, the concept of an interface value is not helpful, and it is more meaningful to think of direct transfer between bulk phase values.

*Settings for PHINT(f)*

The PIL variable controlling f_{i} is PHINT(f). There are many options for setting PHINT values.

€ *Bulk-to-bulk-transfer interface values*

In this case, the interface value of current phase is made equal to the bulk value of another phase, i.e.

f_{i}^{int} = f_{j}

Then, interphase source becomes

S_{IP}=(f_{f,i}+ <m_{ji}>) (f_{j} -f_{i})

where:

f_{i} is current phase bulk value of f and

f_{j} is another phase bulk value of f.

- Usage
- Recommended practice for momentum transfer;
- For energy transfer via overall interphase heat transfer coefficients.

- Activation
- No actions required, i.e.
- PHINT(f
_{1}) = default and PHINT(f_{2}) = default.

€ *Known interface values*

If the values of f at the interface between the phases are known, say,

f_{i}^{int} = Value_{1} and

f_{j}^{int} = Value_{2}

interphase source becomes

S_{IP}=(f_{f,i}+ <m_{ji}>)
(Value_{i}-f_{i})

- Usage
- Whenever the interface values are thought to be the given constants, as, for example,
- saturation interface enthalpies of steam and water at the prevailing pressure.

- Activation
- PHINT(f
_{1}) = Value_{1}and - PHINT(f
_{2}) = Value_{2}

- PHINT(f

Notes:

- A few built-in option are provided for H1 and H2 to deal with pressure-dependent saturation enthalpies at the interface. They are described in PHINT entry of Encyclopaedia;
- User-specified formulae for interface value calculations are accepted via GROUND
provided that PHINT(f
_{i})=GRNDn, as appropriate.

€ *Known difference in interfacial values.*

If the difference in values of f is known at the interface between the phases, say,

diff = f_{i}^{int} - f_{j}^{int}

the actual interface values can then be expressed via phase bulk values as

f_{i}^{int} = f_{j}
+ diff

f_{j}^{int} = f_{i}
- diff

and resulting interphase source becomes

S_{IP}=(f_{f,i}+ <m_{ji}>) (f_{j} - f_{i} + diff)

- Usage
- Whenever there is a "jump" conditions at the interface as, for example,
- a constant heat of evaporation linking the steam and water saturation enthalpies.

- Activation
- PHINT(f
_{1}) = default and - PHINT(f
_{2}) = diff

- PHINT(f

€ *Known ratio of interfacial values.*

If the ratio of f values is known at the interface between the phases, say,

Value_{1} = f_{i}^{int}/ f_{j}^{int}

the actual interface values can then be expressed via phase bulk values as

f_{i}^{int} = Value_{1} × f_{j}

f_{j}^{int} = f_{i}/Value_{1}

and resulting interphase sources become

where:

C_{1} and C_{2} are bulk-to interface transfer coefficients to be
defined in the following considerations.

- Usage
- Whenever there is a "floating-interface" condition as, for example,
- the ratio of specific heats for introduction of interphase heat transfer via enthalpy difference.

- Activation
- PHINT(f
_{1}) = Value_{1}and - PHINT(f
_{2}) = default

- PHINT(f

f_{f,i}, kg/s, featured in S_{IP} is
actually the interphase flux per cell per unit conserved variable difference.

For built-in treatment, the specific interphase fluxes are made proportional to f_{IP},
which is termed the reference interphase transfer coefficient.

**f _{IP}** or

It appears as multiplier, in the built-in formulations of all specific f_{f,i} and m_{ji}. This multiplier is introduced because
of the analogy between friction, heat and mass transfer.

The coefficient f_{IP} is controlled by the PIL variable CFIPS.

If CFIPS is any constant other than a GRND flag, f_{IP} is set as follows:

f_{IP} = CFIPS×RHO1×R1× R2^{*}×Vol, for CFIPS >0.0

f_{IP} = |CFIPS|×RHO2×R2× R1^{*}×Vol, for CFIPS <0.0

In the above, Ri^{*} is the maximum of Ri and RLOLIM. This ensures that, as
long as RLOLIM > 0.0 , the reference interphase transport coefficient remains finite,
even though the volume fraction of either phase falls to zero. Vol is the cell volume.

This has reasonable consequences when friction is concerned, as it ensures that very small bubbles travel at the same velocity as the carrying liquid.

Ignoring the interphase mass transfer, for bulk-to-bulk transport of momentum, the interphase source term becomes:

S_{IP} =f_{IP} × ( vel_{j} - vel_{i})

where S_{IP} has units of Newton and f_{IP} is now the interphase drag
coefficient between phases in units of Ns/m or kg/s. CFIPS of the above options,in units
1/s, is now proportional of the product of phase slip velocity and projected area per unit
volume.

If STORE(CFIP) is set in the Q1 file, then f_{IP} is stored for output
purposes, and also for under-relaxation purposes via the RELAX command.

The built-in options provided are described fully under the Encyclopaedia entries CFIPS and INTERPHASE DRAG MODELS , and so only a brief description of the major ones will be provided here.

- CFIPS=GRND7 selects the DISPERSED-FLOW drag model:
- For CFIPD=4 or 6 the surface tension must also be defined via CFIPC.
- If CFIPB = -Dp, the minus sign activates the removal of R1
^{*}from the f_{IP}formula; a form which is often used for a dilute suspension of particles. - CFIPD=7 selects the particle-fluidization drag model. This model uses a f
_{IP}formula of the above form only when R1 > 0.8, otherwise it uses a formula based on the well-known Ergun correlation for packed beds. - CFIPS=GRND8 selects the same drag models as for CFIPS=GRND7, but with phase 2 as carrier
and phase 1 dispersed, i.e.:
f

_{IP}= 0.75×RHO2×R1×R2^{*}×Vol× <V_{slip},CFIPA> / D_{p} - The commands STORE(SIZE,REYN,VREL,CD,APRJ,EOTV,WEB) will allocate 3D storage for the particulate phase diameter, its Reynolds number, relative slip velocity, drag coefficient, projected area, Eotvos and Weber numbers, when appropriate.
- The options are restricted to rigid particles; they do not allow for particles, droplets or bubbles to change shape, and hence influence the drag.
- The options do not account entirely for multi-particle systems, e.g., the continuous-phase viscosity is used in evaluating Re rather than the apparent mixture viscosity.
- The options are restricted to spherical shapes, but can be extended to irregular-shaped particulates by interpreting Dp as an effective sphere diameter or by using a different drag correlation.
- The options do not account for compressibility effects, in which case Cd becomes a function of both Mach number and Reynolds number.

f_{IP} =0.75 ×C_{d}×RHO1×R2×R1^{*} ×Vol×<V_{slip},CFIPA>
/ D_{p}

wherein phase 1 is the carrier and phase 2 is dispersed.

Here, CFIPA is a parameter used to define the minimum value allowed for slip velocity ,
V_{slip}, between the phases, Cd is a dimensionless particle drag coefficient and
Dp is the particle diameter (defined by CFIPB).

The PIL variable CFIPD selects the Cd correlation:

= 0 - Standard drag curve (default)

= 1 - Stokes-regime drag correlation

= 2 - Turbulent-regime drag correlation

= 3 - Subcritical-regime drag correlation

= 4 - Distorted-bubble "dirty-water" drag correlation

= 5 - Spherical-bubble "dirty-water" drag correlation

= 6 - Ellipsoidal-bubble "clean-water" drag correlation

Notes:

Notes:

In many cases, the driving force is the temperature difference and available heat transfer correlations are based upon it.

With solving for enthalpies, H1 and H2 and ignoring mass transfer the interphase sources become:

S_{H1} = h_{12} × A_{S} × (T_{1}^{int} - T_{1})

S_{H2} = h_{21} × A_{S} × (T_{2}^{int} - T_{2})

where:

h_{12} = bulk1-to-interface heat transfer coefficient, W/m^{2}K

h_{21} = bulk2-to-interface heat transfer coefficient, W/m^{2}K

A_{S} = total interface area, m^{2}

The interface temperatures, T_{1}^{int} and T_{2}^{int},
can be eliminated via overall heat transfer coefficient and bulk-to-bulk temperature
difference as follows:

S_{Hi} = 1 × A_{S} × (T_{j} - T_{i})/(1/C_{1}
- 1/C_{2})

where: C_{1} = h_{ij} and C_{2} = h_{ji}

This formulation is generalised so that the interphase transfer coefficient, f_{IP},
is taken to be harmonic average of the bulk-to-interface coefficients, C_{1} and C_{2}.

Thus,

f_{IP} = 2 / (1/C_{1} - 1/C_{2})

Note: if C_{1} = C_{2} = C, then f_{IP} = C

Built-in provisions have been made to calculate C_{1} and C_{2} via
their ratios to the reference transfer coefficients, as follows

C_{1} = CINT(f_{1}) × FIP

C_{2} = CINT(f_{2}) × FIP

Notes:

- The default values of CINT() are 1.
- If either CINT() is set to a GRND flag, there is no multiplication by FIP for that
phase, and the GROUND set value is used for C
_{i}directly. - Two options are provided internally for CINT(H1) and CINT(H2), see Encyclopaedia entry CINT for details.

In PHOENICS two-phase mode, it is common to solve for enthalpy. It is required a source, in Watts, in the form:

S_{Hi} = h_{ij} × A_{S} × (T_{j} - T_{i})

where:

- h
_{ij}= the heat transfer coefficient, W/m^{2}K - A
_{S}= total surface area of particulates, m^{2}and - the subscripts i and j refer o the current and other phase.

For the particulates of spherical shapes of diameter D_{p}, their surface area
is given by

A_{S} = 6 × R2 × Cell Volume / D_{p}

The heat transfer coefficient is then computed from the local Nusselt Number, Nu = h_{ij}
× D_{p}/k, thus

h_{ij} = k × Nu / D_{p}, where

k = heat conductivity of the carrier phase, W/mK.

Two options for computing Nu have been provided in PHOENICS.

- CINT(H
_{i}) = GRND7

This computes Nu from "fit" to the experimental data on heat transfer from spheres, and it is valid over the complete Re and Prandtl number range. - CINT(H
_{i}) = GRND8

This computes Nu from the so-called Ranz-Marshall correlation, which is valid for 0 <Re < 200.

Notes:

- In both cases, the particle diameter is set through CINH2C.
- However, if CFIPS=GRND7 or GRND8, and SIZE is STOREd, then this will be used instead.
- The command STORE (NUSS) will store the calculated Nusselt Number in a 3-D storage for plotting.

Two variants of CINT(H_{i}) =GRND7 or GRND8 deal with phase specific heats as
follows:

- For constant and equal phase specific heats the following settings should be made:
- If phase 1 is taken as continuous phase:

CINT(H1)=GRND7 or GRND8

CINT(H2)=1.E20

CFIPS=GRND7

CFIPA=minimum slip velocity

CFIPB=particle size - If phase 2 is taken as continuous phase:

CINT(H2)=GRND7 or GRND8

CINT(H1)=1.E20

CFIPS=GRND8

CFIPA=minimum slip velocity

CFIPB=particle size

- If phase 1 is taken as continuous phase:
- For constant and different equal phase specific heats:
- If phase 1 is taken as continuous phase only:

CINT(H1)=GRND7 or GRND8

CINT(H2)=1.E20

CFIPS=GRND7

CFIPA=minimum slip velocity

CFIPB=particle size

PHINT(H1)=C_{P1}/C_{P2}

- If phase 1 is taken as continuous phase only:

The interphase mass transfer rate, m_{ji} in kg/s, is controlled by the PIL
variable CMDOT, as follows:

m_{ji} = CMDOT × f_{IP}

where f_{IP} or FIP = the reference interphase transfer coefficient, kg/s.

Notes:

- The convention is that mass transfer from phase 2 to phase 1 is positive.
- If CMDOT is set to one of the GRND flags, there is no internal multiplication by FIP. The value to be provided is the mass transfer rate per cell.

Two built-in options provided for the calculation of m_{ji}:

- CMDOT = HEATBL
- The setting switches on the calculation of the mass transfer rate from the balance of heat through the interface between phases.
- The option is appropriate for liquid-vapour (eg, water-steam) mixtures in which evaporation or condensation occurs.

- CMDOT = GRND7
- The setting switches on the calculation of the mass transfer rate from the difference of current and saturation values of mixture fraction of coal/oil-derived gases.
- The option is appropriate for the simulation of the combustion of coal, wood or oil in pulverised state.

The virtual- (also called "added-") mass term in the momentum equations for dispersed two-phase flow represents the force required to accelerate the mass of the surrounding continuous phase, in the immediate vicinity of a dispersed-phase fragment, such as a bubble or droplet, when the relative velocity of the phases changes.

It is significant only if the continuous-phase density is of the same order of magnitude as ,or much greater than, that of the dispersed phase, as for water droplets in oil or for steam bubbles in water.

Otherwise, as for sand particles in air or for water droplets in steam (at pressure well below the critical), the virtual-mass effect may be neglected.

*Formulation*

PHOENICS employs the Drew and Lahey (1987) formulation whereby the following source terms are introduced on the right-hand sides of the two momentum equations:

S_{c} = r_{C} × C_{vm} × r_{d}
× a_{vm}

S_{d} = - S_{c}

wherein:

r_{c} = the local density of the continuous phase,
kg/m^{3};

S_{d} and S_{c} = the volumetric force vectors for continuous, c, and
dispersed, d, respectively, N/m^{3};

C_{vm} = the virtual mass coefficient;

r_{d} = the volume fraction of phase d and

a_{vm} = the virtual-mass acceleration vector, m/s^{2};

*Activation*

The virtual-mass sources for the momentum equations are activated by assigning a non-zero value to the PIL variable CVM. The value ascribed to CVM determines how the virtual-mass coefficient is to be calculated, as follows:

- CVM = 0.0 (the default) cuts out the virtual-mass terms entirely.
- CVM = positive constant, K say, sets the virtual-mass coefficient to C
_{vm}=K. - CVM = GRND1 selects: C
_{vm}= CVMA*r_{c},

where CVMA is a positive constant (=0.5 by default). - CVM = GRND2 selects: C
_{vm}= CVMA × ( 1 - 2.78 × min(0.2, r_{d}) )

, where CVMA is a positive constant (=0.5 by default). - CVM = GRND permits the user to supply his own formula for Cvm in Group 10 Section 5 of GROUND.

*Notes*

- The virtual-mass coding presumes phase 1 to be the continuous phase and phase 2 the dispersed phase.
- However, if CVM is set to a negative value (other than GRND), the reverse is presumed, so that phase 2 is taken as the continuous phase and phase 1 as dispersed.
- When STORE(VMSU,VMSV,VMSW) appears in the Q1 file, the virtual-mass forces for each cell of the continuous phase are placed in the 3D stores VMSU, etc, and may be printed in the RESULT file or viewed via PHOTON and AUTOPLOT in the usual way
- The coefficient C
_{vm}having been computed in the user-accessible subroutine GXVMCF. - This subroutine is to be found in the file GX2PHS.FOR, which is located in the directory d_earth/d_opt/d_twophs.
- The more details can be found in "Virtual Mass" Entry of Encyclopaedia.

*Interphase friction, heat transfer and mass transfer rate*

The recommended method of introducing non-standard correlations for dilute ( R2<<R1 ) suspensions is:

- Set CFIPS = GRND.
- In GROUND, compute:
FIP = 0.75×C

_{d}×RHO1×R2×Volume× V_{slip}/ D_{p}Where:

C_{d}= dimensionless drag coefficient

RHO1 = density of continuous phase, kg/m^{3}

R2 = volume fraction of dispersed phase

Volume = unblocked volume of the cell, m^{3}

V_{slip}= relative velocity of the phases, m/s

D_{p}= particle diameter, m - For the momentum equations, leave CINT and PHINT at the default.
- If the transfer rate is driven by bulk phase values, leave both PHINTs at default.
- For the enthalpy equations leave PHINT(H2) at the default and set:
CINT(H1) = GRND and CINT(H2)=GRND

PHINT(H1)=C_{P1}/C_{P2} - In GROUND compute:
CINT(H1)=CINT(H2)=C

_{i}with

C_{i}= (k·Nu/D_{p})·(6·R2/D_{p})·Volume·(C_{P1}+C_{P2})/(2·C_{P1}·C_{P1})

or

CINT(H2) = 1.e20 and

CINT(H1) = (k·Nu/D_{p})·(6·R2/D_{p})·Volume/(2·C_{P1})

Where:

k = heat conductivity of continuous phase, W/(mK)

Nu = dimensionless local particle Nusselt number

C_{P1}= continuous phase specific heat, J/(kgK)

C_{P2}= dispersed phase specific heat, J/(kgK)

- Set CMDOT = GRND.
- In GROUND, calculate the interphase mass transfer rate per cell, in kg/s, from an appropriate expression.

Note:

- Under these circumstances, there is NO internal multiplication by f
_{IP}, and the values set are used directly.

*Virtual-mass coefficient*

If CVM is set equal to GRND in Q1 the user can supply his/her own formula for C_{vm}.

It should be done in Group 10 Section 5 of GROUND.

For example, the FORTRAN lines

CALL SUB2(L0CVM,L0F(LD12),L0R2,L0F(R2)) DO 1052 I-1,NXNY F(L0CVM+I)=CVMA*(1.-2.78*AMIN1(0.2,F(L0R2+I))) 1052 CONTINUE

provide the instructions to compute C_{vm} identically to what is provided in
option CVM=GRND2.

- In many practical cases, the particle size will vary throughout the domain, as the result of combustion or evaporation/condensation.

- This can be dealt with by use of the SHADOW technique. This uses a third phase (the SHADOW phase), which behaves like the disperse phase (usually phase 2), but without interphase mass transfer.
- On the picture, the solid particles are phase 2, and the open ones the shadow phase.
- Then changes in particle size can be calculated from local volume fraction ratios:
D

_{p}/ D_{p,in}= ( R2 / RS )^{1/3} - Solution of the SHADOW phase is activated by the PIL command: SOLVE(RS)
- For the built-in CFIPS options, f
_{IP}will be automatically divided by ( D_{p}/ D_{p,in})^{2}, if SOLVE(RS) and LSG4=T appear in Q1.

- Relevant "help-file" entries are:
- Relevant "Encyclopedia" entries include:

Examples of two-phase flow can be found in the following sections of the Library.

- Two-Phase Flow: IPSA cases W001-W005
- Advanced Multi-Phase Flow: cases P200-P222

On the SHADOW volume-fraction technique:

On interphase drag laws and heat/mass transfer coefficients:

On turbulence modulation by particles:

Chen & Wood, AIChE Journal, Vol 32, No. 1, pp 163-166, (1986).

About multiphase flow in general:

On Basic Equations of Two-Phase Flow:

- D B Spalding, 'Mathematical Modelling of Fluid Mechanics, Heat-transfer and Chemical Reaction Processes: A Lecture Course', Imperial College CFDU Report HTS/80/1.
- M. Ishii, 'Thermo-fluid dynamic theory of two-phase flow', Eyrolles, Book Publication, (1975).

On IPSA:

- D B Spalding, 'Numerical Computation of Multiphase Flow and Heat-transfer', Contribution to 'Recent Advances in Numerical Methods in Fluids', pp 139-167 Eds. C Taylor & K Morgan, Pineridge Press, Swansea, 1980
- D B Spalding, 'IPSA 1981: New Developments and Computed Results', Imperial College CFDU Report HTS/81/2

svz/331/0201

jcl/351/0903