Benchmark on the Dynamics of Liquid Draining Inside a Tank

. Immense information and details observation of flow physics inside a draining tank can be achieved by adopting reliable numerical simulations. Yet the accuracy of numerical results has been always debatable and it is mainly affected by the grid convergence error and computational modeling approaches. Hence, this study is divided into two stages. In the first stage, this paper determines a systematic method of refining a computational grid for a liquid draining inside a tank using OpenFOAM software. The sensitivity of the computed flow field on different mesh resolutions is also examined. In order to study the effect of grid dependency, three different grid refinements are investigated: fine, medium and coarse grids. By using a form of Richardson extrapolation and Grid Convergence Index (GCI), the level of grid independence is attained. In this paper, a monotonic convergence criteria is reached when the fine grid has the GCI value below 10% for each parameter. In the second stage, different computational modeling approaches (DNS, RANS k-ε, RANS k-ω and LES turbulence models) are investigated using the finer grid from the first stage. The results for the draining time and flow visualization of the generation of an air-core are in a good agreement with the available published data. The Direct Numerical Simulation (DNS) seems most reasonably satisfactory for VOF studies relating air-core compared to other different turbulence modeling approaches .


Introduction
The use of cylindrical tank has grown over the past three decades.The liquid draining tank is one of the most command tanks that is frequently used in various engineering applications.The draining tank is important to conserve a system run continuously without any shortage in the supply of the working liquids.The use of liquid draining tanks can be found on the space vehicle for its propellant storage, in the water tank for the domestic use and in the oil industries for its short and long term storages at the refineries.Usually the process of liquid draining from the tank is ordinarily made by pressured system and/or drain by gravity, downward.
The generation of an air-core vortex is a flow phenomenon that frequently occurs inside the liquid draining tanks.Air-core vortex is one of the rotational motions of the liquid at which air entraining the vortex through its core [1] and it can be detected along the draining process inside the tank.As the level reaches a certain critical height, Hc, the dip forms into the air-core on the surface which subsequently enters the outlet [2].The air-core vortex phenomenon continuously been generated at which the dip grows into a vortex with an air-core and it creates a long and slender string shape lengthens to the bottom of the tank [3].This air-core vortex formation is deepened by the intensification of rotational flow along the draining process [2].The rate of liquid draining is decreased and the flow at the outlet nozzle is highly rotational when the core of the vortex extents to the bottom of the tank [1].
In the scope of benchmarking the generation of aircore vortex inside the tank, there are mainly two factors that may lead to inconsistencies of the numerical results.Firstly, is due to the mesh size sensitivity.According to Mathew et al. [4], a grid independence study is implemented to find the most appropriate mesh size.The mesh must be a convergent solution, stable and consistent.A higher mesh density is selected in order to resolve the sharp gradients from the draining effects.Hence, when the air-core is formed, the interface between liquid and gas is better resolved.Son et al. [5] have studied the grid resolution with three different grid systems.In the coarse (21x20x200) and medium (84x20x200) grid system, the air-core is not reproduced.Only after the grid system is doubled from medium in all directions, the air-core finally is reproduced in the finest grid (42x40x400).Hyun et al. [3] in their studies have proposed various grid refinement for the grid dependency.In order to produce a realistic solution, they have discovered that increasing the number of nodes in the radial direction is more efficient than increasing the number of nodes in the axial direction.
Secondly, the difference in computational modelling approaches has also contributed to the deviation in the results.Madsen et al. [6] have studied the applicability of two-phase CFD modelling of internal flow in a largescale pressure-swirl atomizer by applying three methods: 1) a VOF method employing a DNS model, 2) a two-fluid Euler/Euler employing a DNS model, and 3) a VOF method employing LES turbulence model.Similar results to the experiments are achieved and all cases are proven to be able to reproduce the characteristic of aircore vortex.However, the two-fluid and VOF method with DNS appears to give the best agreement with the experimental results.In addition, Hyun et al. [3] have compared the simulation of flow inside the liquid draining tank using turbulence models and DNS.Referring to draining time data, a large difference between turbulence and DNS is observed.A much better agreement is achieved with the DNS results than using the turbulence model.
CFD community has been well defined and agreed that the error from the numerical simulation is not just from the grid convergence error and the selection of computational modelling.But, reducing the error due to grid dependency and properly construct the computational setting, one can minimize the total error.This brings to the objective of the study, which will concentrate on the assessment of the grid dependence and computational modelling approaches.These must be done in a much systematic method in order to reproduce the generation of air-core inside the tank.The fundamental physics flow of the air-core is also reexamined in this study.
This paper is divided into two stages.In the first stage, three computational meshes based on reasonable approximations of cell sizes are created.The grid sensitivity to the parameter drain time,   and flow rate, Q is assessed.Richardson extrapolation and Grid Convergence Index (GCI) have been used to estimate the grid convergence error due to the sensitivity of the resolution [7][8][9].In the second stage, using the finer grid of the first stage, different computational modelling approaches (DNS, RANS k-ε, RANS k-ω and LES turbulence models) are evaluated.
In this study, a model is intentionally made the same as the experimental and numerical investigations of Park and Sohn [10] so that a comparison study can be made.Additionally, this model is also a general model for a draining tank, where a benchmark study can be adopted.The cylindrical tank of length (L) 450mm and diameter (D) of 90mm is filled partially with water.The initial height of the water measured from the bottom of the tank (ho) is 350mm.A drain nozzle of diameter (d) 6mm and length (l) of 15mm is located at the center of the bottom surface of the tank.The bottom and top of the tank are in the atmospheric condition and the fluid is drained naturally by gravity, g, downward.In order to reproduce the generation of the air-core, initial rotation with the speed of 120 RPM is imparted to the wall of the tank before the liquid is started to drain [10].Figure 1 shows the schematic diagram of the draining tank.

Solution methodology 2.1 Governing equations
The continuity and momentum equations are used to model the incompressible draining flow under the influence of gravity, g.Here   is the velocity vector,   is the position vector of coordinate system,  is the pressure,  is the time,  is the dynamic viscosity,  is the density,   is the gravitational acceleration,   is the volumetric surface tension force.
The Euler scheme (1 st order accurate) has been adopted in the first stage.In second stage, it is changed to backward scheme (2 nd order accurate) to reduce the discretization error.

Mesh descriptions
In this grid refinement study, three grid resolutions are tested.Case A has the finest grid, Case B has medium resolution and Case C has the coarsest resolution.Table 1 shows the grid parameter for the Case A, B and C. Referring to Celik et al. [11], the refinement ratio must be greater than the minimum value of 1.3.In this present study, the grid refinement ratio,  for the uniformed meshed is given as below:

Richardson Extrapolation (RE)
Richardson Extrapolation by Roache [8] also known as "the deferred approach to the limit (ℎ → 0)".It defines a higher-order estimate of flow fields from a series of lower-order discrete values ( 1 ,  2 , . .,   ).The equation of Richardson Extrapolation is given as below: Where,  is discrete solutions,  1,  2 are functions and ℎ is grid spacing.If  1 = 0, the quantity  is considered "second-order".The  ℎ=0 is the continuum value at zero grid spacing.As shown in equation ( 4), Richardson Extrapolation can also be generalized to  ℎ order methods and -value of grid ratio without considering the non-appearance of odd power as in equation ( 5): According to Stern et al. [9], equation ( 5) can be measured for order-of-accuracy by applying equation ( 6 1.As the value of convergence ratio (R) is more than zero and less than one, the convergence conditions for parameter draining time,  and flow rate,  are monotonic.Hence, these two parameters are relevant in the Grid Convergence study.

Grid Convergence Index (GCI)
Grid Convergence Index (GCI) method by Roache [8] is based on grid refinement error estimator derived from the generalized Richardson Extrapolation (RE).By using this method, the percentage of differences between computed value and asymptotic value are calculated.It also determines how far the errors of the computed value with the asymptotic value and how much the solution of computed values would transform with further refinement.A small value of GCI percentage displays the computed value is approaching asymptotic range.The GCI for fine grid can be explained as given equation ( 9): Here,   is safety factor.The safety factor should be   = 1.25 [12] since three grids resolutions are used in the present study.
Table 2 shows the Grid Convergence Index (GCI) and order of accuracy for parameter draining time,  and flow rate,  from three different meshes.The GCI value is reduced when the mesh is refined.When the GCI for the finer grid ( 21 ) is lower than the coarser grid ( 32 ), it shows that the dependency of numerical method on the mesh size has been decreased.According to the results listed in Table 2, the changes in the simulation results when the cell size is successively refined, are shown in figure 2 and figure 3 for parameter draining time,   and flow rate, , respectively.Both parameters show a convergence solution to a value near the Richardson extrapolated value.Additionally, the Richardson extrapolated value is within the range of calculated GCI value of 8.75% and 3.98% for,   and flow rate, , respectively.

Comparison of drainage time
In the first stage of the study, the results from the finest grid resolution have shown a satisfactory level of grid independence.The second stage in the study will be using the finer grid resolution for the assessment of various turbulence models in simulating the liquid draining tank.Three turbulence model ( −  ,  −  and LES) are assessed in addition to the Direct Numerical Simulation (DNS).The sensitivity of the results of the time discretization scheme is also assessed for all turbulence models.Table 3 compares the result of the current study with the similar study by Park and Sohn [10].The theoretical value that is calculated from equation ( 10) is also compared.According to Park and Sohn [10], the draining time can be expressed theoretically with the following equation: Here,  is the draining time, ℎ  the initial water level, ℎ the water level at time ,   the tank diameter,   the nozzle diameter, and  gravitational acceleration.Equation ( 10) is derived with the assumption that the flow is irrotational and inviscid.
In the case of 1 st order of time discretization scheme, the draining time completion obtained from the current simulation using DNS, RANS  − , RANS  −  and LES are 91s, 72s, 70s and 101.5, respectively.The draining times for RANS  −  and RANS  −  are 8.26s and 10.26s earlier than the result obtained from the experiment by Park and Sohn [10].Meanwhile, the draining time completion obtained from DNS and LES are 10.74s and 21.24s slower than the result obtained from the experiment by Park and Sohn [10].Not much changes are observed when the time discretization is changed to 2 nd order scheme.

Flow visualization of air-core formation
Figures 4-7 show the progression of liquid draining obtained from DNS, RANS  −  , RANS  −  and LES, respectively.At the beginning of the draining process (t=0), the top surface of the liquid is in parabolic shape.This is due to the centripetal force from the initial wall rotation and the density difference between liquid and air.The parabolic shape for RANS  −  and RANS  −  are more obvious than DNS and LES due to the Reynolds stress term from the convective acceleration which effects on the mean flow [13].However, at time 15s of draining process, the top surface of the liquid is flat for the case RANS  −  and RANS  −  .The shape remains the same until the end of the draining process.A different shape is observed for DNS and LES where a dip is observed near the centre of the tank.As the draining process continues, the dip extends till the outlet of the tank (happens at t=21s) at which the air-core generation is fully completed.At this moment, the dip raises into a vortex with an air-core and the free surface creates a long and slender string shape lengthens to the bottom of the tank, and it is named as air-core vortex [3].As the level reaches a certain critical height Hc, the dip forms into the air-core on the surface which consequently enters the outlet [2].This phenomenon repetitively continued until the draining finished (except when the reverse jet is occurred).The flow at the outlet nozzle is highly rotational and the rate of liquid draining is decreased when the core of the vortex extents to the bottom of the tank [1].[5].According to Son et al. [5], the liquid in the Taylor vortices cannot be combined with the rotating axially in the centre of the tank since it limits the liquid in the off area.So, at the first, only the liquid in the centre of the tank is drained out.The stack structures of Taylor vortices act like a block and makes the condition where a shallow water drain, even though the water level is considerably lower.Thus, the dimple on the free surface is pulled downward and finally, the air-core is reproduced.In the figure 9 and 10, there are no Taylor vortices have been discovered in order to accelerate the axially rotating vortex to regenerate the air-core.Hence, in the case of RANS  −  and RANS  − , there are no generation of air-core is observed.

Conclusions
For the first stage of the study, the grid refinement study on three different meshes was successfully accomplished using OpenFOAM framework.As the grid was refined, the value of Grid Convergence Index (GCI) for parameter of draining time,   and flow rate,  from coarser to finer grid was decreased.The finer grid (Case A) is applied in the second stage of the study because it has the lower GCI and the GCI value is only 3%.
For the second stage, the current simulations (DNS, RANS  − , RANS  −  and LES models) are able to reproduce the liquid draining process inside the tank in all cases.From the comparison of drain time plots, the current DNS demonstrates a very similar pattern and value with the result obtained from the experimental measurement of Park and Sohn.The ellipsoidal shape of the free surface was successfully recreated in all cases at the beginning of the liquid draining.In the DNS and LES cases, the finer grid from the first stage is successfully reproduced the generation of air-core and it is in a good agreement with the result of Son et al.Thus, based on the results from the first and second stage, DNS is most reasonably satisfactory for VOF studies relating air-core compared to other different modelling approaches.

Fig. 2 .
Fig. 2. Comparison of parameter draining time, t drain between three different meshes and Richardson Extrapolation estimation

Fig. 3 .
Fig. 3. Comparison of parameter flow rate, Q between three different meshes and Richardson Extrapolation estimation

Fig. 8 .
Fig. 8. Close-up velocity vector distribution for the current simulation (DNS) at drain time 21s

Fig. 11 .
Fig. 11.Close-up velocity vector distribution for the current simulation (LES) at drain time 21s

Table 1 .
Grid parameter for case A, B and C

Table 2
below displays the results order of accuracy for draining time,   and flow rate,.The results are calculated from three different meshes as shown previously in Table

Table 2 .
Grid Convergence Index (GCI) and order of accuracy for draining time, t drain and flow rate,Q.Index 1, 2 and 3 signify case A, B and C.

Table 3 .
The comparison of draining time of the swirl case between published data by Park and Sohn (experimental, theoretical and numerical), and current studies (DNS, RANS k − ε and RANS k − ω)