Ying-wei SUN, Hai-gui KANG*

School of Hydraulic Engineering, Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian 116023, P. R. China

Application of CLEAR-VOF method to wave and flow simulations

Ying-wei SUN, Hai-gui KANG*

School of Hydraulic Engineering, Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian 116023, P. R. China

A two-dimensional numerical model based on the Navier-Stokes equations and computational Lagrangian-Eulerian advection remap-volume of fluid (CLEAR-VOF) method was developed to simulate wave and flow problems. The Navier-Stokes equations were discretized with a three-step finite element method that has a third-order accuracy. In the CLEAR-VOF method, the VOF functionFwas calculated in the Lagrangian manner and allowed the complicated free surface to be accurately captured. The propagation of regular waves and solitary waves over a flat bottom, and shoaling and breaking of solitary waves on two different slopes were simulated with this model, and the numerical results agreed with experimental data and theoretical solutions. A benchmark test of dam-collapse flow was also simulated with an unstructured mesh, and the capability of the present model for wave and flow simulations with unstructured meshes, was verified. The results show that the model is effective for numerical simulation of wave and flow problems with both structured and unstructured meshes.

water wave; water flow; numerical simulation; CLEAR-VOF;three-step finite element method

1 Introduction

Modeling unsteady free surface flows is an important step in study of water wave problems. At present, the volume of fluid (VOF) method is the most popular free surface modeling method. But the original VOF method, which is based on the finite difference method, is limited to problems with complicated computation domains. In order to make up for this limitation, novel VOF methods have been developed by many researchers. Jeong and Yang (1998) improved the method for calculating the VOF function and presented a numerical model that combined the VOF method with finite element analysis and could be applied to adaptive meshes. L?hner et al. (2006) developed a new VOF model, which can be operated on adaptive unstructured meshes to simulate the interactions between extreme waves and three-dimensional structures. Yang et al. (2006) simulated interfacial flow with the VOF method on unstructured triangular meshes by combining an adaptive coupled level set.Mencinger and ?un (2011) presented a variant of the VOF model, which is based on the piecewise linear interface calculation method and can be used on general moving meshes. Ashgriz et al. (2004) presented a novel VOF method, namely the computational Lagrangian-Eulerian advection remap-VOF (CLEAR-VOF) method. In this method, the volume fraction of each element is calculated in the Lagrangian manner with geometric tools, and the piecewise linear reconstruction method is used to reconstruct the interface. With the help of the finite element method, this method can be easily applied in unstructured meshes.

In this study, a two-dimensional numerical model was established by combining the Navier-Stokes equations with the CLEAR-VOF method. The three-step finite element method (Jiang and Kawahara 1993), which has a third-order accuracy, was used to discretize the Navier-Stokes equations. Using this model, the propagation of regular waves and solitary waves over a flat bottom, and shoaling and breaking of solitary waves on two different slopes were simulated. Numerical results were compared with experimental data and theoretical solutions. To verify the performance of the model on unstructured meshes, a benchmark test of dam-collapse flow was also carried out.

2 Governing equations and discretization

2.1 Governing equations

The governing equations are the following continuity equation and two-dimensional Navier-Stokes equations for incompressible viscous fluid:

whereuis the velocity; the subscriptsiandjrepresent coordinate directions;ρis the fluid density and equals 1000 kg/m3;pis the pressure;fis the unit body force and equals 0 m/s2in thexdirection and –9.8 m/s2in theydirection, respectively;νeis the effective viscosity coefficient and can be described asνe=ν+νt, whereνis the kinematic viscosity and equals 1.002 × 10-6m2/s, andνtis the turbulent viscosity. Here, the Smagorinsky subgrid-scale turbulence model is employed:

wherecsis the Smagorinsky constant ranging from 0.1 to 0.2, andΔis the characteristic length of the element.

2.2 Three-step finite element method

The Navier-Stokes equations are discretized with the three-step finite element method inwhich the time-splitting idea is employed to discretize the time term. As shown in Eq. (4), the Taylor series expansion of velocityuas a function oftcan be described as

wherendenotes the number of time steps. By approximating Eq. (4) to a third-order accuracy, the three-step formula for velocityucan be written as

By applying the three-step formula to Eq. (2), the following discretized Navier-Stokes equations can be obtained:

It can be seen from Eq. (10) thatpn+1should be obtained in advance to calculateun+1. By applying divergence to both sides of Eq. (10) and introducing the incompressible constraint, the following Poisson equation of pressure is obtained:

By applying the standard Galerkin finite element method, the weak form of the momentum equations and the Poisson equation are derived:

wherewis the weight function,Vis the volume of elements,uiβis the velocity at nodeβin directioni,ujγis the velocity at nodeγin directionj, andφβandφγare the basis functions. It should be noted that the same order basisfunctions are used for the interpolation of velocities and pressures. In the present study, the quadrilateral isoparametric element is used, soβandγrange from 1 to 4, and the basis functions are as follows:

whereζandψrepresent the local coordinates of the standard element and range from –1 to 1.

The explicit momentum equations are solved using the lumped-mass matrix method, and the Poisson equation is solved using a preconditioned bi-conjugate gradient method.

3 CLEAR-VOF method

3.1 Calculation of VOF function

In the original VOF method, the VOF functionFis updated every time step by solving its transport equation. Since the CLEAR-VOF method is also applied in Eulerian fixed meshes, and the volume of fluid advection is calculated in the Lagrangian manner with geometric tools, it is unnecessary to solve the transport equation. The procedure for calculating theFfunction with the CLEAR-VOF method can be divided into three steps:

(1) Construction of fluid polygons for elements: In this step, the vertices of fluid polygons in each element need to be determined, as shown in Fig. 1. IfF= 0, there is no fluid polygon because no fluid exists in this element. If 0 ＜F＜ 1, the element is partially filled with fluid: the nodesa,b, anddof the element are inside the fluid polygon, and nodeseandgare the intersections of the interface and the background grid. The nodesa,b,e,g, anddconstitute the vertices of the fluid polygon. IfF= 1, the element is fully filled with the fluid, and the fluid polygon is identical to the background grid. The nodesa,b,c, anddof the grid are also the vertices of the fluid polygon.

Fig. 1 Vertices of fluid polygon

(2) Polygon movement calculation: Once a polygon is identified, its movement is calculated in the Lagrangian manner.andare the velocity components of the vertices of the fluid polygon at thenth time step. After a time step Δt, the displacements and new locations of the vertices can be obtained by Eqs. (17) and (18):

As shown in Fig. 2, nodesa,b,c,d, andein Fig. 2(a) and nodesa1,b1,c1,d1, ande1in Fig. 2 (b) are the vertices of the fluid polygon atthenth time step and the (n+1)th time step, respectively.

Fig. 2 Sketch of CLEAR-VOF method

(3) VOF function calculation at a new time step: The new fluid polygona1b1c1d1e1intersects with the background grids. A portion of the fluid remains inside the home element, and the rest enters its neighboring elements. In Fig. 2, the home element has eight neighboring elements. The fluid-covered areas in the home element and other neighboring elements are calculated, and the new VOF functionFcan be derived by repeating the procedure for all the fluid polygons.

3.2 Interface reconstruction

The piecewise linear reconstruction method (Ashgriz et al. 2004) is used to reconstruct the interface. For an interface elementIwithJneighbouring elements, the unit normal vectormof the line interface, expressed asm=mxi+myj, can be obtained by

where ? is the differential operator. ?Ffor elementIis solved by the following linear system:

where

where (xcen,ycen) is the location of the centroid of each element. Then the interface equation can be obtained by

whereDis a constant.

4 Numerical results

4.1 Regular wave simulation

A two-dimensional numerical wave flume model was established. The domain was 10 m long and 0.8 m high. The computational domain and coordinate system are shown in Fig. 3. The origin of thex-axis was fixed at the inflow boundary on the left side and the origin of they-axis was fixed at the bottom boundary. A mesh of 800 elements in thexdirection and 150 elements in theydirection was used.

Fig. 3 Sketch of computational domain for traveling wave

On the left side, the numerical wave generation theory (Dong and Huang 2004) was adopted to generate second-order Stokes waves. With a piston-type wave board, the horizontal velocityv(t) of the wave paddle satisfies

whereHis the wave height, andLis the wave length.

On the right side, a damping layer boundary condition (Larsen and Dancy 1983) was adopted to absorb the wave energy gradually. In the damping layer, the velocity components are divided byμ(x). The damping factorμ(x) was described as

whereαis equal to 1.1,Lwis the length of the computational domain, andλis the thickness of the damping layer and is set as the wave length.

Fig. 4 and Fig. 5 show the time histories of the simulated regular wave cases with a still water depth (h0) of 0.35 m, wave heights (H) of 0.10 m and 0.12 m, and wave periods (T) of 1.2 s and 1.3 s, respectively. After the first three or four waves, the temporal curves of the simulated waves are very stable. The phases and amplitudes of the waves generated by the numerical model are in good agreement with the theoretical solutions. The nature of the second-order Stokes wave with the wave shape asymmetrical to the still water level is well demonstrated in the numerical results.

Fig. 4 Temporal curves of wave free surface (H= 0.10 m, andT= 1.2 s)

Fig. 5 Temporal curves of wave free surface (H= 0.12 m, andT= 1.3 s)

4.2 Solitary wave simulation

To generate a solitary wave, the following displacement of a piston-type wave maker (Goring and Raichlen 1980) is given:

whereis the initial heig ht of the solitary wave,cis the wave speed and can be calculated

The velocity of a wave maker is described as follows:

wherexis the position of the wave maker.

4.2.1 Propagation of solitary wave

The configuration of the flume is similar to the previous case shown Fig. 3. The flume is 18 m long and 0.5 m deep. The domain is divided by 1 440 elements in thexdirection and 90 elements in theydirection. A solitary wave for, propagating rightward over a water surface with a constant still water depth ofh0=0.3m was simulated.

Fig. 6 Comparison of solitary wave profile

Fig. 6 shows the solitary wave profile generated by the present numerical model using Boussinesq’s theory att =2.5 s. The numerical result is shown to be very close to the analytical solution.

Fig. 7 shows the solitary wave propagation at four moments in the numerical flume. It can be seen that the solitary wave profile is stable in a permanent form during propagation. The attenuation of the solitary wave height during the propagating process was small. According to Mei’s perturbation theory (Mei 1989), the solitary wave height during propagation can be determined by

Fig. 8 shows thatcalculated by the present model agrees with Mei’s perturbation theory.

Fig. 7 Propagation of solitary waves in numerical flume

Fig. 8 Attenuation of solitary wave height during propagation

4.2.2 Shoaling and breaking of solitary waves on slopes

The shoaling and breaking of solitary waves on different slopes were computed with thepresent model. Fig. 9 shows the computational domain.

Fig. 9 Computational domain for solitary wave shoaling and breaking on slopes

In the first case, the slope was 1:15, andwas 0.3. In the second case, the slope was 1:35 andwas 0.2. Figs. 10 and 11 show the comparison of the numerical results of the present method with those of Grilli et al. (1997), in which an experimentally validated fully nonlinear wave model was used to simulate the shoaling and breaking processes of solitary waves. The coordinates in both the horizontal and vertical directions were normalized by the still water depth, andt1,t2,t3, andt4rep resent four different moments. It can be seen that the shoaling and breaking processes in the present study were quite similar to those of Grilli et al. (1997) in both cases. The breaking points also matched well. The water tongues in the present study were a little thicker than those in Grilli et al. (1997).

Fig. 10 Free surface profiles of solitary wave shoaling and breaking on 1:15 slope (=0.3)

Fig. 11 Free surface profiles of solitary wave shoaling and breaking on 1:35 slope (=0.2)

The maximum wave runup (Rm) of a solitary wave on a 1:1 slope was also investigated and a comparison of results is shown in Fig. 12. It can be seen that the results of the present study lie between those of Synolakis (1987) and Maiti and Sen (1999).

Fig. 12 Maximum runup of solitary waves on 1:1 slope

4.3 Dam-collapse flow simulation

Dam-collapse flow simulation, which is a classic-benchmark test for the verification of a free surface model, was conducted with an unstructured mesh to verify the performance of the CLEAR-VOF method. As shown in Fig. 13, an initial water column with both the height (h1) and width (h2) of 0.055 m is determined on the left of the flume. The computational domain, which is 0.3 m long and 0.06 m high, is divided by an unstructured quadrilateral mesh with 2 172 elements and 2 313 nodes.

Fig. 13 Computational domain with unstructured quadrilateral mesh

Fig. 14 shows the evolution of the free surface at the instants oft= 0.06 s,t= 0.08 s,t= 0.10 s, andt= 0.12 s. It can be seen that the interface can be captured by the present model with an unstructured mesh.

Fig. 14 Evolution of free surfaces

Fig. 15 and Fig. 16 show the comparisons of numerical results and experimental data (Martin and Moyce 1952) for the variation of the water surface elevationYon the left side wall and the position of water frontXalong the bottom, respectively. The dimensionless coordinates are defined as

The comparisons of the numerical results and experimental data are fairly acceptable, demonstrating the validity of the present model on unstructured meshes.

Fig. 15 Variation of free surface elevation on left side wall

Fig. 16 Variation of water front along bottom

5 Conclusions

Based on the Navier-Stokes equations and the CLEAR-VOF method, a two-dimensional numerical model for wave and flow simulations was developed in this study. The three-step finite element method was used to discretize the Navier-Stokes equations.

The model was used to simulate the propagation of regular waves and solitary waves over a flat bottom, and shoaling and breaking of solitary waves on different slopes, and a benchmark test of dam-collapse flow with an unstructured mesh was also carried out. The result shows that, compared with the original VOF method that could be only applied in structured meshes, the CLEAR-VOF method can be applied in both structured and unstructured meshes due to its combination with the finite element method. The present model provides a convenient way for simulating wave and flow problems with complex boundaries.

Ashgriz, N., Barbat, T., and Wang, G. 2004. A computational Lagrangian-Eulerian advection remap for free surface flows.International Journal for Numerical Methods in Fluids, 44(1), 1-32. [doi:10.1002/fld.620]

Dong, C. M., and Huang, C. J. 2004. Generation and propagation of water waves in a two-dimensional numerical viscous wave flume.Journal of Waterway, Port, Coastal and Ocean Engineering, 130(3), 143-153. [doi:10.1061/(ASCE)0733-950X(2004)130:3(143)]

Goring, D., and Raichlen, F. 1980. The generation of long waves in the laboratory.Proceedings of the 17th International Coastal Engineering Conference, 763-783. New York: ASCE.

Grilli, S. T., Svendsen, I. A., and Subramanya, R. 1997. Breaking criterion and characteristic for solitary waves on slopes.Journal of Waterway, Port, Coastal and Ocean Engineering, 123(3), 102-112. [doi:10.1061/ (ASCE)0733-950X(1997)123:3(102)]

Jeong, J. H., and Yang, D. Y. 1998. Finite element analysis of transient fluid flow with free surface using VOF(Volume of Fluid) method and adaptive grid.International Journal for Numerical Methods in Fluids, 26(10), 1127-1154. [doi:10.1002/(SICI)1097-0363(19980615)26:10<1127::AID-FLD644>3.0.CO;2-Q]

Jiang, C. B., and Kawahara, M. 1993. The analysis of unsteady incompressible flows by a three-step finite element method.International Journal for Numerical Methods in Fluids, 16(9), 793-811. [doi:10.1002/ fld.1650160904]

Larsen, J., and Dancy, H. 1983. Open boundaries in short wave simulations: A new approach.Coastal Engineering, 7(3), 285-297. [doi:10.1016/0378-3839(83)90022-4]

L?hner, R., Yang, C., and O?ate, E. 2006. On the simulation of flows with violent free surface motion.Computer Methods in Applied Mechanics and Engineering, 195(41-43), 5597-5620. [doi:10.1016/ j.cma.2005.11.010]

Maiti, S., and Sen, D. 1999. Computation of solitary waves during propagation and runup on a slope.Ocean Engineering, 26(11), 1063-1083. [doi:10.1016/S0029-8018(98)00060-2]

Martin, J. C., and Moyce, W. J. 1952. An experimental study of the collapse of liquid columns on a rigid horizontal plane.Philosophical Transaction of the Royal Society of London,Series A, Mathematical and Physical Sciences, 224(882), 312-324.

Mei, C. C. 1989.The Applied Dynamics of Ocean Surface Waves. Singapore City: World Scientific.

Mencinger, J., and ?un, I. 2011. A PLIC–VOF method suited for adaptive moving grids.Journal of Computational Physics, 230(3), 644-663. [doi:10.1016/j.jcp.2010.10.010]

Synolakis, C. E. 1987. The runup of solitary waves.Journal of Fluid Mechanics, 185, 523-545. [doi:10.1017/ S002211208700329X]

Yang, X. F., James, A. J., Lowengrub, J., Zheng, X. M., and Cristini, V. 2006. An adaptive coupled level-set/volume-of-fluid interface capturing method for unstructured triangular grids.Journal of Computational Physics, 217(2), 364-394. [doi:10.1016/j.jcp.2006.01.007]

This work was supported by the National Natural Science Foundation of China (Grant No. 50679008).

*Corresponding author (e-mail:hgkang@dlut.edu.cn)

Received Oct. 24, 2011; accepted Jan. 11, 2012