Chapter 3 Numerical Method for the VOF Equation
3.2 The approximation for the face value of α
The distribution of volume fraction is governed by the equation (2.13):
=0
⋅
∇
∂ +
∂ V
t
αr
α (2.13)
in which the first term on the left-hand-side is an unsteady term and the second term on the left-hand-side is the convection term. The computation of convective fluxes is a complex topic and there are many former researchers devoted to this issue. In the following contents, a simple and efficient treatment will be deduced.
The face value αf can be estimated as a function of neighbor cells. Figure 3.1 shows a control volume and its neighbor cells which determined by flow direction on the control surface. The subscripts U, D and A denote upwind cell, donor cell and accepter cell.
Figure 3.1 the relationship of a control volume and its neighbor cells
For the unstructured grid calculations, only two neighboring cell values can be obtained
directly for a consider face. To discretize the VOF equation, as will be seen in section 3.5, the value of α at each cell face need to be estimated by some approximations and then the value of upwind cell αU will be obtained. The most convenient way to estimate the face value
αf is either the upwind difference scheme or the central difference scheme. These two schemes can be combined as:
( )
The second term of the right-hand-side stands for an anti-diffusion correction to the upwind differencing. This approximation becomes central differencing for γ =1, which may cause oscillations in the region of large gradients. To ensure boundedness of solution γ must be limited. The families of schemes based on the Total Variation Diminishing (TVD) [16] flux limiters are examples of this approach. Furthermore, these schemes are implemented in the context of the normalized variables formulation (NVF) which is proposed by Leonard [16] for the development of normalized variables diagram (NVD) schemes originally. Sweby [22]
presented the limiter γ as a function of the ratio of two consecutive gradients r :
D
The volume fraction value can be normalized as:
U
U
Introducing normalized variables into equation (3.1) leads to
2 ) 3.3 Linear difference schemes
Figure 3.2 shows the Normalized Variable Diagram (NVD) for a number of linear difference schemes. In the Figure the normalized face value is plotted as a linear function of the normalized donor cell value.
Figure 3.2 Normalized Variable Diagram (NVD).
UDS (Upwind Difference Scheme) is unconditionally stable and always produces a bounded solution, but has high numerical diffusion and loses accuracy. DDS (Downwind Difference Scheme) is a first-order scheme that introduces negative numerical diffusion. CDS (Central Difference Scheme), LUS (Linear Upwind Scheme), Fromm’s scheme and QUICK (Quadratic Upwind Interpolation for Convective Kinematics) scheme are second-order
accuracy schemes. The shortcomings of high order schemes are the causes of solution oscillations in sharp gradient regions. In the above diagram, the line of high order schemes would pass the point (0.5, 0.75). As given in equation (3.6), the different scheme at face α~ f can be written as the combination of the UDS and an anti-diffusive term. Different schemes can be represented by specifying a specific flux limiter γ which is a function of the gradient ratio r . A list of the linear difference schemes is given in the following table. Their relationships are also plotted in Figure 3.3.
NVD Flux limiter function
Table 3.1 The normalized variable (NV) and flux limiter function of linear schemes.
Similarly, linear difference schemes can be plotted in TVD diagram (see Figure 3.3).
Figure 3.3 Relationship between γ and rfor linear difference schemes.
3.4 Non-linear difference schemes
There are many high-order linear approaches introduced in section 3.1. However these the high-order schemes do not satisfied the Convective Boundedness Criterion (CBC) which was proposed by Gaskell and Lau [23]. The boundedness criterion of CBC can be shown in the NVD (see the shadow area of Figure (3.2)) :
To satisfy the CBC, a high-order scheme has to be a non-linear framework. Hence, high resolution schemes (HRS) [20] are developed for solving this problem which have the following properties: (i) solutions are free from oscillations or wiggles, (ii) high accuracy is obtained around shocks and discontinuities and (iii) CBC is satisfied. Non-linear approaches get a flux-limiter function to ensure a bounded value. The schemes of SMART and STOIC are classified as NVD category.
Sweby [22] has shown that the following constraint on the limiter guarantees that the scheme satisfies the TVD condition.
0 ( f( )r , f( )) 2 r r
γ γ
≤ ≤ (3.9)
Figure 3.4 shows TVD constraints in NVD and TVD diagram. The hatched region is known as the second-order regime. Therefore, high resolution schemes which pass through the
shadow area of Figure 3.2 satisfy CBC, but not fit in TVD constraints completely. It is apparent that the TVD schemes are more constrained than the NVD schemes. The MUSCL, SUPERBEE, OSHER and Van Leer schemes are developed based on the TVD theory.
(a) (b)
Figure 3.4 TVD Constraints in (a) NVD and (b) TVD diagram.
The following table lists the NVD equation and flux limiter equation in TVD of non-linear high resolution schemes.
NVD Flux limiter equation
SUPERBEE
Table 3.2 The NVD equation and flux limiter function of non-linear difference schemes.
Figure 3.5 to 3.11 show the NVD and TVD diagrams of a number of high resolution schemes.
Figure 3.5 SMART
Figure 3.6 MUSCL
Figure 3.7 SUPERBEE
Figure 3.8 STOIC
Figure 3.9 OSHER
Figure 3.10 BDS
Figure 3.11 Van Leer
Figure 3.12 CHARM
3.5 The normalized variables for unstructured grids
For the unstructured grids, the upwind value αU is required for variable normalization, but is not readily availiable. A estimation of αU can be obtained with a pseudo-upwind node U located at a distance of negative δvDA away from donor node D (see Figure 3.13). Thus,
DA D A
U α α δ
α = −2∇ ⋅ v (3.10)
The above approximation is necessary to be bounded as:
{ }
{
max ,0 ,1}
min U
U α
α = (3.11)
Figure 3.13. The prediction of upwind value for an arbitrary cell arrangement
3.6 Discretization of the indicator equation
Finite volume method is a popular method in CFD and will be applied in this study. This method uses the integral form of the conservation equations. The solution domain is divided into a finite number of control volumes, and the conservation equations are applied to each control volumes. The finite volume method is suitable for arbitrary geometries and can accommodate any type of grid. The finite volume approach is perhaps the simplest to understand and to program. All terms that need be approximated have physical meaning which is why it is popular with engineers.
The finite volume discretization of the equation (2.13) is performed by taking integration of the indicator equation over a control volume :
=0
and then use the Gauss’s divergence theory :
=0
The first term at the left-hand-side is an unsteady term and can be discretized by:
(
oD)
where Δ∀ is the volume of the cell and the superscript n and o denoting new-time level and old-time level.
The second term at left hand side of equation (3.14) is the convection term and can be approximated as
f
where the subscript f denoting that the properties on the surface of a control volume.
For the discretization, the first-order implicit scheme is computationally robust, stable and efficient, but suffers from numerical diffusion. On the other hand, the first-order explicit scheme is stable only for limited time steps. In the presenting work, we will adopt the second-order Crank-Nicolson scheme with large time steps.
In the equation (3.15), the value αf is approximated by
Substituing above equation and using the apprximation given in equation (3.1) in which the anti-diffusion term is treated as a fully explicity manner, the discretization of indicator equation can be presented as:
∑
+where F& is the volume flux at each face, the subscript P denote the primary node and nb V the neighbor node..
3.7 Blending strategy
To further enhance accuracy and reduce numerical diffusion and dispersion, composite schemes with blending strategy is introduced. There is a weighting factor based on the angle between interface and direction of motion. If the cell has started to be filled with fluid from the upwind side of the interface and the interface is parallel to the cell face, then only fluid present at the downwind cell should be convected through the cell face. In this case a compressive scheme should be used (see Figure 3.14(a)). However, if the interface is perpendicular to the cell face, a HR scheme would be appropriate (see Figure 3.14(b)).
(a)
(b)
Figure 3.14 (a) the interface is parallel to the flow direction and the compressive scheme will be used (b) the interface is perpendicular to the flow direction and high resolution schemes are appropriate.
The above-mentioned situations represent extreme cases in which the interface is either parallel or perpendicular to the control volume face. In general, the angle between the interface and control volume face is between these two cases. The normalized face value can be written as
( )
f f HR( ( )
f)
e compressiv f
f α f θ α f θ
α~ = ~ + ~ 1−
) ( )
( (3.21)
where °≤ ≤ °
The following will introduce the HRIC, CICSAM and STACS schemes.
HRIC
The high-resolution volume tracking scheme (HRIC) [23] is base on a blending of the bounded downwind (BD in table 3.2) and upwind differencing schemes. The blending factor
( )
ff θ is given by
( )
f cos( f )f θ = θ (3.22)
The normalized form of HRIC is:
( )
The Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) [25] scheme is also based on the blending strategy. The HYPER-C scheme and the ULTIMATE-QUICKEST (UQ) scheme are utilized when the cell face is perpendicular to the interface and the UQ employed when the face and the interface are parallel.
The HYPER-C scheme is
The ULTIMATE-QUICKEST scheme is
( ) ( )
Therefore, the CICSAM scheme can be written mathematically as
(CICSAM) f HYPER C
( )
f f UQ( ( )
f)
The Switching Technique for Advection and Capturing of Surfaces (STACS) [24] scheme selects downwind differencing as compressive scheme and modified STOIC as high resolution scheme.
The ~ (mod )
STOIC ified
f
α is obtained from
⎪⎪
The blending between the two schemes is preformed using
( )
f[
cos4( f)]
f θ = θ (3.29)
The normalized variables relationship for the STACS scheme is given by
(STACS) f DD
( )
f f ified STOIC( ( )
f)
The following Figures 3.16(a), 3.16(b) and 3.16(c) show the relationship between compressive schemes and high resolution schemes given by equations (3.22), (3.27) and (3.28) for the above composite schemes. In Chapter 5, these different blending factors will be compared in test cases.
Figure 3.16a Figure 3.16b
Figure 3.16c
Figure 3.16a, 3.16b and 3.16c show the relationship between two schemes for HRIC, CICSAM and STACS.
The x-axis stand for θf
and y-axis means the value of f
( )
θf.
Chapter 4 Numerical Method for the Velocity Field
4.1 Introduction
This chapter uses finite volume method to discretize momentum equation and pressure equation. All details of each term will be described. In this study, PISO algorithm is adopted to solve unsteady flow field.
4.2 Discretization
Take volume integral of the the momentum equation (2.5) to yield :
( )
∫∫∫( )
∫∫∫ ∫∫∫ ∫∫∫where φ designates each velocity component. The volume integral of the convection, diffusion and pressure terms can be transform into surface integral by divergence theory as
∫∫∫
In the following the unsteady equation is treated with the use of fully implicit scheme, which is unconditionally stable in the linear analysis and allows larger time steps for unsteady calculation.
Unsteady term
The volume integral of the unsteady term can be discretized as
(
Po)
The old time level term would be absorbed into the source term.Convection fluxes
The surface integral for the convection term can be approximated by
∑
where m&f is the mass flux through the consider face.
As mentioned in the chapter 3, the convective flux is calculated as:
( )
⎟⎟where the subscript D denotes the donor node and A the acceptor node.
The convetive fluxes at donor node and acceptor node can be calculated as:
⎪⎩
In the above the subsript P stands for the primary node, nb the neighbor node (see Figure 4.1) In the momentum calculation the Van Leer scheme is adopted, which is given by
) 1
In equation (4.6), the first term on the right-hand-side is treated implicit by while and the second term is moved into the source term.
Figure 4.1 Illustration of the primary cell P and the neighbor cell nb with a consider face f.
Diffusion fluxes
The surface integral of the diffusion term is approximated by
( )
is a vector pointing in the direction from primary node to the neighbor node (Figure 4.1) defined by :
Pnb
The diffusion flux can then be expressed as:
( ) (S d)
in which the second term at the right-hand-side is also treated explicitly and is moved into the
source term.
Pressure term
The surface integral is approximated as
∫∫ =∑ =∇ Δ∀ Gravitational acceleration term
The gravitational acceleration term can be approximated as
∀
Surface tension term
The surface tension term in the equation (3.2) can be discreted as:
( )
∇ Δ∀As mentioned in the section 2.4, a smooth κ exists in the transitional region where P
( )
∇α P ≠0.where αf is the face value obtained via the interpolation from the two neighboring node values.
The curvature value which defined by equation (2.16) can be calculated by :
∑ ( )
⎥⎥The face value
( )
∇α f is calculated by interpolating from the neighboring node values.Arrangement of the difference transport equation
The discretized momentum transport equation was approximated can be arranged as
∀
where the source term S includes old-time level unsteady term, the anti-diffusion term, the cross diffusion term, the gravitational term and the surface tension term.
4.3 Pressure-velocity coupling with PISO
A non-iterative method named PISO which stands for Pressure-Implicit with Splitting of Operators proposed by Issa [23] for unsteady incompressible flow is used in the present study.
This breakthrough of PISO is attributed to the few more correctors at each time step while SIMPLE has only one corrector proceeded at each iteration step. PISO involve one predictor step and two corrector steps. It is described in the following.
Predictor step
The predictor step is to solve the momentum equation using the prevailing pressure field.
( )
The above equation solves the velocity field but the conservation of volume flux has notbeen satisfied yet. For the purpose of satisfying the continuity equation (2.12), the velocity at face needs to be calculated. The above equation can be rewritten as:
o
Similarly, the relationship between velocity and pressure on the control surface can be deduced as:
o In the equation (4.19), the superscript “─” stands for the value interpolated from the primary cell P and the neighbor cell C with a weighting factor w . f
Substituting equation (4.21) into equation (4.20) to obtain:
o
is an average value of primary cell and neighbor cell.
⎥⎥
The volume flux at face can be calculate by :
d
First corrector step
After the predictor step a velocity field denoted * is obtained. In this corrector step both the velocity and pressure fields are adjusted such that the continuity equation is satisfied. The corrected velocity is denoted by Vv**
and the corrected pressure by P*. They satisfy the following momentum equation.
( )
Dividing the equation byA yields: p* Substituting equation(4.19) from this equation yields a velocity correction eqution.
P
The velocity correction at cell face can be written in a similar manner.
f
The volume flow rate corrector can be obtained from the velocity correction by
f
which can be further rearranged into the following form:
( )
Hence, the corrected velocity can be presented as:
( ) ( )
The discretization of the continuity equation leads to:
* 0
Substituting equation (4.33) into above equation, a pressure correction equation is obtained :
2
Second corrector step
In the second corrector step the velocity and pressure fields are further adjusted Vv*** Substituting equation (4.28) for this equation, a second velocity correction is obtained as :
f The corrected volume flux at a consider face is
*
Similarly, the volume flow rate corrector is
f
In the first correction step ∑∀
f f
*
* is obtained and equal to zero. Enforcing the velocity Vv***
to satisfy the continuity equation yields:
* 0
Substituting equation (4. 39) into (4.40), the second corrector step can be arranged as:
2
( )
Although more corrector steps are needed completed satisfy for the conservation law, two corrector steps are sufficient to have the accuracy of solution within temporal truncation error.
The source term of each algebraic pressure correction equation has been divided into two parts. The first part is accounting for the orthogonality of structured grids, whereas the second one is due to the non-orthogonality of unstructured grids.
4.4 Solution procedure for the two-fluid system
All equations have been discretized and the coupling velocity and pressure has been addressed. It is now to describe the solution procedure of the two-fluid system.
The solution procedure for the two-fluid system is:
1. Initialize all variables at initial time t 0
2. Solve the VOF equation for volume fraction α by using the old time level volumetric fluxes.
3. Update the coefficients of the momentum equations. Use the new α values to obtain an estimate for new viscosity and density.
4. Solve the momentum equation and continue with PISO algorithm.
5. If the final time step has not been reached, advance to the next time step and return to step 2.
Chapter 5 Results and Discussion
5.1 Introduction
A CFD methodology for the volume tracking methods was introduced in the previous chapter. In this chapter, the methodology will be tested against analytical solution and experimental data. The test cases selected are:
Uniform density flow calculation free from gravity and surface tension force. The velocity field in specified without solving momentum equation. First, linear schemes and non-linear schemes are tested. Furthermore, different composite schemes are tests and compare with well-know schemes which based on blending strategy. In the numerical results the combination of SUPERBEE as compressive scheme and MUSCL as high resolution scheme is selected and is adopted in the following cases.
Surface tension forces apply to a square bubble. In this case, gravity is neglected.
The simulation of droplet splashing demonstrates the merging of two independent regions
An air bubble in water is another air/water computational example.
Inviscid and viscous, low amplitude sloshing.
Rayleigh-Taylor instability. This case is density-driven and only influenced by gravity.
Column collapse with and without a small obstacle
A hydraulic bore.
5.2 The VOF schemes test in an oblique velocity field
Two different hollow shapes are convected in an oblique velocity field and Figure 5.1 presents the initial condition of these two shapes. The following describe the information about these two shapes as:
1. A hollow square (Figure 5.1(a)) aligned with the coordinate axes of an outer side length 0.8 and inner side length 0.4.
2. A hollow circle (Figure 5.1(b)) with an outer diameter of length 0.4 and inner diameter of 0.8. Because the volume fraction of each cell can not be obtained easily, there is a special
treatment for the initial setting. In the transition region, the physical surface intersects the cell faces in two points. An approximated surface can be obtained by connecting these two points with a straight line. This line and cell faces build an approximated area (see the hitched area in Figure 5.2). Thus, the portion of an approximated area of a consider cell is the volume fraction value.
The computational domain is a 4× square, subdivided into 4 100×100 grids. All shapes are initially at (0.8, 0.8) with their exact positions centered at (2.8, 1.8) after 1.0 seconds.
Three different time steps, Δt=0.00333, 0.00667 and 0.01, yield Courant number of 0.25, 0.5, 0.75 which will be denoted by low, medium and high, over the uniform grid systems.
Contours are display for
0 . 05 ≤ α ≤ 0 . 95
with 10 levels.For the purpose of comparison the solution errors, the solution error defined as [36]:
∑ the initial condition.
(a)
(b)
Figure 5.1 The flow field domain and initial condition of (a) hollow square, (b) hollow circle
Figure 5.2 Method for catching the value of volume fraction.
The comparison of different schemes
There are many linear and non-linear schemes are considered for solving indicator equation. Linear schemes are upwind difference scheme, central difference scheme, downwind difference scheme and cubic upwind scheme (CUS in Table 3.1). These linear schemes are revised and satisfy the convective boundedness criteria (CBC). In other words, the normalized volume fraction value at a consider face α~ tracks the upwind difference line f in NVD while the donor cell value α~D less than zero or greater than unity. Some famous non-linear schemes are chosen in the test and these schemes (see Figure 5.3 and 5.4) are Bounded Downwind scheme (BD), SMART, MUSCL, OSHER, Van Leer, CHARM and SUPERBEE. The NVD equation and flux limiter function of these non-linear schemes are detailed in Table 3.2.
Among these schemes, upwind difference scheme is a diffusive scheme in each case.
Downwind difference scheme produces highly numerical dispersion and compress shapes, especially at high Courant number (CN=0.75). In hollow circle, downwind difference scheme
Downwind difference scheme produces highly numerical dispersion and compress shapes, especially at high Courant number (CN=0.75). In hollow circle, downwind difference scheme