src/two-phase.h

Two-phase interfacial flows

This file helps setup simulations for flows of two fluids separated by an interface (i.e. immiscible fluids). It is typically used in combination with a Navier–Stokes solver.

The interface between the fluids is tracked with a Volume-Of-Fluid method. The volume fraction in fluid 1 is f=1 and f=0 in fluid 2. The densities and dynamic viscosities for fluid 1 and 2 are rho1, mu1, rho2, mu2, respectively.

#include "vof.h"

scalar f[], * interfaces = {f};
double rho1 = 1., mu1 = 0., rho2 = 1., mu2 = 0.;

Auxilliary fields are necessary to define the (variable) specific volume α=1/ρ and average viscosity μ (on faces) as well as the cell-centered density.

face vector alphav[], muv[];
scalar rhov[];

event defaults (i = 0) {
  α = alphav;
  ρ = rhov;
  μ = muv;
}

The density and viscosity are defined using arithmetic averages by default. The user can overload these definitions to use other types of averages (i.e. harmonic).

#ifndef ρ
# define ρ(f) (clamp(f,0,1)*(rho1 - rho2) + rho2)
#endif
#ifndef μ
# define μ(f)  (clamp(f,0,1)*(mu1 - mu2) + mu2)
#endif

We have the option of using some “smearing” of the density/viscosity jump.

#ifdef FILTERED
scalar sf[];
#else
# define sf f
#endif

event properties (i++) {

When using smearing of the density jump, we initialise sf with the vertex-average of f.

#ifndef sf
#if dimension <= 2
  foreach()
    sf[] = (4.*f[] + 
	    2.*(f[0,1] + f[0,-1] + f[1,0] + f[-1,0]) +
	    f[-1,-1] + f[1,-1] + f[1,1] + f[-1,1])/16.;
#else // dimension == 3
  foreach()
    sf[] = (8.*f[] +
	    4.*(f[-1] + f[1] + f[0,1] + f[0,-1] + f[0,0,1] + f[0,0,-1]) +
	    2.*(f[-1,1] + f[-1,0,1] + f[-1,0,-1] + f[-1,-1] + 
		f[0,1,1] + f[0,1,-1] + f[0,-1,1] + f[0,-1,-1] +
		f[1,1] + f[1,0,1] + f[1,-1] + f[1,0,-1]) +
	    f[1,-1,1] + f[-1,1,1] + f[-1,1,-1] + f[1,1,1] +
	    f[1,1,-1] + f[-1,-1,-1] + f[1,-1,-1] + f[-1,-1,1])/64.;
#endif
#endif 

#if TREE
  sf.prolongation = refine_bilinear;
  boundary ({sf});
#endif
  
  foreach_face() {
    double ff = (sf[] + sf[-1])/2.;
    alphav.x[] = fm.x[]/ρ(ff);
    muv.x[] = fm.x[]*μ(ff);
  }
  foreach()
    rhov[] = cm[]*ρ(sf[]);

#if TREE  
  sf.prolongation = fraction_refine;
  boundary ({sf});
#endif
}

Usage

Examples

Tests