Numerical approach for generic three‐phase flow based on cut‐cell and ghost fluid methods

In this paper, we introduce numerical methods that can simulate complex multiphase flows. The finite volume method, applying Cartesian cut‐cell is used in the computational domain, containing fluid and solid, to conserve mass and momentum. With this method, flows in and around any geometry can be simulated without complex and time consuming meshing. For the fluid region, which involves liquid and gas, the ghost fluid method is employed to handle the stiffness of the interface discontinuity problem. The interaction between each phase is treated simply by wall function models or jump conditions of pressure, velocity and shear stress at the interface. The sharp interface method “coupled level set (LS) and volume of fluid (VOF)” is used to represent the interface between the two fluid phases. This approach will combine some advantages of both interface tracking/capturing methods, such as the excellent mass conservation from the VOF method and good accuracy of interface normal computation from the LS function. The first coupled LS and VOF will be generated to reconstruct the interface between solid and the other materials. The second will represent the interface between liquid and gas.


INTRODUCTION
Simulating multiphase and multimaterial flows is some of the most challenging problems in computational fluid dynamics due to the presence of numerous phases or materials and also to the difficulty of interface treatment. In order to model accurately the physical phenomena between phases or materials, it is crucial to well predict flow fields close to the interfacial layer. The first challenge is the interaction between the solid and fluid phases, when solid is considered as undeformable material. In recent years, applications of Cartesian grids with immersed boundaries between solid and liquid have become increasingly popular. Unlike body-fitting methods where solid boundary is treated explicitly, the immersed boundary method handles the boundary implicitly. The boundary condition may be represented by inserting body force terms into the cell containing solid so that the nonslip condition will be satisfied. [1][2][3][4] Other methods like the ghost cell approach 5-7 define a virtual layer of nodes, which are located inside the body and have at least one adjacent cell in the fluid computational domain. The flow field variables at the ghost nodes are computed based on these values at the neighboring cells and the boundary condition applying on the body surface. The ghost cells then are used for spatial discretization

GOVERNING EQUATIONS
The model equations are derived based on formal volume and ensemble averaging. [24][25][26] An important element is that, when based on volume fractions, accurate boundary positions can be located and correct boundary conditions can be applied at internal and external boundaries. The transport equation for the mass is given as follows: where wall area A w , fluid volume V F , fluid surface S F , wall velocity u w , surface normal vector n F , and wall normal vector n F,w are explained in Figure 1. We will use F as fluid fraction of the control volume. After evaluating the volume and surface integrals, Equation 1 can be rewritten as follows: where ΔV is size of the control volume V and A F is the fluid surface area, which can be determined based on the LS function. The term (u − u w ) · n F,w A w is the change of mass in the cell due to the mass flow across the wall face. For the impermeable inert and steady wall, this term will be negligible. Similarly, the momentum equation reads The volume integrals are first evaluated, d dt ∫ V F udV = ΔV t ( F u) and ∫ V gdV = ΔV F g, where g is the gravity acceleration vector. Next, we do the surface integrals where w is wall shear stress. In this study, a simple formulation is used to compute this term, which is where Δh is the distance from velocity nodes to the wall and is the value of the LS function in the corresponding cells. u (u − u w ) · n F,w A w is the momentum exchange term due to mass transfer. In the case of an inert wall surface, moving through space, we will have zero contribution from this term. It can be applied for a typical fluid-structure interaction problem.

NUMERICAL METHODS
In Figure 2, we depict a typical staggered grid layout in 2D. The LS function and volume fraction are stored at the pressure node. The boundary and solid cells are highlighted. While the location of pressure in unchanged for both standard and

FIGURE 2
The staggered grid layout in 2D. The pressure node is marked by a circle(⊕). Horizontal and vertical arrows represent the location of u and v nodes, respectively boundary cells, the velocity is located at the face center of the open pressure cell surface. In this paper, we assume that the fluids are incompressible and there is no mass transfer across the interface. Therefore, Equation 2 can be discretized in time as follows: where n and n + 1 denote the current and next time steps; Δt is the time step size. Note that the n+1 F is not computed from this equation. For a stationary wall, n+1 F = n F . For a moving wall, we assign n S as the solid volume fraction in the control volume. It can be calculated by where S S is the solid surface. In addition, we also have As a result, The predicted velocity u * can be evaluated from the Navier-Stokes equations as follows: In Equation 10, the convective term is treated explicitly by the second-order Adams-Bashforth scheme. The semi-implicit Crank-Nicolson scheme is employed for discretizing the diffusive terms. For the wall shear stress, a fully implicit scheme is applied to eliminate the stability restriction from this term. The pressure contribution is absent in this formulation due to the stiffness of the moving boundary problem. The velocity at the next time step u n+1 is calculated as follows: By inserting Equation 11 to Equation 6, a Poisson-like pressure equation reads After solving Equation 12 to obtain p n+ 1 2 , the new velocity u n+1 can be computed by applying Equation 6.

The convective flux
In this section, we present a method to estimate the convective flux for the u-momentum equation. This method can be used for the v-momentum equation as well. Based on the cut-cell method introduced by Kirkpatrick et al, 11 the convective flux in the cut-cell is calculated for the u-momentum equation as follows: In the x direction, for the standard cell, a typical central interpolation is used to compute the velocity at the center of the cell face . For a boundary cell, the interpolated velocity is slight off the center of the cell face, as shown in Figures 3 and 4. Therefore, a modification is needed to correct the velocity at this position . In y direction, v n at the north face is calculated as where e = Δx w Δx we . Velocity at the north face is given as The correct velocity at the cell face center reads

The diffusive flux
The diffusive flux for the u-momentum equation is given as follows: for direction.
As seen in Figure 5, the vector connecting new velocity locations E and P may not be perpendicular to the cell face. Therefore, a modification from a convectional central difference is needed to compute the derivative u x and u at the cell face. Taking the derivative of u along the vector S gives where s x and s y are component of S. Using the central difference to approximate u s yields Therefore, where n y is a component of normal vector N at the surface, which passes through e. The velocity u e is evaluated by Equation 14. As a result, Equation 22 becomes

Small cell problem
The presence of an interface creates several velocity cells, which connect to only one pressure cell. Those cells are marked as small cells (slave cells) and are linked to master cells, as shown in Figure 6. The details of this method was presented in the work of Kirkpatrick et al. 11

FIGURE 6
Linking between a slave cell and a master cell

Extension of the cut-cell method for two-fluid flow
As reported in chapter 4 in the work of Tryggvason et al, 27 the cut-cell method can be applied in fluid-fluid interaction directly. In line with this suggestion, we apply the present cut-cell approach for two-fluid flow. We adopt the embedded boundary method 12 to solve the Poisson-like pressure equation. We assign separate flow variables for each fluid, as shown in Figure 7. The predicted velocities for each phase are obtained from the momentum equations. The jump conditions across the interface are used to find a relation between fluid pressures. The mass equations are solved for each fluid to compute the corrected value of pressures and velocities. However, the resulted matrix from discretizing the mass equation is nonsymmetric and the huge variation in fluid densities causes the matrix coefficients to vary significantly. This means that the linear system is very stiff and difficult to solve by well-known methods such as BiCGSTAB or GMRES. Furthermore, it is nontrivial to construct a polynomial for evaluating the pressure gradient at the free surface. We could not find a solution for these challenges, so we choose a ghost fluid method as an alternative approach.

The ghost fluid method
The presence of gas phase in V F can create some problems, because each pressure cell carries only one pressure value. When we try to solve the Poisson-like pressure equation for a cell containing gas and liquid, the same pressure gradient will be applied for both phases. However, the huge difference between liquid and gas densities will cause large variation in the corrected velocity. To tackle this problem, a ghost fluid method is adopted. The idea behind this approach is that we give the distinct pressure gradient for each phase. This gradient is computed basing on a ghost pressure value at the other phase pressure cell. By using the interface jump conditions, the ghost value can be derived from the actual pressure locating at that cell. For simplicity, we present the method as 1D, while an extension for 2D and 3D can be done by a similar technique. We define the fluid LS function F , which denotes the free surface as F = 0. We assume that the centroid of cell i is inside the gas field ( F > 0) and the centroid of cell i + 1 is inside the liquid field ( F < 0) p ghost l,i and p ghost g,i+1 are ghost values of liquid and gas pressure at adjacent cells. Without the present of mass transfer between two phases, the velocity is continuous across free surface. Therefore, the velocity jump condition at the interface yields In the homogeneous model, the predicted gas and liquid velocities are identical (u * g = u * l ). As a result, inserting Equation 26 to Equation 25 gives The liquid and gas pressure at interface are calculated as follows: In our study, we neglect the surface tension and the possible discontinuity of shear stress. Therefore, the pressure is continuous across the interface, which gives the pressure jump condition From Equations 28 and 29, a relation between two ghost values is expressed as follows: From Equations 27 and 30, the gas ghost value can be computed as Inserting the gas ghost value from Equation 31 to Equation 25, the updated velocity is calculated as In this study, we assume that liquid and gas share the same velocity (u i,g = u i,l = u i ). Therefore, we can use Equation 32 for solving the pressure equation in Equation 12.

The density-based computation of the convective flux
There are several proposed methods, which solve the problem of the high density ratio in the momentum equations such as the ones reported in the works of Raessi, 28 Bussmann et al, 29 and Rudman. 30 In order to fit the Adams-Bashforth scheme for the convective term, the flux correction method developed by Desjardins and Moureau 31 is adopted. We denote C n as the convective term at time step n. From Equation 10, C n would be expressed as The new convective flux is given as follows:

FIGURE 8
The volume of fluid in a u cell wherên and̂n +1 are multiphase densities in the velocity cell at time step n and n + 1, respectively. The previous density is calculated aŝn = n l l + n g g with n g and n l as volume fractions of gas and liquid inside the velocity cell. The new densitŷn +1 is updated based on the continuity equation Note that the primary volume fraction is stored at the pressure cell. Therefore, an interpolation technique is needed to compute this parameter in the velocity cell. An example for how to compute the volume fraction in the u cell is given in Figure 8. It is computed by combining the amount of fluid in the left haft of the pressure cell i and right half of the pressure cell i + 1. At the cell face, the density is interpolated based on the first-order upwind scheme. More complex scheme can be found in the work of Rudman. 30

The coupled LS and VOF method
The transport equation for liquid volume fraction ( ) and LS function ( ) is given as Moreover, it is equivalent to The second-order operator split in the works of Puckett et al 32 and Strang 33 for the LS function and the VOF is employed in this study. First, we solve for the x direction and then the y direction. In the next iteration, we switch calculation order. Time and space discretization in the x direction * and in the y direction, The LS function at a cell face is calculated by taking an average of this parameter over the adjacent cells which is proposed by Sussman. 34 It is simpler and less time-consuming than the high-order Essentially Non-Oscillatory (ENO) or Weighted Essentially Non-Oscillatory (WENO) scheme. In order to compute the VOF that passes a cell face at a specific time, an where n(n x , n y ) is interface normal vector. d is the distance from free surface to the cell center. The interface orientation could be calculated by discretization of n = ∇ |∇ | or n = ∇ |∇ | . However, due to the discontinuous nature of the VOF, the -based computation may give an incorrect normal approximation. Therefore, the LS function, which is continuous across the interface, is chosen to estimate the interface normal vector. Usually, a central second-order difference discretization ∇ is used for normal calculation. However, irregular shapes of the free surface can sometimes degenerate the solution. Figure 9 gives an example of when the central difference of n = Δ |Δ | fails to compute the normal vector accurately at (x i , y j ). As we can see from the figure, the LS function at (x i+1 , y j ) gives a distance with respect to interface II, while the one at (x i , y j ) gives a distance with respect to interface I. Therefore, a central difference in the x direction will give a wrong estimation of x . We set "extrema band" as name of the points. In this study, direction difference 35 is employed to minimize the error in normal computing. In the work of Macklin and Lowengrub, 35 "normal quality function" is defined as Q(v) = |1 − |v||, where v is the central second-order difference of ∇ . The Q i, j denotes the evaluation of ∇ at (x i , y j ).
is the threshold value, which is used to determine whether the point (x i , y j ) belongs to extrema band (Q i, j ≥ ). In the present work, = 0.1 and = 0.075 are used for 2D and 3D applications. The direction vector D = (D x , D y ) is defined based on Q(v) as follows: D y is determined similarly. The discretization of x is calculated by For the case when D x is undetermined, the central difference is used instead of the method proposed by Macklin and Lowengrub. 35 For the y direction, y is computed accordingly. After obtaining the normal vector, the distance d is determined based on volume fraction . A detail about this technique was given by Griebel and Klitz. 23 After solving Equations 38 and 39, the LS function will be reinitialized as shortest distance by using the interface configuration in Equation 40. For details of the algorithm, the reader is referred to the work of Wang et al 22 and Griebel and Klitz. 23

Taylor-Couette flow
This test is performed to check the order of accuracy of the cut-cell method. A diagram of Taylor-Couette flow is represented in Figure 10. It consists of two cylinders with a gap between them, which allows the fluid to move freely. The inner cylinder rotates with angular velocity , while the outer cylinder is stationary. The Taylor number Ta, which presents characterization of the Taylor-Couette flow is defined by where R 1 and R 2 are the inner and outer cylindrical radius, respectively. In this test, we set R 1 = 1 and R 2 = 4. As reported by Dou et al, 36 the flow fields are stable with Ta smaller than 1708. According to Cheny and Botella, 37 the exact flow fields in steady state are given as follows: where The infinity norm (L ∞ ) and 2-norm (L 2 ) of the primitive variables' error are computed by where PV n and PV e are numerical and exact solution of flow field variables. In this paper, the Ta is equal to 1000 and the center of the cylinder(x c , y c ) is at (0.023, 0.013). The computational domain is from −5 to 5 in each direction. The grid spacing h is approximated by 10∕N, which N is grid size. The Figure 11 shows the convergence rate of the numerical

FIGURE 10
The schematic of Taylor-Couette flow scheme for L 2 and L ∞ of the velocity field error. While the current method shows second order of accuracy for the L 2 of the u and v velocity errors, the L ∞ is slightly off from the second-order slope. The quadratic convergence is also achieved for L 2 of the pressure error, as indicated in Figure 12A. However, its L ∞ only shows a superlinear convergence.

Flow past a circular cylinder
The second test is aimed at examining the capability of the numerical method in terms of predicting physical quantities. The circular cylinder is chosen in this case due to numerous numerical and experimental research studies reported in the literature. Furthermore, its curved surface creates a wide range of cell shapes, which is notable for investigating the robustness of the cut-cell method. For cylinder flow, its solution strongly depends on Reynolds number (Re). The flow will transform from steady to unsteady state when the Reynolds number is greater than 40 − 50. 11 A diagram of the computational domain is given in Figure 13. A cylinder of diameter D is placed at a distances 20D and 40D from inlet and outlet to minimize the effect of the boundary conditions on the flow field around the cylinder. Slip walls are located 10D from the body center, at the sides of computational domain. For all our computations, the inlet velocity U inlet is calculated based on Reynolds number, object diameter, and fluid kinetic viscosity as U inlet = Re D . Due to the big ratio of computational domain to cylinder, the mesh is generated nonuniformly except for region surrounding the object.   The details of grid sizes using for our simulation is given in Table 1. The pressure, drag, and maximum lift coefficients are calculated by this following formulation: Here, p 0 and U ∞ are inlet pressure and reference velocity, and F D and F L are drag and maximum lift force acting on the cylinder. The present results of the pressure coefficient for Re = 40 achieved on the group of meshes are compared with experiment data 38 in Figure 14. As seen from the figure, good agreement with these reference results is observed. In addition, the current method can predict the pressure distribution quite accurately, even for the coarse grid A. At Re = 40, the flow is steady and creates a vortex ring at the wake. Several characteristics of the wake such as separation angle and length are shown in Table 2

Transversely oscillating cylinder in a free stream
To validate our numerical method in case of a moving boundary problem, we performed a computational simulation for flow over the cylinder with prescribed cylinder oscillation in the vertical direction. The computational configurations are chosen similarly to the previous test in Section 4.2. In this section, all our computations are implemented on the grid C. The translational motion of the cylinder center (x c , y c ) is given by where A = 0.2D is an oscillation amplitude, and f e is an excitation frequency. The free-stream U ∞ is estimated so that the Reynolds number is equal to 185. First, we run a simulation for the stationary cylinder to obtain the natural vortex shedding frequency, f o = 0.195. Afterwards, the excitation frequency can be computed by f e = k f o , which k ∈ [0.8, 1.2]. Figure 15 shows the vortex shedding formed behind the cylinder when the cylinder is at its maximum vertical position. The shed vortices behind the cylinder become numerically dispersed due to the coarser grid downstream. The time evolution of computed drag and lift coefficients for a variety of forcing frequencies is presented in Figure 16. For f e smaller than f o , both C D and C L show a uniform oscillation after a certain time. When f e is greater than f o , the force coefficients shows a nonuniform fluctuating pattern due to the impact of the higher excitation frequency. The computed mean value of drag and lift coefficients are in Figure 17 compared with other numerical results from both a cut-cell method 37 and body-fitted grid. 43 While the root mean square of the drag coefficient shows a good agreement with the reference data, the other quantities show a slightly overprediction. Figure 18 shows a comparison of the pressure coefficient over the cylinder surface. Clearly, our numerical result marks well with the body-fitted grid. 43

FIGURE 15
The vorticity contour [Colour figure can be viewed at wileyonlinelibrary.com]

Flow past a sphere
This test is used to validate our cut-cell method for a 3D application. The computational domain and numerical mesh are shown in Figure 19. The sphere is placed 20D from the inlet and 40D from the outlet, where D is sphere diameter. Four slip walls are located symmetrically at the sides of the domain. A nonuniform mesh of 121 × 81 × 81 grid points is employed in this simulation. We cluster a uniform grid in the vicinity of the sphere with 20Cells∕D, so that the physical phenomena around the body can be captured accurately. As reported in the work of Johnson and Patel, 44 for Reynolds number smaller than 210, the flow is steady and symmetric, as can be observed from Figure 20A. A separation length L∕D and vortex position (x c , y c ) for two specific Reynolds number obtained in various studies is given in Table 5. Clearly, our numerical method can predict quite accurately physical quantities for steady flow over a sphere. When the Reynolds number increases, the flow becomes unsteady and nonsymmetric as indicated in Figure 20B. Furthermore, the time evolution of the drag and lift force coefficients, for Re = 300, are given in Figure 21A. After t * = 100, these quantities show a stationary oscillation. As a result, the lift coefficient is chosen to calculate the Strouhal number. The computed result and comparison with well reported studies in the literature is presented in Table 5. It can be seen that our result is   Figure 21 compares the computed drag coefficient with other numerical results 5,44 and theoretical data. 46 In general, our method can predict the drag coefficient well and in agreement with previous studies.

Advection tests
The tests are designed to validate the ability of CLSVOF scheme to deal with complex topological transformations. In these tests, the transport equations for the VOF and the LS fields will be solved with a specified velocity field. Therefore, the resulting accuracy will only depend on the interface tracking method itself rather than the numerical methods applied for the governing equations. A fluid circular disk of radius 0.15 is located with its center at (0.5, 0.75), inside a unit computational domain. A deforming velocity field is prescribed as follows: where T = 8s is the time period that the fluid disk needs to come back to its original location. Figure 22 shows the development of the interface at t = 0s, 4s and 8s with the grid size of 128 × 128. The fluid circle starts moving in the clockwise direction from its initial position at t = 0s and achieves maximum stretching at half of the time period, t = 4s, before returning to its initial position at t = 8s. As shown in Figure 22B, the interface is not continuous and is broken into several small pieces at the tail. This phenomenon is due to the limitation of the VOF field, which only accepts a single interface inside a cell. Therefore, when filaments become thinner their cross section will be smaller than the grid resolution. The result is interface cracking. This problem could be reduced by refining the mesh. Deformations of the circular disc, calculated by two different schemes for the normal vector approximation, at t = 4s and t = 8s, are shown in Figure 23. Clearly, the new proposed scheme can preserve the interface shape better at both time marks. A grid convergence study for the deforming flow is conducted by performing numerical computation for three different mesh sizes, which are 64 × 64, 128 × 128, and 256 × 256, respectively. The LS function and area losses are examined in  function reduces from second order to first order when we refine grid, while it remains second order for area loss. This phenomenon is due to the mathematical technique computing the LS function from the VOF, when the grid point at the farthest distance from the interface would be the less accurate.

Dam-break flow
Dam-break flow describes the instantaneous movement of a water column in a tank caused by the sudden removal of a vertical obstacle. The wave tip is driven by gravity and travels freely downstream to a steep wall at the end of reservoir. After colliding with the wall, the wave rolls back and impacts with underlying water. The complex physical behavior and severe interface deformation of the flow is a good test case for validating our numerical method. An overview over the simulated dam-break flow configuration is given in Figure 24. The tank dimensions and the water level measuring positions used in the experiment are given in Figure 24A. The computational domain is illustrated in Figure 24B; the liquid column is initialized with the height of 0.3m and width of 0.6m at the left side of the tank. The remaining reservoir is occupied by air. The liquid and air densities are 998.2 kg∕m 3 and 1.205 kg∕m 3 , respectively. The grid size of 161 × 60 is employed in this computation. Figure 25 compares the free surface profile obtained by two different numerical methods calculating the convective flux. In general, a typical scheme can produce spurious interface shapes at particular moments of the simulation, as represented in Figure 25A and 25B. This phenomenon could be eliminated by applying a new method in Section 3.3 to compute the convective flux, as depicted in Figure 25C and 25D. Figure 26 shows the shape of interface captured by experimental work on other numerical simulation and the present study at different times. It can be seen from    result agrees well with reference data. Figure 28 shows temporal variation of water level at specific locations H1, H2, H3, and H4 (given at Figure 24A) from current results and measured data 47,51,52 for the case of water height H = 0.3m. It can be seen that our present method can predict the arrival time of the first wave across the water level measurements H2, H3, and H4 quite accurately. Our model also predicts the development of the liquid height, after arrival of the primary wave, relatively precise. When the second wave arrives and passes through the given locations H1, H2, H3, and H4, there is difference between our computations and the results previously published in the literature studies in both its arrival time and its elevation. We may note that the dam-break flow involves several complex phenomena such as turbulence or fluid-structure interaction. These phenomena differ from case to case and it affects considerably the flow field properties. As a result, it creates random variations of the elevation in the experiments and there is a need to perform ensemble averaging of experiments in our studies, which did not investigate the effects of turbulence and fluid-solid interaction.

Water entry flow
The next test is a water entry simulation, which is used to demonstrate the capability of our numerical method to handle three-phase flows. The shape of the object is given in Figure 29A. The characteristics of water entry flow are illustrated in Figure 29B. is the wedge angle, being 60 o , 90 o and 120 o in our computations. V is the entry velocity. For simplicity, the item is kept standstill and the free surface is let moving with a speed being equivalent to a water entry velocity. V so is water jet speed; o is the angle between splash and a horizontal interface. The computational domain is illustrated in Figure 30A. The size of the rectangular domain is 0.51m of width and 0.32m of height, as reported in the work of Vincent et al. 53 The wedge width d = 0.036m is held constant for the whole simulation. The nonuniform grid of 260 × 180 is shown in Figure 30B. The mesh is clustered in the wedge vicinity with 100cells∕d so that it can resolve the thin geometry of the splash. The jet velocity versus the entry velocity for the case of = 90 o and = 120 o is plotted in Figure 31. Generally, our method overpredicts the jet velocity and underpredict the jet angle in comparison with reference data from the work of Vincent et al 53 for = 90. When = 120, our computed jet quantities fit quite well with the experimental and theoretical results. The time evolution of splash captured by experiment and numerical study is presented in Figure 32. The present method predicts the interface moving faster than experiment at t*=0. 28 and t*=0.56, but good agreement between numerical result and reference data at t*=0.84.

Oil boom case test
The next test is an oil boom case test. Oil booms are important devices for collecting oils spills from the ocean surface. In this case, we present no validation with experimental data, but rather try to test the capability of our numerical method to predict physical phenomena caused by a solid object moving freely on a fixed grid. Here, the Cartesian cut-cell method enables an efficient fluid structure interaction (immersed boundary method) simulation, under interaction with waves. The schematic of the boom case is given in Figure 33. The object can move freely in the vertical direction. The horizontal position is fixed to approximately mimic the towing operation of the boom. The towing speed is approximately constant. The simulated channel (flow domain) length L Ch is 25 m, while its height H Ch is 8 m. The water depth H w is around 6 m at the ocean surface a wave is generated with the amplitude 0.05 m and time period 1.13 seconds. The boom shape is given in Figure 33, which contains two parts. The cylinder with diameter D B of 0.8 m is located at top, while the boom bar is 1.5 m long and 0.16 m wide. The cylinder center locates at 6 m from the bottom boundary. The boom density is assumed to be 355 kg∕m 3 . The water velocity U w at inlet is 1 m/s and gas velocity U g is 3 m/s. The wave shape and boom position over time is given in Figure 34. Because the equilibrium of the object is below the cylinder center, the boom will travel upward  from t = 0 second to t = 0.123 second. Then it moves downward and start oscillating around the equilibrium position as in Figure 35. The oscillation amplitude decreases following time until t = 10 due to dissipation, mainly caused by liquid shear stresses. After that, the movement of the boom become irregular because of interference with surface wave.

CONCLUSIONS
In this paper, we have introduced a numerical approach for generic three phase flow, including two fluids and one solid. The cut-cell method showed good mass conservation when it was applied to regions containing fluid and solid. In the fluid field, the ghost fluid method is adopted to resolve the discontinuity across the gas-liquid interface. In addition, the density-based convection flux is employed to minimize the effect of the high density ratio(large effective density gradients) in the momentum equations. The numerical results show that the proposed scheme can reduce the unphysical deformation of a free surface. For tracking the interface, we applied the CLSVOF method, which yields an excellent mass conservation. A new proposed technique for computing the normal vector indicates an improvement in preserving the interface shape.