Developing accurate, stable and thermodynamically consistent numerical methods to simulate two-phase flows is critical for many applications. We develop numerical methods to simulate two-phase flows with deforming interfaces at various density contrasts by solving thermodynamically consistent Cahn-Hilliard Navier-Stokes equations. We develop three essentially unconditionally energy-stable Crank-Nicolson-type time integration schemes. The first two time integration schemes are fully implicit based on pressure-stabilization techniques. The third approach utilizes projection method to decouple the pressure to improve the efficiency of the fully implicit method. We rigorously prove energy stability of the semi-discrete schemes for the approaches with pressure-stabilization approaches. We also prove the existence of solutions of the advective-diffusive Cahn-Hilliard operator. In the first approach we use a pressure coupled approach with both Navier-Stokes and Cahn-Hilliard equations are solved using fully implicit non-linear schemes in a block iterative manner. In the second approach we extend the block iterative method in the first approach, to a fully-coupled, provably second-order accurate scheme in time, while maintaining energy-stability. The fully coupled approach method requires fewer matrix assemblies in each Newton iteration resulting in faster solution time. The first two methods are based on a fully-implicit Crank-Nicolson scheme in time and a pressure stabilization for an equal order Galerkin formulation. We use a conforming continuous Galerkin (cG) finite element method in space equipped with a residual-based variational multiscale (RBVMS) procedure to stabilize the pressure in the first two approaches.
In the third approach we present a projection based framework extending the fully implicit method in the second approach, to a block iterative, hybrid semi-implicit-fully-implicit in time method. We use a semi-implicit time discretization for Navier-Stokes and a fully-implicit time discretization for Cahn-Hilliard equations. This third approach method requires Newton iteration only for Cahn-Hilliard resulting faster solution time at relatively larger time steps. We use a conforming continuous Galerkin (cG) finite element method in space equipped with a residual-based variational multiscale (RBVMS) procedure. However, pressure is decoupled using a projection step resulting in two linear positive semi-definite systems for velocity and pressure instead of the saddle point system of a pressure-stabilized method. We deploy all three approaches on a massively parallel numerical implementation using fast octree-based adaptive meshes. All the linear systems are solved using efficient and scalable algebraic multigrid (AMG) method. A detailed scaling analysis of all three solvers is presented.
We present comprehensive numerical experiments showing detailed comparisons with results from the literature for canonical cases, including the single bubble rise, Rayleigh-Taylor instability, and lid-driven cavity flow problems for a large range of density ratios
