Two-phase free-surface flow interaction with moving bodies using a consistent, momentum preserving method

The numerical prediction of two-phase ﬂows with an interface is challenging, to a considerable extent because of the high density ratio at the interface. Numerical results become affected by momentum losses, diverging spurious interface velocities, free surface distortion, and even numerical instability. To prevent issues like these, consistent momentum and mass transport with an additional continuity equation were introduced. In this article, we describe how a consistent discretization was incorporated into our own method and extended for ﬂuid-structure interaction (FSI) with moving rigid bodies. The new method was tested against benchmark simulations from literature conﬁrming that consistent transport modeling gives a signiﬁcant improvement compared to non-consistent modeling for the dynamics of two-phase ﬂows. Newly devised proof of principle FSI simulations with momentum transfer from ﬂuid to body in the presence of a high-density ratio between ﬂuids are introduced that could serve as a benchmark for future studies. The simulations demonstrate that consistent modeling gives an order of magnitude improvement in terms of momentum conservation compared to non-consistent modeling. Simulations with the new method are also compared to FSI experiments from literature. Results obtained with the consistent method are closer to the measurements than results of the non-consistent method. The merit of consistent modeling with and without FSI becomes especially apparent for two-phase ﬂows with a high-density ratio between ﬂuids.


Introduction
This article is about a consistent method for two-way coupled fluid-structure interaction (FSI) of two-phase flows with a high-density ratio such as that between water and air, impacting with moving rigid bodies.Two-phase flows, like those surrounding floating marine structures, are difficult to predict.The structure needs to be designed for violent weather conditions to ensure crew safety [14].Linear approaches fall short of representing highly non-linear events to an acceptable level during these weather conditions.Traditionally, experiments at scale are employed to model the events, but in specific circumstances numerical simulations are also possible.Sophisticated two-phase flow models are needed to evaluate the highly non-linear interaction between flow and structure accurately.
Our approach is based on a one-fluid formulation for incompressible flows.A unique, continuous velocity field is solved and the pressure field is relaxed to represent both fluids [24].A finite-volume discretization of the governing equations is adopted with a staggered arrangement of variables.This type of discretization can suffer from numerical problems caused by large density ratios of the fluids near the interface.
The density is discontinuous over the interface.The interfacial discontinuities of fluid properties are often avoided or not well implemented in combination with the continuous velocity field.A lack of attention to the numerical implementation of the discontinuity in density can lead to numerical instabilities and non-physical flow features.When the density ratio between fluids increases, the flow becomes more interfacially driven, adding to the numerical error [33].The error does not become apparent in the heavier fluid but rather in the lighter [5].
Errors can be introduced by the staggered arrangement of variables.The arrangement features different positions of the control volumes for mass conservation and momentum conservation.Errors in momentum will appear for non-matching momentum and mass fluxes.Following Rudman [34], such a method is called "non-consistent" in the remainder of the article.The momentum flux becomes dominant over the mass flux in case the density ratio increases, increasing the magnitude of the errors [33].Numerical error accumulation near highly deformed interfaces can lead to failure of the method or can result in non-physical features [25].
Over the years, many have proposed strategies to mitigate these issues within a wide range of numerical contexts [25]; for example, Li et al. [22] for a moment-of-fluid method, Vaudor et al. [39] for a staggered arrangement with a coupled level-set Volume-of-Fluid (CLSVOF) method, Le Chenadec and Pitsch [20] and Owkes and Desjardins [27] for sharp interface Volume-of-Fluid (VOF) methods, and Jemison et al. [17] and Duret et al. [10] for compressible flow solvers.The methods referred to are similar in the way of discretization but different in the reconstruction of the interface and transport.
Earlier numerical schemes based on VOF [34,5] were not sensitive to the numerical instabilities induced by a large density ratio.The method of Rudman [34] uses, together with a staggered arrangement of the variables, grid refinement for the transport and reconstruction of the interface.On the other hand, the method of Bussmann et al. [5] does not need Rudman's grid refinement because of a collocated arrangement.From the authors' point of view, the last two methods mentioned use a "consistent" coupling of the momentum and mass flux together with an unsplit transport scheme.The key point made by Rudman [34] and Bussmann et al. [5] is that in order to obtain a conservative coupling between mass and momentum fluxes, the density of the momentum flux in the convective term needs to be corrected with the mass flux obtained from the VOF transport equation.Bussmann et al. [5] showed that the collocated arrangement results in easier enforcement of consistency between the momentum and mass transport but introduces difficulties with calculating the fluxes.
In the context of level-set methods, Raessi and Pitsch [33] introduced a consistent, stable method.The method makes use of two geometric reconstruction sweeps of the free surface at two different time levels.The disadvantage of the method is that it is limited to one-and two-dimensional problems.The methods of Ghods and Herrmann [15] for collocated unstructured grids and of Desjardins and Moureau [9] and Nangia et al. [25], adapting the method of Rudman [34] to level set methods, use an extra continuity equation.The extra continuity equation results in a auxiliary density field that is coupled with the momentum flux in the convective term in the momentum equation.Both methods [15,9] showed that non-geometric construction of the interface from a level-set field and grid refinement can also result in consistency between momentum and mass transport.However, being based on level-set, the methods were neither mass nor momentum conserving and used first-order upwind for density and velocity transport, producing diffusive flow features.
Patel and Natarajan [29] achieved higher-order accuracy with a consistent scheme.Patel and Natarajan [30] came up with a consistent convective scheme for momentum and algebraic (non-geometric reconstruction) VOF transport.They showed that any non-consistency in the scheme leads to poor accuracy.
Recently, Zuzio et al. [46], using a staggered grid with CLSVOF for the interface between fluids, describe a consistent momentum and mass method with a temporary continuity equation.Compared to Rudman [34], they do not apply grid refinement for transporting the free surface, thereby reducing the computational effort.We will adopt the strategy of Zuzio et al. [46] into our own method [12].
Arbitrarily shaped bodies can be represented on fixed Cartesian grids by means of immersed boundary methods (IBM).An early such method is that of Peskin [32].A distinction can be made between body forcing methods, using source terms in the governing equations to impose the boundary conditions of the body, and cut-cell methods, in which the position of the boundary is reconstructed and part of the discretization.Whereas a body forcing IBM method is straightforward to implement [41,32], the interface between body and fluids can become diffusive.Lagrange et al. [19] combined a body forcing IBM with the consistent momentum and mass method (CMOM) of Zuzio et al. [46] for non-moving bodies.Because the interface between fluids in our method [12] is also reconstructed, we choose to incorporate a cut-cell method for representing the interface between body and fluids.Drawbacks of cut-cell methods are that work is involved in reconstructing the boundary [26,37,38,6,40] and that special treatments of small cut cells are needed [40].On the other hand, the reconstruction adds to the accuracy of where boundary conditions are imposed and the governing equations do not change compared to uncut cells.Cheny and Botella [6] proposed a cut-cell method where small cut cells do not need any special treatment, followed by Xie et al. [45] who have recently described a cut cell method for moving bodies that is also consistent, but which does not include two-way coupled fluid-structure interaction.
Fluid-structure interaction with two-way coupled body motion appears not to have been investigated in the context of consistent momentum and mass methods, especially in terms of the momentum transfer from fluids to body when impacts take place.In this article, our former method in van der Eijk and Wellens [12] is extended with a cut-cell method similar to Xie et al. [45] and made consistent by incorporating elements from Zuzio et al. [46] to study two-way coupled FSI.
The article starts with the governing equations and the discretization.The implementation of the consistent method and the coupling with the structure are discussed, highlighting the addition of the cut-cell method for moving bodies.In Sec. 5 our method is tested by means of fundamental verification cases with a large density ratio at the interface.The verification includes newly formulated simulations for momentum conservation with FSI, that could serve as a benchmark for future studies.In Sec. 6 simulation results are compared with 2D experimental results with FSI. Conclusions are formulated in the final section.

Governing equations
The governing equations for two-phase flows of immiscible Newtonian fluids are given in Eq. ( 1) and Eq. ( 2).These equations are formulated in a conservative vector form.The one-fluid formulation is used [24] with a single velocity field and a single pressure field.The continuity equation reads in which u denotes the fluid velocity vector, ρ is the mixture density.The momentum equation is given by where p is the relaxed pressure, μ is the dynamic viscosity for a mixture and F f represents the body forces.The body force term contains gravity and capillary stresses, so that

κnδ ).
A transport equation is solved for displacing the interface between fluids, using velocity field u under the assumption of a no-slip condition between the two fluids where f(x, t) = 0 gives the position of the interface.When using the Volume-of-Fluid (VOF) method, the distance function f is replaced by volume fraction C f .The volume fraction is a measure of the ratio of fluid volumes in a cell.The body will be assumed rigid.It is displaced with a state-space representation of Newton's second law, in which the position of the body x b is found from and the body's acceleration is found from with m b the body's mass and F b the force of gravity together with the force of the fluids on the body.The fluid force on the body is equal to the integrated normal pressure along the body boundary.Viscous stresses are neglected in the determination of the force on the body, because our main interest is in representing the two-way coupled fluid-structure interaction of impacts -of the body with the fluid or vice versa.Impacts take place over time spans that are too short for viscous effects such as boundary layers to develop.

Grid and solution variables
Before presenting the numerical discretization of the governing equations and the details of including two-way coupled fluid-structure interaction in a consistent momentum and mass method, the grid structure is introduced with the definitions and notations needed for discussing the solution method.
A fixed Cartesian grid is employed.Labeling of grid cells is used to account for the position of the interface between fluids and the interface between body and fluids, so that it is defined where surface tension needs to be applied and where reconstruction of the interface needs to take place.We have adopted the labeling system of Kleefsman et al. [18], see Fig. 1, that uses the label B(ody) for cells completely filled with body, the label E(mpty) for cells filled with gas only (or the lighter of two fluids), the label S(urface) for cells with some liquid (or the heavier of two fluids) directly adjacent to E-volumes.Remaining cells are labeled F(luid).By definition, F-cells are not allowed to connect with E-cells.Note that F-cells do not necessarily have to be completely filled with liquid; a cell like this is indicated in Fig. 1 with a black marker.A cut cell, containing both fluid and body, is indicated in Fig. 1 with a red marker.The blue marker in Fig. 1 is a "special" cut cell.It contains two interfaces, between liquid and gas and between fluids and body.The reconstruction of the interfaces needs extra attention; this is discussed in Sec.4.1.Kleefsman et al. [18]; labels B, F, S, and E. Body is indicated by .Fluid is indicated by .Cut cells, always labeled with F, are necessary for representing moving bodies in the grid.A F-cell that is not completely filled ( ).A cut cell is illustrated by a marker ( ).A special cut cell is illustrated by another marker ( ) which needs extra attention because it contains two interfaces.(For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)Fig. 2. Difference in labeling between F-cells (left) and C-cells (right).The gray area represents the geometrical reconstruction based on the fill ratio of the cells.Fluid is indicated by .The F-cell is not reconstructed and therefore illustrated with a hatched pattern representing a mixture of liquid and gas.New are the C-cells in which the free surface is reconstructed for a smoother interface between fluids.

Fig. 1. Labeling of cells as in
In the original labeling system [18], reconstruction of the interface between fluids only takes place in S-cells and, hence, not in F-cells.In Fig. 2a, therefore, the interface is shown to be discontinuous and the F-cell is shown with a hatched pattern to indicate that it is not possible to distinguish between liquid and air in this cell.
In this article an additional label is introduced, called a C(orner) cell.When a F-cell diagonally connects with no more than one E-cell, it gets the C-label.The interface in a C-cell is reconstructed, resulting in a more continuous interface, as illustrated in Fig. 2b.Reconstructing the interface in C-cells and including it in evaluating the curvature results in a better representation of capillary effects and prevents artificial air entrainment.A disadvantage is the need for an additional calculation step during labeling.The additional computational cost is negligible compared to solving the Poisson equation.The benefits of C-labeling will become clear in Secs. 5 and 6.
The standard MAC configuration of staggered variables is adopted, where scalar variables (pressure p, density ρ, curvature κ, volume fractions C f and C b ) are stored in the center of a grid cell, and velocities are positioned in the faces of a cell.Uncut cells are represented in Fig. 3 by means of continuous black lines, with open circles in cell centers.Velocities in the cell faces are represented as arrows.Continuity control volumes, used for discretizing the continuity equation, coincide with grid cells.Momentum control volumes are shifted in space with respect to continuity control volumes, horizontally for the equation describing (discrete) conservation of horizontal momentum and vertically for vertical momentum conservation, so that the velocities in the respective directions are in the centers of momentum control volumes.Variables in the centers of continuity control volumes are given indices i (horizontal) and j (vertical) to refer to their position in the grid; variables in the faces of continuity control volumes are given indices i − 1 2 and i + 1 2 to show that they are positioned left or right of a continuity control volume center and, similarly, j − 1 2 and j + 1 2 to indicate their position above or below that center.Averaging is employed when the value of a variable is required at a different location from where it is positioned in the grid.For instance, when a velocity is required at the boundary of a horizontal momentum control volume at position (i, j) it is determined as the average of the velocities on either side.Density averaging, however, or rather the discretization of the density, requires special treatment in order to conserve momentum.It is explained in more detail in Sec.4.2 when CMOM is discussed.
A cut-cell method similar to Fekken [13] and Xie et al. [45] accounts for the presence of the body in the grid.Continuity control volumes that are intersected by the contour of the body are scaled by the part of the volume that is taken up by the body.That scaling factor is called volume aperture C b and defined as the ratio of the volume open to fluid (so one minus the volume of the body) and the size of the control volume itself.It is shown in Fig. 4a.The faces of continuity control volumes are also scaled by face apertures a b .They are defined as the area of the volume's face open to fluid, divided by the area of the face itself, see Fig. 4b.Apertures take on values between 0 and 1.For the visualization, velocities represented by Fig. 3. Standard MAC configuration (staggered); pressure p is defined in the cell center ( ), the horizontal velocity u field is sampled on the vertical faces (→), the vertical velocity v is sampled on the horizontal faces (↑).The subscripts (i, j) for the position in the grid are defined.The overlap of the continuity control volume (-), with a vertical momentum control volume (-), and with a horizontal momentum control volume (-) is shown.Control volumes for horizontal momentum V m,h are shown in Fig. 5 together with the lower half of a control volume for vertical momentum V m,v .Necessarily, the size of the momentum control volume is half of the size of the continuity control volumes on either side.The volume aperture C b of a momentum control volume in uncut cells is also determined as half of the sum of volume apertures of the continuity control volumes, and the face apertures of a momentum control volume V m are determined as half of the sum of face apertures of the continuity control volumes on either side.Applying that rigor to the visualization of volume apertures in cut cells leads to the representation in Fig. 5b where the part of the horizontal momentum control volume V m,h that is open to fluid is precisely half of the open parts of the continuity control volumes.Note that in the discretization, again, the volume sizes and face sizes are not changed, but rather scaled by means of apertures.

Discretization and solution algorithm
Our main interest is in the type of violent fluid-structure interaction encountered in slamming, when a ship meets the free surface again after disconnecting from the main body of water it is sailing in.We have studied slamming with our older method for multiphase flow with gas and liquid, reported in van der Eijk and Wellens [12].That method was nonconservative and non-consistent; what we mean by these terms will become clear from the discussion of the new, consistent method in this section.The first sign of trouble with the older non-consistent method was with the case of a vertical plate being impacted by a wedge-shaped body of water [3], which, before impact and therefore incorrectly, started showing a large and growing local deformation of the interface between gas and liquid.Results of the simulation are discussed in Sec. 5.
Inspiration for a solution to prevent erroneous interface deformations came from the literature regarding consistent discretization of mass and momentum.Consider the time discrete version of Eqs. ( 1), (2), and (3), using Forward Euler in time for brevity of notation, where superscript n indicates the time level.The equations are integrated over control volumes to obtain the weak form that is the basis of the discretization in space.Our method combines elements from Kleefsman et al. [18], van der Eijk and Wellens [12] and the consistent modeling in Zuzio et al. [46] to obtain a new algorithm that can account for two-way coupled fluid-structure interaction with moving bodies.The algorithm comprises the following within a time step.First, the fluids and the body are transported using Eq. ( 3) combined with the assumption of incompressibility of the fluids and with volume fractions C f indicating the liquid fill ratio of volume apertures in which S is the boundary of a control volume.Volume fraction C f is equal to one when a continuity control volume is completely filled with liquid (the heavier of the fluids) and equal to zero when occupied by air (the lighter of the fluids).
Because the fluids are considered to be incompressible, the continuity equation reduces to The time discrete momentum equation reads in which p n+1 = δ p + p n and from which auxiliary vector field ũ is solved from A similar form to Eq. ( 8) is solved by Bussmann et al. [5] and Raessi and Pitsch [33] with ρ n+1 on the left hand side to prevent unphysical velocities.Auxiliary vector field ū in Eq. ( 9) is solved from whereas auxiliary density ρ is solved from Eq. ( 1) used as a temporary continuity equation that is integrated over momentum control volumes Here, ρ * is a density for which a consistent discretization in space is used.A consistent space discretization for momentum and mass is obtained when the fluxes along boundaries that are indicated by the word 'consistent' are treated the same with continuity control volumes and with momentum control volumes, and when the spatial discretization of ρ n and ρ n+1 is the same on continuity control volumes and on momentum control volumes.Different notations from ρ n and ρ n+1 are used for densities ρ * and ρ to make clear that their spatial discretization is different.
The discrete velocity field at the new time level is solved from Eq. ( 8) and substituted into the discrete continuity equation to obtain a Poisson equation for the pressure.A system of equations that combines the Poisson equation for the pressure with the equation of motion for the moving body is solved, after which velocity field u n+1 is reconstructed from the pressure field and a new time step commences.
The flow chart below summarizes these items.The numbers in the flow chart refer to the numbers of the subsections in which the steps are discussed, with emphasis on where the approach differs from earlier work.

Transport fluids and body
Eq. ( 6) is the basis for the discretization of the transport equation.A straightforward discretization of the transport equation reads in which subscript f ace refers to the sides of the continuity control volume and the flux δC f , i.e. the amount of fluid that is transported from one continuity control volume to the next, is of the form Fig. 6 shows the flux at the face of the continuity control volume (cut cell) with index (i + 1 2 , j), with body, liquid and gas present at the same time.
Compared to Zuzio et al. [46], the flux now also depends on the face apertures a b .The stability of the time discretization is not affected by the apertures [18]; the Courant number still needs to be smaller or equal to one.
The discretization of the transport equation in Eq. ( 12) is subject to errors, caused to a considerable extent by the fact that the fluid is transported in the axis directions separately, instead of in the direction of the velocity vector at once.These errors are reduced when using the COSMIC transport algorithm [21] combined with the correction of Weymouth and Yue [44] to improve mass conservation.The COSMIC split scheme we incorporate consists of multiple transport steps and reconstruction steps [21].The COSMIC scheme, in 2D, reads as follows The terms with boxes around them are the steps in which fluxes are calculated.Between these two steps, a reconstruction step takes place.The effect of these two steps on the momentum conservation is discussed in Sec.5.2.Something that we have not seen discussed in literature is that we also use COSMIC for transporting the interface between fluids and body.
After transport, the interface between liquid and air in continuity control volumes is reconstructed with PLIC [28].Reconstruction keeps the interface sharper compared to algorithms without reconstruction.Also the contour of the body is reconstructed by means of PLIC as it moves through the grid.An example of a piecewise-linear representation of the interface between fluids and the contour of the body in a continuity control volume is illustrated in Fig. 6.The reconstruction of the interface between structure and fluid is similar to the interface reconstruction between fluids [31].
Special attention is needed for cut cells labeled as S-or C-cell.These cells contain the interface between fluids as well as the interface between fluid and structure, both requiring reconstruction.An iterative process is required to accomplish this.
We define an initial contact angle between the free surface and the structure.The contact angle θ 1 is illustrated in Fig. 7a.
Geometric reconstruction of the interface between fluids in a S-or C-cell with PLIC requires a stencil of 3×3 cells.The reconstruction in cut cells needs volume fractions of other cut cells or B-cells, where there are no volume fractions.The volume fraction in these cells is determined with a virtual fraction field [7].The virtual fraction field is calculated based on the contact line represented by white dots in Fig. 7a, resulting in the virtual fraction field illustrated in Fig. 7b.The interface between fluids, then, is reconstructed based on the virtual fraction field.
In somewhat more violent two-phase flows, multiple, separate fluid bodies can be found attached to the structure.An example is shown in Fig. 8a.It is necessary to make the assumption that every fluid body has either two or zero cells connecting with the rigid body where θ 1 is defined.The connecting cell is found as follows: the S-labeled cut cell having the most surrounding E-cells compared to the neighboring S-labeled cut cells is defined as a contact cell, in which angle θ 1 is then imposed.The remaining S-and C-cells are assumed to be independent of θ 1 and reconstructed based on the virtual fraction field.Separate bodies of fluid that are nearby can lead to intersecting contact lines and overlapping virtual fraction fields.One such overlap is indicated in Fig. 8b by the circle.To prevent incorrect reconstruction in surrounding S-and C-cells, the overlap needs to be considered, because otherwise fluid bodies that are separate will artificially merge.The reverse is also true.If the overlap is not considered properly, fluid bodies may artificially spawn smaller fluid bodies (droplets or bubbles).An example where fluid is kept nicely together is the buoyant cylinder in Fig. 27a where a thin film of liquid flows down along the boundary of the cylinder as it moves upward.In the remainder of the article, the angle θ 1 is chosen equal to 90 • .

Discretization of the density
A consistent discretization of mass and momentum is obtained through the discretization of the density.The density in different terms of the equations, also indicated by the different notations ρ n , ρ * , ρ and ρ n+1 , is treated differently.In general, near the interface between fluids, the density on any control volume is obtained through averaging by means of the volume fractions of the fluids It is straightforward to use Eq. ( 16) for finding densities ρ i, j in the centers of continuity control volumes.The discussion below is about how to determine the aforementioned densities in momentum control volumes.

Discretization of ρ n and ρ n+1
Consider the horizontal momentum control volume indicated by the blue dashed line in Fig. 9a.The continuity control volumes on either side are indicated by continuous orange lines.The term (ρu) n in the momentum equation is evaluated as (ρ i+ 1 2 , j u i+ 1 2 , j ) n .In our former method [12] density ρ n i+ 1 2 , j is obtained from weighted averaging using the volume apertures in continuity control volumes Averaging the density like this for non-conservative formulations of the equations does not necessarily lead to the problems described in the introduction, so it does not require a consistent discretization.Our problems began when starting to work with the conservative form of the governing equations.A discretization of the density consistent with the size of the momentum control volume is obtained by dividing the momentum control volume into four sub-volumes as in Zuzio et al. [46] but now taking into account the volume apertures in cut cells.By choosing four sub-volumes instead of two, one makes sure that the discretization of the density is also consistent between horizontal and vertical momentum control volumes.For each of the sub-volumes, the volume apertures are determined, after which the density in a sub-volume, such as ρ i+ 1 4 , j+ 1   4   in Fig. 9b, is determined according to Eq. ( 16).
The density for the term (ρu) n in a horizontal momentum control volume is then found with The term ρ n+1 in Eq. ( 9) is discretized in the same way as in Eq. (18).Something that we have not seen discussed in literature before is the discretization of the viscosity μ n near the interface between fluids.As the viscosity in liquids can be Fig. 9. Density calculation at the center of momentum control volume; the continuity control volume (-) and the horizontal momentum control volume (-), the area used by the averaging/discretization method (-), the cell center ( ), the center of a quarter of a cell ( ).Body is indicated by .Fluid is indicated by .quite different from the viscosity in gases, it seemed 'consistent' to also apply the discretization in Eq. ( 18) to the viscosity μ n .The evaluation of this decision has been made explicit for the case of the 2D rising bubble in Sec. 6.

Convective term with ρ *
The convective term S m ρ * u(u • n)dS in the momentum equation requires discretization of ρ * on momentum control volumes.Our older method [12] used weighted averaging based on the volume apertures in continuity control volumes for ρ * , see Eq. ( 17).The new, consistent discretization is similar to [46], but now with apertures.The objective of the discretization of ρ * is to keep the fluxes ρudS in the transport equation, discretized on continuity control volumes, consistent with the momentum fluxes ρ * uudS in the convective term of the momentum equation.For this reason, the term ρ * u in the momentum flux is discretized by means of mass fluxes through the faces of the momentum control volume as The density ρ * in Eq. (19) in the left face of the horizontal momentum control volume shown in Fig. 10 can then be rewritten as δC f ,i, j ρ f + a b,i, j u i, j δtδ y j − δC f ,i, j ρ a a b,i, j u i, j δtδ y j , (20) with the face aperture of the momentum control volume at that location found from

Determine ρ through temporary continuity equation
The left-hand side of Eq. ( 10) features the term ρ ū.In order to prevent instability, the auxiliary density ρ requires a consistent discretization.Such a discretization can be obtained when Eq. ( 11) is used as a temporary continuity equation, see Zuzio et al. [46].Similar discretizations can be found in Bussmann et al. [5], Raessi and Pitsch [33] and Nangia et al. [25].
For a horizontal momentum control volume, discretization of the temporary continuity equation with apertures to account for cut cells yields (δx i + δx i+1 ) , (22) in which V m,h is the size of the horizontal momentum control volume and and Note the use of the horizontal body velocity u b in Eq. ( 23).This is different with respect to Zuzio et al. [46].

Determine temporary velocity ū
Apart from density ρ * , the spatial discretization of the convective term is comparable to Kleefsman et al. [18] and therefore referred to as C(ρ * {u n , u n b })u n .An auxiliary vector field ū is solved from the discrete form of Eq. ( 10) that reads The operators, like the convective operator C and momentum control volume operator m , include the cell face apertures a b and body apertures C b .The operators are explained in more detail by Kleefsman et al. [18].
For the sake of the discussion, Forward Euler time integration of the momentum equation was used in the notation up to here.When using a second-order upwind discretization in the convective term of the momentum equation combined with Forward Euler time integration, it comes with a Courant stability limit that is too low for practical simulation.Using a one-step Adams-Bashforth time integration of the momentum equation instead yields a larger stability region when using a second-order upwind discretization [42], leading to the criterion that the Courant number should remain smaller or equal to 0.25.The convective term The use of the Adams-Bashforth scheme for the convective term implies that this scheme also needs to be used for the advective term of the auxiliary density field.When doing so, spurious velocities are introduced around the interface.Because the Adams-Bashforth uses two different time levels in one step, it is expected that we need the same for the densities ρ * .
However, the densities at the old time level n − 1 are based on an older volume fraction field.This means that they do not match the VOF mass fluxes at time level n and can lead to density values in control volumes without an interface at time level n.Time integration as in Eq. ( 26) is not consistent and can lead to interface distortions.To prevent interface distortions, Eq. ( 26) needs to be solved as follows where (ρ * ) n is based on the volume fluxes to find C n+1 f .The same is true for the operators accounting for the body.An example of the effect of using Eq. ( 26) instead of Eq. ( 27) is shown for a high-density horizontal translating bubble in a vacuum.The results in Fig. 11 are shown after 0.04 [s] where the bubble has an initial speed of 10 [m/s].

Auxiliary vector field ũ
The spatial discretization of gravity in the external force term F f is the same as Kleefsman et al. [18], and therefore not specifically discussed.F f also contains the capillary force at the interface between fluids.
A continuum surface force (CSF) model is used for modeling the capillary force.The CSF-model imposes surface tension as a body force [2,4], for which interface curvature is an input.The curvature κ is calculated in every S-cell and C-cell, using a local height function and is defined for every continuity control volume.The direction in which "height" is defined depends on the largest component of the interface's normal vector found from the PLIC reconstruction.The combination of the CSF-model and the height function ensures that the curvature calculation is not the governing error resulting from the imbalance between pressure and surface tension [1].A height value H consists of the sum of the volume fractions in three cells in height direction as illustrated in Fig. 12.The curvature of the continuity control volume with indices (i, j) in Fig. 12 is found from where .
For the horizontal momentum cell in Fig. 13 the discretization of the capillary force becomes in which (κ i, j + κ i+1, j )/2, if κ i, j and κ i+1, j are defined and both are labeled as S κ i, j , if κ i+1, j is not defined and labeled as E κ i+1, j , if κ i, j is not defined and labeled as E 0, if κ i, j and κ i+1, j are not defined and both are labeled as E .
The subscripts indicate the position relative to continuity control volume with indices (i, j), where i + 1 2 , j indicates the center of the horizontal momentum control volume shown in Fig. 13.The density ρ n+1 i+ 1   2 , j in Eq. ( 29) is determined as in Eq. ( 18).The discussion of the capillary force completes the discretization of the external force term F f .The spatial discretization of the viscous term is comparable to Wemmenhove et al. [43] and referred to as D n u n , keeping in mind that here the kinematic viscosity μ n near the interface is discretized 'consistently' as in Eq. ( 18).Using these, where and the term ū is an auxiliary vector field with the contributions of the convective term, see Eq. ( 25) (or Eq. ( 27) when using Adams-Bashforth).

System of pressure Poisson equation and equation of motion
The continuity equation is solved for the continuity control volumes.When the densities away from the interface are not allowed to vary, the discrete representation of the continuity equation can be written as with u the vector of fluid velocities, M 0 the discrete divergence operator working on the fluid velocities, u b the velocity of the body imposing a boundary condition on the flow and M b the discrete divergence operator working on those boundary velocities.
When u n+1 is solved from Eq. ( 8), it is substituted into Eq.( 32), which results in a Poisson equation for the pressure The operators n+1 m , M n+1 .The operators are explained in more detail by Kleefsman et al. [18] and Fekken [13].
The motion of the body is solved from to determine the area of the body that the pressure acts on.Multiplying the operator with the pressure field results in a force vector.
Using an auxiliary body velocity and the Poisson equation in Eq. ( 33), the coupled system of fluid and body can be written in compact form as where A represents the combination of the operators on the left-hand side of Eq. ( 33) associated with δ p. Solving the system will find the pressure change δ p and the body velocities u b .

New velocity field
The new fluid velocity field is solved from the pressure change δ p as

Proof of principle simulations
Our method is verified by means of several proof of principle simulations that were specifically designed for this article.These simulations show the properties of the method with regard to preserving mass, momentum and the shape of the interface.The implementation of the method that is called "consistent", abbreviated with "C", is compared with the nonconsistent implementation ("NC") of van der Eijk and Wellens [12].The non-consistent method makes use of weighted averaging to determine densities at control volume boundaries.Compared with the consistent implementation described above, the former non-consistent method uses no auxiliary density field, no density discretization based on the mass fluxes through the boundaries of a momentum control volume, and the non-conservative momentum equation is solved.

Moving wedge
The first moment "problems" were encountered because of using a non-consistent method was for the case of a moving wedge [3].The case represents a liquid wave impact against a wall in a LNG container undergoing motion.The wave is represented as a wedge, because the objective in the original article was to study the effect of rise time on the dynamic deformation of the container wall.The fluid configuration is illustrated in Fig. 14a.The domain is 1 [m] by 1 [m] with 80 grid cells in each direction.That means that the height of the wedge is captured by 16 cells.The density of the liquid is 10 3 [kg/m 3 ].The air density is set equal to 10 −9 [kg/m 3 ] to exclude effects that air may have on the wave impact.Typically we find that the bigger the density ratio between liquids, the larger the issues associated with non-consistent modeling.The wedge shape should be preserved before the impact takes place.
The simulation result with NC, however, shows a distorted interface between fluids.It is shown in Fig. 14c for time instance 0.04 [s].Furthermore, it is expected that the momentum in vertical direction remains zero before impact.The total vertical momentum I vert is calculated at each time step as where dV m,v,i, j+ 1 2 is the size of the vertical momentum control volume.The evolution of I vert over time is shown in Fig. 14b.It is clear that the vertical momentum for the non-consistent method shows a spurious increase over time.From Zuzio et al. [46] came the inspiration that the shape deformation could be due to the non-consistent discretization.Where the spurious momentum for a low-density ratio scales with the flow velocity, for a high-density ratio the spurious momentum change due to non-consistent modeling becomes dominant over the changes that scale with the velocity.
With the former method (NC.), momentum transport and mass transport are not consistent; they are decoupled.The non-consistency results in spurious velocities around the interface.Fig. 14d shows the result of applying the new method (C.).No strange interface distortion is visible anymore.

Diagonally translating high-density droplet
After the wedge case with grid-aligned translation, the next case considers diagonal translation and the effect of a split scheme.A diagonally translating droplet with a high density compared to the surrounding gas is modeled.The case was proposed by Bussmann et al. [5] and also used in Zuzio et al. [46].This case is designed to be a good indicator for non-consistency between mass and momentum transport, as the high-density ratio quickly leads to momentum losses and interface deformation.
The density ratio, 10 6 over 1 [kg/m 3 ], is the same as used by others [15,33,46].The large ratio will mainly affect the convective term in Eq. (25).The effect of surface tension and viscosity is neglected.The gravity constant is zero.A square domain with sides of length 1 [m] is used with an initial velocity field of 10 [m/s] in x and y directions throughout both the high-and the low-density fluid.The simulated time is 0.05 [s].The radius of the droplet is 0.15 [m].A Courant restriction of 0.2 is used.The setup is illustrated in Fig. 15.With the large density ratio at the interface, we expect that the interaction between fluids is small and that the shape of the droplet is preserved.The kinetic energy should be conserved during the simulation.It is calculated using the following equation where u and v are the horizontal and vertical velocity, respectively.The subscripts indicate the position of the velocity in the staggered arrangement of variables within a cell.The energy losses over time are depicted in Fig. 16a for the nonconsistent and consistent method.The new consistent method is better at preserving kinetic energy in comparison with the original non-consistent method.Even coarse meshes demonstrate good energy-preserving behavior with the consistent method.Fig. 16b shows an enlargement of the graph of energy over time for the consistent method to investigate grid convergence.The energy loss becomes smaller for higher grid resolutions.The (small) spikes in energy are related to the droplet propagating through cells.
All results in this article were obtained with the COSMIC interface advection method.The densities ρ * and ρ, see Eqs. ( 20) and (22), are based on the mass 'fluxes 1' calculated in the first step of the COSMIC scheme, see the box with number 1 in Eq. ( 14).One of the results, called 'fluxes 2', was obtained with densities ρ * and ρ based on the mass fluxes indicated by the box with number 2 in Eq. ( 15) to evaluate the differences.The simulation with these densities based on the fluxes indicated by 2 had the same setup as described above, and a grid resolution of 120×120.Its result in terms of the kinetic energy is shown in Fig. 16b.Result 'fluxes 2' is accompanied by higher velocities near the interface, with, hence, smaller time steps to satisfy the Courant criterion.The increased computational effort involved in using mass fluxes 2 does not lead to improved kinetic energy preservation.On the contrary, Fig. 16b shows that the non-consistency in using mass fluxes 2 leads to more kinetic energy loss than even the coarsest-grid simulation that is based on mass fluxes 1.In the remainder, the densities ρ * and ρ are based on fluxes 1 in Eq. ( 14).The droplet's interface deformation is illustrated in Fig. 17.These figures show that the consistent method with a grid of 80×80 preserves the shape of the droplet while the non-consistent does not.The C-labeling that was discussed improves the results.Without the extra label, the shape of the droplet is preserved less well and fluids with different densities are kept separate less well.The kinetic energy losses also increase when no C-labels are used from 0.017% to 0.090% with a grid resolution of 40×40.

Momentum conservation droplet impact on fixed wall
The next step is to look at momentum conservation after an impact.In this case, one-way interaction is tested between a solid wall and a droplet of fluid.The setup is similar to the case with the droplet translating diagonally through the domain, but, now apertures (C b , a b ) that represent the solid are involved.
The setup is illustrated in Fig. 18 where the dashed line indicates the expected deformation of the fluid droplet, which because of the 2D nature of the setup is actually a fluid cylinder.During the impact of the droplet with the wall, all horizontal momentum is converted to vertical momentum as the droplet exerts a force on the wall.The initial speed of the fluid droplet is 10 [m/s] and it has a density of 1000 [kg/m 3 ].The density of the gas surrounding the fluid is kept low at 10 −9 [kg/m 3 ], for the same reasons as before.The radius of the cylindrical droplet measures 0.15 [m] initially.Gravity is equal to zero.The effect of viscosity and surface tension is neglected.After 0.03 [s], the impact finishes when the fluid jets reach the bottom and top of the domain.The droplet is positioned such that the loss of momentum takes place during the impact and not because of translation.The edge of the wall being impacted is positioned along the center of a cell.The initial horizontal momentum I 0 is equal to 706.9 [kgm/s].After the impact, loss of horizontal momentum of the fluid should be equal to the integrated force in time acting on the wall.The following equation is used for the horizontal momentum transfer at any time where F b,x is the force acting on the wall as a function of time t n , V m,h is the volume of the corresponding momentum control volume with indices (i + 1 2 , j), n the time level, and I 0 the initial momentum.The value indicates the fluidstructure interaction error during the transfer.The case is simulated with the consistent (C.) and the non-consistent (NC.) method for three different grid resolutions.The integral in (40) is solved using the midpoint rule.The value as a function of time is plotted in Fig. 19a for different grid resolutions and both the non-consistent and the consistent method.The results show convergence.The error for the non-consistent method is an order higher than for the consistent method.The error of the non-consistent method first goes up, but then, unexpectedly, it goes down.It was found that the error goes down because the free surface started to form protrusions similar to those in Fig. 14c with a velocity  in opposite direction.Negative velocities reduce the error as it is defined in Eq. (38).The maximum errors are given in Fig. 19b.

Two-way fluid-structure interaction
Two-way fluid-structure interaction is evaluated next.A droplet of fluid with a horizontal velocity will impact a wall that is initially at rest but free to move.Special about this case is that the apertures are not only included, but that their size now also depends on time.A similar configuration is used as in the previous section, but now the change in horizontal momentum of the fluid should be equal to the gain of (horizontal) momentum of the structure after impact.
The setup is illustrated in Fig. 20.The speed of the droplet is 10 [m/s] initially, while the rectangular wall is at rest.The droplet has a smaller initial radius than the previous case with the fixed wall to make sure the wall has reached its limit velocity before the fluid jets hit the top and bottom of the domain.The droplet, with a radius of 0. We expect that all the horizontal momentum of the fluid is converted into horizontal momentum of the structure.The data required for the comparison is straightforward to obtain from the simulation results until the jets created by the interaction with the structure hit the bottom and top of the domain.Horizontal momentum conservation is calculated using Fluid cylinder: fluid-structure interaction error due moving wall impact for consistent method (C.) and non-consistent method (NC.)

Table 1
Effect of the size of the droplet radius in the case of the two-way interaction with the moving wall for a grid of 200×200.The radius is a measure of the number of interface cells (S, C) over internal cells (F).The maximum momentum loss during translation (trans.)and after impact are given for consistent (C.) and non-consistent (NC.) simulations. R

Analyzing the momentum losses
The governing error appears to be related to an inconsistency between the auxiliary density field ρ and the new density field ρ n+1 .The discretization of ρ, see Eq. ( 22), is based on fluxes through the boundary of a momentum control volume, while the discretization of ρ n+1 , see Eq. ( 18), is based on volumes in the four sub-volumes of the momentum control volume.The difference between these two density fields is caused by reconstruction of the interface that changes the distribution of the fluids in the four sub-volumes of a momentum cell to become different from what one would expect based on the fluxes of the momentum control volume alone.Therefore, the source of the error is located near the interface.This is investigated by varying the radius R of the droplet in the two-way interaction case with a grid resolution of 200×200.For the circular shape of the droplet, the number of internal cells varies with the radius squared, while the number of interface cells, which is related to the circumference, varies linearly with the radius.The results of the investigation are in Table 1.It shows that the momentum errors reduce with increasing radius of the droplet and, hence, with an increase of the number of F-cells over the number of S-cells.In smaller droplets, therefore, the ratio of cells that contribute to the error versus cells that contribute to the total momentum is less favorable.Note that the errors of the consistent simulations in Table 1 are approximately two orders of magnitude smaller than those in the non-consistent simulations.This is an argument that using the consistent method becomes more relevant for keeping errors in check in even more violent interfacial flows with larger curvature than in the test cases with a droplet impacting the moving and fixed walls.
Potential solutions to decrease these errors further could be omitting the reconstruction of the interface and taking the diffusion of the interface for granted, using a collocated arrangement of solution variables, or using the method of Rudman [34] with interface transport and reconstruction on a grid that is twice as fine as the base grid at the expense of computational effort.

Computational costs
The difference in computational cost between the consistent method and the non-consistent method is not straightforward.Considering the droplet impact on a fixed wall with grid 60×240, the consistent method needed 14% more time steps than the non-consistent method.The increase in computational costs is caused by the discontinuous handling of the density, resulting in a better prediction of the moment of impact with higher flow velocities immediately after impact and, hence, smaller time steps because of the Courant restriction.Similarly, for the two-way fluid-structure interaction with the moving wall and grid 300×300, the consistent method needed 16% more time steps.However, for the diagonally translating high density droplet the number of computational time steps with grid 80×80 is almost 50% higher for the non-consistent scheme with respect to the consistent scheme.The increase is due to spurious velocities at the distorted interface.

Comparison benchmark & experiment
Results of the consistent method (this article) and of the older non-consistent method [12] are compared with two benchmarks and with experimental results.

Rising bubble
The first case is a 2D rising bubble benchmark.Our results, with the consistent method (C., this article) and the nonconsistent (NC.) method [12], are compared with the results of Hysing et al. [16].Literature has shown that non-consistent methods with CLSVOF for the interface [36,35] give good results for rising bubble simulations.Therefore, when referring to 'non-consistent' we are only referring to van der Eijk and Wellens [12].The proof of principle simulations in the previous section mainly evaluated the convective term in the momentum equation.Simulating the rising bubble also requires accurate modeling of capillary, diffusive, and buoyancy effects.New in this article is the 'consistent' discretization of the viscosity near the interface, where other methods typically make use of harmonic averaging for the viscosity [12,33].The focus for the rising bubble case is to evaluate the discretization of the viscosity and the curvature.
The benchmark of Hysing et al. [16] with the 2D rising bubble was created due to the absence of analytical solutions and is based on the average of the results of high grid resolution simulations with several methods.The initial fluid configuration of the 2D rising bubble case is illustrated in Fig. 22a.The fluid variables per specific instance of the case are given in the table in Fig. 22b.Results of different simulations are compared in terms of the velocity of the bubble over time, which is calculated as the mean of the vertical velocities in the cells contained within the contour of the bubble in which V b and v c are the volume of the bubble and the terminal velocity, respectively.The bubble is tracked based on the labels.Only cells with labels C, S and E contribute to the average velocity.The subscripts of the velocities v n and v s indicate the northern and southern velocity in the staggered arrangement.The first instance in the table in Fig. 22b does not have large fluid property ratios.The viscosity and density of the liquid are both only one order of magnitude higher than for the gas bubble.As discussed in earlier sections, a low ratio between fluid densities is expected to lead to small differences between the consistent (C., this article) approach and the non-consistent (NC.) approach [12].The NC. results are therefore omitted from this part of the benchmark.
The base case for simulating the first benchmark is the consistent method for the density (this article), harmonic averaging of the viscosity near the interface, and computation of the curvature with a stencil of 3×3 cells.One by one, the improvements discussed in this article are evaluated, the consistent discretization of the viscosity (Visc.), as described in Sec.4.2.1, and a larger stencil of 5×5 cells for the curvature (Curv.).The simulation results in terms of the velocity of the bubble as a function of time are shown in Fig. 23a.
The maximum velocity reached by the bubble is considered a good error measure for the quantification of the difference between our results and the benchmark.The maximum velocity of the consistent simulation is 1.6% smaller than the bench-  mark.When using the flux-based viscosity the difference with the benchmark reduces to 1.2%.Changing the stencil for the curvature from 3×3 to 5×5 cells yield a difference with the benchmark of 0.8%.And when everything is combined, the difference with the benchmark in maximum bubble velocity at a grid resolution of 50×150 drops to 0.6%.The convergence of our newest consistent implementation for different grid resolutions is shown in Fig. 23b.Even at coarse grid resolutions, the difference with the benchmark can be considered small.The difference with the benchmark in maximum velocity of the bubble at the highest resolution of 80×240 is 0.02%.
The second instance of the benchmark in Fig. 22b has larger fluid property ratios than the first instance.We now expect that the effect of using the consistent approach for the density near the interface is larger.The results of the consistent method (C., this article) and the older non-consistent (NC.) method [12] are compared with the benchmark and with each other to show the difference between the two methods.
Fig. 24a compares the velocity of the bubble with large ratios in fluid properties over time between simulations.The maximum velocity of the bubble in the NC.simulation has a difference of 4.45% with the benchmark.When the nonconsistent method for determining the density near the interface is replaced by the consistent method, the difference in maximum velocity with the benchmark reduces to 0.45%.Both the consistent discretization of the viscosity (instead of harmonic averaging) and increasing the stencil for the curvature from 3×3 to 5×5 cells each have a marginal influence on the results and bring the difference in maximum velocity with the benchmark down to 0.44%, which is also the difference with the benchmark when all improvements are combined.The grid convergence study is shown in Fig. 24b.The consistent method with improvements is in good agreement with the benchmark, already at low grid resolutions.At the highest resolution of 80×240, the difference in maximum bubble velocity with the benchmark is as low as 0.03%.
The use of C-labels has not been discussed for the comparison with the rising bubble benchmark.C-labels did not significantly affect the velocity of the rising bubble, but they did affect its shape.That becomes most apparent for the rising bubble when C-labels are combined with the older non-consistent method [12].In Fig. 25 the bubble shapes are shown when C-labels are combined with the older method [12] at grid resolution 50×150.Half of the rising bubble shape is extracted from the non-consistent simulation results by filtering for fluid fractions above 0.05.It is shown to the left.The half of the bubble to the right shows the results when the non-consistent method is combined with C-labels.The C-labels have made sure that a lot less distortion of the shape of the bubble has occurred and that the simulation result is closer to the shape presented in Hysing et al. [16].

Buoyant cylinder water exit & entry
Like in van der Eijk and Wellens [11], our ultimate interest is in modeling the fluid-structure interaction involved in the motions (accelerations) of high-speed craft in heavy seas that frequently disconnect from the main water body when moving vertically upward before falling downward with an impact.An experiment that captures fundamental aspects of these phenomena was conducted by Colicchio et al. [8].The experiment consists of two parts.The first is a buoyant cylinder that moves vertically upward before exiting the water, the other is that same cylinder moving vertically downward before impacting with the water.The experiment was designed to be as 2D as possible.The tests were conducted in a prismatic tank with a length of 3000 [mm], a water depth of 1400 [mm] and a width of 400 [mm] in the direction that is assumed not to be relevant to the dynamics of the cylinder.The cylinder itself has a radius of 0.15 [m] and an effective density of 620 [kg/m 3 ].During testing, the cylinder's position and velocity were measured, together with pressures at discrete locations along the contour of the cylinder.In their article, the measured quantities were also compared to results of a numerical method with a level-set free surface displacement algorithm.Details of the numerical simulations, such as grid size and time step were not reported.
The initial fluid configuration for our simulations is illustrated in Fig. 26, with specifics provided in the table in Comparing numerical simulations to measurements from experiments is challenging, because you need to make sure that the differences between the two are sufficiently small not to have an effect on the comparison.Two, what we think are, cardinal moments during the experiments are shown in Fig. 27.These are the moments when, for the water exit case, a free surface jet impinges on the main body of the fluid, and when, for the water entry case, the fluid from either side starts flowing over the cylinder and meets in the middle.These are events that can be compared to long-crested wave breaking, that are inherently 3D and sensitive to small perturbations [23], even when the setup of the experiment is otherwise flawless.This will be discussed in more detail.

Cylinder exit
We will first discuss the case with the buoyant cylinder rising upward and exiting the water.Much of the dynamics of the cylinder takes place away from the interface between water and air.For this reason, the difference between results of the non-consistent and the consistent method should remain small.The cylinder's velocity and position during exit are shown in Fig. 28.The measurements (Exp.) are compared with levelset method (Num.) of Colicchio et al. [8], our older non-consistent method (NC.), and the consistent method (C.) discussed in this article with all improvements at different grid resolutions.Those are for N = 75, 125, or 175 cells in both directions.
A dashed line is plotted around 1 [s].From that moment, according to our interpretations, 3D effects will start to play a role, see Fig. 27a, and the results in the simulations will start to diverge from the measurements.The simulation results in Fig. 28, for NC. and C. alike, match well with the experimental results before t = 1 [s], and grid refinement with the consistent method (C.) leads to convergence.
The pressures from a simulation with the consistent method at a resolution of 125 cells in both directions at discrete locations along the cylinder contour are compared with the measurements (Exp.) and the numerical results (Num.) of Colicchio et al. [8] in Fig. 29.The pressures were measured at angles of 0, 15, 45, 105, and 145 • , where 0 • indicates the bottom of the cylinder.During the cylinder exiting the water, the low-pressure wake formed behind it during its upward movement interacts with the air phase and entrains air bubbles.The first moment of entrainment is illustrated in Fig. 27a.In the experiment, the entrainment caused noisy pressure records, which were filtered out before being presented [8].The results of our simulations are not filtered.The smaller spikes are due to the cylinder moving from cell to cell, the larger spikes are due to the object moving from cell to cell near the interface between water and air.The pressures in the simulation show a fair agreement with the experiment.The probe at 105 • becomes dry after 0.82 [s] while 0, 15, and 45 • remain wetted by a mixture of water and air.It is our interpretation that the simulations overestimate the pressure in the water-air mixture later in the simulation, because the breaking-up of air bubbles is not part of the formulation of the method.

Cylinder entry
Next, the case of the cylinder entering the water with an impact is discussed.Moments after the cylinder has penetrated the interface, two fluid jets are formed that develop in the direction opposite to the downward motion of the cylinder.The differences between results of using the non-consistent or the consistent method are expected to be larger than for the cylinder exiting the water, because of the higher fluid velocities and the larger contour of the interface along the jets.In the wake above the cylinder as it is falling down through the interface, an air cavity is generated.At some moment, the air cavity becomes unstable and collapses.The collapse results in the mixing of air and water.After the collapse, the free surfaces at both sides of the top of the cylinder move towards each other and make contact (approximately) in the middle.In the experiment, it is around 0.35 [s] when the cylinder gets fully surrounded with water.This moment is shown in Fig. 27b.Our interpretation again is that the experiment is 2D until the free surfaces merge above the cylinder.Just like wave breaking the merge is not expected to occur uniformly in the third dimension, causing 3D effects from that moment on.
The velocity and position are plotted in Fig. 30.The position data extracted from Colicchio et al. [8] is shifted along the vertical axis so that at time instance 0 [s] the position h b is 0. In both directions, N = 75, 100, or 125 cells were used.
The results obtained with the consistent method (C.) converge, and there is a clear difference with the results of the nonconsistent method (NC.).Until 0.35 [s], the position h b with the consistent method matches well with the experiment.With the non-consistent method, the penetration depth is somewhat smaller than in the experiment.After 0.35 [s], the simulation results start to deviate from the experiment.Especially the velocity with the consistent method undergoes a sudden change when the cylinder becomes completely enclosed in water.We think that the water above the cylinder in the experiment is actually a mixture of water and air with a lower density than water alone, affecting the velocity of the cylinder to a lesser extent.This interpretation is in agreement with the difference between the consistent method and non-consistent method, because in the latter the density near the interface is also more diffuse, comparable to Fig.The pressures from the simulations with 125 cells in both directions, with the consistent method (C.) as well as the nonconsistent method (NC.), are compared to the pressures from the experiment of Colicchio et al. [8] in Fig. 31.For reference, their numerical results are also added to the plots.The pressure sensors were placed at 40 or 45 • intervals along the contour of the cylinder where 0 • coincides with the bottom of the cylinder.The agreement in pressures from the simulation and those in the experiment is fair, good perhaps when the pressure from the sensor at 145 • would have been disregarded.This sensor is placed in what can be called the wake of the cylinder and experiences the most complicated combination of water and air.The sensor at 0 • in Fig. 31a registers the impact pressure when the cylinder hits the free surface.Note that the time axis for the graph of this sensor is different from the time axes of the other sensors in Fig. 31.The pressure peak at 0 • in the simulation with the consistent method reaches approximately 2/3 of the peak in the experiment.From the pressure in the experiment one can tell that more is going on: the pressure oscillations following the impact are likely due to vibration of the container in which the cylinder is dropped or density waves at the speed of sound in water.Those oscillations are absent from the simulations.
The comparison between consistent simulation results and non-consistent simulation results in terms of the pressure does not allow interpretation, except for the pressure sensor at 0 • .Oddly, the simulation with the non-consistent method does not register a peak in the pressure at all, and produces an elevated pressure level some time before the impact takes place in the experiment.For this sensor, it is clear that the non-consistent method does not represent the expected physical behavior, and that the consistent method does.

Conclusion
In this article, the method of van der Eijk and Wellens [12] is combined with the consistent mass-momentum transport (CMOM) scheme of Zuzio et al. [46] and extended with two-way coupled fluid-structure interaction (FSI).With the new method, we found similar results in terms of kinetic energy conservation as Zuzio et al. [46] for the diagonally translating high-density bubble.
In addition to the translating bubble, results of our method are compared with the rising bubble benchmark of Hysing et al. [16], which is highly relevant for this type of two-phase flow.The difference in maximum velocity of the bubble between the consistent implementation and the benchmark was 1.6% when the density ratio was 10, and 0.45% when the density ratio was 1000 with a grid of 50×150 cells in horizontal and vertical direction, respectively.That needs to be compared with a 1.6% (density ratio 10) and 4.5% (density ratio 1000) difference between the older, non-consistent implementation at the same resolution and the benchmark.Several other improvements contributed to further reducing the difference in velocity: the viscosity discretized using VOF fluxes had a small, but measurable influence, and the combination of a larger stencil for the curvature of the interface and C-labels near the free surface made sure that also the shape of the bubble was in agreement with the benchmark and that fluids were kept separate better.With these improvements, the consistent method converges to a 0.03% difference in maximum velocity with the benchmark at a grid of 80×240.
It was demonstrated that using a consistent method in a setting with FSI is advantageous for momentum conservation.Newly devised proof of principle simulations involving a droplet of water impacting a moving solid wall showed an order of magnitude improvement in momentum conservation between the consistent method and the older, non-consistent method.
Simulations with both the consistent and the non-consistent method compared well with the cylinder exit experiment of Colicchio et al. [8] until, in our opinion, the experiment cannot be considered 2D anymore.The comparison with a similar experiment of the same authors, involving a cylinder falling on a free surface that is initially at rest, demonstrates that the non-consistent method gives an elevated pressure before the time the impact takes place in the experiment, but does not reproduce a peak impact pressure.In contrast, the consistent method captures the moment of impact well and produces a peak impact pressure that is approximately two-thirds of the impact pressure in the experiment.
The benefits of consistent transport modeling become especially apparent for two-phase flows with a high-density ratio between fluids.This was observed both for the simulations with FSI and without.Consistent modeling improves momentum conservation and yields a more realistic shape of the interface.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 4 .
Fig. 4. Continuity control volume (-) (i, j) for a cut cell.Body is indicated by .Fluid is indicated by .u b and v b are the body velocities.Center (i, j) is given by .

Fig. 5 .
Fig. 5. Horizontal momentum control volume (-) in uncut cells and in a cut cell with indices (i, j).The lower half of a vertical momentum control volume (-) is also shown.Body is indicated by .Fluid is indicated by .The center of a continuity control volume in this representation is given by .

Fig. 6 .
Fig. 6.Flux calculated in a cut cell representation of a continuity control volume.Body is indicated by .Fluid (liquid) is indicated by .Center (i, j) is given by .The amount of fluid being transported (fluxed) is hatched with (-).

step: 1 . 2 . 1 . 2 . 4 . 3 . 5 .
Solve transport equation fluids C n+1 f and body C n+1 b and reconstruct interfaces 2. Determine the discretization of the density before establishing auxiliary vector field ū Determine ρ n (and ρ n+1 and μ n ) 2.2.Determine ρ * 2.3.Solve temporary continuity equation for ρ Determine auxiliary vector field ū Determine auxiliary vector field ũ with p n and ũb 4. Solve system of Poisson equation for pressure change δ p and equation of motion for body velocity u n+1 b Solve new velocity field u n+1 from the updated pressure differences.

Fig. 7 .
Fig. 7. Reconstruction of cut cells including the labels.Body is indicated by .Fluid is indicated by .
in which A is the area of the cell face open for flow.The free surface of the fraction fields C * X f The new reconstructed fraction fields are used to calculate

Fig. 8 .
Fig. 8. Reconstruction of cut cells: multiple fluid bodies with an intersection at .The labels are included.Body is indicated by .Fluid is indicated by .
in which C f is the volume fraction of the liquid and (1 − C f ) the volume fraction of the gas.C b is the part of a cut cell that is open to flow.ρ f is the density of the liquid, and ρ a the density of the gas.Refer to Fig. 4c for the definitions of the volume fractions.
which A is the area of the face open for flow depending on a b and the size of the face.Note that the mass flux δC f needs to be determined in a similar way to what is shown in Fig. 6, but then through the boundary of a momentum control volume instead of a continuity control volume.The parameters required for determining ρ * along the faces of a horizontal momentum control volume -velocities, apertures and mass fluxes -are shown in Figs.10a through 10c.

Fig. 10 .
Fig. 10.Horizontal momentum control volume with velocities (a), face apertures (b) and mass fluxes (c) along the faces of the control volume necessary for the consistent discretization of the density ρ * .

Fig. 11 .
Fig. 11.Comparison of the effect of using Eq.(26) and Eq.(27) on the free surface deformation.Black indicates volume fraction C f equal to one.White indicates volume fraction C f equal to zero.

Fig. 12 .
Fig.12.Stencil for curvature κ in case of a horizontally oriented interface using continuity control volumes.Fluid is indicated by .The considered continuity control volume is given by (-).

Fig. 13 .
Fig. 13.Cell face value of κ when using PLIC.The cells are labeled as E-cell and S-cell.Fluid is indicated by .Horizontal momentum control volume is illustrated by (-).

0 and M n+1 b
depend on the time-varying volume apertures C n+1 b and face apertures a n+1 b

using
Crank-Nicolson time integration.Here, m b represents the body's mass and A n+1 f the operator that integrates the pressure over the surface of the body.The operator A n+1 f makes use of the cell face apertures a n+1 b

Fig. 16 .
Fig. 16.Diagonally translating droplet: kinetic energy as a function of time.

Fig. 20 .
Fig. 20.Two-way fluid-structure interaction: setup of simulation domain with fluid droplet impacting a wall that is free to move.
10 [m] is positioned 0.06 [m] in front of the rectangular wall at the start of the simulation.The rectangular wall has a height of 0.9 [m] and a width of 0.19 [m], centered at 0.655 [m].The impact takes place after around 0.006 [s].The gravity constant is set to zero.The density of the fluid and the structure both are 1000 [kg/m 3 ].The density of the gas around the droplet and the wall is equal to 10 −9 [kg/m 3 ].The Courant restriction equals 0.2.

Fig. 23 .
Fig. 23.2D Rising bubble: results of the first instance of the benchmark with low fluid property ratios.

Fig. 24 .
Fig. 24.2D Rising bubble: results of the second instance of the benchmark with high fluid property ratios.
Fig. 26b.The value h b indicates the center position of the structure, and h b,0 the initial position.The following fluid properties are used: densities ρ a = 1 [kg/m 3 ] and ρ l = 1000 [kg/m 3 ], viscosities μ a = 1.48 • 10 −5 [kg/ms] and μ l = 10 −6 [kg/ms], gravitational constant g = 9.81 [m/s 2 ] and surface tension coefficient σ = 0.072 [N/m].The simulations are run until 1.5 [s] with a maximum Courant number of 0.2.The initial speed for the cylinder exiting the water is 0 [m/s]; for the cylinder falling downward before entering the water, it is 2.55 [m/s].

Fig. 29 .
Fig. 29.Pressure for buoyant cylinder exiting the water with N = 125 cells in each direction (Exp.and Num.results taken from Colicchio et al. [8]).0 • indicates the bottom of the cylinder.

Fig. 31 .
Fig. 31.Pressure for cylinder entering the water with N = 125 cells in each direction (Exp.and Num.taken from Colicchio et al. [8]; note that the time axis is different in Fig. 31a. .