LCRM: Hybrid Interface Tracking for Multiphase Flow

• LCRM is a hybrid interface tracking method for 3D multiphase flows that integrates Lagrangian and Eulerian approaches for improved accuracy and mass conservation.
• It leverages a triangulated Lagrangian interface and periodic Eulerian reconstruction to compute sharp curvature and manage topology changes efficiently.
• The method effectively handles surfactant transport and surface-viscous stresses, demonstrating strong convergence properties in various CFD benchmarks.

Searching arXiv for Level Contour Reconstruction Method and related papers to ground the article. Found relevant papers:

• (Panda et al., 29 Sep 2025) — "Direct numerical simulation of two-phase flows with surfactant-induced surface viscous effects"
• (Shin et al., 2017) — "An Interface-Tracking Technique for Multiphase Flow with Soluble Surfactant"
• (0901.1773) — likely Shin et al. IJNMF 2009 LCRM-related preprint/result

The Level Contour Reconstruction Method (LCRM) is an interface–tracking scheme for three-dimensional multiphase flow that combines an explicit Lagrangian interface with an Eulerian distance function used for periodic interface reconstruction. In the cited formulations, LCRM is presented as a hybrid front-tracking/level-set approach that retains sharp interfaces, accurate curvature, and good mass conservation while handling breakup and coalescence without explicit connectivity bookkeeping. It has been used for surface-tension-dominated multiphase flow, insoluble and soluble surfactant transport, Marangoni stresses, and surfactant-induced surface viscous stresses within the Boussinesq–Scriven constitutive framework (Shin et al., 2017, Panda et al., 29 Sep 2025).

1. Hybrid formulation and defining characteristics

LCRM was originally developed by Shin & Juric as a hybrid scheme that keeps an explicit Lagrangian interface but periodically reconstructs that interface from a distance or level-set field defined on the Eulerian grid. In the descriptions given for both soluble-surfactant and surface-viscous flow formulations, the method is neither a classical front-tracking algorithm with persistent global mesh connectivity nor a pure level-set method in which a signed-distance field is advected and reinitialized. Instead, the interface is advected explicitly as discrete Lagrangian triangular elements, and, when needed, a signed distance function is reconstructed from that interface and used to rebuild the interface as the zero level contour of the distance field in each cell (Shin et al., 2017, Panda et al., 29 Sep 2025).

This hybridization addresses two persistent difficulties in three-dimensional interfacial CFD. Classical front tracking gives sharp interfaces and exact geometric control, but mesh-quality operations, remeshing, and connectivity management become difficult under large deformation and especially under breakup or coalescence. Pure level-set methods handle topology change elegantly and provide smooth normals and curvature, but they are susceptible to mass loss and represent the interface over several cells. LCRM is designed to avoid the connectivity management and mesh-quality operations of classical front tracking in 3D while also avoiding the mass loss and curvature errors inherent to pure level-set advection (Shin et al., 2017).

A common misconception is to identify LCRM with a standard level-set method. The cited formulations explicitly distinguish the two: the distance function is computed directly from the Lagrangian interface, and, in the surface-viscous implementation, the level-set role is geometric and topological rather than kinematic. The interface kinematics remain purely Lagrangian, with marker motion governed by

$$\frac{d\mathbf x_f}{dt}=\mathbf u_f=\mathbf u(\mathbf x=\mathbf x_f),$$

integrated with a second-order Runge–Kutta scheme (Panda et al., 29 Sep 2025).

2. Geometric representation, reconstruction, and surface operators

In three dimensions, the interface is represented by a triangulated Lagrangian surface. Each triangular element $e$ is defined by three vertices, with area and normal given by

$$A_e=\frac{1}{2}\,|\mathbf x_{02}\times \mathbf x_{03}|,\qquad \mathbf n_f=\frac{\mathbf x_{02}\times \mathbf x_{03}}{|\mathbf x_{02}\times \mathbf x_{03}|}.$$

In the soluble-surfactant formulation, the interfacial surfactant concentration $\Gamma$ is stored at the center of each element; in the surface-viscous formulation, the surface carries $\Gamma$, $\sigma(\Gamma)$, $\mu_s(\Gamma)$, $\mu_d(\Gamma)$, and related interfacial stresses (Shin et al., 2017, Panda et al., 29 Sep 2025).

The Eulerian component of LCRM is a signed distance field, denoted $\varphi(\mathbf x)$ or $\phi(\mathbf x)$, constructed directly from the piecewise planar Lagrangian interface rather than by solving a transport PDE. Reconstruction is performed cell-by-cell. Each Eulerian cell is subdivided into tetrahedra, and the zero isocontour inside each tetrahedron is intersected to create local triangles. The cited descriptions emphasize a tetra-marching procedure that ensures a unique, non-ambiguous reconstruction in each cell. In the surface-viscous implementation, reconstruction is performed typically every 25 time steps, with the frequency chosen such that the displacement of markers per step does not exceed the smallest grid spacing. This regular rebuilding keeps elements well-shaped, avoids stretching and clustering, and permits topology change to be inherited from the distance field (Shin et al., 2017, Panda et al., 29 Sep 2025).

Although curvature is computed from the Lagrangian surface, normals for several local operations are often obtained from the distance function,

$e$0

At edge centers, this normal is used to define the binormal direction $e$1 and tangential direction $e$2, which enter the computation of surface gradients of $e$3, $e$4, and, in the Boussinesq–Scriven implementation, velocity-derived surface tensors. This geometric split—Lagrangian curvature and area, Eulerian normals and local orientation—is one of the characteristic features of LCRM (Shin et al., 2017).

3. One-fluid formulation, surfactant transport, and Marangoni coupling

The flow is formulated on a fixed Eulerian grid using a one-fluid or single-field incompressible Navier–Stokes system with discontinuous properties. In the surface-viscous formulation,

$e$5

while the soluble-surfactant formulation uses an indicator $e$6 in an equivalent single-field representation. Incompressibility is enforced through

$e$7

Interfacial forces enter the momentum equation through sharply localized force densities distributed from the Lagrangian interface to the Eulerian grid by a discrete Peskin-type kernel (Shin et al., 2017, Panda et al., 29 Sep 2025).

For surface tension and Marangoni stresses, the interfacial force is written as the sum of a normal curvature term and a tangential term,

$e$8

The normal part is implemented through a hybrid curvature formulation that aligns the discrete curvature force with the pressure gradient operator and reduces parasitic currents. The tangential part is assembled edge-by-edge after decomposing $e$9 into binormal and tangential components (Shin et al., 2017).

In the insoluble-surfactant formulation, the surface concentration satisfies

$$A_e=\frac{1}{2}\,|\mathbf x_{02}\times \mathbf x_{03}|,\qquad \mathbf n_f=\frac{\mathbf x_{02}\times \mathbf x_{03}}{|\mathbf x_{02}\times \mathbf x_{03}|}.$$0

with the surface material derivative defined by

$$A_e=\frac{1}{2}\,|\mathbf x_{02}\times \mathbf x_{03}|,\qquad \mathbf n_f=\frac{\mathbf x_{02}\times \mathbf x_{03}}{|\mathbf x_{02}\times \mathbf x_{03}|}.$$1

In the soluble-surfactant extension, the interfacial conservation law is written as

$$A_e=\frac{1}{2}\,|\mathbf x_{02}\times \mathbf x_{03}|,\qquad \mathbf n_f=\frac{\mathbf x_{02}\times \mathbf x_{03}}{|\mathbf x_{02}\times \mathbf x_{03}|}.$$2

where the adsorption/desorption source is

$$A_e=\frac{1}{2}\,|\mathbf x_{02}\times \mathbf x_{03}|,\qquad \mathbf n_f=\frac{\mathbf x_{02}\times \mathbf x_{03}}{|\mathbf x_{02}\times \mathbf x_{03}|}.$$3

Bulk concentration in one phase evolves through

$$A_e=\frac{1}{2}\,|\mathbf x_{02}\times \mathbf x_{03}|,\qquad \mathbf n_f=\frac{\mathbf x_{02}\times \mathbf x_{03}}{|\mathbf x_{02}\times \mathbf x_{03}|}.$$4

and is coupled to the interface by the sharp Neumann condition

$$A_e=\frac{1}{2}\,|\mathbf x_{02}\times \mathbf x_{03}|,\qquad \mathbf n_f=\frac{\mathbf x_{02}\times \mathbf x_{03}}{|\mathbf x_{02}\times \mathbf x_{03}|}.$$5

In practice this boundary condition is enforced through a ghost-cell method adapted from the “sharp energy” method of Shin & Choi (Shin et al., 2017).

The surface tension is surfactant-dependent. The cited formulations use a Langmuir–Szyszkowski equation of state,

$$A_e=\frac{1}{2}\,|\mathbf x_{02}\times \mathbf x_{03}|,\qquad \mathbf n_f=\frac{\mathbf x_{02}\times \mathbf x_{03}}{|\mathbf x_{02}\times \mathbf x_{03}|}.$$6

and, in some soluble-surfactant benchmarks, a linear equation of state,

$$A_e=\frac{1}{2}\,|\mathbf x_{02}\times \mathbf x_{03}|,\qquad \mathbf n_f=\frac{\mathbf x_{02}\times \mathbf x_{03}}{|\mathbf x_{02}\times \mathbf x_{03}|}.$$7

Marangoni stresses then arise from

$$A_e=\frac{1}{2}\,|\mathbf x_{02}\times \mathbf x_{03}|,\qquad \mathbf n_f=\frac{\mathbf x_{02}\times \mathbf x_{03}}{|\mathbf x_{02}\times \mathbf x_{03}|}.$$8

A notable numerical feature is that surface gradients of $$A_e=\frac{1}{2}\,|\mathbf x_{02}\times \mathbf x_{03}|,\qquad \mathbf n_f=\frac{\mathbf x_{02}\times \mathbf x_{03}}{|\mathbf x_{02}\times \mathbf x_{03}|}.$$9 and $\Gamma$0 are computed with the same probing technique: the distance function is used to place edge-centered probe points along $\Gamma$1 and $\Gamma$2, and two-point finite differences provide directional derivatives while remaining on the reconstructed interface (Shin et al., 2017, Panda et al., 29 Sep 2025).

4. Boussinesq–Scriven surface rheology within LCRM

The 2025 formulation extends LCRM to interfacial flows with surfactant-induced surface viscous stresses. From an extrinsic viewpoint, the interface is treated as a two-dimensional manifold embedded in three-dimensional space. The surface identity tensor is

$\Gamma$3

and, for an inviscid interface, the surface-excess pressure tensor is

$\Gamma$4

Its surface divergence yields

$\Gamma$5

which recovers the normal capillary force and Marangoni tangential force (Panda et al., 29 Sep 2025).

The Boussinesq–Scriven constitutive law introduces surface shear and dilatational viscosities. The surface viscous stress is written as

$\Gamma$6

where the adopted form of the surface rate-of-deformation tensor is the extrinsic–Lagrangian Scriven form,

$\Gamma$7

The total surface lineal force due to surface viscosity becomes

$\Gamma$8

This formulation is implemented directly on the Lagrangian interface and then transferred to the Eulerian grid through the same discrete delta machinery used for capillary and Marangoni forces (Panda et al., 29 Sep 2025).

For the dilatational contribution, the scalar

$\Gamma$9

is defined as a surface-viscous tension. The corresponding force decomposes into a normal term and a tangential term,

$\Gamma$0

The normal component is cast in the same hybrid form used for capillary forces, while the tangential component uses probing to compute $\Gamma$1 edge-by-edge. The paper explicitly notes that the dilatational part is implemented using exactly the same machinery as capillary forces in LCRM, with $\Gamma$2 replaced by $\Gamma$3 (Panda et al., 29 Sep 2025).

For the shear contribution, the force is converted to an edge integral through the surface Gauss theorem,

$\Gamma$4

so accurate evaluation of $\Gamma$5 at edge midpoints is required. The implementation exploits the fact that the full Eulerian gradient $\Gamma$6 is already available on the grid. At an edge midpoint,

$\Gamma$7

with $\Gamma$8 obtained from the distance function. Then

$\Gamma$9

and the surface divergence required in $\sigma(\Gamma)$0 is obtained on the fly as

$\sigma(\Gamma)$1

This is one of the clearest examples of why the Lagrangian–Eulerian division of labor is useful: the interface remains sharp and explicitly represented, while the distance field supplies robust normals and the Eulerian grid supplies the full velocity gradient needed for surface rheology (Panda et al., 29 Sep 2025).

5. Numerical realization, force spreading, and parallel structure

The cited LCRM implementations use a staggered MAC grid and a projection method for incompressible Navier–Stokes. Convective terms are discretized with a second-order ENO scheme, and the pressure Poisson equation is solved with a parallel hybrid multigrid/GMRES solver that is reported as robust up to density ratios $\sigma(\Gamma)$2. Interface points are advanced by second-order Runge–Kutta, and interfacial quantities are updated on the moving triangulated surface (Shin et al., 2017, Panda et al., 29 Sep 2025).

Transfer between the Lagrangian surface and the Eulerian grid is handled through a compact-support discrete delta kernel. In the surface-viscous formulation, the kernel is written as

$\sigma(\Gamma)$3

and any scalar or vector quantity $\sigma(\Gamma)$4 defined on the interface is distributed by

$\sigma(\Gamma)$5

The soluble-surfactant formulation uses the same immersed-boundary style transfer to construct Eulerian fields $\sigma(\Gamma)$6 and $\sigma(\Gamma)$7, after which

$\sigma(\Gamma)$8

This Eulerian distribution preserves surfactant mass, facilitates coupling to the bulk concentration equation, and provides a field from which $\sigma(\Gamma)$9 can be evaluated and re-interpolated after interface reconstruction (Shin et al., 2017, Panda et al., 29 Sep 2025).

The probing technique is central to the numerical treatment of surface gradients. At an edge center, the normal from the distance function defines $\mu_s(\Gamma)$0 and $\mu_s(\Gamma)$1, and probe points are placed at a distance $\mu_s(\Gamma)$2 in appropriate directions while remaining on the interface through the $\mu_s(\Gamma)$3 constraint. Directional derivatives such as

$\mu_s(\Gamma)$4

are then used for diffusion, Marangoni stress, and, by direct analogy, gradients of surface-viscous tension. Because the method avoids connectivity-based stencils, these calculations remain local and are naturally adapted to distributed-memory parallelism (Shin et al., 2017, Panda et al., 29 Sep 2025).

The code base discussed in the cited papers is BLUE. The surface-viscous study reports that BLUE is massively parallel and tested up to 131,072 cores, while the soluble-surfactant study emphasizes that interface operations such as distance computation, delta distribution, and probing are local plus halo exchanges. This locality is presented as a major advantage over classical front-tracking schemes that rely on global mesh adjacency, which are described as cumbersome to parallelize (Shin et al., 2017, Panda et al., 29 Sep 2025).

6. Verification, applications, limitations, and broader significance

The published LCRM formulations are supported by a broad benchmark suite. For soluble surfactant transport, the uniformly expanding sphere test shows that numerical $\mu_s(\Gamma)$5 matches the analytical solution

$\mu_s(\Gamma)$6

extremely well on a coarse $\mu_s(\Gamma)$7 grid with 12 cells across the diameter and 1000 reconstruction steps, indicating very good total surfactant mass conservation and no degradation due to repeated interface reconstruction. For a stationary droplet with surface diffusion, one-dimensional profiles of $\mu_s(\Gamma)$8 match the exact solution closely even on $\mu_s(\Gamma)$9, and grid refinement from $\mu_d(\Gamma)$0 to $\mu_d(\Gamma)$1 yields between first- and second-order convergence in $\mu_d(\Gamma)$2 norm. Marangoni migration velocities converge toward the YGB prediction, and simple-shear drop deformation compares well with experiments and prior simulations; the reported rotation angle is 0.26 rad with surfactant versus 0.27 rad experimentally, and 0.65 rad without surfactant versus 0.61 rad experimentally. Total surfactant mass is conserved to within a few tenths of a percent on refined grids. The same formulation also handles adsorption, desorption, breakup into main and satellite drops, and large-scale annular film flow with insoluble and soluble surfactant (Shin et al., 2017).

For surface-viscous interfacial rheology, LCRM is validated against several canonical problems. In simple shear, mesh refinement with $\mu_d(\Gamma)$3 shows convergence of the major and minor axes $\mu_d(\Gamma)$4 and deformation $\mu_d(\Gamma)$5, with errors less than a few percent relative to Flumerfelt’s small-deformation theory as the grid is refined. Temporal and steady deformations compare closely with Pozrikidis’ BEM results, Flumerfelt’s analytical expression in the Stokes limit, and previous front-tracking simulations. Increasing surface shear viscosity slows deformation and reduces the steady deformation parameter, while large $\mu_d(\Gamma)$6 produces nearly rigid-like surface circulation and reduced tangential velocity gradients. For a rising surface-viscous drop in quiescent fluid, terminal velocity decreases with increasing $\mu_d(\Gamma)$7, and fine-grid errors are reported as $\mu_d(\Gamma)$8–$\mu_d(\Gamma)$9. In vertically oscillated two-layer systems, threshold acceleration $\varphi(\mathbf x)$0 for Faraday instability increases with the total Boussinesq number and matches two-dimensional ALE results within about $\varphi(\mathbf x)$1 across cases; separate effects of $\varphi(\mathbf x)$2 and $\varphi(\mathbf x)$3 are resolved, with slightly larger errors for shear viscosity. In Faraday-wave-induced atomization, the method reproduces deep craters, jet-like ligaments, and pinch-off events in the clean case, while added surface viscosity damps ligament formation, delays pinch-off, and changes droplet sizes and topology (Panda et al., 29 Sep 2025).

These results delimit the method’s current strengths and its acknowledged challenges. The cited comparisons state that LCRM offers excellent mass conservation, reported as less than $\varphi(\mathbf x)$4 error, improved curvature estimation, and reduced parasitic currents, all of which are especially important when interfacial forces depend sensitively on curvature and surface gradients. At the same time, the surface-viscous study notes that accuracy of surface shear viscosity appears slightly more sensitive than dilatational viscosity, with Faraday-wave thresholds deviating by up to about $\varphi(\mathbf x)$5, and that high resolution is required to reduce underestimation or overestimation in demanding DNS settings. The soluble-surfactant formulation also notes that explicit treatment of diffusion and surfactant can impose time-step constraints, and its bulk diffusion examples are limited to one phase, although extension to both phases is described as straightforward (Shin et al., 2017, Panda et al., 29 Sep 2025).

Within the cited research trajectory, LCRM functions as a general three-dimensional interface framework in which surface PDEs, Marangoni forcing, and full Boussinesq–Scriven interfacial rheology can be implemented in a single Navier–Stokes solver. This suggests a methodological niche between implicit-capture methods and connectivity-heavy front-tracking schemes: surface fields remain naturally defined on a Lagrangian manifold, while topology change, probing, and geometric reconstruction are delegated to a locally reconstructed Eulerian distance field. A plausible implication is that this division of labor is the central reason LCRM has been extended from capillary flow to soluble surfactant transport and then to surfactant-induced surface viscosity without changing its core interface philosophy (Shin et al., 2017, Panda et al., 29 Sep 2025).