Artiﬁcial viscosity model to mitigate numerical artefacts at ﬂuid interfaces with surface tension

The numerical onset of parasitic and spurious artefacts in the vicinity of ﬂuid interfaces with surface ten- sion is an important and well-recognised problem with respect to the accuracy and numerical stability of interfacial ﬂow simulations. Issues of particular interest are spurious capillary waves, which are spatially underresolved by the computational mesh yet impose very restrictive time-step requirements, as well as parasitic currents, typically the result of a numerically unbalanced curvature evaluation. We present an artiﬁcial viscosity model to mitigate numerical artefacts at surface-tension-dominated interfaces without adversely affecting the accuracy of the physical solution. The proposed methodology computes an addi- tional interfacial shear stress term, including an interface viscosity, based on the local ﬂow data and ﬂuid properties that reduces the impact of numerical artefacts and dissipates underresolved small scale inter- face movements. Furthermore, the presented methodology can be readily applied to model surface shear viscosity, for instance to simulate the dissipative effect of surface-active substances adsorbed at the inter- face. The presented analysis of numerical test cases demonstrates the eﬃcacy of the proposed methodology in diminishing the adverse impact of parasitic and spurious interfacial artefacts on the convergence and stability of the numerical solution algorithm as well as on the overall accuracy of the simulation results.


Introduction
The dynamics of an interface separating two immiscible fluids is governed by the behaviour of individual molecules, which typically have a size of less than one nanometer. Computational fluid dynamics (CFD), however, is based on continuum mechanics and treats fluids as continuous media. Molecular scales cannot be resolved using continuum mechanics, which leads to a variety of theoretical and numerical difficulties with respect to the representation of fluid interfaces using CFD. These difficulties manifest as parasitic (of unphysical origin) or spurious (of physical origin but numerically misrepresented) flow features in the vicinity of the interface. Numerical artefacts of particular practical interest in interfacial flow modelling that have attracted significant attention and research efforts are parasitic currents [1][2][3][4][5][6][7] and spurious capillary waves [1,[8][9][10][11][12][13] .
The numerical onset of unphysical flow currents in the vicinity of the interface, so-called parasitic currents , are a common oc- * Corresponding author: E-mail address: fabian.denner@gmail.com (F. Denner). currence in the modelling of surface-tension-dominated flows. Previous studies [3][4][5][6] have identified two distinct origins of parasitic currents: a) a discrete imbalance between pressure "jump" and surface tension and b) a numerically unbalanced evaluation of the interface curvature. The imbalance of the pressure gradient and surface tension can be eliminated by employing so-called balanced-force discretisation methods that assure a discrete balance between surface tension and pressure gradient, which has previously been proposed and successfully demonstrated by Renardy and Renardy [14] and Francois et al. [4] for sharp surface representations such as the ghost-fluid method [15,16] , and by Francois et al. [4] , Mencinger and Žun [5] and Denner and van Wachem [6] for methods using a continuum surface force formulation [1,17] . A numerically unbalanced evaluation of the interface curvature, which is typically associated with spatial aliasing errors resulting from the computation of the second derivative of a discrete indicator function or reconstruction [3,6] , leads to an unphysical contribution to the momentum equations (via surface tension) and, consequently, causes an unphysical acceleration of the flow in the vicinity of the interface [4,7] . Using a balanced-force discretisation of the surface tension, parasitic currents are solely dependent on the evaluation of the interface curvature and can, for instance, be eliminated ( i.e. reduced to machine precision) for certain cases by imposing the geometrically exact curvature [4][5][6][7] . Although previous studies [7,18] have shown that, with a careful discretisation of the interface curvature, parasitic currents can be reduced to machine precision in some cases once the interface has reached a numerical equilibrium, parasitic currents are still an issue of significant practical relevance for applications with evolving interfaces, as evident by the considerable body of literature on this subject published over the past five years alone ( e.g. [7,[19][20][21][22] ).
Another important numerical artefact in interfacial flows are spurious capillary waves . In a numerical framework, the shortest wavelength unambiguously resolved by the computational mesh is λ min = 2 x, with x representing the mesh spacing. Because an adequate spatial resolution of waves requires at least 6 − 10 cells per wavelength [13] , capillary waves with a wavelength of λ min ≤ λ ࣠ 3 λ min are not part of the physical solution and can, therefore, be considered to be spurious capillary waves, meaning that these waves are a response of the discretised governing equations to a perturbation of the interface but are numerically misrepresented. Therefore, spurious capillary waves are, contrary to parasitic currents, not the result of numerical errors but the result of the limitations associated with the finite resolution of the computational mesh. The origin of spurious capillary waves can be physical perturbations of the interface, for instance due to the collision of the interface with an obstacle, as well as numerical perturbations, for instance the finite accuracy of numerical algorithms, discretisation errors or parasitic currents [13] . An adequate temporal resolution of the propagation of all spatially resolved capillary waves is essential for a stable numerical solution [1,13] . The dispersion relation of capillary waves in inviscid fluids is given as [23] from which the phase velocity follows as c σ = ω σ /k, where ω σ is the angular frequency of capillary waves, k is the wavenumber, σ is the surface tension coefficient and ρ is the density of the adjacent fluids a and b. Hence, the phase speed of capillary waves increases with decreasing wavelength. This anomalously dispersive behaviour of capillary waves leads to a very rigid time-step restriction for interfacial flow simulations. For the shortest spatially resolved capillary waves, Denner and van Wachem [13] devised the capillary time-step constraint as where x denotes the mesh spacing and u is the flow velocity tangential to the interface. Denner and van Wachem [13] demonstrated that the shortest spatially resolved capillary waves are in most cases not subject to viscous attenuation, since the computational mesh is usually too coarse for typical values of the fluid viscosities to spatially resolve the vorticity generated by the shortest numerically represented waves. The capillary time-step constraint limits the simulation of interfacial flow applications and has been commonly attributed to the explicit numerical implementation of surface tension. As a result, it is widely postulated that an implicit implementation of surface tension would lift or at least mitigate the capillary time-step constraint [1,[8][9][10]18] . Recent research efforts inspired by this assumption aimed at finding an implicit or semi-implicit treatment of the surface tension to lift the time-step restrictions in interfacial flows [8,9,11,12] . Hysing [8] proposed a semi-implicit formulation of surface tension based on the CSF method [1] for a two-dimensional finite element method. The semi-implicit formulation of Hysing includes an additional implicit term which represents artificial shear stresses tangential to the interface. Raessi et al. [9] translated this methodology for finite volume methods and reported that this semi-implicit formulation of surface tension allows to exceed the capillary time-step constraint by up to factor five without destabilising the solution of the presented two-dimensional test cases. Schroeder et al. [11] proposed a two-dimensional method with a semi-implicit implementation of surface tension on a Lagrangian interface mesh ( i.e. using an explicit representation of the interface) coupled to a Eulerian mesh for the flow, presenting stable results for time-steps up to t = 3 t σ . In a similar fashion, Zheng et al. [12] recently proposed a fully-implicit coupling of a Lagrangian interface mesh to a MAC grid and showed that the method can yield stable results for time-steps up to t ≈ 10 3 t σ . However, with respect to methods that rely on an implicit interface representation, such as Volume-of-Fluid methods [24] or Level-Set methods [25,26] , and the CSF method of Brackbill et al. [1] , Denner and van Wachem [13] demonstrated that the temporal resolution requirements associated with the propagation of capillary waves is a result of the spatiotemporal sampling of capillary waves and is independent of whether surface tension is implemented explicit or implicit. Simulating the thermocapillary migration of a spherical drop, Denner and van Wachem [13] demonstrated that without external perturbations acting at the interface (such as parasitic currents), the capillary time-step constraint can be exceeded by several orders of magnitude without destabilising the solution, presenting stable results for t = 10 4 t σ using an explicit implementation of surface tension. Thus, in order to mitigate or lift the capillary time-step constraint for numerical methods that rely on an implicit interface representation, the shortest capillary waves spatially resolved by the computational mesh have to be either filtered out or damped with an appropriate method to mitigate their impact.
Issues regarding numerical artefacts are, however, not limited to interfacial flows. For instance, numerical oscillations induced by a high-order discretisation of advection terms is a longstanding issue in CFD as well as numerical heat and mass transfer, and has been the topic of extensive studies, e.g. [27][28][29][30][31] . Artificial viscosity is a well-established concept to mitigate or eliminate highfrequency oscillations in the solution and improve the stability of the numerical methodology, in particular for shock capturing and in transsonic flows, and numerical models that incorporate artificial viscosity span a wide range of explicit and implicit methods, see e.g. [32][33][34][35][36][37] . Cook and Cabot [35] suggested that the artificial grid-dependent viscosity should be chosen as to only damp wavenumbers close to the Nyquist wavenumber π / x , which in a numerical simulation is the wavenumber of the shortest spatially resolved waves. Discretisation schemes which introduce numerical diffusion to avoid oscillatory solutions, such as TVD schemes [27][28][29] , are often considered to be part of the artificial viscosity models as well. In fact, TVD schemes can be readily translated into an explicit artificial viscosity term, as for instance shown by Davis [38] .
In this study we propose an artificial viscosity model to mitigate numerical artefacts at fluid interfaces, expanding on the work of Raessi et al. [9] . The proposed artificial viscosity model can accommodate arbitrary interface viscosities and two methods to dynamically compute the interface viscosity based on the local flow conditions are presented. We present and discuss the results for a range of numerical experiments, which allow a comprehensive assessment of the efficacy of the methodology, highlight the acting physical mechanisms and provide best practice guidelines for future interfacial flow simulations using the proposed artificial viscosity model. Furthermore, our study demonstrates that the success of the proposed methodology with regards to mitigating the capillary time-step constraint is solely based on the dissipation of surface energy, irrespective of the type of implementation, contrary to previous suggestions [8,9] .
The remainder of this article is organised as follows. Section 2 introduces the governing equations and Section 3 presents the applied numerical framework. The artificial viscosity model is presented and appropriate choices of the interface viscosity as well as the acting physical mechanisms are discussed in Section 4 . The results for a variety of numerical test cases are presented in Section 5 , in order to scrutinise the proposed methodology, show its efficacy and provide a benchmark for future modelling effort s. The article is summarised and concluded in Section 6 .

Governing equations
The considered isothermal, incompressible flow of Newtonian fluids is governed by the continuity equation and the momentum equations, defined as where u is the velocity, p stands for the pressure, ρ is the density, μ is the viscosity of the fluid, t represents time, g is the gravitational acceleration and f s is the volumetric force due to surface tension. Based on hydrodynamic principles, the forces acting on the interface must balance in both phases. Neglecting gradients in surface tension coefficient, the force-balance in the direction normal to the interface is [39] (5) where subscripts a and b denote the two fluids, σ is the surface tension coefficient, κ is the interface curvature and ˆ m is the unit normal vector of the interface (pointing into fluid b).
The Volume of Fluid (VOF) method [24] is adopted to capture the interface between two immiscible fluids. In the VOF method the local volume fraction in each cell is represented by the colour function γ , defined as The interface is, therefore, situated in every mesh cell with a colour function value of 0 < γ < 1. The local density ρ and viscosity μ are defined based on the colour function γ , given as where subscripts a and b denote the two fluids. The advection of the colour function γ is defined based on the underlying flow field as ∂γ ∂t

Numerical methodology
The applied numerical framework is based on a finite volume method with a collocated variable arrangement [6,40] . The transient term of the momentum equations, Eq. (4) , is discretised using the Second-Order Backward Euler scheme and convection is discretised using a central differencing scheme. The continuity equation is discretised using a specifically constructed momentum-weighted interpolation method proposed by Denner and van Wachem [6] , that couples pressure and velocity and assures a discrete balance between pressure gradient, gravity and surface force.
The solution procedure follows a coupled, implicit implementation of the primitive variables, following the work of Denner and van Wachem [6] , where the governing equations of the flow are solved in a single linear system of equations, given as where φ is the solution vector, constituted by the solution subvectors of velocity u = (u, v , w ) T and pressure p . Each time-step consists of a finite number of non-linear iterations to account for the non-linearity of the governing equations. At the end of each nonlinear iteration the deferred terms of the equation system are updated, based on the result of the most recent non-linear iteration. This iterative procedure continues until the non-linear problem has converged to a sufficiently small tolerance. A BiCGSTAB method incorporated in the freely-available PETSc library [41,42] is utilised to solve the preconditioned linear equation system. In this study a compressive VOF methodology [43] is applied to discretise the temporal evolution of the interface, governed by Eq. (9) , using algebraic discretisation schemes to discretise Eq. (9) and transporting the colour function in a time-marching fashion. This compressive VOF methodology inherently conserves mass within the limits of the solver tolerance [43] and is able to capture evolving interfaces with similar accuracy as VOF-based interface reconstruction methods [43,44] . The interface advection and the flow solver are implemented in a segregated fashion and coupled explicitly, see Ref. [43] for details.
Based on the colour function distribution, the surface force per unit volume is discretised using the CSF model [1] as where ˆ m = ∇ γ / |∇ γ | is the interface normal vector. In order to ensure a balanced-force discretisation of the surface force, Eq. (11) is discretised on the same computational stencil as the pressure gradient [6] and no convolution is applied to smooth the surface force [45] .

Artificial viscosity model for interfacial flows
The assumption that an implicit treatment of surface tension can mitigate the capillary time-step constraint, first proposed by Brackbill et al. [1] , has motivated recent effort s to formulate a (semi-)implicit implementation of the CSF model. Following the work of Bänsch [46] , Hysing [8] proposed an additional implicit term for finite element methods that is dependent on the surface tension, which Raessi et al. [9] extended to finite volume methods, given as where ∇ 2 s is the Laplace-Beltrami operator tangential to the interface and δ = |∇γ | is the so-called interface density. The coefficient σ t δ of the tangential Laplacian of the velocity has the same unit (Pa s) as the dynamic viscosity μ and can, therefore, be seen as constituting an artificial interface viscosity. Consequently, the additional term given by Eq. (12) represents artificial shear stresses acting tangential to the interface in the interface region in addition to the shear stresses imposed by the bulk phases, increasing the dissipation of interface movement.

Artificial interfacial shear stresses
Following the analysis of Eq. (12) in the previous paragraph, the additional shear stresses acting tangential to the interface can be generalised as where μ is the interface viscosity, representing a numerical diffusion coefficient applied to smooth the topology of the fluid interface. For an isothermal, incompressible flow the momentum equations, Eq. (4) , including the additional interfacial shear stress term of the artificial viscosity model become The prevailing effect of this additional interfacial shear stress term, Eq. (13) , is an increased dissipation of surface energy through a local increase of the acting shear stresses. The underlying mechanism closely resembles the additional shear stresses imposed at the interface by the adsorption of surface-active substances, which in the literature is often described using a surface viscosity term, following the work of Scriven [47] , that is identical to Eq. (13) .
With respect to surface-active substances the surface viscosity itself is, strictly speaking, not a fluid or interface property but should merely be regarded as a phenomenological constant [39] . The additional interfacial shear stress term, Eq. (13) , also bears similarity to other types of artificial viscosity models, such as the original artificial viscosity model of von Neumann and Richtmyer [32] or the spectral-like artificial viscosity model of Cook and Cabot [35] . With respect to parasitic currents, the additional interfacial shear stresses play the traditional role of increasing the dissipation of parasitic currents, in particular since parasitic currents are spatially very localised and the associated second spatial derivatives are typically high. The effect of the additional interfacial shear stresses on spurious capillary waves is more versatile, as they reduce the frequency and phase velocity of capillary waves [48] and also increase the penetration depth of the vorticity induced by capillary waves [49] , leading to (increased) attenuation of these waves. Furthermore, the increased viscosity acting on capillary waves steepens the energy cascade of capillary wave turbulence [50] and, thus, capillary waves with small wavelength that result from wave interaction are less energetic.

Choice of interface viscosity μ
Defining an appropriate interface viscosity for a given problem is essential for the efficacy of the artificial viscosity model and the overall accuracy of the conducted simulation. If the interface viscosity is too small, parasitic and spurious artefacts cannot be counteracted with maximum effect. On the other hand, if the interface viscosity is too large, the artificial viscosity model compromises the predictive quality of the numerical framework and changes the outcome of the simulation. As a general template it is proposed to define the interface viscosity as μ = μ * ˆ δ (15) where μ * is the base value of the interface viscosity, ˆ δ = δ x is the normalised interface density and x is the mesh spacing. This formulation assures that the interface viscosity is non-zero only in the interface region as well as that the shear stress term of the artificial viscosity model is distributed in the same way as the force due to surface tension and, hence, only acts at the interface. As shown schmematically in Fig. 1 , the normalised interface density is 0 ≤ˆ δ ≤ 0 . 5 and, hence, the applied interface viscosity in any given mesh cell is μ ≤ 0 . 5 μ * . The base value of the interface viscosity can either be fixed to a predefined value or can be calculated based on the fluid properties and the flow field, such as the two specific examples discussed below.
The original formulation of the additional interfacial shear stress term proposed by Raessi et al. [9] is retained by defining the base value of the interface viscosity as This interface viscosity is proportional to the surface tension coefficient σ and the time-step t but inversely proportional to the mesh spacing x , similar to the magnitude of parasitic currents [4,6] .
Denner and van Wachem [13] demonstrated that the viscous attenuation of interface waves in numerical simulations is critically influenced by the mesh resolution, concluding that the attenuation of spurious capillary waves cannot be represented in a physically accurate manner, since the vorticity generated by these waves is typically not discretely resolved by the computational mesh. Thus, an alternative definition of the interface viscosity is proposed based on the length scale of the penetration depth of the vorticity generated by capillary waves, given as [49] where ν = μ/ρ is the kinematic viscosity. With respect to the dispersion of capillary waves, it can be assumed that the fluid properties of both phases act collectively on capillary waves [48] , so that for instance the effective kinematic viscosity is ν = ν a + ν b .
Following this rational for the purpose of defining a length scale that takes into account the natural attenuation of capillary waves, for fluids with different fluid properties and including interface viscosity, this viscous length scale can be reformulated as where ρ = (ρ a + ρ b ) / 2 is the reference density at the interface.
Hence, the base value of the interface viscosity follows as This formulation of the interface viscosity takes into account the damping provided by the bulk phases as well as the increased attenuation of interface features with decreasing mesh spacing, since flow structures associated with parasitic currents and spurious capillary waves equally reduce in size. These flow structures experience a larger natural attenuation because viscosity acts preferably at smaller scales [23] . Furthermore, since underresolved interface waves are underdamped due to an inadequate resolution of the vorticity induced by these waves [13] , a reduced mesh spacing ( i.e. higher mesh resolution) increases the viscous dissipation of those waves. If the bulk phases provide sufficient damping, e.g. in the case of overdamped capillary waves [48] , the interface viscosity becomes zero. Following the argument of Prosperetti [49] , the viscous length scale in Eq. (19) is defined as l ν = x . Thus, the resulting interface viscosity leads to a penetration depth of the vorticity generated by an interfacial wave that is of the order of magnitude of the mesh spacing. The angular frequency ω σ should be based on a representative wavelength, such as the minimum spatially resolved wavelength λ min = 2 x or the smallest wavelength that is spatially adequately resolved, e.g. λ = 3 λ min = 6 x . Note that this represents only one particular choice for the viscous length scale l ν and the angular frequency ω σ for Eq. (19) and that other choices may be more suitable for a given problem.
The particular choice of interface viscosity also dominates the spatiotemporal convergence behaviour of the artificial viscosity model. The interface viscosity defined in Eq. (16) Thus, the interface viscosity diverges with first-order under mesh refinement, with μ * = ∞ for x → 0, and converges with firstorder for decreasing time-steps, with μ * = 0 for t → 0. The interface viscosity defined in Eq. (19) with l ν = x, on the other hand, converges with order 0.5 for spatial refinement, since μ * ∝ √ x , and is independent of the time-step t .

Implementation
Hysing [8] and Raessi et al. [9] implemented Eq. (12) implicitly, assuming that an implicit implementation is required to lift the capillary time-step constraint. This widely advocated assumption has been corrected by Denner and van Wachem [13] , demonstrating that restrictions imposed by the capillary time-step constraint also hold for the implicit implementation of surface tension. Contrary to previous studies [8,9] , the additional interfacial shear stress term of the artificial viscosity model, Eq. (13) , is implemented as an explicit contribution to the momentum equations, defined as where superscripts j and j − 1 denote the current and the previous non-linear iteration, respectively. The Laplace-Beltrami operator of velocity with respect to the tangential vector of the interface is where is the velocity gradient tangential to the interface. Because the shear stress term of the artificial viscosity model is implemented explicitly, the applied time-step has to fulfil the viscous time-step constraint [9] Since the Cauchy stress tensor of the momentum equations, Eq. (4) , is implemented implicity in the applied numerical framework, only the interface viscosity μ of the artificial viscosity model has to be considered in determining the viscous time-step constraint by Eq. (23) . However, strict adherence to the viscous time-step constraint is in many cases not necessary, since the solution algorithm accounts iteratively for non-linearities of the governing equations and, hence, the additional interfacial shear stress term, Eq. (20) , is changing in every iteration. Preliminary tests have shown that, with the numerical framework presented in Section 3 , stability issues associated with the viscous time-step constraint only arise for very large interface viscosities. Dependent on the mesh spacing, the fluid properties and the applied interface viscosity, the viscous time-step constraint t μ can be more restrictive than the capillary time-step constraint t σ , since t μ ∝ x 2 and t σ ∝ x 3/2 . Comparing t μ and t σ , the critical mesh spacing follows as and the largest stable time-step becomes However, in many cases of practical interest x x c with t σ t μ and, hence, the viscous time-step constraint does not impose an additional limit to the applied time-step for such cases. The magnitude of capillary and viscous time-step constraints is further discussed for each individual test case in Section 5 . An implicit implementation of ∇ 2 s u and ∇ s u [see Eqs. (21) and (22) , respectively], as for instance described in detail by Raessi et al. [9] , can remedy this issue for cases in which t μ t σ .

Results
The results of four representative test-cases are presented, scrutinising the proposed artificial viscosity model for interfacial flows by highlighting its acting mechanisms and discussing its impact on realistic applications. In what follows, μ * ,t denotes the interface viscosity based on the work of Hysing [8] and Raessi et al. [9] , see Eq. (16) , with superscript t indicating its time-step dependency, and μ * ,λ denotes the interface viscosity determined by Eq. (19) , where superscript λ stands for its dependency on the chosen reference wavelength. For all presented simulations, the applied time-step satisfies the capillary time-step constraint, Eq. (2) , the viscous time-step constraint, Eq. (23) , as well as a Courant number Co = | u | t/ x < 1 , unless explicitly stated otherwise.

Dispersion of capillary waves
The motion of the interface is induced by surface tension only, gravity is excluded. The domain is λ in width ( x -direction) and 3 λ in height ( y -direction), identical to the domain used by Popinet [18] [51] , which is shown in the graphs as a reference. The temporal evolution of the amplitude of the capillary wave deviates increasingly from the analytical solution as the interface viscosity increases, due to the enhanced attenuation through the additional shear stresses imposed by the artificial viscosity model, given by the damping coefficient where A is the wave amplitude, t stands for time and superscripts 0 and 1 denote two extrema with respect to the temporal evolution of the wave amplitude. Fig. 3 shows the damping ratio ζ = /ω σ of both waves as a function of interface viscosity. The damping ratio increases rapidly with increasing interface viscosity and, as previously indicated, the damping effect of the artificial viscosity model on the shorter wave ( λ = 0 . 1 m ) is significantly stronger.
The temporal evolution of the capillary wave amplitude A , with and without interface viscosity is shown in Fig. 4

Interface instabilities on falling liquid films
Falling liquid films are convectively unstable to long-wave perturbations ( i.e. the wavelength is much longer than the film height) which leads to the formation of periodic or quasi-periodic wave trains [52] . Waves with a frequency higher than the neutral stability frequency f crit are attenuated by surface tension, whereas waves with a frequency lower than f crit evolve as a result of the long-wave instability mechanism [53] . The resulting solitary waves are governed by complex hydrodynamic phenomena and exhibit a dominant elevation with a long tail and steep front, typically with capillary ripples preceding the main wave hump. In inertiadominated film flows, solitary waves exhibit a separation of scales between the front of the main wave hump, where gravity, viscous drag and surface tension balance, and the tail of the wave, characterised by a balance between gravity, viscous drag and inertia [52] .
Two different film flows are considered to study the effect of the artificial viscosity model: a) the attenuation of numerical artefacts at the interface, and b) the spatiotemporal aliasing of spurious capillary waves. Both cases are simulated in a domain of size L x × L y × 0.1 h N , schematically illustrated in Fig. 5 , represented by an equidistant Cartesian mesh with 10 cells per film height h N , where h N is the height of the unperturbed film (Nusselt height) for a given flow rate according to the Nusselt flat film solution [52] . For all cases considered in this study the domain height is L y = 4 h N . The inclined substrate has an inclination angle β to the horizontal plane and is modelled as a no-slip boundary. The gasside (top) boundary is assumed a free-slip wall. A monochromatic perturbation is imposed at the domain-inlet, periodically changing  the mass flow at frequency f and amplitude A from the mean, with a semi-parabolic velocity profile prescribed for the liquid phase at the inlet and a spatially-invariant velocity prescribed at the inlet for the gas-phase where u N is mean flow velocity (Nusselt velocity) based on the Nusselt flat film solution [52] , given as The film height h (x = 0) = h N at the inlet is constant and an open boundary condition is applied at the domain-outlet [54,55] . The film is initially flat and the velocity field is fully developed.

Mitigating parasitic currents
Following the work of Denner et al. [56] and Charogiannis et al.   [56] are shown as a reference. The impact of numerical artefacts at the front of the solitary wave is clearly visible in Fig. 6 , manifesting as wiggles of the film height h (see inset of Fig. 6 a), particularly between the wave crest and the preceding trough, and as strong fluctuations of the interface velocity u (see inset of Fig. 6 b). Since the applied time-step satisfies t σ , no external perturbations of high frequency are applied and given that the numerical artefacts occur in a region of relatively high curvature and low shear stress, the observed numerical artefacts are most likely parasitic currents. No such wiggles and fluctuations are observed if the artificial viscosity model is applied with either of the two considered definitions of the interface viscosity. However, the artificial viscosity model noticeably influences the hydrodynamics of the capillary ripple preceding the solitary wave. With increasing interface viscosity ( μ * ,t = 0 . 0334 Pa s and μ * ,λ = 0 . 277 Pa s ) the amplitude of the capillary ripple preceding the solitary wave is considerably reduced and its wavelength is increased, leading to smaller variations in interface velocity. Capillary ripples form in order for the surface energy of the interface to balance the inertia of the flow and an increasing number of capillary ripples is generally observed for increasing inertia [52,58] . As a consequence of the artificial viscosity model increasing the dissipation in the vicinity of the interface (which is acting particularly at short wavelengths [48] ) and, therefore, dissipating some of the inertia of the flow, the resulting capillary ripple is smaller for increasing interface viscosity.

Spatiotemporal aliasing of spurious capillary waves
Since the interface and the flow field are clearly oriented, falling liquid films are well suited to study the onset of spatiotemporal aliasing of spurious capillary waves as a result of breaching the capillary time-step constraint given in Eq. (2) , as demonstrated by Denner and van Wachem [13] . Spurious capillary waves have typically a frequency that is significantly larger than the neutral stability frequency f crit and, hence, from a purely physical viewpoint, spurious capillary waves should be naturally attenuated. Raessi et al. [9] reported that an implicit implementation of the additional interfacial shear stress term, Eq. (12) , allowed them to conduct numerically stable simulation with time-steps up to 5 times larger than the capillary time-step constraint of Brackbill et al. [1] . As previously discussed in the introduction of this article, Hysing [8] as well as Raessi et al. [9] attributed this behaviour to the implicit implementation of this term, a hypothesis which has later been disputed by Denner and van Wachem [13] . In the current study the additional shear stress term of the artificial viscosity model, Eq. (13) , is implemented explicitly, as detailed in Section 4.3 .
A falling water film on a vertical ( β = 90 • ) substrate is simulated, as previously considered experimentally and numerically by Nosoko and co-workers [54,59] . The liquid phase has a density of ρ l = 998 kg m −3 and a viscosity of μ l = 10 −3 Pa s , whereas the gas phase is taken to have a density of ρ g = 1 . 205 kg m −3 and a viscosity of μ g = 10 −5 Pa s . The surface tension coefficient is σ = 7 . 2205 × 10 −2 N m −1 . The liquid film flow has a Reynolds number of Re = 51 . 1 with h N = 2 . 51 × 10 −4 m . The computational domain has a length of L x = 300 h N and the perturbation frequency is f = 25 s −1 with an amplitude of A = 3% . Given the neutral stability frequency for this case is f crit = 169 . 3 s −1 , the shortest spatially resolved capillary waves with λ min = 2 x and f (λ min ) = 59 . 9 × 10 3 s −1 are naturally attenuated and, hence, any observed growth of these waves is a numerical artefact. Fig. 7 a and b shows the film height as a function of downstream distance of the falling water film at two time instants for different numerical time-steps. The capillary time-step constraint following Eq. (2) is t σ = 5 . 52 × 10 −6 s , assuming u · k = 1 . 5 u N (which is the streamwise interface velocity of the unperturbed film). In Fig. 7 a the onset of spatial aliasing for the case with t = 8 . 31 × 10 −6 s and μ * = 0 is starting to be visible, even in the flat section of the film for x ࣡ 160 h N , while the other two cases do not exhibit any aliasing. The aliasing severely influences the evolution of the falling liquid film shortly after its onset, see Fig. 7 b, in particular at the crest and the preceding trough of the largest interfacial wave. No aliasing is observed if the time-step is reduced to satisfy the capillary time-step constraint ( t < t σ ), exemplified by the case with t = 3 . 92 × 10 −6 s ≈ 0 . 71 t σ in Fig. 7 . Applying the artificial viscosity model with μ * = μ * ,t = 0 . 0239 Pa s , for which t μ = 2 . 634 × 10 −5 s applies, increases the dissipation of the surface energy of spurious capillary waves and, even with t = 8 . 31 × 10 −6 s ≈ 1 . 5 t σ , results in a smooth and accurate evolution of the interface without aliasing. Note that the interface viscosity is 23.9 times larger than the viscosity of the liquid and is, therefore, increasing the penetration depth of the vorticity generated by capillary waves, see Eq. (17) , by almost factor 5. This leads to a significantly increased attenuation of spurious capillary waves, which is essential for the ability to breach the capillary time-step constraint t σ . It is possible to apply a larger time-step than in the shown example but this would generally require a larger interface viscosity and would, therefore, have a more pronounced impact on the accuracy of the simulation results. It is worth noting that despite the explicit implementation of the artificial viscosity model, the numerical algorithm is stable and spurious capillary waves are   not amplified, allowing to breach the capillary time-step constraint without consequences.

Buoyancy-driven rise of a bubble
The rise of a bubble due to the sole action of buoyancy  , on the other hand, results in a noticeably affected Froude number, which is particularly visible in the magnification of Fig. 8 a. Despite the observed differences as a result of the artificial viscosity model and corresponding interface viscosity, the Froude numbers for all three cases are well within the expected range of Froude numbers as given by empirical studies and numerical results discussed in the previous paragraph. The total kinetic energy E kin of the flow in the entire computational do- Table 1 Values of μ * and μ * /μ for different simulation time-steps t . Based on Eq. (2) , t σ 3 . 4 × 10 −6 s .   main exhibits a decreasing total kinetic energy in the domain for increasing interface viscosity as a results of the increased viscous dissipation, as seen in Fig. 8 b.

Capillary instability of a liquid jet
The capillary instability of a water jet, schematically illustrated in Fig. 9 , is simulated and analysed in order to assess the influence of the proposed artifical viscosity model on the prediction of the breakup length of the jet when the capillary time-step contraint is breached. Following the properties used by Moallemi et al. [62] , the liquid phase has a density of ρ = 9 . 97 × 10 2 kg m −3 and a viscosity of μ = 9 . 0 × 10 −4 Pa s , whereas the gas phase has a den-  Considering the equations of motion for an inviscid fluid in a cylindrical coordinate system ( r, θ , x ), Rayleigh [63] conducted a linear stability analysis on the jet radius r under the assumption that with α n r o , and showed 1 that the original equilibrium becomes unstable for n = 0 and for a reduced wavenumber κ = kr o such that | κ| < 1. In that case, the jet radius can be rewritten as and the growth rate of the unstable radius disturbance q is given by 1 In fact, Plateau [64] first proved that a liquid jet is stable for all purely non-axisymmetric deformations, but unstable to axisymmetric modes with | κ| < 1 [ i.e. modes whose wavelength λ(= 2 π r o /κ ) > 2 π r o ]. Rayleigh [63] added to Plateau's theory by describing the dynamics of the instability.
where I n is the modified Bessel's function of the n -th order. The dimensionless disturbance growth rate reaches its maximum ω max 0.3433 for κ 0.697 [63] .
Assuming that the linear approximation remains valid until the jet breaks up, the theoretical breakup time t b is obtained by calculating the time necessary for the radial disturbance = o e qt defined in Eq. (32) to reach = r o , which gives In the cases considered in this study, the inlet jet radius remains constant and equal to r o , while the inlet jet velocity periodically Based on the mechanical energy of the perturbation, Moallemi et al. [62] derived a relationship between a radius perturbation η o and a velocity perturbation ξ o , given as Consequently, by inserting Eq. (38) in Eq. (36) , the breakup length for a jet with a velocity perturbation ξ o becomes The liquid jet in the presented simulations is pulsed by a monochromatic perturbation imposed at the domain-inlet, periodically changing the mass flow with reduced wavenumber κ = 0 . 7 (to trigger the most unstable mode) and amplitude ξ o = 0 . 08 from the mean. Based on the linear stability analysis of Rayleigh [63] , the theoretical breakup length following from Eq. (39) with an ini- Simulations are conducted for μ * ∈ { 0 , μ * ,t , 2 μ * ,t } using different numerical time-steps t , see Table 1 . The interface viscosity μ * ,t (or a multiple thereof), given in Eq. (16) , is deemed to be the most adequate choice for the evaluation of the interface viscosity for the considered case, because of its dependency on the timestep t . Alternatively, the values for μ * given in Table 1 could, of course, be applied explicitly. higher than for the same case without interface viscosity. This observation stands in agreement to previous findings [65] , which suggest that an increase in viscosity increases the breakup time and length. The artificial viscosity model dissipates energy at the interface, thereby reducing the capillary-driven instability. Note that a perfect agreement with the result of the linear stability analysis is, of course, not expected, due to the limiting assumptions of the linear stability analysis as well as the limitations imposed by the finite spatiotemporal resolution of the conducted simulations and the associated difficulties in predicting what is a singular breakup event in reality [66] . However, the presented results with μ * = 0 and t < t σ are in very good agreement with the results reported by Moallemi et al. [62] for a water jet with the same fluid properties and only marginally smaller jet velocity ( W e = 14 . 8 and Oh = 9 . 44 × 10 −3 ) as considered here. Figs. 11 -13 show a slice of the computational domain at t = 3 . 06 × 10 −3 s with the interface contour and the underlying velocity field for each of the considered time-steps and interface viscosities. As previously seen in Fig. 10 , the breakup is delayed by the additional interface viscosity imposed through the artificial viscosity model. This can also be concluded from comparing Figs. 11 d and 13 d. In the case of μ * = 0 depicted in Fig. 11 d, breaching the capillary time-step limit results in visible oscillations of the interface as well as the velocity field. These oscillations are absent if the artificial viscosity model is applied, i.e. μ * > 0 , of which the results are shown in Figs. 12 d and 13 d. This is due to the artificial viscosity model damping the (spurious) capillary waves with the shortest wavelengths. In general, the artificial viscosity has a negligible influence on the overall accuracy of the solution for μ * 3 μ, as for instance observed by comparing the velocity fields in Figs. 11 a-b and 12 a-b, and based on the results plotted in Fig. 10 .

Conclusions
An artificial viscosity model to mitigate the impact of numerical artefacts at fluid interfaces has been presented and validated. The artificial viscosity model is based on the tangential Laplacian of the interfacial velocity and an appropriately chosen interface viscosity. To this effect, the proposed method belongs to the broad class of artificial viscosity models, which are already widely used in other branches of CFD and numerical heat and mass transfer, such as turbulent transsonic flows.
A comprehensive analysis of results for representative test cases have been presented to validate and scrutinise the proposed methodology. The presented results demonstrate that the major physical mechanism associated with the artificial viscosity model is an increased viscous dissipation in the vicinity of fluid interfaces and an ensuing reduction of surface energy, which acts preferably at small length scales. Applying the artificial viscosity model has been shown to allow to breach the capillary time-step constraint and reduce the adverse impact of parasitic currents. As long as the applied interface viscosity is chosen carefully, the artificial viscosity model is found to not affect the accuracy and overall predictive quality of the numerical solution in a noticeable way, while in many cases improving the accuracy of the results due to the mitigation of numerical artefacts. Although the capillary time-step constraint is more restrictive than the viscous time-step constraint for all presented cases, for certain cases, if for instance the applied mesh resolution is very high, the viscous time-step constraint can be more restrictive than the capillary time-step constraint. In such cases, an implicit implementation of the artifical viscosity model is advisable. Nevertheless, the presented results obtained with an explicit implementation of the artificial viscosity model demonstrate the overall efficacy, which is not dependent on the type of implementation, of the proposed model.
Based on the presented results, the interface viscosity derived from the previous work of Raessi et al. [9] , Eq. (16) , performs generally better than the proposed interface viscosity based on the damping of capillary waves. This finding corresponds well with the convenient properties of the interface viscosity as defined by Eq. (16) , as the interface viscosity generally increases under condition which typically lead to larger parasitic currents or faster capillary waves. The proposed interface viscosity derived from the damping of capillary waves, see Eq. (19) , was found to provide typically too much damping in its current form. Nevertheless, it provides an alternative framework for the definition of the interface viscosity.
Following the extensive analysis presented for the dispersion of capillary waves, the evolution of long-wave instabilities on falling liquid films, the buoyancy-driven rise of a bubble and the capillarydriven breakup of a liquid jet, it can be concluded that for an interface viscosity of μ 3 μ the impact on the accuracy of the results is negligible. However, dependent on the particular case, higher interfacial viscosities may be applicable without distorting the result in a significant manner. For instance, in the presented case of the long-wave instabilities on falling liquid films the applied interface viscosity μ ≈ 23 μ prevents the onset of aliasing if the capillary time-step constraint is breached, while the influence of the artificial viscosity model on the evolution of the long-wave instability is insignificant despite the high interfacial viscosity. This can be explained by the long wave length of the examined interfacial waves in this case.