# A mesoscale Eulerian numerical method for short-range repulsion in interfacial dynamics

Théo Rajim, Yedhir Mezache, Bastien Di Pierro, Pierre Trontin

► To cite this version:

Théo Rajim, Yedhir Mezache, Bastien Di Pierro, Pierre Trontin. A mesoscale Eulerian numerical method for short-range repulsion in interfacial dynamics. International Journal of Multiphase Flow, 2026, 194, pp.105495. 10.1016/j.ijmultiphaseflow.2025.105495. hal-05329939

HAL Id: hal-05329939

https://hal.science/hal-05329939v1

Submitted on 24 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A mesoscale Eulerian numerical method for short-range repulsion in interfacial dynamics.

Théo Rajim<sup>a</sup>, Yedhir Mezache<sup>a</sup>, Bastien Di Pierro<sup>a</sup>, Pierre Trontin<sup>a,\*</sup>

<sup>a</sup>Universite Claude Bernard Lyon 1, LMFA, UMR5509, CNRS, Ecole Centrale de Lyon, INSA Lyon, 69622 Villeurbanne France,

**Abstract**

The dynamics of bubble and droplet coalescence and repulsion are critical in a variety of natural and industrial multiphase systems. Their direct numerical simulation (DNS) remains challenging due to the multiscale nature of interfacial interactions. This study presents a mesoscale numerical approach that incorporates a repulsive disjoining pressure into an Eulerian framework to model near-contact interactions without explicitly resolving nanoscale forces. The method computes interfacial distances through a transport equation, enabling efficient and scalable implementation across parallel computing architectures. Unlike traditional marker-based or interface-tracking methods, the approach is independent of the number of bubbles or droplets in the domain. The method is validated on representative test cases, including droplet collisions and static foam networks, and is shown to accurately capture both dynamic inertial effects and quasi-static configurations. Notably, it enables DNS of wet foams with reduced liquid fractions compared to conventional level-set methods where repulsive interactions are not accounted for. This approach lays the foundation for simulating more realistic foam structures and analyzing destabilization phenomena such as drainage and coarsening, thus contributing to the predictive modeling of complex gas-liquid systems.

*Keywords:* short-range repulsive force, near-contact interactions, mesoscale modeling, complex flowing interfaces, foams.

<footnote>*Corresponding author (pierre.trontin@univ-lyon1.fr)</footnote>
Introduction

The dynamics of bubble coalescence and repulsion play a fundamental
role across various scientific and technological fields. Governed by complex
interfacial interactions, these phenomena influence a wide range of natural
and industrial processes. In chemical engineering, bubble coalescence plays
a critical role in gas-liquid reactors, governing the transition from the homo-
geneous regime, characterized by a monomodal bubble size distribution, to
heterogeneous regimes, where larger bubbles dominate [4, 28]. In environ-
mental sciences, it plays a key role in wastewater treatment as in the flotation
process [33], and gas exchange between the ocean and atmosphere [11]. In
medicine, these dynamics are leveraged in applications such as acoustically
driven microbubbles and their interactions with soft and biological matter
[12].

In multiphase flows containing bubbles, a foam network is defined by
a low liquid fraction, where the bubbles are separated solely by thin liquid
films of a soapy solution. Their insulating ability, excellent coverage capacity,
lightness, and strong stabilizing properties make them a natural solution in
many natural and industrial processes [5, 16] with applications in building
design, environment, aerospace, cosmetics and food industries [38].

In two-phase flows containing bubbles, direct numerical simulations (DNS), facilitated by continuously increasing computing capabilities, allow for a detailed characterization of bubble interactions and the underlying physical mechanisms. Methods for interface reconstruction between two phases include diffuse interface approaches, such as the phase-field method [50, 52], and sharp interface methods, which utilize either Lagrangian or Eulerian formulations. Lagrangian or front-tracking methods [48, 47, 36] rely on a marker-based description of the interface, which can make interfacial topological changes such as coalescence and break-up cumbersome. Eulerian or front-capturing, methods describe the interface implicitly, either through the zero isosurface of a signed distance function in the level-set method [29, 41], via the volume fraction of the liquid phase in the volume-of-fluid (VOF) method [49, 32, 34], or by adding an extra force for the lattice Boltzmann method (LBM) to represent the interaction between the two phases [35].

Interactions between bubbles, whether attractive or repulsive, arise from
the competition between nanoscale near-contact forces, including electro-
static, van der Waals and steric forces [5, 2]. Therefore, DNS of bubbly flows
is inherently multiscale, spanning a wide range of scales—from the largest
bulk flow structures (on the order of 100 µm for bubble diameter) to the
smallest scales associated with the thin liquid film between bubbles, which
becomes increasingly thinner (on the order of a few nanometers) during po-
tential bubble collisions. Furthermore, the drainage time required for the
liquid film before bubble coalescence requires smaller time steps and longer
computation times. Thus, with the exception of specific theoretical configu-
rations, such as in [23], where the interaction between two planar interfaces
is calculated using a moving mesh, capturing all scales by DNS is out of reach
for realistic 3D or even 2D bubbly flows [18] with a large number of bubbles
likely to interact by collision, even when local mesh refinement techniques are
used. In configurations where bubble collisions are rare, such as at low void
fractions, or when the focus is on global statistics, such as turbulence induced
in bubble-laden flows rather than individual bubble-bubble interactions, the
aforementioned methods have proven to be highly effective [13, 17, 8]. When
bubble interactions become significant, as in dense bubbly flows, or even
drive the overall equilibrium of a two-phase mixture, as observed in foams,
the impact dynamics governing coalescence or repulsion must be analyzed
for each bubble pair. In the framework of DNS, only a mesoscale representa-
tion of near-contact forces enables the description of bubble interactions at
a scale compatible with current computational capabilities. Existing meth-
ods, as described in [10] and based on the CLSVOF method [40], assign each
bubble its own level-set function and VOF color function. As bubble coales-
cence is not allowed by the method, issues arise where insufficient numerical
resolution prevents accurate representation of fluid motion in the thin film
between the interfaces. The method is therefore only relevant for bubbly
flows with low gas volume fractions. A modification is proposed in [21] to
take collisions into account, but the numerical cost of the method increases
with the number of bubbles and the resulting number of level-set and VOF
fields to be managed.

In addition to bubble interactions, short-range repulsive forces in interfa-
cial dynamics pose significant challenges in the simulation of binary droplet
collisions [37]. To address these effects, a hybrid method was proposed in
[22], coupling the lubrication equation for the gas film with the Navier-Stokes
equations for the droplets, while incorporating nanoscale van der Waals forces
to prevent bouncing. In [24], a VOF-based numerical framework was devel-
oped for droplet tracking, enabling the simulation of both bouncing and
coalescence events. However, these studies impose symmetry boundary con-
ditions, allowing the simulation of only a single droplet and thereby sim-
plifying the binary collision scenario. Despite their accuracy, such methods
demand extremely fine computational grids, and extending them to systems
involving large droplet populations remains computationally intensive and
technically complex.

Among the proposed mesoscale formulations, a common approach is to
incorporate a disjoining pressure term of the order of the capillary pressure
within the surface tension force [51, 27]. This suggests a certain universal-
ity of the physical process, allowing the interaction between two interfaces
to be described globally by modeling its effect on large-scale interface dy-
namics, without the need to consider the precise small-scale nature of the
forces involved. The main challenge here is calculating the distance between
the interfaces for two surface elements facing each other. In [51], it is cal-
culated using a ray-tracing method and is defined as the minimum distance
between pairs of opposing donor/acceptor markers describing facing inter-
faces. Similar to the techniques used in Lagrangian front-tracking methods,
this approach is accurate but computationally expensive due to the necessity
of describing the interface with markers. For parallelized HPC computations,
the fact that an donor/acceptor pair of markers is not processed by the same
processor requires additional inter-processor communications, which can lead
to a significant increase in computation time, especially for a high number
of bubbles in the domain. For non-planar facing interfaces, the very notion
of facing interfaces is not defined, and the concept of reference surface is
introduced in [20] to help define them. The concept of a repulsive force has
previously been introduced within the framework of the LBM method to sim-
ulate the dispersion of one phase into another [1, 30], thereby enabling the
simulation of phenomena such as coarsening in foam-like systems. By appro-
priately tuning the phase-separating interactions, it is possible to obtain a
positive disjoining pressure.

The objective of this paper is to characterize the macroscopic dynamics
governing the interaction between two opposing interfaces, without explicitly
resolving near-contact phenomena such as nanoscale short-range forces. We
propose an approach based on a mesoscale formulation, incorporating a re-
pulsive force between the interfaces based on a disjoining pressure. The inter-
facial distance is calculated in an Eulerian framework by solving a straightfor-
ward transport equation, similar to the one used in the level-set redistancing
step [39]. The method therefore does not depend on the number of bubbles
or droplets in the domain. Its 3D formulation is straightforward, and the
numerical schemes used are the same as those for transporting the level-set
function or for its redistancing. Parallelization is efficient due to the Eulerian nature of the method. Although the ultimate objective is the direct numerical simulation (DNS) of a foam network, the method is sufficiently general to capture interactions between interfaces, such as the collision and rebound of two droplets. This will enable the simulation of drier, and therefore more realistic foams, overcoming a major limitation of current DNS approaches that do not consider repulsive forces and thus cannot simulate low liquid fractions without a high risk of bubble coalescence (see [26] where simulated liquid fractions are greater than 0.3). Unlike LBM-type methods, the present approach employs a sharp-interface formulation that enables the resolution of thin-film drainage without being constrained by interface thickness. This formulation also accommodates large liquid-gas density ratios and improves curvature evaluation through the level-set method, thereby reducing spurious currents. However, in contrast to LBM, which can readily handle realistic geometries and complex meshes, DNS with sharp interfaces are generally limited to academic configurations with simplified grids.

In the first section, the numerical method for reconstructing the distance between two facing interfaces is presented. The disjoining pressure model is also formulated as a function of the reconstructed interfacial distance. Section 2 presents a numerical validation of the proposed approach through four representative test cases, encompassing both droplet collision scenarios and the dynamics of a foam bubble network. Finally, conclusions are drawn in Sec. 3.

# 1. Numerical methods

## 1.1. *Fluid flow solver*

The simulation of liquid-gas flows at moderate velocities is considered assuming an isothermal behaviour. The incompressible Navier-Stokes equations for the one-fluid formulation of immiscible two-phase flows are solved. A level-set method [29, 41] is used to capture the interface between the two phases. It is implicitly given by the zero of a smooth function $\phi(\mathbf{x}, t)$ defined as the signed distance function to the interface (Fig. 1). The level-set function $\phi$ is chosen positive in the liquid and negative in the gas. The interface normal vector $\mathbf{n}$ points towards the positive values of $\phi$. In this study, the level-set method is applied to particle interactions and collisions. More broadly, it represents a general framework for two-phase flows with interfacial dynamics. Surface forces at the interface like capillary and repulsive forces
are taken into account by the CSF (Continuum Surface Force) method [3].
The equation for the velocity field **u** is then

<figure><img src="image_2.png" /><figcaption>Figure 1: The interface is implicitly given by the zero of the level-set function $\phi$. For any point $M(\mathbf{x}, t)$ in the domain, $\phi(\mathbf{x}, t)$ represents the signed distance from $M$ to the interface.</figcaption></figure>

$$ \rho(\phi) \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) = - \nabla p + \nabla \cdot \left[ \mu(\phi) \left( \nabla \mathbf{u} + \nabla \mathbf{u}^T \right) \right] \\ - \left[ \sigma\kappa(\phi) + \Pi_d(\phi) \right] \delta(\phi) \mathbf{n} \quad (1) $$

with the incompressibility condition

$$ \nabla \cdot \mathbf{u} = 0 \quad (2) $$

where $\rho$, $\mu$, $p$, $\sigma$, $\kappa$ and $\delta$ are the density, dynamic viscosity, pressure, surface tension coefficient, interface curvature and Dirac function indicating the interface. The chosen form of $\delta(\phi)$ follows the formulation proposed by [31]:

$$ \delta(\phi) = \begin{cases} \frac{1}{2\epsilon_{\delta}} \left[ 1 + \cos\left(\frac{\pi\phi}{\epsilon_{\delta}}\right) \right] & \text{if } |\phi| \le \epsilon_{\delta} \\ 0 & \text{otherwise} \end{cases} \quad (3) $$

where $\epsilon_{\delta} = 1.5\Delta x$, with $\Delta x$ denoting the uniform grid spacing in all three directions $\mathbf{x}$, $\mathbf{y}$ and $\mathbf{z}$ of the Cartesian grid. Thus, in the context of the CSF method, the capillary and repulsive forces (Eq. (1)) are regularized in the vicinity of the interface over a distance corresponding to 1.5 grid cells on either side of the front. $\Pi_d$ is the disjoining pressure and represents repulsive forces. It will be detailed below. Curvature $\kappa$ is defined by $\kappa = \nabla \cdot \mathbf{n}$.
The evolution of the interface is obtained by solving the following transport equation for $\phi$:

$$ \frac{\partial \phi}{\partial t} + (\mathbf{u} \cdot \nabla)\phi = 0 \qquad (4) $$

Equation (4) ensures that the level-set $\phi = 0$ propagates at the correct velocity. However, as $\phi$ is advected, it loses its intrinsic property of being a distance function ($||\nabla\phi|| \neq 1$). This can lead to significant mass loss or gain. To mitigate this issue, a level-set redistancing step is required. This is accomplished by solving a Hamilton–Jacobi equation for $\phi$ until a steady-state solution is reached [41].

Incompressibility constraint (Eq. (2)) is ensured by classical projection methods [7, 43]. The scalar variables are defined cell-centered, while a staggered MAC [15] uniform Cartesian grid is used for the velocity components. Convective terms in Eqs. (1) and (4) are discretized with a 5th order accurate WENO scheme algorithm [25] while a 2nd order accurate finite difference method is used for the diffusive terms in Eq. (1). Time integration is performed with a 3rd-order accurate TVD Runge–Kutta scheme. The Poisson equation for the pressure is solved by a fast multigrid preconditioned conjugate gradient method [42, 46], the operator being symmetric and negative-definite. All these numerical methods are integrated into a software [9, 53, 54] that has demonstrated its ability to carry out direct numerical simulations in two-phase flow conditions with a highly deformed interface [44, 45].

## 1.2. Mesoscale modeling of repulsive forces between interfaces.

To mimic short-range repulsive effects at the mesoscale, an exponential law is used for the disjoining pressure $\Pi_d$:

$$ \Pi_d = P_0 \exp(-h(\phi)/h^*) \qquad (5) $$

where $P_0$ is the amplitude of the repulsive pressure, $h$ is the distance between two facing interfaces and $h^*$ is a characteristic length scale. It is important to note that the numerical method developed in this study is not constrained to an exponential law but can be applied to any continuously differentiable and monotonically decreasing function for $\Pi_d$. However, the formulation of Eq. (5) avoids singular behaviour, especially for $h \to 0$. Calibration for $P_0$ and $h^*$ is proposed here in the context of foams, i.e. gaseous inclusions within a liquid film. It remains valid for the case of liquid droplets dispersed in a gaseous medium.
Cantat et al. [5] demonstrated that the equality between the Laplace pressure jump and the disjoining pressure is a necessary condition for both the stability and the existence of a wet foam. According to their analysis, $P_0$ is chosen such that $P_0 \sim \sigma/R_{PB}$, where $R_{PB}$ denotes the curvature radius of the Plateau border, thereby ensuring that the capillary and disjoining pressures are of the same order of magnitude.

The parameter $h^*$ characterizes the effective range of the repulsive forces. If $h_c$ is defined as the cutoff scale beyond which $\Pi_d(h) \approx 0$ for $h > h_c$, then, based on the definition given in Eq. (5), a value of $h_c = 5h^*$ is proposed. A calibration for $h_c$ will be presented in Sec. 2.

The distance $h$ between two interfaces (Eq. (5)) can be defined using the level-set function $\phi$ only. Consider the case shown in Fig. 2a, where two bubbles are separated by a liquid film. A 2D configuration has been chosen for the explanation, but the 3D extension is straightforward. We are interested in defining the distance $h$ between the two interfaces facing each other at points A and B. Figure 2b shows the evolution of the signed distance

<figure><img src="image_3.png" /><figcaption>Figure 2: Ridge line and definition of the distance between two facing interfaces.</figcaption></figure>

$\phi$ to the interface as a function of $s$, the curvilinear abscissa along segment $[AB]$. The two bubbles correspond to $s < s_A$ and $s > s_B$ ($\phi < 0$) whereas the liquid film between the bubbles corresponds to $s_A < s < s_B$ ($\phi > 0$). Points A and B are on the interface so that $\phi(s_A) = \phi(s_B) = 0$. Between
A and B, $\phi$ is maximum for $s = s_{RL}$ at the ridge line<fnref n="1" />, and is equal to $h/2$, where $h$ is the distance between the two facing interfaces (Fig. 2a). This means that the distance between the interfaces facing each other at A and B is simply twice the value of the level-set function evaluated at the associated ridge line. Since the level-set function $\phi$ is known at all times and throughout the entire computational domain by solving Eq. (4), all that remains now is to define the ridge lines in the domain and evaluate $\phi$ along these ridge lines. Once the distance between two facing interfaces is known at the ridge line, it must be propagated to the liquid/gas interfaces where the repulsive force is applied according to the CSF formulation of the disjoining pressure. This is described in the next section. In the rest of the paper, we will refer to ridge surfaces as an immediate 3D generalization of 2D ridge lines.

### 1.3. Ridge surface detection and distance propagation between facing interfaces.

The complete process, from ridge surface detection to the propagation of the distance between two facing interfaces, is described in Fig. 3 for the 2D case of three interacting bubbles. Figure 3a shows the level-set field around the three bubbles in the liquid phase.

The first step is the detection of the ridge surfaces within the flow. Due to its distance characteristics, the $||\nabla\phi|| = 1$ property for the level-set function is verified throughout the entire computational domain for the level-set function, except at ridge surfaces, lines or points which are equidistant from two interfaces. The propagation of distance between two facing interfaces for disjoining pressure calculation occurs exclusively within the liquid phase. This can be effectively disregarded by restricting ridge surface detection to the liquid phase ($\phi > 0$). The following simple criterion is used:

$$ | ||\nabla\phi|| - 1 | > \varepsilon_{ridge} \qquad (6) $$

where $\nabla\phi$ is estimated with a second order finite difference scheme. A posteriori numerical tests show that $\varepsilon_{ridge} = 10^{-2}$ is sufficient to reliably detect ridge surfaces. The detection step is illustrated in Fig. 3b, where the ridge lines are shown. Figure 4 illustrates the case of two facing particles of equal diameter, with the ridge region (gray band) captured by the numerical method (Eq. (6)). The iso-contours of the level-set function $\phi$ are displayed,

<footnote n="1">¹The rigde line is defined as the set of points within the liquid phase which are equidistant to two interfaces</footnote>
<figure><img src="image_4.png" /><figcaption>(a) Level-set function field around the three bubbles in the liquid phase.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Ridge line detection and representation of $\psi_0$ along them.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) Propagation of $\psi$ at $\tau = 0.044$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d) Propagation of $\psi$ at $\tau = 0.088$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e) Propagation of $\psi$ at $\tau = 0.177$ (steady state).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(f) Projection of stationary field $\psi$ at $\tau = \tau_{max}$ on the liquid/gas interfaces.</figcaption></figure>

Figure 3: Ridge surface detection algorithm and distance propagation for the 2D case of three interacting bubbles. The dashed lines represent the liquid/gas interfaces.
<figure><img src="image_10.png" /><figcaption>Figure 4: Comparison between the numerically captured ridge region (gray band) and the theoretical ridge line (dashed white line). Iso-contours of the level-set function $\phi$ are shown, with the interfaces of the two facing particles ($\phi = 0$) highlighted in red. The computational mesh is also displayed.</figcaption></figure>
where the particle interfaces (corresponding to $\phi = 0$) are highlighted in red. The theoretical ridge line is indicated by a dashed white line, and the computational mesh is also displayed. Consistent with the discretization scheme applied to Eq. (6), the numerically computed ridge region is localized within a band extending $\pm \Delta x$ around the theoretical ridge line.

The next step is to transport the value of the level-set function estimated at the ridge surfaces (denoted $\phi_{RL}$) to the liquid/gas interfaces. To do this, the following transport equation is solved for the function $\psi(\tau, \mathbf{x})$:

$$ \left\{ \begin{array}{l} \frac{\partial \psi}{\partial \tau} + \mathbf{V} \cdot \nabla \psi = 0 \\ \psi(0, \mathbf{x}) = \psi_0(\mathbf{x}) \end{array} \right. \qquad (7) $$

where $\tau$ is a pseudo-time and the initial condition $\psi_0(\mathbf{x})$ is given by:

$$ \psi_0(\mathbf{x}) = \begin{cases} \phi_{RL} & \text{if } \mathbf{x} \in \text{ridge surface} \\ \psi_0^{max} & \text{otherwise} \end{cases} \qquad (8) $$

$\psi_0^{max}$ denotes an upper bound for $\psi_0$, the magnitude of which is discussed in detail below. Initial condition $\psi_0$ is shown at the ridge lines in Fig. 3b. In Eq. (7), the information is transported downhill with:

$$ \mathbf{V} = \begin{cases} -\nabla\phi, & \text{if } \phi \ge -\epsilon_{\delta} \text{ (liquid phase and interfacial region)} \\ 0, & \text{otherwise (gas phase)} \end{cases} \qquad (9) $$

In this way, $\psi$ propagates along the vector normal to the isosurfaces of the level-set function $\phi$, up to a distance $\epsilon_{\delta}$ inside the bubble, in agreement with the CSF formulation which regularizes the disjoining pressure over a region of thickness $2\epsilon_{\delta}$ surrounding the interface (see Eq. (3)). Section 1.4 discusses the significance of transport direction in defining the distance between two facing interfaces. Figures 3c - 3d illustrate two increasing values of $\tau$ during the propagation step, wherein information progressively propagates from the ridge surfaces to the liquid-gas interfaces. Stationary state is reached for $\tau = \tau_{max}$ (Fig. 3e). Figure 3f shows the values of the stationary $\psi$ field at the centers of the liquid-side cells situated within a distance of 1.5 grid cells from the interface, where the distance information between facing interfaces is required for the CSF formulation of the repulsive forces. Thus, at each interface located at position $\mathbf{x}$, the distance to the facing interface is equal to $2\psi(\tau_{max}, \mathbf{x})$ so that the disjoining pressure in Eq. (5) can be written:

$$ \Pi_d(\mathbf{x}) = P_0 \exp(-2\psi(\tau_{max}, \mathbf{x})/h^*) \qquad (10) $$
In practice, $\tau_{max}$ corresponds to the time needed to transport $\phi_{RL}$ over a distance of at least $5h^*$ at speed $V$, $\Pi_d$ being close to zero for these values of $h$ (Eq. (5)). In Eq. (8), $\psi_0^{max}$ is chosen to be large enough ($\psi_0^{max}/h^* \gtrsim 10$) to ensure that in regions not concerned by the transport of $\phi_{RL}$, the corresponding value of the disjoining pressure $\Pi_d$ is almost zero (Eq. (10)). To avoid diffusion errors during the transport process, Eq. (7) is numerically solved with a global second-order accurate Lax-Wendroff/Runge-Kutta 2 scheme.

The computational cost of propagating interfacial distance information from ridge surfaces to water/gas interfaces is comparable to that of the level-set redistancing step [41]. Importantly, it remains independent of the number of interfacial inclusions (bubbles or droplets) present within the computational domain.

### 1.4. Impact of the definition of interfacial distance

To our knowledge, there is no consensus on the definition of the distance $h$ used in Eq. (5) for the repulsive force. For instance, in [51], it is defined as the closest distance between a pair of donor/acceptor points, each located on two facing interfaces. In [27], it is simply defined as the distance along the normal vector $\mathbf{n}$ to the interface, with no further details on how to calculate it. For planar interfaces, all definitions of the distance between two facing interfaces are equivalent. Differences may arise when the radii of curvature of facing interfaces are comparable to, or smaller than the distance between the interfaces. In this context, we propose to compare different methods of calculating the distances for the disjoining pressure in the case of two spherical particles of diameter $D$ facing each other, and being separated by a distance $h_0$. It is assumed that only the two opposing hemispheres of the facing particles contribute to the force calculation. A description of the 3D half-problem is shown in Fig. 5, the second particle being obtained by symmetry of the first with respect to the shaded ridge plane. The $\mathbf{z}$-axis, which connects the centers of the two opposing particles, is henceforth referred to as the inter-particle axis. Consider a point $M$ on the left sphere identified by $\theta$ and $\phi$, respectively the polar and azimuthal angles. The distance from $M$ to the right sphere is to be calculated. Due to rotational invariance about the $\mathbf{z}$-axis, it depends only on $h_0$ and $\theta$ and is denoted $h(h_0, \theta)$.

We consider the definitions of two distinct distances. The first one, denoted as $h_{ridge}$, is the distance described in this manuscript. Through straightforward geometrical considerations (Fig. 5), the expression for $h_{ridge}$
<figure><img src="image_11.png" /><figcaption>Figure 5: Two particles facing each other. Only one particle is represented, the other being symmetrical to the first with respect to the ridge plane (shaded plane). Distance $h_{ridge}$ is also shown.</figcaption></figure>

is derived as follows:

$$ \frac{1 + h_0/D}{1 + h_{ridge}(h_0, \theta)/D} = \cos(\theta) \qquad (11) $$

For each value of $h_0$, one can define a corresponding angle $\theta_{max}$ such that if $|\theta| > \theta_{max}$ then $h_{ridge} > h_c$ and $\Pi_d(h_{ridge}) = 0$. The expression for $\theta_{max}$ can be written as:

$$ \cos(\theta_{max}) = \frac{1 + h_0/D}{1 + h_c/D} \qquad (12) $$

The parameter $\alpha$ is defined as $\alpha = \theta/\theta_{max}$, with $\alpha \in [0; 1]$ specifying the interval over which $\Pi_d(h_{ridge})$ is nonzero.

The second distance, denoted $h_{min}$, represents the closest distance between $M$ and the right sphere. We now show that, although the distances $h_{ridge}$ and $h_{min}$ differ, the resulting discrepancies in the calculated disjoining pressures and repulsive forces remain limited. Let $F_{rep}$ denote the z-component of the repulsive force exerted by the left sphere on the right one. This force depends solely on $h_0$ and can be expressed as $F_{rep}(h_0) =$
$\int_{0}^{\pi/2} f_{rep}(h_0, \theta) d\theta$ where the force density $f$ is given by:

$$f_{rep}(h_0, \theta) = 2\pi(D/2)^2 P_0 \exp \left[ -\frac{h(h_0, \theta)}{h^*} \right] \sin(\theta) \cos(\theta) \quad (13)$$

The elementary force $f_{rep}d\theta$ corresponds to the force applied by the corona at polar angle $\theta$ with thickness $(D/2)d\theta$ (illustrated in light blue in Fig. 5). Figure 6a shows the variation of the dimensionless force density $f_{rep}/(P_0\pi D^2/4)$ as a function of $\alpha$ for different values of $h_0/h_c$, and for the two formulations $h_{ridge}$ and $h_{min}$. Here, $\alpha$ is defined from Eq. (12) where $h_0 = 0$. For any given

<figure><img src="image_12.png" /><figcaption>Figure 6: Impact of the definition of inter-particle distance on the repulsive force $F_{rep}$.</figcaption></figure>

inter-particle distance $h_0$, the dominant contribution to the repulsive force originates from the corona regions located near $\alpha \approx 0.3$. For coronas closer to the inter-particle axis ($\alpha < 0.3$), the disjoining pressure is higher. However, the limited surface area over which the force is applied results in a reduced net contribution to the overall repulsive interaction. In this region, the choice of the distance definition between $h_{ridge}$ and $h_{min}$ has minimal influence on the force magnitude. For coronas farther from the axis ($\alpha > 0.3$), the area of force application increases, but the disjoining pressure decays exponentially, thereby limiting the overall contribution to the total repulsive force. In this region, noticeable differences arise between $h_{ridge}$ and $h_{min}$ for the computation of $F_{rep}$. In accordance with the definition of $\Pi_d$, the repulsive force is a decreasing function of the inter-particle distance $h_0$.

Regarding the global calculation of $F_{rep}$ over the entire left hemisphere, Fig. 6b shows the variation of $F_{rep}/(P_0\pi D^2/4)$ as a function of the inter-particle distance $h_0$. For $h_0 = 0$, corresponding to the limiting case where
the particles are in contact, the relative error between the two values of $F_{rep}$ calculated from $h_{ridge}$ and $h_{min}$ is 12%. This gap decreases as the inter-particle distance increases, with a concurrent reduction in the magnitude of the repulsive force. Thus, the differences observed between $h_{ridge}$ and $h_{min}$ in the corona regions farthest from the inter-particle axis have only a limited influence on the overall calculation of the repulsive force $F_{rep}$.

We now examine the case of two facing particles that are not necessarily spherical and may undergo deformation due to repulsive forces. To facilitate an analytical formulation of $h_{ridge}$, we limit our analysis to a pair of facing non-deformable oblate spheroids. Each spheroid is generated by rotating an ellipse $\mathcal{E}$ about its minor axis, which is aligned with the z-axis. The semi-major axis $a$ of the ellipse $\mathcal{E}$ is defined such that $2a = D$, where $D$ corresponds to the diameter of the spheres previously considered. The semi-minor axis $b$ of the ellipse $\mathcal{E}$ is determined from its semi-major axis $a$ and eccentricity $e$ using the relation $b = a\sqrt{1 - e^2}$. Figure 7a illustrates the configuration for the case $e = 0.9$. The expressions for $h_{ridge}$ and $f_{rep}$ in the case of two facing

<figure><img src="image_13.png" /><figcaption>Figure 7: Impact of the definition of inter-particle distance on the repulsive force. Case of two facing oblate spheroids.</figcaption></figure>

oblate spheroids are given by:

$$ \frac{1 + \frac{h_0}{2b}}{1 + \frac{h_{ridge}(h_0, \theta)}{2b} \frac{1}{\sqrt{1 - e^2 \sin^2(\theta)}}} = \cos(\theta) \quad (14) $$
and

$$f_{rep}(h_0, \theta) = 2\pi a P_0 \exp \left[ -\frac{h(h_0, \theta)}{h^*} \right] \sin(\theta) n_z(\theta) \sqrt{a^2 \cos^2(\theta) + b^2 \sin^2(\theta)} \quad (15)$$

where parameters $\theta$, $h_0$, and $h_{ridge}$ are derived from Fig. 5, with the z-axis corresponding to the minor axis of the spheroid (Fig. 7a). The component $n_z$ denotes the projection of the outward normal vector to the spheroid at point M onto the z-axis. Equations (14) and (15) reduce appropriately to Eqs. (11) and (13) in the special case of spheres, where $a = b = D/2$ and $e = 0$. Figure 7b shows the variation of the dimensionless force density $f_{rep}/(P_0\pi a^2)$ as a function of $\alpha$ for various spheroidal shapes, or equivalently, for different values of the eccentricity $e$. The two formulations $h_{ridge}$ and $h_{min}$ are studied. Only the ratio $h_0/h_c = 0.1$ is analyzed. As in the spherical case (Fig. 6a), conclusions for higher ratios can be inferred analogously. To enable meaningful comparisons across different geometries, the same value of $\theta_{max}$ as in the spherical case is used when computing $\alpha$ in Fig. 7b. The black curves correspond to the spherical reference case from Fig. 6a, and are included here for comparison. For spheroidal particles, an increase in eccentricity, corresponding to greater flattening of the surfaces facing each other, leads to a divergence between the repulsive forces $f_{rep}$ predicted by the $h_{ridge}$ and $h_{min}$ formulations. This divergence occurs at higher values of $\alpha$ ($\alpha = 0.3$ for $e = 0$ and $\alpha = 0.5$ for $e = 0.9$). This indicates that, during particle collision, the influence of the distance definition for computing disjoining pressure is most pronounced at the initial stage of contact, when the particles are still undeformed. As the collision progresses and the particles approach each other and flatten, the effect of the distance definition gradually diminishes.

When simulating a foam network, the choice of a distance for calculating the repulsive force arises specifically for non-planar interfaces where the local radius of curvature is comparable to the interfacial distance. However, such singularities represent only a minor portion of the liquid-gas interfaces observed near the Plateau borders (see Sec. 2.3). Thus, in bubble networks, the errors arising from the definition of distance are of a lower order compared to those associated with the definition of the constants in the disjoining pressure model, which can exhibit an order-of-magnitude discrepancy, as reported in [51].
## 2. Numerical results

We propose to investigate four test cases. The first aims to validate the integration of disjoining pressure and repulsive forces into the Navier-Stokes equations in an idealized configuration of a head-on collision between two rigid spherical particles. The second, more physically realistic case involves head-on and off-axis bouncing collisions between two particles. The results demonstrate that the implemented repulsive forces successfully capture the bouncing dynamics and effectively prevent coalescence between the particles. The third one is designed to validate the method on a monodispersed, ordered bubble network, with a focus on its capability to accurately reproduce characteristic polyhedral structures. The fourth case extends to a polydisperse foam network.

### 2.1. Head-on collision between two rigid spherical particles

The case of a head-on collision between two facing rigid particles of density $\rho_p = 997 \text{ kg.m}^{-3}$ is studied (Fig. 8). In a parallelepiped box of dimension $(L_x \times L_y \times L_z) = (20 \text{ mm} \times 20 \text{ mm} \times 40 \text{ mm})$, two particles of diameter $D = 10 \text{ mm}$ initially separated by a distance $h_0/D = 1.6$ are injected head-on in the z-direction at velocities $\mathbf{u}_0$ and $-\mathbf{u}_0$ respectively. The density and

<figure><img src="image_14.png" /><figcaption>Figure 8: Head-on collision between two particles: initial conditions.</figcaption></figure>

dynamic viscosity of the fluid surrounding the particles are $\rho_f = 1 \text{ kg.m}^{-3}$ and $\mu_f = 1.85 \text{ mPa.s}$. Initial velocity magnitude is $u_0 = 6 \text{ mm.s}^{-1}$. The
corresponding particle Reynolds number $Re_p = \rho_f u_0 D / \mu_f$ is approximately 0.03. The surface tension coefficient $\sigma$ is set to 750 mN.m$^{-1}$, which induces a Weber number $We = \rho_f u_0^2 D / \sigma$ well below 1 (here $We \approx 5 \times 10^{-7}$), thus preventing deformation of the particles during movement, particularly when subjected to the action of repulsive forces. To avoid resolving the boundary layer at the walls, a slip-type boundary condition is applied to each of the six faces of the parallelepiped. For this analytical test case, the cut-off scale $h_c$ is chosen to be of the same order of magnitude as the initial inter-particle distance $h_0$, so that the two particles experience the repulsive force throughout the entire interaction, from the stage when they approach each other to the stage when they move apart. Thus, $h_c/h_0 = 1.125$ is chosen. To minimize travel distances and reduce the size of the computational domain, the disjoining pressure magnitude is calibrated so that the two particles travel approximately 10% of their initial separation toward each other. This is subsequently verified for $P_0 = 50$ Pa. The maximum simulation time is $t_{max} = 0.6$ s.

We first assess the numerical accuracy of transporting the inter-particle distance from the ridge surfaces to the liquid/gas interfaces of the two facing particles. To this end, the solution in the plane normal to the y-axis through point O (Fig. 5) at the initial condition ($t = 0$ s) is selected. The theoretical distance from any point M on the left sphere to the right sphere can be determined analytically using Eq. (11). It is denoted as $h_{ridge}^{theo}$ and is used as a reference. The distance obtained by DNS through the solution of a transport equation (Eq. (7)), as described in Sec. 1.3, is denoted $h_{ridge}^{num}$. Figure 9 presents the relative error $err_h$ between $h_{ridge}^{num}$ and $h_{ridge}^{theo}$ as a function of $\phi$, ranging from $\phi = h_{ridge}^{theo}/2$ at the ridge surface to $\phi = 0$ at the particle interface(see Fig. 8). Results are presented for different values of $\alpha$ (Eq. (12)), with $h_0/D = 1.6$ corresponding to the initial condition, and for several spatial resolutions, ranging from $D/\Delta x = 16$ to $D/\Delta x = 64$. The configuration with $\alpha = 0.9$ is illustrated in Fig. 8. For coronas near the inter-particle axis ($\alpha \to 0$), the numerical error in transporting the inter-particle distance to the liquid/gas interface is below 1%. Even for the coronas farthest from the axis ($\alpha \to 1$), the error does not exceed 3%, even with the coarsest grid $D/\Delta x = 16$. Thus, the second-order accurate Lax-Wendroff/Runge-Kutta 2 numerical scheme used to transport the inter-particle distance from the ridge surface to the liquid/gas interface introduces an error of only a few percent relative to the theoretical solution. This error is negligible compared to the
<figure><img src="image_15.png" /><figcaption>Figure 9: Relative error (in %) of $h_{ridge}^{num}$ with respect to the theoretical solution $h_{ridge}^{theo}$.</figcaption></figure>

uncertainties in estimating $P_0$ (Eq. (5)), which can sometimes vary by an order of magnitude, as observed in [51]. In addition to assessing the accuracy of the numerical scheme used for the transport of $\psi$, this test case also serves to validate the quality of the initial condition $\psi_0$ (Eq. (7)), obtained from the ridge detection algorithm, as well as the accuracy of the reconstructed transport velocity (Eq. (9)). The latter may be less accurate in the vicinity of the ridge, where the level-set function $\phi$ is no longer continuously differentiable.

Next, we analyze the impact dynamics. Figure 10 illustrates the time evolution of the DNS results for the head-on collision of two particles. The magnitude of the repulsive force is also presented. In the initial step ($t/t_{max} < 0.5$), as the particles move closer together, the repulsive force is triggered, causing the particles to move apart ($t/t_{max} > 0.5$). The iso-value lines for the level-set function $\phi$ and those for the distance $h = 2\psi$ between facing interfaces are orthogonal to each other, confirming the transport of $\psi$ along the gradient of the level-set function $\phi$ (Eqs (7) and (9)). Figures 10a-10c reveal artifacts aligned with the diagonal directions of the grid. This phenomenon arises because the numerical solution of the transport equation for $\psi$ is not based on a fully multidimensional scheme. Instead, it employs the more conventional approach of sequentially applying a one-dimensional algorithm
<figure><img src="image_16.png" /><figcaption>Figure 10: Time evolution of head-on collision. The magnitude of the repulsive force is shown. Black lines: iso-values of the level-set function $\phi$. White lines: iso-values of the distance $h = 2\psi$ between facing interfaces. Black and white lines are orthogonal. Iso-lines are drawn in the median y-normal vertical plane.</figcaption></figure>
along the **x**, **y**, and **z** directions.

To enable a more detailed analysis of the interaction between the two
particles, we first develop a quasi-analytical solution to predict their trajec-
tories. This solution, referred to as *Analytical* and presented in Appendix A,
describes the motion of two rigid spherical particles driven solely by repulsive
and drag forces. Figure 11a shows the time evolution of the inter-particle
distance $h_0$ as computed by DNS. The solid black curve represents the *Ana-

<figure><img src="image_17.png" /><figcaption>Figure 11: Time evolution of the positions and velocities of two rebounding particles. Dashed lines correspond to the case where $F_{rep} = 0$.</figcaption></figure>

lytical solution (see Appendix A), while the dashed black curve illustrates the
solution obtained with the repulsive force set to zero. The direct numerical
solution shows good convergence with the quasi-analytical solution as spatial
resolution increases. At least 64 cells are required to discretize the particle
diameter and accurately capture a symmetric trajectory before and after the
bounce. Beyond the results obtained for the original parallelepiped box, the
effect of computational domain size on particle dynamics is examined by
doubling the lateral dimensions in the **x** and **y** directions. The correspond-
ing results are presented in Fig. 11a under the label "expanded domain".
One spatial resolution $D/\Delta x = 16$ is investigated. During the initial phase
($t/t_{max} \le 0.5$), when the particles approach each other, the domain size has
negligible influence on the dynamics and on the estimation of their relative
positions. Slight deviations appear at later times ($t/t_{max} \ge 0.5$) during the
repulsive phase. Thus, the domain size induces minimal influence on particle
dynamics. It should be noted that, although the size of the computational
domain was doubled, the wall correction factor *K* in the drag model (Eq. (A.2) in Appendix A), which accounts for the confinement effect of the domain, was not adjusted. This confirms the robustness of the drag model presented in Appendix A for the test case under consideration. Figure 11b shows the time evolution of particle velocities. The blue color palette refers to the left particle while the red one concerns the right particle. As with the trajectory, a large number of cells in the particle diameters (at least 64) is required to capture the accurate trend of the quasi-analytical solution. On the other hand, whatever the space resolution, the symmetry between the left and right particles is well observed.

Figure 12 illustrates the time evolution of the amplitude of the repulsive force **F**<sub>rep</sub>. A mesh convergence study is proposed. For the *Analytical* solution

<figure><img src="image_18.png" /><figcaption>Figure 12: Time evolution of the repulsive force magnitude $F_{rep}$ between the two particles.</figcaption></figure>

(black solid curve), **F**<sub>rep</sub> is computed based on the surface definition of the disjoining pressure (Eq. (13)), whereas in DNS, the CSF formulation is applied with a volume regularization of the disjoining pressure. The conclusions regarding the repulsive force, including mesh convergence and the required number of cells per particle diameter, are consistent with those obtained from the trajectory analysis (particle position and velocity). This confirms the validity of the CSF formulation for the disjoining pressure in the DNS. It should be noted that in the most refined calculation ($D/\Delta x = 64$), the loss of smoothness in $F_{rep}$ around $t/t_{max} = 0.35$ and $t/t_{max} = 0.65$ arises
from artifacts in the repulsive force magnitude aligned with the grid diago-
nals (visible at $t/t_{max} = 0.6$ in Fig. 10c), which slightly perturb its evaluation
over the particle.

This test case, which allows comparison with a quasi-analytical solution,
validates the discretization of the repulsive force in the DNS code in particu-
lar through a CSF formulation for the disjoining pressure. The second-order
accurate Lax-Wendroff/Runge-Kutta 2 numerical scheme used to transport
the inter-particle distance from the ridge surfaces to the liquid-gas interfaces
demonstrates sufficient accuracy for the precise calculation of the disjoining
pressure at these interfaces. The following sections present two more physi-
cally realistic test cases: the bouncing between two deformable droplets and
the static equilibrium of a foam network. The calibration of parameters $P_0$
and $h^*$ will be addressed.

## 2.2. Dynamics of bouncing deformable particles

The model is further validated for more inertial regimes, specifically
through simulations of head-on and off-axis bouncing collisions between two
immiscible deformable droplets. In the original experiments by [6], the
droplets consisted of two different fluids: one water and the other diesel
oil. However, the level-set front-capturing method employed in the present
study is limited to representing a single liquid phase in addition to the gas
phase. As all impact parameters in [6] are defined with respect to diesel
oil, both droplets are modeled as diesel oil in the numerical simulations for
consistency. A similar assumption is also adopted in [27].

The impact between the two droplets is characterized by the Reynolds
number $Re = \rho_f U_{rel} D / \mu_f$, Weber number $We = \rho_f U_{rel}^2 D / \sigma$, and the im-
pact parameter $X$. Here, $D$, $\rho_f$, $\mu_f$ and $\sigma$ denote the diameter, density,
dynamic viscosity, and surface tension coefficient of the diesel oil droplet,
respectively. The corresponding values are $D = 750 \ \mu\text{m}$, $\rho_f = 817 \ \text{kg.m}^{-3}$,
$\mu_f = 3.16 \ \text{mPa.s}$ and $\sigma = 28.25 \ \text{mN/m}$. The relative impact velocity be-
tween the two droplets is denoted by $U_{rel}$. The impact parameter $X$, which
ranges from 0 to 1, is defined following [6] as the normalized perpendicular
distance from the center of one droplet to the line of relative velocity passing
through the center of the other droplet. In the frame of reference moving
with one of the droplets, $X = 0$ corresponds to a head-on collision, while
$X = 1$ represents the limiting case in which the droplets pass by each other
without making contact.
As described in [51], the cut-off length scale $h_c$ is chosen such that $h_c/D = 0.2$. The amplitude $P_0$ of the disjoining pressure is selected to satisfy the condition $P_0(\sigma/R)^{-1} = 22$, which enables droplet deformation due to repulsive forces during impact, despite the counteracting effect of surface tension. With respect to spatial resolution, a grid density of $D/\Delta x \approx 100$ is used, ensuring that $h_c/\Delta x \approx 20$. As discussed in previous sections, this resolution is sufficient to accurately capture the repulsive interactions.

Figures 13, 14 and 15 show configurations ($We = 9.6$, $Re = 129$, $X = 0.0$), ($We = 9.9$, $Re = 131$, $X = 0.33$) and ($We = 7.3$, $Re = 112$, $X = 0.82$) respectively. Numerical simulations (bottom) are compared with experiments (top and extracted from [6]) at different times. For numerical simulations,

<figure><img src="image_19.png" /><figcaption>Figure 13: $We = 9.6$, $Re = 129$, $X = 0.0$. Top sequence: experimental results [6]. Bottom sequence: numerical results.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 14: $We = 9.9$, $Re = 131$, $X = 0.33$. Top sequence: experimental results [6]. Bottom sequence: numerical results.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 15: $We = 7.3$, $Re = 112$, $X = 0.82$. Top sequence: experimental results [6]. Bottom sequence: numerical results.</figcaption></figure>

the interface is colored according to the magnitude of the repulsive force. In blue, the force is zero, in red, it is activated. For all three configurations, despite the high Weber numbers, the activation of the repulsive force prevents the droplets from collapsing. The dynamics, particularly the deformation and trajectories of drops after impact, are accurately captured by DNS. To characterize the spatial extent of the interfacial regions contributing to the repulsive forces, Fig. 16 presents the magnitude of the repulsive force along the liquid-gas interfaces within the plane defined by the velocity vectors of the two particles. The three cases in Figs. 13, 14 and 15 are shown at the

<figure><img src="image_22.png" /><figcaption>Figure 16: Magnitude of the repulsive force along the interface, shown in the plane defined by the velocity vectors of the two particles. (a): $We = 9.6$, $Re = 129$, $X = 0.0$ ($t = 1.1$ ms); (b): $We = 9.9$, $Re = 131$, $X = 0.33$ ($t = 1.1$ ms); (c): $We = 7.3$, $Re = 112$, $X = 0.82$ ($t = 1.2$ ms).</figcaption></figure>
same instant ~ 1.1 ms. As expected, repulsive forces arise solely at the facing interfaces.

The influence of the disjoining pressure magnitude $P_0$ is now examined. In [51], which investigates interactions between two bubbles, dividing the multiplicative constant in the repulsion model by a factor of 10 results in a transition from numerical bouncing to coalescence. A similar trend is observed in Fig. 17 for the configuration ($We = 9.9$, $Re = 131$, $X = 0.33$). The bottom sequence shows DNS results obtained with $P_0(\sigma/R)^{-1} = 22$, consistent with the results presented in Fig. 14, while the top sequence corresponds to simulations performed with a reduced disjoining pressure magnitude, $\tilde{P}_0 = P_0/10$. At early times ($t \le 0.3$ ms), when the two droplets remain

<figure><img src="image_23.png" /><figcaption>Figure 17: $We = 9.9$, $Re = 131$, $X = 0.33$. Top sequence (a): $\tilde{P}_0(\sigma/R)^{-1} = 2.2$ with droplet coalescence and the formation of a 48 µm liquid bridge discretized by 6 grid cells. Bottom sequence (b): $P_0(\sigma/R)^{-1} = 22$ with droplet bouncing.</figcaption></figure>

sufficiently separated, the value of $P_0$ has no influence on their dynamics. For longer times ($t > 0.3$ ms), as the droplets get closer, the reduction in $P_0$ leads to a transition from a bouncing-dominated regime to coalescence, characterized by the formation of a liquid bridge at $t = 0.4$ ms, in agreement with [51].

The disjoining pressure model offers an effective framework for capturing droplet impact dynamics, particularly the bouncing phenomenon, which
remains challenging for interface-capturing methods based on the level-set approach. The conclusions derived from binary droplet collisions can be readily extended to droplet populations without additional computational cost, owing to the Eulerian formulation used in reconstructing inter-particle distances. In accordance with previous studies, the parameter $h^*$ is set such that $5h^*/D = 0.2$. The pressure scaling parameter $P_0$ is chosen so that $P_0(\sigma/R)^{-1} \sim 20$. Decreasing $P_0$ by an order of magnitude shifts the impact behavior from bouncing to coalescence, in agreement with the findings reported in [51].

## 2.3. Polyhedral structure of a monodisperse bubble network in foam

The third test case focuses on validating the interface shape under the competing effects of repulsive forces and surface tension. A foam-like configuration is considered, consisting of a monodisperse, ordered bubble network with a low liquid volume fraction (Fig. 18, left). For simplicity, the analysis is restricted to a two-dimensional configuration, though it can be readily extended to three dimensions. Taking the distance between bubble centers $L_c$ as the reference length scale, the computational domain is defined as a rectangular box of dimensions $[L_x, L_y] = [1, \sqrt{3}]$, discretized using a grid of $n_x = 296$ and $n_y = 512$ cells, with periodic boundary conditions.

Four circular half-bubbles, each with an initial radius of $R_0 = 0.48$, are positioned at the center of each domain boundary (Fig. 18, left). The liquid volume fraction is set to $\phi_l = 0.164$. The gas density and dynamic viscosity are given by $\rho_{gas} = 1 \text{ kg.m}^{-3}$ and $\mu_{gas} = 1.85 \ 10^{-5} \text{ Pa.s.}$, respectively. The density and viscosity ratios are $\rho_{liquid}/\rho_{gas} = 2$ and $\mu_{liquid}/\mu_{gas} = 2$, where $\rho_{liquid}$ and $\mu_{liquid}$ represent the density and viscosity of the liquid phase, respectively. The surface tension coefficient is $\sigma = 10^{-2} \text{ N/m}$, and both the gas and liquid phases are initially at rest. The repulsive interaction parameters are set to $h^* = 0.05$, resulting in $h_c/(2R_0) = 0.26$ which is comparable to values commonly used in previous studies [51], and $P_0(\sigma/R)^{-1} = 4.8$, ensuring that the disjoining pressure and surface tension are of comparable magnitude. The resulting interface shape (Fig. 18, right) arises from the local balance between these two competing effects.

By defining $R_0$ as the characteristic length scale and $P_0$ as the scaling factor for the disjoining pressure and surface tension, the corresponding time scale is introduced as $\tau = \sqrt{\rho_g R_0^2 / P_0}$. The temporal evolution of the foam network is obtained through direct numerical simulation (DNS) and illustrated in Fig. 18. Liquid regions are shown in black, while the gas phase is
<figure><img src="image_24.png" /><figcaption>Figure 18: Density contour at $t/\tau = 0, 0.5, 1.0, 1.5$. Black: high density, lightgrey: low density.</figcaption></figure>

indicated in white. The interface gradually evolves into a hexagonal configuration, reaching a steady state at $t \gtrsim 1$, thereby validating the chosen scaling for the characteristic time scale.

To validate the numerical solution, a quasi-exact theoretical solution is constructed. A central stationary bubble (i) is surrounded by six uniformly distributed neighboring bubbles (ii-vii), as illustrated in Fig. 19. In a polar coordinate system $(r, \theta)$ centered on bubble (i), the interface is described by the radial distance $R(\theta)$ from the center of bubble (i). The local curvature $\kappa(\theta)$ is given by:

$$ \kappa(\theta) = \frac{2 \left(\frac{dR}{d\theta}\right)^2 + R(\theta)^2 - R(\theta) \frac{d^2 R}{d\theta^2}}{\left[\left(\frac{dR}{d\theta}\right)^2 + R(\theta)^2\right]^{3/2}}, \quad (16) $$

while the inter-particle distance $h(R, \theta)$ between bubbles (i) and (ii) can be written as:

$$ h(\theta) = L_c - 2R(\theta) \cos(\theta). \quad (17) $$

Up to a phase shift in the angle $\theta$, the inter-particle distances between bubble (i) and the six neighboring bubbles can be computed in the same way.

Since a stationary solution is sought, the pressure is assumed to be constant in each of the phases, while local equilibrium at the interface is given by:

$$ P_{gas}(\theta) - P_{liq}(\theta) = \sigma\kappa(\theta) + P_0 \exp[-h(\theta)/h^*] \quad (18) $$
<figure><img src="image_25.png" /><figcaption>Figure 19: Setup of the quasi-exact solution.</figcaption></figure>

where $P_{gas}(\theta)$ and $P_{liq}(\theta)$ are independent of $\theta$ and are denoted simply as $P_{gas}$ and $P_{liq}$, respectively. The shape of the bubbles is then governed by the constant sum of capillary and disjoining pressures at each angular position, given that both the liquid and gas pressures ($P_{liq}$ and $P_{gas}$ respectively) remain constant. The resulting forces act in the same direction, namely toward the interior of the bubbles. By integrating Eq. (18) with respect to $\theta$ over the interval $[0; 2\pi]$, it can be reformulated as:

$$ \sigma\kappa(\theta) + P_0 \exp[-h(\theta)/h^*] = \frac{1}{2\pi} \int_{0}^{2\pi} (\sigma\kappa(\theta) + P_0 \exp[-h(\theta)/h^*]) d\theta \quad (19) $$

thereby eliminating any reference to the pressure jump across the interface. The system composed of Eqs. (16), (17) and (19), with $R(\theta)$ as the unknown, has no analytical solution and must therefore be solved numerically. To this end, the derivatives in Eq. (16) are approximated using a fourth-order finite difference scheme, while Eq. (19) is solved using Powell's hybrid method (Dog Leg algorithm).

Figure 20 shows a comparison between Direct Numerical Simulation (DNS) results and the quasi-exact stationary solution for a monodispersed, ordered bubble network. Figure 20a illustrates that the interface position and its polyhedral geometry are accurately captured by the DNS. Figure 20b presents the relative error in the radius computed from DNS compared to the analytical solution. The error remains on the order of 1%, demonstrat-
<figure><img src="image_26.png" /><figcaption>Figure 20: Comparison between Direct Numerical Simulation (DNS) results and the quasi-exact stationary solution for a monodispersed, ordered bubble network.</figcaption></figure>

ing the method’s accuracy and its capability to realistically represent a foam network. The asymmetry of the $\varepsilon$ error in Fig. 20(b) between $\theta = -2.6$ and $\theta = 2.6$ results from reconstructing the $\phi = 0$ isocontour by interpolation from the DNS field of the level-set function, defined at the mesh cell centers. The radius $R_{DNS}$ is obtained for each $\theta$ from this isocontour. This asymmetry occurs only for small error values and does not affect the upper bound of the error.

The proposed numerical method, which incorporates a repulsive force calibrated using $h_c/D \sim 0.2$ and $P_0(\sigma/R)^{-1} \sim 5$, successfully captures the characteristic polygonal structures observed in static foam networks. This approach enables DNS of wet foams at reduced liquid fractions (from 0.30 in [26] to 0.16 in the present study), thereby permitting the simulation of more realistic wet foam conditions.

## 2.4. Polydisperse bubble arrangement

The final test case focuses on the interactions between bubbles of varying sizes within a polydisperse foam. In a square periodic computational domain of size $L = 1$, discretized into $n_x = n_y = 256$ grid cells, six bubbles interact solely under the effects of surface tension and disjoining pressure. The initial configuration is shown in Fig. 21(a), while the initial radii and positions of the six bubbles are listed in Tab 1. The bubble index is defined in Fig. 21(a). The physical parameters are identical to those of the monodisperse test case (Sec. 2.3), and the disjoining pressure amplitude $P_0$ is chosen such
<table><caption>Table 1: Initial positions $(x_i, y_i)$ and radii $R_i$ of the bubbles in the polydisperse test case. Bubble numbering (*index*) is shown in Fig. 21(a). The ratio $P_0 R_i / \sigma$ is also given for each bubble.</caption><thead><tr><td>index</td><td>x<sub>i</sub>-coordinate</td><td>y<sub>i</sub>-coordinate</td><td>radius R<sub>i</sub></td><td>P<sub>0</sub>R<sub>i</sub>/σ</td></tr></thead><tbody><tr><td>1</td><td>0.39222</td><td>0.39222</td><td>0.38</td><td>19</td></tr><tr><td>2</td><td>0.80778</td><td>0.80778</td><td>0.18</td><td>9</td></tr><tr><td>3</td><td>0.855</td><td>0.145</td><td>0.12</td><td>6</td></tr><tr><td>4</td><td>0.145</td><td>0.855</td><td>0.12</td><td>6</td></tr><tr><td>5</td><td>0.875</td><td>0.518</td><td>0.1</td><td>5</td></tr><tr><td>6</td><td>0.518</td><td>0.875</td><td>0.1</td><td>5</td></tr></tbody></table>

that the ratios $P_0 R_i / \sigma > 5$ for each bubble (see Tab. 1). As shown in the monodisperse test case, this condition ensures that the repulsive forces are sufficient to prevent the collapse of two opposing bubbles.

Figure 21 presents the temporal evolution of the six bubbles. As in the monodisperse case, the characteristic time $\tau = \sqrt{\rho_g R_1^2 / P_0}$ is defined with respect to the initial radius of the largest bubble $R_1$, which represents the longest and most demanding timescale among the bubble sizes considered. Although no theoretical solution is available for the stationary state of this polydisperse foam network, the arrangement of bubbles of different sizes observed in Fig. 21 shows similarities to the wet foams reported in [19]. The smallest bubbles initially deform into a triangular-like shape under the repulsive force induced by disjoining pressure, but evolve into a quasi-spherical state by the end of the simulation as surface tension dominates over repulsion. The largest bubble first develops a dodecahedral-like configuration due to the repulsive action of its twelve neighbors, and eventually adopts a quasi-square shape. Throughout the temporal evolution, opposing interfaces tend to remain parallel. While the bubble shapes appear qualitatively consistent, future work will be required to quantitatively validate the morphology of the polydisperse network, for instance through comparison with experimental data or analytical solutions.

### 3. Conclusions

The proposed mesoscale approach is based on a disjoining pressure conceived as an effective coarse-grained representation of nanoscale near-contact
<figure><img src="image_27.png" /><figcaption>Figure 21: Density contour at different times. Black: high density, lightgrey: low density. Polydisperse test case.</figcaption></figure>
interactions like short-range forces. This formulation enables the consistent description of macroscopic interfacial dynamics without explicit resolution of the underlying microscopic forces. The disjoining pressure is modeled as a decreasing exponential function of the interfacial distance, defined by two constants. The maximum pressure is estimated from capillary forces at a characteristic length scale, such as the droplet or bubble radius, or ideally the Plateau border radius in foams. The ratio $P_0(\sigma/R)^{-1}$ varies between 5 and 20, covering a range of configurations including droplet interaction and static foam networks. The cut-off length, representing the range of influence of the repulsive force, is set to 20% of the inclusion diameter, as recommended in the literature.

The method is formulated within an Eulerian framework, wherein the distance between facing interfaces is initially computed at the ridge surfaces and subsequently advected to the liquid-gas boundaries. This approach is inherently well-suited for parallel computing architectures. Moreover, its computational cost remains moderate—comparable to that of a level-set reinitialization step—enabling efficient implementation in large-scale simulations. However, the main limitations of the method are as follows. First, unlike Lagrangian approaches, it does not guarantee the prevention of coalescence between two facing interfaces when the relative velocity of interacting particles is high. Second, the current implementation is not capable of capturing particle interactions with boundary conditions such as solid walls, and, in particular, does not account for the presence of the three-phase contact line.

Beyond its computational advantages, the proposed method exhibits a wide range of applicability. It accurately captures both dynamic inertial interactions such as droplet bouncing and head-on collisions (which classical level-set methods cannot resolve without incorporating repulsive forces), as well as quasi-static bubble configurations governed by the balance between capillary and repulsive forces. In particular, the method is able to reproduce the polyhedral geometries characteristic of foam networks in equilibrium, where disjoining pressure plays a crucial role in maintaining structural integrity at low liquid fractions. The method allows reduced liquid fractions from 0.3 to 0.16 in level-set based DNS of wet foams, thereby enabling the computation of increasingly drier foam networks. Upcoming work will focus on optimizing $P_0$ and $h^*$ to simulate even drier foams while preventing interface collapse.

Future studies will focus on investigating foam destabilization mechanisms
nisms, such as coarsening, drainage, and rupture, with the aim of improving the predictive modeling of complex multiphase systems in both natural and industrial settings.

### CRediT authorship contribution statement

**Théo Rajim:** Software, Validation, Visualization. **Yedhir Mezache:** Conceptualization, Methodology. **Bastien Di Pierro:** Conceptualization, Methodology, Software, Validation, Formal analysis, Writing - Original Draft, Writing - Review & Editing, Visualization, Supervision. **Pierre Trontin:** Methodology, Software, Validation, Formal analysis, Writing - Original Draft, Writing - Review & Editing, Visualization, Supervision.

### Acknowledgement

This work was partially supported by the LABEX iMUST of the University of Lyon (ANR-10-LABX-0064), created within the Plan France 2030 set up by the French government and managed by the French National Research Agency (ANR). Access was granted to the HPC resources of PMCS2I (Pôle de Modélisation et de Calcul en Sciences de l'Ingénieur et de l'Information) at Ecole Centrale de Lyon, P2CHPD (Pôle de Compétence en Calcul Haute Performance Dédié) at Université Lyon 1, and GENCI-IDRIS (Grant 2022-A0132B13793).

## Appendix A. Quasi-analytical solution for particle trajectories

The geometrical configuration is as shown in Fig. 8. Since particle density ($\rho_p = 997$ kg.m$^{-3}$) is much greater than that of the surrounding fluid, the effects of virtual mass, Basset, pressure gradient (buoyancy) and lift forces can be neglected compared to those of the drag force. So, during their movement, gravity being assumed equal to zero, the particles are subjected to the drag force $F_d$ and the repulsive force $F_{rep}$, depending on the distance between the two particles. The position $X$ and velocity $u$ of each particle are then solutions of the following system of equations:

$$ \frac{dX}{dt} = u \qquad (A.1a) $$

$$ \frac{\pi}{6} \rho_p D^3 \frac{du}{dt} = F_d + F_{rep} \qquad (A.1b) $$
By symmetry of the problem, only the projection along **z** of Eq. (A.1) is solved. The estimate of $F_{rep} = \mathbf{F}_{rep} \cdot \mathbf{z}$ is obtained from Eq. (13).

Since the particle Reynolds and Weber numbers are much smaller than 1, Stokes' law for viscous drag on a non-deformable spherical particle in a quiescent fluid is used to define $\mathbf{F}_d$:

$$ \mathbf{F}_d = -3K\pi\mu_f D\mathbf{u} \qquad (A.2) $$

$K$ is a dimensionless wall correction factor for rigid spheres moving in a still liquid inside a bounded domain [14]. $K = 1$ corresponds to the case with an unbounded fluid. To determine $K$, the motion of a single particle subject to the only drag force $\mathbf{F}_d$ is investigated using DNS. The geometrical domain and physical parameters are those shown in Fig. 8. The time decay of the particle's velocity is then written:

$$ \ln(\frac{u}{u_0}) = -\frac{t}{\tau} \qquad (A.3) $$

where $\tau = \frac{\rho_v D^2}{18K\mu_f}$. Figure A.22 presents the time evolution of $\ln(\frac{u}{u_0})$ for both a single particle (blue shading) and two particles facing each other (red shading), allowing for an assessment of potential interactions affecting the drag force calculation. A grid convergence study is conducted for different $D/\Delta x$ ratios, with $\Delta x$ the space step. The results indicate that a converged solution is achieved for $D/\Delta x = 32$ in the single-particle case and $D/\Delta = 64$ for two head-on particles. Considering either a single particle or two head-on particles has no significant effect on particle dynamics or drag force estimation. From Fig. A.22 (dashed black line), $t_{max}/\tau \approx 0.39$, yielding $K \approx 2$.

Having defined the models for $\mathbf{F}_d$ and $\mathbf{F}_{rep}$, Eqs. (A.1) are numerically integrated using a 2nd-order accurate Runge-Kutta (RK2) method to derive particle position and velocity.

## References

[1] Benzi, R., Chibbaro, S., Succi, S., 2009. Mesoscopic lattice boltzmann modeling of flowing soft systems. Physical review letters 102, 026002.

[2] Bergeron, V., 1999. Forces and structure in thin liquid soap films. Journal of Physics: Condensed Matter 11, R215.
<figure><img src="image_28.png" /><figcaption>Figure A.22: Time decay of the velocity of a particle subjected to a drag force (*Single particle*). The case of two particles facing each other and subjected to drag force alone is also studied (*Two head-on particles*).</figcaption></figure>

[3] Brackbill, J., Kothe, D., Zemach, C., 1992. A continuum method for modeling surface tension. Journal of Computational Physics 100, 335-354.

[4] Camarasa, E., Vial, C., Poncin, S., Wild, G., Midoux, N., Bouillard, J., 1999. Influence of coalescence behaviour of the liquid and of gas sparging on hydrodynamics and bubble characteristics in a bubble column. Chemical Engineering and Processing: Process Intensification 38, 329-344.

[5] Cantat, I., Cohen-Addad, S., Elias, F., Graner, F., Höhler, R., Pitois, O., Rouyer, F., Saint-Jalmes, A., 2013. Foams: Structure and Dynamics. Oxford.

[6] Chen, R.H., Chen, C.T., 2006. Collision between immiscible drops with large surface tension difference: diesel oil and water. Experiments in fluids 41, 453-461.

[7] Chorin, A.J., 1968. Numerical solution of the navier-stokes equations. Mathematics of computation 22, 745-762.
[8] Cifani, P., Kuerten, J., Geurts, B., 2020. Flow and bubble statistics of turbulent bubble-laden downflow channel. International journal of multiphase flow 126, 103244.

[9] Couderc, F., 2007. Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles. Application à la désintégration assistée d'un jet liquide par un courant gazeux. Thèses. Ecole nationale supérieure de l'aeronautique et de l'espace. URL: https://theses.hal.science/tel-00143709.

[10] Coyajee, E., Boersma, B.J., 2009. Numerical simulation of drop impact on a liquid-liquid interface with a multiple marker front-capturing method. Journal of Computational Physics 228, 4444-4467.

[11] Deike, L., 2022. Mass transfer at the ocean-atmosphere interface: the role of wave breaking, droplets, and bubbles. Annual Review of Fluid Mechanics 54, 191-224.

[12] Dollet, B., Marmottant, P., Garbin, V., 2019. Bubble dynamics in soft and biological matter. Annual Review of Fluid Mechanics 51, 331-355.

[13] Elghobashi, S., 2019. Direct numerical simulation of turbulent flows laden with droplets or bubbles. Annual Review of Fluid Mechanics 51, 217-244.

[14] Haberman, W., Sayre, R., 1958. Motion of rigid and fluid spheres in stationary and moving liquids inside cylindrical tubes. Technical Report 1143. Department of the Navy. David Taylor Model Basin.

[15] Harlow, F.H., Welch, J.E., et al., 1965. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Physics of fluids 8, 2182.

[16] Hill, C., Eastoe, J., 2017. Foams: From nature to industry. Advances in colloid and interface science 247, 496-513.

[17] Hua, J., Stene, J.F., Lin, P., 2008. Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method. Journal of Computational Physics 227, 3358-3382.
[18] Innocenti, A., Jaccod, A., Popinet, S., Chibbaro, S., 2021. Direct numerical simulation of bubble-induced turbulence. Journal of Fluid Mechanics 918, A23.

[19] Jing, Z., Feng, C., 2023. Effect of polydispersity on the structural characteristics of two-phase foam. International Journal of Multiphase Flow 164, 104465.

[20] Kralchevsky, P.A., Ivanov, I.B., 1990. Micromechanical description of curved interfaces, thin films, and membranes: ii. film surface tensions, disjoining pressure and interfacial stress balances. Journal of colloid and interface science 137, 234-252.

[21] Kwakkel, M., Breugem, W.P., Boersma, B.J., 2013. Extension of a CLSVOF method for droplet-laden flows with a coalescence/breakup model. Journal of Computational Physics 253, 166-188.

[22] Lewin-Jones, P., Lockerby, D.A., Sprittles, J.E., 2024. Collision of liquid drops: bounce or merge? Journal of Fluid Mechanics 995, A1.

[23] Li, J., 2016. Macroscopic model for head-on binary droplet collisions in a gaseous medium. Physical review letters 117, 214502.

[24] Liu, M., Bothe, D., 2019. Toward the predictive simulation of bouncing versus coalescence in binary droplet collisions. Acta Mechanica 230, 623-644.

[25] Liu, X.D., Osher, S., Chan, T., 1994. Weighted essentially non-oscillatory schemes. Journal of computational physics 115, 200-212.

[26] Mezache, Y., Detcheverry, F., Di Pierro, B., Spelt, P.D., Biance, A.L., Le Merrer, M., 2024. Surfactant-induced dissipation in sheared foams: mechanics and thermodynamics. Physical Review Fluids 9, 094008.

[27] Montessori, A., Lauricella, M., Tirelli, N., Succi, S., 2019. Mesoscale modelling of near-contact interactions for complex flowing interfaces. Journal of Fluid Mechanics 872, 327-347.

[28] Olmos, E., Gentric, C., Vial, C., Wild, G., Midoux, N., 2001. Numerical simulation of multiphase flow in bubble column reactors. Influence of bubble coalescence and break-up. Chemical engineering science 56, 6359-6365.
[29] Osher, S., Sethian, J.A., 1988. Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations. Journal of computational physics 79, 12–49.

[30] Pelusi, F., Sbragaglia, M., Benzi, R., 2019. Avalanche statistics during coarsening dynamics. Soft matter 15, 4518–4524.

[31] Peskin, C.S., 2002. The immersed boundary method. Acta numerica 11, 479–517.

[32] Puckett, E.G., Almgren, A.S., Bell, J.B., Marcus, D.L., Rider, W.J., 1997. A high-order projection method for tracking fluid interfaces in variable density incompressible flows. Journal of computational physics 130, 269–282.

[33] Rubio, J., Souza, M., Smith, R., 2002. Overview of flotation as a wastewater treatment technique. Minerals engineering 15, 139–155.

[34] Scardovelli, R., Zaleski, S., 1999. Direct numerical simulation of free-surface and interfacial flow. Annual review of fluid mechanics 31, 567–603.

[35] Shan, X., Chen, H., 1993. Lattice boltzmann model for simulating flows with multiple phases and components. Physical review E 47, 1815.

[36] Shin, S., Juric, D., 2002. Modeling three-dimensional multiphase flow using a level contour reconstruction method for front tracking without connectivity. Journal of Computational Physics 180, 427–470.

[37] Sprittles, J.E., 2024. Gas microfilms in droplet dynamics: When do drops bounce? Annual Review of Fluid Mechanics 56, 91–118.

[38] Stevenson, P., 2012. Foam engineering: fundamentals and applications. John Wiley & Sons.

[39] Sussman, M., Fatemi, E., 1999. An efficient, interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow. SIAM Journal on scientific computing 20, 1165–1191.

[40] Sussman, M., Puckett, E.G., 2000. A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows. Journal of computational physics 162, 301–337.
[41] Sussman, M., Smereka, P., Osher, S., 1994. A level set approach for computing solutions to incompressible two-phase flow. Journal of Computational physics 114, 146–159.

[42] Tatebe, O., 1996. MGCG method: a robust and highly parallel iterative method. Ph.D. thesis. University of Tokyo.

[43] Temam, R., 1969. Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (II). Archive for rational mechanics and analysis 33, 377–385.

[44] Trontin, P., Vincent, S., Estivalezes, J., Caltagirone, J., 2010. Direct numerical simulation of a freely decaying turbulent interfacial flow. International Journal of Multiphase Flow 36, 891–907.

[45] Trontin, P., Vincent, S., Estivalezes, J.L., Caltagirone, J.P., 2008. Detailed comparisons of front-capturing methods for turbulent two-phase flow simulations. International journal for numerical methods in fluids 56, 1543–1549.

[46] Trottenberg, U., Oosterlee, C.W., Schuller, A., 2001. Multigrid methods. Academic press.

[47] Tryggvason, G., Bunner, B., Esmaeeli, A., Juric, D., Al-Rawahi, N., Tauber, W., Han, J., Nas, S., Jan, Y.J., 2001. A front-tracking method for the computations of multiphase flow. Journal of computational physics 169, 708–759.

[48] Unverdi, S.O., Tryggvason, G., 1992. A front-tracking method for viscous, incompressible, multi-fluid flows. Journal of computational physics 100, 25–37.

[49] Youngs, D., 1982. Time-Dependent Multi-Material Flow with Large Fluid Distortion. Numerical Methods for Fluid Dynamics/Academic Press, pp. 273-285.

[50] Zhang, A., Guo, Z., Wang, Q., Xiong, S., 2019. Three-dimensional numerical simulation of bubble rising in viscous liquids: A conservative phase-field lattice-boltzmann study. Physics of Fluids 31.
[51] Zhang, L., Peng, K., Shao, X., Deng, J., 2022a. Short-range repulsive force model for near-contact interactions of bubbles. Physical Review E 106, 045110.

[52] Zhang, S., Tang, J., Wu, H., 2022b. Phase-field lattice boltzmann model for two-phase flows with large density ratio. Physical Review E 105, 015304.

[53] Zuzio, D., Estivalèzes, J.L., Di Pierro, B., 2018. An improved multi-scale eulerian-lagrangian method for simulation of atomization process. Computers & Fluids 176, 285-301.

[54] Zuzio, D., Orazzo, A., Estivalèzes, J.L., Lagrange, I., 2020. A new efficient momentum preserving level-set/vof method for high density and momentum ratio incompressible two-phase flows. Journal of Computational Physics 410, 109342.