1. Introduction
Particle-laden flows are prevalent in both natural phenomena and various industrial processes. Despite extensive research efforts, our comprehension of the collective behaviour of particles within a fluid flow remains incomplete. Alongside experimental studies, numerical modelling of particle-laden flows has gained substantial importance over the past few decades, driven by advancements in computational power and the development of more efficient and precise numerical methods.
The array of numerical methods available for simulating particle-laden flows spans from ‘particle-resolved’ methods, such as the immersed-boundary method (IBM), which resolves the flow around each individual particle on a fine computational mesh (e.g. Peskin Reference Peskin1972), to fully interpenetrating Eulerian approaches, in which both the fluid and particle phases are treated as continuous media (e.g. Anderson & Jackson Reference Anderson and Jackson1967). Positioned between these two extremes is the Euler–Lagrange (EL) approach, which treats the fluid phase as a continuum, while resolving the trajectories of individual particles.
Euler–Lagrange point-particle methods are particularly effective for simulating fluid flows laden with up to several millions of particles, providing an accurate, straightforward and cost-efficient solution. In EL approaches, the fluid phase dynamics is solved using a classical Eulerian framework, whereas the positions of the particles, treated as Lagrangian point masses, are evolved based on the computed fluid flow field. The forces acting on the particles, such as drag, are typically estimated using semi-empirical models (e.g. Schiller & Naumann Reference Schiller and Naumann1933; Ergun Reference Ergun1952; Wen & Yu Reference Wen and Yu1966).
Different levels of coupling between the fluid and particulate phases can be considered, each suitable for different particle volume fraction and mass loading regimes. For very dilute particle-laden flows with very low particle volume fractions and/or mass loading, one-way coupling is often assumed. In this scenario, the momentum transfer from the particles to the fluid phase is negligible, meaning the flow is assumed to be unaffected by the presence of the particles. However, when the particle volume fraction exceeds approximately
$10^{-5}$
, the momentum transfer becomes significant and cannot be ignored (Michaelides, Sommerfeld & van Wachem Reference Michaelides, Sommerfeld and van Wachem2022). In such cases, where two-way coupling is applied, the particles influence the flow through source terms in the governing fluid momentum equations.
From a computational perspective, the momentum transfer between the particles and the fluid in the EL framework is commonly addressed using the particle source-in-cell (PSIC) model proposed by Crowe, Sharma & Stock (Reference Crowe, Sharma and Stock1977). This model has become a cornerstone in simulating particle-laden flows due to its ability to incorporate the interactions between discrete particles and the continuous fluid phase effectively and has been extensively utilised over the past decades to model particle-laden flows (e.g. van Wachem et al. Reference van Wachem, Van der Schaaf, Schouten, Krishna and van den Bleek2001; Marchioli et al. Reference Marchioli, Soldati, Kuerten, Arcen, Taniere, Goldensoph, Squires, Cargnelutti and Portela2008; Eaton Reference Eaton2009). However, one of the important assumptions of this approach is that the ratio
$d_{{p}}/\Delta x$
is very small, where
$d_{{p}}$
is the particle diameter, and
$\Delta x$
represents the computational mesh spacing. In our recent work (Evrard, Denner & van Wachem Reference Evrard, Denner and van Wachem2021), we show that the error of the PSIC-EL method in the Stokes regime is linearly proportional to the ratio
$d_{{p}}/\Delta x$
, and can be as large as 10 %, even if following the recommendations of commonly used best-practice guidelines (Sommerfeld, van Wachem & Oliemans Reference Sommerfeld, van Wachem and Oliemans2008), which advises
$d_{{p}}/\Delta x \lt 0.1$
. These errors have contributed to poor results in several detailed validation studies conducted using the PSIC-EL framework, and one of them has concluded ‘Two-way coupled Eulerian–Lagrangian simulations using the point-force technique have not fulfilled their early promise.’ (Eaton Reference Eaton2009).
Recent research efforts have focused on improving the accuracy and extending the PSIC-EL method to handle particulate flow simulations which do not satisfy the condition
$d_{{p}} \ll \Delta x$
, by modifying the momentum transfer using convolution with smooth kernels of a certain scale, and by estimating the drag force based on the undisturbed fluid velocity from the disturbed velocity field available on the Eulerian mesh, along with other flow parameters. Applying a smooth kernel to the momentum transfer between the particles and the fluid flow can mitigate some of the errors associated with two-way coupling between the fluid and particles (e.g. Capecelatro & Desjardins Reference Capecelatro and Desjardins2013; Evrard, Denner & van Wachem Reference Evrard, Denner and van Wachem2019; Poustis et al. Reference Poustis, Senoner, Zuzio and Villedieu2019). Notably, the magnitude of the flow disturbance due to momentum transfer reaches a plateau as the ratio
$d_{{p}}/ \Delta x$
increases, instead of continuing to grow proportionally with this ratio (Evrard, Denner & van Wachem Reference Evrard, Denner and van Wachem2020). The value of this plateau is directly related to the length scale of the regularisation kernel, which spreads the transferred momentum over a broader region, resulting in smaller errors.
For a given regularisation length scale, further error reduction in estimating the drag force necessitates a strategy to improve the estimate of the hydrodynamic forces between the fluid and the particles. Eulerian–Lagrangian frameworks require an estimate of the relative undisturbed velocity between the fluid and the particle, where the undisturbed fluid velocity is the velocity which the fluid would have at the location of the particle, if the particle under consideration were not present, see figure 1. For applying this concept, there are a number of recent research works (e.g. Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2015; Evrard et al. Reference Evrard, Denner and van Wachem2020; Pakseresht & Apte Reference Pakseresht and Apte2021; Balachandar & Liu Reference Balachandar and Liu2022; Horwitz et al. Reference Horwitz, Iaccarino, Eaton and Mani2022; Kim & Balachandar Reference Kim and Balachandar2024; Chandran, Evrard & van Wachem Reference Chandran, Evrard and van Wachem2025; Srinivas & Tomar Reference Srinivas and Tomar2025) which focus on recovering the undisturbed fluid velocity from the actual (disturbed) fluid velocity field available, along with other flow parameters. Determining the correct hydrodynamic forces in a two-way coupled EL framework by estimating the undisturbed fluid velocity at each particle and subsequently using this undisturbed velocity to accurately determine the drag force on the particle using an existing drag model has achieved some success, especially in very dilute flows with low particle Reynolds number. However, there still exist several fundamental problems to overcome before achieving a general solution following this route.
Figure 1. The disturbed flow (a) and undisturbed flow (b) for particle
$\mathcal{P}_i$
under consideration.
Firstly, the concept of the undisturbed fluid velocity, as introduced by Maxey and Riley, and Gatignol (Maxey & Riley Reference Maxey and Riley1983; Gatignol Reference Gatignol1983) may not be the correct quantity to accurately determine the hydrodynamic forces. The undisturbed fluid velocity is the fluid velocity in which the contribution of the particle under consideration is removed and replaced by the fluid, see figure 1. This means that the self-induced flow perturbation by the presence of the particle due to the momentum fed back to the fluid by this particle is not present in the undisturbed flow. For some cases, this undisturbed fluid velocity is easily determined and can be used to compute the hydrodynamic forces on the particle. For instance, when a single particle is falling in a quiescent fluid, the complete motion of the fluid around the particle is due to the disturbance caused by that particle, and the undisturbed fluid velocity should be zero. However, not all possible cases are that easy to analyse. For instance, when a particle in a flow crosses its own trajectory at a later instance in time, the disturbance from the first instance of the presence of the particle should be considered when determining the drag on the particle. The situation gets even more complicated when multiple particles are present. When estimating the undisturbed velocity from a simulation in the absence of the particle, by ‘removing’ the particle under consideration, the secondary effects of this particle, such as the indirect influences occurring on the neighbouring particles, are also removed, which is likely to lead to an incorrect drag force prediction. This is visualised in figure 1, where the fluid velocity field is shown with the particle
$\mathcal{P}_i$
and without it, i.e. the undisturbed velocity field. In the absence of the particle
$\mathcal{P}_i$
, the flow around the neighbouring particles is altered, which has an effect on the force experienced by the particle
$\mathcal{P}_i$
. Typically, models to determine the undisturbed velocity only determine the self-induced velocity disturbance of the particle under consideration, which does not consider the secondary effects on neighbouring particles. Neglecting these secondary effects is only acceptable in the very dilute regime. In the dense regime, force closure models for the estimation of the drag on the particle should account for the particle-induced self-disturbance, the effect this disturbance has on the neighbouring particles, and finally the neighbour-induced disturbance accounting for these secondary effects.
Secondly, accurately determining the undisturbed fluid velocity in practice is generally very computationally expensive because the velocity disturbance is generally a result of non-linear interactions with the background flow. The general governing equation for the disturbance velocity caused by a particle, from which the undisturbed velocity can be obtained, is very similar to the Navier–Stokes equations (Evrard et al. Reference Evrard, Denner and van Wachem2020), with a similar cost to solve. These equations would need to be solved for every particle in the flow, which would be tremendously computationally expensive. Therefore, most approaches to determine the undisturbed fluid velocity assume Stokes or Oseen flow (e.g. Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2015; Ireland & Desjardins Reference Ireland and Desjardins2017; Evrard et al. Reference Evrard, Denner and van Wachem2020; Horwitz et al. Reference Horwitz, Iaccarino, Eaton and Mani2022; Chandran et al. Reference Chandran, Evrard and van Wachem2025), or solve an auxiliary set trying to obtain all the particle-induced velocity perturbations as one, in which particle–particle effects are neglected (e.g. Pakseresht & Apte Reference Pakseresht and Apte2021). However, these are all approximations of the actual disturbance velocity, and a general solution to this problem is likely to be far too computationally expensive in practice.
In the present paper, we propose an alternative framework to accurately obtain the hydrodynamical forces acting on a particle in a flow. This proposed framework does not rely on the undisturbed fluid velocity, but on quantities directly available in an EL framework. A consistent coupling between fluid and particles can be obtained by volume-filtering the Navier–Stokes equations (NSE) (Hausmann et al. Reference Hausmann, Chéron, Evrard and van Wachem2024a ), a concept introduced by Anderson & Jackson (Reference Anderson and Jackson1967). It is important to note that, in EL simulations, the accuracy of the filtered velocity field at the particle location depends on the fidelity of the underlying volume-filtered equations and their associated closures. Next to the hydrodynamic force model, inaccuracies in the subfilter stress closure or the viscous stress closure can also lead to deviations in the recovered filtered velocity, and thus introduce additional error in the predicted hydrodynamic force.
The framework proposed in the present paper links the forces acting on the particles to the local volume-filtered fluid quantities and the ratio of the filter length and the particle diameter, both well defined in the volume-filtered framework and independent of the particle volume fraction or the fluid mesh resolution. We will first illustrate the novel framework by considering a single particle in a Stokes flow, then extend the framework to a single particle in higher-Reynolds-number flows, and finally to an assembly of particles in a flow.
2. Governing equations
We consider the framework of volume-filtering because it allows the derivation of the concept of EL point-particle simulations directly from the NSE. Volume-filtering a flow quantity
$\varPhi$
is defined as the following convolution operation over the whole domain
$\varOmega$
, the union of the fluid domain
$\varOmega _{\!{f}}$
and the particle domain
$\varOmega _{{p}}$
:
\begin{align} \varepsilon _{\!{f}}(\boldsymbol{x}) \overline {\varPhi }(\boldsymbol{x}) = \int \limits _{\varOmega } I_{\!{f}}(\boldsymbol{y}) \varPhi (\boldsymbol{y})g(|\boldsymbol{x}-\boldsymbol{y}|) \mathrm{d}V_y , \end{align}
where
${\rm d} V_y$
indicates an element of volume in the neighborhood of point
$y$
, and
$\overline{\varPhi}$
indicates the filtered quantity, and with the radially symmetrical filter kernel, with a length scale called the filter width,
$\sigma$
, satisfying
\begin{align} \int \limits _{\varOmega }g(|\boldsymbol{x}|) \mathrm{d}V_x = 1, \end{align}
where
${\rm d} V_x$
indicates an element of volume in the neighborhood of
$x$
. Although, in principle, any filter that obeys (2.2) can be used, in this work, we take the filter kernel to be a Gaussian, as it will enable an analytical treatment of some of the closure terms. The fluid indicator function,
$I_{\!{f}}$
, is defined as
\begin{align} I_{\!{f}}(\boldsymbol{x})=\begin{cases} 1 & \text{if } \boldsymbol{x} \in \varOmega _{\!{f}}\\ 0 & \text{else } \end{cases} \end{align}
and
$\varepsilon _{\!{f}}$
is the fluid volume fraction. The volume-filtered NSE of an incompressible flow with constant density
$\rho _{\!{f}}$
and constant dynamic viscosity
$\mu _{\!{f}}$
may be written as (Hausmann et al. Reference Hausmann, Chéron, Evrard and van Wachem2024a
)
\begin{align} \dfrac {\partial \varepsilon _{\!{f}}}{ \partial t} + \dfrac {\partial }{\partial x_i}(\varepsilon _{\!{f}}\bar {u}_i) &= 0, \end{align}
\begin{align} \rho _{\!{f}}\dfrac {\partial \varepsilon _{\!{f}}\bar {u}_i}{\partial t} + \rho _{\!{f}}\dfrac {\partial }{\partial x_j}(\varepsilon _{\!{f}}\bar {u}_i\varepsilon _{\!{f}}\bar {u}_j) &= -\dfrac {\partial \varepsilon _{\!{f}}\bar {p}}{\partial x_i}+\mu _{\!{f}} \dfrac {\partial ^2\varepsilon _{\!{f}}\bar {u}_i}{\partial x_j \partial x_j} - \sum _q s_{q,i}+ \mu _{\!{f}} \mathcal{E}_i -\rho _{\!{f}}\dfrac {\partial }{\partial x_j}\tau _{{\textit{sfs}},{\textit{ij}}}, \end{align}
where
$u_i$
is the fluid velocity,
$p$
is the fluid pressure,
$\mathcal{E}_i$
represents the viscous closure, and
$\tau _{{\textit{sfs}},{\textit{ij}}}$
the subfilter stress tensor. Here,
$\mathcal{E}_i$
can be expressed analytically, and
$\tau _{{\textit{sfs}},{\textit{ij}}}$
requires modelling and is important to take into account for larger Reynolds numbers (Hausmann et al. Reference Hausmann, Chéron, Evrard and van Wachem2024a
). These closures have been derived in Hausmann et al. (Reference Hausmann, Chéron, Evrard and van Wachem2024a
), where their origin and potential impact are also discussed. The particle momentum source is defined as the following sum of integrals over the surfaces of the particles with index
$q$
:
\begin{align} s_i = \sum _q \int \limits _{\partial \varOmega _{{p},q}}g(|\boldsymbol{x}-\boldsymbol{y}|)\left (-p \delta _{ij} + \mu _{\!{f}} \left (\dfrac {\partial u_i}{\partial y_j}+\dfrac {\partial u_j}{\partial y_i}\right )\right )n_j\mathrm{d}A_y, \end{align}
where
$n_j$
is the normal vector at the surface of the particle. Since
$s_i$
depends on the unfiltered fluid velocity and pressure fields, it typically requires modelling. For large filter widths,
$\sigma$
, to particle diameter,
$d_{{p}}$
, ratios, as it is common for point-particle simulations, the particle momentum source is typically approximated as (see e.g. Capecelatro & Desjardins Reference Capecelatro and Desjardins2013)
\begin{align} s_i \approx \sum _q F_{{h},q,i} \, g(|\boldsymbol{x} - \boldsymbol{x}_{{p},q}|), \end{align}
which has been shown in Hausmann et al. (Reference Hausmann, Chéron, Evrard and van Wachem2024a
) to be a reasonable assumption for
$\sigma /d_{{p}}\geqslant 1$
. Here,
$F_{{h},i}$
is the hydrodynamic force on the particle. In flows with large particle-to-fluid density ratios, the fluid adjusts rapidly to the particle motion, allowing the hydrodynamic force to be well approximated by a quasi-steady drag formulation. In this work, we assume that the dominant hydrodynamic force on the particle is the drag force,
$\boldsymbol{F}_{{h}} \approx \boldsymbol{F}_{{d}}$
, although the framework can be extended to account for other hydrodynamic forces as well. The drag model is a crucial aspect when predicting the motion of particles in a fluid flow. It encompasses a wide range of complexities, from the basic linear drag law for an isolated particle in Stokes flow, to more sophisticated formulations that consider various flow regimes and particle interactions.
3. A novel class of hydrodynamic force correlations based on the volume-filtered fluid velocity
In the hydrodynamic force models used in EL frameworks, the free-stream fluid velocity,
$\boldsymbol{u}_{\infty }$
, is typically interpreted as the undisturbed fluid velocity. Additionally, in some cases, self-induced velocity disturbance correction models are used to determine
$\boldsymbol{u}_{\infty }$
from the volume-filtered fluid velocity which is available in the EL point-particle simulation. The drawbacks of this procedure have been discussed in § 1. We propose a novel class of hydrodynamic force correlations, in which the hydrodynamic force depends directly on the volume-filtered fluid velocity at the particle position and the filter width as inputs. This requires changing existing force correlations, or even deriving new force correlations.
The principle will be first shown for a single particle in Stokes flow and a single particle in a finite-Reynolds-number flow, and will also be derived for a flow past particles in assemblies with varying volume fractions.
3.1. A sphere in Stokes flow
Stokes drag on a sphere is the simplest form of drag force determination, applicable to small spherical particles moving at low particle Reynolds numbers,
${{\textit Re}_p} \ll 1$
, where inertial effects are negligible compared with viscous forces. The particle Reynolds number is the Reynolds number based on the particle diameter and the relative velocity between the fluid and the particle. The expression for Stokes drag for a uniform flow in an infinite domain is analytically given as
\begin{equation} \boldsymbol{F}_{{d}} = 3 \unicode{x03C0} \mu _{\!{f}} d_{{p}} \boldsymbol{U}_{\!{\textit rel}}, \end{equation}
where
$\boldsymbol{U}_{\! \textit{rel}} =(\boldsymbol{u}_{\infty }-\boldsymbol{v})$
is the relative velocity between the particle and the uniform fluid velocity very far away from the particle,
$\boldsymbol{u}_{\infty }$
, herein referred to as the undisturbed velocity. The particle velocity is written as
$\boldsymbol{v}$
. This linear relationship between the drag force and undisturbed velocity indicates that the drag force is directly proportional to the particle velocity. Stokes drag is derived under the assumption of steady, laminar flow with no significant effects from the particle wake.
In order to obtain the drag force dependence on the volume-filtered fluid velocity at the particle centre,
$\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p}$
, and the non-dimensional relative filter width,
$\sigma ^\prime = \sigma /d_{{p}}$
, the analytical velocity,
$\boldsymbol{u}$
, of the Stokes flow past a sphere moving with a velocity
$\boldsymbol{v}$
is volume-filtered. This can be written as follows:
\begin{equation} \varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p} = 2 \unicode{x03C0} \int _0^{\unicode{x03C0} } \int _{d_{{p}}/2}^{\infty } (\boldsymbol{u} + \boldsymbol{v}) \ g \ r^2 \sin \theta \ \mathrm{d}r \ \mathrm{d}\theta , \end{equation}
where
$g$
is the Gaussian distribution function with standard deviation
$\sigma$
. It should be noted that
$\boldsymbol{u}$
is defined in the reference frame of the moving particle in spherical coordinates (Batchelor Reference Batchelor1967). Utilising the symmetry of
$\boldsymbol{u}$
in the azimuthal direction, and integrating
$\boldsymbol{u}$
along the polar (
$\theta$
) and radial (
$r$
) directions, we obtain
\begin{equation} \varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p} = (\boldsymbol{u}_{\infty } - \boldsymbol{v}) {\textit{erfc}}{ \left ( \frac {1}{2 \sqrt {2} \sigma ^\prime } \right ) } + \varepsilon _{\!{f}} \boldsymbol{v}. \end{equation}
Substituting (
$\boldsymbol{u}_{\infty } - \boldsymbol{v}$
) from (3.3) above into (3.1), the drag force on the particle as a function of the local filtered velocity is as follows:
\begin{equation} \boldsymbol{\tilde {F}}_{{d}} = \dfrac {3 \unicode{x03C0} \mu_{f} d_{{p}} (\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p}-\varepsilon _{\!{f}}\boldsymbol{v})}{{\textit{erfc}}{\left ( \dfrac {1}{2 \sqrt {2} \sigma ^\prime } \right ) } }. \end{equation}
The resulting expression for the drag force shares a similar functional form with the correlation proposed by Ireland & Desjardins (Reference Ireland and Desjardins2017), and could, in principle, also be derived within the framework of undisturbed velocity-based models. However, the conceptual foundation of our approach is fundamentally different. Rather than relying on an estimate of the undisturbed fluid velocity, we formulate the drag law directly in terms of the volume-filtered fluid velocity that is readily available in EL simulations.
In the limit of
$\sigma ^\prime \to \infty$
,
$\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p} \to \boldsymbol{u}_{\infty }$
and (3.4) converges to the standard Stokes drag defined in (3.1). The drag force on the particle,
$\boldsymbol{\tilde {F}}_{{d}}$
, is a function of the volume-filtered velocity at the particle centre, a quantity that can be directly interpolated from the computational grid in a volume-filtered EL simulation, and the relative filter width.
3.2. A sphere in finite-Reynolds-number flow
As the particle Reynolds number increases beyond the Stokes regime, the drag force no longer remains proportional to the relative velocity. To determine the drag force for particles with intermediate particle Reynolds numbers, 1
$\lt$
${\textit{Re}}_{{p}}$
$\lt$
1000, empirical correlations are often used to account for the increased nonlinearity in the drag force. One common empirical formula is the correlation of Schiller & Naumann (Reference Schiller and Naumann1933), which proposes a drag coefficient
$C_{{D}}$
of
\begin{equation} C_{{D}} = \frac {24}{{\textit{Re}_p}} \big(1 + 0.15 \, {\textit{Re}_p}^{0.687}\big), \end{equation}
and the drag force is then given by
\begin{equation} \boldsymbol{F}_{{d}} = \frac {1}{2} C_{{D}} \, \rho _{\!{f}} \, {\frac {\unicode{x03C0} }{4}} d^2_{{p}} \left | \boldsymbol{U}_{\! \textit{rel}} \right | \boldsymbol{U}_{\! \textit{rel}}, \end{equation}
where
$\rho _{\!{f}}$
is the density of the fluid.
Similar to the model proposed for a sphere in Stokes flow, we now seek a relation between
$\boldsymbol{u}_{\infty }$
and the volume-filtered velocity at the centre of the particle,
$\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p}$
, for a uniform flow around a sphere at particle Reynolds numbers larger than zero. In this case, it is impossible to derive an analytical relation between
$\boldsymbol{u}_{\infty }$
and
$\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p}$
, as the fluid velocity field around a particle is not known analytically. In order to derive an empirical relation between
$\boldsymbol{u}_{\infty }$
and
$\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p}$
for the uniform flow around a sphere at finite Reynolds numbers, we carry out particle-resolved direct numerical simulations (PR-DNS) of a sphere at Reynolds numbers in the range
$1\leqslant {\textit{Re}_p}\leqslant 200$
.
Our PR-DNS framework employs a finite-volume approach to solve the incompressible Navier–Stokes equations with second-order accurate spatiotemporal discretisation over a mesh. The flow is driven by a body force in the direction of the primary flow, while no-slip and no-penetration boundary conditions at the particle surface are enforced using a momentum source term computed via the hybrid IBM (Chéron et al. Reference Chéron, Evrard and van Wachem2023). The fluid governing equations are solved numerically on a grid which is refined near the surface of the particle, achieving a resolution of
$d_{{p}}/\Delta x \approx 36$
, where
$\Delta x$
is the cell size of the grid near the particle surface. In figure 2, a visualisation of the flow past an isolated particle is shown for the case of particle Reynolds number
${\textit{Re}_p}=100$
.
Figure 2. Flow past a single particle at
${\textit{Re}_p}=100$
. Top half is coloured by the flow velocity normalised by the free-stream velocity,
$u_\infty$
, along with streamlines. Bottom half is coloured by the pressure normalised by
$\rho _{\!{f}} u_\infty ^2 / 2$
.
The volume-filtered fluid velocity at the centre of the particle is obtained by explicitly volume-filtering the velocity field obtained from the PR-DNS, which is done in two steps. Firstly, the velocity field is multiplied with the fluid indicator function, i.e. the fluid velocities in the mesh cells inside the particle are multiplied with zero. Secondly, the convolution integral for the volume-filtered velocity at the centre of the particle is computed discretely as
\begin{align} \varepsilon _{\!{f}}(\boldsymbol{x}_{{c}}) \overline {\boldsymbol{u}}(\boldsymbol{x}_{{c}}) = \int \limits _{\varOmega } I_{\!{f}}(\boldsymbol{y}) \boldsymbol{u}(\boldsymbol{y})g(|\boldsymbol{x}_{{c}}-\boldsymbol{y}|) \mathrm{d}V_y \approx \sum _{n=1}^{N_{\!{f}}} I_{{f},n} \boldsymbol{u}_n\tilde {g}_n\Delta V_n, \end{align}
where
$ I_{{f},n}$
and
$\boldsymbol{u}_n$
are the values of the fluid indicator function and the velocity in mesh cell
$n$
,
$\boldsymbol{x}_c$
indicates the centre of the particle,
$\Delta V_n$
is the volume of the mesh cell,
$\tilde {g}_n$
is the integral contribution of the Gaussian to the mesh cell, and
$N_{\!{f}}$
is the total number of mesh cells in a simulation.
From explicitly volume-filtering the PR-DNS results, pairs of
$\boldsymbol{u}_{\infty }$
and
$\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p}$
for Reynolds numbers in the range of
${{{Re}}_p}\in [1,10,50,100,200]$
and for the relative filter widths
$\sigma ^\prime \in [0.5,1,2,3,4,5]$
are obtained. Note that for values of relative filter widths of
$\sigma ^\prime \gt 5$
,
$\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p}$
no longer changes significantly compared with its value obtained with a relative filter width of
$5$
.
The relation between
$\boldsymbol{u}_{\infty }$
and
$\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p}$
must satisfy two asymptotic limits: (i) (3.3) must be recovered as
${{\overline {Re}}_p}\rightarrow 0$
and (ii)
$\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p} \rightarrow \boldsymbol{u}_{\infty }$
as
$\sigma ^\prime \rightarrow \infty$
. We propose the following empirical correlation that satisfies these asymptotic limits:
\begin{align} \boldsymbol{U}_{\! \textit{rel}} = \boldsymbol{u}_\infty - \boldsymbol{v} = \dfrac {\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p} - \varepsilon _{\!{f}}\boldsymbol{v}}{{\textit{erfc}}{ \left(\dfrac{1}{2 \sqrt {2} \sigma ^\prime } \right)} }(1+k_\sigma k_{{Re}}), \end{align}
with
\begin{align} k_\sigma = \frac {1}{2}\left (\frac {a_0 (\sigma ^{\prime }-0.5)^{a_1}}{1+a_0(\sigma ^{\prime }-0.5)^{a_1}} -1 \right ) \end{align}
and
\begin{align} k_{{Re}} = \frac {1}{2}\left (1 + {erf}(a_2 \log _{10}({{\overline {Re}}_{{p}}}) - a_3) \right )\!, \end{align}
where
${{\overline {Re}}_{{p}}}=\rho _{\!{f}} |\varepsilon _{\!{f}}\bar {\boldsymbol{u}}_{f@p} - \varepsilon _{\!{f}}\boldsymbol{v}|d_{{p}}/\mu _{\!{f}}$
is the filtered particle Reynolds number. The coefficients of the empirical correlation are given in table 1. After obtaining
$\boldsymbol{U}_{\! \textit{rel}}$
, the drag force on the particle can then be computed with (3.6).
Table 1. Coefficients for the empirical correlation to determine
$\boldsymbol{U}_{\! \textit{rel}}$
for finite
${\textit{Re}_p}$
, which is given in (3.8).
Figure 3 shows the target undisturbed Reynolds number,
${\textit{Re}}_{{p}}^{(\textit{targ})}$
, divided by the predicted undisturbed Reynolds number,
${\textit{Re}}_{{p}}^{(\textit{pred})}$
, computed with the empirical correlation for
$\boldsymbol{U}_{\! \textit{rel}}$
for different Reynolds numbers. The maximum deviation of approximately 3 % occurs at a Reynolds number of one for a relative filter width of three. We assume this deviation to be significantly smaller than the modelling and discretisation errors that typically arise in EL point-particle simulations.
Figure 3. Target undisturbed Reynolds number,
${\textit{Re}}_{{p}}^{(\textit{targ})}$
, divided by the predicted undisturbed Reynolds number,
${\textit{Re}}_{{p}}^{(\textit{pred})}$
, with the empirical correlation for
$\boldsymbol{U}_{\! \textit{rel}}$
as given in (3.8) for different normalised filter widths
$\sigma ^\prime$
.
3.3. Suspension of monodisperse spheres
In particle-laden flows of practical relevance, the individual particles do not experience a uniform flow field in an unbounded domain; instead, the presence of surrounding particles alters the local flow so it becomes non-uniform, thereby also influencing the hydrodynamic forces acting on each particle. There are a number of proposed empirical models (e.g. Wen & Yu Reference Wen and Yu1966; Gidaspow Reference Gidaspow1986; Beetstra, Vanderhoef & Kuipers Reference Beetstra, Vanderhoef and Kuipers2007; Tenneti, Garg & Subramaniam Reference Tenneti, Garg and Subramaniam2011) to predict the drag force on a particle in a suspension. These models are typically based on the spatial mean of the relative velocity between the fluid and the particle, the mean particle volume fraction, and the properties of the particle. For instance, in Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011), PR-DNS of various particle assemblies in periodic domains are carried out, and the corresponding forces on the particles in the assemblies are averaged and are used to propose an empirical expression for the drag force, which depends on the ‘global’ fluid velocity and the ‘global’ volume fraction only, where ‘global’ refers to the average over the periodic domain size in which the simulations are carried out.
However, because the flow around each particle, particularly when surrounded by other particles, can be highly complex and vary significantly from one particle to another, the drag force can also exhibit strong particle-to-particle variability. As a result, conventional drag correlations based solely on global parameters cannot provide accurate predictions at the individual particle level. More recently, models have been proposed that aim to predict the deviation of the hydrodynamic force acting on individual particles (e.g. Akiki, Jackson & Balachandar Reference Akiki, Jackson and Balachandar2017; Hardy et al. Reference Hardy, Simonin, De Wilde and Winckelmans2022; Siddani et al. Reference Siddani, Balachandar, Zhou and Subramaniam2024; van Wachem, Elmestikawy & Chéron Reference van Wachem, Elmestikawy and Chéron2024). These models incorporate, in some form, information about the relative position of neighbouring particles. When such local structural information is accounted for, improved accuracy in drag force prediction is possible. However, the definitions of the superficial velocity and global particle volume fraction, on which these models also rely, become ambiguous in inhomogeneous flows, such as flows with particle clustering, since the value of both quantities depends on how large the averaging volume is, which is an arbitrary choice. Additionally, in particle assemblies, particularly at higher concentrations, the notion of an undisturbed velocity becomes increasingly ambiguous due to the complex and overlapping flow disturbances generated by neighbouring particles. As a result, accurately determining this velocity at the particle location becomes highly challenging, both conceptually and computationally.
The presently proposed framework, in which the force correlations depend on volume-filtered quantities, can be extended to particle suspensions. Since the flow is more complicated than the uniform flow around a single particle, additional uniquely defined volume-filtered quantities that are accessible in EL simulations have to be considered to predict the hydrodynamic force accurately. Independent of the complexity of the particle-laden flow, the hydrodynamic force on particle
$q$
is expressed as a function of the pressure and velocity field as
\begin{align} F_{{h},q,i} = \int \limits _{\partial \varOmega _{{p},q}} \left [-p \delta _{ij} + \mu _{\!{f}} \left (\dfrac {\partial u_i}{\partial x_j}+\dfrac {\partial u_j}{\partial x_i}\right )\right ] n_j\mathrm{d}A_x, \end{align}
where
$\partial \Omega_{p,q}$
is the surface of particle
$q$
,
$\delta_{ij}$
is the Kronecker delta function, and
$d A_x$
indicates an element of the surface in the neighborhood of point
$x$
on the particle. If the values of
$u_i$
and
$p$
are known on all points on the surface of the particle. Formally, the process of volume-filtering does not remove any information, but the discretisation of the filtered solution on a finite fluid mesh does. Therefore, an exact relation between the volume-filtered flow quantities and the hydrodynamic force exists, but since the volume-filtered flow quantities are approximated in EL simulations, the predicted hydrodynamic force is also an approximation. Since less information is removed for small filter widths, it is expected that the estimation of the hydrodynamic force is also more accurate for small filter widths.
As a conceptual proof of the proposed framework for force correlations, the functional approach of the existing mean-drag-force correlation of Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011) is adapted to the framework proposed in this paper. The adaption consists of two essential modifications. (i) The input parameters are the volume-filtered velocity and the volume fraction at the particle position as well as the filter width, instead of the superficial velocity and the global volume fraction, as in the original model. (ii) The coefficients of the expression are obtained by minimising the deviation of the predicted force from the actual force of each individual particle, instead of the deviation of the predicted force from the mean of the forces acting on all the considered particles in each PR-DNS of a particle assembly.
To determine the coefficients of the proposed expression, we have carried out PR-DNS of assemblies of particles. The PR-DNS are performed in a periodic cubic domain of volume
$V_{\textit{tot}}$
, containing
$N_{{p}}$
non-overlapping, fixed, monodisperse spherical particles. The domain is discretised into
$N_{\!{f}}$
Eulerian fluid cells, in which the governing equations for fluid flow are solved to obtain the local fluid properties. Each fluid cell may be fully occupied by fluid, fully occupied by a particle or partially occupied by both. Flow through the periodic domain is driven by imposing a pressure drop,
$\delta P$
, which is introduced as an additional body force in the
$x$
-direction in the fluid momentum equations, so the correct superficial fluid velocity is obtained. For further details, we refer to van Wachem et al. (Reference van Wachem, Elmestikawy and Chéron2024).
In this study, a total of six global particle volume fractions are considered, namely
$\langle \varepsilon _{{p}}\rangle$
=
$0.1$
,
$0.2$
,
$0.3$
,
$0.4$
,
$0.5$
and
$0.6$
. For each global particle volume fraction, multiple flow conditions are investigated by varying the superficial particle Reynolds number as
$\textit{Re}_{{s}} = 0.1$
,
$1$
,
$10$
,
$50$
,
$100$
and
$300$
. The superficial particle Reynolds number is defined as
$\textit{Re}_{{s}}=\rho _{\!{f}}d_{{p}} u_{{s}}/\mu _{\!{f}}$
, where the superficial velocity is given as the fluid domain average of the fluid velocity:
\begin{align} u_{{s}} = \dfrac {1}{V_{\textit{tot}}} \int \limits _{\varOmega _{\!{f}}} u_x(\boldsymbol{x})\mathrm{d}V_x. \end{align}
The total volume of the domain,
$V_{\textit{tot}}$
, is the sum of the volume occupied by the fluid and the volume occupied by the particles. To ensure statistical convergence, three independent realisations are simulated for each combination of solid volume fraction and superficial particle Reynolds number. Each simulation is carried out until a statistically steady state is achieved. In total, 108 PR-DNS are performed, with an average of 136.33 particles per simulation. This results in a total of 14 724 data points across the entire parameter space, which serve as the basis for the development of the hydrodynamic force correlations. In figure 4, a two-dimensional section of the normalised streamwise velocity field is shown for
$\langle \varepsilon _{{p}} \rangle =0.2$
and
$0.5$
, both at
${\textit{Re}_s}=50$
.
Figure 4. Examples of (
$u_x$
) stream-wise velocity fields normalised by the superficial velocity as predicted by PR-DNS, for two different volume fractions.
From the results of the PR-DNS, we determine the volume-filtered flow quantities for various filter widths, for each simulation configuration. To obtain the volume-filtered flow quantities, we exploit the fact that the domain is periodic and that a convolution becomes a multiplication in spectral space. Therefore, the volume-filtered velocity is given as
\begin{align} \varepsilon _{\!{f}}\overline {\boldsymbol{u}} = g\ast (I_{\!{f}}\boldsymbol{u}) = \mathcal{F}^{-1}\{ \mathcal{F}\{ g\}\mathcal{F}\{I_{\!{f}}\boldsymbol{u}\} \}, \end{align}
where
$\ast$
is the convolution operator, and
$\mathcal{F}$
corresponds to the Fourier transform.
In order to derive a drag correlation as a function of the local volume-filtered quantities, we propose a drag correlation of a similar form as that in Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011), except that the input parameters here are the local fluid volume fraction, the volume-filtered velocity, and the ratio of the filter width to the particle diameter,
$\sigma ^{\prime }$
. The general form of the normalised drag correlation is as follows:
\begin{align} \boldsymbol{\tilde {F}}_{{d}} &= 3 \unicode{x03C0} \mu_{f} d_{{p}} (\varepsilon _{\!{f}} \bar {\boldsymbol{u}}_{{f@p}} - \varepsilon _{\!{f}}\boldsymbol{v}) \left [ \frac {C_{\mathrm{D}} \, {{\overline {Re}}_{{p}}}}{24} \frac {1}{(1 - \delta _{\varepsilon })^3} + a_0 \frac {\delta _{\varepsilon }}{(1 - \delta _{\varepsilon })^3} \right . \nonumber \\ &\quad \left . +\, a_1 \frac {\delta _{\varepsilon }^{1/3}}{(1 - \delta _{\varepsilon })^4} + \delta _{\varepsilon }^{a_4} {{\overline {Re}}_{{p}}} \left (a_2 + a_3 \frac {\delta _{\varepsilon }^{a_5}}{(1-\delta _{\varepsilon })^2} \right ) \right ]\!, \end{align}
where
$\delta _{\varepsilon }$
is the difference between the fluid volume fraction of an isolated particle and
$\varepsilon _{\!{f}}$
, both of which are evaluated at the particle centre. Here,
$C_{{D}}$
is the Schiller–Naumann correlation based on
${\textit{Re}_p}$
. It is to be noted that, for the correlation to completely depend on volume-filtered quantities, we convert
${{\overline {Re}}_{{p}}}$
to
${\textit{Re}_p}$
using (3.8). For a Gaussian kernel,
$\delta _{\varepsilon }$
can be mathematically expressed as
\begin{equation} \delta _{\varepsilon } = {\textit{erfc}} \left (\frac {1}{2 \sqrt {2} \sigma ^\prime } \right ) + \frac {\exp \left (-\dfrac {1}{8 {\sigma ^\prime }^2}\right )}{\sigma ^\prime \sqrt {2 \unicode{x03C0} }} - \varepsilon _{\!{f}}. \end{equation}
Equation (3.15) is obtained by using the fluid volume fraction for an isolated particle (Balachandar & Liu Reference Balachandar and Liu2022), and using L’Hôpital’s rule to retrieve the volume fraction at the particle centre. In the limit of a very dilute volume fraction, i.e.
$\delta _{\varepsilon } \to 0$
, and a large filter width, (3.14) reduces to the standard Schiller–Naumann correlation for a single particle. Mean-drag-force models, such as that proposed by Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011), use the ‘global’ superficial fluid velocity, and the global particle volume fraction to predict a mean drag force in an assembly. When
$\sigma ^{\prime } \to \infty$
,
$\delta _{\varepsilon }$
simplifies to the global particle volume fraction,
$\langle \varepsilon _{{p}} \rangle$
, and the volume-filtered velocity at the particle centre simplifies to the superficial fluid velocity,
$u_{{s}}$
. In this limit, (3.14) uses the same inputs as in a mean-drag-force model, such as that proposed by Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011).
To evaluate the coefficients in (3.14), we define an error measure that minimises the deviation of the drag force prediction from the actual drag experienced by each individual particle in the random assembly. The mean relative deviation of the predicted drag force in the streamwise direction acting on particle
$q$
,
$\tilde {F}_{{d},q,x}$
, from the actual drag force acting on particle
$q$
,
$F_{{d},q,x}$
, is therefore defined as
\begin{align} {E}_{{F}} = \dfrac {1}{N_{{p}}} \sum _{q=1}^{N_{{p}}} \bigg | \dfrac {F_{{d},q,x} - \tilde {F}_{{d},q,x}}{F_{{d},q,x}} \bigg | . \end{align}
Table 2 presents the fitted coefficients of (3.14) obtained by minimising the mean deviation in the predicted individual particle force against the corresponding particle force from the PR-DNS data using (3.16). The datasets are grouped based on
$\sigma ^\prime$
, and use the Limited-memory Broyden-Fletcher-Goldfarb-Shanno for Box contraints algorithm for limited memory bounded-constraint optimisation to evaluate the coefficients (Byrd, Peihuang & Nocedal Reference Byrd, Peihuang and Nocedal1996). A maximum mean error of approximately 27 % in predicting the individual particle forces across the six different values of
$\sigma ^\prime$
is observed.
Table 2. Coefficients of the empirical correlation for
$\boldsymbol{\tilde {F}}_{{d}}$
given in (3.14).
4. Results and discussion
4.1. A sphere falling in a fluid in the Stokes regime
We demonstrate the advantages of the newly proposed framework of force correlations by first considering a single falling spherical particle in a large domain filled with an initially quiescent fluid under gravity in the Stokes regime. The particle, which is initially at rest, accelerates in the fluid until it reaches its terminal velocity, which can be determined analytically in this regime. We perform three configurations of two-way coupled EL point-particle simulations of this case by:
-
(i) solving the volume-filtered NSE with the classical Stokes drag;
-
(ii) solving the volume-filtered NSE with the newly proposed filtered Stokes drag, (3.4); and
-
(iii) solving the commonly used PSIC method (Crowe et al. Reference Crowe, Sharma and Stock1977) and neglecting
$\mathcal{E}_i$
and
$\tau _{{\textit{sfs}},{\textit{ij}}}$
.
For each configuration, the domain is fully periodic and has a size of
$L_x\times L_y\times L_z = 100d_{{p}}\times 100d_{{p}}\times 100d_{{p}}$
and various mesh resolutions and, where applicable, filter widths are simulated. The particle Reynolds number based on the terminal velocity is
${\textit{Re}_p}=1.11\times 10^{-4}$
and the density ratio is
$\rho _{{p}}/\rho _{\!{f}}=2000$
, where
$\rho_p$
is the particle density.
Figure 5 shows the particle velocities of the three different simulation configurations for three different resolutions,
$d_{{p}}/\Delta x = 0.25$
,
$d_{{p}}/\Delta x = 1$
and
$d_{{p}}/\Delta x = 2$
, and in the figure the simulation results are compared with the analytical result. The filter widths,
$\sigma$
, are chosen according to the guidelines provided in Hausmann et al. (Reference Hausmann, Chéron, Evrard and van Wachem2024a
), such that
$\sigma ^\prime \geqslant 1$
and
$\sigma /\Delta x \geqslant 1$
are always satisfied. The time is normalised by
$t_n = (\rho _{{p}}-\rho _{\!{f}})d_{{p}}^2/(18 \mu _{\!{f}})$
.
The simulations using the classical Stokes drag overpredict the magnitude of the terminal velocity considerably, although the simulations do reach a stable velocity after some time. The simulations adopting the PSIC framework show an even larger error, and do not even converge for
$d_{{p}}/\Delta x = 2$
. When using the newly proposed drag force correlation based on the volume-filtered velocity and the ratio of the filter width and the particle diameter, the analytical particle velocity is reproduced accurately in all cases.
Figure 5. Relative particle velocities of a single falling sphere over time for the volume-filtered simulations using classical Stokes drag, the volume-filtered simulations with the filtered Stokes drag, and the PSIC simulation frameworks, compared with the analytical solution. The simulations are performed with three different resolutions: (a)
$d_{{p}}/\Delta x = 0.25$
, (b)
$d_{{p}}/\Delta x = 1$
and (c)
$d_{{p}}/\Delta x = 2$
.
$v_{\infty}$
is the terminal velocity of the particle.
4.1.1. Comparison with undisturbed velocity models
A falling sphere in a finite domain is often used to validate models based on the undisturbed fluid velocity principle. The most recent velocity disturbance corrections for transient flow (Evrard et al. Reference Evrard, Chandran, Cortez and van Wachem2025; Chandran et al. Reference Chandran, Evrard and van Wachem2025) can predict the trajectory of a single falling sphere in a finite domain accurately because this is one of the few cases in which the undisturbed velocity is well defined. We investigate a single falling particle under gravity in a small periodic domain, a simple case where the concept of subtracting the velocity disturbance to achieve the undisturbed fluid velocity fails when the particle is confronted with its historic effect on the fluid.
To simulate this case, two-way coupling is used, i.e. the drag force acting on the particle is fed back to the fluid momentum. This results in a constant acceleration of the mean fluid velocity. Since in steady state the drag force on the particle must be equal to the sum of the gravitational force,
$\boldsymbol{F}_{\!{G}}$
, and the buoyancy force,
$\boldsymbol{F}_{\!{B}}$
, the mean flow accelerates with a magnitude of
$|\boldsymbol{F}_{\!{G}}+\boldsymbol{F}_{\!{B}}|/m_{\!{f}}$
, where the mass of the fluid,
$m_{\!{f}} = \rho _{\!{f}}L_xL_yL_z$
. The periodic domain has a size of
$L_x\times L_y\times L_z = 50d_{{p}}\times 50d_{{p}}\times 50d_{{p}}$
. A steady state is achieved only in the sense of the relative velocity between the particle and the fluid, as the mean fluid flow continues to accelerate, just as the particle. However, the particle velocity in the non-moving Eulerian frame of reference continues to increase in time. The results of the EL simulation using the concept of the undisturbed fluid velocity with the EL simulation using the volume-filtered drag model, as well as the theoretically expected trend, is shown in figure 6. It is clearly evident that, in this case, the concept of reconstructing the undisturbed velocity to compute the drag force fails to produce accurate results, whereas the simulation employing the volume-filtered drag force framework yields predictions in good agreement with the reference solution.
One could argue that a single particle which is falling in a periodic domain does not have much practical relevance. However, the failure of the concept of velocity disturbance correction translates directly to configurations with more than one particle, which is frequently studied in the literature (see e.g. Uhlmann & Doychev Reference Uhlmann and Doychev2014; Capecelatro, Desjardins & Fox Reference Capecelatro, Desjardins and Fox2015; Willen & Prosperetti Reference Willen and Prosperetti2019; Hausmann et al. Reference Hausmann, Elmestikawy and van Wachem2024b ; Xia et al. Reference Xia, Lin, Guo and Yu2024). If these cases were simulated with the EL point-particle approach with an accurate velocity disturbance correction, incorrect results would be obtained.
Figure 6. A falling isolated particle in a periodic domain. Simulations of the volume-filtered NSE using the proposed filtered Stokes drag are compared with the theoretical acceleration of the particle and the ideal correction of the flow disturbance of the particle.
4.2. A sphere falling in a fluid at higher Reynolds numbers
In order to test the novel force correlation that we propose for finite
${\textit{Re}_p}$
, the study of the falling particle discussed in the previous section is extended to larger
${\textit{Re}_p}$
. Similarly to the configuration of the falling particle in the Stokes regime, the domain size remains
$L_x\times L_y\times L_z = 100d_{{p}}\times 100d_{{p}}\times 100d_{{p}}$
for this validation. The density ratio is kept at
$\rho _{{p}}/\rho _{\!{f}}=2000$
, and the particle Reynolds number is varied by varying the fluid viscosity. At finite
${\textit{Re}_p}$
, the results cannot be compared with an analytical solution, but for a single particle, a reference solution can be obtained with a one-way coupled simulation. In the one-way coupled simulation, the fluid velocity remains zero at all times and the classical Schiller–Naumann correlation, as given in (3.5), provides an accurate drag force on the particle. In addition to the one-way coupled simulations, we perform two two-way coupled EL point-particle simulations of this case by:
-
(i) solving the volume-filtered NSE with the classical Schiller–Naumann drag correlation; and
-
(ii) solving the volume-filtered NSE with the newly proposed filtered version of the Schiller–Naumann drag correlation using (3.8).
The simulations are performed with three different particle Reynolds numbers,
${\textit{Re}_p}\in [1.11\times 10^{-4},0.9646,38.7]$
and two relative filter widths
$\sigma ^\prime \in [1,4]$
, whereas
$\sigma /\Delta x=1$
for all simulations.
In figure 7, the particle settling velocities are shown for the three
${\textit{Re}_p}$
and for a relative filter width of
$\sigma ^\prime =1$
. The volume-filtered simulations with the newly proposed filtered Schiller–Naumann drag force model predict the particle settling velocities accurately for all of the three
${\textit{Re}_p}$
. With the classical Schiller–Naumann drag, significant deviations from the correct particle settling velocities are observed, whereby the deviations decrease with increasing particle Reynolds number. For large particle Reynolds numbers, the flow disturbance induced by the particle becomes small, and the volume-filtered velocity at the particle position is relatively close to the undisturbed velocity.
Figure 7. Relative particle velocities of a single falling sphere over time for the volume-filtered simulations using classical Schiller–Naumann drag and the volume-filtered simulations with the filtered Schiller–Naumann drag. The corresponding one-way coupled simulation with classical Schiller–Naumann drag is shown as reference. The simulations are performed with three different values for
${\textit{Re}_p}$
and with a relative filter width of
$\sigma ^\prime =1$
.
Figure 8. Particle velocities of a single falling sphere over time for the volume-filtered simulations using classical Schiller–Naumann drag and the volume-filtered simulations with the filtered Schiller–Naumann drag. The corresponding one-way coupled simulation with classical Schiller–Naumann drag is shown as reference. The simulations are performed for three different values of
${\textit{Re}_p}$
and with a relative filter width of
$\sigma ^\prime =4$
.
Figure 8 shows the particle settling velocities for a larger filter width,
$\sigma ^\prime =4$
. Ideally, the predicted particle settling velocities should not depend on the filter width or the spatial resolution and the same results should be obtained as with
$\sigma ^\prime =1$
. This is the case for the volume-filtered simulations using the filtered Schiller–Naumann drag. With the classical Schiller–Naumann drag, the particle settling velocity is still inaccurate for the smaller
${\textit{Re}_p}$
, but for the larger
${\textit{Re}_p}$
the particle settling velocities are predicted accurately.
This means that, when using the classical Schiller–Naumann drag with the volume-filtering framework, the results deteriorate as the relative spatial resolution increases or the local
${\textit{Re}_p}$
decreases. For large filter widths, the flow disturbance by the particle is spread over a wider region, which leads to a volume-filtered velocity at the particle position closer to the undisturbed velocity than for small filter widths. Therefore, the disturbance is negligible for the largest
${\textit{Re}_p}$
and
$\sigma ^\prime \geqslant 4$
.
Although accurate particle settling velocities are predicted with the volume-filtering framework and the newly proposed filtered Schiller–Naumann drag, some small deviations to the one-way coupled simulations remain, which have multiple origins. (i) The volume-filtered NSE are solved on a discrete fluid mesh, which leads to a discretisation error. Furthermore, the volume fraction, the mass source and the closures, including the analytical viscous closure, are also represented on a discrete fluid mesh. The resulting discretisation errors lead to a deviation of the numerical solution from the actual volume-filtered flow field, which has been used to fit the coefficients of the force correlation. (ii) The model for the subfilter stress tensor is accurate for small filter widths but the modelling error increases as the filter width increases (Hausmann et al. Reference Hausmann, Chéron, Evrard and van Wachem2024a
), and (iii) a falling particle is a transient process but the filtered Schiller–Naumann correlation is obtained with data from the stationary flow around a sphere. In the acceleration phase, the flow field around the particle can be different from the steady-state flow field at the same
${\textit{Re}_p}$
. Unlike undisturbed velocity models, which attempt to remove transient effects only to reintroduce them through corrective terms such as the history force, our approach directly models the total hydrodynamic force using the filtered flow quantities. This avoids ambiguities in defining the undisturbed velocity, and offers a consistent, potentially unsteady-aware framework. The extent to which transient effects are retained remains an open question for future research.
These and other error sources are common in EL point-particle simulations, but the present results indicate that they are small with the volume-filtering framework, at least for the falling sphere configurations investigated.
4.3. Suspension of monodisperse spheres
As a final test case, the novel force correlation framework is evaluated for a suspension of monodisperse spheres. The hydrodynamical force on each individual particle predicted by the force correlation derived in § 3.3 is compared with the actual force given in the PR-DNS of the flow through the particle arrangements. In figures 9 and 10, the mean relative force error,
$E_{{F}}$
, of the correlation given in (3.14) is shown for different filter widths, different global particle volume fractions,
$\langle \varepsilon _{{p}} \rangle$
, and different superficial particle Reynolds numbers,
$\textit{Re}_{{s}}$
. It should be noted that
$\textit{Re}_{{s}}$
and
$\langle \varepsilon _{{p}} \rangle$
are not the filtered quantities at the particle positions but the averages over the entire domain of the respective case, which corresponds to filtered quantities with an infinite filter width. Therefore, the values on the
$x$
-axes are not the direct values used in the correlation, but they characterise the simulation case, which is why some values on the
$x$
-axis are the same, although the filter width is different. This means that the mean relative error
$E_{{F}}$
is the error observed for a specific superficial Reynolds number and a global particle volume fraction, which is somewhat arbitrary. It may very well be that particles belonging to cases with different superficial velocities and different global particle volume fractions have similar volume-filtered velocities or local volume fractions at the particle positions.
Figure 9. Mean relative force error as a function of the superficial Reynolds number for different filter widths and different global particle volume fractions compared with the correlation of Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011).
Figure 10. Mean relative force error as a function of the global particle volume fraction for different filter widths and different superficial Reynolds numbers compared with the correlation of Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011).
The relative errors observed in figures 9 and 10 are of the order of
$20\,\%$
. The proposed correlation cannot be expected to be much more accurate because the volume-filtered velocity and volume fraction are not sufficient to predict the force in such a complex flow with high accuracy. Furthermore, the functional approach used is designed to accurately predict the mean force and not the force on each individual particle. However, more relevant for the proposed framework for force correlations than the magnitude of the error is how the error changes with the filter width.
With the filter width
$\sigma ^\prime =5$
, the volume-filtered velocity and volume fraction at the particle positions are almost identical to the superficial velocity and the global particle volume fraction for all particles. Therefore, the correlations of all particles in the same simulation case have almost identical input values, which leads to predicted forces by the correlation that are also almost identical for all particles. For such large filter widths and if only the volume-filtered velocity and the volume fraction are considered as input parameters, the proposed framework yields a correlation that is similar to existing mean force correlations. However, the proposed framework is also suitable for smaller filter widths. As observed in figures 9 and 10, the mean relative force error decreases for almost every
$\textit{Re}_{{s}}$
and
$\langle \varepsilon_p\rangle$
as the filter width decreases to
$\sigma ^\prime = 0.5$
. For the smaller filter widths, the volume-filtered velocity and volume fraction at the different particle positions varies and the correlation can predict the different drag forces on the particles. Figures 9 and 10 also show the predictions of the correlation of Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011) evaluated at the same filtered Reynolds numbers and volume fractions. For the majority of
${\textit{Re}_s}$
and
$\langle \varepsilon_p\rangle$
, the newly proposed correlation based on the filtered flow quantities shows a smaller error, especially for relative filter widths less than five. To improve the overall accuracy of the force correlation, additional volume-filtered flow quantities need to be included in the correlation, such as the volume-filtered velocity gradient or the volume fraction gradient, which becomes particularly evident for large global particle volume fractions.
As already discussed for the configurations with isolated particles, the proposed framework of drag force correlations does not require the determination of the undisturbed velocity, but is fully based on filtered flow quantities. Moreover, the accuracy of the force prediction does not deteriorate as the filter width or mesh spacing is decreased, as with existing drag force prediction frameworks, but the accuracy of the drag force prediction is improved.
5. Conclusions
This study presents a novel framework for modelling the hydrodynamic forces in EL point-particle simulations that does not rely on the classical concept of the undisturbed fluid velocity. By leveraging the volume-filtered NSE, the proposed approach directly relates the drag force on each of the particles to volume-filtered fluid quantities which are readily available in EL simulations.
The accuracy of the framework has been demonstrated through validation against analytical solutions, one-way coupled simulations, and PR-DNS across a wide range of particle Reynolds numbers and filter widths. In particular, the filtered drag force correlations accurately recover the correct terminal velocity for a falling sphere in both the Stokes and finite-Reynolds-number regimes.
Moreover, the framework has been extended to particle assemblies, where a generalised drag correlation has been derived based on the volume-filtered velocity and the volume fraction. The new correlation shows an improved predictive performance over existing mean-force models, especially at smaller filter widths and mesh resolutions, and at lower volume fractions. To increase the accuracy of the newly proposed framework, apart from the volume-filtered fluid velocity, additional volume-filtered quantities, such as the volume-filtered fluid velocity gradients, are required.
Unlike traditional approaches that depend on estimating undisturbed velocities, the proposed methodology remains robust in both dilute and moderately dense regimes, avoids high computational cost, and enables the formulation of force models that are inherently consistent with the volume-filtered flow description. The proposed framework provides a theoretically sound and computationally efficient basis for the development of next-generation force correlations in EL point-particle simulations of multiphase flows.
Funding
This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation): project-ID 457515061, project-ID 457509672, project-ID 422037413 - TRR 287 and project-ID 466092867 - SPP 2331.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are reproducible and are openly available in the repository https://doi.org/10.5281/zenodo.15364320.
Open access
