A Novel Algorithm of Advection Procedure in Volume of Fluid Method to Model Free Surface Flows

Abstract

In this study, the developed procedure of advection in volume of fluid (VOF) method is presented for free surface modeling. The fluid is assumed to be incompressible and viscous and therefore, Navier-Stokes and continuity are considered as governing equations. Applying Youngs’ algorithm in staggered grids, it is assumed that fluid particles in the cell have the same velocity of the cell faces. Therefore, fluxes to neighboring cells are estimated based on cell face velocities. However, these particles can show different velocities between two adjacent cell faces. In developed model, the velocity in mass center of fluid cell is evaluated to calculate fluxes from cell faces. The performance of the model is evaluated using some alternative schemes such as translation, rotation, shear test, and dam break test. These tests showed that the developed procedure improves the results when using coarse grids. Therefore, the Modified Youngs-VOF (MYV) method is suggested as a new VOF algorithm which models the free surface problems more accurately.

1. Introduction

In the numerical computations of free surface flows such as water waves and splashing droplets, accurate representation of the interface is very important. The volume of fluid (VOF) method is a convenient and powerful tool for modeling the free surface flows, where the fluid location is determined using related function. In this method, the VOF function is averaged over each computational cell. It is set as unity, zero, and a value between unity and zero for full fluid, empty, and free surface cells respectively. Using this function, the VOF method is capable of modeling flows with complex free surface geometries such as rising bubbles [1], the merging and fragmentation of the drop [2]. In addition, in comparison with other methods, yet it is remarkably economical in computational point of view. It is due to requiring only one array for storing the VOF function and a simple algorithm to advect the function during each computational time step. Several volume advection techniques have been developed with the aim of maintaining sharp interface. The more famous ones are the simplified line interface calculation (SLIC) method of Noh and Woodward [3] such as SOLA-VOF [4] or its corrected form as in NASA-VOF2D [5], the VOF method of Hirt and Nichols [6], and the method of Youngs [7, 8]. VOF advection algorithm can be classified based on free surfaces reconstructing technique in each cell and the method of computing boundary flux integration. Some VOF methods represent free surface interfaces as a line parallel to one of the grid coordinates which are referred as piecewise constant scheme. Some of them are the methods used by Nichols et al. [4], Hirt and Nichols [6], Torrey et al. [5], and Duff [9] where free surface interfaces are constructed in a stair-shaped profile. The alternative methods are known as piecewise linear schemes. They are developed by Rider and Kothe [10], Harvie and Fletcher [11, 12], Geuyffier et al. [13], and Scardovelli and Zaleski [14, 15]. In these methods, oriented free surface interface is in a direction perpendicular to the locally evaluated VOF gradient. These schemes are complex but more accurate than their piecewise constant counterparts associated with more computational costs. In this paper, a new advection model as modified Youngs’ VOF (MYV) method is presented.

2. Governing Equations

The fluid is considered to be Newtonian and incompressible. Therefore, 2D continuity and Navier- Stokes equations (NSEs) are used as 𝜕𝑉∇⋅𝑉=0𝜕𝑡+(𝑉⋅∇)𝑉=−∇𝜌𝑝+∇⋅𝜈∇𝑉+∇𝑉^𝑇+𝐵,(2.1) where 𝑡 is time, 𝑉 velocity vector, 𝑝 hydrodynamic pressure, 𝜈 kinematic fluid viscosity, 𝜌 fluid density, ∇ gradient operator, and 𝐵 is gravity based body force. In turbulent flow, the effect of turbulence can be considered using eddy viscosity models [16]. Researchers have used different models such as 𝑘-𝜀 and 𝑘-𝑤 [17–20] to model turbulent flow. In this paper, the standard two equation 𝑘-𝜀 model is used, where the first equation involves turbulence kinetic energy (𝑘) and represents the velocity scale. The second one takes into account turbulent dissipation rate (𝜀) and represents the length scale. The two-equation 𝑘-𝜀 turbulence model accounts for the effect of turbulence as follows: 𝜕𝑘𝜕𝑡+𝜕𝜕𝑢𝑥𝑘+𝜕𝜕𝑣𝑦𝑘=𝜕𝜕𝑥𝜈^𝑘𝜕𝑘𝜕𝑥+𝜕𝜕𝑦𝜈^𝑘𝜕𝑘𝜕𝑦−𝐺^𝑆𝜕𝜀−𝜀,𝜕𝑡+𝜕𝜕𝑢𝑥𝜀+𝜕𝜕𝑣𝑦𝜀=𝜕𝜕𝑥𝜈^𝜀𝜕𝜀𝜕𝑥+𝜕𝜕𝑦𝜈^𝜀𝜕𝜀𝜕𝑦−𝐶¹𝐺_𝑆−𝐶₂𝜀²𝑘,(2.2) where 𝜈=𝜇𝜈_𝑡𝜌,=𝐶^𝜇𝑘²𝐺_𝑆𝜀,=𝜈_𝑡2𝜕𝜕𝑥𝑢²+2𝜕𝑣𝜕𝑦²+𝜕𝑢𝜈_𝑘𝜕𝑦+𝜕𝜕𝑥𝑣²=𝜈+𝜈^𝑡𝜎_𝑘𝜀𝜈_𝜀,=𝜈+𝜈^𝑡𝜎_𝜀.,(2.3) Values of 𝐶_𝜇, 𝐶₁, 𝐶₂, 𝜎_𝜀 and 𝜎_𝑘 are selected as in Table 1.

Table 1 

Coefficients for standard 𝑘-𝜀 turbulence model.

Finally, the Reynolds-averaged Navier-Stokes equations (RANSEs) are used to model the turbulent flows as 𝜕𝑉𝜕𝑡+(𝑉⋅∇)𝑉=−∇𝜌𝑝+∇⋅𝜈+𝜈^𝑡∇𝑉+∇𝑉^𝑇+𝐵.(2.4) To model the free surface by VOF method, a step function of 𝑓(𝑥,𝑦,𝑡) is used. This function is expressed as [3] 𝜕𝑓𝜕𝑡+(𝜈⋅∇)𝑓=0.(2.5) Based on this equation, 𝑓 varies versus both time and space. Derivations of 𝑓 (and its continuous spectrum 𝐹) are used to determine the direction and curvature of the free surface. As direction, position and curvature of the free surface is determined, the surface tension can be computed. The continuous spectrum of 𝑓 is able to represented as 𝐹_𝑖(𝑡)=1𝐴_𝑖_𝐴𝑖𝑓_𝑖(𝑥,𝑡)𝑑𝑥.(2.6) Combining (2.5) and (2.6), the colour function is obtained as 𝜕𝐹𝜕𝑡+(𝜈⋅∇)𝐹=0,(2.7) where, 𝐹 has a specific value in each scalar cell as 𝐹=⎧⎪⎨⎪⎩10between0and1fiinnresseiiddseeuwarfiaartceercells.(2.8) Specific technique must be used to discrete (2.7).

3. Solution Algorithm and Stability Criteria

In this study, two-step projection method is used to solve the time-dependent NSE. The governing equations are discretized on a Cartesian-staggered grid system [21]. In this method, the two-step advancement algorithms are used as𝑈−𝑈^𝑛Δ𝑡=−conv^𝑛𝑈+Diff^𝑛+𝐵^𝑛,(3.1)_𝑛+1−𝑈^𝑛Δ𝑡=−∇𝑃^𝑛+1,(3.2) where 𝑈^𝑛 is the velocity field in old time level, 𝑈 intermediate velocity field, 𝑈^𝑛+1 the new velocity field, conv^𝑛 convection term, Diff^𝑛 diffusion term, and 𝐵^𝑛 the body force. The present method is explicit and first order in time, whereas the order of the space discretisation depends on the scheme used for the convective terms. The velocity field must satisfy the continuity equation. Therefore taking the divergence of (3.2), the Poisson’s equation for the pressure is obtained: ∇²𝑃^𝑛+1=1Δ𝑡∇𝑈.(3.3) So, at first, the intermediate velocity field is obtained by (3.1). Then, the pressure distribution in the new time level is obtained using (3.3). In this step, the new velocity field is calculated using (3.2). Finite difference and three-diagonal matrix algorithm (TDMA) methods are used to discretize and solve the Poisson’s equation.

The position of the free surface is then updated using VOF method expressed by (2.7). In this procedure, proper selection of the time step (Δ𝑡) affects the accuracy of final results. Therefore, time step was selected based on two stability criteria [22, 23] as Courant and diffusion conditions, respectively as follows: Δ𝑡_𝑐=minminΔ𝑥^𝑖Δ𝑡_𝜈||𝑢_𝑖,𝑗||,min||Δ𝑣𝑦^𝑗_𝑖,𝑗||,=12⎡⎢⎢⎣𝜈_𝑒11/Δ𝑥^𝑖²+1/Δ𝑦_𝑗²⎤⎥⎥⎦.(3.4) Using the above-mentioned relationships, the time step must be calculated to consider the smaller one in the numerical simulation.

4. Definition of the New Advection Method

In solution procedure of governing equations, it needs to decide where to store the velocity components. staggered and collocated grids can be used to evaluate this problem. On staggered grids, the velocity components are stored at the cell faces, and the scalar variables such as pressure are stored at the central nodes. However, on collocated grids, all parameters are defined at the same location at the central nodes. The staggered grids method gives more accurate pressure gradient estimation. However, collocated grids method is more simpler for solving the equations [24].

In Youngs’ VOF (Y-VOF) which is a volume of fluid method, free surface is tracked based on fluxes between two neighboring cells [25, 26]. In the staggered grids, the velocity components are defined in the cell faces, and fluxes are estimated using those velocity components. Flux translation is estimated as volume of fluid passed from the faces with constant velocities of faces. While the particles of fluid between two adjacent faces have different velocities. In this paper, a new advection method based on velocity in the mass gravity of fluid in a cell is introduced as MYV method. At first, estimation is made for the interface orientation 𝛽. The interface within a cell is then approximated by a straight line segment with orientation 𝛽 as shown in Figure 1.

Figure 1 

Interface orientation in a free surface (𝑖,𝑗) cell.

Free surface cuts the cell in such a way that the fractional fluid volume is given by 𝐹(𝑖,𝑗). The geometry of the fluid resulting from this reconstruction is then used to determine the fluxes through any side on which the velocity is directed out of the cell. For example, flux from right cell face (𝐹_𝑟) for the cell shown in Figure 2(a) can be estimated as 𝐹_𝑟=⎧⎪⎨⎪⎩12𝑈(𝑖,𝑗)𝑑𝑡2−𝑈(𝑆𝑖,_𝑏𝑗𝑑)𝑥∗𝑑𝑡𝑆^𝑟𝐹(𝑖,𝑗)∗𝑑𝑥∗𝑑𝑦𝑑𝑦iiff𝑈𝑈((𝑖𝑖,,𝑗𝑗))𝑑𝑑𝑡𝑡<>𝑆𝑆_𝑏𝑏𝑑𝑥𝑑𝑥.(4.1) Similarly, it is assumed that fluid is passed from cell face with constant velocity equal to cell face velocity. So, in the developed model, fluxes are calculated based on horizontal and vertical velocity of mass center of fluid cells. Figure 2 shows new arrangement of velocity in MYV method.

(a)

(b)

(c)

(d)

Figure 2 

Definition of velocity arrangement in new advection algorithm in MYV method.

The new equations of horizontal and vertical velocity (𝑈_𝐶𝐺 and 𝑉_𝐶𝐺) in the mass center are estimated using equations summarized in Table 2.

Table 2 

Velocity calculation for new arrangement in MYV method.

For example, fluxes from cell faces in (4.1) are estimated using these new velocities as 𝐹_𝑟=⎧⎪⎨⎪⎩12𝑈^𝐶𝐺(𝑖,𝑗)𝑑𝑡2−𝑈^𝐶𝐺(𝑖,𝑗)∗𝑑𝑡𝑆_𝑏𝑑𝑥𝑆^𝑟𝐹(𝑖,𝑗)∗𝑑𝑥∗𝑑𝑦𝑑𝑦iiff𝑈𝑈_{𝐶𝐺𝐶𝐺}(𝑖,𝑗)𝑑𝑡<𝑆_𝑏𝑑𝑥(𝑖,𝑗)𝑑𝑡>𝑆_𝑏𝑑𝑥.(4.2)

5. Model Validation

To validate the modified model, a series of standard tests such as lid-driven cavity, sloshing problem, constant unidirectional velocity field, shear test, and dam break over a dry bed was carried out.

5.1. Lid-Driven Cavity

Lid-driven cavity is the fluid flow in a rectangular container which moves tangentially to itself and parallel to one of the sidewalls. Due to the simplicity of the cavity geometry, applying a numerical method on this flow problem in terms of coding is quite easy and straightforward. Despite its simple geometry, the driven cavity flow retains a rich fluid flow physics manifested by multiple counter rotating recirculation regions on the corners of the cavity depending on the Reynolds number. The boundary conditions and induced eddies are shown in Figure 3.

(a)

(b)

Figure 3 

Lid-driven cavity test: (a) boundary conditions; (b) induced eddies.

A sensitivity analysis is performed on the problem mesh size. For example, horizontal velocities for Re = 100 and 3200 and for different mesh sizes are shown in Figure 4. As can be seen in this figure, by increasing the mesh size, the results become independent to it.

(a)

(b)

Figure 4 

Lid-driven cavity test; mesh independency for (a) Re = 100; (b) Re = 3200.

The model is then performed for various Reynolds numbers at a range of 100 to 10000 and results were compared with those of Ghia et al. [27] which is regarded as the true solution of this problem. For example, the results for Re = 1000 and 10000, are presented in Figures 5 and 6. The existence of a good agreement between results can be seen in these figures.

(a)

(b)

Figure 5 

Lid-driven cavity test, comparison between results of horizontal velocity of present model and those of Ghia et al. [27] for (a) Re = 1000; (b) Re = 10000.

(a)

(b)

Figure 6 

Lid-driven cavity test, comparison between results of vertical velocity of present model and those of Ghia et al. [27] for (a) Re = 1000; (b) Re = 10000.

5.2. Sloshing Problem

Sloshing of a low-amplitude liquid wave under forced movement is another problem to test the interfacial flow solver. In this test, a rectangular tank with a width of 0.9 m and a water depth of 0.6 m was exposed to a horizontal periodic sway motion as 𝑋=0.002sin(5.5𝑡). Therefore, 𝑎_𝑥=−0.0605sin(5.5𝑡) was considered as exciting acceleration. The displacement of a node on the free surface in contact with right-hand sidewall (𝜂) was calculated by the model and compared with those of Nakayama and Washizu [28] in Figure 7. The existence of a good agreement between results can be seen in this figure.

Figure 7 

Comparison between new models result for sloshing problem and that of Nakayama and Washizu [28].

5.3. Simple Advection Test

Before attempting to couple the advection 𝐹 to solution of the momentum equations, a comparison of the volume-tracking schemes with analytical solution is performed to validate VOF solver. In this regard some standard tests are used.

5.3.1. Constant, Unidirectional Velocity Field

The simplest test involves advection of a geometric shape in the computational domain. In this test, the geometric shape should remain intact, and total amount of the fluid within the region should be conserved. The test examined here is a hollow box being translated by a uniform constant velocity field. This is chosen to highlight the problem existing in some recent methods [13]. The 2D Cartesian region as shown in Figure 9 has dimension of 1m×1m and is composed equally sized cells. A square fluid blocks of dimensions 0.1m×0.1m moves with equal horizontal and vertical velocities of 1 m/s towards the top right-hand corner of the computational domain. Figure 8 shows the fluid position computed using the Hirt and Nichols, YV and MYV methods, every 0.1 sec until 0.7 sec. The computational time step in these calculations is 1×10⁻³sec leading to a courant number of 0.1.

(a)

(b)

(c)

(d)

Figure 8 

Numerical results of advection test for validation of VOF model: (a) exact solution; (b) Hirt- Nichols; (c) YV; (d) MYV.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9 

Final shape of circle after 2000 steps forward and then 2000 steps backward for different mesh sizes using YV (left) and MYV (right) methods: (a) and (b) 𝑁𝑋=𝑁𝑌=50; (c) and (d) 𝑁𝑋=𝑁𝑌=100; (e) and (f) 𝑁𝑋=𝑁𝑌=200.

It is evident in these figures that the geometric shape of the translated hollow square in modified model has not been improved comparing with original one. It is due to this fact that the velocity field is constant in this physical problem. Therefore, the velocity of mass center and cell faces are the same.

5.3.2. Shear Test, the Rudman Vortex

The final advection test examined in this study employs a nonuniform vorticity field which stretch and shear free surface interfaces as fluid is translated through the computational domain. A 2D computational domain with dimensions of 3.1m×3.1m composed of 𝑁𝑥×𝑁𝑦 uniformly sized square cells is used in the shear test. The velocity field is specified as: 𝑉𝑈(𝑥(𝑥,,𝑦𝑦))==−𝐴𝐴scinos(𝜋(𝜋𝑥𝑥)c)osisn(𝜋(𝜋𝑦𝑦)),,(5.1) in which 𝐴 equals to 1 for first 𝑁 computational time step, and −1 for the second 𝑁 time steps. The initial fluid geometry is a circle of radius 0.2 m. Time step is selected by courant number of 0.25 based on the maximum velocity within the computational domain. A sensitivity analysis is performed on mesh sizes. The final shape of circle after 2000 steps forward and then 2000 steps backward for different mesh numbers are estimated using different VOF methods. The results are presented in Figure 9. The results show that the present methods improve the results comparing with original ones for coarse grids, but they have the same accuracy for fine grids. The results in three steps for 𝑁𝑋=𝑁𝑌=100 and for YV and MYV methods are presented in Figure 10 as a sample.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 10 

Results for shearing field using YV (left) and MYV(right) methods and 𝑁𝑋=𝑁𝑌=100: (a) and (b) after 1000 steps forward; (c) and (d) after 2000 steps forward; (e) and (f) after 2000 steps forward and then 2000 steps backward.

5.4. Dam Break over a Dry Bed

Another test problem used for free surface case is the collapse of water column over a dry bed. This problem was first studied and used as benchmark by the developers of SOLA-VOF (Nichols et al. [4]). It is a very useful benchmark providing extreme conditions to assess the numerical stability as well as the capability of the model to treat the free surface problem. In this test, a square computational domain with a length and height of 22.8 cm is set. A water column with the width of 𝐿=5.7 cm and height of 2𝐿 is assumed at the left-hand side of the computational domain surrounded by walls with no-slip boundary condition. The spatial step sizes in the horizontal and vertical direction are Δ𝑥=Δ𝑦=0.1×𝐿. The comparison of the computed results using present methods with the experimental data given by Martin and Moyce [29] is shown in Figure 11. In this figure, 𝑋 is the location of wave front in time 𝑡. It is seen that the developed models passes this test successfully as the calculated result are well agreed with experimental data.

Figure 11 

Dam breaking test, comparison of models’ results with that of Martin and Moyce.

To compare the results of different VOF algorithms, two error criteria as sum square error (SSE) and sum absolute error (SAE) are used. The results are summarized in Table 3.

Table 3 

Comparison of SSE and SAE errors of different VOF algorithms in Dam break test.

6. Discussion and Conclusion

In this paper, modified Youngs’ VOF algorithm denoted by MYV method is presented. In Youngs’ VOF algorithm, for staggered grids fluxes to neighboring cells are estimated based on cell face velocities. However in practice, these particles have a variable velocity between velocities of two adjacent cell faces. In the developed model, the velocity in mass center of fluid cell is estimated and used to calculate fluxes from cell faces. To validate the modified model, a serious of standard tests such as Lid-driven cavity, Sloshing problem, constant unidirectional velocity field, Shear test, and Dam break over a dry bed was carried out. The results showed that in some cases such as Constant, unidirectional velocity field, geometric shape of the translated hollow square in modified model has not been improved comparing with original one. It is due to the constant velocity field leading to equal velocities of mass center and cell faces in this physical problem. In some other cases such as Shear test, the modified model improves the results especially for coarse grids. In Dam break test, present method improves the results dramatically. Therefore, based on these results, the MYV method is suggested as a new VOF algorithm which models the free surface problems more accurately.