Free-surface flows involve the interaction of two or more distinctly different fluids, separated by sharply defined interfaces.
The position of the interface is not known a priori.
The mathematical model will need to:
Two special techniques for treating free surface flows are available in PHOENICS.
The Scalar-Equation Method deduces the interface position from the solution of the conservation equation for a scalar "fluid-marker" variable. This works well, but requires special techniques to combat numerical diffusion.
The Height-of-Liquid Method, when it is applicable, i.e. when the interface is not convoluted, needs no anti dispersion device, and is simple, effective and economical.
The two methods are similar, in that a "fluid-marker" is used to determine the physical properties. They differ in the method used to determine the distribution of the marker variable.
This lecture concentrates on the Scalar-Equation Method.
Other ways in which a free surface can be represented in PHOENICS include:
These methods suffer from various drawbacks, such as smearing of the interface, computational expense, or providing solutions on only one side of the interface.
The Scalar Equation (and Height of Liquid) methods are both essentially single phase treatments.
The two fluids separated by the interface are treated as both belonging to a single phase, i.e. there is only one value of each velocity component, temperature, concentration, etc., for each computational cell.
The relevant governing equations need to be solved, accounting for the discontinuities in the physical properties at fluid interfaces.
The location of the interface is deduced from some marker variable.
We will now look at the set of equations that need to be solved.
The single phase continuity equation can be written as:
dr/dt + d(rui )/ dxi = 0
This can also be written as:
D(lnr)/Dt + dui/dxi = 0
Assuming incompressible flow:
dui/dxi = 0
This implies that knowledge of the density is immaterial for the solution of the continuity equation if the flow is incompressible.
The continuity condition in terms of volumetric conservation is valid even when the density changes from point to point at the interface.
This is embodied in GALA (Gas And Liquid Algorithm), available in PHOENICS.
The conservation equations for other variables (i.e. momentum, energy, etc.) are solved in the conventional manner.
However, density and viscosity fields which are dependent on the local fluid type are required for the solution of these equations.
The remainder of this lecture concentrates on the Scalar Equation Method, which is a means of determining the property fields.
This method deduces the fluid interfaces from the solution of a conservation equation for a scalar "fluid-marker" variable.
The local physical properties such as density and viscosity are set to those of the appropriate fluid according to the value of the scalar marker variable.
The transient convection of a scalar variable is described by the following equation:
The scalar variable F is used as a marker according to which the fluid properties are set:
F = 0.0 - fluid 1
F = 1.0 - fluid 2
The governing hydrodynamic equations are solved in tandem with the transport equation for the scalar-marker variable.
The marker variable is used to update the fluid properties.
It is well known that the numerical discretisation of Equation (1) creates unphysical smearing of the discontinuity at the interface, known as numerical diffusion.
Several remedies exist to reduce this undesirable effect.
Special care in the specification of the properties:
At the interface, the properties are assumed to vary linearly from those of one fluid to those of the other.
The gradient of the property interface can be varied by limiting the range of PHI within which properties can vary.
Van Leer discretisation of the scalar-convection terms.
The "fully implicit upwind" (PHOENICS default) scheme considers the transport of a variable, F, across a cell face (e.g. the east face) according to the following formula:
Fe = FP for ue > 0 and
Fe = FE for ue < 0
where FP and FE are the values of F at the grid nodes at the end of the time step.
Figure 2: Higher-Order Differencing Template
The van Leer approach is an explicit finite-difference scheme which modifies the first order upwind formulation using the theory of characteristics:
Fe= FP + (dF/dx)P (dx - uedt )/2 for ue>0
Fe= FE - (dF/dx)E (dx + uedt )/2 for ue<0
where FP and FE are the values of F at the end of the time step.
The dF/dx is based on values at the start of the time step, which implies that the scheme is explicit in nature.
This gradient is specified such that the scheme reduces to the purely upwind form in the presence of local extrema.
This Total Variation Diminishing (TVD) approach offers the advantages of a higher-order scheme while avoiding under- or overshoots.
Due to the explicit formulation, the Courant criterion places a maximum limit on the time increment for the stability of the solution:
dt = min ( dx/u, dy/v, dz/w )
where this minimum is with respect to every cell in the grid, not just at the interface.
The SEM is applicable to:
Examples of applications:
The method cannot be used to produce steady-state solutions directly. However, steady-state conditions can be achieved asymptotically from a suitable transient simulation.
Surface tension effects at the interface are not included in the model. This is however a feasible extension to the model.
Numerical diffusion, despite treatments for its limitation, still imposes a restriction on minimum grid requirements.
Due to the explicit nature of the van Leer scheme, restrictions on the size of the time step increment can make long transient processes computationally expensive.
The Scalar Equation method has been implemented in the subroutine GXSURF which is available with PHOENICS.
It can be activated from PHOENICS-VR, by selecting Scalar Equation Method from the Free Surface Models menu of the Models panel of the Main menu.
Should the user wish to activate SEM 'by hand' from the Q1, the following commands are required:
In Group 7
STORE ( DEN1 , PRPS )
SOLVE ( SURN , VFOL )
SURN stores the scalar-marker variable, the value of which is unity in cells completely filled with the heavier fluid. VFOL stores the value of SURN at the old time step and provides information on the inflow volumetric fluxes.
Note that STORE(VFOL) will suffice if there is no inflow of fluid (i.e. for a sealed container). If inflows occur, VFOL should be SOLVED as it is required for the correct setting of the inflow boundary conditions.
In Group 8
GALA = T
Activates a solution procedure based on volume continuity.
TERMS(SURN,N,N,N,N,P,P); TERMS(VFOL,N,N,N,N,P,P)
TERMS ensures that the solution of SURN and VFOL is performed completely in GXSURF. This is required because the equation has no density in it.
RLOLIM is the value of SURN below which the cell will be regarded as being full of the lighter fluid. RUPLIM is the value of SURN above which the cell will be regarded as being full of the heavier fluid.
These two parameters are used to sharpen the fluid interface.
Typical values might be RLOLIM = 0.4 and RUPLIM = 0.6
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
Conditions are required for SURN, VFOL, DEN1 and PRPS to determine the initial position of the interface.
FIINIT, PATCH and INIT commands can be used to set the initial distributions of the above variables.
Typical settings to place a volume of water in surrounding air would be:
FIINIT(SURN)=0; FIINIT(VFOL)=0
FIINIT(PRPS)=0; FIINIT(DEN1) = 1.189
INIADD=F
PATCH (WATER, INIVAL, IXF,IXL, IYF,IYL,IZF,IZL, 1,1)
INIT (WATER, SURN, 0, 1)
INIT (WATER, VFOL, 0, 1)
INIT (WATER, PRPS, 0, 67)
INIT (WATER, DEN1, 0, 1000.5)
In the VR-Editor, the material for a Blockage object can be set to Light fluid or Heavy fluid.
In Group 13
PATCH ( NAME , CELL ,1 ,NX ,1 ,NY ,1 ,NZ ,1 ,LSTEP)
COVAL ( NAME , SURN , GRND , GRND )
This source reintroduces the transient term into the equation for SURN which was removed by the TERMS command.
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 )
COVAL ( NAME , SURN , FIXFLU , VEL * SURNin )
RHOM is the incoming fluid mixture density, and SURNin is the scalar value carried by the fluid.
The setting of an inflow value for VFOL is a mechanism which will satisfy the continuity equation, as in GALA the default source for inflows in the overall continuity equation is:
Sv=mass inflow/(cell density)
If the cell density is different from that of the incoming flow, the volumetric source will not be correct. However GALA will recognise the inlet value of VFOL and set the source in the overall continuity equation as:
Sv= mass flow × incoming VFOL =
(RHOM*Vel)
×(1./RHOM) = Velocity
Outflow Conditions
PATCH (NAME, AREA , IF, IL, JF, JL, KF, KL, TF, TL)
COVAL (NAME, P1, GRND1, Pext)
COVAL (NAME, SURN, FIXFLU, GRND1)
Fixed pressure conditions are generally applied. If the interface crosses the boundary, the above settings ensure that the continuity and scalar marker equations are treated correctly.
In Group 16
LITER(SURN)=1
Ensure solver does only one iteration for SURN.
In Group 18
VARMIN(SURN)=0; VARMAX(SURN)=1
Ensure SURN does not go outside physical limits.
In Group 19
SURF = T
This activates the special ground subroutine.
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.
An account has been given of the Scalar Equation Method, and its implementation in PHOENICS.
The method is provided in the subroutine GXSURF.
The method has been applied successfully to a variety of free-surface flows, especially those involving multiple and overturning interfaces such as mould filling.
Limitations on the application of the method have been indicated.
The PHOENICS Encyclopedia contains further details under the headings:
Free-Surface-Flow; Scalar-Equation 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
About the Scalar Equation method:
Jun L, Spalding DB 1988 "Numerical simulation of flows with moving interfaces." PHOENICS Journal of Computational Fluid Dynamics Vol 10, No 5/6, pp 625-637 1988. Published by Pergamon Press
Hamill IS, Jun L, Waterson N 1991 "A model for the simulation of three-dimensional mould-filling processes with complex geometries." Presented at the International Conference on Mathematical Modelling of Materials Processing, Bristol, 23-25 September 1991
Spalding DB, Jun L 1988 "Numerical simulation of flows with moving interfaces." Published in PCH PhysioChemical Hydrodynamics Vol 10, No. 5/6, pp 625-637 1988
