On sharp surface force model: Effect of sharpening coefficient

Amongst the multitude of approaches available in literature to reduce spurious velocities in Volume of Fluid approach, the Sharp Surface Force (SSF) model is increasingly being used due to its relative ease to implement. The SSF approach relies on a user-defined parameter, the sharpening coefficient, which determines the extent of the smeared nature of interface used to determine the surface tension force. In this paper, we use the SSF model implemented in OpenFOAM® to investigate the effect of this sharpening coefficient on spurious velocities and accuracy of dynamic, i.e., capillary rise, and static bubble simulations. Results show that increasing the sharpening coefficient generally reduces the spurious velocities in both static and dynamic cases. Although static millimeter sized bubbles were simulated with the whole range of sharpening coefficients, sub-millimeter sized bubbles show nonphysical behavior for values larger than 0.3. The accuracy of the capillary rise simulations has been observed to change non-linearly with the sharpening coefficient. This work illustrates the importance of using an optimized value of the sharpening coefficient with respect to spurious velocities and accuracy of the simulation.


Introduction
Modelling surface tension dominant multiphase flows is relevant in a multitude of industrial processes like lab-on-chip, atomization, and boiling. One of the main approach to capture the interface dynamics is the Volume of Fluid method which uses the advection of scalar volume fraction based on algebraic (interface compression) or geometric (piecewise linear interface calculation or PLIC) reconstruction algorithms in order to preserve the sharpness of interfacial region (Cifani et al., 2016). The VOF based solver available in OpenFOAM ® , interFoam, which generates an interface which is smeared over a few computational cells uses the interface compression method due to its relative simplicity (Deshpande et al., 2012).
In the VOF approach used in interFoam, the volume fraction field is used to determine curvature and corresponding surface tension force based on models like the widely used Continuum Surface Force (CSF) approach (Brackbill et al., 1992). Due to the smeared nature of the interface, the curvature and the pressure jump across the interface obtained from the simulations do not match the theoretical value which generates spurious velocities (Deshpande et al., 2012). These spurious velocities introduce nonphysical flows near the interface which may cause the bubble to numerically drift as well as alter the heat/mass transfer coefficients in supersaturation and temperature driven phase change processes (Samkhaniani and Ansari, 2016;Saufi et al., 2019;Vachaparambil and Einarsrud, 2020). The works by Popinet (2018) and Deshpande et al. (2012) have reviewed the various approaches reported to mitigate these effects, namely: improved curvature estimation, force balance between surface tension and pressure gradient (for static cases), time step constraint when surface tension is calculated explicitly and temporally implicit approach to estimate surface tension.
Amongst the approaches proposed, the Sharp Surface Force model, developed by Raeini et al. (2012), estimates surface tension based on a smoothed interface curvature and a sharpened interface region defined using a user defined sharpening coefficient (C sh ). This model, which is relatively simple to implement compared to height function based approach (Pavuluri et al., 2018), has been shown to reduce spurious velocities in comparison to commonly used CSF model (Pavuluri et al., 2018;Vachaparambil and Einarsrud, 2019a). The SSF model, VVol. 3, No. 3, 2021, 226-232 Experimental andComputational Multiphase Flow https:// doi.org/10.1007/s42757-020-0063-5 using C sh [0, 0.5] Î , has been used to simulate dynamic cases like rising bubbles (Vachaparambil and Einarsrud, 2019a), microfluidic T-junction (Soh et al., 2016), microchannels (Pavuluri et al., 2018), capillary rise (Raeini et al., 2012;Vachaparambil and Einarsrud, 2019a), interfacial mass transfer (Maes and Soulaine, 2018), and bubble growth (Vachaparambil and Einarsrud, 2020) whereas when modelling static cases, like stationary millimeter sized bubble, C sh is set equal to 0.98 (Vachaparambil and Einarsrud, 2019a). The choice of the value C sh used in the simulations is often heuristic and to the best knowledge of the authors there has not been a systematic attempt to quantify the effect of this user-defined parameter.
In this paper, we investigate the effect of the sharpening coefficient used in the SSF model, as developed on OpenFOAM ® 6 by Vachaparambil and Einarsrud (2019a), to model two dimensional dynamic cases like capillary rise and static cases like millimeter and sub-millimeter bubbles. All the simulations discussed in this work use the sharpening coefficient typically used in simulating practical flow scenarios, i.e.,  (2020)).

Governing equations
The volume fraction ( 1 α ) used in VOF method is a scalar field that is zero in gas phase, unity in the liquid phase, and 1 0 1 α < < at the interface. The interface dynamics is captured based on the advection of 1 α as where U is the velocity in both phases and the third term is the interface compression term that acts in the interfacial region to prevent excessive smearing using U r which is defined as where α C , φ , f S , and n represent adjustable compression factor (set equal to unity as recommended by Greenshields (2019)), volumetric flux across the cell face, cell face surface area, and unit normal to interface respectively (Deshpande et al., 2012). The fluid properties like density ( ρ ) and viscosity ( μ ) are calculated as where 2 1 1 α α = -. The mass conservation of the incom-pressible phases is described using continuity equation as The momentum equation is written based on a modified pressure ( rgh p ), defined as rgh where F ST is the surface tension modelled based on SSF described in Vachaparambil and Einarsrud (2019a) and Raeini et al. (2012). Initially, a smoothed volume fraction is obtained using a three consecutive smoothing steps ( i = 1, 2, 3) as where 1 s α = 1 α , C is equal to 0.5, and 1 c f α  < > represents the interpolation of 1 α from cell center to face. The unit normal to the smoothed interface is calculated and corrected for the effects of contact angle, see Vachaparambil and Einarsrud (2019a). Subsequently an initial estimate of curvature is calculated as where , δ defined as is used to prevent denominator from becoming zero. The curvature is smoothed using a two step procedure (i=1, 2) as The final curvature is calculated as The surface tension is estimated based on ST where C sh is the sharpening coefficient which when equal to zero produces sh α that is equivalent to 1 α . Due to the coupled nature of Eq. (4) and Eq. (5), these equations are solved by Pressure-Implicit with Splitting of Operators (PISO) algorithm (Deshpande et al., 2012). PISO algorithm involves estimation of a predicted velocity that is used to calculate pressure, using pressure correction equation, which is used to update the velocity in an iterative manner (Deshpande et al., 2012). In order to reduce spurious velocities, the force balance between pressure gradient, surface tension, and gravitational force due to discretization is ensured by calculating the gradients at cell faces as described in Deshpande et al. (2012). However the iterative procedure used to solve rgh , p i.e., the PISO algorithm, converges based on a user defined tolerance (Deshpande et al., 2012). This tolerance, required to calculate rgh p , introduces a force imbalance between surface tension, gravitational force, and pressure gradient which can be reduced by setting a very low convergence criterion, like 10 20 used in Table 1, as recommended by Deshpande et al. (2012).

Computational domain and solver settings
The governing equations are discretized using first and second order methods in time and space respectively, see Vachaparambil and Einarsrud (2019a), and solved based on methods described in Table 1. Other numerical settings like the sub-cycling of volume fraction equation and momentum predictor, which are relevant in solving the governing equations, are set based on OpenFOAM ® default settings/recommendations for simulating multiphase flows which has also been used in Vachaparambil and Einarsrud (2019a). The simulations are run with no under-relaxation factor and maximum time step is calculated as ( ) m/s 2 whereas stationary bubble simulations neglect gravity (Yamamoto et al., 2017;Vachaparambil and Einarsrud, 2019a). The computational domain used for the capillary rise simulations is 20 mm × 1 mm and meshed with a hexahedral grid of 400 × 20. This mesh is chosen based on the work by Yamamoto et al. (2017) that investigated the α is zero gradient at the outlet, Dirichlet condition equal to one at inlet, and zero gradient with a constant contact angle of 45° at the side boundaries. The modified pressure ( rgh p ) uses Dirichlet condition, equal to zero, at inlet and outlet but the side walls are assigned the zero gradient condition. The boundary conditions for U at side boundaries are set as no slip whereas the inlet and outlet are assigned a pressureinlet outlet velocity condition (Greenshields, 2019). The simulations are initialized with liquid column at a height of 8 mm (from the inlet) in the computational domain. These simulations are run until 1.5 s which is enough to reach steady capillary rise height with maximum time step, calculated based on Eq. (11), equal to 3.5 μs.
In order to model a stationary bubble of radius R, which is initialized at the center of a square computational domain of dimensions 4 4 R Ŕ , gravity is neglected. The four boundaries are assigned zero gradient condition for U and 1 α but the rgh p employs a Dirichlet condition equal to the operating pressure (equal to 101,325 Pa). The bubbles modelled in this work are a millimeter sized bubble of radius equal to 2.5 mm and a sub-millimeter bubble of radius equal to 0.25 mm. These simulations are run until 0.05 s and the corresponding time step constraints based on the mesh resolution are discussed in Section 4.2 and Section 4.3.

Results and discussions
In order compare the results from the dynamic and static simulations, spurious velocities, denoted by U sc , are calculated as max(| | U ). The time averaging of an arbitrary parameter Φ and spurious velocities are represented as with an over bar as Φ and sc U respectively.

Capillary rise
For 2D capillary rise, the equilibrium height (h T ) at which when gravitational force balance the vertical component of surface tension force for a liquid column rising between two parallel plates can be theoretical calculated as where , θ t, and ρ D are the contact angle, distance between parallel plates (equal to 1 mm), and difference between densities of the phases respectively (Bullard and Garboczi, 2009). h T , based on Eq. (12), is equal to 9.91 mm and the capillary rise height from the simulations is calculated as where the numerator represents the area of the liquid in the computational domain. The capillary rise from the simulations is compared to Eq. (12) in Table 2. The temporal evolution of the capillary rise heights and spurious velocities (U sc ) are plotted in Fig.  1 and Fig. 2 respectively. Although S h obtained from the simulations stabilize after the initial transients, the capillary rise height obtained with sh 0 3 C = . oscillates slightly, by approximately 9.30 ± 0.009 mm, as shown in Fig. 1. This oscillation in the interface position, using sh 0 3 C = . , also cause the periodic variation of sc U which is shown in Fig.  2. As the oscillations in capillary rise height are lower than ± 0.1% of the capillary rise height, we assume that the simulations have converged reasonably. h S obtained when using C sh = 0.5 matches the capillary rise height reported by Vachaparambil and Einarsrud (2019a) using SSF model. It is also worth pointing out that sc U obtained using sh 0 0 C = . is an order larger than spurious velocities obtained with other values of sharpening coefficients, see Table 2 and Fig. 2.

Millimeter sized stationary bubble
The Laplace pressure in a 2D bubble can be theoretically calculated using the Young-Laplace equation as  which for the bubble radius of 2.5 mm in the simulation is equal to 28 Pa. For simulations, the Laplace pressure in the bubble is calculated as where p 0 is the operating pressure (equal to 101,325 Pa). The mesh resolution and the time step constraints (calculated based on Eq. (11)) used in the simulations are summarized in Table 3. The stationary millimeter bubble has been modelled with the three meshes as well as for a range of sharpening coefficients between 0 and 0.5, see Table 4 and Fig. 3. Spurious velocities are observed on both sides of the interface for all the cases modelled, as illustrated in Fig. 3. The use of larger sharpening coefficients seems to reduce the error in calculating Laplace pressure as well as spurious velocities in the simulations, see Table 4. Decreasing the mesh size does not always exacerbate spurious velocities which is contrast to the increasing sc U observed with CSF model in the work by Deshpande et al. (2012) and Vachaparambil and Einarsrud (2019a). The variation between sc U reported in Table 4  (c) sh 0 2 C = .
(f) sh 0 5 C = .  Einarsrud (2019a) is due to the difference in the sh C and solver setting, in Table 1, used for the simulations.

Sub-millimeter sized stationary bubble
around 50-60 is typically used in thermal and supersaturation driven phase change processes (Samkhaniani and Ansari, 2016). Consequently, a sub-millimeter bubble, of radius equal to 0.25 mm, is initialized in a 1 mm 2 domain that is meshed by 120×120 cells and the corresponding maximum time step, calculated  (11), is set at 0.6 μs. The Laplace pressure in the sub-millimeter bubble, equal to 280 Pa, is compared to the corresponding value obtained from simulations in Table 5. The sub-millimeter bubble could be modelled with sh 0 3 C . ≤ , see Fig. 4 and Fig. 5, but for sharpening coefficient of 0.4 and 0.5 the bubble numerically drifts. For sh 0 2 C . ≤ , the interface undergoes slight periodic deformation which is reflected in the oscillations in U sc observed in Fig. 6. This deformation is not substantial enough to observe a noticeable deviation from the circular bubble shape, see Fig. 4. At t = 0.05 s, the simulations using C sh = 0.3 seems to have very low spurious velocities on both sides of the interface when compared to other sharpening coefficients, see Fig. 5.

Conclusions
The effect of sharpening coefficient used in Sharp Surface Force model, developed in the work by Vachaparambil and Einarsrud (2019a), is investigated for capillary rise and stationary bubbles of radii equal to 0.25 and 2.5 mm. The solver ensures force balance between pressure gradient, surface tension, and gravitational force due to discretization    (c) sh 0 2 C = .
(d) sh 0 3 C = .  and iterative procedure used to solve for p rgh as recommended by Deshpande et al. (2012). In order to prevent the growth of spurious velocities, time step constraint based on fluid viscosity and density as well as mesh size, in Eq. (11), proposed by Deshpande et al. (2012) is used. The simulations for a range of value of sharpening coefficients, sh 0 05 C . ≤ ≤ , shows that  The use of a larger value of C sh generally reduces the spurious velocities in capillary rise and stationary bubble simulations.  The mesh refinement does not always exacerbate spurious velocities, see Table 4, unlike while using CSF model (Deshpande et al., 2012;Vachaparambil and Einarsrud, 2019a).  The millimeter sized bubble can be modelled with sh 0 05 C . ≤ ≤ and the three meshes. Using the finest M3 mesh and C sh equal to 0.5 provides the lowest spurious velocities as well as the most accurate prediction of Laplace pressure.  The sub-millimeter bubble can be modelled with sh 0 3 C . ≤ but the lowest spurious velocities and error in Laplace pressure are observed when C sh = 0.3.  The capillary rise simulations show a non-linear variation of h S with C sh albeit the decrease in spurious velocities. The reduced spurious velocities and error in the capillary rise height is obtained when using C sh = 0.2. Although this paper investigated the effect of C sh for a few flow scenarios, the results show the importance of choosing an optimized value of the sharpening coefficient for future applications of SSF model to simulate two-phase flow phenomena.

Open Access This article is licensed under a Creative Commons
Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.