Numerical method

Interface capturing

To begin with, I consider an indicator function \(H\), which works as a marker taking \(0\) or \(1\) when the infinitesimal control volume is occupied with the primary or the secondary liquids, respectively. By this definition, as an equation governing the evolution of \(H\), I have

\[\pder{H}{t} + \frac{u_j}{h_{\gcs^j}} \pder{H}{\gcs^j} = 0.\]

Since I assume the liquids are incompressible, this advection equation is equal to

\[\pder{H}{t} + \frac{1}{J} \pder{}{\gcs^j} \left( \frac{J}{h_{\gcs^j}} u_j H \right) = 0,\]

whose volume-integrated form leads to

\[J \pder{\phi}{t} + \int_{\partial V} \frac{J}{h_{\gcs^j}} u_j H n_j dS = 0,\]

where \(\phi\) is the rate of the primary phase volume inside the control volume, or known as volume-of-fluid. Note that I assume the coordinate systems do not change in time. To integrate this equation, I employ the THINC/QQ method (Xie and Xiao, J. Comput. Phys. (349), 2017), whose numerical treatment is extensively discussed in the other project.

Consistency

Note that the following idea is widely known in the phase-field community (e.g. Mirjalili and Mani, J. Comput. Phys. (426), 2021), which is slightly customised and applied here.

As discussed above, how the volume fraction is evolved is totally determined by the volume-of-fluid method. Since \(H\) and \(\rho\) are related by

\[\rho = 1 + \left( \hat{\rho} - 1 \right) H,\]

the evolution of the density (i.e. mass per unit volume) is also governed by the VOF method. Since the momentum is tightly linked to the flux of the mass, it is natural to apply this information (mass flux) to update the momentum field. To this end, I focus on the evolution of \(\phi\) again:

\[J \pder{\phi}{t} + \pder{}{\gcs^j} \left( \frac{J}{h_{\gcs^j}} H u_j \right) = 0.\]

As the flux of the indicator function \(H u_j\) is already known, the density flux \(\rho u_j\) is obtained by

\[\rho u_j = u_j + \left( \hat{\rho} - 1 \right) H u_j.\]

src/interface/mass_flux.c

for(int j = 0; j <= jsize + 1; j++){
  for(int i = 1; i <= isize + 1; i++){
    FLUXX(i, j) = UX(i, j) + (denr - 1.) * FLUXX(i, j);
  }
}

src/interface/mass_flux.c

for(int k = 0; k <= ksize + 1; k++){
  for(int j = 0; j <= jsize + 1; j++){
    for(int i = 1; i <= isize + 1; i++){
      FLUXX(i, j, k) = UX(i, j, k) + (denr - 1.) * FLUXX(i, j, k);
    }
  }
}

src/interface/mass_flux.c

for(int j = 0; j <= jsize + 1; j++){
  for(int i = 0; i <= isize + 1; i++){
    FLUXY(i, j) = UY(i, j) + (denr - 1.) * FLUXY(i, j);
  }
}

src/interface/mass_flux.c

for(int k = 0; k <= ksize + 1; k++){
  for(int j = 0; j <= jsize + 1; j++){
    for(int i = 0; i <= isize + 1; i++){
      FLUXY(i, j, k) = UY(i, j, k) + (denr - 1.) * FLUXY(i, j, k);
    }
  }
}

src/interface/mass_flux.c

for(int k = 0; k <= ksize + 1; k++){
  for(int j = 0; j <= jsize + 1; j++){
    for(int i = 0; i <= isize + 1; i++){
      FLUXZ(i, j, k) = UZ(i, j, k) + (denr - 1.) * FLUXZ(i, j, k);
    }
  }
}

Surface tension force

For simplicity, the continuum surface force model (Brackbill et al., J. Comput. Phys. (100), 1992) is adopted to model the interfacial tension:

\[f_i \approx \frac{2 \rho}{1 + \hat{\rho}} \sigma \kappa \frac{1}{h_{\gcs^i}} \pder{\phi}{\gcs^i}.\]

src/interface/force.c

for(int j = 1; j <= jsize; j++){
  for(int i = 2; i <= isize; i++){
    const double den_x0 = (
        + 0.5 * DEN(i-1, j  )
        + 0.5 * DEN(i  , j  )
    );
    const double kappa = (
        + 0.5 * CURV(i-1, j  )
        + 0.5 * CURV(i  , j  )
    );
    const double delta = 1. / HXXF(i  ) * (
        - VOF(i-1, j  )
        + VOF(i  , j  )
    );
    IFRCX(i, j) = den_x0 * refdeninv * tension * kappa * delta;
  }
}

src/interface/force.c

for(int k = 1; k <= ksize; k++){
  for(int j = 1; j <= jsize; j++){
    for(int i = 2; i <= isize; i++){
      const double den_x0 = (
          + 0.5 * DEN(i-1, j  , k  )
          + 0.5 * DEN(i  , j  , k  )
      );
      const double kappa = (
          + 0.5 * CURV(i-1, j  , k  )
          + 0.5 * CURV(i  , j  , k  )
      );
      const double delta = 1. / HXXF(i  ) * (
          - VOF(i-1, j  , k  )
          + VOF(i  , j  , k  )
      );
      IFRCX(i, j, k) = den_x0 * refdeninv * tension * kappa * delta;
    }
  }
}

src/interface/force.c

for(int j = 1; j <= jsize; j++){
  for(int i = 1; i <= isize; i++){
    const double den_y0 = (
        + 0.5 * DEN(i  , j-1)
        + 0.5 * DEN(i  , j  )
    );
    const double kappa = (
        + 0.5 * CURV(i  , j-1)
        + 0.5 * CURV(i  , j  )
    );
    const double delta = 1. / hy * (
        - VOF(i  , j-1)
        + VOF(i  , j  )
    );
    IFRCY(i, j) = den_y0 * refdeninv * tension * kappa * delta;
  }
}

src/interface/force.c

for(int k = 1; k <= ksize; k++){
  for(int j = 1; j <= jsize; j++){
    for(int i = 1; i <= isize; i++){
      const double den_y0 = (
          + 0.5 * DEN(i  , j-1, k  )
          + 0.5 * DEN(i  , j  , k  )
      );
      const double kappa = (
          + 0.5 * CURV(i  , j-1, k  )
          + 0.5 * CURV(i  , j  , k  )
      );
      const double delta = 1. / hy * (
          - VOF(i  , j-1, k  )
          + VOF(i  , j  , k  )
      );
      IFRCY(i, j, k) = den_y0 * refdeninv * tension * kappa * delta;
    }
  }
}

src/interface/force.c

for(int k = 1; k <= ksize; k++){
  for(int j = 1; j <= jsize; j++){
    for(int i = 1; i <= isize; i++){
      const double den_z0 = (
          + 0.5 * DEN(i  , j  , k-1)
          + 0.5 * DEN(i  , j  , k  )
      );
      const double kappa = (
          + 0.5 * CURV(i  , j  , k-1)
          + 0.5 * CURV(i  , j  , k  )
      );
      const double delta = 1. / hz * (
          - VOF(i  , j  , k-1)
          + VOF(i  , j  , k  )
      );
      IFRCZ(i, j, k) = den_z0 * refdeninv * tension * kappa * delta;
    }
  }
}

The pre-factor is computed here:

src/interface/force.c

static double compute_refdeninv(
    const double denr
){
  return 2. / (1. + denr);
}