Encyclopaedia Index

## TWO-PHASE FLOWS

### Contents:

1. General introduction
2. IPSA in PHOENICS
3. Interphase transport processes
4. Particle-Size Calculation
5. Encyclopaedia and Help Entries
6. Library Examples
7. References

### Practical Importance

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.

### Definition Of Phase

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.

### Methods Of Solution

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

• up to 3 velocity components for each phase, ui
• 1 volume fraction for each phase, ri

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 i at the point and instant in question;
• 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 ri, ui, vi, wi, Ti, Ci, 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.

### Phase continuity

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

d(riri)/dt + div( ririVi - Grigrad(ri) ) = �ji
Transient Convection Phase Diffusion Mass Source

Where:

• ri = phase volume fraction, m3/m3
• ri = phase density, kg/m3
• Vi = phase velocity vector, m/s
• Gri = phase diffusion coefficient, Ns/m2
• ji = net rate of mass entering phase i from phase j, kg/(m3s)

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:

r1 + r2 = 1

### Phase conservation equations

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

Transient Convection Within-phase and Phase Diffusion Sources

Where:

• ri = phase volume fraction, m3/m3
• ri = phase density, kg/m3
• Vi = phase velocity vector, m/s
• Gfi = within-phase diffusion coefficient, Ns/m2
• Gri = phase diffusion coefficient, Ns/m2
• Si = within-phase volumetric sources, kg�f/(m3s)
• Sip = interphase volumetric sources, kg�f/(m3s)

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.

### Implementation and Activation

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.

### Pairing of Variables

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

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.

### Boundary Conditions

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.

### Further Details Of Mass Flow Conditions

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

m' = Type � Cp1 � ( Vp1 - 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 � Cp1 � ( Vp1 - P1 ) for P1 < VPi {inflows}
m' = Type � R1 � Cp1 � ( Vp1 - P1 ) for P1 > VPi {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:

Gfi = RHOi � ( 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:

Gri = RHOi � ( ENUL / PRNDTL(Ri) + ENUT / PRT(Ri) )

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.

### Interphase Source Term

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

SIP=(ff,i+ <mji>) (fiint -fi)

where:

ff,i = interphase (diffusive) transfer coefficient, kg/s;
mji = net mass transfer rate between the phases, kg/s;
<> = maximum of 0.0 and the quantity enclosed;
fiint = value of f at the interface between the phases and
fi is a bulk phase value of conserved variable f.

Notes:

• The unit of SIP is (kg/s) � (unit of f),i.e.
• if f is velocity, m/s, the units of SIP are Newtons and
• if f is enthalpy, J/kg, the units of SIP are Watts.
• SIP, obviously, appears in the equations for each of paired variables, f1 and f2 ,i.e

SIP,1=(ff,1+ <m21>) (f1int -f1)
SIP,2=(ff,2+ <m12>) (f2int -f2)

wherein m21 = - m12;

• 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 Interface Value: PHINT(f)

The concept of the interface value fiint is the value of fi 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 fi 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.

fiint = fj

Then, interphase source becomes

SIP=(ff,i+ <mji>) (fj -fi)

where:
fi is current phase bulk value of f and
fj 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(f1) = default and PHINT(f2) = default.

Known interface values

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

fiint = Value1 and
fjint = Value2

interphase source becomes

SIP=(ff,i+ <mji>) (Valuei-fi)

• 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(f1) = Value1 and
• PHINT(f2) = Value2

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(fi)=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 = fiint - fjint

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

fiint = fj + diff
fjint = fi - diff

and resulting interphase source becomes

SIP=(ff,i+ <mji>) (fj - fi + 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(f1) = default and
• PHINT(f2) = diff

Known ratio of interfacial values.

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

Value1 = fiint/ fjint

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

fiint = Value1fj
fjint = fi/Value1

and resulting interphase sources become where:

C1 and C2 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(f1) = Value1 and
• PHINT(f2) = default

### Interphase Transfer Coefficient - FIP

ff,i, kg/s, featured in SIP is actually the interphase flux per cell per unit conserved variable difference.

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

fIP or FIP, kg/s=Ns/m, is the force per cell per unit velocity difference.

It appears as multiplier, in the built-in formulations of all specific ff,i and mji. This multiplier is introduced because of the analogy between friction, heat and mass transfer.

### Reference Interphase Transport Coefficient - CFIPS

The coefficient fIP is controlled by the PIL variable CFIPS.

If CFIPS is any constant other than a GRND flag, fIP is set as follows:

fIP = CFIPS�RHO1�R1� R2*�Vol, for CFIPS >0.0
fIP = |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.

### Interphase Friction - Built-in Options

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

SIP =fIP � ( velj - veli)

where SIP has units of Newton and fIP 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 fIP 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:
• fIP =0.75 �Cd�RHO1�R2�R1* �Vol�<Vslip,CFIPA> / Dp

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 , Vslip, 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:

• 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 fIP 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 fIP 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.:

fIP = 0.75�RHO2�R1�R2*�Vol� <Vslip,CFIPA> / Dp

Notes:

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

### Interphase Heat Transfer - CINT

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:

SH1 = h12 � AS � (T1int - T1)
SH2 = h21 � AS � (T2int - T2)

where:

h12 = bulk1-to-interface heat transfer coefficient, W/m2K
h21 = bulk2-to-interface heat transfer coefficient, W/m2K
AS = total interface area, m2

The interface temperatures, T1int and T2int, can be eliminated via overall heat transfer coefficient and bulk-to-bulk temperature difference as follows:

SHi = 1 � AS � (Tj - Ti)/(1/C1 - 1/C2)

where: C1 = hij and C2 = hji

This formulation is generalised so that the interphase transfer coefficient, fIP, is taken to be harmonic average of the bulk-to-interface coefficients, C1 and C2.

Thus,

fIP = 2 / (1/C1 - 1/C2)

Note: if C1 = C2 = C, then fIP = C

Built-in provisions have been made to calculate C1 and C2 via their ratios to the reference transfer coefficients, as follows

C1 = CINT(f1) � FIP
C2 = CINT(f2) � 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 Ci directly.
• Two options are provided internally for CINT(H1) and CINT(H2), see Encyclopaedia entry CINT for details.

### Interphase Heat Transfer - Built-in Options

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

SHi = hij � AS � (Tj - Ti)

where:

• hij = the heat transfer coefficient, W/m2K
• AS = total surface area of particulates, m2 and
• the subscripts i and j refer o the current and other phase.

For the particulates of spherical shapes of diameter Dp, their surface area is given by

AS = 6 � R2 � Cell Volume / Dp

The heat transfer coefficient is then computed from the local Nusselt Number, Nu = hij � Dp/k, thus

hij = k � Nu / Dp, where

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

Two options for computing Nu have been provided in PHOENICS.

• CINT(Hi) = 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(Hi) = 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(Hi) =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
• 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)=CP1/CP2

### Interphase Mass Transfer - CMDOT

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

mji = CMDOT � fIP

where fIP 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.

### Interphase Mass Transfer - Built-in Options

Two built-in options provided for the calculation of mji:

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

### "Virtual mass" effects

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:

Sc = rC � Cvm � rd � avm
Sd = - Sc

wherein:

rc = the local density of the continuous phase, kg/m3;
Sd and Sc = the volumetric force vectors for continuous, c, and dispersed, d, respectively, N/m3;
Cvm = the virtual mass coefficient;
rd = the volume fraction of phase d and
avm = the virtual-mass acceleration vector, m/s2;

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 Cvm =K.
• CVM = GRND1 selects: Cvm = CVMA*rc,
where CVMA is a positive constant (=0.5 by default).
• CVM = GRND2 selects: Cvm = CVMA � ( 1 - 2.78 � min(0.2, rd) )
, 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 Cvm 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.

### Introduction of Non-Standard Correlations

Interphase friction, heat transfer and mass transfer rate

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

1. Set CFIPS = GRND.
2. In GROUND, compute:

FIP = 0.75�Cd�RHO1�R2�Volume� Vslip / Dp

Where:
Cd = dimensionless drag coefficient
RHO1 = density of continuous phase, kg/m3
R2 = volume fraction of dispersed phase
Volume = unblocked volume of the cell, m3
Vslip = relative velocity of the phases, m/s
Dp = particle diameter, m

3. For the momentum equations, leave CINT and PHINT at the default.
4. If the transfer rate is driven by bulk phase values, leave both PHINTs at default.
5. For the enthalpy equations leave PHINT(H2) at the default and set:

CINT(H1) = GRND and CINT(H2)=GRND
PHINT(H1)=CP1/CP2

6. In GROUND compute:

CINT(H1)=CINT(H2)=Ci with
Ci = (k�Nu/Dp)�(6�R2/Dp)�Volume�(CP1+CP2)/(2�CP1�CP1)
or
CINT(H2) = 1.e20 and
CINT(H1) = (k�Nu/Dp)�(6�R2/Dp)�Volume/(2�CP1)
Where:
k = heat conductivity of continuous phase, W/(mK)
Nu = dimensionless local particle Nusselt number
CP1 = continuous phase specific heat, J/(kgK)
CP2 = dispersed phase specific heat, J/(kgK)

7. Set CMDOT = GRND.
8. 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 fIP, 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 Cvm.
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 Cvm identically to what is provided in option CVM=GRND2.

### Particle Size Calculation

• 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:

Dp / Dp,in = ( R2 / RS )1/3

• Solution of the SHADOW phase is activated by the PIL command: SOLVE(RS)
• For the built-in CFIPS options, fIP will be automatically divided by ( Dp / Dp,in )2 , if SOLVE(RS) and LSG4=T appear in Q1.

### Encyclopedia And Help File Entries

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

### Library Examples

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

### References

• D B Spalding, 'The SHADOW method of particle size calculation in two phase flow combustion', 19th Symposium on Combustion, pp941-951. The Combustion Institute, Pittsburgh, 1982.
• On interphase drag laws and heat/mass transfer coefficients:

• R Clift, J R Grace & M E Weber, 'Bubbles, Drops and Particles' Academic Press, new York, 1978.
• J. T. Kuo and G. B. Wallis, 'Flow of bubbles through nozzles', Int. J. Multiphase Flow, Vol. 14, No. 5, p547, (1988).
• M. Ishii and N. Zuber, 'Drag coefficient and relative velocity in bubbly, droplet or particulate flows', AIChE Journal, Vol. 25, No. 5, p843, (1979).
• On turbulence modulation by particles:

• Mostafa & Mongia, Int J Heat Mass Transfer, 31,10, pp 2063-2075, (1988)
Chen & Wood, AIChE Journal, Vol 32, No. 1, pp 163-166, (1986).
• About multiphase flow in general:

• G Hetsroni (Ed.), 'Handbook of Multiphase Systems', Hemisphere, Washington, (1982).
• P. B. Whalley, 'Boiling, condensation and gas-liquid flow', Clarendon Press, Oxford, (1990).
• 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