Total volume conservation in numerical simulation of liquid sloshing phenomena in partially filled containers

A computational model is developed for incompressible, 2-dimensional, unsteady free surface flows to study the conditions of total volume conservation. Free surface tracking methods used in the numerical simulation of unsteady free surface flows may introduce sources or sinks resulting in changes in total fluid volume in the computational domain. The model is based on finite volume discretization of the Navier-Stokes equations coupling momentum and mass conservation. Free surface position is tracked using VOF method. Possible free surface cell configurations and a solution procedure for continuity are described. Liquid sloshing in a partially filled rectangular containers are simulated. Numerical solutions preserving total volume are presented. Computed free surface profiles are verified by experimental. The Navier-Stokes equations are assumed to hold in the liquid domain and have been solved numerically in developed two-dimensional grids, while the dynamical effects in the ideal gas are disregarded. Also all of the vortexes generated in the flow have been modelled and analyzed numerically. Comparisons of the computed results with exact solutions showed that the method is capable of achieving high accuracy and efficiency with minimal computational effort.DOI: http://dx.doi.org/10.5755/j01.mech.18.4.2326


Introduction
Flow with a free surface includes a wide step of flows in the nature and industry.Among them are the flow through channels and rivers feeding the oceans and dam break flow.Kačeniauskas [1] performed the dam break flow simulation by the finite elements and the pseudoconcentration method.One of the issues related to the flow with the free surface is the "slosh" phenomenon.Liquid oscillation due to forced excitation is called sloshing.Liquid oscillation in a partially filled tank which is associated with various engineering problems.The violent sloshing of liquid creates highly localized impact pressure on the container walls, which may in turn cause structural damages and may even create sufficient moment to destabilize the vessel that carries the liquid container.In order to simulate the fluid flow phenomena with satisfactory resolution in space and time domains, the Navier-Stokes equations are solved numerically using a suitable discretization technique such as finite difference, finite volume, or finite element methods.The finite volume method is widely used in computational fluid dynamics (CFD) applications due to the advantages of integrated equations to enforce conservation of physical quantities in arbitrary control volumes.Free surface simulation became a popular research area with the rapid development of computing power.There are several methods to track the free surface position such as Marker-And-Cell (MAC) method [2], volume of fluid (VOF) [3], and level set methods [4], smoothed particle hydrodynamics (SPH) [5].The MAC method described by Harlow and Welch [2] was the first to track the free surface on a discretized domain.This method uses mass less mark-er particles, which are used to indicate the fluid configuration showing which region is occupied by fluid and which region is empty.The marker particles are moved to new positions using local fluid velocities.Chan and Street developed the SUMMAC method [6], which is a modification of the MAC method, using an interpolation technique for free surface velocity instead of solving the continuity equation on the free surface.Level-set methods, originally introduced by [4], have been applied to a wide variety of immiscible interfacial problems.Kačeniauskas [7] described the development of the mass correction technique for viscous incompressible moving interface flows by interfaces capturing approach.These methods use a level function  with the 0 contour level defining the material interface, much the same as the fractional volume function c in the VOF method, to indicate the shortest distance to the interface.The VOF method [3,8] is a popular interface tracking algorithm that has proven to be a useful and robust tool since its development decades ago.However, these methods are convenient for large domain computations since they are easy to implement, computationally efficient, and require minimum memory storage.Exact volume conservation is not guaranteed in the above methods and it can be difficult to preserve the total fluid volume in the computational domain for long simulation durations.Free surface tracking methods used in the simulation of free surface flows may produce sources or sinks during the computations, causing the total fluid volume to change over time Such observations have been reported by Wang et al. [9].After some numerical tests it was recognized that there may be two sources of error causing variations in total fluid volume in the computational domain.One source of the error is due to inappropriate boundary conditions applied at the far field of a large domain extending to infinity.The second source of error may originate from the free surface tracking algorithms.The motivation behind the present study was to investigate the conditions of satisfying exact volume conservation in a numerical numerical simulation of incompressible free surface flows so that the total fluid volume in the computational domain remains the same throughout the simulation time.The main aim of the study is to develop a computer code for unsteady free surface flows suitable for large spatial domains such as dam reservoirs and sloshing liquid in the container.The numerical methods adopted in the present simulations are based on Hirt and Nichols VOF method [3] coupled with Youngs piecewise linear interface calculation (PLIC) scheme [10], Brack_bills continuum surface force (CSF) model [11], and solved by algebraic multigrid (AMG) solver [12], as well as k-e turbulence model [13], and the pressureimplicit with splitting of operators (PISO) scheme for pressure-velocity coupling [14].

Governing equations
Mass conservation is expressed as volume conservation for incompressible flows considered in this study.The equations of motion for 2-dimensional in compressible flows in a vertical plane are given as 2 1 () where x and z are coordinate axes in horizontal and vertical directions, respectively, x a and z a are ground accelera- tions, V is the velocity vector, p is pressure,  is fluid density.
Ground accelerations are included to represent earthquake excitations or accelerations due to shaking of experimental tanks.The computational domain is assumed to move with the ground.To obtain finite volume formulation of free surface flow, Eqs. ( 1), (2), and (3) are integrated applying Gauss divergence theorem.

Mathematical formulation
where Δt is the time step, Δx and Δz are mesh sizes and Dif and Con represent diffusive and convective fluxes.
Respectively to utilize the advantage of the staggered grid system, convenient control volumes are selected for each equation as shown in Figs. 1 and 2. First order derivatives in diffusive fluxes are discretized using second order polynomial approximation on a variable mesh.Convective fluxes are evaluated by first order upwind (FOU) and by quadratic upstream interpolation for convective kinetics (QUICK) approximation using a 3-point upstream-weighted quadratic interpolation for cell face values.The code can switch between FOU and QUICK, depending on user preferences.Detailed descriptions on temporal and spatial l discretizations can be found in Li and Baldacchino [15].Pressure solution is obtained from the Poisson equation for pressure.The discretized form of the Poisson equation for pressure is obtained by substituting Eqs. ( 7) and ( 8) into (11).
The momentum Eqs. ( 7) and ( 8) and the pressure Poisson Eq. ( 12) are solved by sequential iterations.A computer code in Matlab software developed by the authors to perform the computations.In the numerical solution, location of the free surface must be followed by an appropriate algorithm identifying the computational cells as fluid cell (F) when it is completely filled by fluid, empty cell (E) when there is no fluid, surface cell (S) when the cell is partially filled, and boundary cell (B) on the solid boundaries (Fig. 3).

Computational cells and labelling
When applying any numerical procedure on a computational cell, 4 neighbours, namely, east, west, north, (12) and south contiguous cells, must be identified.Integration and interpolation practices depend on the type of neighbouring cells.A standard solution procedure is applicable for F cells when the neighbours are also composed of F or S cells.Special procedures are required to obtain velocity and pressure when E cells appear as neighbours.A common approach used to calculate the velocity components at the interfaces of E and F cells is extrapolation of the velocity components from the closest available velocities obtained from momentum solutiotailed description of the extrapolation procedure can be found in Armenio [16].
The pressure equation is solved in F cells, the horizontal momentum equation is solved at F-F, F-S, and S-S interfaces, and the vertical momentum equation is solved at F-F and F-S interfaces.Pressure on the free surface in S cells is computed from the free surface stress conditions given by Batchelor [17].A nonzero pressure just on the free surface is evaluated from Then, the pressure at the centre of the surface cell is calculated by linear interpolation (Fig. 4) where A detailed description of free surface boundary conditions can be found in Chen et al [18].A detailed description of free surface boundary conditions can be found in Chen et al., Kleefsman [7], Griebel [19].The momentum equations, pressure equation, and free surface equation are solved by sequential iterations using an explicit procedure.In order to ensure computational stability of the numerical algorithm, a combined stability condition is imposed on the time step based on the convection and diffusion processes (Chan and Street) [20]     were CFL is the Courant-Friedrichs-Lewy number, which is fixed as 0.5 throughout this study, and c is the surface wave celerity.
The pressure Eq. ( 12) is solved by the Point Successive-Over-Relaxation (PSOR) method.

Review of VOF method
The primary works conducted on the VOF method [3] are related to the early 1970.Three methods of tracking the volume were introduced which included Debar method [21], and SLIC method (simple-line-interface) presented by Woodward and Noah [22], and VOF method which became well-known by Donar-acceptor method of Hirt and Nicholas [3].The main development in surface tracing was resulted under the Youngs method as a method entitled Piece Linear Schemes of Youngs [23].Much of the improvements and developments in the work of Young have been made after his work.These versions are known as PLIC methods.In this paper, the PLIC method has been used for constructing and tracing the free surface because of if its higher precision and applicability in complex flows compared to other methods such as donar-acceptor and the Euler-explicit plan and implicit plan.

2. Governing equations
In VOF method, the movement of the interface surface between several immiscible fluids with different density and viscosity is tracked by function of volume of fraction and interface surface is obtained through the following three conditions With respect to the local values of k C , the appro- priate and relevant variables and properties for each control volume inside the domain are determined.The volume of fraction function is obtained by the volume fraction equation where U is the velocity vector of the fluid.The two-phase flows are modeled through the Navier-Stokes equation where g is the gravity vector and P describes pressure.The velocity field of the incompressible flow follows the following limitation In two-phase systems, the properties which exist in the momentum equation are controlled through determining the components of the phases in each cell volume.The average values of density and viscosity are obtained through the following formula [24] ) ( The sum of the fraction volumes , ij C in each cell is equal to the unit.So by substituting the average values of above properties in the Navier-Stokes and Volume fraction Equations the magnitude of the volume fraction , ij C at each cell in the interface were defined and by applying this property in the free surface reconstruction algorithm we can predict the free surface shape.A VOF geometric reconstruction scheme is divided into two parts a construction step and a propagation step.The key part of the reconstruction step is the determination of the orientation of the segment.This is equivalent to the determination of the unit normal vector n to the segment.Then the normal vector  ( ,0) Once the interface is reconstructed, the velocity at the interface is interpolated linearly and the new position of the interface is calculated by the following formula

Numerical analysis
In this paper, we have examined the fluid flow in a reservoir with motive bottom surface in two-dimensional case.The reservoir used in this study is 1  1 meters and the numbers of meshes are 200  200 and 400×400.We have studied different regimes, laminar and turbulent states in a reservoir and then analyzed the stream line contours in different time steps.We have specified the position and location of the vortexes.Then we have compared the results obtained from excited the reservoir in harmonic condition with code and experimental results, which indicate the high precision of the analysis.

Liquid sloshing in an oscillating rectangular
In this state we consider a reservoir with 1m×1m dimensions and with 100  100 and 300  300 grids.Liquid sloshing inside a partially filled rectangular tank is one of the well known test cases for unsteady free surface flows.The tank (Fig. 5) is oscillated along the horizontal axis by a sinusoidal excitation.Horizontal displacement α of the container is given by (2 / ) where A is the amplitude of horizontal displacement and T is the period of oscillation Eq. ( 25) is differentiated twice to obtain the horizontal acceleration.The parameters of oscillation selected are the same as those in the experiments reported by Okamoto and Kawahara [25] (A= 0.93 cm, T =1.183 s, b = m, H = 0.5 m) to compare the free surface computations with the experimental data.This test case was defined such that the frequency of oscillation was equal to the natural frequency of the tank to observe possible resonating free surface oscillations.The boundary conditions defined as stationary asymptotic walls and motive bottom surface having variable velocity and upper surface with pressure outlet condition, and the reservoir is filled in half by liquid which has viscosity  = 0.001 kg/ms and density  =1000 .The free surface profile obtained from computer code by VOF method at t = 1.2 sec compared to experimental measurements given by Okamoto and Kawahara [25] in Fig. 6 and shows high precision of the computer code.

Fig. 5 Definition of geometric parameters in an oscillating
Different grid distributions are used to investigate the mesh size effect on the free surface computations.Grid clustering in the vertical direction is applied around the free surface to reduce possible truncation errors due to partially filled surface cells.Two uniform grids arrangements, 100  100, 300  300, are used to obtain the free surface profiles at t = 1.2 s.Free surface profiles obtained by VOF compared to experimental data of Okamoto and Kawahara [25] are shown in Fig. 6.VOF solutions exhibit small oscillations but volume conservation is satisfied exactly at all time steps.The spatial and temporal agreement between the computations and the experiment clearly indicates the accuracy achieved by computational developments and the adequacy of the grid arrangement (300  300) selected in the solution.Another test for volume conservation was performed by observing the diminishing of surface deformations when the tank oscillation is stopped.The tank is oscillated for the first 12 s and then stopped.It is observed that the steady-state water depth at the end of computations is the same as the initial water depth, indicating no leaks or sources of fluid volume after 1000 seconds of simulation.We consider a cavity as 1×1meters with 300×300 meshes and boundary conditions defined as having stationary asymptotic walls and motive bottom wall with 1 m/s velocity and upper surface with Pressure Outlet conditions, and then fill half of it with a liquid having kinematical viscosity (μ = 0.09 kg/ms) and density (ρ = 450 3 kg/m ) or in other words, the Reynolds number Re = 5000.A time step = 0.1 was applied to solve the unsteady state laminar flow.After passing a time of about 30.01 seconds the flow reaches to the steady state Fig. 10, a. Numerical simulation was applied for different times Figs.8 and 9.In order to achieve this aim a CPU time = 186500 seconds was passed.In the simulation of multiphase flows the VOF method for tracking the surface and the PISO solution control method has been used.For solving the flow in turbulent state, we consider the physics of the problem just like the laminar state except that the liquid used in this state is a liquid with μ = 0.001 and ρ = 998, or in other words, Re = 6  10 .It must be noted that the size number of grids in this state is appropriate for the turbulent state, and the (κ-ε) turbulent model has been used.The momentum equations, the volume of fraction equations, and the disturbing loss and kinetic energy equations have been solved by second-order up wind method.We see that the unsteady state turns to the steady state after passing 127.7 seconds Fig. 11, b  For laminar state at 30.1 sec and for turbulent state at 127.7 sec the flow becomes steady and the free surface profile never changes.

Conclusion
A computational model for simulating 2-dimensional unsteady free surface flows is developed.The model solves Navier-Stokes equations using a finite volume technique.The free surface tracking is accomplished by VOF method to investigate total volume conservativeness of the formulations.The model has been tested in sloshing in an oscillating tank.The following conclusions are drawn from this study: 1.The VOF method produces free surface profiles satisfying mass conservation exactly.The algorithm developed here can deal with steep surface slopes, such as dambreak, without computational problems.Accuracy in predicting local water depths may deteriorate due to extrapolations in the free surface computation when a coarse grid is used.However, total volume conservation is always satisfied.
2. Total volume conserving property and computational accuracy are independent of the dimensions of the computational domain.This feature of the model makes it suitable for computations in large spatial domains.
3. Use of FOU or QUICK for convection terms has no effect on volume conservation in the VOF solution.By comparing the answers of the steady states we find that in laminar state (Fig. 10, a) the resulted circulation is closer to the right wall and the flow in the bottom right corner is sharper.The resulted circulation in the laminar flow is more symmetric, this is resulted from the velocity horizontal component distribution graphs in Fig. 12.
4. As it was expected, the laminar flow reaches the steady state sooner (30.1 sec) (Fig. 10, a), while the turbulent flow, at the same time (30.1 sec), is in the middle part.By solving the laminar flow in unsteady state and with time step = 0.1, after passing a time of about 30.01 seconds the flow reaches the steady state and no change is observed in the circulation.In order to achieve to this purpose, a time equal to CPU time = 186500 s should be passed.Turbulence state after a time of 127.7 seconds the flow reaches the steady state and no change is observed in the circulation it obvious that the problem in laminar state takes shorter time to reach the steady state condition than turbulent state.

Fig. 4
Fig. 4 Pressure interpolation in a surface cell

C 3 . 3 .
uniquely determine a straight line.Once the interface has been reconstructed, its motion by the underlying flow field must be modelled by a suitable algorithm.The fluid advection algorithm During the advection step, the volume of fraction ij C is truncated by the formula at the (n+1) time step ,, 1,

CC
field is obtained according to the local velocity field, and fluxes DC at each cell are determined by algebraic or geometric approaches.Here, the simplest operator split advection (geometric) algorithm is used as with a y sweep.

Fig. 6
Fig. 6 Comparison of the free surface profiles obtained by VOF method with experimental [25] in different grid arrangements (100  100 and 300  300) at time = 1.2 Computations are continued with no excitation until the kinetic energy of the oscillating fluid volume is dissipated totally and the free surface becomes nearly horizontal.Water level at the left wall is shown as a function of simulation time in Fig. 7.During the first 16 s, the surface waves are amplified, reaching a maximum wave height of approximately 0.4 m, and then start to decrease when the tank oscillations are stopped.It takes about 300 s for the surface waves to be reduced to negligible amplitude.It is observed that the steady-state water depth at the end of computations is the same as the initial water depth, indicating no leaks or sources of fluid volume after 1000 seconds of simulation.

Fig. 7
Fig. 7 Water level at the left wall of the oscillating tank as a function of time 4.2.Simulating the laminar flow in the rectangular reservoir with motive bottom surface

4. 3 .
Simulating the turbulent flow in the liquid reservoir with motive bottom surface

Fig. 14 Fig. 15 Free
Fig. 14 Illustrates the skin friction coefficient on the reservoir moving wall in term of x-position of reservoir

A
. Kh. Poorfar, M. M. Shahmardan, M. Sedighi TOTAL VOLUME CONSERVATION IN NUMERICAL SIMULATION OF LIQUID SLOSHING PHENOMENA IN PARTIALLY FILLED CONTAINERS S u m m a r y