2.1. Governing equations

We consider the motion of a single or multiple rigid and non-rotating particles immersed in an incompressible Newtonian fluid flow with constant dynamic viscosity $\mu _{{f}}$ and density $\rho _{{f}}$. Outside of the particles, the fluid motion is governed by the NSE

(2.1)$$\begin{gather} \dfrac{\partial u_i}{\partial x_i} = 0, \end{gather}$$

(2.2)$$\begin{gather}\rho_{{f}}\dfrac{\partial u_i}{\partial t} + \rho_{{f}} \dfrac{\partial u_i u_j}{\partial x_j} ={-} \dfrac{\partial p}{\partial x_i} + \mu_{{f}} \dfrac{\partial}{\partial x_j}\left(\dfrac{\partial u_i}{\partial x_j}+\dfrac{\partial u_j}{\partial x_i}\right), \end{gather}$$

where $u_i$ is the fluid velocity vector and $p$ the pressure. At the boundary of the particle with the index $q$, $\partial \varOmega _{{p},q}$, the no-slip boundary condition is assumed:

(2.3)\begin{equation} u_i|_{\partial \varOmega_{{p},q}} = v_{q,i}, \end{equation}

where $v_{q,i}$ is the velocity of the particle with the index $q$. The particle motion is governed by Newton's second law:

(2.4)\begin{equation} \rho_{{p},q} V_{{p},q} \dfrac{\mathrm{d}v_{q,i}}{\mathrm{d} t}=F_{{h},i} , \end{equation}

where $F_{{h},i}$ is the total hydrodynamic force acting on the particle resulting from the pressure and viscous stresses at the particle surface, $V_{{p},q}$ is the volume of the particle, and $\rho _{{p},q}$ is the density of the particle. Note that gravity and other body force are not considered, although these do not change the analysis from a fundamental point of view. The hydrodynamic force is obtained by integrating the fluid stresses over the particle surface:

(2.5)\begin{equation} F_{{h},i} = \int_{\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)\hat{n}_j\,\mathrm{d}A, \end{equation}

where $\delta _{ij}$ is the Kronecker delta. We assume that the particle surface normal vector $\hat {n}_j$ points towards the outside of the particle.

From a computational point of view, it is desirable to consistently smooth the high gradients arising near the particle to allow for a coarser resolution of the flow. The smoothing can be realized by volume filtering the equations governing the flow using a filter kernel $g$, as originally proposed by Anderson & Jackson (Reference Anderson and Jackson1967) and later extended to the Euler–Lagrange framework by Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013). An arbitrary quantity $\varPhi$ that is defined in the fluid domain $\varOmega _{{f}}$ can be volume filtered according to

(2.6)\begin{equation} \epsilon_{{f}}(\boldsymbol{x}) \bar{\varPhi}(\boldsymbol{x}) = \int_{\varOmega} I_{{f}}(\boldsymbol{y}) \varPhi(\boldsymbol{y})g(|\boldsymbol{x}-\boldsymbol{y}|)\,\mathrm{d}V_y , \end{equation}

where $\bar {\varPhi }(\boldsymbol {x})$ is the filtered quantity, $\varOmega =\varOmega _{{f}}\cup \varOmega _{{p}}$ is the union of the fluid and the particle domain, and $\epsilon _{{f}}(\boldsymbol {x})$ is the fluid volume fraction. The fluid indicator function is defined as

(2.7)\begin{equation} I_{f}(\boldsymbol{x})=\begin{cases} 1 & \text{if } \boldsymbol{x} \in \varOmega_f,\\ 0 & \text{else }. \end{cases} \end{equation}

A uniform and symmetric filter kernel is assumed that satisfies

(2.8)\begin{equation} \int_{\varOmega}g(|\boldsymbol{x}|) \,\mathrm{d}V_x = 1. \end{equation}

The kernel possesses a characteristic width $\delta$ that we refer to as the filter width. Applying the volume-filtering operation to the NSE gives (see Anderson & Jackson Reference Anderson and Jackson1967)

(2.9)$$\begin{gather} \dfrac{\partial \epsilon_{{f}}}{ \partial t} + \dfrac{\partial}{\partial x_i}(\epsilon_{{f}}\bar{u}_i) = 0, \end{gather}$$

(2.10)$$\begin{gather}\rho_{{f}}\dfrac{\partial \epsilon_{{f}}\bar{u}_i}{\partial t} + \rho_{{f}}\dfrac{\partial}{\partial x_j}(\epsilon_{{f}}\overline{u_i u_j}) ={-}\dfrac{\partial \epsilon_{{f}}\bar{p}}{\partial x_i}+\mu_{{f}}\dfrac{\partial}{\partial x_j}\left[ \epsilon_{{f}}\left( \overline{\dfrac{\partial u_i}{\partial x_j}} + \overline{\dfrac{\partial u_j}{\partial x_i}} \right) \right] - s_i \end{gather}$$

with the momentum source term expressed as a sum over all particles $q$:

(2.11)\begin{equation} s_i = \sum_q \int_{\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)\hat{n}_j\,\mathrm{d}A_y, \end{equation}

where we dropped the spatial dependence for simplicity. Note that volume-filtered quantities are fundamentally different than their respective unfiltered relatives. A kinetic energy derived from the volume-filtered velocity, e.g. $K_{{VF}}=\epsilon _{{f}}\bar {u}_i\bar {u}_i/2$, does not necessarily obey the same conservation laws as $K=u_i u_i/2$. The volume-filtered scales only represent a portion of the total scales of interest and kinetic energy may be exchanged between the flow present at the volume-filtered scales and the flow present at the smaller subfilter scales. A direct projection of the physical principles to the volume-filtered quantities can lead to wrong conclusions about the consistency of numerical methods for particle-laden flows.

Equations (2.9) and (2.10) contain quantities that are not filtered, namely the fluid velocity and pressure in the particle momentum source term $s_i$, the fluid velocity in the advective term, and the fluid velocity gradient in the viscous term. Therefore, a general solution of these governing equations is not possible without providing closures for the unfiltered terms. The advective term, the second term on the left-hand side of (2.10), is typically replaced with

(2.12)\begin{equation} \dfrac{\partial}{\partial x_j}(\epsilon_{{f}}\bar{u}_i \bar{u}_j), \end{equation}

which we refer to as single volume fraction advective term, as it is commonly used in volume-filtered Euler–Lagrange simulations (see, e.g. Capecelatro & Desjardins Reference Capecelatro and Desjardins2013; Mallouppas & van Wachem Reference Mallouppas and van Wachem2013; Arolla & Desjardins Reference Arolla and Desjardins2015; Ireland & Desjardins Reference Ireland and Desjardins2017). This replacement of the advective term is only valid if the single volume fraction advective term is subtracted and the original advective term of (2.10) is added, which can be realized by introducing a subfilter stress tensor, $\tau _{{sfs},ij}$. Therefore, the subfilter stress tensor accounts for volume filtering the fluid velocity instead of its dyadic product in the advective term. Since the solution of the volume-filtered momentum equation provides an expression for $\epsilon _{{f}}\bar {u}_i$ and $\epsilon _{{f}}\bar {p}$, a division by $\epsilon _{{f}}$ is required to compute the single volume fraction advective term. In cases where the fluid volume fraction is locally very small, small errors of the numerical solution of the volume-filtered momentum equation may lead to instabilities, as is reported by Dave, Herrmann & Kasbaoui (Reference Dave, Herrmann and Kasbaoui2023). A small fluid volume fraction arises for filter widths equal to or smaller than the particle size in the centre of the particle or even for larger filter width to particle size ratios in the case of many closely packed particles. Numerical interventions such as an artificial increase of the filter width are required in these cases (Link et al. Reference Link, Cuypers, Deen and Kuipers2005; Jing et al. Reference Jing, Kwok, Leung and Sobral2016; Chen & Zhang Reference Chen and Zhang2022). Using the identity that the sum of the fluid volume fraction and the particle volume fractions $\epsilon _{{p},q}$ for every particle $q$ equals one,

(2.13)\begin{equation} \epsilon_{{f}} +\sum_q \epsilon_{{p},q} = 1, \end{equation}

we can decompose the single volume fraction advective term as follows:

(2.14)\begin{equation} \dfrac{\partial \epsilon_{{f}}\bar{u}_i\bar{u}_j}{\partial x_j} = \dfrac{\partial \epsilon_{{f}} \bar{u}_i \epsilon_{{f}} \bar{u}_j}{\partial x_j} + \sum_q \dfrac{\partial \epsilon_{{f}}\bar{u}_i \epsilon_{{p},q} \bar{u}_j}{\partial x_j}. \end{equation}

The first term on the right-hand side turns out to be numerically much more suitable as advective term in the volume-filtered NSE than the single volume fraction advective term, since no division by the (potentially very small) fluid volume fraction is required to obtain the filtered fluid flow velocity. We can rewrite the volume-filtered momentum equation (2.10) such that terms that are closed and terms that require closure are separated:

(2.15) $$\begin{align} \rho_{{f}}\dfrac{\partial \epsilon_{{f}}\bar{u}_i}{\partial t} + \rho_{{f}}\dfrac{\partial}{\partial x_j}(\epsilon_{{f}}\bar{u}_i\epsilon_{{f}}\bar{u}_j) &={-}\dfrac{\partial \epsilon_{{f}}\bar{p}}{\partial x_i}+\mu_{{f}} \dfrac{\partial^2\epsilon_{{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_{{sfs},ij}, \end{align}$$

where $\mathcal {E}_i$ is the closure originating from switching the filter and spatial derivative in the viscous term

(2.16)\begin{equation} \mathcal{E}_i = \dfrac{\partial}{\partial x_j}\left[ \epsilon_{{f}}\left( \overline{\dfrac{\partial u_i}{\partial x_j}} + \overline{\dfrac{\partial u_j}{\partial x_i}} \right) \right] - \dfrac{\partial^2\epsilon_{{f}}\bar{u}_i}{\partial x_j \partial x_j}, \end{equation}

and $\tau _{{sfs},ij}$ is the subfilter stress tensor, which arises from the advective term

(2.17)\begin{equation} \tau_{{sfs},ij} = \epsilon_{{f}} \overline{u_i u_j} - \epsilon_{{f}}\bar{u}_i \epsilon_{{f}}\bar{u}_j = \epsilon_{{f}} \overline{u_i u_j} - \epsilon_{{f}}\bar{u}_i\bar{u}_j + \sum_q\epsilon_{{f}}\bar{u}_i \epsilon_{{p},q}\bar{u}_j, \end{equation}

which can be written as a function of $\epsilon _{{f}}\bar {u}_i \epsilon _{{f}}\bar {u}_j$, which we refer to as double volume fraction advective term, or the single volume fraction advective term, $\epsilon _{{f}}\bar {u}_i\bar {u}_j$. Furthermore, $s_{q,i}$ is the particle momentum source originating from the particle with the index $q$ and given as

(2.18)\begin{equation} s_{q,i}= \int_{\partial\varOmega_{{p},q}}g(|\boldsymbol{x}-\boldsymbol{y}|)\left({-}p \delta_{i j} + \mu_{{f}} \left(\dfrac{\partial u_{i}}{\partial y_j}+\dfrac{\partial u_j}{\partial y_{i}}\right)\right)\hat{n}_j\,\mathrm{d}A_y. \end{equation}

Note that (2.18) depends on unfiltered quantities. The total hydrodynamic force acting on each particle, $|\boldsymbol {F}_{{h},q}|=|\int _{\varOmega }\boldsymbol {s}_q\,\mathrm {d}V|$, is typically modelled using empirical correlations. What remains to be closed is the distribution of the hydrodynamic force in space, which is expressed by the regularization kernel $\mathcal {K}_q$. In the direction of the resultant hydrodynamic force, which is indicated by the index $\alpha$, the regularization kernel is defined as

(2.19)\begin{equation} s_{q,\alpha}=\mathcal{K}_q(|\boldsymbol{x}-\boldsymbol{x}_{{p},q}|)F_{{h},q,\alpha}, \end{equation}

where $\boldsymbol {x}_{{p},q}$ is the position of the particle centre of the particle with the index $q$. Note that the particle momentum source may be non-zero in other directions than the direction of the resultant force, although its spatial integral is zero. The regularization kernel $\mathcal {K}_{q}$ is influenced by, e.g. the Reynolds number, velocity gradients, other particles or walls in the environment of the particle. At very large filter width to particle size ratios, however, the filter kernel in (2.18) may be assumed to vary little over the particle surface such that it can be taken outside the integral. Consequently, $\mathcal {K}_q \approx g$ for very large filter widths, which suggests that the aforementioned influences become less important with increasing filter width.

Note that the second term on the right-hand side of (2.14) is absorbed in the subfilter stress tensor, such that (2.15) is equivalent to (2.10). The double volume fraction advective term, in combination with providing the closures in (2.15), allows to use a single-phase flow solver if the solution variables are $u_{\epsilon,i} = \epsilon _{{f}}\bar {u}_i$ and $p_{\epsilon } = \epsilon _{{f}}\bar {p}$, which gives the following volume-filtered NSE:

(2.20)$$\begin{gather} \dfrac{\partial \epsilon_{{f}}}{ \partial t} + \dfrac{\partial u_{\epsilon,i}}{\partial x_i} = 0, \end{gather}$$

(2.21)$$\begin{gather}\rho_{{f}}\dfrac{\partial u_{\epsilon,i}}{\partial t} + \rho_{{f}}\dfrac{\partial}{\partial x_j}(u_{\epsilon,i}u_{\epsilon,j}) ={-}\dfrac{\partial p_{\epsilon}}{\partial x_i} +\mu_{{f}} \dfrac{\partial^2u_{\epsilon,i}}{\partial x_j \partial x_j} -s_i + \mu_{{f}} \mathcal{E}_i- \rho_{{f}}\dfrac{\partial }{\partial x_j}\tau_{{sfs},ij}. \end{gather}$$

Some expressions of the volume-filtered equations in the literature are derived using assumptions, e.g. that filtering an already filtered quantity gives approximately the filtered quantity (Anderson & Jackson Reference Anderson and Jackson1967; Capecelatro & Desjardins Reference Capecelatro and Desjardins2013). Equations (2.20)–(2.21) are derived without a priori assumptions and are exact (given that the closures are provided). Different postulated expressions for the volume-filtered viscous term with different definitions of the viscous closure are discussed in Appendix A.

2.2. Galilean invariance

Specific care has to be taken to ensure Galilean invariance of the momentum equation when closures are derived, which is much less trivial than in the single-phase flow case. Galilean invariance is satisfied if the governing equations are invariant with respect to a constant velocity of the frame of reference $\boldsymbol {u}_{{ref}}$, such that

(2.22)$$\begin{gather} \mathcal{U}_i(\boldsymbol{x},t) = u_i(\boldsymbol{x} - \boldsymbol{u}_{{ref}}t,t) + u_{{ref},i}, \end{gather}$$

(2.23)$$\begin{gather}\mathcal{V}_i(\boldsymbol{x},t) = v_i(\boldsymbol{x} - \boldsymbol{u}_{{ref}}t,t) + u_{{ref},i}, \end{gather}$$

where $\mathcal {U}_i$ and $\mathcal {V}_i$ are the fluid and particle velocity in the fixed frame of reference, and $u_i$ and $v_i$ are the fluid and particle velocity in the frame of reference moving with the velocity $\boldsymbol {u}_{{ref}}$. It can be easily verified that the right-hand side of the volume-filtered momentum equation (2.10) satisfies the Galilean invariance. For the left-hand side, we find

(2.24)\begin{equation} \dfrac{\partial \epsilon_{{f}} \bar{\mathcal{U}}_i}{\partial t} = \dfrac{\partial \epsilon_{{f}} \bar{u}_i}{\partial t} - u_{{ref},j}\dfrac{\partial \epsilon_{{f}}\bar{u}_i}{\partial x_j} - u_{{ref},i}\dfrac{\partial \epsilon_{{f}}\bar{u}_j}{\partial x_j} - u_{{ref},i}u_{{ref},j}\dfrac{\partial \epsilon_{{f}}}{\partial x_j}, \end{equation}

using the continuity equation (2.9), and

(2.25)\begin{equation} \dfrac{\partial \epsilon_{{f}}\overline{\mathcal{U}_i \mathcal{U}_j}}{\partial x_j} = \dfrac{\partial \epsilon_{{f}}\overline{u_i u_j}}{\partial x_j} + u_{{ref},j}\dfrac{\partial \epsilon_{{f}}\bar{u}_i}{\partial x_j} + u_{{ref},i}\dfrac{\partial \epsilon_{{f}}\bar{u}_j}{\partial x_j} + u_{{ref},i}u_{{ref},j}\dfrac{\partial \epsilon_{{f}}}{\partial x_j}. \end{equation}

The sum of both terms leads to the original left-hand side of (2.10), which verifies the Galilean invariance. In the same manner, it can be shown that if we replace the actual advective term, $({\partial }/{\partial x_j})(\epsilon _{{f}}\overline {u_i u_j})$, with the single volume fraction advective term, $({\partial }/{\partial x_j})(\epsilon _{{f}}\bar {u}_i\bar {u}_j)$, Galilean invariance of the momentum equations is achieved even if the subfilter stress tensor is not modelled and fully neglected. With the double volume fraction advective term, $({\partial }/{\partial x_j})(\epsilon _{{f}}\bar {u}_i\epsilon _{{f}}\bar {u}_j)$, the subfilter stress tensor is not Galilean invariant. Consequently, at least one term that cancels out the terms induced by a change of frame of reference has to be included as part of the model for the subfilter stress tensor.

Changing the frame of reference of the double volume fraction advective term causes desired and undesired additional terms:

(2.26) \begin{equation} \dfrac{\partial \epsilon_{{f}} \bar{u}_i \epsilon_{{f}} \bar{u}_j}{\partial x_j} =\underbrace{\dfrac{\partial \epsilon_{{f}}\bar{u}_i\bar{u}_j}{\partial x_j}}_{\text{terms compensated by (2.24)}\ }- \underbrace{\sum_q \dfrac{\partial \epsilon_{{f}}\bar{u}_i \epsilon_{{p},q} \bar{u}_j}{\partial x_j}}_{\text{additional terms}}. \end{equation}

The terms that arise by changing the frame of reference of the first term on the right-hand side are required to compensate the terms arising from the temporal velocity derivative equation (2.24). However, the second term on the right-hand side causes terms that are not compensated without further interventions, i.e. the volume-filtered momentum equation is not Galilean invariant. Therefore, a compensating term has to be added. The minimal expression for compensating the additional terms is

(2.27)\begin{equation} \tau_{{sfs},ij}^{{G}} = \sum_{q}\left[ \epsilon_{{p},q} v_{q,j}\epsilon_{{f}} \bar{u}_i + \epsilon_{{p},q} v_{q,i}\epsilon_{{f}} \bar{u}_j - \epsilon_{{p},q} v_{q,i} \epsilon_{{f}} v_{q,j}\right]. \end{equation}

We consider $\tau _{{sfs},ij}^{{G}}$ as the Galilean invariance part of the subfilter stress tensor. This term does nothing else but compensate the additional terms arising from a change of the frame of reference. It is not a model for the physics in the particle fixed frame, but it ensures that the interscale energy and momentum transfer is the same in all frames of reference. Note that the expression for $\tau _{{sfs},ij}^{{G}}$ given in (2.27) is not unique, but only one possibility to form a Galilean invariant term when adding it to the second term on the right-hand side of (2.26).

The behaviour of the volume-filtered equations when the frame of reference is changed can be intuitively understood when the definition of the volume-filtering equation (2.6) is considered. If the flow quantity $\varPhi$ in (2.6) has the dimension of a velocity, its value changes when the frame of reference is moved. Due to the indicator function, however, the value of the integrand is always zero outside of the fluid region. Consequently, gradients of flow quantities close to the interface are dependent on the frame of reference, which affects all terms in the continuity and momentum equation that contain spatial derivatives of volume-filtered velocities.

A main outcome of the Galilean invariance is that a temporal change of the volume fraction, caused by a change of the frame of reference or a moving particle, leads to additional terms if a term contains spatial gradients of $\epsilon _{{f}}\bar {u}_i$. This has to be considered in the closure modelling for the volume-filtered momentum equation (2.15).

3.1. Viscous closure $\mathcal {E}_i$

The term $\mathcal {E}_i$ arises because the volume-filtered velocity gradients in the viscous stresses of (2.10) are replaced with the gradients of the volume-filtered velocity. Formally, we can write $\mathcal {E}_i$ as the difference between the volume-filtered viscous term and the Laplacian of $\epsilon _{{f}}\bar {u}_i$ (see (2.16)). The effect of $\mathcal {E}_i$ is typically modelled by an additional viscosity (Zhang & Prosperetti Reference Zhang and Prosperetti1997; Patankar & Joseph Reference Patankar and Joseph2001; Capecelatro & Desjardins Reference Capecelatro and Desjardins2013) or neglected (Subramaniam & Balachandar Reference Subramaniam and Balachandar2022). To the best of the authors’ knowledge, however, no study exists that proves the validity of either of the two assumptions. In the present work, we derive an expression for $\mathcal {E}_i$ analytically that is valid without any further assumptions.

We can approach the problem of finding an expression for $\mathcal {E}_i$ to search for a closure for $\epsilon _{{f}}\overline {{\partial u_i}/{\partial x_j}}-({\partial }/{\partial x_j})(\epsilon _{{f}} \bar {u}_i)$. With the definition of the volume-filtering equation (2.6), we obtain

(3.1)\begin{align} \epsilon_{{f}} \overline{\dfrac{\partial u_i}{\partial x_j}}-\dfrac{\partial}{\partial x_j}(\epsilon_{{f}} \bar{u}_i) &= \int_{\varOmega_f}g(|\boldsymbol{x-y}|)\dfrac{\partial u_i(\boldsymbol{y})}{\partial y_j} - \dfrac{\partial g(|\boldsymbol{x-y}|)}{\partial x_j}u_i(\boldsymbol{y})\,{\rm d}V_y \nonumber\\ &= \int_{\varOmega_f}g(|\boldsymbol{x-y}|)\dfrac{\partial u_i(\boldsymbol{y})}{\partial y_j} + \dfrac{\partial g(|\boldsymbol{x-y}|)}{\partial y_j}u_i(\boldsymbol{y})\,{\rm d}V_y \nonumber\\ &= \int_{\varOmega_f} \dfrac{\partial }{\partial y_j}\left( g(|\boldsymbol{x-y}|) u_i(\boldsymbol{y}) \right)\,{\rm d}V_y. \end{align}

For the following analysis, a distinction between the volume fraction of each particle is required if more than one particle is considered. We consider multiple particles by assigning a particle volume fraction $\epsilon _{{p},q}$ to every particle such that (2.13) is satisfied. By applying the divergence theorem, considering that the normal vector $n_j$ points towards the inside of the particle, and the fluid velocity equals the particle velocity at the particle surface, we can further simplify

(3.2)\begin{align} \epsilon_{{f}} \overline{\dfrac{\partial u_i}{\partial x_j}}-\dfrac{\partial}{\partial x_j}(\epsilon_{{f}} \bar{u}_i) &=\int_{\partial \varOmega_{f}}g(|\boldsymbol{x-y}|) u_i(\boldsymbol{y}) n_j \,{\rm d}A_y \nonumber\\ &={-}\sum_q\int_{\partial \varOmega_{{p},q}}g(|\boldsymbol{x-y}|) u_i(\boldsymbol{y}) \hat{n}_j \,{\rm d}A_y \nonumber\\ &={-}\sum_q v_{q,i}\int_{\partial \varOmega_{{p},q}}g(|\boldsymbol{x-y}|)\hat{n}_j \,{\rm d}A_y \nonumber\\ &= \sum_q v_{q,i}\int_{\varOmega_{{p},q}}\dfrac{\partial g(|\boldsymbol{x-y}|)}{\partial x_j}\,{\rm d}V_y \nonumber\\ &= \sum_q v_{q,i}\dfrac{\partial \epsilon_{{p},q}}{\partial x_j}, \end{align}

where $\hat {n}_j$ represents the normal vector pointing into the fluid. Although particle rotation is excluded in the present study, it is worth noting that a pure rotation of the spherical particle does not contribute to the closure (see Appendix B). From (3.2), we can conclude that if the volume fraction does not change in time, the filtered velocity gradient and the gradient of the filtered velocity are equivalent.

Since the volume of each particle stays constant, we can write

(3.3)\begin{equation} \dfrac{\mathrm{D}\epsilon_{{p},q}}{\mathrm{D}t} = 0 = \dfrac{\partial \epsilon_{{p},q}}{\partial t} + v_{q,i}\dfrac{\partial \epsilon_{{p},q}}{\partial x_i}. \end{equation}

Using the conservation of volume fraction, (2.13), we get

(3.4)\begin{equation} \dfrac{\partial \epsilon_{{f}}}{\partial t} = \dfrac{\partial (1-\sum_q \epsilon_{{p},q})}{\partial t} = \sum_q v_{q,i}\dfrac{\partial \epsilon_{{p},q}}{\partial x_i}, \end{equation}

and by incorporating the continuity equation (2.9),

(3.5)\begin{equation} \dfrac{\partial \epsilon_{{f}} \bar{u}_i}{\partial x_i} ={-}\sum_q v_{q,i}\dfrac{\partial \epsilon_{{p},q}}{\partial x_i}. \end{equation}

By inserting (3.2) into the definition of the viscous closure equation (2.16), we obtain

(3.6)\begin{align} \mathcal{E}_i &= \dfrac{\partial}{\partial x_j}\left[ \dfrac{\partial \epsilon_{{f}} \bar{u}_i}{\partial x_j} + \sum_q v_{q,i}\dfrac{\partial \epsilon_{{p},q}}{\partial x_j} + \dfrac{ \partial \epsilon_{{f}} \bar{u}_j}{\partial x_i} +\sum_q v_{q,j}\dfrac{\partial \epsilon_{{p},q}}{\partial x_i}\right] - \dfrac{\partial^2\epsilon_{{f}}\bar{u}_i}{\partial x_j \partial x_j} \nonumber\\ &= \dfrac{\partial}{\partial x_j}\left[ \sum_q v_{q,i}\dfrac{\partial \epsilon_{{p},q}}{\partial x_j} + \dfrac{ \partial \epsilon_{{f}} \bar{u}_j}{\partial x_i} +\sum_q v_{q,j}\dfrac{\partial \epsilon_{{p},q}}{\partial x_i}\right] \nonumber\\ &= \sum_q v_{q,i}\dfrac{\partial^2 \epsilon_{{p},q}}{\partial x_j \partial x_j} + \dfrac{\partial}{\partial x_i}\left( \dfrac{ \partial \epsilon_{{f}} \bar{u}_j}{\partial x_j} \right)+\sum_q v_{q,j}\dfrac{\partial^2 \epsilon_{{p},q}}{\partial x_i \partial x_j} \nonumber\\ &= \sum_q v_{q,i}\dfrac{\partial^2 \epsilon_{{p},q}}{\partial x_j \partial x_j}, \end{align}

where we used (3.5). In numerical simulations, the particle volume fraction is an easily accessible quantity with a known analytical expression for spherical particles that can be found, for instance, in Balachandar & Liu (Reference Balachandar and Liu2022), using a Gaussian filter kernel $g$. Its second spatial derivative can be computed numerically or analytically and $\mathcal {E}_i$ can be added as an explicit source term in the fluid momentum equation (2.15). Note that the expression for the viscous closure is also valid for different filter kernels but the particle volume-fraction needs to be adjusted accordingly.

Together with the closure $\mathcal {E}_i$, the viscous contribution of (2.15) is Galilean invariant. Consequently, the closure must be included to ensure a physical representation of the volume-filtered fluid–particle system when the volume fraction changes in time. In the case of fixed particles, the viscous closure is zero.

3.2. Particle momentum source $s_i$

The particle momentum source term, $s_i$, represents the momentum exchange between the particles and the filtered flow scales. We exclude the discussion of the modelling of the hydrodynamic force on the particle, $F_{{h},i}$. Instead, we discuss the regularization kernel, i.e. how $s_i$ varies in space, which is typically approximated with a Gaussian centred at the particle centre (referred to as Gaussian regularization kernel), a polynomial approximation of a Gaussian or a top-hat function (Crowe et al. Reference Crowe, Sharma and Stock1977; Maxey & Patel Reference Maxey and Patel2001; Maxey Reference Maxey2017; Poustis et al. Reference Poustis, Senoner, Zuzio and Villedieu2019; Evrard et al. Reference Evrard, Denner and van Wachem2020; Keane et al. Reference Keane, Apte, Jain and Khanwale2023). It is commonly argued that if the filter width is much larger than the size of the particle, i.e. $\delta \gg a$, the filter kernel $g$ varies insignificantly across the volume of the particle and can be treated as a constant for the integration (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013):

(3.7)\begin{align} s_i &= \sum_q \int_{\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 \nonumber\\ &\approx \sum_q F_{{h},q,i}g(|\boldsymbol{x}-\boldsymbol{x}_{{p},q}|), \end{align}

where $\boldsymbol {x}_{{p},q}$ is the centre of the respective particle. It is not clear, however, at which filter width this is a valid assumption and what the effect of this approximation is on the velocity field. This is discussed in more detail in § 4.

For a general expression for $s_i$, the exact pressure and velocity field at the surface of the particle is required, and its convolution with the filter kernel $g$ must be integrable. We find such an expression in the limit of ${Re}\rightarrow 0$, where the flow is governed by the Stokes flow equations. By exploiting that the fluid stresses in Stokes flow are constant over the whole spherical particle surface with radius $a$, we get

(3.8)\begin{align} s_i &= \int_{\partial\varOmega_{{p}}}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 \nonumber\\ &= \dfrac{3\mu_{{f}}(u_{\infty,i}-v_{i})}{2a}\int_{\partial\varOmega_{{p}}}g(|\boldsymbol{x}-\boldsymbol{y}|)\,\mathrm{d}A_y \nonumber\\ &= F_{{h},i} \mathcal{K}_{{St}}(|\boldsymbol{x}-\boldsymbol{x}_{{p}}|), \end{align}

where $u_{\infty,i}$ is the fluid velocity far away from the particle. If $g$ is assumed to be Gaussian,

(3.9)\begin{equation} g(\boldsymbol{x}) = \dfrac{1}{(2{\rm \pi} \sigma^2)^{3/2}}\exp\left( -\dfrac{|\boldsymbol{x}|^2}{2 \sigma^2} \right), \end{equation}

the explicit evaluation of the integral in (3.8) gives the regularization kernel for an isolated spherical particle in Stokes flow

(3.10)\begin{equation} \mathcal{K}_{{St}}(r) = \dfrac{\exp{\left(-\dfrac{(a+r)^2}{2\sigma^2}\right)}}{4{\rm \pi}\sqrt{2 {\rm \pi}\sigma^2}ra} \left( \exp{\left(\dfrac{2ar}{\sigma^2}\right)} - 1\right), \end{equation}

where $\sigma$ is the standard deviation of the Gaussian. For the remainder of this work, we assume the relation $\sigma = \delta \sqrt {2/9{\rm \pi} }$. This somewhat arbitrary choice of the filter width is equal to the support of the Wendland kernel, a practical relevant polynomial approximation of a Gaussian, which is used, e.g. by Evrard et al. (Reference Evrard, Denner and van Wachem2020) to regularize the particle momentum source.

Note that an expression for the regularization kernel is not limited to a Gaussian filter kernel, but can be obtained for every filter kernel that can be integrated according to (3.8). It has to be highlighted that the kernel $\mathcal {K}_{{St}}$ is only valid in uniform Stokes flow over a single sphere in an infinite domain. At higher Reynolds number, in shear flow, in the vicinity of a wall, or if multiple particles are present, the particle momentum source is distributed differently in space. Note that although the Stokes flow equations are linear, that does not allow for a superposition of the kernels $\mathcal {K}_{{St}}$ for multiple spheres, as the fluid stresses are no longer constant over the particle surface.

To assess if the commonly applied Gaussian kernel is a suitable choice for the regularization of the momentum source, it is compared with the analytical regularization kernel in Stokes flow. In figure 1, the analytical regularization kernel of Stokes flow around a sphere $\mathcal {K}_{{St}}$ is shown for different filter width to particle size ratios along with the Gaussian kernel $g$ that is commonly used to regularize the force in point-particle simulations. The Gaussian and the analytical kernel are normalized with the volume of the spherical particle, $V_{{p}}=4{\rm \pi} a^3/3$. The largest deviations are observed for the smallest filter widths, where the Gaussian has much larger values than the analytical kernel. An interesting property of the analytical kernel is that for $\delta /a<3$, it is centred at the particle surface, whereas the Gaussian kernel is centred at the particle centre. As the filter width increases, the magnitudes and shapes of the two kernels approach each other and at filter widths of approximately $\delta /a\ge 8$, only marginal differences between the kernels can be observed. For smaller filter widths, however, significant deviations are observed.

Figure 1. Gaussian kernel (dashed lines) and the analytical kernel in the Stokes flow (solid lines) for different filter widths. The colour indicates the different filter widths. (a) For $1\leq \delta /a \leq 5$ and (b) for $5 \leq \delta /a \leq 10$.

As the fluid stresses become non-uniformly distributed over the particle surface because of a higher Reynolds number or the presence of other particles, the kernel looses its spherical symmetry. In such a case, $\mathcal {K}_{{St}}$ becomes an approximation. The Gaussian kernel, however, is an approximation for all ${Re}$. As shown in figure 1, the Gaussian kernel is a very poor approximation of the real regularization kernel at small ${Re}$ and small filter width to particle size ratios ($\delta /a=O(1)$).

3.3. Subfilter stress tensor $\tau _{{sfs},ij}$

The subfilter stress tensor, $\tau _{{sfs},ij}$, originates from the difference of the volume-filtered tensor product of the fluid velocity vectors and the tensor product of the volume-filtered fluid velocity vectors. With the choice of the advective term in the volume-filtered fluid momentum equation (2.10), a modelled subfilter stress tensor needs to at least contain a contribution to satisfy Galilean invariance $\tau _{{sfs}, ij}^{{G}}$, which is discussed in § 2.2. Furthermore, the nonlinear subfilter interactions are modelled with $\tau _{{sfs}, ij}^{{inter}}$, which is Galilean invariant itself, such that

(3.11)\begin{equation} \tau_{{sfs}, ij}^{{mod}} = \tau_{{sfs}, ij}^{{G}} + \tau_{{sfs}, ij}^{{inter}}. \end{equation}

While the expression to compensate for additional terms arising from the change of the frame of reference $\tau _{{sfs}, ij}^{{G}}$ is known, the nonlinear subfilter interactions $\tau _{{sfs}, ij}^{{inter}}$ require modelling.

A major difficulty for understanding and modelling the subfilter stress tensor is that, unlike other closures, it can be non-zero also in the absence of particles. When a turbulent single-phase flow with a filter width larger than the Kolmogorov length scale is considered, as is the case in classical LES, $\tau _{{sfs},ij}$ accounts for the mainly dissipative effect of the unresolved turbulent velocity on the filtered scales. The effect of the subfilter stresses on the filtered scales is relatively well understood in single-phase flow cases such as homogeneous isotropic turbulence (see, e.g. Sagaut Reference Sagaut2006). However, if particles are present, one does not even fully understand the effect of a single particle subject to laminar flow on the subfilter stress tensor, and especially not the effects of an additional turbulent background flow or of fluid velocity disturbances due to neighbouring particles. Mehrabadi et al. (Reference Mehrabadi, Tenneti, Garg and Subramaniam2015) studied the behaviour of the pseudo-turbulent Reynolds stresses, $R_{ij}$, with different configurations of particle-resolved direct numerical simulations. The pseudo-turbulent Reynolds stresses, $R_{ij}$, are defined as

(3.12)\begin{equation} R_{ij} = \epsilon_{{f}} \overline{u_i^{\prime} u_j^{\prime}}, \end{equation}

where the subfilter velocity is given by

(3.13)\begin{equation} u_i^{\prime} = u_i - \bar{u}_i. \end{equation}

Note that Mehrabadi et al. (Reference Mehrabadi, Tenneti, Garg and Subramaniam2015) define the pseudo-turbulent Reynolds stresses based on ensemble averaging, which is related, but not equivalent to, the volume-filtering studied in the present work. There have been attempts at modelling the trace of $R_{ij}$, the pseudo-turbulent kinetic energy, and reconstruct $R_{ij}$ based on it (Mehrabadi et al. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015; Shallcross et al. Reference Shallcross, Fox and Capecelatro2020). However, finding a model for the pseudo-turbulent Reynolds stresses in the context of volume filtering is not equivalent to finding a closure for the subfilter stress tensor. The first term of $\tau _{{sfs},ij}$ can be decomposed as follows:

(3.14)\begin{equation} \epsilon_{{f}}\overline{u_i u_j} = \epsilon_{{f}}\overline{ \bar{u}_i \bar{u}_j} + \epsilon_{{f}}\overline{ \bar{u}_i u^{\prime}_j}+ \epsilon_{{f}}\overline{ u^{\prime}_i\bar{u}_j } + R_{ij}. \end{equation}

Consequently, to model $\tau _{{sfs},ij}$, a model for $R_{ij}$ is not sufficient. Although the importance of the remaining three terms is unknown for a particle-laden flow, it is known from single-phase turbulence that the contribution of these three terms is not negligible and, in some configurations, even dominate the contribution of $R_{ij}$ (Hunt & Carruthers Reference Hunt and Carruthers1990; Mann Reference Mann1994; Laval, Dubrulle & Nazarenko Reference Laval, Dubrulle and Nazarenko2001; Ghate & Lele Reference Ghate and Lele2020).

There have been suggestions to model the subfilter stress tensor with an additional viscosity according to the Boussinesq hypothesis in turbulence (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013; Subramaniam & Balachandar Reference Subramaniam and Balachandar2022). Hausmann et al. (Reference Hausmann, Evrard and van Wachem2023) derived an eddy viscosity model for LES of particle-laden turbulent flows that accounts for the turbulence modulation by the particles, but with the assumption that the effect of the sub-Kolmogorov particles on the subfilter stress tensor is negligible if the turbulence is fully resolved.

To indicate the filter width of a volume-filtered flow quantity weighted with the fluid indicator function, $I_{{f}}\varPhi$, we introduce the following notation:

(3.15)\begin{equation} \overline{I_{{f}}\varPhi}\vert_{\sigma} = \int_{\varOmega} I_{{f}}(\boldsymbol{y}) \varPhi(\boldsymbol{y})g(|\boldsymbol{x}-\boldsymbol{y}|)\,\mathrm{d}V_y, \end{equation}

where $(\cdot )\vert _{\sigma }$ indicates that the flow quantity is filtered using a Gaussian with a standard deviation $\sigma$. Assuming that the filtering kernel $g$ is Gaussian, it can be verified that the volume-filtering operation with (2.6) is equivalent to solving the diffusion equation for a flow quantity weighted with the fluid indicator function, $I_{{f}}\varPhi$, (see, e.g. Johnson Reference Johnson2021)

(3.16)\begin{equation} \dfrac{\partial \overline{I_{{f}}\varPhi}\vert_{\sigma}}{\partial (\sigma^2)} = \dfrac{1}{2}\nabla^2\overline{I_{{f}}\varPhi}\vert_{\sigma}, \end{equation}

with the initial condition

(3.17)\begin{equation} \overline{I_{{f}}\varPhi}\vert_{\sigma=0} = I_{{f}}\varPhi. \end{equation}

The filtered result is

(3.18)\begin{equation} \overline{I_{{f}}\varPhi}\vert_{\sigma} = \epsilon_{{f}} \bar{\varPhi}. \end{equation}

With the definition of the nonlinear advective closure equation (2.17), we can find the exact relation (see Appendix C)

(3.19)\begin{equation} \dfrac{\partial \tau_{{sfs},ij}\vert_{\sigma}}{\partial (\sigma^2)} = \dfrac{1}{2}\nabla^2\tau_{{sfs},ij}\vert_{\sigma} + \dfrac{\partial \overline{I_{{f}}u_i}\vert_{\sigma}}{\partial x_k}\dfrac{\partial \overline{I_{{f}}u_j}\vert_{\sigma}}{\partial x_k}, \end{equation}

with the initial condition

(3.20)\begin{equation} \tau_{{sfs},ij}\vert_{\sigma=0} = 0. \end{equation}

Equation (3.19) is analogous to the expression obtained from Johnson (Reference Johnson2021) for the subfilter stress tensor in a single-phase flow. Since the Green's function $G$ of (3.19) is given as

(3.21)\begin{equation} G\vert_{\sigma}(|\boldsymbol{x}|) = \begin{cases} g\vert_{\sigma}(|\boldsymbol{x}|), & \sigma^2\geq 0,\\ 0, & \sigma^2<0, \end{cases} \end{equation}

where $g\vert _{\sigma }$ is the Gaussian of standard deviation $\sigma$, its formal solution reads as the convolution integral

(3.22)\begin{equation} \tau_{{sfs},ij}\vert_{\sigma}(\boldsymbol{x}) = \int_0^{\sigma^2} \int_{\varOmega} G\vert_{\sqrt{\sigma^2-\theta}}(|\boldsymbol{x}-\boldsymbol{s}|)\dfrac{\partial \overline{I_{{f}}u_i(\boldsymbol{s})}\vert_{\sqrt{\theta}}}{\partial x_k}\dfrac{\partial \overline{I_{{f}}u_j(\boldsymbol{s})}\vert_{\sqrt{\theta}}}{\partial x_k}\,\mathrm{d}\boldsymbol{s}\,\mathrm{d}\theta. \end{equation}

After evaluating the inner integral, which is essentially another filtering operation, we obtain

(3.23)\begin{align} \tau_{{sfs},ij}\vert_{\sigma} &= \int_0^{\sigma^2} \left.\overline{ \dfrac{\partial\overline{ I_{{f}}u_i}\vert_{\sqrt{\theta}}}{\partial x_k} \dfrac{\partial\overline{ I_{{f}}u_j}\vert_{\sqrt{\theta}}}{\partial x_k}}\right\vert_{\sqrt{\sigma^2-\theta}}\,\mathrm{d}\theta \nonumber\\ &= \sigma^2 \left.\overline{I_{{f}}\dfrac{\partial u_i}{\partial x_k}}\right\vert_{\sigma}\left.\overline{I_{{f}}\dfrac{\partial u_j}{\partial x_k}}\right\vert_{\sigma} + \int_0^{\sigma^2} \left.\overline{ \dfrac{\partial\overline{ I_{{f}}u_i}\vert_{\sqrt{\theta}}}{\partial x_k} \dfrac{\partial\overline{ I_{{f}}u_j}\vert_{\sqrt{\theta}}}{\partial x_k}}\right\vert_{\sqrt{\sigma^2-\theta}}\,\mathrm{d}\theta \nonumber\\ &\quad -\int_0^{\sigma^2} \left.\overline{ \left.\overline{I_{{f}}\dfrac{\partial u_i}{\partial x_k}}\right\vert_{\sqrt{\theta}} }\right\vert_{\sqrt{\sigma^2-\theta}} \left.\overline{ \left.\overline{I_{{f}}\dfrac{\partial u_j}{\partial x_k}}\right\vert_{\sqrt{\theta}} }\right\vert_{\sqrt{\sigma^2-\theta}}\,\mathrm{d}\theta, \end{align}

where we use the relation (see, e.g. Johnson Reference Johnson2021)

(3.24)\begin{equation} \overline{\overline{I_{{f}}u}_i\vert_{\sigma_1}}\vert_{\sigma_2} = \overline{I_{{f}}u}_i\vert_{\sqrt{\sigma_1^2 + \sigma_2^2}}. \end{equation}

The first term on the right-hand side of (3.23) contains only filtered quantities, whereas the remaining terms also consist of the fluid velocity at scales that are smaller than $\sigma$. The first and the last term on the right-hand side of equation (3.23) are Galilean invariant, but the second term is not. The subfilter stress tensor, as defined in (2.17), is not Galilean invariant. For single-phase turbulence, the contribution of the integral terms to the kinetic energy transfer are minor compared with the contribution of the first term for small filter widths. The first term can be computed from knowledge of volume-filtered quantities only, which essentially leads to a well-known LES model for single-phase flow turbulence, sometimes referred to as the nonlinear model (Liu, Meneveau & Katz Reference Liu, Meneveau and Katz1994; Borue & Orszag Reference Borue and Orszag1998). This model was first derived by Leonard (Reference Leonard1975) with a Taylor expansion of the filtered velocity in space. It has been shown by Johnson (Reference Johnson2021) that the resolved contribution of (3.23) is dominant for small filter widths in single-phase flow turbulence. Borue & Orszag (Reference Borue and Orszag1998) showed that the nonlinear model strongly correlates with the subfilter stress tensor. Therefore, the subfilter stress tensor is proposed to be modelled as the resolved contribution of (3.23), i.e. the first term on the right-hand side of (3.23),

(3.25)\begin{equation} \tau_{{sfs},ij}^{{mod,NL}} = \tau_{{sfs},ij}^{{G}} + \sigma^2 \left[ \dfrac{\partial \epsilon_{{f}} \bar{u}_i }{\partial x_k} + \sum_q v_{q,i} \dfrac{\partial \epsilon_{{p},q}}{\partial x_k} \right]\left[ \dfrac{\partial \epsilon_{{f}} \bar{u}_j }{\partial x_k} + \sum_q v_{q,j}\dfrac{\partial \epsilon_{{p},q}}{\partial x_k} \right]. \end{equation}

We refer to this model as the nonlinear model since it is a generalization of the model proposed by Leonard (Reference Leonard1975). The nonlinear model can be formulated as a tensorial particle viscosity model in analogy to a turbulent viscosity in turbulent single-phase flow. For comparison, we adapt the most common scalar turbulent viscosity model, the Smagorinsky model, to a volume-filtered particle-laden flow. The adaption is required to ensure Galilean invariance of the model. Using the identity equation (3.2), we can write the modelled subfilter stress tensor analogously to the Smagorinsky model in turbulent single-phase flow,

(3.26)\begin{equation} \tau_{{sfs},ij}^{{inter,Sm}} - \dfrac{1}{3}\tau_{{sfs},kk}^{{inter,Sm}}\delta_{ij} ={-} 2\nu_{{p}}\bar{S}_{ij}, \end{equation}

with the additional viscosity that we refer to as particle viscosity,

(3.27)\begin{equation} \nu_{{p}} = (C_{{s,p}}\sigma)^2\sqrt{2\bar{S}_{ij}\bar{S}_{ij}}, \end{equation}

where $C_{{s,p}}$ is a constant. The Galilean invariant volume-filtered strain-rate tensor is given by

(3.28)\begin{equation} \bar{S}_{ij} = \dfrac{1}{2}\left[ \dfrac{\partial \epsilon_{{f}} \bar{u}_i }{\partial x_j} + \sum_q v_{q,i} \dfrac{\partial \epsilon_{{p},q}}{\partial x_j} + \dfrac{\partial \epsilon_{{f}} \bar{u}_j }{\partial x_i} + \sum_q v_{q,j} \dfrac{\partial \epsilon_{{p},q}}{\partial x_i}\right], \end{equation}

where the volume-filtered velocity gradients are replaced using (3.2). Note that the Smagorinsky model, such as any model for the subfilter stress tensor relying on the Boussinesq hypothesis, assumes an alignment of the eigenvectors of the subfilter stress tensor with the eigenvectors of the strain-rate tensor. At least in single-phase turbulence, it is known that this alignment is wrong (Horiuti Reference Horiuti2003). The Smagorinsky model adapted to the volume-filtered framework is, therefore, not expected to capture the correct directional dependencies in anisotropic flows.

4.1. Considered cases

To test the implications of the previous section and to assess the importance of the different closures in the volume-filtered momentum equation, we study two different configurations of particles immersed in different flows. The first configuration is an isolated sphere in uniform flow at different Reynolds numbers. In Stokes flow, the existence of an analytical solution for this configuration allows us to study the effect of volume filtering in an infinite domain without a discretization error. We extend this configuration to higher Reynolds numbers (${Re}=20$ and ${Re}=100$) and refer to this first configuration as a single sphere configuration. With the second configuration, we study the influence of neighbouring particles by simulating a sphere in a fully periodic computational domain, which we refer to as a periodic array of spheres.

The cases of finite ${Re}$ require numerical simulations, which are carried out with our in-house flow solver. The flow is solved with a finite-volume method, which employs momentum-weighted interpolation for a coupled solution of continuity and momentum that is second order in time and space (Bartholomew et al. Reference Bartholomew, Denner, Abdol-Azis, Marquis and van Wachem2018; Denner, Evrard & van Wachem Reference Denner, Evrard and van Wachem2020). The no-slip boundary condition at the particle surface is enforced with an immersed-boundary method (IBM) as described by Abdol Azis, Evrard & van Wachem (Reference Abdol Azis, Evrard and van Wachem2019) and Chéron, Evrard & van Wachem (Reference Chéron, Evrard and van Wachem2023b). Chéron et al. (Reference Chéron, Evrard and van Wachem2024) perform validations to obtain temporal and spatial resolutions that lead to converged results. In the present work, we adapt these resolutions.

Figure 2(a) shows a sketch of the simulation domain of the single sphere configuration. The simulation domain is a cube of size $40a$ with the particle placed in the centre and the flow is driven by a uniform inflow with the velocity $u_{\infty }$ at the lower face of the steamwise direction, $x_1$. The other velocities and the pressure gradient are zero. The same boundary conditions are applied at the faces in the two normal directions, $x_2$ and $x_3$. At the upper streamwise face, the velocity gradients and the pressure are set to zero.

Figure 2. Sketch of the simulation domain of (a) the single sphere configuration and (b) the periodic array of spheres.

The simulation domain of the periodic array of spheres is depicted in figure 2(b). The cubic periodic simulation domain has a size of $4a$, which corresponds to a global particle volume fraction of $\epsilon _{{p,glob}}=0.065$. The flow is driven by a constant volume force, $\delta f$, applied in the streamwise direction. The Reynolds number of the periodic array of spheres is defined as

(4.1)\begin{equation} {Re}_{{periodic}} = \dfrac{\rho_{{f}}|\boldsymbol{u}_{{sf}}|2a}{\mu_{{f}}}, \end{equation}

where the superficial velocity is given by

(4.2)\begin{equation} u_{{sf},i} = \dfrac{1}{V_{{periodic}}}\int _{\varOmega_{{periodic}}} I_{{f}}u_i \,\mathrm{d}V. \end{equation}

The volume of the periodic simulation domain $\varOmega _{{periodic}}$, which includes the volume occupied by the fluid and the particle, is indicated with $V_{{periodic}}$.

For the single sphere configuration, adaptive mesh refinement is used to obtain a spatial resolution corresponding to 60 fluid mesh cells per particle diameter. Far away from the particle, the fluid mesh is significantly coarser, because the flow is essentially uniform there. In the periodic array of spheres, the fluid mesh is equidistant with 64 fluid mesh cells per particle diameter. The time step, $\Delta t$, for the single sphere configuration is chosen such that $\mathrm {CFL}=u_{{max}}\Delta t/\Delta x_{{min}}<0.02$ and for the periodic arrays of spheres, $\mathrm {CFL}=0.14$, where $u_{{max}}$ is the maximum velocity in the domain and $\Delta x_{{min}}$ is the smallest fluid mesh spacing.

The volume-filtered results are obtained by explicit volume filtering the flow fields from the particle-resolved simulations, which is detailed in Appendix D.

The filter widths studied in the remainder of the present section range from $\delta /a = O(1)$ to $\delta /a = O(10)$. Filter widths of the order $\delta /a = O(1)$ are not applicable to point-particle simulations because of major modelling errors with different origins. However, $\delta /a = O(1)$ is still an interesting regime to study to guide new simulation methods that serve as a bridge between particle-resolved simulations and point-particle simulations, and to gain a general understanding of the effect of volume-filtering on the flow quantities.

4.2. Stokes flow around a single sphere

We consider the Stokes flow around a sphere in an infinite domain. In one configuration, the particle steadily moves with a constant velocity in the negative streamwise direction in a quiescent fluid and in the other configuration, the particle is fixed but there is a constant positive fluid velocity. Although these cases are analogous, by considering a different reference frame in which the sphere moves, the viscous closure, $\mathcal {E}$, may be investigated, which vanishes if the sphere does not move with respect to the reference frame.

We can define a kinetic energy of the volume-filtered fluid velocity $K_{{VF}}=\epsilon _{{f}}\bar {u}_i\bar {u}_i/2$. By taking the dot product of the volume-filtered fluid momentum equation (2.15) with $\bar {u}_i$, we obtain a transport equation for the kinetic energy of the volume-filtered fluid velocity

(4.3)$$\begin{gather} 2\rho_{{f}}\dfrac{\partial K_{{VF}}}{\partial t} - \rho_{{f}}\epsilon_{{f}}\bar{u}_i\dfrac{\partial \bar{u}_i}{\partial t} + \bar{u}_i\rho_{{f}}\dfrac{\partial}{\partial x_j}(\epsilon_{{f}}\bar{u}_i\epsilon_{{f}}\bar{u}_j) ={-} \bar{u}_i\dfrac{\partial \epsilon_{{f}}\bar{p}}{\partial x_i}+\mu_{{f}} \bar{u}_i\dfrac{\partial^2\epsilon_{{f}}\bar{u}_i}{\partial x_j \partial x_j} - \bar{u}_i s_i \nonumber\\ +\mu_{{f}} \bar{u}_i\mathcal{E}_i - \rho_{{f}} \bar{u}_i\dfrac{\partial }{\partial x_j}\tau_{{sfs},ij}. \end{gather}$$

In the limit of steady Stokes flow, the transport equation reduces to

(4.4)\begin{equation} 0 ={-} \bar{u}_i\dfrac{\partial \epsilon_{{f}}\bar{p}}{\partial x_i}+\mu_{{f}} \bar{u}_i\dfrac{\partial^2\epsilon_{{f}}\bar{u}_i}{\partial x_j \partial x_j} - \bar{u}_i F_{{h},i}\mathcal{K}_{{St}}(|\boldsymbol{x}-\boldsymbol{x}_{{p}}|)+ \mu_{{f}} \bar{u}_i\mathcal{E}_i. \end{equation}

For conciseness, we abbreviate the pressure term $\mathcal {P}_i={\partial \epsilon _{{f}}\bar {p}}/{\partial x_i}$ and the viscous term $\mathcal {V}_i={\partial ^2\epsilon _{{f}}\bar {u}_i}/{\partial x_j \partial x_j}$. In figure 3, the different volume integrated energy transfer rates $\varPsi$ that occur in the kinetic energy transport equation are plotted for the flow around a fixed sphere. The energy transfer rate of a quantity $\varPhi _i$ in the volume-filtered fluid momentum equation (2.15) is given as

(4.5)\begin{equation} \varPsi = \int _{\varOmega_{{int}}}\bar{u}_i \varPhi_i \,\mathrm{d}V, \end{equation}

where $\varOmega _{{int}}$ is the integration region, which corresponds to a cube with an edge length of ten particle diameters, i.e. $20a$. The trends with respect to the filter width of the presented results are verified to be insensitive to $\varOmega _{{int}}$, although the magnitudes of some energy transfer rates change. Therefore, the results must be interpreted as the energetic effect of particles in an influence region $\varOmega _{{int}}$, which is the most interesting region because the energy transfer rates possess the largest magnitudes there.

Figure 3. Volume integrated energy transfer rates of the flow over a fixed sphere in Stokes flow. The energy transfer rates of the particle momentum source, the viscous term and the pressure term are plotted.

When the particle is fixed and the filter width $\delta =0$, no kinetic energy is exchanged between fluid and particle because the velocity at the particle surface is zero. The energy transfer rate of the particle momentum source $\varPsi =-F_{{h},i}v_i$ is zero. With increased filter width, the energy transfer rate is negative. Consequently, the energy transfer rate of the particle momentum source removes kinetic energy from the filtered scales and adds kinetic energy to the subfilter scales, such that the total energy transfer rate of all scales equals zero. This is in alignment with the observation that particles in turbulence remove kinetic energy from large scales and add them to smaller scales (Squires & Eaton Reference Squires and Eaton1994; Sundaram & Collins Reference Sundaram and Collins1999; Ferrante & Elghobashi Reference Ferrante and Elghobashi2003; Mallouppas, George & van Wachem Reference Mallouppas, George and van Wachem2017).

We further observe that the energy transfer rate induced by the pressure term remains relatively constant as the filter width increases indicating that the pressure energy exchange takes place at large scales. The energy transfer rate of the viscous term possesses its minimum at $\delta =0$. Most of the viscous dissipation takes place at a small filter width. For larger filter widths, the viscous source term even becomes positive, which seems counter-intuitive. The viscous energy source term can be decomposed as

(4.6)\begin{equation} \int_{\varOmega}\mu_{{f}} \bar{u}_i\dfrac{\partial^2\epsilon_{{f}}\bar{u}_i}{\partial x_j \partial x_j}\,\mathrm{d}V = \int_{\varOmega}\dfrac{\mu_{{f}}}{2\epsilon_{{f}}} \dfrac{\partial^2 (\epsilon_{{f}}\bar{u}_i)^2}{\partial x_j \partial x_j}\,\mathrm{d}V - \int_{\varOmega}\dfrac{\mu_{{f}}}{\epsilon_{{f}}}\dfrac{\partial \epsilon_{{f}}\bar{u}_i}{\partial x_j}\dfrac{\partial \epsilon_{{f}}\bar{u}_i}{\partial x_j}\,\mathrm{d}V. \end{equation}

The second term on the right-hand side is always negative and analogous to the viscous fluid flow dissipation which is also present in a single-phase flow. In the fixed particle configuration, the curvature of the squared volume-filtered velocity is positive and the local fluid volume fraction is smaller than one. Consequently, positive curvatures of the squared volume-filtered velocity receive a larger weight than negative curvatures, and the first term on the right-hand side is positive and even outweighs the second term for large filter widths.

In figure 4, the frame of reference is shifted, such that the particle moves in a quiescent fluid. The energy transfer rate induced by the viscous closure is compared with the energy transfer rate of the total viscous term. The largest contribution of the viscous closure to the energy transfer rate is at $\delta \approx a$, where it makes up approximately 20 % of the energy transfer rate of the viscous term. For larger filter widths, the energy transfer rate induced by the viscous closure approaches zero. In the fixed sphere configuration, the closure is zero because the particle translational velocity is zero, i.e. the volume fraction does not change.

Figure 4. Volume integrated energy transfer rates of a sphere moving in quiescent fluid in the Stokes regime. The energy transfer rates of the viscous closure term, the viscous term and their sum are plotted.

We consider the moving particle in the Stokes flow to evaluate the spatial distribution of the viscous closure $\mathcal {E}_i$. Figure 5 compares the contours in a plane normal to the streamwise direction of the viscous closure computed by subtracting the viscous stresses as in (2.16) using quantities not accessible in the volume-filtered framework and computed using the exact relation of (3.6) containing only quantities that are known in the volume-filtered framework. The results are normalized with $f_{{ref}}=|\boldsymbol {F}_{{h}}|/V_{{p}}$. The agreement of the two contours confirms the derivation in § 3.1. Note that $\mathcal {E}_i$ is spherically symmetric, which is always guaranteed for spherical particles. Since the particle volume fraction corresponding to a spherical particle is spherical symmetric, its Laplacian possesses spherical symmetry too. The derivation of the expression for $\mathcal {E}_i$ shows that it is valid also outside the Stokes regime and in the vicinity of other particles.

Figure 5. Contours of the viscous closure at filter width $\delta /a= 1$. (a) Closure by explicitly subtracting the stresses and (b) closure that is derived analytically in § 3.1. The black circle indicates the surface of the particle.

A very common simplification of the volume-filtered Euler–Lagrange framework is to use the following fluid momentum equation instead of (2.15) (Maxey & Patel Reference Maxey and Patel2001; Maxey Reference Maxey2017; Poustis et al. Reference Poustis, Senoner, Zuzio and Villedieu2019; Evrard et al. Reference Evrard, Denner and van Wachem2020):

(4.7)\begin{equation} \rho_{{f}}\dfrac{\partial u_{{pp},i}}{\partial t} + \rho_{{f}}\dfrac{\partial u_{{pp},i}u_{{pp},j}}{\partial x_j} ={-}\dfrac{\partial p_{{pp}}}{\partial x_i} + \mu_{{f}} \dfrac{\partial^2 u_{{pp},i}}{\partial x_j \partial x_j} - F_{{h},i}g(|\boldsymbol{x}-\boldsymbol{x}_p|), \end{equation}

where the index ${pp}$ indicates a quantity in the framework that we refer to as the simplified point-particle framework and $g$ is a Gaussian regularization kernel representing the fluid–particle momentum exchange. The simplification that $\epsilon _{{f}}=1$ is typically used for dilute regimes. We neglect the assumption of a constant fluid volume fraction $\epsilon _{{f}}$ and focus on the last term on the right-hand side of (4.7), the particle momentum exchange. We simplify our discussion by considering the Stokes flow around an isolated spherical particle

(4.8)\begin{equation} 0 ={-}\dfrac{\partial p_{{pp}}}{\partial x_i} + \mu_{{f}} \dfrac{\partial^2 u_{{pp},i}}{\partial x_j \partial x_j} - F_{{h},i}g(|\boldsymbol{x}-\boldsymbol{x}_p|). \end{equation}

Since the volume integrals of the Gaussian $g$ and the analytical kernel in the Stokes limit $\mathcal {K}_{{St}}$ equal one, the total momentum together with the particle is conserved. However, it has been argued in previous studies that the simplified momentum equation (3.1) is not energetically consistent if $\int u_{{pp},i}F_{{h},i}g(|\boldsymbol {x}-\boldsymbol {x}_p|)\,\mathrm {d}V\not = v_i F_{{h},i}$ (Xu & Subramaniam Reference Xu and Subramaniam2007; Fröhlich et al. Reference Fröhlich, Schneiders, Meinke and Schröder2018; Evrard et al. Reference Evrard, Denner and van Wachem2020), where energetical consistency is satisfied if viscous dissipation or inelastic collisions are the only mechanisms to reduce the total kinetic energy of the fluid–particle mixture. The simplified framework is indeed not energy consistent, but for a different reason. The reason for the inconsistency is that $\int u_{{pp},i}F_{{h},i}g(|\boldsymbol {x}-\boldsymbol {x}_p|)\,\mathrm {d}V\not =\int \bar {u}_is_i\,\mathrm {d}V$, i.e. a Gaussian is, strictly speaking, not the right kernel to regularize the particle momentum source.

Figure 6 shows the energy transfer rate of the particle momentum exchange term for the filtered momentum equation and the simplified momentum equation in Stokes flow. The results are shown for the moving sphere and the fixed sphere in the Stokes regime. The results of the simplified momentum equation are obtained analytically using the solution of a regularized Stokeslet (see, e.g. Cortez Reference Cortez2001).

Figure 6. Volume integrated energy transfer rates of the analytical particle momentum source and the particle momentum source in the simplified point-particle framework in Stokes flow. Black indicates configuration of the moving sphere in a quiescent fluid and red indicates the flow over a fixed sphere.

As the filter width approaches zero, the explicitly volume-filtered energy transfer terms converge towards $v_i F_{{h},i}$. In the simplified point-particle framework, the energy transfer term becomes infinitely large. Evrard et al. (Reference Evrard, Denner and van Wachem2020) show that, for a filter width of $\delta /a\approx 2{.}148$, the energy transfer equals the energy that is transferred to (or removed from) the particles $v_i F_{{h},i}$ in the Stokes regime if a Wendland kernel is employed to regularize the momentum source, which is a polynomial approximation of the Gaussian kernel. In contrast to what is claimed by Evrard et al. (Reference Evrard, Denner and van Wachem2020), $\delta /a\approx 2{.}148$ does not correspond to a consistent energy transfer because this would assume that the complete energy exchange with the particle only takes place at the filtered scales. In fact, the energy transfer is consistent if, and only if, the energy transfer term equals the the energy transfer rate from the explicit filtering. However, it is observed in figure 6 that the energy transfer of the simplified point-particle framework approaches the explicitly filtered solution with increasing filter width. This is supported by the finding that $\mathcal {K}_{{St}}\approx g$ for $\delta \gg 1$.

It is worth highlighting that instabilities in the simplified point-particle framework can occur when the filter width approaches zero, since the Gaussian approaches a singular Dirac delta distribution as it is observed by Chéron et al. (Reference Chéron, Brändle de Motta, Ménard, Poux and Berlemont2023a).

4.3. Single sphere at higher Reynolds number

We split the analysis of the single sphere in a flow with finite Reynolds number into an a priori analysis where the explicitly volume-filtered flow fields of particle-resolved DNS, i.e. a simulation that resolves the small flow structures around the particle, is investigated and an a posteriori analysis where the volume-filtered Navier–Stokes equations are solved with the proposed models for the closures and compared with the explicitly volume-filtered flow fields of particle-resolved DNS.

4.3.1. A priori analysis of the energy transfer

Similar to the studies in Stokes flow outlined earlier, we consider the flow around a single sphere at ${Re}=100$ in two frames of reference, a fixed sphere in a uniform flow and a moving sphere in a quiescent fluid. The purpose of the moving sphere configuration is to investigate the influence of the viscous closure, $\mathcal {E}_i$, at higher ${Re}$, which is zero for a fixed particle.

Figure 7 shows the volume integrated terms of the kinetic energy equation (4.3) for the fixed sphere configuration. The volume for the integration is a cube with an edge length of ten particle diameters. Note that $\mathcal {A}_i={\partial \epsilon _{{f}}\overline {u_i u_j}}/{\partial x_j}$ is an abbreviation for the explicitly volume-filtered advective term. Similar to the Stokes flow around an isolated spherical particle, the energy transfer rate of the particle momentum source is zero for $\delta =0$ and becomes increasingly negative with increasing filter width. The particle momentum source removes energy from the filtered scales and adds energy to the subfilter scales. For the larger ${Re}$ than in the Stokes limit, however, the energy transfer takes place at much smaller filter widths. This can be concluded because the slope of the energy transfer rate of the particle momentum source in figure 7 is very steep for small filter widths and flat for larger filter widths compared with the Stokes flow around the sphere. The energy transfer rate of the viscous term possesses the minimum at $\delta =0$ and approaches zero as the filter width increases. Consequently, the viscous dissipation predominates at the small scales. The energy transfer rate of the pressure term is positive for all filter widths. That means that at all filter widths, the pressure term adds energy to the filtered scales. The same is observed for the energy transfer rate of the advective term. However, it is observed that both energy transfer rates are not monotonically increasing or decreasing. The pressure term possesses a maximum and the advective term a minimum of energy transfer rate in the region of $1< \delta /a < 2$. Filtering out the small scale contribution of the energy transfer rate of the pressure term or the advective term in regions with positive slope of the pressure term or the advective term in figure 7 increases the respective energy transfer rate with respect to a smaller filter width. A negative slope of the pressure term or the advective term decreases the energy transfer rate with respect to a smaller filter width. Therefore, the pressure term and the advective term remove energy in regions of positive slope and add energy in regions of negative slope. The advective term removes energy from the large scales and adds energy to the small scales, which is qualitatively similar to what is known for single-phase turbulence (see, e.g. Pope Reference Pope2000). The results for ${Re}=20$ are qualitatively similar and, therefore, not shown.

Figure 7. Volume integrated energy transfer rates of the flow over a fixed sphere at $Re=100$. The energy transfer rates of the particle momentum source, the viscous term, the pressure term and the advective term are plotted.

4.3.2. A priori analysis of the viscous closure

In figure 8, the energy transfer rate of the viscous closure, $\mu _{{f}}\bar {u}_i\mathcal {E}_i$, is shown for the case of ${Re}=100$. Note that the frame of reference is shifted such that the particle moves with the velocity of $v_i$ in a quiescent fluid to obtain a non-zero viscous closure. Along with the energy transfer rate of the viscous closure, in figure 8, we show the energy transfer rate of the viscous term, as it occurs in (4.3), and the sum of both terms, which corresponds to the energy transfer rate of the explicitly volume-filtered viscous term in the NSE. In alignment with the observations in Stokes flow, the energy transfer rate of the viscous closure reaches a maximum at filter widths $\delta <1$ and approaches zero as the filter width increases. Consequently, the energy transfer rate of the viscous term deviates from the energy transfer rate of the explicitly filtered viscous term at small filter widths, but both approach each other for increasing filter widths. Although the viscous closure has negligible impact on the energy transfer for a large filter width in the investigated cases, common configurations exist where the viscous closure has a significant influence, even for large filter widths. In a channel flow, for instance, the fluid and the particles move with a large velocity relative to the frame of reference. The large particle velocity in the frame of reference increases the magnitude of the viscous closure compared with the present configuration. Since we provide an analytical expression for the viscous closure, its closure is recommended especially in such configurations.

Figure 8. Volume integrated energy transfer rates of a sphere moving in quiescent fluid at $Re=100$. The energy transfer rates of the viscous closure term, the viscous term and their sum are plotted.

4.3.3. A priori analysis of the subfilter stresses

For the study of the subfilter stresses, we consider the frame of reference where the particle is fixed, since the subfilter stress tensor is non-zero also in this configuration. In figure 9, the contributions of the subfilter stress tensor to the energy transfer are shown. In addition to the subfilter stress tensor, we also plot the energy transfer rate of the nonlinear model that is proposed in § 3.3, the Smagorinsky model adapted to a particle-laden flow and the subfilter stress tensor as it would be defined if the single volume fraction advective term $\tau _{{sfs},ij}^{{SVF}}$ would be used. Furthermore, the energy transfer rate of the explicitly filtered advective term is shown for comparison. In figure 9(a), the results of ${Re}=20$ are shown and in figure 9(b), the results of ${Re}=100$ are shown. For $\delta =0$ and both ${Re}$, the energy transfer rate of the subfilter stress tensor and the subfilter stress tensor with the single volume fraction definition of the advective term are zero. With increasing filter width, both terms yield increasingly negative energy transfer rates. With both definitions, $\tau _{{sfs},ij}$ and $\tau _{{sfs},ij}^{{SVF}}$, the subfilter stress tensor removes energy from the filtered scales and adds energy to the subfilter scales. This is analogous to the behaviour of the subfilter stress tensor in LES of single-phase turbulence, where the unresolved turbulent subfilter scales have a mainly dissipative effect on the filtered scales. This motivates the use of an additional (typically positive) turbulent viscosity as a model for the subfilter stress tensor. The energy transfer rates of both subfilter stress tensor definitions and both investigated ${Re}$ do not approach zero as the filter width increases. At the largest investigated filter width, $\delta /a=5$, the energy transfer rates of the subfilter stress tensors make up approximately 50 % of the energy transfer rate of the advective term for ${Re}=100$. For ${Re}=20$, however, both definitions of the subfilter stress tensors contribute only approximately 20 % to the energy transfer rate of the advective term. The significant contribution of the subfilter stress tensor to the energy transfer suggest that it requires modelling for the investigated configuration of an isolated sphere, even at larger filter widths. Note that even at a small energy transfer rate of the subfilter stress tensor, its contribution to the momentum can be significant. Figure 10 compares the contributions to the fluid momentum equation of the particle momentum source and the subfilter stress term in the streamwise direction for the flow over a fixed sphere at ${Re}=100$ and filter width $\delta /a=4$. The magnitude of the subfilter stress term is of similar magnitude and even exceeds the magnitude of the particle momentum source, which clearly highlights that the subfilter stress tensor requires modelling. In this particular configuration, neglecting the subfilter stress tensor would be as severe as omitting the particle momentum source.

Figure 9. Volume integrated energy transfer rates of the flow over a fixed sphere at (a) ${Re}=20$ and (b) ${Re}=100$. The energy transfer rates of the advective term, $\rho _{{f}}\bar {u}_i \mathcal {A}_i$, the subfilter stress tensor, $\rho _{{f}}\bar {u}_i \partial \tau _{{sfs},ij}/\partial x_j$, the subfilter stress tensor with the single volume fraction definition, $\rho _{{f}}\bar {u}_i \partial \tau _{{sfs},ij}^{{SVF}}/\partial x_j$, the modelled subfilter stress tensor with the nonlinear model, $\rho _{{f}}\bar {u}_i \partial \tau _{{sfs},ij}^{{mod,NL}}/\partial x_j$, and the subfilter stress tensor with the adapted Smagorinsky model, $\rho _{{f}}\bar {u}_i \partial \tau _{{sfs},ij}^{{mod,Sm}}/\partial x_j$, are plotted.

Figure 10. Contours of the momentum contributions in the streamwise direction of (a) the particle momentum source and (b) the divergence of the subfilter stress tensor for the fixed sphere at ${Re}=100$ and filter width $\delta /a=4$.

We also show the energy transfer rate of the nonlinear model for the subfilter stress tensor, which we propose in § 3.3, and the Smagorinsky model adapted to a particle-laden flow in figure 9. Up to a filter width of $\delta /a\approx 1$, the energy transfer rate of the nonlinear model almost perfectly agrees with the energy transfer rate of the subfilter stress tensor for ${Re}=20$ and ${Re}=100$. As the filter width increases further, however, the agreement deteriorates. The same behaviour is observed in single-phase turbulence (Borue & Orszag Reference Borue and Orszag1998; Johnson Reference Johnson2021). For larger filter widths, the neglected integral terms in (3.23) predominate. The Smagorinsky model yields an energy transfer rate of a too large magnitude for $\delta /a<4$. At larger filter widths, the magnitude of the energy transfer rate is too small. Although the Smagorinsky model contains a constant, which we assumed to be $C_{{s,p}}=1$ in the present case, no constant exists that leads to a matching energy transfer rate for a wide range of filter widths, which is clearly observed in figure 9 as the energy transfer rate by the Smagorinsky model behaves differently with respect to the filter width than the energy transfer rate of the subfilter stress tensor. Therefore, the Smagorinsky model can generally not be recommended as a model for the subfilter stress tensor; at least not for the flow around an isolated sphere. A potential improvement of the Smagorinsky model could be to vary $C_{{s,p}}$ as a function of the filter width, which could be achieved with a dynamic procedure.

Figure 11 shows contours of the three components of the divergence of the subfilter stress tensor for the subfilter stress tensor and the nonlinear model proposed in § 3.3 for a filter width of $\delta /a=1$ and ${Re}=100$. The results are normalized with $f_{{ref}}=|\boldsymbol {F}_{{h}}|/V_{{p}}$. The shape of the contours of the modelled subfilter stress tensor show similar structures as the contours of the subfilter stress tensor. The first component, which is the streamwise component, possesses the largest magnitude. The contours are very similar in shape but the magnitudes deviate. The second component possesses the smallest magnitude. Although the overall shape of the contours is similar, minor deviations in shape occur near the upstream surface of the sphere. The third component agrees well in shape but shows the same discrepancies in magnitude as the other components. For smaller filter widths, the agreement in shape and magnitude improves, but for larger filter widths, the modelled subfilter stress tensor shows larger deviations, which supports the observations of the energy transfer rates.

Figure 11. Divergence of the subfilter stress tensor of a flow around a fixed sphere at ${Re}=100$. The plots (a,c,e) are obtained from explicitly volume filtering the results of the particle-resolved simulation and the plots (b,d,f) are computed with the nonlinear model. The filter width of the shown case is $\delta /a=1$.

4.3.4. A priori analysis of the particle momentum source

In figure 12, the spatial distribution of the particle momentum source in a flow of ${Re}=100$ is compared with a Gaussian and the analytical kernel in Stokes flow $\mathcal {K}_{{St}}$ for a filter width $\delta /a=1$. The spatial distribution of the particle momentum source appears as a ring in the two-dimensional slice, and is a spherical hull in three dimensions. In contrast to the spherical symmetry of the kernel in the Stokes limit, an upstream–downstream asymmetry is observed for larger ${Re}$. Near the stagnation point, the greatest magnitude of the kernel is observed. In figure 12(b), the kernel is plotted along lines through the centre of the sphere in the streamwise and normal directions together with the analytical kernel in the Stokes limit and the commonly used Gaussian regularization of the same filter width. The magnitude of the kernel near the stagnation point is approximately six times as large as the magnitude on the other side of the sphere. In the normal direction, the kernel is symmetric with a maximum magnitude slightly larger than the magnitude downstream. The analytical kernel in the Stokes limit possesses spherical symmetry and a maximum magnitude of approximately one third of the actual kernel near the stagnation point. However, the positions of the maxima are approximately similar to those of the actual kernel for ${Re}=100$. The maximum of the Gaussian is in the centre of the sphere and has a much too large magnitude. It can be concluded that at such small filter widths, the Gaussian is a poor approximation of the regularization of the particle momentum source in the Stokes limit and at larger ${Re}$. Note that at finite ${Re}$, the particle momentum source does not only have a contribution in the streamwise direction, but also in the normal, or perpendicular, directions. The normal momentum source components have a much smaller magnitude than the streamwise component and their integrals over the entire volume are zero, i.e. there is no resultant momentum source in the normal directions.

Figure 12. Shape of the particle momentum source for a flow over a fixed sphere at ${Re}=100$ and filter width $\delta /a=1$. (a) Contour of the regularization kernel. (b) Regularization kernel in the streamwise and normal direction together with the analytical regularization kernel for Stokes flow and a Gaussian of the same filter width.

Figure 13 compares the shapes of the regularization of the particle momentum source at $\delta /a=4$. The shape of the kernel is relatively close to a Gaussian that is shifted slightly upstream. Except for this upstream shift, $\mathcal {K}_{{St}}$ approximates the shape and the magnitude of the kernel well. Although the Gaussian is much closer to the kernel than for $\delta /a=1$, it still overpredicts the magnitude of the kernel by a factor of approximately $1.5$.

Figure 13. Shape of the particle momentum source for a flow over a fixed sphere at ${Re}=100$ and filter width $\delta /a=4$. (a) Contour of the regularization kernel. (b) Regularization kernel in the streamwise and normal direction together with the analytical regularization kernel for Stokes flow and a Gaussian of the same filter width.

The studies of the flow relative to a single isolated sphere show that all of the three closures, the viscous closure, the subfilter stress tensor and the regularization kernel of the particle momentum source, can significantly contribute to the energy transfer and require appropriate modelling in the investigated case. The viscous closure can be obtained analytically and is shown to match the explicit evaluation of $\mathcal {E}_i$. With the nonlinear model, the subfilter stresses can be approximated well for filter widths up to approximately the particle size, but the Smagorinsky model is shown to not be suitable for modelling the subfilter stress tensor. For larger filter widths, $\mathcal {K}_{{St}}$ can recover the shape of the regularization of the particle momentum source well. For filter widths $\delta /a>20$, which are relevant for simulations using the PSIC method, conclusions cannot be drawn with certainty from the present results. The observed trends suggest, however, that the shape of the regularization kernel approaches a very wide Gaussian and that the subfilter stress tensor still contributes significantly to the energy and momentum transfer, such that it requires modelling.

4.3.5. A posteriori analysis

When the volume-filtered equations are solved in a reduced order simulation, such as an Euler–Lagrange or an Euler–Euler simulation, it is particularly relevant what the impact of the closures is on the resulting fluid velocity field. Therefore, we solve the volume-filtered equations with different models for the closures and compare it with the explicitly filtered flow fields from particle-resolved simulations. To isolate the effect of the closure models on the resulting fluid velocity field, the volume-filtered equations are solved on a very fine mesh, i.e. 15 cells per particle diameter, which is much finer than any expected application of volume filtering. However, this high resolution minimizes the numerical error and the deviations can be traced back to the modelling errors. The simulation set-up is similar to the interphase-resolving simulation of the isolated sphere.

In figure 14, we show the contours of the volume-filtered streamwise velocity of the flow past a fixed sphere of ${Re}=20$ and ${Re}=100$ for three different cases: the explicitly volume-filtered particle-resolved velocity, the Gaussian regularization of the particle momentum source without a model for the subfilter stress tensor, and the regularization with the analytical kernel in the Stokes limit with the nonlinear model for the subfilter stress tensor. The results are shown for a filter width of $\delta /a=2$. The case with the Gaussian regularization corresponds to the simplified point-particle framework, which is introduced in § 4.2. Note that the viscous closure is zero in the present case because the particles are not moving in the considered frame of reference.

Figure 14. Contours of the filtered streamwise fluid velocity of a flow over a fixed sphere at (a,c,e) ${Re}=20$ and (b,d,f) ${Re}=100$. Compared are (a,b) the explicitly volume-filtered resolved simulation, (c,d) the simplified point-particle framework with a Gaussian regularization and without a model for the subfilter stress tensor, and (e,f) the volume-filtering framework with the analytical kernel of the Stokes regime and the nonlinear model for the subfilter stress tensor. The investigated case corresponds to a filter width $\delta /a=2$.

In figure 14, it can be observed for both ${Re}$ that the explicitly volume-filtered fluid velocity possesses a small positive minimal velocity slightly downstream from the particle centre. For ${Re}=100$, the minimal fluid velocity is further downstream than for ${Re}=20$. With the simplified point-particle framework, where a Gaussian is used for regularization of the particle momentum source, a large negative (upstream) velocity is induced for both ${Re}$. At the considered filter width, the Gaussian is much too localized or concentrated in the particle centre. For ${Re}=20$, the minimum streamwise velocity is present at the particle centre and for ${Re}=100$, even upstream of the particle. A large spurious velocity is induced that reduces the accuracy of the undisturbed velocity prediction, which is required for the empirical correlations for the drag force and other forces on the particle. A variety of different methods have been developed to subtract this partially spurious self-induced fluid velocity (Ireland & Desjardins Reference Ireland and Desjardins2017; Balachandar et al. Reference Balachandar, Liu and Lakhote2019; Evrard et al. Reference Evrard, Denner and van Wachem2020; Pakseresht & Apte Reference Pakseresht and Apte2021; Balachandar & Liu Reference Balachandar and Liu2022). With the analytical regularization kernel of the Stokes regime and the nonlinear model for the subfilter stress tensor, the minimal velocities are larger positive (downstream) values than for the explicitly volume-filtered velocity and the minima are too far downstream. However, the velocity magnitudes are much closer to the explicitly filtered solution than the velocity of the Gaussian regularization. At ${Re}=20$, a good approximation of the explicitly filtered fluid velocity is obtained. For ${Re}=100$, we show in § 4.3.4 that the shape of the particle momentum source is strongly asymmetric. The assumption of a spherical symmetric kernel of the Stokes regime becomes decreasingly valid for increasing ${Re}$.

It is worth highlighting that even at the small filter width of $\delta /a=2$, a stable solution of the volume-filtered momentum equation (2.15) together with the volume-filtered continuity equation (2.9) is obtained. With the single volume fraction advective term, even at larger filter widths, numerical interventions are typically required (Link et al. Reference Link, Cuypers, Deen and Kuipers2005; Jing et al. Reference Jing, Kwok, Leung and Sobral2016).

The results of the same configurations, but with a filter width $\delta /a=4$, are shown in figure 15. At this filter width, the contribution of the subfilter stress tensor is significant and the nonlinear model for the subfilter stress tensor predicts too small values. It is shown in § 4.3.4 that $\mathcal {K}_{{St}}$ is a good approximation of the shape of the particle momentum source, even at ${Re}=100$, and that the Gaussian is still too localized and possesses a too large magnitude at $\delta /a=4$.

Figure 15. Contours of the filtered streamwise fluid velocity of a flow over a fixed sphere at (a,c,e) ${Re}=20$ and (b,d,f) ${Re}=100$. Compared are (a,b) the explicitly volume-filtered resolved simulation, (c,d) the simplified point-particle framework with a Gaussian regularization and without a model for the subfilter stress tensor, and (e,f) the volume-filtering framework with the analytical kernel of the Stokes regime and the nonlinear model for the subfilter stress tensor. The investigated case corresponds to a filter width $\delta /a=4$.

The streamwise fluid velocities in figure 15 show fewer deviations than for $\delta /a=2$. The Gaussian regularization captures the magnitude of the minimal velocity for ${Re}=20$ well and slightly overestimates it for ${Re}=100$. With $\mathcal {K}_{{St}}$ and the nonlinear model, the minimal streamwise velocity is overestimated more than with the Gaussian regularization for both ${Re}$. Although $\mathcal {K}_{{St}}$ is a better approximation of the shape of the particle momentum source than the Gaussian, the streamwise fluid velocity field is in better agreement with the explicitly filtered velocity. The reason is that the neglected subfilter stress tensor is at least partially compensated by the overestimation of the regularization kernel. This compensation of errors causes the better agreement when a Gaussian is used at filter width $\delta /a=4$. Instead of relying on compensation of errors, however, the more physical kernel $\mathcal {K}_{{St}}$ should be used and a model for the subfilter stress tensor that is also valid at large filter widths needs to be derived.

4.4. Periodic array of spheres

To investigate the presence of multiple spheres on the closures, a periodic array of fixed spheres is considered. The flow is driven by a constant momentum source.

Figure 16 shows the volume integrated energy transfer rates for the fluid flow through the periodic array of spheres at ${Re}_{{periodic}}=50$. Since the energy is introduced by the additional driving momentum source, the energy transfer rate of the pressure term is zero at $\delta =0$. The conclusions on the energy transfer are similar to the single sphere configurations, with the exception of the advective term. Because of the periodic domain and the resulting repetition of the flow patterns, all existing flow structures are captured by integrating over a finite volume, whereas in the single sphere configuration, a finite integration volume always omits some of the existing flow patterns. The behaviour of the advective term with increasing filter width is very similar to the single sphere configuration, but at $\delta =0$, the energy transfer rate of the advective term in the periodic case is zero. In the periodic configuration, the advective term removes energy from large scales and adds it to small scales. However, the contribution of the advective term to the energy transfer is relatively small for all filter widths.

Figure 16. Volume integrated energy transfer rates of the flow over a periodic array of spheres at ${Re_{periodic}=50}$. The energy transfer rates of the particle momentum source, the viscous term, the pressure term and the advective term are plotted.

In figure 17, the energy transfer rates of the subfilter stress tensor, the single volume fraction definition of the subfilter stress tensor and the modelled subfilter stress tensor using the nonlinear model are plotted. The energy transfer rate of the advective term is shown as a reference. The energy transfer rate of the present definition and the single volume fraction definition of the subfilter stress tensor are relatively similar and deviate only slightly from each other. The energy transfer rate contribution of the subfilter stress tensor is significant compared with the energy transfer caused by the advective term, especially at large filter widths. Similar to the single sphere cases, the nonlinear model for the subfilter stress tensor approximates the energy transfer rate of the subfilter stress tensor well for small filter widths but increasingly deviates from the energy transfer rate of the subfilter stress tensor as the filter width increases.

Figure 17. Volume integrated energy transfer rates of the flow over a periodic array of spheres at ${Re_{periodic}=50}$. The energy transfer rates of the advective term, $\rho _{{f}}\bar {u}_i \mathcal {A}_i$, the subfilter stress tensor, $\rho _{{f}}\bar {u}_i \partial \tau _{{sfs},ij}/\partial x_j$, the subfilter stress tensor with the single volume fraction definition, $\rho _{{f}}\bar {u}_i \partial \tau _{{sfs},ij}^{{SVF}}/\partial x_j$, and the modelled subfilter stress tensor with the nonlinear model, $\rho _{{f}}\bar {u}_i \partial \tau _{{sfs},ij}^{{mod,NL}}/\partial x_j$, are plotted.

Although the energy transfer rate of the subfilter stress tensor is small compared with the energy transfer rate of the particle momentum source, the subfilter stress tensor can still significantly contribute to the momentum. Figure 18 shows the particle momentum source and the subfilter stress momentum contributions in the streamwise direction for a filter width of $\delta /a=4$. At this filter width, the energy transfer rate of the subfilter stresses is almost zero compared with the energy transfer rate of the particle momentum source. The magnitude of the momentum contribution of the subfilter stresses is approximately one third of the magnitude of the particle momentum source. The significant relative contribution of the subfilter stresses to the fluid momentum confirms that they require modelling even at such large filter widths.

Figure 18. (a) Contours of the momentum contributions in the streamwise direction of the particle momentum source and (b) the divergence of the subfilter stress tensor for a periodic array of spheres at ${Re}_{{periodic}}=50$ and filter width $\delta /a=4$.

As already discussed in the single sphere configuration, another error source in the volume-filtered momentum equation is the regularization of the particle momentum source. For the periodic array of spheres at ${Re}_{{periodic}}=50$, the regularization kernel is shown for $\delta /a=1$ in figure 19 together with the analytical kernel of the Stokes regime and a Gaussian. Along the streamwise direction, the kernel is asymmetric. The largest value of the kernel is at the stagnation point of the sphere. Downstream of the particle, the kernel is even negative. In the normal direction, the kernel is approximated well by the analytical regularization kernel of the Stokes flow around an isolated sphere, $\mathcal {K}_{{St}}$, but the Gaussian of the same filter width possesses a much too large magnitude and the wrong shape.

Figure 19. Shape of the particle momentum source for a flow over a fixed array of spheres in the periodic configuration at ${Re}_{{periodic}}=50$ and filter width $\delta /a=1$. (a) Contour of the regularization kernel. (b) Regularization kernel in the streamwise and normal direction together with the analytical regularization kernel for Stokes flow and a Gaussian of the same filter width.

The regularization kernel is plotted for $\delta /a=4$ in figure 20. The shape of the kernel is relatively similar to a Gaussian in the streamwise and the normal directions but slightly shifted upstream. Due to the periodic domain, the kernel is always significantly larger than zero. At this larger filter width, $\mathcal {K}_{{St}}$ is a good approximation of the kernel in terms of shape and magnitude. However, the Gaussian overpredicts the magnitude of the kernel.

Figure 20. Shape of the particle momentum source for a flow over a fixed array of spheres in the periodic configuration at ${Re}_{{periodic}}=50$ and filter width $\delta /a=4$. (a) Contour of the regularization kernel. (b) Regularization kernel in the streamwise and normal directions together with the analytical regularization kernel for Stokes flow and a Gaussian of the same filter width.

The main conclusions of the single sphere configuration also apply to the periodic array of spheres. The subfilter stress tensor requires modelling, because it contributes significantly in the momentum equation of the fluid relative to other momentum sources, such as the particle momentum source. For small filter widths, the nonlinear model captures the energetic influence of the subfilter stresses. The shape of the regularization of the particle momentum source is poorly approximated by the commonly used Gaussian. At larger filter widths that are of practical importance for Euler–Lagrange simulations, $\mathcal {K}_{{St}}$ approximates the kernel well and is to be preferred over a Gaussian.

5.1. Choosing an appropriate filter width

The volume-filtered NSE, including the derived closures, can be solved in the scope of a point-particle simulation, where the influence of the particles on the fluid is considered by the regularized particle-momentum source. Point-particle simulations are typically carried out on a computational fluid mesh that is coarse compared with the size of the particles and, in addition to the modelling errors discussed in the previous sections, discretization errors occur.

There are three length scales in point-particle simulation that have to be considered, the particle diameter, $d_{{p}}$, the size of the computational fluid mesh cells, $\Delta x$, and the filter width, $\delta$. Typically, the particle diameter is predetermined and the remaining parameters have to be chosen accordingly. It is demonstrated in the previous sections that a Gaussian kernel approximates the particle momentum source well, even for Reynolds numbers up to ${Re}=100$, if the filter width is large enough. In Stokes flow, a filter width of $\delta /d_{{p}} \gtrapprox 3.5$ yields an analytical kernel that deviates insignificantly from a Gaussian. A Gaussian kernel is advantageous, as it can be easily integrated over a cell of the fluid mesh. As shown by Balachandar & Liu (Reference Balachandar and Liu2022), at large filter widths, the following approximation becomes valid:

(5.1)\begin{equation} \epsilon_{{p}}(\boldsymbol{x}) = \int_{\varOmega} I_{{p}}(\boldsymbol{y})g(|\boldsymbol{x}-\boldsymbol{y}|)\,\mathrm{d}V_y \approx V_{{p}}g(|\boldsymbol{x}|), \quad \sigma\gg d_{{p}}, \end{equation}

where $I_{{p}}$ is the particle indicator function and $\sigma$ is the standard deviation of the Gaussian. For $\sigma /d_{{p}}=1$, which corresponds to $\delta /d_{{p}}\approx 3.76$, this approximation leads to a maximum deviation of $7.7\,\%$ from the analytical volume fraction. To justify the Gaussian approximation of the particle momentum source and the volume fraction, we suggest to choose the filter width of at least $\delta /d_{{p}}>3.76$ or $\sigma /d_{{p}}>1$ independent of the fluid mesh size.

The filter width cannot be chosen independently of the mesh size. The purpose of solving the volume-filtered NSE instead of the NSE is that large gradients that cannot be resolved by a coarse fluid mesh are avoided. In the volume-filtered NSE, large gradients are mainly avoided by choosing a widely distributed particle momentum source. If the size of the fluid mesh cell is too coarse for the chosen filter width, the discretization error imposes a spatially varying unknown filter, and the closures, which are derived for a spatially constant Gaussian filter of a specific filter width, are not valid. If the filter width is chosen as $\sigma /\Delta x \gtrapprox 1$, a Gaussian can be well resolved by the mesh. Sampling a Gaussian on a mesh with a significantly smaller ratio of $\sigma /\Delta x$ may capture only a fraction of its energy.

For point-particle simulations, we recommend to choose the filter width based on what is larger, the particle diameter, $d_{{p}}$, or the fluid mesh size, $\Delta x$. The filter width should satisfy both $\sigma /d_{{p}} \gtrapprox 1$ and $\sigma /\Delta x \gtrapprox 1$ to justify the Gaussian approximation of the volume fraction and the particle momentum source, and to avoid the solution being altered by the filter imposed by a too coarse fluid mesh. Note that point-particle simulations with filter widths larger than the grid spacing may become computationally expensive, which is competing with the criterion $\sigma /\Delta x \gtrapprox 1$. In practice, a compromise must be found between sufficient resolution and computational expenses.

The guidelines also hold in a turbulent flow independent of the size of the Kolmogorov length scale. In an LES, even without particles, $\sigma /\Delta x \gtrapprox 1$ should be satisfied to avoid flow structures that are too small to be resolved by the mesh.

5.2. Implementation of the closures

In the volume-filtered continuity equation (3.5), in the Galilean invariance part of the subfilter stress tensor $\tau _{{sfs},ij}^{{G}}$ and in the viscous closure $\mathcal {E}_i$, spatial gradients of the volume fraction are required. In a finite volume framework, the spatial gradients are required as integrals over a fluid mesh cell. If the volume fraction can be approximated with a Gaussian, the following integrals over a fluid mesh cell $\varOmega _{\varDelta }$ can be computed analytically as

(5.2)$$\begin{gather} \int_{\varOmega_{\Delta}}\dfrac{\partial \epsilon_{{p}}}{\partial x_i}\,\mathrm{d}V \approx V_{{p}} \int_{\varOmega_{\Delta}}\dfrac{\partial g}{\partial x_i}\,\mathrm{d}V, \end{gather}$$

(5.3)$$\begin{gather}\int_{\varOmega_{\Delta}}\dfrac{\partial^2 \epsilon_{{p}}}{\partial x_i\partial x_i}\,\mathrm{d}V \approx V_{{p}} \int_{\varOmega_{\Delta}}\dfrac{\partial^2 g}{\partial x_i\partial x_i}\,\mathrm{d}V, \end{gather}$$

(5.4)$$\begin{gather}\int_{\varOmega_{\Delta}}\dfrac{\partial \epsilon_{{p}}^2}{\partial x_i}\,\mathrm{d}V \approx V_{{p}}^2 \int_{\varOmega_{\Delta}}\dfrac{\partial g^2}{\partial x_i}\,\mathrm{d}V. \end{gather}$$

Note that with these analytical integrals, there is no significant difference in computational costs between the discrete computation of the spatial gradient of the sum of all particle volume fractions, i.e.

(5.5)\begin{equation} \dfrac{\partial }{\partial x_i}\sum_q \int_{\varOmega_{\Delta}}\epsilon_{{p},q}\,\mathrm{d}V, \end{equation}

and the sum of the analytical gradient of the particle volume fractions, i.e.

(5.6)\begin{equation} \sum_q \int_{\varOmega_{\Delta}}\dfrac{\partial \epsilon_{{p},q}}{\partial x_i}\,\mathrm{d}V. \end{equation}

However, the former causes a discretization error and the latter is exact up to machine precision.

An accurate estimate of the spatial gradients of the volume fraction in the continuity equation, in the closure $\tau _{{sfs},ij}^{{G}}$ and in the viscous closure $\mathcal {E}_i$ is crucial to achieve Galilean invariance in the numerical simulations. Galilean invariance is satisfied if the following relation for the volume-filtered velocity in the new frame of reference $\epsilon _{{f}}\bar {\mathcal {U}}_i$ holds:

(5.7)\begin{align} \epsilon_{{f}}\bar{\mathcal{U}}_i(\boldsymbol{x},t) &= \int_{\varOmega} I_{{f}}(\boldsymbol{y}) (u_i(\boldsymbol{y}-\boldsymbol{u}_{{ref},i}t,t) + u_{{ref},i})g(|\boldsymbol{x}-\boldsymbol{y}|)\,\mathrm{d}V_y \nonumber\\ &= \epsilon_{{f}}\bar{u}_i(\boldsymbol{x}-\boldsymbol{u}_{{ref},i}t,t) + \epsilon_{{f}}u_{{ref},i}(\boldsymbol{x},t). \end{align}

To isolate the effect of determining the terms involving spatial gradients of the volume fraction, a fixed particle in a quiescent fluid is considered. The particle-momentum source is zero in this case. Transforming the frame of reference with a velocity $u_{{ref}}$ yields a fluid velocity and a particle velocity equal to $u_{{ref}}$ in the new frame. The volume-filtered fluid velocity is given according to (5.7). Any deviation from $\epsilon _{{f}}\bar {\mathcal {U}}_i$ of the volume-filtered velocity in the new frame is considered a spurious current. Without the terms $\tau _{{sfs},ij}^{{G}}$ and $\mathcal {E}_i$ in the volume-filtered momentum equation, spurious currents occur because $\epsilon _{{f}}\bar {u}_i$ is not divergence-free. Based on the velocity of the moving frame of reference, we consider two Reynolds numbers, ${Re}_{{ref}}=u_{{ref}}d_{{p}}\rho _{{f}}/\mu _{{f}}$. For ${Re}_{{ref}}=100$, the closure $\tau _{{sfs},ij}^{{G}}$ is dominant and for ${Re}_{{ref}}=1$, the viscous closure $\mathcal {E}_i$ dominates. Figure 21 shows the spurious currents when the frame of reference is transformed using the proposed closures computed as described and without the closures and a resolution of $d_{{p}}/\Delta x = 2$. The filter width is chosen according to the guidelines $\sigma /d_{{p}}=1$. Without the closures $\tau _{{sfs},ij}^{{G}}$ and $\mathcal {E}_i$, a spurious velocity of up to $8\,\%$ of $u_{{ref}}$ occurs. With the proposed discretization of the volume fraction gradients, the closures drastically reduce the spurious currents at this resolution. Note that the closures are exact and solely the discretization causes the observed small spurious currents. Figure 22 shows the spurious currents at a coarser resolution of $d_{{p}}/\Delta x = 0.5$ and a filter width $\sigma /\Delta x = 1$, according to the guidelines. The maximum magnitude of the spurious velocity is relatively small even without closures. The closures, however, reduce the spurious velocity even further. The spurious currents reduce when the filter width is increased but they can reach a significant magnitude even if the filter width is large. In configurations with a large mean flow velocity and relatively small relative velocity between fluid and particles, such as channel flows, $u_{{ref}}$ is large and can induce spurious currents that are large compared with the velocity induced by the particle momentum source.

Figure 21. Contours of the spurious streamwise fluid velocity in a configuration without relative velocity between fluid and particle in a frame of velocity $u_{{ref}}$. The spurious currents resulting from the volume-filtered NSE with and without the closures $\tau _{{sfs},ij}^{{G}}$ and $\mathcal {E}_i$ are compared for ${Re}_{{ref}}=100$ and ${Re}_{{ref}}=1$ with a resolution $d_{{p}}/\Delta x = 2$ and a filter width $\sigma /d_{{p}}=1$.

Figure 22. Contours of the spurious streamwise fluid velocity in a configuration without relative velocity between fluid and particle in a frame of velocity $u_{{ref}}$. The spurious currents resulting from the volume-filtered NSE with and without the closures $\tau _{{sfs},ij}^{{G}}$ and $\mathcal {E}_i$ are compared for ${Re}_{{ref}}=100$ and ${Re}_{{ref}}=1$ with a resolution $d_{{p}}/\Delta x = 0.5$ and a filter width $\sigma /d_{{p}}=2$.

The present test case demonstrates that an additional viscosity, as applied by Zhang & Prosperetti (Reference Zhang and Prosperetti1997), Patankar & Joseph (Reference Patankar and Joseph2001) and Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013), is not suitable to replace the viscous closure, as it solely decreases ${Re}_{{ref}}$, which does not avoid the spurious currents.

It is shown in the previous sections that the subfilter stress tensor can significantly contribute to the energy transfer and the fluid momentum balance for small and large filter widths and, therefore, requires appropriate modelling. The proposed nonlinear model for the subfilter stress tensor predicts excellent agreement with the subfilter stress tensor for filter widths up to approximately the particle size, but for large filter widths, large deviations are observed. A Smagorinsky type model for the subfilter stress tensor cannot predict the energy transfer of the subfilter stress tensor, independent of the model constant. For filter widths according to the guidelines, neither of the discussed models for $\tau _{{sfs},ij}$ can be recommended for point-particle simulation and further research is required to derive models that are valid for large filter widths.

The closures that are discussed in the present paper, including their expressions, their implementation and the conditions under which they are valid, are summarized in table 1.

Table 1. Summary of the closures discussed in the present paper, how they can be implemented in point-particle simulations and under which conditions they are valid.

5.3. Remark on the PSIC method

With the very popular PSIC method (Crowe et al. Reference Crowe, Sharma and Stock1977), the particle momentum source is applied to the computational cell in which the particle lies, which may satisfy $\sigma /d_{{p}} \gtrapprox 1$ for small particles, but always violates $\sigma /\Delta x \gtrapprox 1$. Therefore, the PSIC method is not in accordance with the guidelines for the filter width. This can have at least the following potentially problematic consequences. (i) If particles move from one fluid mesh cell to another, sudden jumps in the particle momentum source may occur. (ii) A fluid mesh cell may contain significantly more particles than a neighbouring mesh cell, which can cause large gradients of the particle momentum source and, consequently, large fluid velocity gradients that cannot be resolved by the fluid mesh. The filter cannot be assumed to be a spatially uniform Gaussian and additional unknown closures arise. (iii) The PSIC method proposed by Crowe et al. (Reference Crowe, Sharma and Stock1977) leads to a regularization of the particle momentum source and, consequently, an induced local fluid velocity that depends on the ratio $d_{{p}}/\Delta x$. In point-particle simulations, the local fluid velocity interpolated to the particle position is used to compute the drag force on the particle. As shown by Evrard, Denner & van Wachem (Reference Evrard, Denner and van Wachem2021), the PSIC method can cause significant local fluid velocity disturbances even for $d_{{p}}/\Delta x\approx 0.1$, which leads to an erroneous prediction of the forces acting on the particle. (iv) Since the local filter width implied by the PSIC method depends on the local mesh resolution, local mesh refinement, e.g. near a solid wall, leads to a non-uniform filter and additional closures that are neglected.

Eaton (Reference Eaton2009) summarizes commonly observed discrepancies between simulations using the PSIC method and experiments, which are potentially caused by these problems. Regularizing the particle momentum source over multiple cells using a Gaussian with a sufficiently large spatially uniform standard deviation can mitigate the mentioned problems of the PSIC method. However, the complexity of the implementation is increased and some computational overhead is added.