Coupled Level-Set LBM on Adaptive Grids

• The Coupled Level-Set Lattice Boltzmann Method is a computational framework that uses separate solvers for liquid and gas phases with a high-order level-set for sharp interface tracking.
• It rigorously enforces jump conditions at the interface, achieving high-fidelity results in benchmarks such as stratified flows and rising bubble tests.
• Adaptive mesh refinement via parallel quad/octree grids concentrates computational effort near interfaces, enhancing accuracy and efficiency without extra smoothing.

The coupled level-set lattice Boltzmann method (LS-LBM) on adaptive Cartesian grids is a computational framework for simulating liquid–gas multiphase flows characterized by sharp interfaces and large density and viscosity ratios. This approach employs separate lattice Boltzmann solvers for each fluid phase, coupled via sharp-interface boundary conditions, and couples them with a high-order level-set method to track the interface position. Adaptive mesh refinement (AMR) using a parallel quad/octree structure concentrates computational effort near interfacial regions, thereby improving accuracy and efficiency. The method facilitates the rigorous enforcement of interface jump conditions, treatment of complex interfacial dynamics, and has demonstrated high-fidelity results across standard test cases involving stratified flows and rising bubbles (Vorspohl et al., 9 Jan 2026).

1. Governing Equations and Discretization

The LS-LBM framework employs distinct lattice Boltzmann solvers for each fluid phase, liquid and gas, each solving the discrete Boltzmann equation using a collision-streaming splitting. In each time step $\Delta t$, the particle distribution function $f_i(\mathbf{x}, t)$ undergoes collision

$$f^{c}_i(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\{f\}) + F_i$$

and then streaming

$$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f^c_i(\mathbf{x}, t).$$

Macroscopic quantities are computed by

$$\rho = \sum_i f_i,\quad \mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i.$$

For the liquid phase (low Reynolds number), the BGK collision operator is used:

$\Omega_i(f) = \omega_{BGK}(f^{eq}_i - f_i),$

where the Maxwell equilibrium is

$$f_i^{eq} = w_i \rho\Big[1 + \frac{\mathbf{c}_i\cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2c_s^2}\Big],\quad \omega_{BGK} = \frac{c_s^2}{\nu + \frac{1}{2}c_s^2}.$$

For the gas phase (high Reynolds number, 3D), the cumulant collision operator is employed:

• Transform $f_i$ to cumulants $c_\alpha$.
• Relax each cumulant:

$$c^c_\alpha = c_\alpha + \omega_\alpha (c_\alpha^{eq} - c_\alpha),$$

with $f_i(\mathbf{x}, t)$0 for $f_i(\mathbf{x}, t)$1, which enhances stability at high Re.

External forces such as gravity are incorporated via

$f_i(\mathbf{x}, t)$2

The interface is represented by a signed–distance function $f_i(\mathbf{x}, t)$3, with $f_i(\mathbf{x}, t)$4 in the liquid, $f_i(\mathbf{x}, t)$5 in the gas, and $f_i(\mathbf{x}, t)$6 at the interface $f_i(\mathbf{x}, t)$7. Its advection follows:

$f_i(\mathbf{x}, t)$8

discretized using a 5th-order upwind-central scheme in space and a 3rd-order TVD–RK in time. The velocity $f_i(\mathbf{x}, t)$9 at the interface is extended off $$f^{c}_i(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\{f\}) + F_i$$0 by solving a hyperbolic equation in artificial time $$f^{c}_i(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\{f\}) + F_i$$1:

$$f^{c}_i(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\{f\}) + F_i$$2

where $$f^{c}_i(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\{f\}) + F_i$$3.

To maintain $$f^{c}_i(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\{f\}) + F_i$$4, $$f^{c}_i(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\{f\}) + F_i$$5 is reinitialized using a constrained Hamilton–Jacobi reinitialization, which preserves the zero level set exactly.

2. Coupling Strategy at the Sharp Interface

Enforcement of velocity and stress jump conditions at $$f^{c}_i(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\{f\}) + F_i$$6 uses a modified two-fluid bounce–back method analogous to Bouzidi et al. In cut cells (intersected by the interface), streaming into the opposite phase is disallowed; instead, the missing post-collision $$f^{c}_i(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\{f\}) + F_i$$7 is constructed as:

$$f^{c}_i(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\{f\}) + F_i$$8

with $$f^{c}_i(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\{f\}) + F_i$$9 the linearly interpolated interface velocity (ensuring $$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f^c_i(\mathbf{x}, t).$$0), $$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f^c_i(\mathbf{x}, t).$$1 indicating the interface position along $$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f^c_i(\mathbf{x}, t).$$2, $$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f^c_i(\mathbf{x}, t).$$3, and $$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f^c_i(\mathbf{x}, t).$$4 the total jump in viscous stress, pressure, and capillary stress.

The stress jump $$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f^c_i(\mathbf{x}, t).$$5 decomposes as: \begin{align*} [S] : \mathbf{n}\otimes\mathbf{n} &= \frac{[p] + 2\sigma\kappa}{2\bar{\eta}} - \frac{[\eta]}{\bar{\eta} : \mathbf{n}\otimes\mathbf{n}}\ [S] : \mathbf{n}\otimes\mathbf{t}_j &= - \frac{[\eta]}{\bar{\eta} : \mathbf{n}\otimes\mathbf{t}_j} \end{align*} where $$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f^c_i(\mathbf{x}, t).$$6 is the pressure jump, $$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f^c_i(\mathbf{x}, t).$$7 the viscosity jump, $$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f^c_i(\mathbf{x}, t).$$8, and $$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f^c_i(\mathbf{x}, t).$$9 tangential directions.

Surface tension can be included via a continuum surface force:

$$\rho = \sum_i f_i,\quad \mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i.$$0

where $$\rho = \sum_i f_i,\quad \mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i.$$1 is surface tension, $$\rho = \sum_i f_i,\quad \mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i.$$2 the curvature, and $$\rho = \sum_i f_i,\quad \mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i.$$3 a regularized Heaviside function; $$\rho = \sum_i f_i,\quad \mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i.$$4 enters the forcing term during collision.

3. Adaptive Cartesian Grid Infrastructure

The computational domain is implemented as a parallel quad-tree or octree ("forest-of-trees"), shared by the phase-resolved solvers (liquid LBM, gas LBM, level-set). This eliminates the need for inter-solver data communication, as all fields reside in the same memory.

Refinement and coarsening are governed by solver-specific sensor functions $$\rho = \sum_i f_i,\quad \mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i.$$5. For instance, near the interface, refinement is triggered if $$\rho = \sum_i f_i,\quad \mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i.$$6, with coarsening if $$\rho = \sum_i f_i,\quad \mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i.$$7 and the local level equals $$\rho = \sum_i f_i,\quad \mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i.$$8.

Interpolation of missing $$\rho = \sum_i f_i,\quad \mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i.$$9 values across grid refinement jumps uses the Dupuis & Chopard strategy: spatial and temporal interpolation from adjacent levels, scaling $\Omega_i(f) = \omega_{BGK}(f^{eq}_i - f_i),$0 to ensure $\Omega_i(f) = \omega_{BGK}(f^{eq}_i - f_i),$1 across all levels. When the phase changes (as $\Omega_i(f) = \omega_{BGK}(f^{eq}_i - f_i),$2 sweeps a cell), the cell's $\Omega_i(f) = \omega_{BGK}(f^{eq}_i - f_i),$3 are refilled from neighbors along a lattice direction $\Omega_i(f) = \omega_{BGK}(f^{eq}_i - f_i),$4 closest to $\Omega_i(f) = \omega_{BGK}(f^{eq}_i - f_i),$5, and the equilibrium is reconstructed by three-point interpolation using interface and neighbor values.

4. Computational Workflow and Implementation

The main computational loop per global time step $\Omega_i(f) = \omega_{BGK}(f^{eq}_i - f_i),$6 proceeds as follows:

1. For each phase-LBM (liquid and gas): a) Collision: $\Omega_i(f) = \omega_{BGK}(f^{eq}_i - f_i),$7 (including body forces and, optionally, surface tension). b) Apply bounce–back for cut cells using the $\Omega_i(f) = \omega_{BGK}(f^{eq}_i - f_i),$8 formula. c) Streaming: propagate $\Omega_i(f) = \omega_{BGK}(f^{eq}_i - f_i),$9. d) Update macroscopic moments ($$f_i^{eq} = w_i \rho\Big[1 + \frac{\mathbf{c}_i\cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2c_s^2}\Big],\quad \omega_{BGK} = \frac{c_s^2}{\nu + \frac{1}{2}c_s^2}.$$0, $$f_i^{eq} = w_i \rho\Big[1 + \frac{\mathbf{c}_i\cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2c_s^2}\Big],\quad \omega_{BGK} = \frac{c_s^2}{\nu + \frac{1}{2}c_s^2}.$$1).
2. Extend interface velocity off $$f_i^{eq} = w_i \rho\Big[1 + \frac{\mathbf{c}_i\cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2c_s^2}\Big],\quad \omega_{BGK} = \frac{c_s^2}{\nu + \frac{1}{2}c_s^2}.$$2 by solving the hyperbolic extension equation in $$f_i^{eq} = w_i \rho\Big[1 + \frac{\mathbf{c}_i\cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2c_s^2}\Big],\quad \omega_{BGK} = \frac{c_s^2}{\nu + \frac{1}{2}c_s^2}.$$3.
3. Advect the level-set function: update $$f_i^{eq} = w_i \rho\Big[1 + \frac{\mathbf{c}_i\cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2c_s^2}\Big],\quad \omega_{BGK} = \frac{c_s^2}{\nu + \frac{1}{2}c_s^2}.$$4 via the advection equation (5th-order WENO in space, 3rd-order RK in time).
4. Reinitialize $$f_i^{eq} = w_i \rho\Big[1 + \frac{\mathbf{c}_i\cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2c_s^2}\Big],\quad \omega_{BGK} = \frac{c_s^2}{\nu + \frac{1}{2}c_s^2}.$$5 to maintain the signed-distance property (constrained Hamilton–Jacobi).
5. Compute interface normals $$f_i^{eq} = w_i \rho\Big[1 + \frac{\mathbf{c}_i\cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2c_s^2}\Big],\quad \omega_{BGK} = \frac{c_s^2}{\nu + \frac{1}{2}c_s^2}.$$6 and curvatures $$f_i^{eq} = w_i \rho\Big[1 + \frac{\mathbf{c}_i\cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2c_s^2}\Big],\quad \omega_{BGK} = \frac{c_s^2}{\nu + \frac{1}{2}c_s^2}.$$7 in cut cells.
6. Optionally, assemble surface tension force $$f_i^{eq} = w_i \rho\Big[1 + \frac{\mathbf{c}_i\cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2c_s^2}\Big],\quad \omega_{BGK} = \frac{c_s^2}{\nu + \frac{1}{2}c_s^2}.$$8 for the next LBM forcing step.
7. Evaluate sensor functions and tag cells for refinement or coarsening.
8. Execute AMR: refine or coarsen the quad/octree structure, redistribute data, and refill newly created phase cells using the described non-equilibrium refilling strategy.
9. Advance to the next time step.

Special handling includes the use of cumulant collision for high-Re gas flows, modified bounce–back for sharp two-fluid coupling, and the lack of additional smoothing or filtering required for large density ratios.

5. Validation and Test Cases

Validation covers canonical problems and benchmarks:

• Laminar Stratified Couette and Poiseuille Flows (2D): Two immiscible fluids occupy $$f_i^{eq} = w_i \rho\Big[1 + \frac{\mathbf{c}_i\cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2c_s^2}\Big],\quad \omega_{BGK} = \frac{c_s^2}{\nu + \frac{1}{2}c_s^2}.$$9 and $f_i$0, respectively. For viscosity ratio $f_i$1, the velocity profile’s $f_i$2 error is $f_i$3, matching analytical solutions.
• Single Rising Bubble in 3D (Safi et al. benchmark):
  • Case I: $f_i$4, $f_i$5, $f_i$6, $f_i$7
  • Case II: $f_i$8, $f_i$9, $c_\alpha$0, $c_\alpha$1
  • On a uniform grid with $c_\alpha$2, the liquid phase uses BGK with $c_\alpha$3; the gas uses the cumulant operator. The bubble center-of-mass velocity deviates by $c_\alpha$4 from reference; shape and interface deformation are well captured, including pressure jumps and vorticity structure.
• Adaptive Mesh Refinement Study (Case II):

Starting from a coarse base ($c_\alpha$5), local refinement ($c_\alpha$6) within $c_\alpha$7. With a narrow refinement band ($c_\alpha$8), rise velocity error is $c_\alpha$9; for $$c^c_\alpha = c_\alpha + \omega_\alpha (c_\alpha^{eq} - c_\alpha),$$0, error falls below $$c^c_\alpha = c_\alpha + \omega_\alpha (c_\alpha^{eq} - c_\alpha),$$1. Fine cells are confined to $$c^c_\alpha = c_\alpha + \omega_\alpha (c_\alpha^{eq} - c_\alpha),$$2–$$c^c_\alpha = c_\alpha + \omega_\alpha (c_\alpha^{eq} - c_\alpha),$$3 of the domain, yielding substantial resource savings.

• Cluster of Nine Rising Bubbles (3D):

Two staggered layers, domain size $$c^c_\alpha = c_\alpha + \omega_\alpha (c_\alpha^{eq} - c_\alpha),$$4. Wake interaction, downwash, and bubble–bubble coupling are faithfully resolved, with no observed artifacts from mesh modification.

In summary, the method achieves second-order spatial accuracy in stratified flow tests, robustly enforces interface jump conditions, resolves interfacial structure and flow features without the need for smoothing at large density ratios, and realizes significant efficiency gains from AMR (Vorspohl et al., 9 Jan 2026).

6. Significance and Context

The coupled level-set LBM on adaptive Cartesian grids enables accurate direct simulation of sharp liquid–gas interfaces with arbitrary density and viscosity ratios, leveraging separate solution strategies for each phase and coupling via sharp-interface boundary conditions. The composite quad-/octree AMR infrastructure focuses computational effort, providing benchmark-level accuracy for classical multifluid flow problems while minimizing computational cost. No additional smoothing or filtering is required to stabilize large-density ratio flows, suggesting robustness for a broad class of high-fidelity multiphase flow simulations.

A plausible implication is the method’s suitability for industrial or geophysical multiphase flows where sharp interface dynamics, large density contrasts, and complex domain topologies are present. The hybrid approach provides a compelling alternative to diffuse interface or single-fluid methods when precise handling of interface physics is necessary (Vorspohl et al., 9 Jan 2026).