Two-scale isothermal atomization

An air jet shears a liquid column that atomizes, with inter-scale mass transfer feeding a disperse phase.

Two-scale isothermal atomization

A liquid column is suddenly exposed to a fast air stream. The shear tears the column apart, stretching it into ligaments that break into a fine spray. This is the separate-to-disperse transition at the heart of primary atomization.

The case is the air-blasted liquid column of Section 5.2 of the reference paper, computed with a unified two-scale two-phase model: the same set of equations describes both the resolved (large-scale) interface and an unresolved (small-scale) disperse phase, and an inter-scale mass transfer moves liquid from one representation to the other when the interface curvature exceeds a threshold. This transfer is the main novelty of the model.

Model

Both phases are isothermal and governed by a barotropic (linearized) equation of state. The large-scale liquid is tracked through its volume fraction $\alpha_\ell$ , transported with the mixture and coupled to surface tension through the capillary term $\sigma\,\kappa\,\mathbf{n}$ , where the curvature $\kappa$ and interface normal $\mathbf{n}$ are reconstructed from $\nabla\alpha_\ell$ :

\partial_t (\rho \mathbf{u}) + \nabla\cdot(\rho\,\mathbf{u}\otimes\mathbf{u}) + \nabla p = \nabla\cdot\big(\sigma\,(\|\nabla\alpha_\ell\|\,\mathbf{I} - \tfrac{\nabla\alpha_\ell\otimes\nabla\alpha_\ell}{\|\nabla\alpha_\ell\|})\big).

Where the interface becomes too curved to resolve, a bound-preserving relaxation transfers mass into a small-scale disperse volume fraction $\alpha_d$ , closing the separate-to-disperse cycle.

Numerical method

Space: second-order finite volumes; an HLLC solver for the hyperbolic subsystem and a dedicated flux for the surface-tension subsystem.
Time: an operator-splitting, second-order (two-stage) scheme.
Relaxation: a Newton-based, bound-preserving update of $\alpha_\ell$ that enforces the admissible-state bounds and drives the inter-scale mass transfer.
Adaptation: samurai multiresolution refines the mesh along the interface and the freshly created ligaments.

What to look for

Watch the liquid column (bright) flatten, wrap into a horseshoe and shed ligaments downstream. The multiresolution mesh (grid overlay) chases the interface, concentrating cells exactly where atomization happens and coarsening in the quiescent gas.

References

See case.yaml for the reference list.

References

Source code two_scale_capillarity.cpp

Powered by the 2026_06_two_scale_isothermal solver engine - the code below is this case's scenario (initial & boundary conditions); the flux, time integration and mesh adaptation come from the shared engine.

two_scale_capillarity.cpp 180 lines

// Copyright 2021 SAMURAI TEAM. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
//
// Author: Giuseppe Orlando, 2026
//
#include <CLI/CLI.hpp>

#include <nlohmann/json.hpp>

#include "include/two_scale_capillarity.hpp"

// Main function to run the program
//
int main(int argc, char* argv[]) {
  using json = nlohmann::json;

  auto& app = samurai::initialize("Finite volume example for the air-blasted liquid column configuration", argc, argv);

  json input;
  try {
    std::ifstream ifs("input.json"); // Read a JSON file
    input = json::parse(ifs);
  }
  catch(const json::parse_error& e) {
    throw std::runtime_error("Cannot parse parameter file 'input.json'. Please verify that the file is present");
  }

  /*--- Set and declare simulation parameters ---*/
  using Number = TwoScaleCapillarity<EquationData::dim>::Number;
  Simulation_Parameters<Number> sim_param;

  // Physical parameters
  sim_param.xL = input.value("xL", static_cast<double>(0.0));
  sim_param.xR = input.value("xR", static_cast<double>(4.0));
  sim_param.yL = input.value("yL", static_cast<double>(0.0));
  sim_param.yR = input.value("yR", static_cast<double>(2.0));

  sim_param.t0 = input.value("t0", static_cast<Number>(0.0));
  sim_param.Tf = input.value("Tf", static_cast<Number>(2.5));

  sim_param.sigma = input.value("sigma", static_cast<Number>(1e-2));

  sim_param.apply_relaxation = input.value("apply_relaxation", true);
  sim_param.mass_transfer    = input.value("mass_transfer", false);
  sim_param.Hmax             = input.value("Hmax", static_cast<Number>(40.0));
  sim_param.kappa            = input.value("kappa", static_cast<Number>(1.0));
  sim_param.alpha_d_max      = input.value("alpha_d_max", static_cast<Number>(0.5));
  sim_param.alpha_l_min      = input.value("alpha_l_min", static_cast<Number>(0.01));
  sim_param.alpha_l_max      = input.value("alpha_l_max", static_cast<Number>(0.1));

  sim_param.x0         = input.value("x0", static_cast<Number>(1.0));
  sim_param.y0         = input.value("y0", static_cast<Number>(1.0));
  sim_param.U0         = input.value("U0", static_cast<Number>(6.66));
  sim_param.U1         = input.value("U1", static_cast<Number>(0.0));
  sim_param.V0         = input.value("V0", static_cast<Number>(0.0));
  sim_param.R          = input.value("R", static_cast<Number>(0.15));
  sim_param.eps_over_R = input.value("eps_over_R", static_cast<Number>(0.2));

  // Numerical parameters
  sim_param.Courant = input.value("cfl", static_cast<Number>(0.4));

  sim_param.alpha_residual       = input.value("alpha_residual", static_cast<Number>(1e-8));
  sim_param.mod_grad_alpha_l_min = input.value("mod_grad_alpha_l_min", static_cast<Number>(0.0));

  sim_param.lambda           = input.value("lambda", static_cast<Number>(0.9));
  sim_param.atol_Newton      = input.value("atol_Newton", static_cast<Number>(1e-12));
  sim_param.rtol_Newton      = input.value("rtol_Newton", static_cast<Number>(1e-10));
  sim_param.max_Newton_iters = input.value("max_Newton_iters", static_cast<std::size_t>(60));

  // MR paramters
  sim_param.min_level     = input.value("min-level", static_cast<std::size_t>(8));
  sim_param.max_level     = input.value("max-level", static_cast<std::size_t>(8));
  sim_param.MR_param      = input.value("MR_param", static_cast<double>(1e-1));
  sim_param.MR_regularity = input.value("MR_regularity", static_cast<double>(0));

  // Output parameters
  sim_param.save_dir = input.value("save-dir", fs::current_path());
  sim_param.nfiles   = input.value("nfiles", static_cast<std::size_t>(10));

  // Restart file
  sim_param.restart_file = input.value("restart_file","");

  /*--- Allow for parsing from command line ---*/
  // Physical parameters
  app.add_option("--xL", sim_param.xL, "x Left-end of the domain")->capture_default_str()->group("Physical parameters");
  app.add_option("--xR", sim_param.xR, "x Right-end of the domain")->capture_default_str()->group("Physical parameters");
  app.add_option("--yL", sim_param.yL, "y Bottom-end of the domain")->capture_default_str()->group("Physical parameters");
  app.add_option("--yR", sim_param.yR, "y Top-end of the domain")->capture_default_str()->group("Physical parameters");

  app.add_option("--t0", sim_param.t0, "Initial time")->capture_default_str()->group("Physical parameters");
  app.add_option("--Tf", sim_param.Tf, "Final time")->capture_default_str()->group("Physical parameters");

  app.add_option("--sigma", sim_param.sigma, "Surface tension coefficient")->capture_default_str()->group("Physical parameters");

  app.add_option("--apply_relaxation", sim_param.apply_relaxation, "Apply or not relaxation")->capture_default_str()->group("Physical parameters");
  app.add_option("--mass_transfer", sim_param.mass_transfer,
                 "Choose whether to perform or not the mass transfer")->capture_default_str()->group("Physical parameters");
  app.add_option("--kappa", sim_param.kappa,
                 "Small-scale disperse phase raidus with rispect to maximum curvature")->capture_default_str()->group("Physical parameters");
  app.add_option("--Hmax", sim_param.Hmax,
                 "Maximum curvature before activating atomization")->capture_default_str()->group("Physical parameters");
  app.add_option("--alpha_d_max", sim_param.alpha_d_max,
                 "Maximum admitted small-scale volume fraction")->capture_default_str()->group("Physical parameters");
  app.add_option("--alpha_l_min", sim_param.alpha_l_min,
                 "Maximum effective volume fraction for the mixture region")->capture_default_str()->group("Physical parameters");
  app.add_option("--alpha_l_max", sim_param.alpha_l_max,
                 "Maximum effective volume fraction for the mixture region")->capture_default_str()->group("Physical parameters");

  app.add_option("--x0", sim_param.x0, "Liquid column x-center")->capture_default_str()->group("Physical parameters");
  app.add_option("--y0", sim_param.y0, "Liquid column y-center")->capture_default_str()->group("Physical parameters");
  app.add_option("--U0", sim_param.U0, "Parameter for initial horizontal velocity")->capture_default_str()->group("Physical parameters");
  app.add_option("--U1", sim_param.U1, "Parameter for initial horizontal velocity")->capture_default_str()->group("Physical parameters");
  app.add_option("--V0", sim_param.V0, "Initial vertical velocity")->capture_default_str()->group("Physical parameters");
  app.add_option("--R", sim_param.R, "Initial radius of the liquid column")->capture_default_str()->group("Physical parameters");
  app.add_option("--eps_over_R", sim_param.eps_over_R,
                 "Initial interface thickness with respect to the radius")->capture_default_str()->group("Physical parameters");

  // Numerical parameters
  app.add_option("--cfl", sim_param.Courant, "The Courant number")->capture_default_str()->group("Numerical parameters");

  app.add_option("--alpha_residual", sim_param.alpha_residual, "Residual large scale volume fraction")->capture_default_str()->group("Numerical parameters");
  app.add_option("--mod_grad_alpha_l_min", sim_param.mod_grad_alpha_l_min,
                 "Tolerance for zero gradient volume fraction")->capture_default_str()->group("Numerical parameters");

  app.add_option("--lambda", sim_param.lambda,
                 "Parameter for bound-preserving strategy")->capture_default_str()->group("Numerical parameters");
  app.add_option("--atol_Newton", sim_param.atol_Newton,
                 "Absolute tolerance of Newton method for the relaxation")->capture_default_str()->group("Numerical parameters");
  app.add_option("--rtol_Newton", sim_param.rtol_Newton,
                 "Relative tolerance of Newton method for the relaxation")->capture_default_str()->group("Numerical parameters");
  app.add_option("--max_Newton_iters", sim_param.max_Newton_iters,
                 "Maximum number of Newton iterations")->capture_default_str()->group("Numerical parameters");

  // MR parameters
  app.add_option("--min-level", sim_param.min_level, "Minimum level of the AMR")->capture_default_str()->group("AMR parameter");
  app.add_option("--max-level", sim_param.max_level, "Maximum level of the AMR")->capture_default_str()->group("AMR parameter");
  app.add_option("--MR_param", sim_param.MR_param, "Multiresolution parameter")->capture_default_str()->group("AMR parameter");
  app.add_option("--MR_regularity", sim_param.MR_regularity, "Multiresolution regularity")->capture_default_str()->group("AMR parameter");

  // Output parameters
  app.add_option("--save-dir", sim_param.save_dir, "Output directory")->capture_default_str()->group("Output parameters");
  app.add_option("--nfiles", sim_param.nfiles, "Number of output files")->capture_default_str()->group("Output parameters");

  // Restart file
  app.add_option("--restart_file", sim_param.restart_file, "Name of the restart file")->capture_default_str()->group("Restart");

  /*--- Set and declare simulation parameters related to EOS ---*/
  EOS_Parameters<Number> eos_param;

  eos_param.p0_phase_liq   = input.value("p0_phase_liq", static_cast<Number>(1e5));
  eos_param.rho0_phase_liq = input.value("rho0_phase_liq", static_cast<Number>(1e3));
  eos_param.c0_phase_liq   = input.value("c0_phase_liq", static_cast<Number>(1e1));

  eos_param.p0_phase_gas   = input.value("p0_phase_gas", static_cast<Number>(1e5));
  eos_param.rho0_phase_gas = input.value("rho0_phase_gas", static_cast<Number>(1.0));
  eos_param.c0_phase_gas   = input.value("c0_phase_gas", static_cast<Number>(1e1));

  app.add_option("--p0_phase_liq", eos_param.p0_phase_liq, "p0_phase_liq")->capture_default_str()->group("EOS parameters");
  app.add_option("--rho0_phase_liq", eos_param.rho0_phase_liq, "rho0_phase_liq")->capture_default_str()->group("EOS parameters");
  app.add_option("--c0_phase_liq", eos_param.c0_phase_liq, "c0_phase_liq")->capture_default_str()->group("EOS parameters");
  app.add_option("--p0_phase_gas", eos_param.p0_phase_gas, "p0_phase_gas")->capture_default_str()->group("EOS parameters");
  app.add_option("--rho0_phase_gas", eos_param.rho0_phase_gas, "rho0_phase_gas")->capture_default_str()->group("EOS parameters");
  app.add_option("--c0_phase_gas", eos_param.c0_phase_gas, "c0_phase_gas")->capture_default_str()->group("EOS parameters");

  /*--- Create the instance of the class to perform the simulation ---*/
  CLI11_PARSE(app, argc, argv);

  xt::xtensor_fixed<double, xt::xshape<EquationData::dim>> min_corner = {sim_param.xL, sim_param.yL};
  xt::xtensor_fixed<double, xt::xshape<EquationData::dim>> max_corner = {sim_param.xR, sim_param.yR};
  auto TwoScaleCapillarity_Sim = TwoScaleCapillarity(min_corner, max_corner,
                                                     sim_param, eos_param);

  TwoScaleCapillarity_Sim.run(sim_param.nfiles);

  samurai::finalize();

  return 0;
}