Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 46, 2, pp. 325-345, Warsaw 2008 COMPARISON OF CICSAM AND HRIC HIGH-RESOLUTION SCHEMES FOR INTERFACE CAPTURING Tomasz Wacławczyk Tadeusz Koronowicz The Szewalski Institute for Fluid Flow Machinery, Polish Academy of Sciences, Gdańsk, Poland e-mail: twacl@imp.gda.pl; ttk@imp.gda.pl Thepaper concerns themodelling of aflowwithhighdensity andviscosi- ty ratios using theVolume of Fluidmethod (VOF). Two high resolution schemes are used for discretization of the nonlinear convective term in the equation for transport of the volume fraction. Properties of the sche- mes are compared by taking into consideration their ability to capture the interface subjected to strong deformation. We show that diffusive properties of the high-resolution schemes affect the obtained solution. Verification of the schemes is performed in two test cases: rotation of a solid body with a slot and sloshing of water in a rectangular tank. Key words: free surface flow, Volume of Fluid method, high-resolution schemes 1. Introduction Numerical modelling of multiphase flows is a developing branch of the Com- putational Fluid Dynamics (CFD). In this paper, incompressible viscous two- phasefluidswithawell defined interfacewill be considered.Themaindifficulty here is tracking of the interface position during integration in time since this information is needed for proper arraignment of boundary conditions. To cap- ture the interface position, the Volume of Fluid (VOF) method is employed (Hirt andNicholls, 1981). In the framework of theVOFmethod, the evolution of the free surface in time is described by an additional transport equation for the volume fraction of the background fluid. The value of the volume fraction indicates presence or absence of the background fluid in the control volume. Control volumes that are partially filledwith the background fluid contain the interface.When coupledwith theNavier-Stokes equations, the volume fraction 326 T. Wacławczyk, T. Koronowicz is treated as an active scalar (it has influence on the velocity field). Themain numerical difficulties connected with discretization of the transport equation for the volume fraction are: keeping constant width of the interface, i.e. avo- iding artificial diffusion of the step interface profile and assuring a monotonic change of the variables. This last condition is also known as the boundedness criterion (Ferziger and Perič, 2002; Hirt andNicholls, 1981; Versteeg andMa- lalasekera, 1998). In order to overcome the aforementioned problems different methods were proposed. For instance, in Hirt and Nicholls (1981) a Donor- Acceptor Scheme (DAS), based on the availability criterion, was introduced. Problems that arose when using this scheme, see Noh andWoodward (1976), Youngs (1982), provoked other proposals that follow the idea of geometric interface reconstruction; examples are SLIC (Simple Line Interface Calcula- tion) method (Noh and Woodward, 1976), PLIC (Piecewise Linear Interface Construction) method (Youngs, 1982) or more recent methods that use the least-square procedure (Pillod and Puckett, 1991, 2004) or splines (Lopez et al., 2004). Methods that employ these ideas give good approximation of the shape of the interface and they allow for proper calculation of the fluxes thro- ugh faces of the control volumes. However, their application is often restricted to structured grids with simple shapes of the control volumes. Moreover, sin- ce estimation of a spatial orientation of the interface from the distribution of the volume fraction needs a substantial number of numerical operations, interface reconstruction methods increase the computational effort (Zaleski, 2002). In this work, to discretize the convective term in the equation for trans- port of the volume fraction, high-resolution schemes: CICSAM (Compressive Interface Capturing Scheme for Arbitrary Meshes) (Ubbink and Issa, 1999) and HRIC (High Resolution Interface Capturing scheme) (Muzaferija et al., 1998), are used. Unlike geometric interface reconstruction methods, the high- resolution schemes donot introduce geometrical representation of the interface but try to satisfy the aforementioned conditions by properly chosen discreti- zation scheme. Tomodel dynamics of a set of immiscible fluids, a solution to the equation of transport of the volume fraction needs to be coupled with the solver of the Navier-Stokes equations. In this paper, an in-house, finite volume second order accurate code is used where the SIMPLE procedure is employed for pressure- velocity coupling, see Ferziger and Perič (2002), Versteeg and Malalasekera (1998). The developed algorithm is used to model the flow in a system of immiscible fluids and for comparison of the aforementioned high-resolution schemes. Comparison of CICSAM and HRIC high-resolution schemes... 327 2. One-fluid formulation for multiphase flow In this paper, one-fluid formulation for a multiphase flow is used. According to thismodel, dynamics of viscous, incompressible and immiscible fluids is go- verned by singleNavier-Stokes equations (2.1)1 and continuity equation (2.1)2 ∂(ρui) ∂t + ∂(ρujui) ∂xj =− ∂p ∂xi + ∂ ∂xj [ µ (∂ui ∂xj + ∂uj ∂xi )] +ρgi+σκni (2.1) ∂ui ∂xi =0 where ui is the i-th velocity component, ρ denotes density, µ is the dynamic viscosity, p is the pressure, gi is the i-th component of the gravitational ac- celeration. The additional term that appears in Navier-Stokes equation (2.1)1 is inherited from one-fluid formulation. This term represents surface tension acting locally at the interface between fluids. In the case of constant surface tension coefficient σ, the force resulting from surface tension acts in the di- rection normal to the interface n = ∇φ/|∇φ|, where φ denotes the volume fraction of the background fluid. One needs to notice that ∇φ has a non-zero value only at the interface, indicating a local character of the surface tension term. In this paper, the term that models surface tension effects will be ne- glected, since inertial forces are dominant in the case of water sloshing in a rectangular tank. To model dynamics of a set of immiscible fluids, one needs to add con- stitutive relations that convolve material properties of fluids with the volume fraction of the background fluid (here water) in the control volume. Since in this paper the set of two fluids is considered, the constitutive relations write ρ=(1−φ)ρ1+φρ2 µ=(1−φ)µ1+φµ2 (2.2) where ρj, µj represent, respectively, the density and dynamic viscosity of the j-th fluid. Substituting constitutive relation (2.2)1 to the conservative formof the continuity equation and taking into account the assumption about incompressibility of the fluids, one can derive equation for transport of the volume fraction ∂φ ∂t +uj ∂φ ∂xj =0 (2.3) To model a multiphase flow and capture the interface position, one needs to solve Eqs. (2.1)-(2.3). It should be noticed that both fluids are treated as single continuous since they share the same velocity in each control-volume. 328 T. Wacławczyk, T. Koronowicz This simplification allows for relatively easy numerical solution of the problem alsowhenmore than two fluids are considered.On the other hand, sharing the same velocity is a disadvantage of this formulation, because one is not able to model the slip phenomenon. 3. Numerical methods for solution of conservation equations To solve numerically the set of equations (2.1)-(2.3) it is necessary to discretize it in order to obtain a system of algebraic equations. In the present work, an in-house, three-dimensional, finite volume, second order accurate code was used to model the dynamics of a multiphase flow. Next, a brief description of the numerical methods employed for discretization of the Navier-Stokes equations will be given. Two high-resolution schemes used for discretization of the convective term in Eq. (2.3) will be presented in more details. 3.1. Discretization of the Navier-Stokes equation The first step to solve Navier-Stokes equations (2.1)1 is to write them in a conservative, integral form (finite volume method) and to apply the Gauss theorem to convective and diffusive terms (Ferziger and Perič, 2002; Versteeg andMalalasekera, 1998). In the case of the second-order accuratemethod, the assumption about a constant value of the variables, e.g. velocity component, at the cell face centres of the control volumes, leads to the conservative formu- lation. One needs to remember that in a collocated-mesh arrangement, used here, the solution is searched in centres of the control volumes. Therefore, va- lues of variables at the cell face centres have to be connected with their cell centre neighbours through interpolation. In this work, both convective and diffusive terms (fluxes through the faces of the control volumes) in Eq. (2.1)1 are discretized using a deferred correction approach (Ferziger andPerič, 2002; Versteeg and Malalasekera, 1998). The pressure gradient and mass force are discretized in a conservative and non-conservative manner, respectively. The second step is integration in time. It is assumed that the surface and volume force terms are constant over one time step. However, the component of velocity in the nonlinear convective term is assumed to be taken from the new time level. Integration of the local derivative over one time step, together with application of the first order implicit Euler scheme, leads to an impli- cit, first order accurate algorithm. Additionally, one needs to correct velocity, which does not satisfy continuity equation (2.1)2 after segregated solution to Comparison of CICSAM and HRIC high-resolution schemes... 329 the momentum equations. In our code, the SIMPLE algorithm is used for pressure-velocity coupling (Ferziger and Perič, 2002; Versteeg and Malalase- kera, 1998). In this method, the Poisson equation for pressure correction is obtained. Moreover, it can be shown that the velocity correction is propor- tional to the gradient of the pressure correction. Both velocity and pressure corrections are used to obtain a divergence free velocity field as well as cor- rected pressure field. In the case of unsteady flows, it is usually necessary to perform more than one outer iteration per time step to accurately satisfy continuity equation (2.1)2. Boundaryconditions forEqs. (2.1)-(2.3) and thePoisson equationobtained in the SIMPLEprocedure are test case dependent. In thiswork, the symmetry plane and no-slip wall boundary conditions were used, see Ferziger and Perič (2002), Versteeg and Malalasekera (1998) for Eq. (2.1)1 and Eq. (2.3). Since boundary conditions for the Navier-Stokes equation determine velocity at the borders of the computational domain, the Neumann boundary condition is applied for the Poisson equation at each side of the domain. 3.2. Discretization of the equation for transport of the volume fraction, high-resolution schemes Proper discretization of Eq. (2.3) is crucial for simulation of a multiphase flow. It can be shown (Wacławczyk and Koronowicz, 2005), that numerical schemes, commonly used for discretization of the convection term, introduce numerical diffusion or numerical dispersion phenomena. For this reason, some additional techniques are needed, i.e., high-resolution schemes. Themain task for a high-resolution scheme is to discretize the convective term in Eq. (2.3) in a way that will prevent artificial smearing of the step interface profile due to numerical diffusion.Moreover, it should preservemonotonic distribution of the variable, i.e., it should satisfy the boundedness criterion. Discretization of Eq.(2.3) in time is performed with the Crank-Nicolson method, shown to be consistent with an unsplit time integration procedure (Ubbink and Issa, 1999). Since in this paper two high resolution schemesCICSAMandHRICare compared, abrief description of schemeswill be given, emphasising similarities and differences between both formulations. 3.2.1. Normalised Variable Diagram TheNormalised Variable Diagram (NVD) (Leonard, 1991) provides the foun- dation to CICSAM and HRIC high-resolution schemes. Let us consider the control volumes denoted by U, D, A for upwind, donor and acceptor faces respectively; f denotes the face of the control volume, see Fig.1a. Since the 330 T. Wacławczyk, T. Koronowicz Fig. 1. (a) CBC criterion on the three control volumes U upwind, D donor, A acceptor and face f, (b) Normalised Variable DiagramNVD, UDS upwind scheme, CDS central scheme, DDS downwind scheme value of the variable at the face of the control volume should satisfy the Co- nvective Boundedness Criterion (CBC), see Leonard (1991), Ubbink and Issa (1999), Versteeg and Malalasekera (1998), one can write: φD ¬ φf ¬ φA. Using this constraint and information about the value of the variable in the upwind control volume φU, normalised variables are introduced φ̃f = φf −φU φA−φU φ̃D = φD−φU φA−φU (3.1) where φD, φA, and φU are the volume fractions in the donor, acceptor and upwind cells, respectively. It could be noticed that when φ̃f = φ̃D, which means φf = φD, the cell face value is calculated with the first order accu- rate Upwind Scheme (UDS). On the other hand, when φ̃f = 1 which entails φf = φA, the cell face value is calculated with the first order accurate Down- wind Scheme (DDS).When φ̃D < 0 or φ̃D > 1, the only scheme that satisfies CBC is UDS. Thus, when 0 < φ̃D < 1, only schemes that satisfy Eq. (3.2), i.e. lie in the shaded region on Fig.1b, preserve CBC. φ̃f = { φ̃D for φ̃D < 0, φ̃D > 1 φ̃D ¬ φ̃f ¬ 1 for 0¬ φ̃D ¬ 1 (3.2) One needs to notice that when using Eqs. (3.1) it is possible to find the value of the volume fraction at the control volume face φf φf =(1− β̃f)φD + β̃fφA β̃f = φ̃f − φ̃D 1− φ̃D (3.3) which is used for calculation of the volume fraction flux in Eq. (2.3). However, to calculate φf, one needs to find φ̃f first. Comparison of CICSAM and HRIC high-resolution schemes... 331 3.2.2. Compressive Interface Capturing Scheme for Arbitrary Meshes – CICSAM In the case of the CICSAM scheme, an additional assumption about the de- pendence of the region where the CBC is satisfied on the CFL condition is used (Leonard, 1991). The local value of the Courant number, defined at the face Sf of the control volume Cf = ufSf∆t/V together with the CBC gives the following formula φ̃D ¬ φ̃f ¬min{1, φ̃D/Cf}, which can be plotted at a normalised variable diagram, see Fig.2. One needs to notice that for explicit schemes, if a local value of the Courant number is equal to Cf =1, only the UDS satisfies the CBC criterion. On the other hand, the smaller the value of the local Cf the lager the domain where the CBC criterion is satisfied. Fig. 2. Dependence of the CBC region on the local CFL condition Coupling of the donor-acceptor scheme (DAS) (Hirt and Nicholls, 1981), with NVD formulation dependent on the CFL condition, results in the first component of CICSAMknown as the HYPER-C scheme φ̃fCBC =    φ̃f = φ̃D for φ̃D < 0, φ̃D > 1 min { 1, φ̃D Cf } for 0¬ φ̃D ¬ 1 (3.4) The above HYPER-C scheme satisfies the CBC criterion and is shown to be compressive, which means that it changes any gradient of φ to a step profile due to DDS employed. However, the compressive character of theHYPER-C is not always desira- ble.When the interface is tangential to the flow direction, it is shown that the aforementioned scheme tends to artificially deform its shape.For this reason, it is found to be necessary to switch between scheme Eq. (3.4) and an other less 332 T. Wacławczyk, T. Koronowicz compressive formulation. In the case of the CICSAM, the Ultimate-Quickest (UQ) scheme is employed (Leonard, 1991) φ̃fUQ =    φ̃D for φ̃D < 0, φ̃D > 1 min {8Cfφ̃D+(1−Cf)(6φ̃D +3) 8 , φ̃fCBC } for 0¬ φ̃D ¬ 1 (3.5) To switch smoothly between both schemes, a linear blending is used with the blending factor 0¬ γf ¬ 1. The value of the γf depends on the angle θf, see Eq. (3.6)2 between the unit vector normal to the interface n = ∇φD/|∇φD| and the unit vector parallel to the line between centres of the donor D and acceptor A cells d = −−→ DA/| −−→ DA|, see Fig.3a. When the interface position is normal to the direction of the flow, γf =1 and Eq. (3.4) is used. In the case of tangential orientation of the interface, γf =0 and Eq. (3.5) is employed φ̃f = γfφ̃fCBC +(1−γf)φ̃fUQ 0¬ γf ¬ 1 θf =arccos |dn| (3.6) γf =min {1+cos2θf 2 ,1 } It should be noticed that the above derivation of the high-resolution CICSAM scheme was carried out only in one-dimension. To extend above formulation to multiply dimensions, the cell Courant number is introduced CD = ∑n f=1max{Cf,0}, where n denotes the number of control volume fa- ces. Fig. 3. (a) Definition of the vector normal to the interface n=∇φ/|∇φ| and vector parallel to the line between the centres of D, A control volumes d= −−→ DA/| −−→ DA|, (b) continuity of the HRIC scheme in the time domain Comparison of CICSAM and HRIC high-resolution schemes... 333 3.3. High Resolution Interface Capturing (HRIC) scheme To simplify the above procedure and to get rid of the explicit dependence on theCFLcondition, theHRIC (Muzaferija et al., 1998) schemewas introduced. As itwasmentioned, this scheme also relies on theNVDandnormalised varia- bles. Application of the HRIC scheme can be divided into three steps. Firstly, the normalised cell face value will be estimated from a scheme that conti- nuously connects the upwind and downwind schemes on the NVD diagram, cf. Fig.1b φ̃f =    φ̃D for φ̃D < 0, φ̃D > 1 2φ̃D for 0¬ φ̃D < 1 2 1 for 1 2 ¬ φ̃D ¬ 1 (3.7) Secondly, since DDS can cause alignment of the interface with the mesh and its artificial deformation (as in the case of the HYPER-C scheme), one needs an other scheme that satisfies the CBC. In the case of the HRIC, first order UDS is used as themost straightforward choice. Again, the blending factor γf connected with the angle θf is introduced, see Eq. (3.8)2, to switch smoothly between the schemes φ̃∗f = γfφ̃f + φ̃D(1−γf) γf = √ |cosθ| (3.8) The blending of the UDS and DDS schemes is dynamic and takes into acco- unt the local distribution of the volume fraction. In the case when the CFL condition is not satisfied, the dynamic character of this scheme can cause sta- bility problems. Therefore, φ̃∗f is corrected with respect to the local Courant number Cf, see Eq. (3.9). The goal of this correction is to force continuous switching between the schemes also in the time domain, see Fig.3b φ̃∗∗f =    φ̃∗f for Cf < 0.3 φ̃D for Cf > 0.7 φ̃D+(φ̃ ∗ f − φ̃D) 0.7−Cf 0.7−0.3 for 0.3¬Cf ¬ 0.7 (3.9) When using this scheme inmultiple dimensions, the local Courant number Cf is again replaced by its cell definition CD. One can notice that the main difference between the CICSAMand HRIC is the order of accuracy of the component schemes and the differently included CFLcondition.Next, two testswill beperformedto asses howthosedifferences influence the results obtainedwith the two presented high-resolution schemes. 334 T. Wacławczyk, T. Koronowicz 4. Advection test case The aim of the advection test case is to check how the high-resolution schemes preserve the initial shape of the interface in a given velocity field. Since it is known that when convection is the only transport phenomenon, the distance between two points advected in the uniform flow should remain constant. InWacławczyk andKoronowicz (2006) a circular spot advection test case was carried on to show that both the CICSAM and HRIC schemes are first order accurate. Themain observation concerning properties of the considered high-resolution schemes was better, i.e. sharper and less diffusive, interface reconstructionof theadvected shapeby theCICSAMscheme.However, in that previously performed test case, we analysed only a two dimensional contour of the volume fractiondistribution φ=0.5,whichdidnot givewhole information about reconstruction of the advected shape. Fig. 4. Computational domain for rotation of a solid body with a slot Another advection test case known in literature (Ubbink and Issa, 1999) is rotation of a solid body with a slot placed in a rotational velocity field. This test case allows one to assess how the shape of the interface is reconstructed in amore complicated velocity field, comparing thewhole range of the volume fraction distribution φ= [0,1]. Here, a circular solid bodywith a slot is revolved with a constant angular velocity ω that is chosen to satisfy CD ≈ 0.15 at the body edges. The final solutionwas obtained after 1360 time steps,with ∆t=0.005s on 64×64 grid. Figure 4 depicts the computational domainwith the initial distribution of the volume fraction and used boundary conditions. Inside of the two dimensional solid body the volume fraction is set to one, outside the body is set to zero indicating the presence of the background fluid, see Fig.5a. Comparison of CICSAM and HRIC high-resolution schemes... 335 Fig. 5. Final shape of the solid body after one revolution in rotational velocity field, (a) initial condition, (b) solution with CICSAM, (c) solution with HRIC Comparing the results obtained with the CICSAM and HRIC schemes, one can notice that, again, the first one preserves the shape of the interface better (volume fraction shape is closer to the initial condition) than the latter, compare Fig.5a and Fig.5b, what is visible in the whole range of volume fraction distribution. The interface calculated with HRIC is more smeared, non-zero volume fraction values are visible in the slot after one revolution, see Fig.5c. Conclusion from this test case is that CICSAMworks better, which is in agreement with the results obtained in the first circular spot advection for small CD values (see Wacławczyk and Koronowicz, 2006). This behaviour of both schemes is caused by the existence of numerical diffusion introduced by the differently included first order upwind scheme.Notice the explicit presence of theUDS inHRIC formulation during blendingwith respect to the interface position, see Eq. (3.8)1, for more details see alsoWacławczyk et al. (2007). 5. Sloshing of water in an oscillating tank To compare the properties of the schemes in more realistic situations, we carried out simulation of water sloshing in a quadrangle container. This issue is important for vehicle safety, e.g. ship or rig, when the cargo transported is liquid. Themain goal here is to find resonance periods for swayingmotions of the tank because the sloshing liquid can damage it. In the next sections of the paper, experimental data from Hinatsu et al. (2001) are compared with the results of our numerical simulation. The experimental setup for this test case consists of a tank filled with a liquid in two different levels. At the south, east and north walls of the tank, pressure transducers are installed, see Fig.6, to record the pressure history during motion of the tank. The parameters governing behaviour of the fluid inside the tank are: filling level of the tank Φ [%], period T [s] and ampli- 336 T. Wacławczyk, T. Koronowicz Fig. 6. Experimental setup for measurements of pressure during sloshing of water in a quadrangle container, P1,P2, . . . ,P9 depict position of transducers for pressure measurement tude A [m] of the sway motion. Several cases were investigated in Hinatsu et al. (2001) with two different filling levels; details are listed in Table 1. In this paper, we will consider the filling level Φ = 20% and two different time periods: T =1.74s when resonance of the water/tank system is observed and T =1.94swithout resonance, to compare dynamic behaviour of the fluids, i.e. water and air, in both cases. Themeasurement data consists of time history of the pressure, i.e., the signal registered by the set of transducers during motion of the tank. Additionally, the position and acceleration of the tank were measured as well, see Fig.7. Table 1. Sloshing conditions Φ [%] A [m] T [s] Observed events Active transducers 20 0.06 1.74 resonance P1,P2,P3 20 0.06 1.94 – P1,P2,P3 60 0.015 1.40 - P4,P6,P8,P9 60 0.015 1.47 resonance P4,P6,P8,P9 In the present study, the position and acceleration provided by experimen- tal data are used to tune sway motion parameters. For the sake of the higher accuracy amplitude, the period and phase were changed: A1 = 0.0605m, T1 = 1.74s, φ1 = 0.41 in the first case with resonance and A2 = 0.064m, T2 =1.74s, φ2 =0 in the second case without resonance, compare Table 1. Comparison of CICSAM and HRIC high-resolution schemes... 337 Fig. 7. Measured (solid line) and theoretical (dashed line) history of swayingmotion and acceleration of the tank for two periods of swaying: A=0.06, T =1.74 (top), A=0.06, T =1.94 (bottom), in both cases the tank filling was Φ=20% 5.1. Numerical simulation – sway of the tank The numerical simulation was carried out using a two dimensional compu- tational domain in a real scale, discretized with a uniform orthogonal grid consisting of 256×128 control volumes. At all borders of the domain, no-slip wall boundary condition was imposed. The time step size used during simu- lation was set to ∆t = 1−4 s to make the maximal cell Courant number CD smaller than 0.2. During both simulations, the whole set of conservation equ- ations Eqs. (2.1)-(2.3) was solved, using the following material properties for air: ρ1 = 1.205kg/m 3, µ1 = 1.8 −5kg/ms and for water: ρ2 = 998kg/m 3, µ2 =1 −3kg/ms.Toobtain the effect of acceleration acting on the tankandflu- id, an external force Fe was added to the right hand side of the x-momentum equation. Sincemotion of the tank in time is describedbyaharmonic function, the source term that models the external force Fe can be written Fe =−ẍρV =Aω 2cos(ωt+φ) (5.1) 338 T. Wacławczyk, T. Koronowicz where A,ω=2π/T ,φ are swaymotion parameters: amplitude, frequency and phase, respectively. The external force is addedwith the opposite sign because of the reference frame associated with the fluid. 5.2. Discussion of the results The results obtained during the numerical simulation are presented in Figs. 8, 9, 10. Let us follow thewater behaviour inmore details. In the first case when T =1.74s, the orientation of the air/water interface is directly comparedwith the photographs fromHinatsu et al. (2001), see Fig.8. Fig. 8. Comparison of the water/air interface position from experiment (rightmost column) with numerical results obtained using CICSAM (leftmost column) and HRIC (middle column) at different times t: A1 =0.0605m, T1 =1.74s, φ1 =0.41 (first period of sloshing) Comparison of CICSAM and HRIC high-resolution schemes... 339 Fig. 9. Pressure history obtained from simulation with CICSAM (left column) and HRIC (right column) schemes compared with experimental data (dashed line) recorded by transducers P1, P2, P3,A1 =0.0605m, T1 =1.74s, φ1 =0.41 At thebeginningof the simulation, due to the initialmovement of the tank, water accelerates. This can be observed as a smooth variation of the pressure about t=1.5s, see Figs.9, 10. In the period of time t from 1.5s to approxi- mately 2.8s, in both considered cases, thewater wave travels from thewest to the east wall of the tank. The firstmaximum of the pressure is observed close to t=2.9s.At this time, the travelling wave hits thewestwall which causes a rapid pressure increase that is recorded by all active transducers P1,P2,P3, see Figs.8c, 9, 10. Afterwards, splashed water climbs up along the wall and then collapses due to the gravitational force, however, not immediately. There 340 T. Wacławczyk, T. Koronowicz is a time interval from 2.9s till about 3.3s, which is needed for collapsing and forming of the water wave that travels back to the west wall of the tank, see Fig.8d. Fig. 10. Pressure history obtained from simulation with CICSAM (left column) and HRIC (right column) schemes compared with experimental data (dashed line) recorded by transducers P1, P2, P3,A2 =0.064m, T2 =1.94s, φ2 =0 This event can be observed in Figs.9, 10 as a platou between the pressure peak (i.e. wave impact) and the region where the pressure starts to decrease. Observed pressure drop is caused by the wave travelling back to the west wall of the tank (pressure transducers aremounted on the east side,Fig.6). Finally, the whole phenomenon repeats, following the described scenario. The pressure history obtained during simulations compares well with the experimental evidencewhen one considers the shape of the solution and dura- Comparison of CICSAM and HRIC high-resolution schemes... 341 tion of characteristic events in this flow, e.g., prediction of the wave collision with the tank wall or time interval of the water column collapsing. However, the maximal values of the pressure during the collision of the travelling wave with the wall are overestimated or underestimated when compared with the experiment. This effect is particularly evident in the case of pressure transdu- cer P3 and interestinglymore significantwhen the resonance occurs, compare the pressure recorded by probe P3 in Figs.9, 10. This transducer records the highest variations of the pressure signal because it is placed directly in front of the hitting wave. The pressure variation reconstructed at the positions of transducers P1, P2 follow more accurately the experimental data, especially in the case without resonance, see Fig.10. The pressure maxima that result from the collision of the approaching wave with thewest wall are larger in the first case T =1.74s,which indicates thepresence of resonance.This is thema- in difference between the no-resonance and resonance frequency of the swaying tank, compare Figs.9 and 10, since in the case of resonance, a larger amount of energy is transported from the moving tank to the air/water system. This last fact has also influence on the flow, and hence the accuracy of numerical reconstruction. 5.3. Influence of properties of the high-resolution scheme on results Aqualitative comparison of the air/water interface position shows that simu- lations carried out with the CICSAM and HRIC schemes recover its shape accurately during the first period of oscillations, see Fig.8. Themain features of the flow are reconstructed in the sameway by both schemes, but the inter- face representation with the CICSAM scheme is sharper despite equal mesh resolution in both cases. The difference in the interface representation can be noticed while the wave collides with the east wall of the tank, see Fig.8c, and the water column collapses, see Fig.8d. When the fluid inside the tank experiences violent behaviour due to mo- tion with the resonance frequency, the air/water system obtains more energy from the external force, which considerably influences the behaviour of the flow.When the resonance occurs, a mixture of the fluid inside the tank beco- mes dispersed and heterogeneous (droplets of water mixed with air). For this reason, as can be noticed in Fig.9, every next sloshing period is reconstructed less accurately by the numerical simulation. This can be compared with the no-resonance case, see Fig.10, where there is no visible shift in time between thenumerical simulations (usingboth schemes) anddata fromthe experiment. The influence of schemes properties on the solution are compared closely in Fig.11. Pressure behaviour during the first period of sloshing in both the reso- 342 T. Wacławczyk, T. Koronowicz nance and no-resonance case is predicted reasonably by both high-resolution schemes (compare the top and bottom rows in the case of the first period of sloshing registered byprobes P1 or P3). Interesting conclusions can bedrawn when one compares the last period of sloshing in the resonance (top row) and no-resonance case (bottom row). In the resonance case, it is visible that lar- ger numerical diffusivity of the HRIC scheme introduces larger shift in time, moreover, the pressure variation becomes smoother. This is connected with a less accurate representation of the interface in the case of a heterogeneous multi-dispersedflow.TheCICSAMbetter represents the shape of the pressure signal but it also introduces smearing and, thus, a time shift in the pressure re- construction. Unlike in the resonance case, when the flow remains calmer and the number of interface reconnections is limited, both schemes reconstruct the pressure history accurately in the first and last period of sloshing, see Fig.11 (bottom row). Fig. 11. Comparison of pressure histories from probes P1, P3 (first and last period of sloshing) in the case with resonance T =1.74s (top) and without resonance T =1.94s (bottom) during simulations with the CICSAM andHRIC schemes It can be also noticed that between succeeding wave impacts, when the wave travels in the direction of the west wall, overshoots and undershoots in the pressure history obtained from the numerical simulation are visible. Comparison of CICSAM and HRIC high-resolution schemes... 343 For instance, in Fig.9 for time t = 5s till t = 6s. This effect is connected with the choking phenomenon that appears, when movement of the water column is so rapid that the air is entrapped in the water near the south-east corner of the tank, which results in a series of oscillations. This effect was not recorded during the experiment, however, other numerical simulations, e.g. Nielsen (2003), also show this deviation from the experimental evidence. 6. Conclusions In the present paper, the VOF method with high-resolution scheme metho- dology was successfully applied and tested. Two high-resolution schemes the CICSAM and HRIC were employed to discretize the convective term in the equation for transport of the volume fraction, Eq. (2.3), in order to capture the interface position without introduction of numerical diffusion, keeping values of the variables bounded. Results of computations were compared with the experimental evidence, showing good quantitative and qualitative agreement. The main features of the flow in a rectangular tank are properly recovered, see Fig.8, by both high- resolution schemes. However, one needs to notice that when the air/water system becomes dispersed, both the CICSAM and HRIC introduce an artifi- cial time shift in the pressure history, see Fig.9. This artificial effect is smaller in theCICSAMscheme.Moreover, the results obtainedwith theCICSAMare less affected by diffusivity and preserve the shape of the pressure signal mo- re accurately. It also important to notice that in the case without resonance, see Fig.10, the results obtained with both schemes follow experimental data withoutavisibledecrease of accuracy in time.Not enoughaccurate reconstruc- tion of the dispersed flows was partially caused by the fact of neglecting the three dimensional effects and usage of the first order upwind scheme in both formulations. Therefore, themain disadvantage of the high-resolution schemes presented here is their low order of accuracy, and thus the furtherwork should be connected with developing of higher order high-resolution schemes. References 1. Ferziger J.H., PeričM., 2002,Computational Methods for Fluid Dynamics, Springer, Berlin 344 T. Wacławczyk, T. Koronowicz 2. Hinatsu M., Tsukada Y., Fukasawa R., Tnanaka Y., 2001, Two-Phase flows for joint research,Proceedings of SRI-TUHHminiWorkshop onNumerical Simulation of Two-Phase Flows, edited byM.Hinats, 12-19,NationalMaritime Research Institute Japan 3. Hirt C.W., Nicholls B.D., 1981, Volume of fluid method for dynamics of free boundaries, J. Comput. Phys., 39, 201-221 4. Leonard B.P., 1991, The ULTIMATE conservative difference scheme ap- plied to unsteady one-dimensional advection, Comput. Methods Appl. Mech. and Eng., 88, 17-74 5. Lopez J.,Hernandez J.,GourezP.,FauraF., 2004,Avolumeoffluidme- thod based onmultidimensional advection and spline interface reconstruction, J. Comput. Phys., 195, 718-742 6. Muzaferija S., PericM., SamesP., SchelinT., 1998,A two-fluidNavier- Stokes solver to simulate water entry,Proc. Twenty-Second Symposium on Na- val Hydrodynamics 7. NielsenK.B., 2003,Numerical prediction of greenwater loads on ships,Tech- nical University of Denmark, PhDThesis 8. Noh W.F., Woodward P.R., 1976, SLICK (Simple Line Interface Calcula- tion), Lecture Notes in Physics, 59, 330-360 9. PillodJ.E., PuckettE.G., 1991,AVolume ofFluid interface tracking algo- rithmwith applications to computing shockwave refraction, 4-th International Symposium on Computational Fluid Dynamics 10. Pillod J.E., Puckett E.G., 2004, Second order-accurate volume-of-fluid al- gorithms for trackingmaterial interfaces, J. Comput. Phys, 199, 718-742 11. Ubbink O., Issa R.I., 1999, Method for capturing sharp fluid interfaces on arbitrarymeshes, J. Comput. Phys., 153, 26-50 12. Versteeg, H.K., Malalasekera W., 1998, An Introduction to Computa- tional Fluid Dynamics, Longman, London 13. Wacławczyk T., Gemici Ö.C., Schäfer M., 2007, Novel high-resolution schemes for interface capturing inmultiphaseflows,Proceedings of International Conference on Multiphase Flows, Leipzig, 82 14. Wacławczyk T., Koronowicz T., 2005,Modeling of the flow in systems of immiscible fluids usingVolume of Fluidmethodwith CICSAM scheme,Turbu- lence, 8, 267-276 15. Wacławczyk T., Koronowicz T., 2006, Modelling of the free surface flow with high-resoluion schemes,Chemical and Process Engineering, 27, 783-802 16. Youngs D.L., 1982, Time-dependent multi-material flow with large fluid di- stortion,Numerical Methods for Fluid Dynamics, 1982, 273-285 Comparison of CICSAM and HRIC high-resolution schemes... 345 17. Zaleski S., 1996-2002,Simulations of highReynoldsnumberbreakupof liquid- gas interfaces, von Karman Institute for Fluid Dynamics, Lecture Series Porównanie własciwości schematów o wysokiej rozdzielczości CICSAM i HRIC Streszczenie Tematem niniejeszej pracy jest modelowanie przepływuw układzie płynów o du- żym stosunku gęstości i lepkości z dobrze określoną powierzchnią rozdziału: wody i powietrza. Symulacje komputeroweprzeprowadzono za pomocąwłasnego kodu z za- implementowanąmetodą Volume of Fluid (VOF), w której pozycja powierzchni roz- działu zostaje określona poprzez śledzenie rozkładu objętościowego współczynnika wypełnienia w czasie. Dwa schematy o wysokiej rozdzielczości CICSAM i HRIC słu- żące do dyskretyzacji członu konwekcyjnego w równaniu transportu objętościowego współczynnikawypełnienia zostałyporównanewdwóchprzypadkach testowych: rota- cji ciała ze szczelinąoraz ruchuwodywzbiornikuwywołanegodziałaniemzewnętrznej siły.Głównymwnioskiem zwykonanych symulacji oraz porównaniawynikówobliczeń z danymi doświadczalnymi jest fakt zależności otrzymanychwyników od właściwości numerycznychanalizowanychschematów, aw szczególności sposobu rekonstrukcji po- wierzchni rozdziału.Wprzypadku niejednorodnejmieszaniny płynów obawspomnia- ne schematywprowadzajądodatkoweprzesunięciew czasie do historii zmian ciśnienia oraz tłumienie ekstremalnych wartości amplitudy ciśnienia wskazujące na nieprawi- dłową rekonstrukcję dynamikimieszaniny. SchematCICSAM rekonstruje zachowanie płynóww sposób bardziej zbliżony do danych doświadczalnychw prównaniu ze sche- matem HRIC. Dowodem na to jest mniejsze przesunięcie w czasie względem danych doświadczlnych orazmniej wyraźne tłumienie ekstremalnych wartości ciśnienia. Manuscript received January 8, 2007; accepted for print February 7, 2008