Wave breaking is one of the most complicated phenomena that place at the ocean-atmosphere interface, and it is still poorly understood. The breaking of free-surface waves is characterized by a very wide range of scales.
Large-scale breaking waves are characterized by strong turbulence with a significant
amount of drops, spray and bubbles about the breaker front (whitecaps) while at the shortest scales, the stabilizing actions of gravity and surface tension dominate over the disrupting effect of the turbulence.
Breaking of ocean waves is of relevant interest because of its implications in many physical, chemical and biological processes that take place at the air-water interface, in particular, the air entrainment generated by this phenomenon is responsible for the bulk of mass transfer that occurs at the free surface of oceans, rivers, and streams, providing the transport of oxygen and carbon dioxide (which are critical to the survival of these ecosystems), as well as playing a central role in the global climate evolution.
From a mathematical point of view, for multiphase flow problems, the exact analytical solutions exist only for the simplest problems, and experimental studies of multiphase flows are not easy either. The need for numerical solutions of the governing equations has, therefore, been felt by the multiphase research community since the origin of computational fluid dynamics, in the late fifties and early sixties. Although much has been accomplished, simulations of multiphase flows have remained far behind homogeneous flows where direct numerical simulations (DNS) have become a standard tool in turbulence research.
For these reasons, in this work, we will develop an efficient direct numerical solver for the Navier Stokes equation in multi-phase flow, our goal is to produce faster and more accurate code than those available to study the wave-breaking phenomena with three-dimensional direct numerical simulations.
A common solution technique for the incompressible Navier-Stokes equations is the projection method (A.J.Chorin,1968). in which a Poisson equation must be solved numerically at each time step. Consequently, much work has gone into developing so-called "fast Poisson solvers", which use a combination of fast Fourier transforms (FFT) and Gauss elimination to solve the Poisson equation directly in Fourier space. However, fast Poisson solvers cannot solve the Poisson equation for the pressure that arises in two-fluid flows. Fast Poisson solvers require constant coefficient on the left-hand side of the Poisson equation, whereas the coefficient on the left-hand side of the pressure equation varies in space and time. To solve the variable coefficient equation, a standard practice has been to use iterative methods, e.g., multigrid methods, Krylov methods or Krylov methods preconditioned with multigrid. The downside is that iterative methods are often slower than fast Poisson solvers. In fact, they can be up to ten times slower. Additionally, an iterative methods operation count depends on problem parameters (e.g., density ratio) and convergence tolerance, whereas fast Poisson solvers have the advantage of fixed operation counts. To reduce the computational time in solving variable density incompressible flows, Guermond and Salgado (2009) adopted a penalty formulation, whereby only a constant coefficient Poisson equation must be solved at each time step. More recently, for solving the two-fluid coupled Navier-Stokes Cahn-Hilliard (phase-field) equations, Dong and Shen (2012) developed a velocity-correction method that transforms the Poisson equation from a variable to a constant coefficient equation. The underlying idea is to split the variable-coefficient pressure- gradient term into a constant term and a variable term, and then treat the constant term implicitly and the variable term explicitly leads to a constant coefficient Poisson equation for the pressure that can be solved directly using a fast Poisson solver. Therefore, by using the splitting technique in the pressure equation, there is potential for significant speedup in solving the Navier-Stokes equations for incompressible two-fluid flows. Another bottleneck in the VOF method is the interface reconstruction and update. The algorithm for geometrically reconstructing the interface affects the quality of the numerical solutions, and it remains so far an active topic. Recent signs of progress of the VOF schemes are achieved along the path of the PLIC (piecewise linear interface calculation) reconstruction. However, a geometrical reconstruction makes the computer coding complicated and filled with "if" logic operations, which usually prevents them from being immediately adaptable to the users. A VOF algorithm that does not need any geometrical reconstruction, referred to as algebraic type VOF method in our context, was developed by F.Xiao et al. (2005). The original THINC scheme effectively maintains the interface thickness but tends to ruffle the surface when the flow direction is almost parallel to the interface. A significant improvement in the interface shape faithfulness can be made by blending the THINC scheme and the first-order upwind advection scheme (Xiao et al., 2011). In conclusion, in this work, we will develop an efficient direct numerical solver for the Navier Stokes equation in multi-phase flow. Our goal is to produce a code faster and more accurate than those available using the splitting technique for the pressure equation, so that we can use a "fast Poisson solver", whereas for the advection equation of marker function we will use an algebraic technique, like to the THINC scheme (Xiao et al., 2011).