The lecture on Scalar Equation Method introduced free-surface flows, which involve the interaction of two or more distinctly different fluids, separated by sharply defined interfaces.
A mathematical model, the Scalar Equation Method (SEM), was described. The main limitation of the SEM is that it requires special techniques to combat numerical diffusion.
This lecture will introduce a method, namely the Height-of-Liquid Method, which, when it is applicable, needs no anti dispersion device, and is simple, effective and economical.
As for the SEM method, continuity is expressed from a volume conservation equation instead of mass conservation. This is activated by setting GALA=T in PHOENICS.
The density field must be known in advance, and at all times for the solution of the other conservation equations, i.e. momentum, energy, concentrations, etc. These are then solved in the conventional manner.
The remainder of this lecture will concentrate on the Height-of-Liquid (HOL) Method, which is a means of determining the density field.
The HOL method is applicable to:
The only restriction is that the flow considered should exhibit no "overturning" of the interface.
This means that there must exist one direction, designated the "up" direction, along which only ONE intersection of the interface exists.
Figure 1: Overturning Interface
In the HOL method, the location of the interface is determined from the solution of liquid-balance equations. To facilitate the understanding of the description which follows, several main variables of the HOL method are defined:
RHOL: liquid density;
RHOG: gas density
VFOL: volume fraction of liquid; CVOL: cell volume
TMOL: total mass of liquid in a vertical column of cells.
LMOL(IZ): lower mass of liquid, or TMOL below the interface cell. e.g. if IZ=4, the sum of liquid masses in cells 1, 2, and 3.
MOLO: total mass of liquid in a vertical column at the old time step. (i.e. the old equivalent of TMOL)
Mass balance definition
INOU: sum of inflows of liquid to the cell column minus the sum of outflows of liquid from the cell column. The inflows and outflows are computed from the volumetric flows obtained as convective fluxes, multiplied by the upwind values of VFOL*RHOL.
Definition of the Interface
The position of the interface is defined by calculating the volume fraction of liquid, VFOL in every cell, from:
VFOL = (TMOL-LMOL)/(RHOL*CVOL)
The interface is defined to lie in the cells with a VFOL in the range 0.0 <VFOL < 1.0
The unknowns to be determined are TMOL & LMOL.
TMOL: always known at the beginning of the calculation.
For all subsequent times, it is determined from:
TMOL = MOLO + INOU
LMOL(IZ): is in fact the mass of liquid in a column of cells below the current IZ, and is related to the cell volume CVOL by:
LMOL(IZ) = LMOL(IZ-1) + RHOL*CVOL(IZ-1)
with LMOL(1) = 0
LMOL is greater than TMOL when the interface lies below the currently-considered IZ, producing a negative VFOL, which is of course unrealistic, and which should be increased to zero.
Similarly, when TMOL is greater than LMOL, VFOL could be greater than one, which should be reduced to one.
Therefore:
VFOL = max (VFOL, 0.0)
VFOL = min (VFOL, 1.0)
NB: The solution of VFOL is performed entirely within the above equations.
The density field can be calculated as follows:
RHO = RHOG * (1-VFOL) + RHOL * VFOL
and the viscosity field is calculated as follows:
ENUL = ENULG * (1-VFOL) + ENULL * VFOL
In practice, VFOL is encoded into the PRPS variable at the start of each timestep. which is then used to obtain the cell-by-cell properties.
The HOL method has been implemented in the GX routine GXHOL, which is available with PHOENICS.
It can be activated from PHOENICS-VR, by selecting Height of Liquid from the Free Surface Models menu of the Models panel of the Main menu.
Should the user wish to activate HOL 'by hand' from the Q1, the following commands are required:
In Group 7
STORE ( DEN1, PRPS ) SOLVE ( VFOL )
VFOL stores the volume fraction of liquid in each cell.
If there are no inflows (i.e. for a sealed container) it is sufficient to STORE(VFOL) since the solution of VFOL is performed entirely within GXHOL.
If inflows occur, VFOL should be SOLVEd as inlet boundary conditions must be obtained.
In Group 8
GALA = T
Activates a volume-continuity-satisfying flow field.
TERMS(VFOL,N,N,N,N,P,P)
TERMS ensures that the solution of VFOL is performed in GXHOL.
In Group 9
The flow properties are calculated in GXPRPS, based on the local property marker, PRPS. PRPS is updated from the VFOL field at the start of every time step.
No specific settings are needed in Group 9 - the PRPS values for the heavy and light fluids are set in Group 19.
In Group 11
Initial conditions are required for VFOL and DEN1 for the initial position of the interface.
Typical settings to place a volume of water in surrounding air would be:
FIINIT(VFOL)=0; FIINIT(DEN1)=1.189
INIADD = F
PATCH (WATER, INIVAL, IXF,IXL, IYF,IYL, IZF,IZL,ITF,ITL)
INIT (WATER, VFOL, 0, 1)
INIT (WATER, DEN1, 0, 1000.5)
It is not necessary to explicitly set PRPS, as this is updated internally from the VFOL field.
In the VR-Editor, the material for a Blockage object can be set to Light fluid or Heavy fluid.
In Group 13
Inflow conditions
INLET ( NAME , AREA , IF,IL,JF,JL,KF,KL,TF,TL)
VALUE ( NAME , P1 , VEL * RHOM )
VALUE ( NAME , U1 , UIN )
VALUE ( NAME , V1 , VIN )
VALUE ( NAME , W1 , WIN )
VALUE ( NAME , VFOL, 1. / RHOM )
RHOM is the incoming fluid density, computed from the rule given in panel 11.
In GALA the default source for inflows in the continuity equation is:
Sv=mass inflow/(cell density)
If the cell density is different from that of the incoming flow, this is not correct. GALA will, however, recognise the inlet value of VFOL and set the source in the continuity equation as:
Sv= mass flow * incoming VFOL = (RHOM*Vel) * (1./RHOM) = Vel
Outflow conditions
PATCH (NAME, AREA, IF,IL, JF,JL, KF,KL, TF,TL)
COVAL (NAME, P1, GRND, Pext)
Fixed-pressure conditions are generally applied. If the interface crosses the boundary, the pressure Coefficient must be replaced with the local value of the density, which is done automatically by using GRND as the Coefficient.
In Group 16
RELAX(P2, LINRLX, factor)
Applies a linear relaxation factor to the rate of change of TMOL with sweep. For steady calculations, a very small number is common.
For transient calculations, no relaxation is necessary or possible, as the relaxation factor is set internally to 1.0 .
Note that the use of P2 is merely a device to pass a relaxation factor to a variable solved wholly within HOL. There is no need to SOLVE or STORE P2.
In Group 19
HOL = T (activates the special ground subroutine).
The first sweep at which TMOL is calculated is set using:
SPEDAT(SET,HOL,IFIRST,I,ifirst).
where ifirst is the number of the first adjustment sweep. The default value is 5. For transient cases this is always 1 and cannot be changed.
The frequency in terms of sweeps of the calculation of TMOL is set using:
SPEDAT(SET,HOL,IFREQU,I,ifrequ).
where ifrequ is the sweep frequency for adjustments. The default value is 5. For transient cases this is always 1 and cannot be changed.
IHOLA - This determines the "up" direction.
IHOLA = 1 : selects X = 2 : selects Y = 3 : selects Z
Negative values of IHOLA signify that the height of liquid is calculated with reference to the "high" end of the coordinate direction.
IPRPSA is the index selecting the heavier fluid in the PROPS file.
IPRPSB is the index selecting the lighter fluid in the PROPS file.
For example, if the heavier fluid is water, and the lighter fluid is air at constant density, IPRPSA should be set to 67, and IPRPSB to 0.0.
HOL cannot cope with "overturning" motions of fluid.
When NONORT is set to TRUE for highly non-orthogonal grids, the way that cell face mass fluxes are calculated may cause convergence problems for the HOL method.
The surface tension effects at the interface are not included in the present model. This could be implemented in the model.
HOL is not immediately amenable to the simulation of flows where the domain consists of inter-connected regions of open spaces and solid blocks.
Any blockages present must be at least partially immersed at the start of the calculation, or located on the edge of the domain. The density of the blocks must be the same as that of the liquid.
Conjugate heat transfer calculations are possible, but great care must be taken with the properties of the block bearing in mind that the density must be that of the fluid.
Hence, for transient cases, the specific heat must be multiplied by the density ratio to ensure the correct thermal capacity of the block.
The following cases cannot be performed with HOL.
Figure 3: Cases which cannot be treated with HOL
Advantage of the HOL compared to the SEM:
Advantage of the SEM compared to the HOL:
An account has been given of the HOL method, and its implementation in PHOENICS.
The method is provided in the subroutine GXHOL.
The method is very successful in cutting down numerical diffusion on the evolution of an interface.
Limitations on the application of the method have been indicated.
The PHOENICS Encyclopedia contains further details under the headings:
Free-Surface-Flow; Height-of-Liquid Method
Library examples can be found in the multi-Phase flow library, which can be accessed from the library option of the top menu, or from command mode by SEELIB(P).
Active demonstrations can be found in Active Demos, under
Extra multi phase features
Free surface flow
