# Free Surface Flow Simulation by a Viscous Numerical Cylindrical Tank
## 3 Numerical Method
### VOF Method  
As the volume fraction used in the VOF method is a step function, the problem arises of how to advect without numerical diffusing, dispersing or wrinkling accompanying with its convection across a mesh. To avoid those problems, various VOF methods that can be categorized as Lagrangian and algebraic methods have been proposed. With the Lagrangian methods, the free surface should be firstly rebuilt at each computational time step based on the volume fraction with the interface reconstruction methods, like the SLIC method (Simple Line Interface Calculation, presented by [PERSON] and [PERSON] in [19]) and the PLIC method. Then, the reconstructed surface was advected explicitly. With the algebraic methods, such as the Thinc method reported in [20], the Flux-Corrected Transport method presented in [21] and the CICSAM method, the high order convection schemes and mathmatical limitor was used to constrain the flux across the cell faces. To avoid that complexity, algebraic methods are more suitable methods to track the free surface in the cylindrical coordinates.  
In order to solve the VOF convection equation, the CICSAM method based on the donor and acceptor method in the original VOF method by [PERSON] and [PERSON] in [13], and the high-resolution scheme was applied in this work. In the CICSAM method, volume fraction was first approximated with the implicit Crank-Nicholson method, and then an iterative correction procedure was applied to correct the unphysical values. To simplify the computation, [PERSON] et al. implemented a Jacobi  
Figure 2: Staggered grid arrangement in control cells: (**a**) grid arrangement at normal control cell and (**b**) grid arrangement at the control cells at the interface of the coarse and fine blocks. (The subscripts of \(P\) and \(u\) represent the indexes of the control cells).  
type dual-time stepping formulation approach in [22] by reformulating the original VOF convection equation as  
\[\frac{F^{\tau+1}-F^{\tau}}{\Delta\tau}+\frac{F^{\tau}-F^{n}}{\Delta t}=-\frac{1}{ 2}(\frac{1}{r}\frac{\partial ru^{\tau}}{\partial\tau}^{\tau}+\frac{1}{r}\frac{ \partial vE^{\tau}}{\partial\theta}^{\tau}+\frac{\partial wF^{\tau}}{\partial z }^{\tau}+\frac{1}{r}\frac{\partial ruF^{n}}{\partial r}^{n}+\frac{1}{r}\frac{ \partial vF^{n}}{\partial\theta}^{n}+\frac{\partial wF^{n}}{\partial z}) \tag{13}\]  
by approaching the implicit solution with an explicit method as it converges in the pseudo time \(\tau\). With the finite volume method, the discretized form of Equation (13) could be written as:  
\[\frac{F^{\tau+1}-F^{\tau}}{\Delta\tau}+\frac{F^{\tau}-F^{n}}{\Delta t}=-\frac{ 1}{2}\frac{1}{\Delta V}\int{(F^{\tau}_{f}+F^{n}_{f})U^{n+1}n_{f}\mathrm{d}A} \tag{14}\]  
where \(F^{\tau}_{f}\) and \(F^{n}_{f}\) are volume fractions defined at the faces of the control cell, which can be calculated by the CICSAM scheme, and \(\Delta\tau\) is the pseudo time step set to be \(2\Delta t/3\).  
To avoid numerical oscillations at the interfaces between two blocks, the flux-constraint technique used in the velocity computation was adopted to keep the flux of the volume fraction the same in this study.