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 criical), the virtual-mass effect may be neglected.

In the following discussion, rhoc stands for the local density of the continuous phase (called phase c), and rhod stands for that of the dispersed phase (called phase d).

Formulations of the virtual-mass term have been proposed by a number of workers, including Cook and Harlow [1984] and Drew and Lahey [1987]. PHOENICS employs the latter formulation whereby the following source terms are introduced on the right-hand sides of the two momentum equations:

Sc = rhoc * Cvm * rd * avm ; Sd = - Sc (1)

wherein: Sc and Sd are the volumetric force vectors for phases c

and d, respectively;

Cvm is the virtual mass coefficient;

rd is the volume fraction of phase d; and

avm is the virtual-mass acceleration vector, given by:

avm = D(Ud,i)/Dt - D(Uc,i)/Dt

= d(Ud,i-Uc,i)/dt + Ud,i d(Ud,i)/dxi - Uc,i d(Uc,i)/dxi (2)

Here D/Dt is the "substantial" derivative, and Ud,i and Uc,i are the velocity vectors of phases d and c, respectively.

The dimensionless coefficient Cvm describes the volume of displaced fluid that contributes to the effective mass of unit volume of the dispersed phase.

In general, this parameter is likely to be a function of the Reynolds number of the relative motion, and of other parameters such as fragment shape, but it is often taken as a constant. For spherical particles at high Reynolds numbers, Cvm=0.5 is appropriate.

PHOENICS provides not only for a constant value of Cvm, but also for values which are a function of rd, as follows:

Cvm = Cvma * rc ; and (3)

Cvm = Cvma * [1.-2.78 * min( 0.2,rd )] (4)

where Cvma is a constant, usually 0.5, and rc is the volume fraction of phase c. Equation (3) is based on the proposal of Watanabe et al [1990], while equation (4) is taken from Huang [1989] and Kowe et al [1988].

The virtual-mass source is computed within EARTH, 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.

Subroutine GXVMCF is called from Group 1 Section 1 and Group 10 Section 5 of GREX3.FOR.

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) and rc is the volume fraction of the continuous phase.

CVM = GRND2 selects: Cvm = CVMA*[ 1 - 2.78 * min( 0.2,rd ) ]

where CVMA is a positive constant (=0.5 by default) and rd is the volume fraction of the dispersed phase.

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.

CVM=GRND permits the user to supply his own formula for Cvm in Group 10 Section 5 of GROUND. For example:

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

computes Cvm from equation (4) above.

When STORE(VMSU,VMSV,VMSW) appears in the Q1 file, the virtual-mass forces for each cell of the continuous phase, as given by equations (1) and (2) and integrated over the control volume, are placed in the 3D stores VMSU, etc, and may be printed in the RESULT file or viewed via PHOTON and AUTOPLOT in the ususal way.

The virtual-mass momentum source term to be introduced is:

T = Cvm * integral{ rd*rhoc*[D(Ud,i)/dt - D(Uc,i)/dt] dVol } (5)

where the integral is over the cell volume. The term T, which must be subtracted from the dispersed-phase momentum equation and added to the continuous-phase momentum equation, is expressed approximately as:

T = Cvm*[ - (rhoc/rhod)*sum( Md,n*[Ud,n-Ud] )

+ (rd/rc)*sum( Mc,n*[Uc,n-Uc] ) ] (6)

where Md,n and Mc,n are the phase-specific coefficients of the finite-volume equations representing the effects of spatial and temporal convection. The subscript n stands for 'neighbour' in space and time.

The word 'approximately' is appropriate because the volume fractions and densities used in formulating the Md's and Mc's are neighbour-cell rather than in-cell values. This approximation is insignificant in comparison with the uncertainty regarding the proper value of the virtual-mass coefficient Cvm.

The finite-volume form of the indivicual-phase momentum equations can be expressed in correction form as:

Ud' = [sd + f*(Uc-Ud)-T]/apd (7)

Uc' = [sc + f*(Ud-Uc)+T]/apc (8)

where f is the interphase-friction coefficient, apd and apc are the central coefficients, sd and sc are all other terms in the corresponding momentum equation, and the ' denotes the correction to be applied.

In the interests of convergence, T is expressed as (9) T = T* + T'

where the * denotes the value based upon current values, and the correction T' is given by

T' = g*Ud' - G*Uc' (10)

where: g = Cvm*(rhoc/rhod)*sum(Md,n), and G = Cvm*(rd/rc)*sum(Mc,n).

If equations (9) and (10) are substituted in equations (7) and (8), the following linearised form of the momentum equations is obtained:

Ud' = [sd + f*(Uc-Ud) - T* + (f+G)*Uc']/(apd + f + g) (11)

Uc' = [sc + f*(Ud-Uc) + T* + (f+g)*Ud']/(apc + f + G) (12)

In PHOENICS, equations (11) and (12) are combined and then rearranged for solution by the Partial Elimination Algorithm (PEA), as follows. First, equations (11) and (12) are re-written as:

Ud' = [ Bd + (f+G)*Uc' ]/Ad (13)

Uc' = [ Bc + (f+g)*Ud' ]/Ac (14)

and after successive substitution one obtains:

Uc' = [Bc*Ad + (f+g)*Bd]/[Ad*Ac - (f+g)*(f+G)] (15)

and

Ud' = [Bd*Ac + (f+G)*Bc]/[Ad*Ac - (f+g)*(f+G)] (16)

as the final form of the momentum equations.

Examples of the use of the virtual-mass feature may be found in the advanced multi-phase section of the Input File Libraries.

T.L.Cook and F.H.Harlow, 'Virtual mass in multiphase flow', Int. J.Multiphase Flow, Vol.10, No.5, p691, (1984).

D.A.Drew and T.J.Lahey, 'The virtual mass and lift force on a sphere in rotating and straining inviscid flow',Int.J.Multiphase Flow, Vol.13, No.1, p113, (1987).

G.F.Hewitt, J.M.Delhaye and N.Zuber,' Test 2.4, Sedimentation', in Multiphase Science and Technology, Vol.3, Ed. Book publication, Hemisphere Pub. Corp., p474, (1987).

B.Huang, 'Modelisation numerique d'ecoulements diphasiques a bulles dans des reacteurs chimiques', PhD Thesis, L'Universite Claude Bernard - Lyon, (1989).

R.Kowe, J.C.R.Hunt, A.Hunt, B.Couet and L.J.S.Bradbury, 'The effects of bubbles on the volume fluxes and the pressure gradients in unsteady and non-uniform flow of liquids', Int. J.Multiphase Flow, Vol.14, No.5, p587, (1988).

S.A.Morsi and A.J.Alexander, 'An investigation of particle trajectories in two-phase flow systems', J.Fluid Mech., Vol.55, p193, (1972).

D.B.Spalding, 'The virtual-mass force in two-phase flow', CHAM Technical file note dated 7/5/97, (1997).

T.Watanabe, M.Hirano, F.Tanabe and H.Kamo, 'The effect of the virtual mass force term on the numerical efficiency of system calculations', Nuclear Engng. & Design, Vol.120, p181, (1990).

wbs