Comparison of Different Numerical Interface Capturing Methods for the Simulation of Faraday Waves

Abstract

Faraday instability is a classic problem that occurs due to the relative displacement of the interface that separates two immiscible fluids placed in a closed container under oscillating acceleration parallel to gravity. The interface deformation and the induced flow patterns of this two-phase flow are very complex and numerical simulations could allow a deeper understanding of the dynamics of these systems. Some tests have been performed to establish a reference solution, but further validation is needed in order to ensure the validity of these solutions. In this work, we compare some numerical solutions for the linear and nonlinear regimes using the phase field scheme with predictions obtained using different numerical schemes such as Front Tracking, Volume of Fluid, and Element-based Finite Volume Method. The results show that, in both linear and nonlinear regimes, some important differences in the prediction of the interface dynamics between the methods are observed, and the need to provide a reference numerical solution for future benchmarks is highlighted.

1. Introduction

Instabilities in mechanical systems are present in many areas of engineering processes [1]. In general, the response of different variables of a mechanical system can be especially sensitive to changes in one or several parameters [2]. For example, in the well-known Melde’s experiments [3] in which a horizontal string is subjected to a tension at one end and made to vibrate at the other end, the formation of standing waves is possible for specific values of the excitation frequency that depends on geometrical and dynamic variables such as the length and tension of the string. In this case, the parameter controlling the wave pattern is the excitation frequency. Determining the value of the parameter range or the phase space of the parameter set, where the system is stable or not, is crucial in many applications.

The dynamics of some fundamental and industrial applications of mechanical systems involving fluids is based on the same principles. Applications include spray formation [4], liquid atomization [5], and the creation of liquid-based templates to generate patterns to print micro-structures [6]. In particular, many natural phenomena and industrial processes originate due to the effects caused by the acceleration of a container or vessel containing two fluids (both liquids or a liquid and a gas) separated by an interface. The occurrence of an instability at the interface is taken advantage of, which may result in the appearance of standing wave patterns or the detachment of small droplets of the denser fluid. This instability will be of a different nature depending on how the acceleration is applied [7]. Indeed, depending on whether the applied acceleration is constant, impulsive, or harmonic, instabilities such as Rayleigh–Taylor [8], Richtmyer–Meshkov [9], or Faraday [10] occur, respectively. Often, the solutions of the governing equations are very sensitive to the parameters that trigger the instability and describing their phase space can be a very difficult task [11]. In addition, analytic solutions for the governing equations are often not available and numerical methods have to be employed to find approximate solutions. Special care must be taken when considering different numerical approaches to ensure that the solutions are independent of the discretization scheme chosen to represent the governing equations.

In this work, we are interested in the numerical simulation of Faraday instability [10]. This instability is produced by the relative displacement of the interface separating two non-miscible fluids placed in a closed container, subjected to an oscillating acceleration parallel to gravity. If both fluids and the container are at rest, the interface that separates both fluids is flat. Faraday [10] found that, for a given value of the frequency of the oscillating acceleration, when values of the amplitude of the oscillating acceleration, $a_f$, are below a critical value, $a_c$, no appreciable motion of the interface is observed, and any perturbation of the interface is dissipated. Otherwise, if $a_f$ is greater than $a_c$, the interface becomes unstable showing complex behavior. Under certain conditions, for the threshold value $a_c$, the interface separating both fluids shows a regular standing waves pattern. Faraday [10] observed that, when $a_f$ value corresponds to the threshold value $a_c$, the frequency of the standing waves was half the frequency of the acceleration of the container while Matthiessen [12] reported that the frequency was equal to that of the container. The value of $a_c$ is a function of the dimensions of the container, the physical properties of fluids, and the oscillating frequency.

Several efforts have been made in the past to determine the values of critical parameters for the Faraday instability and to analyze the interface behavior for values of $a_f$ greater than $a_c$. The study of Faraday instability has been approached by considering two different regimes, linear and nonlinear. In the linear regime, the terms of the equations of motion involving products between velocity components, or with its derivatives, are neglected as Kumar [2] does. A variation of this approach is called weakly nonlinear and consists of expanding the nonlinear terms and neglecting the fourth- and higher-order terms [12]. Consequently, in the linear regime, the oscillation amplitudes of the interface must be relatively small compared to the wave length, so the convective term in the equations of motion can be neglected in the governing equations.

The Faraday instability has been analyzed by many authors. Kumar and Tuckerman [2] determine the critical values of the oscillating acceleration magnitude as a function of the wavenumber for a two-dimensional system consisting of a rectangular container with two ideal fluids separated by a planar interface using a Floquet analysis. They found that, for the threshold value, standing waves could be harmonic or sub-harmonic with frequencies that are equal to or half the oscillating acceleration, respectively.

Subsequently, Wright et al. [7] compare, for non-viscous fluids, the solutions from a boundary-integral method with the vortex-sheet method. For small interface deformations, the predictions of both methods are very similar. They found that plumes, droplet formation, and ejection could appear for moderate values of the Atwood number.

The Faraday instability with viscous fluids is studied by Kumar [13] who presents a linear analysis. He found that, for the same value of frequency and amplitude of the oscillating acceleration, it is possible to have both harmonic and sub-harmonics’ standing waves patterns. Later, Périnet et al. [14] use a finite-difference projection method coupled with a front-tracking method to track the interface. They perform analyses in two and three dimensions and their results for the threshold values of the oscillating acceleration amplitude are in good agreement with the Floquet analysis of Kumar and Tuckerman [2]. In order to validate the numerical results of Périnet et al. [14] with a different numerical method, Takagi and Matsumoto [15] use a numerical scheme based on the phase-field method. They perform an analysis for both linear and nonlinear regimes. Their results for the magnitude of the critical acceleration in the linear regimes are in good agreement with those of Périnet et al. [14], but, in the nonlinear regime, the comparison of the interface dynamics in their test case, taken from Wright et al. [7], is only qualitative due to differences between properties of fluids (densities), the size of the physical domain, and the consideration of viscous effects. Takagi and Matsumoto [15] show that the trends of the interface dynamics is similar in both methods and conclude that their model is able to provide accurate numerical results of the Faraday instability. Other different nonlinear cases are tested and show good agreement with experimental results of Jiang et al. [16], in particular, concerning the occurrence of the period tripling of standing waves. Takagi and Matsumoto perform another study [17] of the same mechanical system considering a two-frequency forcing, but they use a different numerical scheme, the particle-level method, to model the interface dynamics. Accurate modeling of the forces due to interfacial tension is very important to adequately describe the interface dynamics. Different methodologies are used to represent the effect of surface tension, depending on how it is represented. The region that separates fluids can be represented, in grid-based methods, as a diffuse or a sharp interface. The natural choice for representing interfacial forces, based on the definition of the forces themselves, is to express them as tangential forces to the interface. This method then requires defining the interface using adaptive meshes so that the interface coincides with a surface (line in 2D flows) of the grid, or, in fixed meshes, by introducing a set of markers that represent the interface and move with the flow. Methods using adaptive moving meshes allow interfacial conditions on the tangential and normal stresses of the surface to be imposed in a natural way, but have limitations in their application when the interface deforms considerably [18]. On the other hand, there are methods that use markers that are carried by the flow and allow for tracking the interface. These methods allow for considering large curvatures but have the disadvantage that can lead to interfaces that deform at spatial scales that cannot be represented in the grid [19]. Among the methods belonging to this family is Front Tracking [20].

On the other hand, another family of methods uses diffuse interfaces for which the term representing the interfacial forces is expressed as a volumetric distribution of forces. In this case, this originally local force (appearing only at the interface) is distributed over a small set of cells. For this purpose, the function representing an interface concentrated in a region of space is approximated by a Heaviside function that extends the interface in a direction perpendicular to itself, an ad-hoc distance (usually the width of a domain cell). The interfacial forces are then calculated from an integration over the entire fluid domain. Volume-of-fluid type methods use this approach [19]. Popinet [21] presents an excellent review of methodologies to calculate these interfacial forces.

To consider Takagi and Matsumoto [15] results as reference solutions for the 2D Faraday instability, a few points must be well established. First, comparisons in linear regime must go beyond checking in a number of cases the computed values for the critical acceleration. The dynamics of the interface must be analyzed to be sure that the actual motion is accurately represented. In addition, the comparison with the nonlinear case must also be verified quantitatively, since Takagi and Matsumoto [15] only compare qualitative trends of the interface movement at the midpoint of the container and the formation of plumes. In their simulation, Takagi and Matsumoto [15] use fluid density values different than those of Wright et al. [7] due to limitations in their implementation of the phase-field scheme to consider high Atwood numbers. In addition, the results of Wright et al. [7] were obtained in a semi-infinite domain; consequently, the influence of fluid confinement on the interface dynamics could be important.

In this work, different numerical schemes are considered for the simulation of the two-dimensional Faraday instability. We use three different codes, two in-house codes based on Front-tracking (sharp interface representation) and Volume of Fluid (diffuse interface representation) methods, respectively, to shift the interface, and the commercial software ANSYS-CFX (Element-based Finite Volume Method). We consider cases in the linear and nonlinear regimes in the same manner as Takagi and Matsumoto [15]. Comparison of the threshold accelerations and interface dynamics for the linear harmonic and sub-harmonic cases and a nonlinear case resulting in the formation of a plume have been considered.

This paper is organized as follows: first, we describe the physical domain, the mathematical framework, and the particular features of the numerical schemes for each code. In Section 3, we present the details of the verification tests for the in-house codes. Section 4 contains the numerical results for the Faraday instability while Section 5 discusses the implications of comparing the results of all codes with the phase-field code-based scheme. Finally, some consequences of the comparisons made and their implications for further developments for better modeling of the Faraday instability are presented.

2. Mathematical and Numerical Models

We consider two immiscible fluids in a two-dimensional container, separated by a sharp interface. Initially, both fluids are at rest, and, due to buoyancy effects, the lighter (heavier) fluid is naturally placed on top (bottom) of the container as shown in Figure 1.

Both fluids are Newtonian, incompressible, and immiscible, and their physical properties are constant. The fluid at the top side of the container has density and dynamic viscosity ${\rho }_t$ and ${\mu }_t$, respectively, while the fluid in the bottom side has density and dynamic viscosity ${\rho }_b$ and ${\mu }_b$, respectively. The vertical axis z is aligned with the gravity acceleration g. The container is submitted to an oscillating vertical external force resulting in a total variable acceleration a in the vertical direction as shown in Equation (1),

$$a=-g+a_{f}cos\left(\omega t\right),$$

(1)

where $a_f$ and $\omega$ are the amplitude and angular frequency of the variable component of the acceleration, respectively.

To describe the dynamics of these fluids, we consider the Navier–Stokes equations for both fluids in laminar flow regime:

$$\nabla ·\mathbf{u}=0$$

(2)

$$\rho \left[\frac{\partial \mathbf{u}}{\partial t}+\nabla ·\left(\mathbf{u}\mathbf{u}\right)\right]=-\nabla p+\nabla ·\left(\mu \left[\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^T\right]\right)+\rho \mathbf{g}+\rho {\mathbf{f}}_s$$

(3)

where $\rho$ and $\mu$ are the local fluid density and dynamic viscosity, respectively, p is the pressure, $\mathbf{u}$ is the velocity field, $\mathbf{g}$ is the gravitational acceleration, and ${\mathbf{f}}_s$ takes into account the surface tension effects between the two fluids.

We compare the results obtained with the Phase-Field, Front-Tracking, Volume-Of-Fluid, and Element-based-Volume-of-Fluid schemes after numerically solving Equations (2) and (3). In all the considered schemes, the domain is decomposed into a number of elementary cells. The two fluids are separated by an interface that is determined by either: (a) a function $\varphi$ that changes its values between two extreme values (for example −1 and 1 or between 0 and 1), the value representing the location of the interface being a specific value (e.g., $\varphi$ = 0 or $\varphi$ = 0.5), or (b) the location of a discrete set of points or markers that allow for reconstructing the interface when they are connected. In either case, a subset of cells are traversed by the interface and there the physical properties, density, and viscosity are determined so that their value is proportional to the fluid volumetric fraction of each fluid in the cell.

The boundary conditions for the system of Equations (2) and (3) are: (a) the upper and the lower boundaries are non-slip boundaries and (b) the left and right boundaries are periodic boundaries. Depending on the selected numerical scheme, a different approach is followed to account for interfacial effects. For the sake of completeness, a simplified version of each scheme is presented.

2.1. Front Tracking Scheme

Oliva [22] developed an in-house code, hereafter referred to as ${FT}_S$, based on the Front Tracking method. His implementation follows the single-field formulation proposed by Tryggvasson et al. [20,23], where a single set of conservation equations is used with the addition of a marker function that identifies the fluids. The discretization of the Navier Stokes equations is performed using the finite volume method in a staggered mesh. The linear momentum conservation equations are solved by the Semi-Implicit Method for Pressure-Linked Equations Revised, SIMPLER, (Patankar [24]). The marker points are advected from the velocity field and then used to reconstruct the marker function and estimate the surface tension force, in the same way as Perinet et al. [14]. When the simulation is performed in a two-dimensional domain, the interface separating the fluids is two-dimensional, and the structure of the front can be handled as a succession of ordered markers connected one after the other. The transfer of information between the front and the staggered grid is done by the smooth weight function proposed by Brackbill and Ruppel [25]. The reconstructing of the marker function is done by first computing its gradient and then integrating it. To maintain accuracy and efficiency, Oliva’s implementation [22] allows for restructuring of the front, involving, inserting, or removing markers in the interface when they get too far away or too close to each other while moving.

2.2. Phase-Field Based Scheme

Phase-field methods model the sharp interface as a diffuse interface. A field variable representing the local presence of each fluid changes in a controlled manner in a finite region of the space. The equation describing this field variable is the Cahn and Hilliard equation [26]. This equation is obtained by considering the time evolution of the fluid volume fraction of each phase, $FVF$, of the free energy density function of van der Waals. This equation includes the energy gradient (proportional to the local variation of $FVF$) and the second bulk energy density (proportional to a function of $FVF$ that considers the local presence rate of each phase, with minimum values at each phase and a maximum at the interface). The effect of surface tension is obtained from the variational derivation of the energy field (Jacmin [27]) and is included in the linear momentum conservation equations as an external force. The conservation equations are then solved simultaneously with the Cahn–Hilliard equation to determine the values of the flow variables and the location of the interface separating the fluids.

The numerical scheme consists of using a staggered grid to represent the variables. The advection term of the Cahn–Hilliard equation is discretized by standard second-order finite difference and aliasing is not eliminated. Concerning the discretization of the chemical potential term, it is split in three terms: the fourth order term which is fully implicit discretized while the other two terms (linear and the third order for free energy) are semi-implicit discretized. On the other hand, the Navier–Stokes equations are discretized using a finite difference projection method [28]. For more details, see Takagi and Matsumoto [15]. This scheme is hereafter called ${PF}_S$.

2.3. Volume of Fluid Scheme

In order to have a different approach than Front Tracking and Phase fields methods to numerically solve the interface motion, we developed another in-house code combining the Finite Volume method to discretize the Navier–Stokes equations and the Volume of Fluid (VOF) method to model dynamics of the two-phase flow [29].

In the implementation of Machado’s [29] VOF scheme, hereafter called ${VOF}_S$, the governing equations are solved numerically in a staggered grid, using a system of equations for the entire flow field. The different fluids are identified by a marker function that takes different values for each fluid. In the cells where the interface is located, a value proportional to the relative amount of each fluid must be calculated. As the fluids move, the interface separating them changes location. The new locations are calculated by solving a transport equation for the interface that is identified by the gradient of the marker function. Therefore, an accurate calculation of the relative values of the volume of each fluid in each cell where the interface is located and the normal vector to the interface is needed to solve its transport equation with the precision to ensure proper modeling of the flow dynamics. In our combined Finite Volume discretization-Volume of Fluid implementation, we used the Piecewise Linear Interface Calculation, PLIC, to estimate and reconstruct the interface and estimate the volume fraction of each fluid in the cells and the Efficient Least-square Volume-of-fluid Interface Reconstruction Algorithm, ELVIRA [30], to determine the local normal vector in each cell. To model the effects of surface tension, we use the Continuum Surface Force algorithm, or the CSF model, developed by Brackbill et al. [31], in particular the Marker and Cell (MAC) version in which the zero-thickness interface is replaced by a discrete-thickness interface with a finite thickness, which allows for a smooth variation of the marker function. This smoothed function is an interpretation of the interface as a transition zone, where it changes from a sharp jump in physical properties to a smoother or degraded transition, which helps to obtain more accurate results, favoring the stability of the different numerical methods. The time-advance of the algorithm follows these steps: starting from a velocity field, the interface is advected, then the position of the interface is computed by using PLIC reconstruction and the ELVIRA algorithm [30], the surface tension forces are computed by the CSF method, a new field of velocity and pressure is computed for both fluids, and so on.

2.4. ANSYS CFX Scheme

ANSYS-CFX is a commercial software package with several years on the market. Its solver is based on the finite volume method (Element-based Finite Volume Method, EbFVM). Elemental volumes are obtained from the mesh nodes. Each elemental volume is divided into small subvolumes and field variables are interpolated within each subvolume using interpolating functions. Once the equations are integrated, the contribution of all subvolumes is summed to have the elemental volume contribution to the overall solution. It uses the volume of fluid formulation to model flows with several fluids simultaneously. It has been proven in many applications and today is one of the references in commercial software for fluid flows analysis in a variety of applications.

In this work, we use the second order backward Euler scheme for temporal discretization and the High Resolution scheme for spatial discretization. This scheme is the second order of spatial accuracy in general and first order near discontinuities. To model the interface dynamics, we use the homogeneous model. This model computes common fields of all variables except the volume fraction of each phase. The effects of surface tension interactions are considered using the CSF model of Brackbill et al. [31]. More information can be found in the documentation of the software [32]. This scheme is hereafter referred to as ${EbFVM}_S$.

2.5. Validation and Verification of In-House Codes

The in house codes developed with the Front Tracking and Volume of Fluid schemes were carefully verified to ensure their ability to accurately model the deformable interfaces as well as the motion of the interfaces. This analysis is presented in Appendix A.

3. Results

The comparison between the different methodologies was performed for the linear and nonlinear regimes. The physical domain is rectangular, as shown in Figure 1, with dimensions $L_x$ × $L_z$. The boundary conditions considered were: (a) the lower ($z=0$) and upper (z = $L_z$) boundaries are non-slip walls, and (b) the left ($x=0$) and right (x = $L_x$) vertical walls are periodic boundaries (see Figure 1). The numerical domain is discretized by a uniform mesh with $N_x$ × $N_z$ cells, the size of each cell is $\Delta x$ × $\Delta z$ with $\Delta x$ = $L_x/N_x$, and $\Delta z$ = $L_z/N_z$.

Initially, to properly compare with the results presented by Takagi and Matsumoto [15], in all cases, it was considered that, at $t=0$, both fluids were at rest and separated by a slightly deformed interface. The initial configuration of the interface is:

$$z=\frac{L_z}{2}+bcos\left(kx\right)$$

(4)

where b = $L_z/\Delta z$ and $k=2\pi /L_x$.

3.1. Linear Regime

To study the linear regime, we consider the same four cases as Takagi and Matsumoto [15]. The physical properties of the fluids, gravity acceleration, vertical length, and angular frequency are given in

Table 1. Physical parameters defining linear test cases.

${\mathit{\rho }}_{\mathit{b}}$	${\mathit{\rho }}_{\mathit{t}}$	${\mathit{\mu }}_{\mathit{b}}$	${\mathit{\mu }}_{\mathit{t}}$	$\mathit{\sigma }$	g	${\mathit{L}}_{\mathit{z}}$	$\mathit{\omega }$
(kg/m3)	(kg/m3)	(Pa·s)	(Pa·s)	(N/m)	(m/s2)	(m)	(s−1)
519.933	415.667	3.908 × 10−5	3.124 × 10−5	2.181 × 10−6	9.8066	2.31 × 10−4	200$\pi$

.

A careful and detailed sensitivity analysis of mesh size and time steps was performed for all schemes for both linear and nonlinear regime; details are presented in Appendix A. The horizontal length was defined as $L_x=2\pi /k$, and the mesh discretization was $N_x$ × $N_z=128\times 128$ for all codes. Time-step was set to $\Delta t$ = 2.5 × ${10}^{-6}$ s for all methods.

Table 2. Comparison of critical amplitude values $a_c$/g in the linear regime.

k (mm−1)	Linear Theory	${\mathit{P}\mathit{F}}_{\mathit{S}}$	${\mathit{F}\mathit{T}}_{\mathit{S}}$	${\mathit{V}\mathit{O}\mathit{F}}_{\mathit{S}}$	${\mathit{E}\mathit{b}\mathit{F}\mathit{V}\mathit{M}}_{\mathit{S}}$
28	4.37	4.65	4.36	4.43	4.35
48	12.50	12.50	12.17	11.75	13.55
73	28.50	28.40	27.65	25.25	29.20
94	51.00	50.00	48.06	48.00	49.90

shows the values of the critical acceleration $a_c$ for each considered case.

The values of the critical acceleration are very close. The relative differences between these values and those predicted by the linear theory are shown in

Table 3. Relative difference between computed critical amplitude values $a_c$/g and linear theory predictions.

k (mm−1)	${\mathit{P}\mathit{F}}_{\mathit{S}}$	${\mathit{F}\mathit{T}}_{\mathit{S}}$	${\mathit{V}\mathit{O}\mathit{F}}_{\mathit{S}}$	${\mathit{E}\mathit{b}\mathit{F}\mathit{V}\mathit{M}}_{\mathit{S}}$
28	6.41	−0.2	1.37	−0.46
48	0.00	−2.64	−6.00	8.40
73	−0.35	−2.98	−11.40	2.46
94	−1.96	−5.76	−5.88	−2.16

.

In all the cases analyzed, some of the schemes show differences with the reference solution of more than 5%. The first one corresponds to the ${PF}_S$ scheme for $k=28{\mathrm{mm}}^{-1}$. Takagi and Matsumoto [15] argue that this difference would be a consequence of the relatively large value of $L_x/N_x$. However, values of $a_c$ for all the other schemes are below 1.5% using the same discretization. The cases correspond to $k=$ 48 ${\mathrm{mm}}^{-1}$, and $k=$ 73 ${\mathrm{mm}}^{-1}$ show differences of more than 8%.

Figure 2 shows the time evolution of the interface at x = $L_x/2$ for three different values of k. Although all schemes are able to describe the motion of the interface, only the ${PF}_S$ scheme shows discrepancies in the first period. In particular, the difference between the numerical predictions of the phase field scheme and the other schemes for $k=$ 94 ${\mathrm{mm}}^{-1}$ in the first period is remarkable. However, all schemes are capable of accurately describing the dynamics of the interface.

3.2. Nonlinear Regime

Once all the schemes have been validated, we proceed to analyze a case in the nonlinear regime. In the nonlinear regime, the size of the oscillation amplitude of the interface may be larger compared to the linear regime and the possibility of having plume formation and drop ejection if the external forces operate for several periods. We choose the same test case of Takagi and Matsumoto [15], a modification of a test case considered by Wright et al. [7]. The physical and geometrical parameters are described in

Table 4. Physical parameters defining the nonlinear test case.

${\mathit{\rho }}_{\mathit{b}}$	${\mathit{\rho }}_{\mathit{t}}$	${\mathit{\mu }}_{\mathit{b}}$	${\mathit{\mu }}_{\mathit{t}}$	$\mathit{\sigma }$	${\mathit{L}}_{\mathit{x}}$	${\mathit{L}}_{\mathit{z}}$	${\mathit{a}}_{\mathit{f}}$	$\mathit{\omega }$
(kg/m3)	(kg/m3)	(Pa·s)	(Pa·s)	(N/m)	(m)	(m)	(m/s2)	(s−1)
1000	818.18	0.2	0.2	0.072	1	1	0.09	0.26728

.

The interface is slightly deformed at $t=0$ according to Equation (4) with b = −0.01 m. The mean gravity acceleration is 0, and the grid discretization is $N_x\times N_z=256\times 256$ for all schemes. The time steps for each scheme, 1.18 × ${10}^{-2}$ s, 3.36 × ${10}^{-2}$ s and 1.0 × ${10}^{-2}$ s for the ${EbFVM}_S$, ${FT}_S$ and ${VOF}_S$, respectively, were chosen to ensure that the numerical results were independent of the time discretization. This analysis is shown in Appendix A.2.2. There are two important differences between the values of the physical and geometrical parameters between situations considered by Takagi and Matsumoto [15] and Wright et al. [7]. First, the original case of Wright et al. [7] considers an Atwood number, $At=({\rho }_b-{\rho }_t$)/(${\rho }_b+{\rho }_t$) = 0.65, while the case proposed by Takagi and Matsumoto [15] corresponds to $At=0.1$. Takagi and Matsumoto [15] argue that their ${PF}_S$ scheme cannot cope with Atwood numbers greater than 0.4. The other difference is the size of the physical domain. Wright et al. [7] consider an infinite physical domain in the vertical direction. For these reasons, comparisons between our results and those of Takagi and Matsumoto [15] with numerical results of Wright et al. [7] can only be qualitative. Figure 3 shows the motion of the interface point located at the midpoint of the container, i.e., at $L_x/2$. Although the trends of all numerical simulations are similar, there are some important differences between them. First, there is a very close agreement between the results of the ${FT}_S$, ${VOF}_S$ and ${EbFVM}_S$ schemes with those of Wright et al. [7] in the first two cycles (the curves corresponding to the three schemes are superimposed on each other). For $t/T_v>2$, the predictions of these three schemes show the same trends with relative maximum and minimum values at the same time as those of Wright et al. [7].

On the other hand, the results of the phase-field scheme follow the trend of those of Wright et al. [7] but not as closely as the other schemes. In the first oscillation cycle, the predicted motion of the interface midpoint of ${PF}_S$ is in an opposite (downward) direction to the behavior of the other three schemes (upward). Subsequently, the relative maximum and minimum values are delayed compared to Wright et al. [7]. In all three schemes tested in this work, we found the droplet ejection limit very close to $t/T_v=4$.

4. Discussion

Comparison between the numerical results of the three numerical schemes Volume of Fluid (${VOF}_S$), Front-Tracking (${FT}_S$), and Element-based Finite Volume Method (${EbFVM}_S$), and those of the phase field scheme [15], show several important differences. First, in the linear regime, although the acceleration threshold values, $a_c$, determined by all schemes are relatively close for all the cases analyzed. However, none of the schemes analyzed were capable of modeling the four cases considered, with relative differences of less than 5% with respect to the reference solutions. For the case $k=28{\mathrm{mm}}^{-1}$, the ${PF}_S$ scheme shows important relative differences compared to those presented by the other schemes. In the case $k=48{\mathrm{mm}}^{-1}$, the ${EbFVM}_S$ scheme presents the most important deviations while, in the other two cases analyzed, the ${VOF}_S$ scheme has the worst performance. The ${FT}_S$, ${PF}_S$ and ${EbFVM}_S$ schemes are below 3% in at least three of the four cases analyzed.

In all cases, simulations were started for fluids at rest. Only the ${PF}_S$ scheme takes time to synchronize with the reference solutions provide by the linear analysis [2]. In particular, the behavior in the subharmonic case $k=94{\mathrm{mm}}^{-1}$ as it is shown in Figure 2 is remarkable.

Concerning the analysis of the nonlinear regime, for the first two periods, the numerical results of ${FT}_S$, ${VOF}_S$ and ${EbFVM}_S$ are almost identical to those of Wright et al. [7]. Then, for $t/T_v>2.3$, the Wright et al. [7] solution separates from the other three schemes (which continue to agree). This behavior is a consequence of the confinement of fluids by the container walls which are located at dimensionless surface elevation $h/L_x=\pm$0.5. This effect is not appreciable in Wright et al. [7] because their physical domain is infinite in the vertical direction. However, the differences between the prediction of the ${PF}_S$ scheme and the other schemes are appreciable. The surface elevations predicted by the ${PF}_S$ scheme are lower in the first period than those of the vortex-sheet method of Wright et al. [7] and the three schemes developed in this work. Next, for ${PF}_S$, the interface position at the middle of the tank grows while it goes down with our three schemes (for $1<t/T_v<2.3$). After $t/T_v>2.3$${PF}_S$, minimum and maximum values are significantly different from all the other considered schemes and the solutions of Wright et al. [7].

To check the possible causes of this behavior, we performed a last simulation considering that the oscillating acceleration phase was opposite to that considered by Wright et al. [7]. The evolution of the interface midpoint position is presented in Figure 4.

Except for the simulations of Wright et al. [7], the trends of all numerical simulations are the same. Moreover, in the first period, for very small amplitudes, the agreement between all simulations is remarkable as it is shown in Figure 5. From a qualitative point of view, the effect of this phase change in acceleration is to invert the direction of the inertial force on the interface which now acts in the same direction as the interfacial force, i.e., the interface moves towards the bottom of the container, in the opposite direction compared to the original case. This difference in the behavior of the interface can be appreciated in Figure 5.

For $t/T_v>1$, the ${PF}_S$ prediction of the surface elevations diverges from the predictions of other schemes, but the trends are similar, the relative minimum and maximum values occur almost at the same time up to $t/T_v\sim 3.5$. However, from a quantitative point of view, numerical predictions of the ${PF}_S$-scheme present important differences with the other three schemes considered. Moreover, Figure 4 shows some important differences in surface elevation of the interface midpoint predicted by the ${FT}_S$ and ${VOF}_S$ schemes are appreciable when $t/T_v>3$. Figure 6 shows the shape of the interface for various time instants (the same reported in [15]), and it is evident that, after $t/T_v>3$, although the shape of the interface follows the same trend, there are some important differences between the ${FT}_S$ scheme and the ${VOF}_S$ and ${EbFVM}_S$ schemes. Figure 6 includes results from [15] even though the direction of acceleration is reversed for illustrative purposes only.

In the search for possible explanations to the difficulties encountered in numerical modeling of Faraday instability, we compare the different time scales involved with the forces present could give some ideas for future work. A dimensional analysis considering the different forces involved (gravitational, viscous, and interfacial tension) leads to the following time scales: $t_g={(\lambda /g)}^{1/2}$, $t_{\mu }=({\lambda }^2\rho /\mu$) and $t_{\sigma }={({\lambda }^3\rho /\sigma )}^{1/2}$, where $\lambda$ is a reference length (wavelength or physical domain length), $\rho$ and $\mu$ are the average densities and viscosities, respectively, $\sigma$ is the surface tension, and g is the maximum value of the modified gravity acceleration.

Table 5. Comparison of time scales in linear and nonlinear cases.

Case	${\mathit{t}}_{\mathit{g}}$ (s)	${\mathit{t}}_{\mathit{\mu }}$ (s)	${\mathit{t}}_{\mathit{\sigma }}$ (s)
$Lineark=28{\mathrm{mm}}^{-1}$	0.0040	0.6679	0.0491
$Lineark=48{\mathrm{mm}}^{-1}$	0.0023	0.1700	0.0176
$Lineark=73{\mathrm{mm}}^{-1}$	0.0015	0.0987	0.0117
$Lineark=94{\mathrm{mm}}^{-1}$	0.0010	0.0594	0.0080
$NonLinear$	3.33	4545	355

shows the values of these scales in all cases analyzed.

The time scale associated with gravity is at least one order of magnitude smaller with respect to the other two time scales in the linear regime simulations and is even three orders of magnitude in the nonlinear regime. Consequently, the difficulty in numerically simulating this problem may lie in solving all these time scales with the necessary accuracy. Consequently, an analysis should be performed using numerical schemes with higher temporal accuracy or based on other formulations such as the use of adaptive meshing at fluid interface.

5. Conclusions

In this work, a comparison is made between different numerical schemes to solve the Navier–Stokes equations for modeling the Faraday instability. Results calculated using the phase field scheme were reviewed from several test cases, in the linear and nonlinear regimes, for the Faraday instability. All cases were modeled using three different codes based on numerical schemes developed from Front Tracking and Volume of Fluid methods and a commercial ANSYS-CFX software developed from a variant of Volume of Fluid methods. The comparison was based on predictions about the evolution of the interface shape or the dynamics of a particular interface point.

Although the numerical predictions in a linear regime are relatively close, some differences between them can be appreciated. In particular, in all cases, for some numerical schemes, the value of the critical acceleration magnitude is not so close to the reference inviscid solution. Moreover, in a linear regime, the ${PF}_S$ scheme requires extra time compared to other schemes to synchronize with the corresponding standing wave in sub-harmonic cases.

In the nonlinear regime, the differences between the numerical solutions for the test case were largest between the ${PF}_S$ scheme and all other schemes considered. Moreover, the excellent agreement between the solutions obtained by the ${FT}_S$ and ${VOF}_S$ schemes and the remarkable quantitative agreement with the reference inviscid solution in the low amplitude time range allow us to conclude that the use of phase-field methods for modeling Faraday instability should be revisited. This is evident in the variant of the nonlinear case (with a reversed phase of the oscillating acceleration), which shows that the trends of all numerical schemes are the same from a qualitative point of view, but some minor but distinguishable differences are observed in the shape of the interfaces by the three different schemes tested.

Consequently, despite the important coincidences between the results obtained between some schemes, some results for both the linear and nonlinear regimes need to be clarified in order to establish a reference solution of Faraday instability. A further study is needed that considers the influence of fluid viscosities, different boundary conditions on the top and bottom walls to make a more reasonable comparison with inviscid solutions, a sensitivity analysis of mesh density and physical domain size, and the use of more accurate temporal schemes. This topic will be considered in future work.