The conservation laws for mass, momentum, and energy are used to describe the physical behavior of the fluid flow. The general form of the conservation equation for a flow property φ in the control volume (C.V.) system shown in figure 2.1, is
C D S flux over the boundary due to convection, Vr
the fluid velocity, FD the flux over the boundary due to diffusion, ∀ the control volume and S the control surface. We can use the Gauss’s theorem to rewrite equation (2.1):
C D S
The above equation can be rewritten to the general conservative differential form when the control volume is contracted to a single point:
S
2.3 Mass and momentum conservation equation
The conservation equations for mass and momentum can be obtained by substituting φ=ρ by neglecting the diffusion term and the source terms, and φ=Vr
with the assumption of a laminar Newtonian working fluid under unsteady and incompressible conditions with body force and surface tension force. The mass and momentum conservation equations can be written as following:
=0
the velocity, P the pressure, μ the viscosity coefficient, g the gravitational acceleration and f the surface tension. σ
The density and viscosity of the effective fluid in the equation (2.4) and (2.5) can be calculated by the volume fraction, as
)
where the subscripts 1 and 2 denote the two fluids, α is the volume fraction. The density and viscosity of the different fluids are considered as variables through the full domain but constants in each kind of fluids. As mentioned above, all properties are piecewise continuous due to the volume fraction.
2.4 VOF equation
In the following, we will define different fluids by using the volume fraction in the volume-of-fluid (VOF) method on the Eulerian grid system (see figure 2.2). The value of volume fraction can be defined as
Volume
Therefore, the fluids through the entire computational domain can be divided into three
categories by the volume fraction as
The two-fluid system is propagated as the Lagrangian invariant and thus has a zero material derivative [19]:
=0 The above equations (2.4) to (2.10) can describe the fluid flow of the two-fluid system.
However, the form of the volume fraction, shown as (2.10), is not a conservative one and is not suitable for numerical solution. Because of the reason, it must to be reformulated [33].
The mass conservation equation (2.4) is a conservation form. It can be rewritten as 0
This non-conservation form of the mass conservation equation is much suitable for the two-fluid system with high density ratio, because the Vr
on the interface is defined as continuous. Figure 2.3 shows the densities at the inlet and outlet are not the same in the closed domain. The velocity Vr
of the fluid of entering and leaving the domain is the same, but the momentum Vρr of the fluid entering and leaving the domain is different. In addition, the fluids of this study are the assumption that they are incompressible. By substituting equation (2.6) into equation (2.11) the non-conservative equation becomes
(
1 2)
2The equation (2.10) can be rearranged into a conservation form with the incompressible
condition by recognizing that ∇⋅αVr =α⋅∇Vr+Vr⋅∇α
In the present study, the continuity equation (2.12), the momentum equation (2.5) and the VOF equation (2.13) together with the equation (2.6) and (2.7) will be employed to model the two-fluid flow.
2.5 Surface tension
As mentioned above, the surface tension will be modeled by the continuum surface force (CSF) model [37]. Surface tension creates a pressure jump which supplies the mean interface curvature with its necessary work on the interface. The surface tension coefficient σ exists for any pair of fluids and its magnitude is determined by the nature of the fluids. The value of σ is always positive for immiscible fluids and negative for miscible fluids [38]. The pressure jump is a function of the mean interface curvature, and it can be shown as [39]:
1 2 tension coefficient and κ is the mean interface curvature. For κ>0 fluid 1 lies on the concave side of the interface and for κ<0 fluid 2 lies on the concave side (figure 2.4). The gradient of α which is zero everywhere except at transient region, gives the normal vector, which always point from fluid 2 toward fluid 1 (figure 2.4):
α
∇
=
nr (2.15)
Thus, the mean interface curvature κ can be rewritten in terms of divergence of the unit normal vector as:
⎟⎟
By substituting equation (2.16) into equation (2.14), the surface tension term in the
momentum equation can be expressed as
2.6 Boundary conditions
Inlet: A velocity distribution is specified at the inlet.
Outlet: The outlet boundary condition uses the fixed pressure boundary condition. The boundary values are obtained from convective boundary condition [40]
=0
where φ represent the transported property and Vvc
is the convective velocity.
Rigid boundary (walls): A rigid boundary is generally defined as a non-slip boundary condition (u=0, v=0).
Chapter3 Numerical Method
3.1 Introduction
In the chapter 2, the mathematical model of the two-fluid flow has been described in detail. It is a necessary to choose a suitable discretization method, such as the finite difference (FD), finite volume (FV) or the finite element (FE) methods. These methods approximate the differential equations by a system of algebraic equations. Finite volume method is the method which uses integral form of conservation equations. The calculated domain can be divided into many several control volumes. In our computation, the VOF equation and the momentum will be discretized by using the finite volume method. The coupling between pressure and velocity will be treated by the PISO algorithm [41].
3.2 Discretization of the VOF equation
The finite volume method for the VOF equation of equation (2.13) is first integrated over a control volume, and then can be transformed by the Gauss divergence theory as:
0
The unsteady term can be discretized as:
(
Dn Do)
where Δ∀ is the volume of the cell, the superscripts n and o denote respectively the new and old time steps, D is the donor cell.
The second term can be discretized as:
( ) ( *)f f
where the subscript f denotes the properties on the surface of a control volume.
In the above equation, the value of α is obtained by using the second-order Crank-Nicolson *f scheme.
The volume fraction on the considering face can be determined by a function of neighbor cells and a flux limiter function γ
( )
r , which will be introduced in next chapter, shown as:( )
2A D
f D r α α
α =α +γ ⎛⎜⎝ − ⎞⎟⎠ (3.6)
Substituting the above equations into equation (3.1) yields
∑
+ denotes the neighbor node.3.3 Discretization of the momentum equation
The momentum equations can be expressed by( )
V( )
Qwhere Qφ is the source term of momentum equation, and φ represents the velocity components. Then, take a volume integral of the above equation and make use of Gauss theorem to yield:
( ) ( )
3.3.1 Unsteady term
The volume integral of the unsteady term can be discretized as
( )
3.3.2 Convection term
The surface integral for the convection term can be approximated by
( )
( )f f fC f fwhere m&f is the mass flux through the considering face.
The convective flux on the considering face can be estimated as a function of neighbor cells by a flux limiter function γ
( )
r , which will be introduced in the next Chapter.( )
⎟⎟where the subscript D denotes the donor node and A the acceptor node. In our calculation, the Van Leer scheme will be employed in the momentum equations.
The form of flux limiter function of Van Leer will be described in Chapter 4. The convetive fluxes at donor node and acceptor node can be calculated as:
⎪⎩
where the subsript P denotes the primary node and nb stands for the neighboring node (see Figure 3.1).
3.3.3 Diffusion term
The surface integral of the diffusion term is approximated by
( )
( o )f f fDis a vector pointing in the direction from the primary volume center to the
neighboring volume center (Figure 3.1). The length dr
was considered to be the factor affecting the diffusion dominancy and numerical stability. Hence, the over-relaxed approach for dr
was introduced
2
The diffusion flux can then be expressed as:
( ) (S d)
The volume integral of the diffusion term is approximated by ( )P
V
Q dφ ∀ ≈ qφΔ∀
∫
(3.21)The source terms in the momentum equation are pressure, gravitational acceleration and surface tension terms. In the following, each term will be introduced.
Pressure term
The surface integral of the pressure term is approximated as
o o
f f S f
Pd S =
∑
P S = ∇ Δ∀P∫
uv r (3.22)Gravitational acceleration term
The gravitational acceleration term can be approximated as
o
V
gd g
ρ ∀ = ρ Δ∀
∫
v v (3.23)Surface tension term
The volume integral of the surface term can be obtained from the approximation
( )
∇ = Δ∀∑
two neighboring nodes. Thus,( )
3.3.5 Arrangement of the difference transport equation
The discretized form of the momentum equation is approximated by the following form.
n n o
3.4 Pressure-velocity coupling of the PISO algorithm
The method of Pressure-Implicit with Splitting of Operators, which is proposed by Issa [41], is called PISO algorithm. In this study, the PISO algorithm will be used to deal with the unsteady problems. In the following, the procedure of PISO algorithm is addressed.
Predictor step
The predictor step is to solve the momentum equation using the prevailing pressure field.
( )
* * o
P P nb nb P
nb
A Vr =
∑
A Vr + S− ∇ Δ∀P (3.27)The above equation solves the velocity field but the mass conservation law has not been satisfied yet. Dividing the above equation by A yields P
*
P
In the equation (3.29), the superscript “─” stands for the value interpolated from the primary cell P and the neighbor cell C with a weighting factor w . f
First corrector step
The corrector steps are taking care of the mass conservation law of the flow field by updating the corresponding pressure. The new velocities and the corresponding new pressure assumed to be obtained from the first corrector step are denoted with superscript ** and *.
( )
** * *
P P nb nb P
nb
A Vr =
∑
A Vr + S− ∇ Δ∀P (3.31)Dividing the above equation by A and using the definition of (3.29) P
*
The first velocity correction equation (V ′vP
) is obtained by V′≡V**−V*, P′≡P* −Po, and
At the cell face, the velocity correction equation can be obtained in a similar manner.
f
Therefore, the volume flux correction equation is obtained as
f
The first volume flux correction equation can be presented as
( ) ( )
The volume flux at the face (∀&*f) of the above equation can be obtained by the following.
The relationship between velocity and the pressure at the face can be written as the form of (3.30) similarly.
o
*
Substituting equation (3.40) into equation (3.39) to obtain:
o
is an average value of primary cell and neighbor cell.
⎥⎥
The continuity equation is discretized as
The pressure correction equation is obtained by substituting (3.38) into (3.46).
where
Second corrector step
To enhance the SIMPLE procedure PISO performs a second corrector step. Similarly, the momentum equation is taken as
( )
where the new velocities and the corresponding new pressure are denoted with superscript*** and **.
Similarly, the second velocity corrector can be deduced as
f Again, the volume flow rate corrector is
f
The second corrector step is
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.
Solution procedure of PISO
Step 1. read the velocities and pressure of the flow field from the old time level.
Step 2. solve the momentum equation (3.27) to get V *
Step 3. compute P′ by solving the first pressure correction equation (3.47) to update velocities and pressure to get V** and P . *
Step 4. compute P ′′ by solving the second pressure correction equation (3.51) to further update velocities and pressure to get V*** and P . **
Step 5. if the required time step is achieved, then stop the calculation and output the data otherwise proceed to the next time step and repeat all over the way from step 1 to step 4.
3.5 Solution procedure
The VOF equation and momentum equation have been discretized. In this section, the solution procedure of two-fluid flow system will be described.
1. Initialize all variables at initial time t 0
2. Solve the VOF equation for volume fraction α by using the old time 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.
Chapter4 High Resolution Differencing Schemes
4.1 Introduction
As mentioned above, an effective scheme adopted to solve the two-phase flow must have some features, such as small diffusion, boundedness, and maintenance of sharp interface. The schemes which are introduced in the first chapter can not include above features at the same time. For example, the first-order upwind difference scheme is bounded but too diffusive. The problem which results from numerical diffusion is very important in the two-fluid system with the interface between two fluids. The strong numerical diffusion will smear the characteristic of the step function on the interface. In this chapter, the linear and non-linear scheme will be presented, and a composite scheme with switching function which can maintain the sharpness and shape of the interface will be proposed.
4.2 Convective flux of volume fraction
As mentioned above, the discretization of the VOF equation has been established. In the calculation process, the value on the considering face of a cell must be approximated. The approximation of the face value is necessary to ensure the accuracy and stable. The face value
αf can be estimation about a function of neighbor cells. Only two neighbor cells should be considered in the unstructured grid system. Figure 4.1 shows a control volume and its neighbor cells including the upwind and accepter. The subscript U, D and A denote upwind cell, donor cell and accepter cell.
In general, the methods, such as the upwind difference scheme (4.1) and the central difference scheme (4.2), are adopted to approximate this value.
f D
α =α (4.1)
1( )
2
f A D
α = α +α (4.2)
where αf is the value of volume fraction on the face, αD the value of volume fraction of the donor cell, and α the value of volume fraction of the accepter cell.
These two schemes can be associated by the variable of γ shown as:
( )
2
f D A D
α =α +γ α −α (4.3)
The term which has γ in the equation (4.3) is called the anti-diffusion correction to the upwind differencing. When γ=1, the equation (4.3) becomes central difference scheme. Then, this approximation will result in oscillations in the regions where the gradients are large.
Because of the above, the variable γ must be limited. The schemes with limitation will present high accuracy and resolution results which guarantee boundedness. The schemes with total variation diminishing (TVD) flux limiters were proposed in [23] in order to ensure the bounded solution. These schemes are implemented in the context of the normalized variables formulation (NVF) [24] for the development of normalized variables diagram (NVD) schemes originally. The limiter γ is defined as a function of the gradient r [26], shown as:
with the normalization, we can get following equations:
0
Introduce this normalized variable into equation (4.3), and it can be rewritten as following equation, and α~ is just the function of f α~D.
( )(1 )
4.3 CBC and TVD constraints
In order to ensure a bounded value, the high resolution schemes (HRs) [23] must satisfy the Convective Boundedness Criterion (CBC) or total variation diminishing (TVD) condition.
The high resolution schemes will prevent the oscillation or wiggles and get more accurate results around shocks and discontinuities in the two-fluid flow simulation. The Convective Boundedness Criterion (CBC) which was proposed by Gaskell and Lau [27] can be shown in the NVD (Fig 4.2):
Sweby [26] has proposed another constraint which makes the scheme satisfy the TVD condition and it can be shown as:
0 ( f( )r , f( )) 2 r r
γ γ
≤ ≤ (4.9)
The above constraints can be illustrated from Fig. 4.2 to Fig. 4.4. Fig. 4.2 and Fig 4.3 present the comparison with TVD constraint and the CBC constraint in the NVD. TVD constraint in the TVD diagram is showed in Fig. 4.4, and the hatched region is known as the second-order regime.
4.4 Linear and non-linear schemes
In this section, several high order schemes will be introduced. Generally speaking, these schemes can be divided into linear and non-linear schemes. First, linear schemes are built and can be explained by using of the combination of UDS and an anti-diffusion term or by the flux limiter function γ( )r . The normalized variable and the flux limiter function of the linear schemes can be found in table 4.1. Furthermore, these schemes will be plotted in the normalized variables diagram (NVD) and total variation diminishing (TVD) diagram (see Fig.
4.5 and 4.6).
UDS and DDS are both the first-order schemes. UDS (Upwind Difference Scheme) is an unconditional stable scheme and provides a stable solution at any time. However, this scheme always causes too much numerical diffusion which decreases the accuracy of the simulation.
DDS (Downwind Difference Scheme) results in less numerical diffusion result but unbounded.
The other schemes, such as CDS (Central Difference Scheme), LUS (Linear Upwind Scheme) and Fromm scheme, are two-order accuracy schemes. The schemes, including QUICK (Quadratic Upwind Interpolation for Convective Kinematics) and CUS (Cubic Upwind Scheme), have higher-order accuracy than above schemes. In the figure 4.5 and 4.6, the lines of second- and third- order schemes pass by the point (0.5, 0.75) in the NVD and the point (1, 1) in the TVD diagram. These linear schemes except the upwind difference scheme are not satisfied Convective Boundedness Criterion (CBC). High resolution schemes are developed by changing the high-order schemes into a non-linear one which satisfy the CBC. These non-linear schemes can be divided two categories. SMART and STOIC schemes are referred to NVD scheme, and MUSCL, SUPERBEE, OSHER and Van Leer schemes can be referred to TVD scheme. In present paper, two modified NVD schemes called M-MUSCL (Modified MUSCL) and M-BDS (Modified Bounded Downwind Scheme) will be developed to simulate the indicator function of the volume fraction. The normalized variable and the flux limiter function of these non-linear schemes can be summarized in table 4.1, and Fig. 4.7 to 4.16 show the NVD and TVD diagrams of these non-linear schemes.
4.5 Composite scheme with switching function
Solving the face value of the volume fraction by using only one high resolution scheme can not give consideration to both less numerical diffusion and non-deformed interface. The high resolution schemes preserve the shape of the interface but can not reduce the numerical diffusion. The compressive schemes are able to reach less numerical diffusion but let the interface deformed. Most composite schemes switch the compressive scheme and high
resolution scheme with a switching function about the slope of the interface in order to overcome the above problem. The method about composite scheme has been used in the former investigations. The key issue is how to switch schemes not when. The switching function has the basic on the angle between interface and direction of motion. It will decide which kind of the scheme to calculate the normalized face value of the volume fraction. If the interface is perpendicular to the cell face, a high resolution scheme would be appropriate. If the interface is parallel to the cell face, a compressive scheme would be appropriate. In general, the angle between the interface and control volume face is between these two cases.
The normalized face value can be written as:
( ) ( ( ) )
The well known composite schemes, such as HRIC of Muzaferija [33] and CICSAM of Ubbink [35] will be introduced. The composite of modified MUSCL and modified BDS will be also developed in this section.
CICSAM
In order to reduce the numerical diffusion and keep the sharpness of the interface, the compressive scheme called Hyper-C has been proposed by Leonard [42]. In general, the most compressive scheme is very suitable to the two-fluid flow with moving interface, but the Hyper-C may sometimes make the interface deformed or wrinkled. Therefore, Ubbink uses the ULTIMATE-QUICKEST (UQ) scheme to preserve the shape of the interface. The UQ (4.12) and Hyper-C (4.13) scheme (Fig. 4.17) can be shown as:
( ) ( )
The Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) scheme is developed by using the composite of the Hyper-C and UQ scheme with switching function
( )
fψ θ (Fig. 4.18) given as:
Then, the CICSAM scheme can be written as:
(CICSAM) f HYPER C
( )
f f UQ( ( )
f)
This composite scheme, like the above method, switches the upwind difference scheme and bounded downwind scheme with the switching function (Fig. 4.19) shown as
Then, the HRIC scheme can be shown as:
( )
Composite of modified MUSCL and modified BDS
⎥⎦
In present paper, the development of a new composite scheme switches the modified MUSCL and modified BDS scheme posed in the above chapter with the switching function (Fig. 4.20) given as:
( )
f[
cos4( f)]
f θ = θ (4.18)
This switching function has been proposed in [34], and this method can be formulated as:
This switching function has been proposed in [34], and this method can be formulated as: