Enforcing accurate volume conservation in VOF‐based long‐term simulations of turbulent bubble‐laden flows on coarse grids

This study proposes two different strategies to enforce accurate volume conservation in volume‐of‐fluid (VOF)‐based simulations of turbulent bubble‐laden flows on coarse grids. It is demonstrated that, without a correction, minimal volume errors on a time‐step level, caused by the under‐resolution of the interface, can accumulate to significant deviations from the intended flow conditions despite the comparably good volume conservation properties of the geometric VOF method. In particular, large volume errors are observed for challenging setups combining coarse grid resolutions and comparably high Reynolds and Eötvös numbers. The problem is reinforced for long‐term simulations in periodic domains, which are often performed to collect flow statistics of bubbly flows. The first proposed volume conservation method simply corrects the volume error of a bubble by uniformly adding or removing the respective amount of gas volume in the interface cells. The second proposed method performs an additional reconstruction and advection step of the VOF field using a non‐divergence‐free velocity field, which can be interpreted as a slight dilatation or contraction of the bubble. A comparison between the global flow statistics as well as the individual bubble dynamics for both volume conservation methods reveals that the results are quasi‐identical for a number of challenging test cases, while the gas volume is accurately conserved. The proposed methods allow to perform numerical simulations of freely deformable bubbles in turbulent flows for setups that have previously been out of reach for this numerical framework.

such as non-intrusive imaging are usually limited to very dilute void fractions, 1,2 numerical investigations play a major role in improving the fundamental understanding of such flows.4][5][6][7][8] While the pioneering studies of Lu and Tryggvason [9][10][11] featured relatively small-scale channel flows with up to 72 bubbles, recently developed massively parallel solvers 5,12 are able to perform DNS with up to  (10 4 ) bubbles in full-scale channels, while maintaining excellent scalability.However, due to the computational cost, DNS will remain limited to reduced-complexity bubbly flows at moderate Reynolds numbers for the foreseeable future.On the other hand, simplified numerical approaches representing bubbles as Lagrangian point particles instead of resolving them are unsuitable for many flow configurations.This is particularly the case for scenarios in which bubbles assume complex shapes, or in which resolving the gas-liquid interface as well as the internal flow field of the bubbles plays an important role.One example in this regard is the process of pool scrubbing, where air contaminated with particles is injected into a water pool in order to transfer the particles from the air bubbles to the surrounding water.Consequently, the combination of large eddy simulation (LES) and an interface capturing method is a promising approach for the design of technical devices. 13Compared to DNS, much less of the interface dynamics can be resolved in LES.This leads to challenges regarding the numerical algorithms for two-phase flows, such as the well-known problem of imprecise curvature computation.Another particularly important issue in this regard is the violation of volume conservation, which is the subject of this study.In general, the volume-of-fluid (VOF) 14 method, which is used in the present work, is known to exhibit significantly better volume conservation properties than, for example, the Level-Set 15 or Front-Tracking 16 methods in their original forms.Nevertheless, as will be demonstrated, small volume errors on a time step level, caused by the under-resolution of the interface, accumulate to severe errors if no correction procedure is carried out.
Recently, Takeuchi and Tryggvason 17 have discussed this issue for the Front-Tracking method and proposed a respective volume conservation approach, which has been validated by simulating a single rising bubble.In the context of VOF-based studies, one of the most popular methods to correct volume errors stemming from over-and undershoots of the VOF variable is to use local redistribution algorithms, such as the ones proposed by Harvie and Fletcher 18 or Baraldi et al. 19 Their general principle is to adjust the VOF value in a selected cell adjacent to the one where the VOF value must be clipped back to the interval [0, 1], such that the volume error due to the clipping is rectified.Another option is to use the VOF advection procedure by Weymouth and Yue, 20 which per construction avoids over-and undershoots and thus conserves volume to machine precision on cartesian grids.However, as will be demonstrated and discussed in Section 3.2, other sources for volume errors can be present besides over-and undershoots of VOF values, which requires alternative approaches to enforce volume conservation.
In many well-known VOF test cases, such as a single rising bubble, Zalesak's disk, the phase inversion benchmark 21,22 or bubble-induced turbulence, 23 volume errors never reach a problematic level due to the short transients that are investigated.However, volume conservation errors have been observed and measured in some previous VOF studies.A very recent example is the work of Di Giorgio et al., 22 where the well-known phase inversion benchmark 21 has been studied for both geometric as well as algebraic VOF, as well as for different momentum advection schemes and grid resolutions.While, due to the short transient investigated, the volume errors are very low, it is clearly demonstrated that deviations are already present at the end of the transient.Significantly larger volume errors have been observed in the work of Weymouth and Yue, 20 where a high-energy breaking air-water wave is investigated using different geometric VOF methods.While the simulation performed with the VOF algorithm proposed by the authors accurately conserves volume, the cases which employ alternative VOF schemes already lead to a total volume change of up to 12.63% after one wave period.
For this work, the frequently studied test case of a periodic channel flow 24 is extended by including a significant number of freely deformable gas bubbles.The resulting bubble-laden downflow channel is a well-established and frequently studied simulation setup. 5,6,9,10To demonstrate the problem and propose solutions to it, the presented bubbly flows are purposefully set up in a way that, compared to previous numerical studies of turbulent bubbly flows, results in very challenging conditions for numerical simulation.In particular, the setups stand out for the following specifications: 1.The grid resolution is significantly coarser than what would be required for a two-phase flow DNS, and rather targeted at typical LES-like resolutions.2. The Reynolds and Eötvös numbers are relatively high compared to previous studies of bubble-laden channel flows.
This also leads to significant deformations of the bubbles.3. A large number of flow-through times (i.e., a long simulation interval) is computed to collect converged statistics.This is especially important for statistical quantities involving the gas volume fraction.
This work proposes two different methods to enforce accurate volume conservation for VOF-based simulations of such challenging flow configurations.Compared to previous approaches, their main advantage is that volume is conserved irrespective of the numerical origin of the volume error.In order to validate the proposed methods, they are tested for a variety of different setups.First, the Eötvös number is varied to demonstrate that the methods are not limited to nearly spherical shapes, but can also handle ellipsoidal and wobbling bubbles.Moreover, it is shown that the methods are reliable for different grid resolutions, which is of great importance for efficient industry-scale simulation.The validation is not limited to global flow statistics, such as the mean velocity or turbulent kinetic energy profiles of the channel flow.Instead, also the individual bubbly dynamics are thoroughly investigated in order to ensure that the interface dynamics are not negatively affected by the volume conservation methods.This is important for certain application scenarios involving turbulent bubbly flows, such as the previously mentioned process of pool scrubbing.By using the proposed volume conservation methods combined with the multiple marker formulation, 25 the simulation of freely deformable bubbles in turbulent flows becomes feasible for setups that have previously been out of reach for the geometric VOF method.The methods can easily be be added to a VOF-based solver, and can thus contribute to the further development of reliable and efficient CFD codes for predicting turbulent bubble-laden flows.
The rest of the article is organized as follows: Section 2 introduces the governing equations, the numerical methods and the simulation setup.Section 3 demonstrates the problem of volume conservation for different test cases and discusses the causes for the observed volume errors.Subsequently, the two volume conservation methods are introduced.The resulting flow and bubble statistics as well as the individual bubble dynamics for both methods are presented in Section 4. Finally, Section 5 summarizes the main findings and draws conclusions.

NUMERICAL METHOD AND FLOW CONFIGURATION
The simulations for this study have been performed using the state-of-the-art open source software "TBFsolver," 5 which is based on the finite volume method.Section 2.1 introduces the governing equations and their numerical implementation.Section 2.2 presents the flow configuration.

Governing equations and numerical framework
The mathematical model describing the immiscible gas-liquid flow uses the one-fluid formulation 26 of the incompressible Navier-Stokes equations.A single set of equations is solved for the two phases and a geometric VOF method 14,27 comprising a piecewise-linear interface calculation (PLIC) accounts for the interfacial terms, as well as the local densities and viscosities.The latter are linearly interpolated from the material properties of the liquid and gas phase on the basis of the local gas volume fraction.To avoid the unphysical merging of colliding bubbles, a phenomenon referred to as numerical coalescence, a multiple-marker formulation 25 is used.It transports an individual VOF marker function f b for each bubble b present in the domain: The split-direction advection method used to solve Equation ( 1) is investigated in more detail in the context of the volume conservation strategies.Thus, further details regarding the two implemented advection procedures by Puckett et al. 28 and Weymouth and Yue 20 are provided in Section 3.1.
In the usual notation (velocity components u i , pressure p, density , average density in the domain  0 , dynamic viscosity , gravitational acceleration g i ), the mass and momentum equation are given as Following the Continuous-Surface-Force method, 29 the surface tension force is determined from the surface tension coefficient , the liquid-phase oriented unit vector normal to the interface n i = (f b ∕x i )∕|∇f b |, the Dirac interface indicator function  S = |∇f b |, and the interface curvature  = (n i ∕x i ), where  is computed with a generalized height function approach. 30he Poisson pressure equation is transformed into a constant-coefficient formulation, as this allows the application of fast Poisson solvers. 31,32The introduced equations are solved with the projection method.Time integration is performed using a fully-explicit second-order Adams-Bashforth scheme, while the spatial discretization of the momentum equation is realized using "Quadratic Upstream Interpolation for Convective Kinematics" (QUICK) 33 for the convective term and central differencing for the diffusive term.A staggered grid arrangement is implemented, that is, the velocity components are stored at the cell faces, whereas the pressure, the gas volume fraction as well as the material properties are stored in the cell centers.

Case definition
Figure 1 illustrates the vertical downflow channel investigated in this study.The stream-wise (x), wall-normal (y) and span-wise (z) dimensions are L x = 8H, L y = 2H and L z = 4H, respectively.The channel half-width H is taken to be unity.No-slip walls are prescribed in the y-direction, whereas the x-and z-direction are periodic.A constant pressure gradient corresponding to a friction Reynolds number of Re  = (u  H)∕ l = 590 drives the flow downward.
Here, the friction velocity is defined as u  = √  w ∕ l .Moreover,  w denotes the average wall shear stress in the statistically steady state, and  l and  l are the density and the kinematic viscosity of the liquid phase, respectively.To initialize the simulation, 780 initially spherical bubbles with a diameter of d = 0.25H are inserted, which corresponds to a technically relevant gas volume fraction of 10%.Due to the implemented numerical method, the bubbles are freely deformable.The density and the dynamic viscosity ratio are both specified as  l ∕ g =  l ∕ g = 20, where the subscripts l and g denote the liquid and gas phase, respectively.Consequently, the bubbles rise relative to the liquid carrier phase.The gravitational acceleration is selected such that the Galilei number is Ga =  l d √ gd∕ l = 417.19.Finally, the time step is adjusted such that the CFL number is 0.2.All specifications provided up to this point are identical for all setups simulated in this work.The next sections will discuss several variations that are performed to investigate the influence of different parameters on the volume error and the proposed volume conservation methods.
To study the influence of bubble deformability, two different Eötvös numbers are considered in this work.The lower value is specified as Eo = ( l gd 2 )∕ = 0.67.In combination with the resulting bubble Reynolds number, which is Re b = (u rel d)∕ l ≈ 478, this leads to nearly spherical, mildly wobbling bubbles. 34The average relative velocity between the two phases in stream-wise direction is defined as u rel = ⟨u⟩ l − ⟨u⟩ g , where the phase-conditional averaging ⟨⋅⟩ l,g is limited to the bulk region 0.5 < y∕H < 1.5 to exclude the relatively bubble-free wall boundary layers. 7The higher Eötvös number of Eo = 3.75 is prescribed by decreasing the surface tension coefficient, while retaining all other parameters.Here, the bubble Reynolds number is approximately Re b ≈ 358, which in combination with Eo leads to ellipsoidal and wobbling bubbles.The bubble shapes are visualized in Figure 1, and a quantitative analysis of the sphericity and the strength of the shape oscillations is presented in Section 4.2.
In addition, the effect of the bubble resolution on the volume error and the outcome of the volume conservation approaches is investigated.For both resolutions, the domain is discretized using an equidistant cubic mesh.The cell size is varied such that the bubbles are either resolved with d∕Δ c = 12.5 or d∕Δ f = 25 cells per bubble diameter.The coarse resolution (subscript c) leads to a total of 8 million cells, whereas 64 million cells are used for the fine resolution (subscript f ).
Besides the different Eötvös numbers and mesh resolutions, an additional variation is introduced by the two distinct volume conservation methods.Finally, as already mentioned in Section 2.1, the VOF advection scheme is varied between the ones of Puckett et al. 28 and Weymouth and Yue 20 for some selected cases.Since the VOF advection schemes and the volume conservation methods are introduced in Sections 3.1 and 3.3, respectively, more detailed overviews of the simulated cases (see Tables 1 and 2) are presented in the respective sections.

VOLUME CONSERVATION: BACKGROUND AND METHODS
For the implemented numerical framework combining the multiple-marker formulation and a directionally split geometric VOF advection method, the three identified causes for volume errors of the tracked gas phase are: 1. Depending on the selected algorithm, the split-direction advection can lead to over-(f b > 1) and undershoots (f b < 0) of local cell VOF values.These VOF values must then be clipped to f b = 1 or f b = 0, respectively.The first case therefore leads to a decrease of gas volume, whereas the second case leads to an increase of gas volume.This is a well known issue 5,27 and it will be investigated in further detail in Section 3.1.As will be discussed there, this problem occurs when the VOF field is advected with the method by Puckett et al., 28 however it is avoided when using the approach of Weymouth and Yue. 202. To avoid smearing of the gas-liquid interface, the VOF values in each bubble's interface cells (0 < f b < 1) are reset to zero if no full gas cell (f b = 1) is detected within a reasonably wide stencil of cells surrounding the considered interface cell.In the present framework, the respective stencil centered at the interface cell consists of 5 × 5 × 5 cells, as shown in Figure 2.This procedure can only lead to a decrease of gas volume.It is important to note that comparable procedures are also used in other state-of-the-art geometric VOF codes, 35 and therefore this issue is not specific for the chosen solver.The necessity of this approach in the context of the geometric VOF method is demonstrated in Appendix A. In principle, the most useful choice of stencil size as well as frequency with which this algorithm is performed can vary to some degree depending on the investigated flow setup.This will influence the corresponding loss of gas volume.However, as long as the algorithm is used, volume conservation violations are unavoidable and must be corrected.It should be noted that, while having well-documented disadvantages compared to PLIC-based VOF methods, 27,36 algebraic VOF methods, such as the ones based on the "tangent of hyperbola for interface capturing" (THINC) scheme, 36 might not suffer from this specific issue.3. Small gas satellites can break off a bubble in regions where the shear forces overcome the restoring surface tension forces.In particular, this issue is observed for coarse grid resolutions, since a certain amount of surface tension forces cannot be resolved any more. 37These small satellites are deleted, leading to a decrease of gas volume.It is important to note that, unlike the smeared interface regions mentioned in the previous point, these satellites have completely separated from the main bubble.This difference is illustrated in Figure 2. Depending on the specific numerical approach, other origins such as slight inaccuracies in the computation of volume fluxes or the reconstruction of the interface could also lead to volume errors.However, the detailed analysis in Section 3.2 has revealed that, compared to the three issues listed above, these are negligible for the present study.Since the volume conservation approaches introduced in Section 3.3 correct volume errors independent of their origin, any additional volume deviations arising from other issues would be corrected automatically.Moreover, it is important to note that using the multiple marker formulation, which assigns an individual VOF function to each bubble, does not cause any additional volume errors compared to tracking all bubbles with one common VOF field.
In the following, this section first introduces the two selected advection methods for the gas volume fraction (see Section 3.1).Subsequently, the order of magnitude of the volume conservation errors is investigated, and the relative importance of the three origins listed above is analyzed (see Section 3.2).Finally, two different volume conservation methods are derived and discussed (see Section 3.3).

Advection of the gas volume fraction
The geometric VOF method used in this study 27 to solve Equation (1) consists of two separate sub-steps.In the first step, the interface is reconstructed using a state-of-the-art PLIC approach.The reconstruction procedure, which is shortly summarized in the following, is not modified in this work and is therefore identical for all cases.Using the PLIC approach, the interface plane in a grid cell is defined as where x pos is its position vector and q denotes the volume fraction enclosed under the cutting plane.A geometrical interpretation of Equation ( 4) is that the plane is moved along its normal vector n (see Section 2.1) until q corresponds to the volume prescribed by the VOF value in the local cell.For further details, the reader is referred to Cifani et al. 27 The second step of the geometric VOF method is the advection of the VOF field.As already mentioned in the beginning of Section 3, two different schemes are tested.Both are based on an operator-split approach, that is, the transport is carried out sequentially for the three dimensions.To avoid asymmetries, the sequence is permutated on a time step level.For simplification, and since the actual sequence is irrelevant, the following explanation assumes that the sequential transport is carried out in the order x → y → z.The first method, which is also the standard approach in the solver, is the one of Puckett et al. 28 As presented by Cifani et al., 27 the discretized three-dimensional algorithm transporting a VOF field f in ) . ( For better clarity, the velocity components are denoted as u, v, and w, and the spatial directions are x, y, and z.In Equation ( 5), f * , f * * , and f * * * are the VOF fields after each one-dimensional update, respectively.The second expression on the right hand side of the first three lines denotes the actual geometric fluxes in the respective direction, whereas the last term is a correction term that limits potential under-and overshoots at each sub-step of the algorithm.This is necessary, since the single sweeps are performed using one-dimensional velocity fields, which are not divergence-free. 20In the last line, these correction terms are subtracted again.It is noteworthy that, for lines one to three of Equation ( 5), the respective updated VOF field on the left hand side is also present in the correction terms on the right hand side.Therefore, for implementation purposes, Equation ( 5) can be written as As discussed in previous studies, 27,39 the algorithm presented in Equation ( 5) does not perfectly conserve volume, as over-and undershoots are still possible in the separate sub-steps of the advection process, and a respective clipping of the VOF values is necessary to ensure 0 ≤ f ≤ 1.This issue was identified as one of the causes for volume errors listed in the beginning of Section 3. Therefore, an alternative advection methodology is considered in the following.
Weymouth and Yue 20 name three requirements for constructing an advection scheme that conserves volume to machine precision.First, the flux terms must be computed in a conservative manner.Second, the velocity fields must be numerically divergence-free.The third requirement is that over-and undershoots are avoided per construction, that is, 0 ≤ f ≤ 1 is not violated at any stage, and thus cutting off too large and too small values is not necessary.Unlike the latter condition, the first and second requirement are satisfied for the implementation of the method of Puckett et al. 28 in the flow solver. 5The flux terms are computed conservatively, and, due to the way the velocity and pressure field are corrected in the solver, 5 the resulting velocity field used to transport the VOF field is numerically divergence free.Using the same notation as in Equation ( 5), the method of Weymouth and Yue, 20 which is applicable on Cartesian grids and additionally satisfies the third requirement, is given as ) . (7) Here, the following binary definition of f bin is used: This binary expression, which remains constant during all sub-steps of Equation ( 7), replaces the VOF values f * , f * * , and f * * * as an estimate of f in the correction term.Since, as discussed above, ∇ ⋅ u = 0 is satisfied to machine precision for the present solver, the last line of Equation ( 7) can be rewritten as f n+1 = f * * * and the sum of the correction terms vanishes automatically.In addition, precise volume conservation is only ensured if the following Courant condition is satisfied: For further details, the reader is referred to the article of Weymouth and Yue. 20

Investigation of the volume errors
This section provides a detailed investigation of the volume errors for a number of test cases.Previous studies with the same solver have demonstrated that, for less challenging flow conditions than the ones investigated in this work, the volume errors are negligible even for long simulation intervals.In their work, Cifani et al. 5 have observed a final volume conservation error of only ∼0.04% for their DNS conducted at Re  ≈ 154, Eo = 3.125 and a resolution of d∕Δ x = d∕Δ z = 15.3 and d∕Δ y = 20, respectively.The setups investigated in this section are listed in Table 1.In order to perform a meaningful analysis that is independent of the initial conditions and allows a comparison between the single cases, the investigation of the volume errors has been conducted in the statistically steady state.For this purpose, one of the volume conservation methods ("dilatation and contraction method") presented in Section 3.3 has been used to conserve volume accurately until a statistically steady state has been reached.In other words, all six cases listed in Table 1 are initialized from equivalent flow conditions, where no volume error is present at the initialization.
Figure 3 shows the global relative volume error of the tracked gas phase  rel , given in percent, for the six setups (see Table 1) over the dimensionless time.In this context, the relative error is defined as  rel (t) = (V g (t) − V g,0 )∕V g,0 , where V g (t) and V g,0 denote the current and the initial total gas volume, respectively.As the figure demonstrates, a decrease of gas volume is observed for all simulations, and the relative volume error can accumulate to a severe deviation from the intended flow conditions.The error is more pronounced for the coarser grid resolution, as well as for the more deformable bubbles at Eo = 3.75.The volume advection method of Weymouth and Yue 20 is only tested for the coarse grid simulations.Despite its ability to avoid over-and undershoots and therefore a clipping of the gas volume fraction at any stage (see Section 3.1), its utilization does neither lead to satisfying volume conservation, nor to consistently improved volume conservation compared to the scheme of Puckett et al.In fact, using the method of Puckett et al. leads to a slightly smaller gas volume reduction at Eo = 0.67, whereas the scheme of Weymouth and Yue shows a clear improvement at Eo = 3.75.Consequently, solely using a VOF advection scheme that avoids over-and undershoots of VOF values does not solve the problem, since the other two issues listed in the beginning of Section 3 remain unsolved.

TA B L E 1
Table 1 additionally provides an overview of the relative importance of the three volume conservation issues listed in the beginning of Section 3. The final (t * = 75, see Figure 3) relative volume error for the gas phase is broken down into the contributions of the three error sources, which have been accumulated individually for the time interval 0 ≤ t * ≤ 75.When offsetting the errors against each other, the final (t * = 75) total errors from Figure 3 are retrieved.As the results demonstrate, the VOF advection scheme of Weymouth and Yue 20 indeed avoids over-and undershoots and a corresponding cut-off of the volume fraction values.It is important to note that the errors listed in Table 1 exceed machine accuracy, since the net error accumulated for around 5 × 10 5 time steps across all bubbles is provided.The scheme of Puckett et al., 28 on the other hand, causes a net increase of gas volume due to over-and undershoots for all setups investigated.For the simulation on the coarse grid (d∕Δ c = 12.5) and Eo = 3.75, a significant increase of 3.74% during the investigated time interval is observed.The by far largest volume errors are, however, generated by the algorithm removing smeared interface regions, illustrated in the top row of Figure 2. Apart from the respective error for the simulation at Eo = 0.67 using the fine grid (d∕Δ f = 25), all observed errors for this category are above 10%.For the most challenging simulation (d∕Δ c = 12.5 and Eo = 3.75), around 50% of the gas volume vanishes due to this error source.On the other hand, the removal of small gas satellites that have completely broken off the respective main bubble (i.e., error 3, illustrated in the bottom row of Figure 2) does not contribute significantly to the volume errors.The strong interplay between the removal of smeared interface regions and the formation of broken-off satellites is demonstrated in detail in Appendix A.
The results presented in Table 1 also explain why the global relative gas volume error for the simulation at Eo = 0.67 and a resolution of d∕Δ c = 12.5 is larger for the Weymouth and Yue advection method compared to the scheme of Puckett et al. (see left plot of Figure 3), despite the fact that over-and undershoots are avoided for the former approach: due to the net increase of gas volume generated by the clipped over-and undershoots of the VOF values for the scheme of Puckett et al., the large decrease of gas volume caused by the prevention of interface smearing is mitigated.The results clearly demonstrate the need for explicit volume conservation methods, such as the ones presented in Section 3.3.The effects of the volume errors on the global flow statistics will be further discussed in Section 4.1.

Volume conservation methods
This section presents two different methods to enforce accurate volume conservation for turbulent bubble-laden flows.
The individual problem of over-and undershoots of VOF values (see Problem 1), which requires a clipping to 0 ≤ f ≤ 1, is well documented in the literature, and different approaches are available to solve it.Some state-of-the-art solvers such as the "Parallel Robust Interface Simulator (PARIS)" 35 or "Basilisk" 40 include the VOF advection method of Weymouth and Yue 20 to avoid the problem from the outset.As discussed by Cifani et al., 5 an alternative solution to mitigate the problem is to apply a redistribution algorithm that locally redistributes over-and undershoots to the neighboring cells.Examples can be found in the works of Harvie and Fletcher, 18 Kwakkel et al., 41 and Baraldi et al. 19 However, as demonstrated in Section 3.2, solely solving the issue of over-and undershoots of the VOF values is not sufficient for the challenging bubbly flow setups investigated in this study.Therefore, to correct the volume error originating from all three error sources listed in the beginning of Section 3, additional correction steps would be required in addition to a redistribution approach.Consequently, no local volume redistribution algorithm is included in the proposed procedure.
The two methods introduced in the following have several common features.First, they are both based on the idea that the volume errors originating from all three error sources are rectified in a single correction step.Second, since every bubble is tracked using its own VOF marker function (see Equation ( 1)), the volume correction is performed for each bubble individually.Third, it is important that the volume conservation methods do not exert a notable influence on the evolution of the flow and bubble behavior.In other words, their only influence on the simulation should be to avoid volume errors.Since it has been observed that the volume error is negligible for a single time step and only accumulates for long simulation intervals, each bubble's volume is corrected in every single time step.Thus, the influence of the correction approaches on the predicted physics can be considered vanishingly small.This last assumption will be reviewed in Section 4.
The first method ("method A") is referred to as the "simple" method in the following, as it uses a straightforward approach to correct the volume error ΔV n b for bubble b in time step n.In a first step, the N interface cells of the bubble are identified.These are defined as cells with a VOF value of  ≤ f b ≤ (1 − ), where  is a tolerance value.Then, the volume error is corrected by modifying the VOF values in these interface cells such that the correct bubble volume is retained.This is achieved by adding the value −ΔV n b ∕(N ⋅ V cell ), which can take positive and negative values, to the VOF values of all interface cells, that is, the VOF correction in all interface cells of one bubble is identical.This approach does not rely on kinematic assumptions regarding the deformation of bubbles and is easy to implement.The default tolerance value  used in this work is  = 1 × 10 −4 .If a very small value is selected, volume might be added to (subtracted from) a large number of cells that are nearly full (empty).This can cause further overshoots (undershoots) of VOF values, which need to be cut-off again.If rather large tolerance values are selected, this problem is avoided, however the correction volume might be unevenly distributed to the interface region.Thus, the "simple" method is not parameter-free.However, as demonstrated in Appendix B, a comparison of the two tolerance values 1 × 10 −4 and 1 × 10 −2 has revealed quasi-identical results.Therefore, it can be assumed that the method is insensitive to  as long as the tolerance value is in a reasonable range.Since each bubble's volume error is very small on a single time step basis and the volume correction is performed in every single time step, the amount of gas volume added or removed in an interface cell in a given time step is very low.Thus, no issues regarding numerical stability have occurred using this approach, and no additional smoothing is required.In theory, the "simple" method could also be used in combination with an adaptive mesh refinement (AMR) functionality, which is, for example, available in the "Basilisk" code. 40In the context of AMR, however, a more general formulation would be required, since the volumes of the interface cells are not necessarily identical.
The second method ("method B") uses a more sophisticated, parameter-free approach.The general transport equation for the VOF field can be rewritten as The integral form of Equation ( 11) is Using the fact that the geometrical fluxes, represented by the second term, vanish at the boundaries of each bubble sub-domain, Equation ( 12) further simplifies to Finally, Equation ( 13) is considered for a single bubble b, which leads to Equation ( 14) describes the evolution of a single bubble's volume, which demonstrates that its volume can be manipulated by imposing an artificial, non-divergence-free velocity field u * .Assuming ∇ ⋅ u * =  = const., a discrete version of Equation ( 14) in any given time step reads where V err.b is the current erroneous bubble volume after the transport of the VOF field and the removal of smeared interface regions and satellites.Moreover, V corr.b is the correct bubble volume to be conserved, and the divergence  is a constant which is valid only in the considered time step and for the investigated bubble.Since all other quantities are known, Equation ( 15) can be solved for .This allows to construct a fictitious velocity field u * with ∇ ⋅ u * = , which is determined as In this context, (x 0 , y 0 , z 0 ) is the center of the bubble.Finally, an additional reconstruction and advection step (see Section 3.1) for the bubble's VOF field is performed, in which the VOF field is transported with the non-divergence-free velocity field u * .If the volume of the bubble is too small before the correction (i.e., V err.b < V corr.

b
), this corresponds to a slight dilatation of the bubble.If, on the other hand, the bubble's volume exceeds the correct volume (i.e., V err.b > V corr.b ), a slight contraction is performed.Therefore, the method will also be referred to as the "dilatation and contraction method" in the following.Figure 4 provides a geometrical interpretation of this procedure.Due to the fact that the geometrical fluxes vanish, the correction procedure in time step n + 1 effectively leads to a correct final bubble volume of As Equation ( 17) demonstrates, the time step width Δt used for the additional reconstruction and advection procedure has no influence: a higher (lower) value of Δt in the correction procedure leads to a lower (higher) value of  (see Equation (15)) and thus to reduced (increased) imposed velocities (see Equation ( 16)) for the additional advection step.However, this is counterbalanced by the fact that the advection procedure in the correction step is performed with a higher (lower) time step value.For this work, the time step is simply taken to be the same as in the regular VOF transport.Finally, it is important to note that the real velocity field u in the framework of the momentum equation remains unchanged, as u * is solely used for performing the volume conservation procedure.In summary, the proposed methods are both based on correcting the volume error of each individual bubble with a single correction step in each time step.Both approaches have limitations regarding the maximum volume error that can be compensated.However, running the respective correction algorithm in each time step ensures that the volume that is added or removed is negligibly small.Therefore, as Section 4 will demonstrate, the volume conservation methods do not exert any notable influence on the flow physics of the present setup (see Section 2.2).The following investigation demonstrates the applicability of the introduced approaches to bubbly channel flows without breakup or coalescence.This flow setup is relatively simple, yet represents an important physical process in a variety of industrial applications such as reactor cooling systems, heat exchangers in power plants, chemical reactors in the process industry or air lift pumps, 7,42 and thus is of great interest.Therefore, in the recent years, first large-scale simulations of freely deformable bubbles in turbulent flow have been performed.Compared to previous studies, [5][6][7]42 the utilization of the volume conservation methods allows to study particularly challenging setups combining high Reynolds and Eötvös numbers and coarse grid resolutions, which have been out of reach for this numerical framework previously. Ths is a first step towards LES of turbulent bubble-laden flows in realistic industrial applications.The two proposed volume conservation methods can both be implemented in any finite-volume-based geometric VOF solver.Due to its low computational cost (see Section 4.3), simple formulation and applicability to arbitrarily complex VOF structures, this holds especially true for the "simple method."While first investigations for setups such as liquid-liquid emulsions show promising results in this regard, a detailed analysis for this and other complex test cases, such as the phase inversion benchmark 21 or Rayleigh-Plateau instability, is beyond the scope of this work.

RESULTS AND DISCUSSION
This section discusses the results of the simulations performed with the two volume conservation methods introduced in Section 3.3.First and foremost, the property of volume conservation, which can be evaluated individually for each simulation, is analyzed.Furthermore, an emphasis is put on ensuring that, apart from accurately conserving the gas volume, the volume conservation methods do not exhibit a notable influence on the flow and bubble behavior.For this purpose, the two methods are validated against each other.It is assumed that, if both methods lead to identical results despite their different approaches, their overall influence on the flow behavior is negligible.While it would be desirable to perform a reference DNS, such a simulation is currently not possible for the challenging parameter range investigated in this work.More importantly, it would not allow a validation of the volume conservation approach, since the results of a coarse grained simulation (e.g., LES) will strongly depend on the mesh, the numerical scheme and several unclosed terms interacting with each other. 13Thus, the most useful practical approach the authors can think of for testing one volume conservation method without a reference case using the alternative method is to make sure that the volume error per bubble and time step is very low, which can be verified without such a reference case.If this condition is satisfied, the proposed methods are expected to accurately conserve the gas volume without exhibiting a notable influence on the flow and bubble behavior.Section 4.1 investigates the global flow statistics, whereas individual bubble dynamics are studied in Section 4.2.Finally, Section 4.3 provides information regarding computational cost of both approaches.
Table 2 gives an overview of the twelve setups investigated in this section.They correspond to the six cases previously simulated without volume conservation method (see Section 3.2), however each simulation is now repeated with both correction approaches.
Figure 5 illustrates the magnitude of the global relative gas phase volume error  rel (t) = |V g (t) − V g,0 |∕V g,0 in percent for three simulations with different volume conservation settings.The investigated time interval is 0.04 (H∕u  ), which corresponds to around 300 time steps.The three simulations have been conducted at Eo = 0.67 and d∕Δ c = 12.5, and the VOF advection scheme by Puckett et al. has been used.While the "dilatation and contraction method" and the "simple method" with  = 1 × 10 −2 lead to volume conservation with machine accuracy, minor volume errors are observed for the "simple method" with  = 1 × 10 −4 .As discussed in Section 3.3, this is due to the fact that small tolerance values  can lead to a certain amount of VOF over-and undershoots during the correction step, which must be cut off again.On the other hand, the volume correction is distributed more evenly over the interface cells for  = 1 × 10 −4 compared to  = 1 × 10 −2 , as more interface cells satisfy  ≤ f ≤ (1 − ).Since the illustrated errors correspond to the volume error for the entire gas phase, they are considered acceptable.Moreover, as is demonstrated in Appendix B, the flow and bubble statistics for the simulations with the "simple method" are quasi identical for  = 1 × 10 −4 and  = 1 × 10 −2 .For the rest of this section, the "simple method" has been used with  = 1 × 10 −4 .Most importantly, Figure 5 reveals that all methods conserve the gas phase volume and the slight errors on a single time step level do not accumulate over time.

Global flow statistics
This section compares the global flow statistics resulting for the two volume conservation methods across the investigated setups listed in Table 2.More precisely, the wall-normal profiles of the average gas fraction ⟨f ⟩, the normalized stream-wise velocity ⟨u⟩∕u  and the normalized turbulent kinetic energy ) are analyzed.The wall-normal distance is expressed in terms of y + = yu  ∕ l .In addition to the profiles for the simulations with enforced volume conservation, an additional profile for the respective simulation without volume conservation is provided.
Figure 6 shows the results for the simulations at Eo = 0.67 and d∕Δ c = 12.5, and for both investigated VOF advection schemes.As the figure demonstrates, the two volume conservation methods lead to identical global flow statistics.Note: The latter two columns provide the time-and bubble-averaged mean bubble deformation and oscillation strength, which are discussed in Section 4.2.
Profiles of the average gas fraction ⟨f ⟩, the normalized stream-wise velocity ⟨u⟩∕u  and the normalized turbulent kinetic energy k∕u 2  over the non-dimensional wall distance y + for the simulations at Eo = 0.67 and d∕Δ c = 12.5.The results for the simulations with the VOF advection scheme of Puckett et al. are shown in the top row (A-C), whereas the results for the simulations using the VOF advection scheme of Weymouth and Yue are illustrated in the bottom row (D-F).The results have been averaged over a time interval of 400 (H∕u  ).
[Colour figure can be viewed at wileyonlinelibrary.com]The simulations without explicit volume conservation, on the other hand, reveal a strong deviation and a significant tendency towards the profiles that would be expected for the corresponding single-phase channel flow.After averaging for a time interval of 400 (H∕u  ), the gas fraction profiles reflect the loss of gas volume (see Section 3.2), the stream-wise velocity in the channel center is observed to rise, and the turbulent kinetic energy increases in the near-wall region, whereas it decreases in the channel center.This clearly demonstrates the importance of explicit volume conservation.
Similar results are obtained for the setups at Eo = 3.75 and d∕Δ c = 12.5, which are shown in Figure 7. Again, simulations with both VOF advection schemes have been carried out.Once more, the profiles of the investigated quantities perfectly match each other for the two volume conservation methods.However, for the more deformable bubbles at Eo = 3.75, conducting the simulations without an explicit volume conservation approach leads to even stronger deviations than for the nearly spherical bubbles at Eo = 0.67.This is reasoned by the fact that the volume error increases with Eo, as demonstrated in Section 3.2.
Figure 8 illustrates the global flow statistics for the setups at Eo = 0.67 and d∕Δ f = 25.The simulations have only been performed with the VOF advection scheme of Puckett et al., and the averaging interval is reduced from 400 (H∕u  ) to 75 (H∕u  ).Due to the higher resolution, which exceeds typical LES-like resolutions, and the nearly spherical bubble shape, these simulations are the least challenging regarding the violation of volume conservation.This is also reflected by the profiles for the simulation conducted without an explicit volume conservation approach, which are nearly identical to the profiles resulting from simulations with accurate volume conservation.Only the gas fraction profile shows a slight deviation after the shorter averaging interval of 75 (H∕u  ).The profiles for the two simulations with different volume conservation methods reveal no differences.
Finally, the global flow statistics for the setups at Eo = 3.75 and d∕Δ f = 25, which are only computed with the VOF advection scheme of Puckett et al., are provided in Figure 9. Once more, identical profiles are observed for the two simulations performed with the two different volume conservation methods.Again, the averaging interval for the simulation without explicit volume conservation is reduced to 75 (H∕u  ).While the gas fraction reveals a visible deviation, the velocity and turbulent kinetic energy profiles have not changed notably during the investigated interval.To conclude, the results provided in this section demonstrate that the global flow statistics are identical for the two volume conservation methods across all investigated setups.This indicates that, besides accurately conserving the gas phase volume, the volume conservation methods do not exert a notable influence on the physics of the flow, which is the ideal scenario.Section 4.2 additionally investigates whether the individual bubble dynamics are modified by the volume conservation methods.On the other hand, the coarse-grid simulations that have been performed without an additional volume conservation approach show an increasing error in the flow statistics, which demonstrates the need for explicit volume conservation.

Individual bubble dynamics
As Section 4.1 demonstrates for a variety of test cases, the two volume conservation methods do not exhibit a notable influence on the global flow statistics.This section investigates whether this is also the case for the individual bubble dynamics, which are studied on the basis of the temporal behavior of the bubble's interface areas.The normalized interface area A b (t) of bubble b at time t is defined as where A S is the interface area of the volume-equivalent sphere.The time-averaged normalized interface area can therefore be computed as with t s and t e denoting the start and end time of statistical averaging.A b can be interpreted as the average deformation of the bubble, with A b = 1 representing a perfect sphere.Finally, the root mean square of the fluctuations of the normalized interface area around its mean value is given as A high value of A ′ b is observed for strongly oscillating bubbles, whereas a low value is indicative of a rather static bubble shape.After separately calculating A b and A ′ b for each bubble, the results are averaged over all 780 bubbles present in the domain, which yields the time-and bubble-averaged deformation and oscillation characteristics A and A ′ .
The results for the mean deformation A and strength of the oscillations A ′ computed from the bubble's interface areas across all simulated setups are provided in Table 2.As the results for A and A ′ show, the bubbles are nearly spherical and only mildly wobbling for the cases at Eo = 0.67, whereas rather ellipsoidal shapes and significant oscillations are observed for Eo = 3.75.It is important to note that, for the coarser grid resolution of d∕Δ c = 12.5 and Eo = 0.67, the values for A are below their theoretical lower limit of one, which is simply due to the fact that the numerical approximation of the interface area given in Equation ( 18) becomes less accurate on coarse grids.This clarifies that this issue is not caused by the implemented volume conservation methods.A comparison of the values for all pairs of setups that only differ regarding the volume conservation method reveals that A and A ′ are quasi identical for the respective two simulations.This indicates that there is no apparent modification of the individual bubble dynamics by the two volume conservation methods.

Computational cost
This section contains a brief investigation of the computational cost related to the two volume conservation methods.Four simulations using the VOF advection scheme of Puckett et al. have been conducted at the lower grid resolution of d∕Δ c = 12.5, as this is the more relevant resolution from an LES perspective.The investigation has been conducted for both Eötvös numbers and both volume conservation methods.For the simulations at Eo = 0.67 (Eo = 3.75), the "dilatation and contraction method" has been observed to consume around 50 (42) times more time than the "simple method".This can be explained by the significantly more sophisticated approach of the "dilatation and contraction method": in each time step, an entire additional reconstruction and advection step of the VOF field is conducted.When using the "simple method," the computational time required for volume correction accounts for 0.38% of the total computing time for Eo = 0.67, and for 0.49% of the total computing time for Eo = 3.75.For the "dilatation and contraction method," on the other hand, 15.73% of the computational time are spent on performing the volume correction for Eo = 0.67, whereas 17.38% are used for the same task at Eo = 3.75.Finally, an investigation of the performance difference for the full simulation reveals that the two simulations at Eo = 0.67 and Eo = 3.75 using the "dilatation and contraction method" have taken around 21% longer than the equivalent setups using the "simple method."

SUMMARY AND CONCLUSIONS
As this article demonstrates, the violation of volume conservation can be a significant obstacle for VOF-based simulations of turbulent bubble-laden flows on LES-like grids.Especially for high Reynolds-and Eötvos numbers, slight volume errors on a time step level can accumulate to severe deviations for long-term simulations in periodic domains.Large volume errors lead to a significant distortion of the global flow statistics and therefore make the results unreliable.To counteract this problem, two different methods for explicit volume conservation are proposed and applied to a variety of test cases.
The study covers different grid resolutions and bubble shapes and investigates whether, besides accurately conserving the tracked gas volume, the implemented approaches exert any notable influence on the global flow statistics and individual bubble dynamics.As the results reveal, the investigated statistical quantities across all studied setups are quasi identical for the two methods.This demonstrates that the physical flow characteristics are not manipulated by the implemented procedures, which justifies their application.
Both procedures can easily be added to a VOF-based solver for freely deformable bubbles and allow to perform numerical studies for flow setups that have previously been out of reach for this numerical framework.The "simple method" is more straightforward, free of kinematic assumptions and reveals an advantage in terms of computational cost.It is not parameter-free, however additional investigations have demonstrated that the results are insensitive to the choice of the single parameter.The "dilatation and contraction method" uses a more sophisticated approach that corrects the volume in a manner consistent with the transport of the VOF field.It is parameter-free and has previously been used to perform LES of turbulent bubble-laden flows. 13Both methods are applicable to a wide range of bubble shapes, including strongly deformed structures.In principle, the "simple method" can also be used to conserve the volume of arbitrarily complex VOF structures, including skirted bubbles, which are not studied in this work.This section demonstrates the necessity of removing smeared interfacial gas regions (see issue 2 in the list of origins of volume errors given in the beginning of Section 3) in order to avoid both the smearing of the interface and the related frequent formation of broken-off satellites in the context of the geometric VOF method.As shown in Table 1, the volume error caused by the removing of smeared interface regions reaches the largest value for the simulation at d∕Δ c = 12.5 and Eo = 3.75, which is computed using the VOF advection scheme of Puckett et al. 28 Therefore, this setup is chosen for demonstration purposes.To simplify the analysis, only one bubble is inserted into a converged single-phase channel flow at Re  = 590.Volume conservation is not enforced in the statistically steady state, as the investigation of a short time interval is sufficient.All other specifications are identical to those listed in Section 2.2.The effect of removing smeared interface regions and satellites is demonstrated in Figure A1.The left plot (A) compares the number of cells with VOF values of 0 < f b < 1 for three simulations, where either both smeared interfacial regions and broken-off satellites, only satellites, or neither of both are removed.If neither smeared interface regions nor satellites are removed ("nothing is removed," blue), the number of cells with VOF values of 0 < f b < 1 does not only exhibit the expected oscillatory behavior due to the bubble deformations at Eo = 3.75 (see Section 4.2), but also a rapid upward trend, in particular compared to the behavior of the other two cases.This can be explained by the fact that, in the very short time interval up to t * = t∕(H∕u  ) ≈ 2.5 simulated for this specific case, already six broken-off satellites have been formed.As these are not deleted, they clearly increase the overall gas-liquid interface area.As the multiple marker formulation assumes that a single bubble is represented by one VOF marker field, the simulation is not continued after reaching t * = t∕(H∕u  ) ≈ 2.5.The very rapid formation of broken-off satellites under the challenging conditions combining a coarse grid and a bubble at Eo = 3.75 and Re  = 590 demonstrates that the removal of broken-off satellites is inevitable, as otherwise the bubble would produce a large number of smaller satellites over time.If, on the other hand, only broken-off satellites are removed ("only satellites are removed," green), the number of cells with VOF values of 0 < f b < 1 does not grow excessively.However, as the right plot (B) shows for this simulation, a considerable number of relatively large satellites breaks off the bubble, which is consistent with the observations for the previous case ("nothing is removed," blue).As discussed previously, these gas satellites, which together exhibit volumes of up to nearly 2.2% of the initial volume of the main bubble in some time steps, must be removed from the simulation.This corresponds to a removal of large amounts of gas volume unevenly distributed in time, which occurs in 129 of the 70,497 time steps computed during the illustrated time interval.
Finally, the simulation which removes both smeared interfacial regions, as well as broken-off satellites in case of their formation ("satellites and smeared interface regions are removed," red), remains stable, while completely avoiding the formation of satellites in the investigated time interval.This demonstrates that the continuous removal of small smeared-out interfacial regions suppresses the formation of broken-off satellites, while keeping the liquid-gas interface sharp.Moreover, from a volume correction point of view, very small volume errors in each time step are preferable compared to less frequent, but much more severe volume losses.Thus, the continuous removal of smeared interfacial gas regions is an important feature for the challenging setups investigated in this work.It is important to note that, despite the fact that no broken-off satellites have been observed in the investigated time interval if the removal of smeared interface regions is activated (see "satellites and smeared interface regions are removed," red, in Figure A1), such satellites can still form very rarely, and then must be removed.However, as Table 1 demonstrates, the removal of these satellites plays a much smaller role compared to the removal of smeared interface regions.

APPENDIX B. VARIATION OF THE TOLERANCE PARAMETER FOR THE "SIMPLE METHOD"
This section investigates the influence of the tolerance parameter  defining the interface cells used for the volume correction in the "simple method" (see Section 3.3).For this purpose, two volume-conserving simulations at Eo = 0.67 and d∕Δ c = 12.5 are carried out using the VOF advection scheme of Puckett et al., and the tolerance is varied between the default value  = 1 × 10 −4 and  = 1 × 10 −2 .Figure B1 demonstrates that the global flow statistics are completely insensitive to a variation of , as long as the value is selected from a reasonable interval.The mean bubble deformation A obtained for both tolerance values is 0.9846, whereas the oscillation parameter A ′ is 1.5600 × 10 −2 for  = 1 × 10 −4 , and 1.5590 × 10 −2 for  = 1 × 10 −2 (see Section 4.2).Consequently, also the individual bubble dynamics are identical.

F I G U R E 1
Investigated downflow channel configuration.The bubbles are visualized as iso-surfaces corresponding to a volume fraction value of f b = 0.5 (gray).The left image shows the entire channel together with a wall-normal slice of the velocity magnitude.The enlarged sections on the right side illustrate the bubble shapes for the two Eötvös numbers Eo = 0.67 and Eo = 3.75.[Colour figure can be viewed at wileyonlinelibrary.com]

2
Exemplaric illustration of the removal of smeared interface regions (top row) and broken-off gas satellites (bottom row), see issues 2 and 3 in the list of volume error origins in the beginning of Section 3. [Colour figure can be viewed at wileyonlinelibrary.com]

3
Puckett et al. +1.89 −10.95 −0.38 Note: The last three columns provide a breakdown of the final (t * = 75) relative volume error for the gas phase (in percent) illustrated in Figure 3 into the three causes of volume errors listed at the beginning of Section 3. When the single errors, which are accumulated for the time interval 0 ≤ t * ≤ 75, are set off against each other, the final total errors illustrated in Figure 3 are retrieved.Global relative volume error of the gas phase  rel = (V g (t) − V g,0 )∕V g,0 given in percent, over the dimensionless time t * = t∕(H∕u  ).The left plot (a) shows the results for Eo = 0.67, whereas the right plot (b) shows the results for Eo = 3.75.The volume advection method used for the respective simulation is indicated by the abbreviations P (Puckett et al.) and WY (Weymouth and Yue) in the legend.The time interval 0 ≤ t * ≤ 75 simulated in the statistically steady state corresponds to approximately 135 flow-through times in vertical direction.[Colour figure can be viewed at wileyonlinelibrary.com]

4
Two-dimensional geometrical interpretation of the "dilatation and contraction" method.The red and green circles visualize the bubble before and after the correction step, respectively.The blue vector field represents the non-divergence-free velocity field u * .If the bubble volume is too low (V err.b < V corr.b , (A) on the left side),  is positive and the bubble volume increases.If the bubble volume is too high (V err.b > V corr.b , (B) on the right side),  is negative and the bubble volume decreases.[Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 5
Magnitude of the global relative gas phase volume error | rel | = |V g (t) − V g,0 |∕V g,0 for three simulations with different volume conservation settings at Eo = 0.67 and d∕Δ c = 12.5.All simulations have been performed with the VOF advection scheme by Puckett et al.The relative error is given in percent, and plotted over the dimensionless time t * = t∕(H∕u  ).The visualized time interval spans around 300 time steps.Missing data points correspond to time steps with | rel | = 0. [Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 7 8 9
Profiles of the average gas fraction ⟨f ⟩, the normalized stream-wise velocity ⟨u⟩∕u  and the normalized turbulent kinetic energy k∕u 2  over the non-dimensional wall distance y + for the simulations at Eo = 3.75 and d∕Δ c = 12.5.The results for the simulations with the VOF advection scheme of Puckett et al. are shown in the top row (A-C), whereas the results for the simulations using the VOF advection scheme of Weymouth and Yue are illustrated in the bottom row (D-F).The results have been averaged over a time interval of 400 (H∕u  ).[Colour figure can be viewed at wileyonlinelibrary.com]Profiles of the average gas fraction ⟨f ⟩, the normalized stream-wise velocity ⟨u⟩∕u  and the normalized turbulent kinetic energy k∕u 2  over the non-dimensional wall distance y + for the simulations at Eo = 0.67 and d∕Δ f = 25.The results have been averaged over a time interval of 75 (H∕u  ).[Colour figure can be viewed at wileyonlinelibrary.com]Profiles of the average gas fraction ⟨f ⟩, the normalized stream-wise velocity ⟨u⟩∕u  and the normalized turbulent kinetic energy k∕u 2  over the non-dimensional wall distance y + for the simulations at Eo = 3.75 and d f ∕Δ = 25.The results have been averaged over a time interval of 75 (H∕u  ).[Colour figure can be viewed at wileyonlinelibrary.com] Number of cells with VOF values of 0 < f b < 1 for the single-bubble test case computed with the VOF advection scheme of Puckett et al. at d∕Δ c = 12.5 and Eo = 3.75 over the dimensionless time t * = t∕(H∕u  ) (A).The legend indicates which structures are removed from the bubble during the simulation.The right plot (B) shows the total volume loss due to the removal of broken-off satellites in single time steps, given in percent of the initial bubble volume, for the case where only satellites are removed ("only satellites are removed," green).Values are only shown for 129 of the 70,497 time steps in the investigated time interval, as in all other time steps no satellites are removed, and thus the volume loss is zero.For the other two cases, a removal of satellites does not occur in the illustrated time interval, either because the removal of satellites is deactivated ("nothing is removed," blue), or since the removal of smeared interface regions suppresses the formation of satellites ("satellites and smeared interface regions are removed," red).[Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E B1 Profiles of the average gas fraction ⟨f ⟩, the normalized stream-wise velocity ⟨u⟩∕u  and the normalized turbulent kinetic energy k∕u 2  over the non-dimensional wall distance y + for the simulations using the "simple method" (A) and two different tolerance values  = 1 × 10 −2 and  = 1 × 10 −4 .Both simulations have been conducted at Eo = 0.67 and d∕Δ c = 12.5, and the VOF advection scheme of Puckett et al. has been used.[Colour figure can be viewed at wileyonlinelibrary.com] Case overview for the simulations performed without a volume conservation method in terms of grid resolution (d∕Δ), Case overview for the simulations performed with the two volume conservation methods in terms of grid resolution (d∕Δ), Eötvös number (Eo), VOF advection scheme and volume conservation method.