1. Introduction
Oil–water emulsions containing solid particles are encountered in numerous industrial processes, including enhanced oil recovery (Adil & Onaizi Reference Adil and Onaizi2022), oil sands processing (Velayati et al. Reference Velayati, Habibi, Nikrityuk, Tang and Zeng2022), food and cosmetic formulations (Yang et al. Reference Yang, Fang, Chen, Zhang, Xie, Chen, Liu and Yuan2017), microplastics separation (Nabi et al. Reference Nabi, Bacha and Zhang2022) and mineral extraction (Kusaka et al. Reference Kusaka, Arimoto, Nakahiro and Wakamatsu1994). In these systems, micron- and sub-micron-sized particles can either adsorb at the oil–water interface – forming Pickering emulsions (Pickering Reference Pickering1907) – or remain dispersed within one of the bulk phases. This behaviour is primarily governed by particle wettability, i.e. the preferential affinity of particles for either the water or oil phase, which classifies them as hydrophilic or lipophilic. Particle wettability is therefore a key factor in determining particle trapping, separation, and transport efficiency in oil–water systems.
Particle trapping at quiescent fluid interfaces is relatively well understood (Martínez-Pedrero et al. Reference Martínez-Pedrero, Carbone, Rubio, Ortega and Guzmán2024). By contrast, under turbulent conditions, the interactions between particles and interfaces remain poorly characterised, particularly when interface deformation, break-up and coalescence occur simultaneously (Outokesh et al. Reference Outokesh, Saeedipour and Hlawitschka2026). This knowledge gap stems from the complexity of multiphase turbulence, in which particles interact with droplets that continuously break up and coalesce, undergoing repeated attachment to and detachment from interfaces in response to hydrodynamic fluctuations (Hajisharifi et al. Reference Hajisharifi, Marchioli and Soldati2021; Hajisharifi et al. Reference Hajisharifi, Marchioli and Soldati2022). As a result, the overall dynamics is intertwined and difficult to disentangle. In this context, numerical simulations represent a powerful tool to investigate these processes, as they provide detailed space- and time-resolved information that enables the identification and quantification of particle–interface interaction mechanisms under turbulent conditions.
Fully resolved simulations, in which all relevant scales are captured from first principles, remain computationally prohibitive, as particle–interface interactions involve physical mechanisms spanning a wide range of scales (Kostoglou et al. Reference Kostoglou, Karapantsios and Matis2006; Soligo et al. Reference Soligo, Roccon and Soldati2021). Consequently, a variety of modelling approaches have been developed, which differ in how interfaces and particles are represented, and in which physical mechanisms are explicitly resolved (Singh & Joseph Reference Singh and Joseph2005; Joshi & Sun Reference Joshi and Sun2010; Aland et al. Reference Aland, Lowengrub and Voigt2011; Jansen & Harting Reference Jansen and Harting2011; Millett & Wang Reference Millett and Wang2011; Gao et al. Reference Gao, Evans, Wanless and Atanasio2017; Zeng et al. Reference Zeng, Lu and Tryggvason2025). In this work, we adopt an interface-resolved approach coupled with a point-particle description. This strategy captures the essential features of turbulent multiphase flows while remaining computationally affordable, and enables simulations with large particle populations, which would be otherwise prohibitive for fully resolved particle–interface interactions requiring extremely fine grids to resolve lubrication films and contact-line dynamics. We consider a turbulent oil–water emulsion laden with sub-Kolmogorov particles in a channel flow. The oil–water interface is captured using a phase-field method, while particle dynamics is described using a Lagrangian framework that accounts for interfacial adhesion and wettability effects. Building on previous approaches (Gu & Botto Reference Gu and Botto2020; Hajisharifi et al. Reference Hajisharifi, Marchioli and Soldati2021; Hajisharifi et al. Reference Hajisharifi, Marchioli and Soldati2022), we introduce a wettability-dependent particle–interface interaction that allows particles to exhibit tunable wettability (hydrophilic or lipophilic). Using this framework, we perform original simulations that resolve the coupled dynamics of turbulence, interface deformation and asymmetric particle capture. In the present work, the flow configuration is intentionally kept fixed in order to isolate the specific role of particle wettability on interfacial trapping dynamics, independently of additional variations in turbulence and flow parameters. We show that the resulting trapping dynamics can be described by a first-order kinetic balance between attachment and detachment processes.
2. Governing equations
We perform direct numerical simulations (DNS) of an incompressible multiphase flow consisting of an oil-in-water emulsion in a channel flow configuration, laden with sub-Kolmogorov particles characterised by different wettability. Following Tryggvason et al. (Reference Tryggvason, Dabiri, Aboulhasanzadeh and Lu2013), the term DNS refers here to the fluid-phase simulations, in which all continuum time and length scales are fully resolved across the complete range of scales present in the unsteady flow. The dynamics of the oil–water interface is captured by a phase-field method (Jacqmin Reference Jacqmin1999; Badalassi et al. Reference Badalassi, Ceniceros and Banerjee2003; Roccon et al. Reference Roccon, Zonta and Soldati2023), while the motion of the sub-Kolmogorov particles is described using a point-particle approach (Brandt & Coletti Reference Brandt and Coletti2022; Marchioli et al. Reference Marchioli2025). The interaction of the Lagrangian particles with the oil–water interface is modelled by introducing additional forces that account for adhesion and wettability effects (Gu & Botto Reference Gu and Botto2020; Hajisharifi et al. Reference Hajisharifi, Marchioli and Soldati2021; Iberl et al. Reference Iberl, Trautner, Klein and Hasslberger2026). The adopted modelling strategy makes it possible to simulate large particle populations while retaining the essential physical mechanisms governing interfacial trapping. A direct resolution of the interface deformation induced by the particles is currently beyond the available computational resources, owing to the large scale separation between particle–interface interactions and turbulent flow scales. The adopted modelling framework has been previously validated and applied to investigate particle capture and interfacial dynamics in turbulent multiphase flows (Gu & Botto Reference Gu and Botto2020; Hajisharifi et al. Reference Hajisharifi, Marchioli and Soldati2021; Hajisharifi et al. Reference Hajisharifi, Marchioli and Soldati2022; Iberl et al. Reference Iberl, Trautner, Klein and Hasslberger2026).
2.1. Phase-field method
The dynamics of the oil–water interface is captured by a phase-field method, which relies on a phase-field variable
$\phi$
(Jacqmin Reference Jacqmin1999; Badalassi et al. Reference Badalassi, Ceniceros and Banerjee2003; Roccon et al. Reference Roccon, Zonta and Soldati2023). This variable has uniform values in the bulk of the two phases –
$\phi =-1$
in the water carrier, and
$\phi =+1$
inside the oil droplets – and transitions smoothly across a thin transition layer. The characteristic length scale of this layer is controlled by the Cahn number
$\textit{Ch}=\xi /h$
, the ratio between its dimensional value and the channel half-height
$h$
. The evolution of the phase-field variable is governed by a Cahn–Hilliard equation, which in dimensionless form reads
\begin{align} \frac {\partial \phi }{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\phi =\frac {1}{\textit{Pe}_{\phi }} {\nabla} ^2\mu _{\phi }+\frac {\lambda }{\textit{Pe}_{\phi }}\left [{\nabla} ^2\phi -\frac {1}{\sqrt {2}\, \textit{Ch}}\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\big (1-\phi ^2\big )\frac {\boldsymbol{\nabla }\phi }{|\boldsymbol{\nabla }\phi |}\right )\right ]\!, \end{align}
where the first term on the right-hand side represents the diffusive flux induced by the phase-field chemical potential
$\mu _\phi =\delta \mathcal{F}/\delta \phi =\phi ^3-\phi -\textit{Ch}^2 {\nabla} ^2 \phi$
, which is obtained from a Ginzburg–Landau free-energy functional (Roccon et al. Reference Roccon, Zonta and Soldati2023). The second term represents a regularisation term that is employed to better preserve the hyperbolic tangent profile of the transition layer during the computation (Li et al. Reference Li, Choi and Kim2016; Soligo et al. Reference Soligo, Roccon and Soldati2019). The Péclet number
$ \textit{Pe}_\phi =u_\tau h/(\mathcal{M}_\phi \beta )$
represents the ratio of convective to diffusive time scales, where
$u_\tau$
is the friction velocity,
$\mathcal{M}_\phi$
is the phase-field mobility, and
$\beta$
is a positive constant used in the dimensionless procedure. The parameter
$\lambda$
is a dimensionless numerical coefficient that sets the strength of the profile-correction term introduced to preserve the hyperbolic tangent shape of
$\phi$
within the interfacial transition layer during the simulation (Soligo et al. Reference Soligo, Roccon and Soldati2019).
2.2. Hydrodynamics
To describe the flow field, a one-fluid approach is employed, and a single set of Navier–Stokes and mass conservation equations is solved in the entire domain (Tryggvason et al. Reference Tryggvason, Scardovelli and Zaleski2011). We assume that the fluids have equal density (
$\rho _o=\rho _w=\rho$
) and viscosity (
$\mu _o=\mu _w=\mu$
), thus mimicking a realistic emulsion composed of water and light oil (e.g. ExxonMobil Isopar H). With these assumptions, the mass and Navier–Stokes equations read as
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{u}}=0 , \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }\left ({\boldsymbol{u}\otimes \boldsymbol{u}}\right )=-\boldsymbol{\nabla }p+\frac {1}{Re_\tau }{\nabla} ^2\boldsymbol{u}+\frac {3}{\sqrt {8}}\frac {\textit{Ch}}{\textit{We}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }_c, \\[0pt] \nonumber \end{align}
where
$\boldsymbol{u}$
is the velocity vector, and
$p$
is the pressure field. The shear Reynolds number is defined as
$Re_\tau =\rho u_{\tau }h/\mu$
, where
$\mu$
is the dynamic viscosity. Surface tension forces are introduced via a continuum surface stress formulation. Specifically, the resulting interfacial forces are computed by taking the divergence of the Korteweg stress tensor
$\boldsymbol{\tau }_c=\lvert \boldsymbol{\nabla }\phi \rvert ^2\,\boldsymbol{I}-\boldsymbol{\nabla }\phi \otimes \boldsymbol{\nabla }\phi$
, where
$\boldsymbol{I}$
is the identity matrix. The Weber number
$\textit{We}=\rho u_{\tau }^2h/\sigma$
quantifies the ratio of inertial to interfacial tension forces, with
$\sigma$
denoting the interfacial tension.
2.3. Solid spherical particles
The solid spherical particles are assumed to be smaller than any dynamically relevant flow scale, and are therefore modelled as material points within a point-particle framework, which also accounts for interface–particle interactions (Gu & Botto Reference Gu and Botto2020; Hajisharifi et al. Reference Hajisharifi, Marchioli and Soldati2021; Iberl et al. Reference Iberl, Trautner, Klein and Hasslberger2026). The motion of the particles is controlled by the balance equations
\begin{align} \frac {\text{d}\boldsymbol{x}_p}{\text{d}t}=\boldsymbol{u}_p,\quad m_p\frac {\text{d}\boldsymbol{u}_p}{\text{d}t}= \begin{cases} \boldsymbol{F}_D, & \left |\phi _p\right |\gt \phi _{\textit{lim}}, \\[5pt]\boldsymbol{F}_D+\boldsymbol{F}_A+\boldsymbol{F}_{\kern-1pt H}, & \left |\phi _p\right |\leq \phi _{\textit{lim}}, \end{cases} \end{align}
where
$\boldsymbol{x}_p$
and
$\boldsymbol{u}_p$
denote the particle position and velocity, respectively, and
$\phi _p$
is the phase-field value sampled at the particle position, used to determine whether the particle is located within the interfacial layer or in one of the bulk phases through the threshold value
$\phi _{\textit{lim}}=0.964$
, following Gu & Botto (Reference Gu and Botto2020). The particle mass is given by
$m_p=\pi d_p^3\rho _p/6$
, with
$d_p$
the particle diameter, and
$\rho _p$
its density. The hydrodynamic drag force
$\boldsymbol{F}_D$
is modelled using the Stokes drag law corrected for finite particle Reynolds number effects (Schiller & Naumann Reference Schiller and Naumann1933):
\begin{align} \boldsymbol{F}_D =3\pi d_p\mu \left (\boldsymbol{u}_{\kern-1pt f,p}-\boldsymbol{u}_p\right ) \left (1+0.15\,Re_p^{0.687}\right )\!, \end{align}
where
$\boldsymbol{u}_{\kern-1pt f,p}$
is the fluid velocity interpolated at the particle position, and
$Re_p=d_p\,|\boldsymbol{u}_{\kern-1pt f,p}-\boldsymbol{u}_p|/\nu$
is the particle Reynolds number. As a particle comes into contact with the fluid–fluid interface, additional interfacial effects influence its dynamics (Binks & Horozov Reference Binks and Horozov2006). In the present work, these effects are modelled through two forces: a capillary adhesion force
$\boldsymbol{F}_A$
, and a wettability-induced force
$\boldsymbol{F}_{\kern-1pt H}$
, accounting for energetic trapping at the interface and preferential phase affinity, respectively. Both forces are non-zero only within the thin interfacial layer, and are expressed in terms of the local phase-field value and its gradient, which respectively encode the signed distance from the
$\phi =0$
iso-surface and the interface normal vector. The capillary adhesion force is modelled following the fast interface particle interaction approach proposed by Gu & Botto (Reference Gu and Botto2020):
\begin{align} \boldsymbol{F}_A=-\mathcal{A}\pi \sigma \mathcal{D}\boldsymbol{n}\,, \end{align}
where the dimensionless parameter
$\mathcal{A}$
controls the adhesion strength,
$\boldsymbol{n}=\boldsymbol{\nabla }\phi /|\boldsymbol{\nabla }\phi |$
is the interface normal, and
$\mathcal{D}$
is the signed distance from the iso-level
$\phi =0$
. In the vicinity of the interface, the signed distance
$\mathcal{D}$
is obtained from the phase-field variable as
\begin{align} \mathcal{D}=\frac {\sqrt {2}}{2} \textit{Ch}\, \log \left (\frac {1+\phi _p}{1-\phi _p} \right )\!. \end{align}
Overview of simulation parameters.
Table 1. Long description
The table presents an overview of simulation parameters for three cases: Hydrophilic, Neutral, and Lipophilic. It consists of 3 rows and 9 columns. The columns are labeled as follows: Case, Reτ, We, Ch, Peφ, St+, ρp/ρ, dp+, A, and kH+. Each row corresponds to a different case with specific values for each parameter. Row 1: Case Hydrophilic, Reτ 300, We 3.0, Ch 0.02, Peφ 50, St+ 0.6, ρp/ρ 10, dp+ 1.04, A 0.1, kH+ -1720. Row 2: Case Neutral, Reτ 300, We 3.0, Ch 0.02, Peφ 50, St+ 0.6, ρp/ρ 10, dp+ 1.04, A 0.1, kH+ 0. Row 3: Case Lipophilic, Reτ 300, We 3.0, Ch 0.02, Peφ 50, St+ 0.6, ρp/ρ 10, dp+ 1.04, A 0.1, kH+ +1720.
Particle wettability is modelled as a differential energetic affinity for the two phases, consistent with the adsorption free-energy description introduced by Pieranski (Reference Pieranski1980), and adapted here to point particles within a diffuse-interface framework. Specifically, wettability is included through a particle–phase-field interaction potential proportional to the local phase-field value, giving rise to a force proportional to the phase-field gradient:
\begin{align} \boldsymbol{F}_{\kern-1pt H} = \pi k_{\kern-1pt H}\,\boldsymbol{\nabla }\phi \,, \end{align}
which acts normal to the interface. Positive values of
$k_{\kern-1pt H}$
correspond to lipophilic particles, negative values to hydrophilic particles, and
$k_{\kern-1pt H}=0$
to neutral wettability, while the magnitude of
$k_{\kern-1pt H}$
quantifies the strength of this preference. Because the adhesion and wettability forces are conservative, their combined effect can be expressed in terms of a scalar potential
$U$
, such that
$\boldsymbol{F}_A + \boldsymbol{F}_{\kern-1pt H} = -\boldsymbol{\nabla }U$
. In viscous units, this potential reads
\begin{align} U^+\left (\phi _p\right )=\frac {1}{2}\pi \mathcal{A}\frac {Re_{\tau }}{\textit{We}} {\mathcal{D}^+}^2 -\pi k_{\kern-1pt H}^+\phi _p+C\,, \end{align}
where
$k_{\kern-1pt H}^+=u_\tau k_{\kern-1pt H}/(\nu ^3\rho )$
is the dimensionless wettability parameter, and
$C$
is an integration constant. In the sharp-interface formulation of Pieranski (Reference Pieranski1980), the adsorption free energy includes a quadratic contribution in the particle distance from the interface, associated with capillary adhesion, and a linear contribution accounting for preferential wettability. In the present diffuse-interface formulation, the adhesion term remains quadratic in the signed distance, whereas the wettability contribution includes a smooth hyperbolic tangent dependence on the distance that recovers the linear behaviour locally near the interface. The competition between adhesion and wettability contributions determines the equilibrium particle position within the interfacial layer (figure 1
$a$
). Specifically, under quiescent conditions (
$\boldsymbol{F}_D=\mathbf{0}$
), the equilibrium particle position corresponds to the minimum of the potential well defined by (2.9). Once
$Re_\tau$
,
$\textit{We}$
,
$\textit{Ch}$
and the adhesion parameter
$\mathcal{A}$
are fixed, this equilibrium position is uniquely determined by the wettability parameter
$k_{\kern-1pt H}^+$
, as illustrated in figure 1(
$b$
).
Model of hydrophilic, neutral and lipophilic particle behaviour at the carrier–droplet interface. (a) Schematic representation of particle positioning within the interfacial layer as a function of the wettability parameter
$k_{\kern-1pt H}^+$
. Hydrophilic (blue) particles (
$k_{\kern-1pt H}^+\lt 0$
) preferentially reside on the carrier side (
$\phi \lt 0$
), and lipophilic (red) particles (
$k_{\kern-1pt H}^+\gt 0$
) on the droplet side (
$\phi \gt 0$
), while neutral (grey) particles (
$k_{\kern-1pt H}^+=0$
) localise symmetrically at the interface (
$\phi =0$
). Dashed contours represent
$\phi$
levels, and
$\mathcal{D}$
indicates the signed distance from the interface. (b) The corresponding free-energy profiles
$U^{+}(\phi )$
for different values of the wettability parameter
$k_{\kern-1pt H}^+$
, obtained at fixed
$Re_\tau = 300$
,
$\textit{We} = 3.0$
,
$\textit{Ch} = 0.02$
, and adhesion parameter
$\mathcal{A} = 0.1$
. Positive
$k_{\kern-1pt H}^+$
stabilises lipophilic particles on the droplet side, negative
$k_{\kern-1pt H}^+$
stabilises hydrophilic particles on the carrier side, and
$k_{\kern-1pt H}^+=0$
corresponds to neutral wettability with a minimum at the interface.
3. Simulations
The governing equations are solved using a pseudo-spectral method in a channel flow configuration, utilising the GPU-accelerated, open-source solver FLOW36 (Roccon et al. Reference Roccon, Soligo and Soldati2025). Simulations are performed at a shear Reynolds number
$Re_\tau =300$
. The domain is a channel of dimensions
$L_x \times L_y \times L_z = 4\pi h \times 2\pi h \times 2h$
, discretised using
$N_x\times N_y\times N_z=512\times 256\times 513$
grid points. To represent a realistic oil–water emulsion, the Weber number is set to
$\textit{We}=3.0$
. The Cahn number is set to
$\textit{Ch} = 0.02$
, compatibly with the grid resolution. The phase-field Péclet number is chosen as
$ \textit{Pe}_\phi =1/\textit{Ch}=50$
, and the parameter
$\lambda$
, which controls the magnitude of the profile-corrected term, is set using the scaling
$\lambda =0.0625/\textit{Ch}=3.125$
(Soligo et al. Reference Soligo, Roccon and Soldati2019). To investigate the role of wettability on particle–interface dynamics, we perform three simulations with different particle sets, characterised by different wettability parameters:
$k_{\kern-1pt H}^+ = - 1720$
(hydrophilic),
$k_{\kern-1pt H}^+=0$
(neutral) and
$k_{\kern-1pt H}^+= 1720$
(lipophilic), as summarized in table 1. Particles are simulated in a one-way coupling regime. The particle Stokes number is fixed at
$\textit{St}^+ = 0.6$
, corresponding to non-tracer, weakly inertial particles. The particle-to-fluid density ratio is set to
$\rho _p/\rho = 10$
, representative of typical metals or minerals suspended in water/oil, resulting in particle diameter
$d^+ = 1.04$
, which is smaller than the Kolmogorov scale in the channel, and therefore consistent with the point-particle modelling approach adopted in this work. The particle volume fraction is sufficiently small (
$\varPhi _p = 6.91 \times 10^{-5}$
) to justify the use of a one-way coupling regime, in which particles do not affect the carrier flow. The adhesion parameter is kept constant across all simulations at
$\mathcal{A} = 0.1$
. The value of
$\mathcal{A}$
is selected based on dynamical considerations of particle motion within the interfacial layer. In particular, the adhesion force acts as an elastic force normal to the interface, while viscous drag provides dissipation. For the particle size and density ratio considered in this work, the resulting dynamics is overdamped, such that particles relax monotonically towards their equilibrium position without oscillations. To obtain statistically stationary droplet–particle configurations, simulations are conducted in three stages. First, a single-phase turbulent channel flow is simulated until a statistically steady state is reached. Subsequently, the dispersed phase is introduced in the form of an array of
$N_d = 256$
spherical droplets. Each droplet has initial diameter
$d = 0.4h$
, corresponding to a dispersed-phase volume fraction
$\varPhi _d = 5.4\,\%$
. After the injection, droplets evolve freely under turbulent stresses, undergoing deformation, break-up and coalescence until reaching a statistically steady configuration characterised by a stable droplet count and interfacial area. At this stage, three particle sets, each containing
$N_0=5\times 10^5$
particles, are injected into the carrier phase. The system is then further evolved to collect particle–interface interaction statistics for the three wettability cases.
4. Results
We begin with a qualitative overview of the system configurations attained at statistical steady state. Figure 2 displays a vis-à-vis comparison of renderings of the system, showing the steady-state configurations observed for hydrophilic particles (figure 2
a) and lipophilic particles (figure 2
b). The droplet interface is identified by the iso-contour
$\phi =0$
. A comparison between the two cases immediately reveals that the particle trapping efficiency at the oil–water interface strongly depends on wettability. For hydrophilic particles (figure 2
a), a larger fraction remains dispersed in the carrier phase, indicating reduced trapping efficiency, whereas lipophilic particles (figure 2
b) are trapped more efficiently. Close-up views further show that hydrophilic particles preferentially accumulate on the carrier-fluid side of the interface (water phase), whereas lipophilic particles accumulate on the droplet side (oil phase). This preferential accumulation on one side of the interface is a direct consequence of the wettability force, which shifts the particle equilibrium position towards positive or negative values of the phase-field variable, and consequently towards the oil or water phase, as shown in figure 1(a).
Three-dimensional rendering of the steady-state configuration obtained for (
$a$
) hydrophilic and (
$b$
) lipophilic particles. The interface of the oil droplets is visualised by the
$\phi =0$
iso-contour. The first close-up view shows particle trapping by oil droplets, whereas the second close-up illustrates the equilibrium particle position at the interface. Lipophilic particles preferentially reside on the oil side of the interface, while hydrophilic particles reside on the water side. For visualisation purposes, Lagrangian particles are not shown to scale.
To quantify this latter observation, we compute the probability density function (PDF) of the phase-field values sampled by the particles. As the phase-field behaves as a signed distance (at least in the vicinity of the interface), this observable provides a direct indication of the instantaneous positions of particles trapped at the interface. The resulting PDFs are shown in figure 3(
$a$
) for the three cases considered: hydrophilic (blue), neutral-wettability (grey) and lipophilic (red). In all cases, the PDFs exhibit a peak at well-defined locations within the interfacial layer, corresponding to the minima of the wettability-dependent interaction potential (2.9), here reported with dashed vertical lines (see also figure 1). For hydrophilic particles (blue), the PDF exhibits a peak at negative values of
$\phi$
, which corresponds to positions located towards the water phase. Differently, the PDF of neutrally wetting particles peaks at the centre of the transition layer (
$\phi =0$
), and finally, for lipophilic particles, the PDF shows a peak in correspondence with positive values of
$\phi$
. The locations of these peaks indicate that particle configuration at the interface is controlled by wettability. Notably, for all cases, the PDFs are rather narrow, and almost all particles trapped at the interface remain confined in a region of the transition layer characterised by the same sign: hydrophilic particles preferentially sample negative values of
$\phi$
, while lipophilic particles sample positive values. Finally, the role of the adhesion force becomes apparent: by attracting the particle towards the centre of the transition layer, this force produces a positively skewed PDF for hydrophilic particles, and a negatively skewed PDF for lipophilic particles.
This wettability-controlled interfacial configuration leads to systematic variations in trapping efficiency, which increases with lipophilicity, i.e. for increasing values of
$k_{\kern-1pt H}^+$
. This trend is quantified in figure 3(
$b$
), which reports the fractions of particles trapped at the interface,
$N_i/N_0$
, and those dispersed in the carrier phase,
$N_c/N_0$
, averaged over the statistically steady regime (
$t^+ \gt 6000$
), as functions of the wettability and normalised by the total number of particles
$N_0$
. Trapping efficiency, which can be quantified by the fraction of particles trapped at the interface
$N_i/N_0$
, is asymmetric with respect to neutral wetting, and reaches its maximum for lipophilic particles (
$k_{\kern-1pt H}^+\gt 0$
). These trends can be interpreted in terms of attachment and detachment mechanisms. Turbulent fluctuations continuously transport particles towards droplet surfaces, promoting attachment to the interface. Once attached, particles detach only when local turbulent stresses overcome the finite energetic barrier associated with removal from the interface. Therefore, particle retention results from the competition between turbulent forcing and interfacial detachment energy. Although the energetic barrier for detachment is symmetric for lipophilic particles detaching towards the droplet interior, and hydrophilic particles released into the carrier phase, the turbulent stresses sampled by the particles on the two sides of the interface differ remarkably. In particular, turbulence intensity is significantly attenuated inside droplets (Scarbolo & Soldati Reference Scarbolo and Soldati2013), with corresponding modifications to the flow topology (Soligo et al. Reference Soligo, Roccon and Soldati2020). Consequently, hydrophilic particles – predominantly exposed to the carrier phase – experience stronger turbulent stresses and higher detachment rates than lipophilic particles, which are partially shielded on the droplet side. This asymmetry leads to a decreased likelihood of detachment, and thus enhances the trapping efficiency for lipophilic particles.
(
$a$
) The PDF of the phase-field values sampled by the particles for the three wettability cases. Vertical dashed lines indicate the corresponding minima of the potential. (
$b$
) The steady-state fractions of particles located at the interface
$\langle N_i/N_0 \rangle$
, and particles located in the carrier phase
$\langle N_c/N_0 \rangle$
, where
$N_0$
is the total number of particles. Colours identify hydrophilic (blue), neutral (grey) and lipophilic particles (red). Error bars represent the standard deviation.
Figure 3. Long description
The image contains two graphs that analyze the behavior of particles in oil-water emulsions based on their wettability. Panel A: A histogram displays the probability density function (PDF) of the phase-field values sampled by the particles for three wettability cases: hydrophilic, neutral, and lipophilic. The x-axis represents the phase-field values (phi_p) ranging from -1 to 1, and the y-axis represents the PDF on a logarithmic scale. Vertical dashed lines indicate the corresponding minima of the potential. The hydrophilic case is shown in blue, the neutral case in grey, and the lipophilic case in red. Panel B: A line graph shows the steady-state fractions of particles located at the interface (N_i/N_0) and particles located in the carrier phase (N_c/N_0), where N_0 is the total number of particles. The x-axis represents the wettability parameter (k_H+) ranging from -2000 to 2000, and the y-axis represents the fractions. The hydrophilic region is marked in blue, the neutral region in grey, and the lipophilic region in red. Error bars represent the standard deviation.
4.1. Kinetic description of trapping efficiency
To assess the validity of the previous observations on the effect of wettability on particle trapping, we propose a kinetic model that describes particle scavenging as a balance between attachment to and detachment from the interface. In particular, we model these events as stochastic jump processes governed by first-order kinetics. We start by classifying particles depending on their state: particles in the carrier phase (state
$c$
), particles in the droplet phase (state
$d$
) or trapped at the interface (state
$i$
). The particle state is determined by the phase-field value sampled by the particle. Particles are classified as trapped at the interface when
$|\phi _p|\le \phi _{\textit{lim}}$
, while values
$\phi _p\lt -\phi _{\textit{lim}}$
and
$\phi _p\gt \phi _{\textit{lim}}$
identify particles in the carrier and droplet bulks, respectively. The transfer of particles between the three states occurs through attachment events from the phase bulks to the interface, and detachment events from the interface back into either phase bulk. Attachment fluxes per unit area are modelled as proportional to the bulk particle concentration (Hajisharifi et al. Reference Hajisharifi, Marchioli and Soldati2021),
$J_{c\rightarrow i}=\alpha _{c\rightarrow i}{N_c}/{V_c}$
and
$J_{d\rightarrow i}=\alpha _{d\rightarrow i}{N_d}/{V_d}$
, where
$V_c$
and
$V_d$
are the carrier and droplet volumes, and
$\alpha _{c\rightarrow i}$
and
$\alpha _{d\rightarrow i}$
are attachment rate coefficients reflecting near-interface transport. Similarly, detachment fluxes are modelled as the result of first-order kinetic processes,
$J_{i\rightarrow c}=\beta _{i\rightarrow c}{N_i}/{A}$
and
$J_{i\rightarrow d}=\beta _{i\rightarrow d}{N_i}/{A}$
, where
$A$
is the total interfacial area, and
$\beta _{i\rightarrow c}$
and
$\beta _{i\rightarrow d}$
are detachment rate coefficients associated with the competition between local turbulent stresses and detachment energy barriers. The resulting evolution equations are given by
\begin{align} \left \{ \begin{aligned} \frac {\textrm{d}N_c}{\textrm{d}t} &= -\frac {A}{V_c}\alpha _{c\rightarrow i} N_c + \beta _{i\rightarrow c} N_i, \\[3pt]\frac {\textrm{d}N_i}{\textrm{d}t} &= \frac {A}{V_c}\alpha _{c\rightarrow i} N_c + \frac {A}{V_d}\alpha _{d\rightarrow i} N_d - (\beta _{i\rightarrow c}+\beta _{i\rightarrow d}) N_i, \\[3pt]\frac {\textrm{d}N_d}{\textrm{d}t} &= -\frac {A}{V_d}\alpha _{d\rightarrow i} N_d + \beta _{i\rightarrow d} N_i, \end{aligned} \right . \end{align}
where
$N_i(t)$
denotes the number of particles residing within the interfacial layer at time
$t$
, and
$N_{c}(t)$
and
$N_{d}(t)$
denote the numbers of particles in the carrier bulk and in the droplet bulk, respectively. The population of particles is conserved, implying constant
$N_0=N_i(t)+N_{c}(t)+N_{d}(t)$
. The transfer rate coefficients are obtained directly from the DNS by tracking particle state changes between consecutive saved times, estimating the corresponding transition probabilities among the carrier bulk, interface and droplet bulk states, and converting these probabilities into continuous-time rates assuming first-order memoryless kinetics. Additional details are provided in Appendix A.
(a) The fraction of particles at the interface, and (b) the ones in the carrier fluid over time. Each plot shows the time evolution computed from the DNS (solid markers) and from the first-order kinetic model (lines). Results are shown for the three different wettability cases. The insets show the transfer rate coefficients obtained from transition probabilities across the simulations. The inset in (a) shows the deposition rates
$\alpha _{c\rightarrow i}$
and
$\alpha _{d\rightarrow i}$
, while the inset in (b) reports the detachment rates
$\beta _{i\rightarrow c}$
and
$\beta _{i\rightarrow d}$
.
Figure 4(
$a$
) shows the temporal evolution of the fraction of particles trapped at the interface,
$N_i/N_0$
, while figure 4(
$b$
) reports the fraction of particles in the carrier phase,
$N_c/N_0$
. The DNS results are shown using solid markers, while predictions of the first-order kinetic model are represented by solid lines. First, we observe that the system reaches a statistically steady state at
$t^+ \simeq 6000$
, indicating that a balance between attachment and detachment processes has been attained. A comparison between the model predictions and DNS data shows that the model closely reproduces the temporal evolution of both the fraction of particles trapped at the interface and those dispersed in the carrier phase. Owing to mass conservation, this implies that the fraction of particles present in the droplets is also well predicted. The good agreement shows that the dominant particle-exchange dynamics can be reduced to a three-state kinetic process governed by effective constant rates, with attachment and detachment events induced by turbulent fluctuations well represented, at the population level, by memoryless Poisson transition processes.
The transfer rate coefficients shown in the insets of figure 4 provide direct insight into the role of wettability in controlling trapping efficiency. The attachment rate coefficient (inset in figure 4
a) remains unchanged across the different wettability cases, since the underlying flow field – and therefore the turbulent transport of particles toward the interface – is the same in all simulations. In contrast, the detachment rate coefficient from the interface to the carrier phase (inset in figure 4(b), filled symbols) is significantly lower for lipophilic particles (
$k_{\kern-1pt H}^+\gt 0$
) than for hydrophilic ones. The detachment rate coefficient from the interface to the droplet phase (inset in figure 4(b), empty symbols) exhibits the opposite trend; however, its magnitude is orders of magnitude smaller than the other coefficients, therefore its contribution is negligible. Overall, the behaviour of the attachment and detachment rates confirms that the enhanced trapping efficiency of lipophilic particles stems from reduced detachment rates, leading to more stable interfacial trapping.
The observed trends in detachment rates can be attributed to differences in turbulent stresses experienced at the interface, which depend on the wettability-controlled equilibrium position of particles. Lipophilic particles, concentrated on the internal side of droplet surfaces, are exposed to weaker fluctuations than hydrophilic ones concentrated on the external side, consistent with the reduction of turbulent kinetic energy across the interface (Scarbolo & Soldati Reference Scarbolo and Soldati2013). This interpretation is further supported by the PDFs of the turbulent kinetic energy
$k^+$
sampled by interfacial particles (figure 5). Lipophilic particles (red) exhibit a less populated high-
$k^+$
tail and a higher probability of low-
$k^+$
events (highlighted in the inset), indicating reduced exposure to intense fluctuations. Since detachment occurs when turbulent stresses overcome the interfacial trapping potential, the lower frequency of such events decreases the likelihood of removal from the interface.
The PDFs of the turbulent kinetic energy
$k^+$
sampled by particles located at the interface for the three wettability conditions. The inset presents the same PDFs on a linear scale, with a magnified view of the low-
$k^+$
region. Colours denote hydrophilic (blue), neutral (grey) and lipophilic (red) particles. All quantities are expressed in wall units.
5. Conclusions
We performed interface-resolved direct numerical simulations (DNS) of turbulent oil–water emulsions laden with particles of different wettability (hydrophilic, neutral and lipophilic). A phase-field description of the deformable interface, coupled with a Lagrangian particle model including capillary adhesion and wettability forces, enables the direct resolution and quantification of particle–interface interactions under fully turbulent conditions. The simulations show a clear asymmetry in trapping efficiency: lipophilic particles are retained at the interface more effectively than hydrophilic ones, consistent with their preferential affinity for the droplet phase. This difference originates from their wettability-controlled equilibrium position within the interfacial layer. Particles residing on the droplet side experience weaker turbulent stresses, which reduces detachment rates and increases interfacial residence times. Trapping is therefore governed not only by interfacial energetics, but also by the asymmetric turbulent environment across the interface (Scarbolo & Soldati Reference Scarbolo and Soldati2013; Soligo et al. Reference Soligo, Roccon and Soldati2020). Finally, a minimal first-order kinetic model quantitatively reproduces the DNS results, and identifies reduced detachment rates as the dominant mechanism underlying enhanced trapping. These findings provide a mechanistic framework for wettability-driven particle segregation in turbulent emulsions, with implications for the design and optimisation of flotation, emulsification and demulsification processes.
Acknowledgements
We acknowledge the Vienna Scientific Cluster (VSC) for generous allowance of computational resources (project ID 72788). The authors acknowledge TU Wien for financial support through its Open Access Funding Programme.
Funding
A.R., D.Z. and A.S. gratefully acknowledge financial support from PRIN 2022, ‘The fluid dynamics of interfaces: mesoscale models for bubbles, droplets, and membranes and their coupling to large scale flows’ (2022R9B2MW – G53C24000810001).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
Data are available on request from the authors.
Appendix A. Computation of transfer rates
The transfer coefficients appearing in (4.1) are inferred directly from the statistics of individual particle transitions observed in the point-particle DNS. Specifically, at every saved time instant, each particle is assigned to one of three states according to the value of the phase field at the particle position, as defined in § 4.1: carrier bulk (
$c$
), interface (
$i$
), or droplet bulk (
$d$
). By comparing the state of the same particle at two consecutive outputs separated by the saving interval
$\Delta T^+ = 2000 \Delta t^+$
(with
$\Delta t^+$
the simulation time step), each event is classified as one of the following transitions: attachment to the interface,
$c\to i$
and
$d\to i$
, or detachment from the interface,
$i\to c$
and
$i\to d$
. All transitions are accumulated over all particles and over all sampling intervals within the statistically stationary regime (
$t^+\gt 6000$
), and we denote with
$\mathcal N_{x\to y}$
the total number of particles that transitioned from state
$x$
to state
$y$
, and with
$\mathcal N_x$
the total number of particles occupying state
$x$
at the beginning of the sampling intervals considered. The corresponding probabilities that a transition between two states occurs are estimated through the pooled estimators
\begin{align} p_{c\to i}=\frac {\mathcal N_{c\to i}}{\mathcal N_c}, \quad p_{d\to i}=\frac {\mathcal N_{d\to i}}{\mathcal N_d}, \quad p_{i\to c}=\frac {\mathcal N_{i\to c}}{\mathcal N_i}, \quad p_{i\to d}=\frac {\mathcal N_{i\to d}}{\mathcal N_i}. \end{align}
Thus
$p_{x\to y}$
represents the probability that a particle initially in state
$x$
has transitioned to state
$y$
after a time
$\Delta T^+$
. Assuming first-order memoryless kinetics, each transition process can be modelled as Poissonian, and each transition probability can be related to the corresponding continuous-time transition rate
$\vartheta$
as
\begin{align} \vartheta =-\frac {1}{\Delta T^+}\ln (1-p). \end{align}
Applying this relation to attachment towards the interface gives the per-particle rates
\begin{align} \vartheta _{c\to i}= -\frac {1}{\Delta T^+}\ln (1-p_{c\to i}), \quad \vartheta _{d\to i}= -\frac {1}{\Delta T^+}\ln (1-p_{d\to i}), \end{align}
which have dimensions of inverse time, and represent the attachment frequency of particles initially located in the carrier and droplet bulk phases, respectively. However, in the population-balance model (4.1), attachment is formulated in terms of interfacial fluxes per unit area, consistent with the model proposed by Hajisharifi et al. (Reference Hajisharifi, Marchioli and Soldati2021):
\begin{align} J_{c\to i}=\alpha _{c\to i}\frac {N_c}{V_c}, \quad J_{d\to i}=\alpha _{d\to i}\frac {N_d}{V_d}. \end{align}
Multiplication of these fluxes by the total interfacial area
$A$
gives the total number of attachment events per unit time. Therefore, equating this representation with the formulation based on the rates (A3), inferred from the pooled probabilities, yields
\begin{align} \vartheta _{c\to i}N_c=A J_{c\to i}, \quad \vartheta _{d\to i}N_d=A J_{d\to i}, \end{align}
from which we compute the attachment rates (
$\alpha _{c \to i} $
and
$\alpha _{d \to i}$
) appearing in the present model:
\begin{align} \alpha _{c\to i}=\frac {V_c}{A}\vartheta _{c\to i}, \quad \alpha _{d\to i}=\frac {V_d}{A}\vartheta _{d\to i}. \end{align}
For particles initially residing at the interface, detachment may occur either towards the carrier phase or towards the droplet phase. Since these two outcomes are mutually exclusive over one sampling interval, the probability of leaving the interface during
$\Delta T^+$
is
$p_{i\to c}+p_{i\to d}$
, and the corresponding total detachment rate is therefore
\begin{align} b_i= -\frac {1}{\Delta T^+} \ln \!\left (1-p_{i\to c}-p_{i\to d}\right )\, . \end{align}
In the reduced model, detachment from the interface is described by two first-order transfer processes with coefficients
$\beta _{i\to c}$
and
$\beta _{i\to d}$
, and fluxes
\begin{align} J_{i\to c}=\beta _{i\to c}\frac {N_i}{A}, \quad J_{i\to d}=\beta _{i\to d}\frac {N_i}{A}, \end{align}
so that the total rate of leaving the interfacial state is simply
\begin{align} b_i=\beta _{i\to c}+\beta _{i\to d}. \end{align}
This total exit rate is then partitioned into detachment towards the carrier phase and detachment towards the droplet phase according to the relative frequency of the two observed detachment events:
\begin{align} \beta _{i\to c} = b_i\frac {p_{i\to c}}{p_{i\to c}+p_{i\to d}}, \quad \beta _{i\to d} = b_i\frac {p_{i\to d}}{p_{i\to c}+p_{i\to d}}. \end{align}
Once the four coefficients
$\alpha _{c \rightarrow i}$
,
$\alpha _{d \rightarrow i}$
,
$\beta _{i \rightarrow c}$
and
$\beta _{i \rightarrow d}$
are determined, the population-balance equations can be solved directly. The good agreement between the predicted evolution of
$N_c(t)$
,
$N_i(t)$
and
$N_d(t)$
and the DNS results demonstrates that although the coefficients are inferred from transition probabilities measured directly from the DNS, the essential mechanisms governing particle–interface exchange can be effectively reduced to and accurately captured by a three-state kinetic model governed by constant transfer rates.
Appendix B. Sensitivity of the results to the phase-field parameters
To assess the robustness of the proposed framework, we perform a sensitivity analysis on the influence of the interface thickness. Specifically, all cases presented above are repeated at the same grid resolution with an increased Cahn number, corresponding to a thicker interfacial layer. When the Cahn number is varied, the adhesion parameter
$\mathcal{A}$
must be rescaled proportionally to
$\textit{Ch}^{-2}$
in order to preserve the same particle interaction potential with the interface,
$U^+(\phi _p)$
, expressed as a function of the phase field variable sampled by the particle,
$\phi _p$
. This follows directly from the interaction potential introduced in (2.9),
\begin{align} U^+(\phi _p)=\frac {1}{2}\pi \mathcal{A}\frac {\textit{Re}_{\tau }}{\textit{We}} {\mathcal{D}^+}^2 -\pi k_{\kern-1pt H}^+\phi _p+C\, , \end{align}
and from the expression for the signed distance given in (2.7). In fact, expressing this relation in wall units gives
\begin{align} \mathcal{D}^+ = \frac {\sqrt {2}}{2}\textit{Re}_\tau \textit{Ch}\, \log \left (\frac {1+\phi _p}{1-\phi _p}\right )\! , \end{align}
and substituting into the potential yields
\begin{align} U^+(\phi _p) = \frac {\pi }{4}\mathcal{A}\frac {\textit{Re}_\tau ^3}{\textit{We}}\textit{Ch}^2 \left [ \log \left (\frac {1+\phi _p}{1-\phi _p}\right ) \right ]^2 -\pi k_{\kern-1pt H}^+\phi _p + C\, , \end{align}
which shows that at fixed flow conditions (
$Re_\tau$
,
$\textit{We}$
) and fixed wettability parameter
$k_{\kern-1pt H}^+$
, the potential depends on
$\mathcal{A}$
and
$\textit{Ch}$
only through the product
$\mathcal{A}\,\textit{Ch}^2$
. Consequently, a variation of the Cahn number from
$\textit{Ch}_1$
to
$\textit{Ch}_2$
must be compensated by rescaling the adhesion parameter such that
$\mathcal{A}_1\textit{Ch}_1^2=\mathcal{A}_2\textit{Ch}_2^2$
, in order to preserve the same potential profile as a function of
$\phi _p$
. This invariance implies that both the equilibrium particle position (corresponding to the minimum of the potential) and the trapping energy barriers associated with detachment towards the carrier and droplet phases remain unchanged.
Sensitivity analysis of the results to the Cahn number. (a) The steady-state PDF of the phase-field values
$\phi _p$
sampled by the particles as a function of the wettability (
$k_{\kern-1pt H}^+=-1720$
hydrophilic,
$k_{\kern-1pt H}^+=0$
neutral,
$k_{\kern-1pt H}^+=+1720$
lipophilic). The inner vertical dashed lines mark the equilibrium positions predicted by the minima of the effective potential, while the outer dashed lines delimit the bulk carrier (
$\phi _p\lt -0.964$
) and droplet (
$\phi _p\gt +0.964$
) regions. (b) The corresponding steady-state fractions of particles at the interface
$\langle N_i/N_0\rangle$
(circles), and in the carrier phase
$\langle N_c/N_0\rangle$
(triangles), as functions of
$k_{\kern-1pt H}^+$
; error bars denote the temporal standard deviation in the steady-state window. In both plots, full-colour curves and solid symbols refer to the simulations at
$\textit{Ch}=0.03$
, whereas faded curves and symbols correspond to
$\textit{Ch}=0.02$
.
In the present analysis, the reference value
$\textit{Ch}=0.02$
is increased to
$\textit{Ch}=0.03$
, and application of the above rescaling to the reference adhesion parameter
$\mathcal{A}=0.1$
yields the corresponding value
$\mathcal{A}=0.044$
for the thicker interfacial layer. The results are reported in figure 6, which shows the steady-state PDF of the phase-field values sampled by the particles (figure 6
a), together with the corresponding steady-state fractions of particles at the interface
$\langle N_i/N_0\rangle$
, and in the carrier phase
$\langle N_c/N_0\rangle$
(figure 6
b). The two sets of simulations, performed at
$\textit{Ch}=0.02$
and
$\textit{Ch}=0.03$
, are denoted by faded and full-colour symbols, respectively.
The distributions of
$\phi _p$
are very similar in the two cases, with the peaks associated with the equilibrium positions of hydrophilic, neutral and lipophilic particles occurring at the same values of
$\phi _p$
and comparable widths, consistent with the preservation of the interaction potential in the
$\phi _p$
coordinate. A small but systematic difference is nevertheless observed in the steady-state particle fractions. When increasing the Cahn number from
$\textit{Ch}=0.02$
to
$\textit{Ch}=0.03$
, the fraction of particles trapped at the interface slightly increases, while the fraction of particles remaining in the carrier phase decreases accordingly. The largest deviations are found for the lipophilic and neutral cases. Although the rescaling preserves the depth and location of the potential well as a function of
$\phi _p$
, the physical thickness of the diffuse interface in wall units increases with
$\textit{Ch}$
, and as a result, the attractive interaction is active over a wider region, leading to an increased trapping efficiency. Overall, the differences between the two datasets remain limited, indicating that the model presented above is robust with respect to the choice of Cahn number within the range considered.
kH+
kH+<0
ϕ<0
kH+>0
ϕ>0
kH+=0
ϕ=0
ϕ
D
U+(ϕ)
kH+
Reτ=300
We=3.0
Ch=0.02
A=0.1
kH+
kH+
kH+=0
a
b
ϕ=0
a
b
⟨Ni/N0⟩
⟨Nc/N0⟩
N0
αc→i
αd→i
βi→c
βi→d
k+
k+
ϕp
kH+=−1720
kH+=0
kH+=+1720
ϕp<−0.964
ϕp>+0.964
⟨Ni/N0⟩
⟨Nc/N0⟩
kH+
Ch=0.03
Ch=0.02
Open access
