Numerical time-step restrictions as a result of capillary waves

Article history: Received 26 August 2014 Received in revised form 9 January 2015 Accepted 14 January 2015 Available online 19 January 2015


Introduction
The numerical simulation of interfaces with surface tension is only stable if the propagation of capillary waves, which are interface waves governed by surface tension, is temporally resolved by the applied numerical time-step. The phase velocity of capillary waves is given as [1] which follows from the dispersion relation (2) where ω σ is the angular frequency of capillary waves, k is the wavenumber, σ is the surface tension coefficient and ρ a and ρ b are the densities of the two fluids. Thus, capillary waves of short wavelength travel faster than capillary waves of long wavelength. For numerical simulations this anomalously dispersive behaviour of capillary waves leads to very rigid time-step restrictions. In a numerical framework, the shortest wavelength unambiguously resolved by the computational mesh and, therefore, the lower limit for the wavelength, is [2][3][4] λ min = 2 x, (3) with x representing the mesh spacing. This definition of the minimum resolved wavelength is best understood as a sampling criterion. Considering a series of standing waves, the amplitude at any given point on a wave is repeated at distances of integer-multiples of half the wavelength. Hence, in order to unambiguously represent the amplitude of the wave as well as its dispersion, the phase shift between two neighbouring sampling points, i.e. computational nodes with respect to numerical simulations, has to be smaller than π [3], which corresponds to λ/2. Because an adequate spatial resolution of waves requires at least 6-10 cells per wavelength [4], the dynamic behaviour of the shortest spatially resolved capillary waves cannot be represented in a physically accurate manner. Hence, capillary waves with a wavelength between 2 x and, approximately, 6 x are not part of the physical solution and can, therefore, be considered spurious capillary waves. The origin of spurious capillary waves are perturbations of the interface induced by external forces, for instance due to the collision of the interface with an obstacle, as well as the finite accuracy of numerical algorithms and discretisation errors, in particular parasitic currents. Parasitic currents, even if they are small in magnitude, cause a local displacement of the interface that initiates capillary waves of short wavelength. Other frequently encountered sources of spurious capillary waves are sudden changes in pressure or the impact of eddies on the interface.
Brackbill et al. [5] recognised the necessity to temporally resolve the propagation of the shortest numerically represented capillary waves and derived the capillary time-step constraint where superscript BKZ stands for the surnames of the original authors, based on the phase velocity of inviscid capillary waves c σ and the shortest spatially resolved wavelength λ min . The common consensus in the research community is that the capillary time-step constraint is required as a result of the explicit treatment of surface tension [5][6][7][8][9][10][11][12][13][14][15][16][17][18]. Brackbill et al. [5] and many consecutive researchers have suggested that this capillary time-step constraint can be mitigated if the surface force is treated implicitly. This assumption has motivated recent efforts to devise an implicit or semi-implicit implementation of the surface force, in order to lift the time-step restrictions in interfacial flow simulations [9,14]. Hysing [9] proposed a semi-implicit formulation of the surface force for finite element methods based on the CSF model of Brackbill et al. [5], introducing an additional implicit term which represents additional shear stresses tangential to the interface. Hysing [9] reported stable results for time-steps up to two orders of magnitude higher than the capillary time-step constraint. The semi-implicit formulation of Hysing [9] was later modified for finite volume methods by Raessi et al. [14], who reported that the additional shear stresses at the interface enables them to exceed the capillary time-step constraint by at least factor five without destabilising the solution. A study recently presented by Denner and van Wachem [19] suggests that the success of this additional interfacial shear stress term is based on the dissipation of the surface energy of capillary waves with short wavelength, independent of whether this additional shear stress term is implemented explicitly or implicitly. A different picture is presented by considering the resolution of capillary waves from a signal processing viewpoint. Sampling a spatiotemporally evolving wave, temporal aliasing occurs if the sampling frequency is below the Nyquist sampling rate, which stipulates that the sampling rate has to be at least twice the maximum frequency of the sampled signal in order to avoid temporal aliasing [20]. This suggests that a maximum time-step with respect to a given mesh spacing exists, similar to the capillary time-step constraint given by Eq. (4), which assures an aliasing-free result. Exceeding this time-step eventually results in aliasing of the capillary waves, leading to an ambiguous representation of the frequency and amplitude of capillary waves [3], an unphysical evolution of the interface and potential divergence of the numerical solution algorithm. Since the temporal sampling rate is independent of the applied numerical technique, the time-step constraint imposed by the Nyquist sampling rate does not depend on the numerical treatment of surface tension.
In recent years, several studies reported coupled numerical frameworks, e.g. [21][22][23][24], where the primitive variables are solved in a single system of linearised equations. A coupled implicit approach provides a strong, implicit pressure-velocity coupling, which is particularly desirable for the simulation of two-phase flows, given the discontinuous pressure jump resulting from surface tension and the potentially large density and viscosity ratios between the interacting fluids. Denner and van Wachem [25] have recently presented a coupled numerical framework for two-phase flows on arbitrary meshes, including a segregated implementation of the interface advection and an explicit treatment of surface tension. The results reported by Denner and van Wachem [25] demonstrate the robustness and accuracy of the coupled framework, even for very high surface tension and large density ratios [25,26]. However, a coupled numerical framework in which the governing equations as well as the interface advection are solved in a single equation system and in which surface tension is implemented implicitly has not been reported yet.
This article investigates the time-step requirements associated with resolving the dynamics of the equations governing capillary waves, focusing on two main objectives: (a) to determine whether explicit and implicit treatments of surface tension have different time-step requirements with respect to the dispersion of capillary waves, and (b) the formulation of an accurate time-step criterion for the propagation of capillary waves based on established numerical principles and signal processing theory. A fully-coupled numerical framework with implicit coupling of the governing equations and the interface advection and an implicit treatment of surface tension is proposed, based on the coupled framework of Denner and van Wachem [25]. This fully-coupled implicit framework is used to study the temporal resolution of capillary waves with explicit and implicit treatment of surface tension. Furthermore, a revised capillary time-step constraint is proposed, based on a detailed analysis of the numerical domain of dependence and spatiotemporal sampling requirements of capillary waves, including the superposition of capillary waves on an underlying fluid motion.
This article is structured as follows. Section 2 outlines the governing equations and the general numerical framework. In Section 3, the explicit and implicit treatments of surface tension are presented and a coupled flow-interface advection methodology is proposed. Subsequently, Section 4 investigates the dispersion of spurious capillary waves in real fluids. In Section 5, the temporal resolution requirements for the simulation of capillary waves are devised and a revised capillary time-step is proposed. In Section 6, the surface tension treatments are compared with respect to the capillary time-step constraint and the revised capillary time-step constraint is validated. The article is concluded in Section 7.

Numerical methods
The considered incompressible Newtonian fluids are governed by the continuity equation and the momentum equations, defined as where u is the velocity, t represents time, p stands for the pressure, ρ is the density, μ is the viscosity of the fluid, g is the gravitational acceleration and f s is the volumetric surface force. In Eq. (6) the transient term is discretised using the Second-Order Backward Euler scheme and convection is discretised using a central differencing scheme. If non-isothermal flows are considered, the fluids are additionally governed by the energy equation where T is the temperature, c p is the specific heat capacity and Λ is the thermal conductivity.
The Volume of Fluid (VOF) method [27] is adopted to capture the interface between two incompressible and immiscible fluids. In the VOF method, the local volume fraction in each cell is represented by the colour function γ , defined as γ (x, t) = 0 fluid a, Thus, the interface is located in every mesh cell holding a colour function value of 0 < γ < 1. The density ρ and the viscosity μ are defined based on the colour function γ as where subscripts a and b denote the two fluids. If heat transfer is considered, the specific heat capacity c p and the thermal conductivity Λ are treated in the same manner as density and viscosity. The colour function γ is advected by the hyperbolic advection equation based on the underlying flow with velocity u. In the present study a compressive VOF methodology is applied as detailed in [28], using algebraic discretisation schemes to discretise Eq. (11) and transport the colour function in a time-marching fashion. The CICSAM scheme by Ubbink and Issa [29] is used to discretise the spatial advection term in Eq. (11) and the transient term of Eq. (11) is discretised using a Crank-Nicolson scheme. The error in mass conservation of this compressive VOF method is of the order of the applied solver tolerance, as reported by Denner and van Wachem [28] for various representative test cases on structured and unstructured meshes. Assuming surface tension is a volume force acting in the interface region, the surface force per unit volume is described by the CSF model [5] as f s = σ κmδ Σ , (12) where σ is the surface tension coefficient, κ represents the interface curvature, m is the unit normal vector of the interface and δ Σ is the, so-called, interface density. Following Brackbill et al. [5], the interface density is defined as δ Σ = |∇γ |.
the volumetric surface force becomes f s = σ κ∇γ . (15) If heat transfer is considered, the temperature-dependency of the surface tension coefficient as well as the resulting thermocapillarity is included in the definition of the surface force and Eq. (15) becomes [30] f s = σ κ∇γ + ∇ s σ |∇γ |, (16) where |∇γ | is the interface density as defined in Eq. (13) and ∇ s σ is the gradient of the surface tension coefficient tangential to the interface, given as A linear relationship between temperature gradient and surface tension coefficient is assumed, with the surface tension coefficient being where σ 0 is the surface tension coefficient at temperature T 0 and σ T is the temperature coefficient of surface tension.
In order to ensure a balanced-force implementation of the surface force, Eq. (12) is discretised on the same computational stencil as the pressure gradient [25]. The pressure gradient term of the momentum equations as well as the gradient of the colour function are evaluated using the Gauss theorem, defined for the colour function gradient of mesh cell P as (19) where V P is the volume of cell P , subscript f denotes the mesh faces bounding cell P , n f is the outward-pointing unit normal vector of face f and A f stands for the area of face f . The surface force is taken as is, and no convolution is applied to smooth the surface force [31].

Solution methodology
The solution procedure follows a coupled, implicit implementation of the primitive variables, following the work of Denner and van Wachem [25], where the governing equations of the flow are solved in a single linear system of equations, given as where A is the coefficient matrix, φ is the solution vector and b is the right-hand side vector. Inside matrix A, A i j represents the coefficient submatrix for primitive variable j of the i-th equation, with superscripts x, y and z denoting the three momentum equations and c denoting the continuity equation. The solution vector φ is constituted by the solution subvectors of the three velocity components u, v and w, and pressure p. The discretised energy equation solved for temperature T is added to the linear system of equations, Eq. (20), in non-isothermal cases.
In order to provide an additional relationship between pressure and velocity to close the linear system of equations, the continuity equation is formulated using a specifically constructed advecting velocity u n and is defined as where f denotes all bounding faces of a given mesh cell P , n f is the outward-pointing unit normal vector of face f and A f is the area of face f . For isothermal two-phase flows on arbitrary meshes, Denner and van Wachem [25] devised the advecting velocity u n f based on a momentum interpolation method as where c f , d f and α f are coefficients based on the flow field and the mesh orientation, see [25] for details, and s f is the vector connecting the cell centres P and Q adjacent to face f , as depicted in Fig. 1. The gravity source term S g is discretised as [25] S g,P = 1 where a f = x f − x P is the co-variant cell-face vector. The advecting velocity resulting from Eq. (22) couples pressure and velocity and assures a discrete balance between pressure gradient, gravity and surface force. The pressure terms in Eq. (22) represent a filter that is nonzero if the pressure profile is locally cubic or higher, thus, suppressing cubic and higher-order pressure fluctuations, which are responsible for pressure-velocity decoupling on collocated meshes. In the continuity equation, Eq. (21), the velocity term u f n f as well as the pressure gradient across the cell face α fd f (p Q − p P )/|s f | of the advecting velocity u n , Eq. (22), are implemented implicitly. The advecting velocity is also applied to define the volume flux F f = u n f A f used for the convection term of the momentum equations, Eq. (6), and the spatial advection term of the interface advection equation, Eq. (11), thus, providing a consistently defined spatial advection of momentum and the VOF colour function [28].
Each time-step consists of a finite number of non-linear iterations to account for the non-linearity of the governing equations (inexact Newton method). At the end of each non-linear iteration the deferred terms of the equation system, such as the face flux F f , are updated. This iterative procedure continues until the non-linear problem has converged to a sufficiently small tolerance. The main convergence criterion considered in the presented simulations is the divergence of the velocity field, representing the conservation of continuity. The chosen tolerance, ranging between 10 −6 and 10 −8 in the simulations presented in this article, should reflect the expected length and time scales of the flow field. For two-phase flows, the required number of non-linear iterations to convergence, typically 3-30, is predominantly governed by the prevailing pressure gradient at the interface.
Two ways of implementing the interface advection and the surface force are considered in this article and explained in the following sections: (a) an explicit implementation of the surface force with a segregated advection of the interface and (b) an implicit implementation of the surface force with a coupled interface advection. Comparing these two surface force implementations will reveal the effect of an explicit and an implicit treatment of surface tension on the capillary time-step constraint.

Explicit treatment of surface tension
Applying an explicit treatment of surface tension, the interface is advected prior to computing the flow field at a given time instant in a segregated fashion, as schematically depicted on the left in Fig. 2. A linear equation system is constructed, containing the discretised VOF advection equation, Eq. (11), for every mesh cell, using the face fluxes F f resulting from the previous time-step. After the colour function field has been advected, the colour function gradient as defined in Eq. (19), the fluid properties, see Eqs. (9) and (10), and the interface curvature κ are updated. In this case, the CSF model is implemented as an explicit source term in the momentum equations. Therefore, the colour function gradient and the interface curvature updated after the interface advection are applied to calculate the surface force as defined in Eq. (15). Hence, the surface force remains constant during the iterative solution of the equation system given in Eq. (20) for any given time-step. The velocity field and the interface advection as well as the surface force and the pressure gradient are coupled explicitly.

Implicit treatment of surface tension
For the implicit treatment of surface tension, the interface advection equation becomes part of the flow equation system presented in Eq. (20). Thus, for isothermal flows, the equation system contains five linearised equations and five unknown variables for each mesh cell and follows as where A i γ represents the submatrix with the coefficients of the VOF colour function γ for equation i, with superscript γ denoting the VOF advection equation. The solution algorithm is schematically illustrated on the right of Fig. 2. In this implementation, the discretisation of the colour function gradient, Eq. (19), applied to determine the surface force f s as well as the colour function gradient across the mesh faces of the continuity equation, α fd f σ (γ Q − γ P )/|s f | as given in Eq. (22), are implemented implicitly. The interface curvature as well as the fluid properties are deferred, i.e. taken based on the result of the previous non-linear iteration. The surface force is, therefore, not a constant throughout a given time-step unlike in the explicit implementation described in the previous section, but is directly coupled and fed-back to the flow field in each iteration.

Dispersion of spurious capillary waves
The dispersion of capillary waves as given by Eq. (2) is, strictly speaking, only valid in ideal, inviscid fluids without gravity. In reality, however, the dispersion of capillary waves is influenced by gravity, molecular viscosity in the bulk phases as well as surface-active substances. Following Eqs. (1) and (2), the phase velocity c σ and angular frequency ω σ of capillary waves in ideal, inviscid fluids rises to infinity for increasing wavenumber. In real fluids, however, viscous stresses attenuate interface waves and become increasingly significant with increasing wavenumber, as the angular frequency becomes [32] where Γ represents the damping rate. The phase velocity of capillary waves reduces accordingly. The presence of surfaceactive substances has a similar effect on the dispersion of capillary waves, increasing the dissipation of surface energy and reducing the angular frequency. At high wavenumber, which is the focus of the work presented in this article, interface waves are dominated by surface tension, whereas gravity plays a negligible role. For decreasing wavenumber, however, the effect of gravity becomes increasingly important and the impact of surface tension gradually diminishes. For wavenumbers [33,34] k g < ρ g σ (26) the dispersion relation of interface waves, so-called capillary-gravity waves, is given as (27) where g is the gravitational acceleration and n Σ is the normal vector of the interface. For interface waves of very long wavelength the effect of surface tension becomes negligible.
Spurious capillary waves as well as the flow field in the vicinity of these waves are not accurately represented by the spatial discretisation and the numerical methodology. Hence, the impact of viscous stresses on spurious capillary waves is not accurately modelled in numerical simulations. The viscous length scale with respect to the frequency of capillary waves is [35,36] l ν = μ ρω , (28) which Prosperetti [35] associated with the penetration depth of vorticity generated by interface waves. The vorticity generated by interface waves decays exponentially in the viscous bulk phases with increasing distance to the interface [1,32] and this vorticity becomes insignificant at a distance d ν = √ 2πl ν from the interface [1]. As a result, the vorticity generated by the shortest spatially resolved waves is not discretely represented by the computational mesh and, thus, viscous stresses cannot act at the corresponding length scales. For instance, for a capillary wave propagating on a water-air interface with wavelength λ = 10 −3 m, d ν = 6.83 × 10 −5 m in the water phase and d ν = 2.65 × 10 −4 m in the air phase. If this wave represents the shortest spatially resolved wave, i.e. λ = 2 x, the mesh spacing is x = 5 × 10 −4 m and, consequently, the vorticity generated by this capillary wave is not discretely resolved. The shortest spatially resolved capillary waves are, therefore, not attenuated by viscous stresses, at least for a computationally significant time span after inception, and the phase velocity relevant for the temporal resolution of these waves is . (29) In order to test the approximation made in Eq. (29), the damped oscillation of a small-amplitude capillary wave in viscous fluids is simulated with various mesh resolutions. The wavelength of the capillary wave is λ = 1 m, and, as in previous studies [8,15,37], the initial amplitude of the capillary wave is h 0 = 0.01λ, with a Laplace number of La = σρλ/μ 2 = 3000.
The non-dimensional viscosity of the fluids is = 4μπ 2 /ρω σ λ 2 = 6.472 × 10 −2 . Both fluids are initially at rest, gravity is neglected and the motion of the interface is induced by surface tension only. The domain is λ in width (x-direction) and 3λ in height ( y-direction), identical to the domain used by Popinet [15]. All domain boundaries are treated as free-slip walls. The periodic motion of the capillary wave is discretised with 1000 time-steps t per period, resulting in a time-step well within the capillary time-step constraint as given in Eq. (4). The damping coefficient of the damped oscillatory motion of the standing wave is calculated as where a is the wave amplitude, t stands for time and superscripts 0 and 1 are two extrema with respect to the temporal evolution of the wave amplitude. This damping coefficient is compared with the damping coefficient Γ exact obtained from the analytical solution for a standing capillary wave in the limit of equal fluid properties, small amplitude and a domain of infinite extend, devised by Prosperetti [35,38]. Fig. 3 shows the relative damping coefficient Γ /Γ exact obtained by varying the resolution in the x-direction and the y-direction. The damping coefficient reduces rapidly for decreasing mesh resolution with respect to both coordinate axes, as observed in Fig. 3. Although the specific case of λ = 2 x = 2 y cannot be considered for a standing wave, the results strongly support the approximation made in Eq. (29) for the phase velocity of spurious capillary waves, in particular for waves with wavelength λ → λ min = 2 x.

Temporal sampling of capillary waves
The simulation of interfacial flows with surface tension is only stable if spatially resolved capillary waves are appropriately resolved by the applied time-step. The applicability of the capillary time-step constraint as presented in Eq. (4) Fig. 3. Computed damping coefficient Γ of a standing wave relative to the damping coefficient Γ exact obtained from the analytical solution of Prosperetti [35] as a function of mesh spacing. In the left graph the mesh resolution along the x-axis is varied and along the y-axis is constant y = λ/2, whereas in the right graph the mesh resolution in the x-direction is x = λ/5 and the mesh resolution along the y-axis is varied. has been demonstrated by a large number of studies since its introduction [5,8,14,15,17,39,40]. Nevertheless, the analysis presented in Sections 5.1-5.4 suggests that the capillary time-step constraint derived by Brackbill and co-workers [5] is incomplete and founded on incorrect assumptions, such as the influence of waves travelling in opposite directions (see Section 5.3). A revised capillary time-step constraint is proposed in Section 5.5 based on the presented analysis.

Numerical domain of dependence
The principle of domain of dependence, initially introduced by Courant et al. [41], stipulates that in each time-step the numerical solution must not advance further in space than the physical solution over the same time increment. Hence, the value at a given computational node at the new time level is a function of its neighbour values at the previous time level, defined as with u being the maximum propagation speed of information, t is the time-step and d is the distance between computational nodes. The condition given by Eq. (31) is widely known as CFL condition and essentially couples the numerical solution to the physical solution. With respect to capillary waves, the maximum speed of propagation is the phase velocity c σ , as defined in Eq. (1), of the shortest spatially resolved capillary wave, with a wavelength of λ min = 2 x. The shortest spatially resolved wave is equivalent to the maximum wavenumber unambiguously resolved by the computational mesh which forms the basis of the derivation of Eq. (4). The CFL condition for capillary waves can, therefore, be expressed as A similar stability condition, or in fact the same stability condition applying central differencing, can be deduced from a von Neumann stability analysis (see e.g. [42]) with respect to the numerically stable solution of capillary waves.
It is important to understand that numerical stability is not a sufficient condition for physically viable results. If the CFL condition is violated, even a numerically unconditionally stable differencing scheme may lead to unphysical results and to unsustainable numerical errors, which can eventually cause the solving algorithm to diverge.

The signal processing viewpoint
Signal processing offers an alternative viewpoint with respect to the numerical resolution of capillary waves. According to the Nyquist-Shannon sampling theorem, the sampling frequency has to be at least twice the maximum frequency of the sampled bandlimited signal, given as also known as Nyquist sampling rate, in order to avoid temporal aliasing [20]. Temporal aliasing leads to abrupt spatiotemporal variations and an ambiguous representation of frequency and phase speed of a given signal [3], with large local gradients. Even small temporal aliasing errors of a spatially propagating signal have been found to lead to substantial errors upon integration or differentiation [4,43]. At the maximum spatially resolved wavenumber, defined in Eq. (32), the sampling frequency is Given the angular frequency of capillary waves, ω σ = c σ k, the sampling frequency becomes f s ≥ c σ x . (36) Hence, the capillary time-step based on the Nyquist-Shannon sampling theorem follows as (37) which is exactly the same time-step condition as given by the CFL condition, see Eq. (33).
Crucially with respect to the numerical representation of interface waves, the time-step requirement associated with the Nyquist-Shannon sampling theorem, Eq. (37), is independent of the numerical implementation, since the Nyquist sampling rate is a result of the discrete spatiotemporal representation of a continuous signal and not the result of the numerical discretisation of the governing equations. Hence, the temporal resolution requirement imposed by the Nyquist sampling rate applies independent of the particular choice of numerical discretisation schemes and numerical implementation methodology.

Waves travelling in opposite directions
Brackbill et al. [5] derived the capillary time-step constraint given in Eq. (4) based on the phase velocity of capillary waves c σ and the assumption arguing that the value of 0.5 on the right-hand side is necessary to account for the case of two capillary waves travelling in opposite directions. Taking this assumption into account, the capillary time-step constraint introduced by Brackbill et al. [5], Eq. (4), corresponds to the time-step condition derived from the CFL condition and the Nyquist sampling rate, as given in Eqs. (33) and (37).
However, the assumption that two oppositely travelling waves divide the time-step requirement in half is incorrect. The time-step limit depends only on the local mesh spacing and the local phase velocity. The numerical differencing scheme treats each wave individually and, thus, two oppositely travelling waves do not affect the numerical stability. Hence, from a numerical viewpoint, the time-step is not directly affected by waves propagating in opposite directions. However, two oppositely travelling waves are superimposed and may decrease the time interval between successive local extrema, thus, locally increasing the effective wavenumber sampled by the numerical algorithm, which becomeŝ k = k I + k II , (39) where subscripts I and II denote the two oppositely travelling waves. The special case of two oppositely travelling waves with λ = 2 x entering the same mesh cell at the same time-step results in a standing wave and has, therefore, no impact on the phase velocity. In fact, a standing wave forms whenever two waves of the same wavelength and frequency interfere [34].
Waves with a wavelength λ < 2 x are not unambiguously represented by the spatial discretisation and can be ignored. The highest local wavenumber is, hence, generated by two oppositely travelling waves, where one wave has the minimum wavelength λ I = 2 x, while the other wave has a wavelength λ II → 2 x and satisfies the condition λ II > 2 x. Thus, the wavenumbers of these waves are k I = k max and k II → k max , respectively, and, based on Eq. (39), the interference of these waves results in an effective wavenumber of The increase in effective wavenumber k caused by oppositely travelling capillary waves leads to a numerical misrepresentation of the dispersion of these waves. Given the dispersion relation of capillary waves, see Eq. (2), this increase in effectively sampled wavenumber corresponds to a maximum increase in phase velocity by factor √ 2 and the effective phase velocity becomeŝ

Doppler shift
Considering the motion of a two-phase flow, waves on the interface may be considered in a reference frame moving with the flow velocity tangential to the interface. In this case a stationary computational mesh has to be regarded as a stationary observer, resulting in a Doppler shift of the capillary wave frequency. The phase velocity of capillary waves relative to the stationary mesh is, hence, where k is the vector of the wavenumber, u Σ is the velocity of the flow at the interface and u k Σ is the velocity of the flow at the interface parallel to the wavenumber vector. Thus, the underlying motion of the flow tangential to the interface increases the maximum phase velocity of capillary waves relative to the stationary computational mesh and, consequently, reduces the maximum time-step that assures a robust numerical solution.

Revised capillary time-step constraint
Based on the analysis presented in Sections 5.1-5.4, a revised capillary time-step constraint is proposed. Assuming a static case (u Σ · k = 0) and given the time-step requirement for capillary waves derived from the CFL condition, Eq. (33), and the Nyquist-Shannon sampling theorem, Eq. (37), as well as the increase in effective phase velocity caused by waves travelling in opposite directions, Eq. (41), the revised capillary time-step constraint is defined as where superscript stat denotes the static case. Comparing this revised capillary time-step constraint with the capillary timestep constraint introduced by Brackbill et al. [5], given in Eq. (4), it is obvious that the revised capillary time-step constraint defined in Eq. (43) is less restrictive. The only difference between these two time-step constraints are the assumptions about oppositely travelling waves. The static case, however, only applies to a very limited number of scenarios, since most flows feature a finite velocity along the interface. Considering fluid motion tangential to the interface by including the Doppler shift as defined in Eq. (42), the revised capillary time-step constraint becomes where superscript dyn denotes the dynamic case. Strictly speaking, the definition of the time-step constraint as given in Eq. (44) is only correct on equidistant meshes. On arbitrary meshes it is more descriptive to define the capillary timestep constraint based on the vector s Σ connecting neighbouring computational nodes in the interface region. The general formulation of the revised capillary time-step constraint can, thus, be expressed as Any time-step that fulfills this requirement should result in a robust numerical solution without aliasing.

Numerical validation
Two representative test cases are considered in this study to scrutinise the numerical stability of the explicit and implicit treatments of surface tension and to validate the proposed capillary time-step constraint. First, the numerical stability of the proposed solution methods is tested without perturbation of capillary waves, focusing only on the implementation of surface tension. Second, the accuracy and applicability of the revised capillary time-step constraint is validated, with explicit and implicit treatment of surface tension.

Numerical stability
The capillary time-step constraint is commonly attributed to numerical stability considerations associated with the explicit treatment of surface tension [5][6][7][8][9][10][11][12][13][14][15][16][17][18], a claim which has never been substantiated by theoretical analysis or numerical experiments. If the capillary time-step constraint is a requirement for numerical stability due to the explicit treatment of surface tension, the capillary time-step constraint should only depend on the magnitude of the force due to surface tension and, therefore, must be valid even if no capillary waves are propagating on the interface. A well-known example of such a stability requirement is the explicit treatment of diffusion, which is only stable for time-steps t ≤ ρ x 2 /2D φ , where D φ is the relevant diffusion coefficient. As elaborately explained by Patankar [44], among others, exceeding this time-step limit changes the coefficients of the primitive variable and, as a result, violates the physical cause-effect relationship. However, with the explicit as well as the implicit treatment of surface tension, the volumetric source term representing the surface tension, Eq. (15), is a function of the surface tension coefficient, the interface curvature and a passive scalar, i.e. the interface indicator function. The source term representing the surface tension is not a function of the primitive variables and the type of implementation does not change the matrix coefficients of the primitive variables. The implementation (explicit or implicit) of the surface force as defined by Brackbill et al. [5], given in Eq. (15), should, therefore, not impact the numerical stability. Consequently, if interface perturbations are absent and assuming all other time-step requirements that are not a result of capillary waves are fulfilled, the stability of the solution algorithms presented in Section 3 should be independent of the applied time-step.
In order to scrutinise the general numerical stability of the explicit and implicit treatments of surface tension, the thermocapillary migration of a fluid particle with Marangoni number Ma = ρc p U r r 0 /Λ = 0.01 and Reynolds number The surface tension coefficient σ is varied to accommodate different capillary numbers Ca = U r μ/σ and the temperature dependence of the surface tension is σ T = −0.0012 N m −1 K −1 . Because Re 1, the interface can be regarded as nondeformable and the exact values for interface curvature and interface normal vector are prescribed at the interface. All domain boundaries are taken to be free-slip walls and the mesh resolution is 10 cells per fluid particle diameter. Because the applied numerical framework maintains an exact discrete balance between pressure gradient and surface force, and because the geometrically exact interface normal vector and curvature are prescribed, parasitic currents are absent and the interface is not perturbed by any other source [25]. Young et al. [48] reported an analytical solution for the migration velocity of a fluid particle at zero Marangoni number in the creeping flow regime in a domain of infinite extend, given as where superscript ∞ denotes the domain of infinite extend and subscripts c and d denote the continuous and the dispersed phase, respectively. Klein and Feuerbach [49] found Eq. (46) to hold for a Marangoni number of Ma ≤ 2 with respect to the continuous phase and Ma ≤ 10 −2 of the dispersed phase. In order to account for the finite size of the domain, the theoretically devised migration velocity is corrected using the semi-empirical correlation of Harmathy [50] and is given as with L being the domain size perpendicular to the direction of migration. Fig. 4 shows the temporal evolution of the migration velocity of the fluid particle for three different time-steps and a capillary number Ca = 0.01, with the explicit as well as the implicit surface tension treatment. The predicted migration velocities are in excellent agreement with each other as well as with the expected value as given by Eq. (47), even if the static capillary time-step constraint, Eq. (43), is breached by three orders of magnitude and irrespective of the surface tension treatment. The acceleration phase of the fluid particle is clearly not accurately captured with t = 10 3 t stat , yet the terminal migration velocity is still predicted accurately. Crucially, the numerical stability is unaffected by the particular implementation of surface tension, even for t = 10 3 t stat ≈ 1.41 × 10 3 t BKZ . Furthermore, the numerical stability and predictive quality of the migration velocity are independent of the capillary number, see Fig. 5, and the density ratio, see Fig. 6. This demonstrates that the capillary time-step constraint is completely independent of the surface tension treatment, especially considering that the treatment of surface tension is irrelevant for the numerical stability even for the case with Ca = 10 −2 , ρ c /ρ d = 10 24 and t = 10 3 t stat as well as the case with Ca = 10 −4 and t = 10 4 t stat . It is worth mentioning that the Courant number for all simulated time-steps is Co = |u| t/ x < 0.81.

Temporal resolution of capillary waves
The discrete representation of capillary waves with a wavelength close to the minimum spatially resolved wavelength, i.e. λ → 2 x + , is inherently inadequate and the interface curvature of such waves is numerically ill-defined, which means the dynamics of the interface are no longer accurately represented and which can lead to parasitic currents that impact  the result. Consequently, it is very difficult to numerically generate such spurious capillary waves in a way that allows to compare different time-steps and numerical methodologies and, at the same time, observe the evolution of the spurious capillary waves. Film flows are ideally suited to study the propagation of capillary waves and, in particular, the time-step associated with capillary waves, since the interface as well as the flow field are clearly oriented. Moreover, waves with a frequency higher than the neutral stability frequency f c are naturally attenuated by the acting surface tension, whereas waves with a frequency lower than f c evolve into long-wave instabilities [51]. Since spurious capillary waves typically have a frequency higher than the neutral stability frequency, assuming an adequate mesh resolution, any observed short wave amplification is a numerical artifact and not part of the physical solution.
A vertically falling water film in contact with air is simulated to determine the impact of the treatment of surface tension on the capillary time-step constraint and to verify the revised capillary time-step constraint proposed in Section 5.5. The where u N is the mean velocity, based on the Nusselt flat film solution [53], The velocity of the gas phase at the inlet is set to u(x = 0, h N < y ≤ 3h N ) = 1.5u N . The film is initially flat, with h(x = 0) = h N prescribed at the inlet, and the velocity field is fully developed. The frequency of the shortest spatially resolved capillary waves, f max ≥ 59.9 × 10 3 s −1 for Re = 51.5 and f max ≥ 49.5 × 10 3 s −1 for Re = 75, are two orders of  respectively. Hence, any observed amplification of spurious capillary waves is purely numerical. In the considered cases, interface waves are perturbed by an explicitly imposed interface curvature of κ in = (A x) −1 at the inlet, where A is a constant depending on the simulation setup, providing a time-invariant perturbation and natural excitation of instabilities. In particular the sudden perturbation of the flat film at the beginning of the simulation triggers interface waves with a wide range of wavelengths. The magnitude of the perturbation, controlled by constant A, has to be chosen carefully, because if the perturbation is too small, capillary waves vanish immediately. On the other hand, if the magnitude of the perturbation is too large, the largest interface waves triggered by the perturbation break and small-scale capillary waves are absorbed. For the presented simulations, constant A is ranging from 9.5 to 22.15. Fig. 7 shows the result of this film flow at Reynolds numbers Re = 51.5 and Re = 75, simulated with t = t BKZ σ and treating surface tension as an explicit source term. In both cases, the interface near the inlet adapts to the explicitly imposed curvature at the inlet and triggers a long-wave instability, clearly visible in Fig. 7 for both Reynolds numbers. This is expected since a vertically falling film is unstable to long-wave perturbations for all finite Reynolds numbers [53]. No unphysically amplified spurious capillary waves are visible in either case. In Figs. 8 and 9, results for the simulated Reynolds numbers are shown, this time simulated with t = 2 t BKZ σ and treating surface tension as an explicit and an implicit source term. For both Reynolds numbers the interface waves generated by the perturbation at the inlet are artificially amplified, regardless of the surface tension treatment. The unbounded amplification of these waves leads, eventually, to the divergence of the solving algorithm. The type of surface tension treatment, i.e. explicit or implicit, has evidently no effect on the capillary time-step constraint with respect to numerical stability. The amplification of spurious capillary waves observed in Figs. 8 and 9 is smaller for the implicit treatment of surface tension for both Reynolds numbers, which can be attributed to the different feedback structure of the explicit and the implicit treatments of surface tension, as illustrated in Fig. 2.    magnitude of the perturbation is increased, see Fig. 12, the surface tension is treated implicit, as shown in Fig. 13, or if the mesh resolution is increased, see Fig. 14. Although small-scale waves are visible in Figs. 11 and 12 applying t = 0.99 t dyn σ , these waves are not amplified and are advected without aliasing, hence, these waves do not cause numerical problems. Interestingly, the spurious capillary waves evolve more quickly using the implicit surface tension treatment than with the explicit surface tension treatment for Re = 75 and t = 0.99 t stat σ , as observed in Figs. 11 and 13, contrary to the behaviour   observed with t = 2 t BKZ σ shown in Fig. 9. This suggests that an a priori statement about the evolution of spurious capillary waves cannot be made with respect to the surface tension treatment.
It is worth pointing out that the Doppler shift relative to the effective phase velocity of the fastest capillary waves in the presented cases is only u k Σ /ĉ σ ≈ 0.073 for Re = 51.5 and, dependent on the mesh spacing, u k Σ /ĉ σ ≈ 0.081-0.099 for Re = 75. This difference is small compared to its profound effect on the numerical simulation of capillary waves and the accurate prediction of this sharp transition provides a strong argument for the revised capillary time-step constraint.

Conclusions
The temporal resolution of spatially resolved capillary waves typically imposes a significant restriction on the applied time-step associated with solving the governing equations and is a major contributor to the required computational resources for two-phase flow simulations. The hypothesised reason for this time-step restriction is the explicit treatment of surface tension and it is widely assumed that an implicit treatment of surface tension would remove or at least mitigate the constraint on the temporal resolution [5,9,14]. To elucidate the actual mechanism for the restriction of the time-step associated with solving the equations governing a two-phase flow dominated by surface tension, this article has revisited the hypotheses for the time-step restriction and has compared the explicit treatment and the implicit treatment of the force resulting from surface tension.
A fully-coupled numerical framework for two-phase flows with an implicit implementation of surface tension has been introduced in this article. This fully-coupled framework has then been used to compare the influence of the surface tension treatment on the time-step restrictions resulting from capillary waves. The conducted study demonstrates that restrictions on the numerical time-step resulting from capillary waves are valid and unchanged regardless of the numerical treatment of surface tension. Since surface tension is not a function of pressure or velocity, the change in implementation does not affect the matrix coefficients of the primitive variables and, thus, numerical stability is independent of the treatment of surface tension. Further analysis shows that the capillary time-step constraint is a requirement imposed by the spatiotemporal sampling of capillary waves, which is independent of the applied numerical methodology.
The presented analysis of temporal resolution requirements associated with capillary waves reveals shortcomings of the original capillary time-step constraint, as introduced by Brackbill et al. [5]. A revised capillary time-step constraint has been proposed in Eq. (45), based on sound numerical and signal processing principles. The revised capillary time-step constraint also includes the Doppler shift resulting from flow relative to a stationary computational mesh, which has been neglected in previous studies. The presented results successfully validate the accuracy of the revised capillary time-step constraint and highlight the impact of the Doppler shift.