# Free Surface Flow Simulation by a Viscous Numerical Cylindrical Tank
## 3 Numerical Method
### Discretization of the Continuity and Momentum Equations  
To solve the control equations, the staggered grid arrangement was applied to store velocities at the cell faces and the scalar variables (pressure, density, viscosity, etc.) at the cell center (Figure 2), and the velocity and pressure field was coupled with the projection method, where an intermediate velocity field \(U^{*}\) was computed and projected into a divergence-free field with Equations (7a) and (7b), as  
\[\frac{U^{*}-U^{n}}{\Delta t}= -\frac{1}{r}\frac{\partial(rtU)^{n}}{\partial r}+\frac{1}{r} \frac{\partial(rtU)^{n}}{\partial\theta}+\frac{\partial(uU)^{n}}{\partial \theta} \tag{7a}\] \[+\frac{1}{r}\frac{\partial}{\partial r}(vr\frac{\partial U}{ \partial r})^{n}+\frac{1}{r}\frac{\partial}{\partial\theta}(v\frac{\partial U }{\partial\theta})^{n}+\frac{\partial}{\partial x}(v\frac{\partial U}{ \partial\theta})^{n}+\frac{\partial}{\partial x}(v\frac{\partial U}{\partial \theta})^{n}+S_{U}^{n}\] \[U^{n+1}=U^{*}-\frac{\Delta t}{\rho^{n+1}}\
abla(P^{n+1}-P^{n})^{ P^{\prime}=P^{n+1}-P^{n}}\
abla\cdot U^{*}=\Delta t\
abla.\frac{1}{\rho^{n+1}} \
abla P^{\prime} \tag{7b}\]  
where \(n\) and \(n+1\) donate the previous and new time layers. The intermediate velocity could be computed by discretizing the convection, diffusion and source terms in Equation (3a) over the control volumes by using the finite volume method. For the orthogonal grids, by integrating the convection term over the control volumes, it yielded  
\[\int(\frac{1}{r}\frac{\partial(ruU)}{\partial r}+\frac{1}{r}\frac{\partial(vU )}{\partial\theta}+\frac{\partial(wU)}{\partial z})\mathrm{d}V=\int U_{f}Un_{ f}\mathrm{d}A \tag{8}\]  
in which \(n_{f}\) is the unit surface vector of cell faces and \(U_{f}\) is the velocity at the cell faces of the control volume, which could be determined by a high-resolution scheme as the Sharp and Monotonic Algorithm for Realistic Transport (SMART) scheme by [PERSON] in [PERSON] and [PERSON] (1986).  
Integrating the diffusion term in Equation (7a) over the control volume, it was given as  
\[\int(\frac{1}{r}\frac{\partial}{\partial r}(vr\frac{\partial U}{\partial r} )+\frac{1}{r}\frac{\partial}{\partial\theta}(\frac{v}{r}\frac{\partial U}{ \partial\theta})+\frac{\partial}{\partial z}(v\frac{\partial U}{\partial z})) \mathrm{d}V=\int v\
abla U_{f}n_{f}\mathrm{d}A, \tag{9}\]  
where the gradient of the velocity at the faces of control volume was computed by the central difference scheme. The source term in Equation (7a) was assumed to be piecewise uniform over each volume, discretized as  
\[\int S_{U}\mathrm{d}V= S_{U}\Delta V \tag{10}\]By invoking the [PERSON]'s theorem, Equation (7b) could be written as  
\[\int U\pi_{f}\mathrm{d}A=\frac{\Delta t}{\rho^{n}}\int\
abla P^{r}n_{f}\mathrm{d}A \tag{11}\]  
where \(\
abla P^{r}\) or \(\
abla P\) at the normal faces is evaluated with a central difference scheme, and special concerns should be paid in the treatment of \(\
abla P^{r}\) or \(\
abla P\) at the interface between two blocks. [PERSON] et al. simulated a two-phase flow motion over an adaptive octree grids system in the three-dimensional Cartesian Coordinates by treating the gradient at the interface to be interpolated over cells that shared the same faces in [18], and the interpolation scheme used in that study was introduced to the present algorithm. Therefore, the pressure gradient at the west faces of the control cells \(EN\) and \(ES\) shown in Figure 2b should be kept the same, which could be formulated as  
\[\frac{\partial P}{\partial r}_{ES}=\frac{\partial P}{\partial r}_{EN}=\frac{(P _{ES}+P_{ES})-2P_{P}}{2\Delta r} \tag{12}\]  
This interpolation yields a symmetric coefficient matrix for the pressure, which enables the use of the iterative solver, such as the preconditioned conjugate gradient (PCG) method. In this paper, the pressure Poisson equation was solved by using an algebraic multi-grid method (AMG).