Modelling of filtered drag force for clustered particle-laden flows based on interface-resolved simulation data

Abstract

Interface-resolved direct numerical simulations (DNS) of clustered settling suspensions in a periodic domain are performed to study the filtered drag force for clustered particle-laden flows. Our results show that, for the homogeneous system, the filtered drag is independent of the filter size, whereas for the clustered particle-laden flows, the averaged drag becomes smaller than the homogeneous drag at the filter size above 4 particle diameters. The drag reduction saturates at the filter size being comparable to the cluster size in the horizontal direction in our simulations. A new correlation is proposed to account for the mesoscale effect on the filtered drag force by using drift velocity and variance of the solid volume fraction, based on the modification of existing subgrid drag models for the inhomogeneous system. The existing models for the drift velocity and the variance of the solid volume fraction are assessed using our DNS data. A new model for the drift velocity and the variance of the solid volume fraction is proposed, based on the combination and modification of the previous models. All mesoscale models considered can predict well the filtered drag with comparable accuracy, and are superior to the homogeneous drag model for the clustered system. Our models with the same parameter values obtained from the large-scale system can also predict well the filtered drag for smaller computational domain sizes.

1. Introduction

Particle-laden flows are ubiquitous in numerous natural and industrial processes, such as fluidized beds, drug delivery, drilling processes and sand storms. A great many works have been dedicated to the development of computational models for particle-laden flows. The two-fluid model (i.e. the Euler–Euler method) (Ding & Gidaspow 1990; Ishii & Hibiki 2010) and the discrete particle (or element) model (i.e. the Euler–Lagrange method) (Tsuji, Kawaguchi & Tanaka 1993; Deen et al. 2007) are two models widely used for simulations of industrial-scale multiphase flows. The drag model is most important for the accuracy of these two methods and is the subject of extensive studies (Tenneti & Subramaniam 2014; Wang 2020). The feature of the two-fluid model is that the particle phase is treated as a continuous phase, like a fluid, and a computational cell could contain a large number of particles whose positions are not tracked. In the present study, we are concerned with modelling of the average drag used in the two-fluid model for dense particle-laden flows, and thus, in the following we only review the works on the drag models for two-fluid simulations, although there are many other state-of-the-art methods for the particle-laden flows based on filtering approaches, such as the multiphase turbulence models (Fox 2014) and the Lagrangian probability-density-function model (Innocenti et al. 2019; Innocenti, Fox & Chibbaro 2021).

The early and widely used drag models were obtained from experimental data (Wen & Yu 1966; Gidaspow 1994; Zhang & Reese 2003). In the past two decades, interface-resolved direct numerical simulations (IR-DNS) (van der Hoef et al. 2008; Balachandar & Eaton 2010; Tenneti & Subramaniam 2014; Maxey 2017) have been employed extensively to develop more accurate drag models based on either fixed particle beds (Hill, Koch & Ladd 2001; van der Hoef, Beetstra & Kuipers 2005; Beetstra, van der Hoef & Kuipers 2007; Tenneti, Garg & Subramaniam 2011; Tang et al. 2015; Zhou & Fan 2015) or dynamic particle-laden flows (Luo et al. 2016; Rubinstein, Derksen & Sundaresan 2016; Tang, Peters & Kuipers 2016; Tavanashad, Passalacqua & Subramaniam 2021; Xia et al. 2022b).

The above-mentioned models were mainly developed for the homogeneous system where the particles are distributed homogeneously within the computational cell. In many industrial-scale applications, particles with random and uniform distributions initially may self-organize into heterogeneous structures (clusters) with a wide range of spatial and temporal scales. The heterogeneous structure of particles is also referred to as the mesoscale structure. When there exists a particle cluster, namely the particle distribution is inhomogeneous in a computational cell in the coarse-grid simulations, the drag models for the homogeneous system are no longer accurate, and often overestimate the drag significantly (Sundaresan, Ozel & Kolehmainen 2018). The filtered method (Igci et al. 2008) and the heterogeneity-based method (including the energy minimization multi-scale model, i.e. EMMS model) (Li 1994) are two typical methods accounting for the effects of unresolved mesoscale structure on the drag force (Wang 2020).

The basic idea of the filtered model is to establish the drag model based on the data from highly resolved simulations, including the fine-grid two-fluid model (TFM) (Igci et al. 2008; Igci & Sundaresan 2011), discrete element method–computational fluid dynamics (DEM-CFD) (Radl & Sundaresan 2014; Ozel et al. 2017; Beetham, Fox & Capecelatro 2021) and interface-resolved simulations (Rubinstein et al. 2017; Liu, Ge & Wang 2020) with a filtering (averaging) procedure. There are two categories for the filtered drag models,depending on whether the drift velocity is used for the modelling. In the first category, drag models are modified by a mesoscale correction coefficient directly using available filtered quantities such as filtered volume fraction, filtered slip velocity and filter size (Milioli et al. 2013; Radl & Sundaresan 2014; Sarkar et al. 2016; Zhu, Liu & Luo 2018; Cloete et al. 2019; Milioli & Milioli 2023).

In the second category, the effects of the mesoscale structures on the drag are correlated with the subgrid quantities, including the drift velocity and the variance of the solid volume fraction, which also need to be modelled with the filtered quantities. The drift velocity is defined as the difference between the mean fluid velocity seen by the particles and the phase-averaged fluid velocity in the cell. Parmentier, Simonin & Delsart (2012) found that it is necessary to take the subgrid drift velocity into account to predict the significant drag decrease obtained from fine-grid TFM results. The variance of the solid volume fraction, which indicates the heterogeneity of particle distribution in the filter cell, has been found to be related with the drag force (Schneiderbauer 2017). Chen et al. (2020) proposed an expression of the filtered drag force using the Taylor polynomial approximation of a microscopic drag coefficient, which reveals that the drift velocity and the variance of the solid volume fraction are two important subgrid quantities for the estimation of the drag force. The DEM-CFD simulations (Ozel et al. 2017) and interface-resolved simulations (Rubinstein et al. 2017; Zhao, Chen & Zhou 2021) demonstrated that the drift velocity and variance of the solid volume fraction are predominant markers in the modelling of local drag force. Schneiderbauer (2017) employed the variance of the solid volume fraction and subgrid velocity fluctuations of the fluid and particle phases as markers in the drag closure. Furthermore, the anisotropy of the filtered drag was found to be crucial in the simulations (Cloete et al. 2018; Rauchenzauner & Schneiderbauer 2020).

Although the subgrid markers can account for the contribution of the mesoscale structure on the drag force, they are not available in coarse-grid simulations. The drift velocity can be modelled via a scale-similarity approach (Parmentier et al. 2012) and a gradient model (Ozel, Fede & Simonin 2013). By performing momentum balance analysis, Jiang, Chen & Zhou (2020) proposed an explicit model of drift velocity by introducing a scaled filtered pressure gradient. Rauchenzauner & Schneiderbauer (2020) proposed a dynamic anisotropic model, where the drift velocity was modelled by components of the filtered fluid-phase velocity fluctuations. Schneiderbauer & Saeedipour (2018, 2019) employed an approximate deconvolution model to reconstruct the unresolved flow field using the information around the local filtered cell, and then the approximate subgrid quantities were directly obtained by a deconvolution computation. Recently, the artificial neural network method has been shown to be a promising tool to construct the drift velocity implicitly (Zhang et al. 2020; Jiang et al. 2021; Zhu et al. 2021). Beetham et al. (2021) proposed models for the drift velocity and volume fraction variance using symbolic regression from Eulerian–Lagrangian simulation data.

In principle, the interface-resolved simulations are able to provide accurate drag forces on the particles, which is the advantage over the fine-grid two-fluid simulations and the Euler–Lagrange simulations for dense particle-laden flows. The drawback of the interface-resolved simulations is the huge computational cost. The interface-resolved simulations have been used to establish drag models for inhomogeneous systems. Wang, Liu & You (2011) proposed a drag model by introducing a non-uniformity coefficient. In some other structure-dependent models (Zhou et al. 2014; Li et al. 2017), the drag force depends on the volume fraction gradient and the slip velocity. Zhao et al. (2022) proposed a drag model that considers the effect of the solid volume fraction in the surrounding cells on the drag force at low Reynolds numbers.

From the perspective of filtered drag modelling using subgrid quantities, Rubinstein et al. (2017) modelled the effects of drift velocity and variance of the solid volume fraction on the local drag on the basis of interface-resolved simulations for the filter size ranging from $3d_p$ to $5d_p$. Following the methodology of Rubinstein et al. (2017), Zhao et al. (2021) proposed a drag model that depends on the solid drift flux and variance of the solid volume fraction that are available in DEM-CFD simulations. We note that these studies are limited to small-scale simulations of macroscopically homogeneous systems under low-Reynolds-number conditions.

Due to huge computational costs, there are very few works on mesoscale drag modelling based on large-scale DNS. Liu et al. (2020) proposed a scale-dependent drag model as a function of the grid Froude number from the results of large-scale interface-resolved simulations of heterogeneous gas–solid flows, but no subgrid markers were introduced. In the present study, we aim to check whether the existing mesoscale drag models can predict the drag in our large-scale interface-resolved simulations of inhomogeneous particle-laden flows, and also attempt to improve the mesoscale drag model with our DNS data. We first study the effects of the subgrid drift velocity and variance of the solid volume fraction on the filtered drag force and propose a model for mesoscale drag correction as a function of both filtered and subgrid quantities. Then the existing subgrid models and a new combined model for drift velocity and variance of the solid volume fraction are evaluated via a priori study. Finally, we assess the filtered drag models using the predicted subgrid quantities. To the best of our knowledge, this is the first work to employ the data of large-scale interface-resolved direct simulations of inhomogeneous particle-laden flows to develop a drag correlation for the two-fluid model via the filtering method or using the drift velocity.

The rest of the paper is organized as follows. In § 2 we describe the numerical method, simulation set-up and filtering method. In § 3 we present the results on the development of mesoscale drag correlations based on subgrid models and our DNS data. In § 4 we summarize the principal contributions of the present study.

2. Interface-resolved simulations of sedimentation suspensions

2.1. Direct-forcing/fictitious domain method

In the present study the simulations are conducted with the direct-forcing/fictitious domain (DF/FD) method (Yu & Shao 2007; Yu et al. 2016). The key idea of this method is to fill the interior domain of particles with fictitious fluids and rigid-body constraints on the fictitious fluid are enforced through a pseudo-body force (i.e. distributed Lagrange multiplier) in the particle domain (Glowinski et al. 1999).

The dimensionless formulation of the fictitious domain method for incompressible fluid (written for simplicity for a system with a single particle) can be written as

(2.1)$$\begin{gather} \frac{\partial {\boldsymbol{u}}}{\partial t}+{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}= \frac{\nabla^2 {\boldsymbol{u}}}{Re}-\boldsymbol{\nabla} p +{\boldsymbol{\lambda}} \quad \text{in} \ \varOmega, \end{gather}$$

(2.2)$$\begin{gather}{\boldsymbol{u}}={\boldsymbol{U}}+{\boldsymbol{\omega}_p} \times {\boldsymbol{r}}\quad \text{in} \ {P}, \end{gather}$$

(2.3)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}} =0 \quad \text{in} \ {\varOmega}, \end{gather}$$

(2.4)$$\begin{gather}(\rho_r -1)V_p^*\left(\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} t}-Fr \frac{\boldsymbol{g}}{ g}\right)={-}\int_P {\boldsymbol{\lambda}}\, {\rm d}{\kern0.8pt}{\boldsymbol{x}}, \end{gather}$$

(2.5)$$\begin{gather}(\rho_r -1)\frac{{\rm d} ({\boldsymbol{J}}^* \boldsymbol{\cdot}{\boldsymbol{\omega }}_p)}{{\rm d} t}={-}\int_P\boldsymbol{r}\times{\boldsymbol{\lambda}}\,{\rm d}{\kern0.8pt}{\boldsymbol{x}}, \end{gather}$$

where $\boldsymbol{u}$, $p$ and $\boldsymbol{\lambda}$, $\boldsymbol{r}$ represent the fluid velocity, pressure, pseudo-body force and position vector with respect to the mass centre of the particle, respectively. The particle domain and entire fluid/particle domain are denoted by $P$ and $\varOmega$. The viscosity and density of the fluid are $\mu$ and $\rho _f$. The particle density, volume and moment of inertia tensor, translational velocity and angular velocity are denoted by $\rho _s$, $V_p$, $\boldsymbol{J}$, $\boldsymbol{U}$ and $\boldsymbol{\omega}_p$, respectively. Here $\rho _r$ is the particle–fluid density ratio defined by $\rho _r=\rho _s/\rho _f$. The characteristic scales for the non-dimensionalization are as follows: $L_c$ for length, $U_c$ for velocity, $L_c/U_c$ for time, ${\rho _f}U_c^2/L_c$ for the body force. Here $Re$ is the Reynolds number defined by $Re=\rho _{f} U_c L_c/\mu$, $Fr$ is the Froude number defined by $Fr=gL_c/U_c^2$, with $g$ being the gravitational acceleration; $V_p^*$ and $\boldsymbol{J}^*$ are the dimensionless particle volume and moment of inertia tensor, defined by $V_p^*=V_p/L_c^3$ and $\boldsymbol{J}^*=\boldsymbol{J}/\rho _s L_c^5$.

A fractional-step temporal scheme is used to decouple the system (2.1)–(2.5) into two subproblems.

1. (1) Fluid subproblem for ${\boldsymbol {u}}^*$ and $p$ is given by

   (2.6)$$\begin{gather} \frac{\boldsymbol{u}^*-\boldsymbol{u}^n}{\Delta t}-\frac{1}{2}\frac{\nabla^2\boldsymbol{u}^*}{Re}={-}\boldsymbol{\nabla} p-\frac{1}{2}[3(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u})^n-(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u})^{n+1}]+\frac{1}{2}\frac{\nabla^2\boldsymbol{u}^n}{Re}+\boldsymbol{\lambda}^n, \end{gather}$$
   (2.7)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}^*=0. \end{gather}$$

   A finite-difference-based projection method on a homogeneous half-staggered grid is used. All spatial derivatives are discretized with a second-order central difference scheme.

2. (2) Particle subproblem for $\boldsymbol{U}^{n+1}$, $\boldsymbol{\omega }_p^{n+1}$, $\boldsymbol{\lambda }^{n+1}$ and $\boldsymbol{u}^{n+1}$ is given by

   (2.8)$$\begin{gather} \rho_r V_p^* \frac{\boldsymbol{U}^{n+1}}{\Delta t}=(\rho_r-1) V_p^* \left(\frac {\boldsymbol{U}^n}{\Delta t}-Fr\frac{\boldsymbol{g}}{g}\right)+\int_P \left({\frac{\boldsymbol{u}^*}{\Delta t}- \boldsymbol{\lambda}^n}\right){\rm d}{\kern0.8pt}{\boldsymbol{x}}, \end{gather}$$
   (2.9)$$\begin{gather}\rho_r \frac{\boldsymbol{J}^{*}\boldsymbol{\cdot} \boldsymbol{\omega }_p^{n+1}}{\Delta t}=(\rho_r-1) \left[\frac{\boldsymbol{J}^{*} \boldsymbol{\cdot} \boldsymbol{\omega }_p^{n}}{\Delta t}-\boldsymbol{\omega }_p^{n} \times (\boldsymbol{J}^{*} \boldsymbol{\cdot} \boldsymbol{\omega }_p^{n})\right]+\int_P \boldsymbol{r} \times \left({\frac{\boldsymbol{u}^*}{\Delta t}- \boldsymbol{\lambda}^n}\right) {\rm d}{\kern0.8pt}\boldsymbol{{x}}. \end{gather}$$
   The pseudo-body force $\boldsymbol{\lambda}$ is then updated from
   (2.10)\begin{equation} \boldsymbol{\lambda}^{n+1}=\frac{\boldsymbol{U}^{n+1}+\boldsymbol{\omega}_p^{n+1}\times \boldsymbol{r}-\boldsymbol{u}^*}{\Delta t}+\boldsymbol{\lambda}^n. \end{equation}
   The fluid–particle force acting on the particle can be directly obtained by integrating the pseudo-body force over the volume of each particle as
   (2.11)\begin{equation} \boldsymbol{F}^H={-}\int_P \boldsymbol{\lambda} \,\mathrm{d}{\kern0.8pt}\boldsymbol{x}+\frac{M}{\rho_{r}} \frac{\mathrm{d} \boldsymbol{U}}{\mathrm{d} t}. \end{equation}
   Finally, the fluid velocities $\boldsymbol{u}^{n+1}$ at the Eulerian nodes are corrected as
   (2.12)\begin{equation} \boldsymbol{u}^{n+1}=\boldsymbol{u}^{*}+\Delta t (\boldsymbol{\lambda}^{n+1}-\boldsymbol{\lambda}^{n}). \end{equation}

In the above manipulations, a discrete $\delta$ function in the form of a trilinear function is employed to transfer quantities between the Eulerian and Lagrangian nodes.

We implement a soft-sphere collision model with the lubrication correction to solve the short-ranged interparticle interaction. The soft-sphere collision model is a discrete element model based on a nonlinear soft-sphere approach. The lubrication force model is used to calculate the unresolved hydrodynamic forces exerted on the particles (Brändle de Motta et al. 2013). The time step for the collision process is set to be one tenth of that for the flow solution (i.e. $\Delta t_c = \Delta t/10$). Our collision model combined with the lubrication correction has been validated in our previous work (Xia et al. 2020), and the reader is referred to this work for the details of our collision model.

2.2. Simulation set-up

As shown in figure 1, particles freely settle in a fully periodic domain of dimensions $( L_{x},L_{y},L_{z}) = (320d_p,10d_p,80d_p)$. The gravitational acceleration $\boldsymbol {g}$ points along the negative $x$ direction. Initially, particles are placed randomly in the computational box, and the mean flux of the suspension is kept to zero by instantaneously adjusting the pressure gradient. The number of particles is given by

(2.13)\begin{equation} {N_p} = \left[ {\frac{{6{\langle \phi \rangle}}}{{\rm \pi} }\frac{{{L_x}{L_y}{L_z}}}{{d_p^3}}} \right], \end{equation}

where $[ ({\cdot } ) ]$ denotes the nearest integer and $\langle \phi \rangle$ denotes the domain-averaged solid volume fraction.

Figure 1. Snapshots for the large-scale IR-DNS of particle-laden flows for (a) case 1, $\langle \phi \rangle =0.1$, (b) case 2, $\langle \phi \rangle =0.15$, (c) case 3, $\langle \phi \rangle =0.2$, (d) case H1, $\langle \phi \rangle =0.1$.

We take the sphere diameter $d_p$ as the characteristic length, and the characteristic velocity $u_c$ is selected as

(2.14)\begin{equation} {u_c} = \sqrt {\frac{{4( {{\rho_r} - 1} )g{d_p}}}{{3}}}. \end{equation}

The Galileo number $Ga = \sqrt {d_p^3( {{\rho _r} - 1} )g} /\nu$ and the single-particle settling Reynolds number $Re_p = {u_t}{d_p}/\nu$ are widely used to measure the particle settling effect, with $u_t$ being the terminal settling velocity of an isolated particle in the form

(2.15)\begin{equation} {u_t} = \sqrt {\frac{{4( {{\rho_r} - 1} )g{d_p}}}{{3{C_{D}}}}}, \end{equation}

where $C_{D}$ is the standard drag coefficient

(2.16)\begin{equation} C_{D} = \frac{{24}}{{R{e_{{{p}}}}}}( {1 + 0.15Re_{{{p}}}^{0.687}} ). \end{equation}

Then, a relation between the single-particle settling Reynolds number and the Galileo number can be obtained as

(2.17)\begin{equation} G{a^2} = 18Re_p( {1 + 0.15Re_p^{0.687}}). \end{equation}

The particle Stokes number, which represents the ratio of the particle relaxation to the fluid relaxation time, is defined as

(2.18)\begin{equation} S{t_p} = \frac{{{\rho_s} {{{u}_t}} {d_p}}}{{18\mu }} = \frac{{{\rho_r}Re_p}}{{18}}. \end{equation}

The parameter settings for our simulations are listed in table 1. A large density ratio of $\rho _r=100$ is chosen to generate the clustered flows. To address the difference between the homogeneous and clustered flows, two additional simulations with $\rho _r=2$ are performed to produce nearly homogeneous flows, as shown in figure 1(d).

Table 1. Parameter settings for the cases studied.

The domain is discretized using a uniform Cartesian grid and the grid resolution with respect to the particle size is $d_p/\Delta x = 12.8$, which was shown to have good accuracy for particle Reynolds numbers below 100 for the sedimentation of a single sphere (Yu & Shao 2007; Xia et al. 2022b). We examine the effects of grid resolution on the drag force of moderately dense suspensions in Appendix A, and show that this grid resolution can produce reasonably accurate results. The dimensionless time step is $\Delta t=0.01$ with ${\rm CFL}=u_c \Delta t / \Delta x$ (where CFL denotes Courant–Friedrichs–Lewy) is about 0.15, and we set the dimensionless Young's modulus $E/(\rho _f u_c^2)=3 \times 10^5$, the Poisson ratio $\sigma$ = 0.33 and the friction coefficient $f = 0.3$ in the collision model. The sampling data are obtained for the duration of around 400 non-dimensional time units ($d_p/u_c$) with the sampling interval being 10 $d_p/u_c$, after the statistically steady state for the domain-averaged particle velocity is reached.

2.3. Filtering procedure

The goal of the present study is to analyse and model the local drag force for the filtered Euler–Euler approach with coarse fluid grids, based on interface-resolved simulation data. An efficient two-step filtering process (Capecelatro & Desjardins 2013) is employed. First of all, the Lagrangian particle variables are transferred onto Eulerian grids with a base cell size of $\varDelta _b = d_p$ through a mollification kernel $g_{\mathcal {M}}(r)$ as

(2.19)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle\phi = \sum_i^{{N_p}} {{g_{\mathcal{M}}}( {\boldsymbol{x} - {\boldsymbol{X}_i}} )V_p}, \\ \displaystyle{\boldsymbol{f}_p} = \sum_i^{{N_p}} {{g_{\mathcal{M}}}( {\boldsymbol{x} - {\boldsymbol{X}_i}} ){\boldsymbol{F}_i}}, \\ \displaystyle{\boldsymbol{u}_s} = \dfrac{1}{\phi}\sum_i^{{N_p}} {{g_{\mathcal{M}}}( {\boldsymbol{x} - {\boldsymbol{X}_i}} ){\boldsymbol{v}_i}}. \end{array}\right\} \end{equation}

Here, $\phi$, $\boldsymbol {u}_s$ and $\boldsymbol {f}_p$ are the microscopic solid volume fraction, particle velocity and fluid–particle force per unit volume on the base grid. The mollification kernel $g_{\mathcal {M}}(r)$ is defined as

(2.20)\begin{equation} g_{\mathcal{M}}(r)=\frac{1}{\sigma \sqrt{2 {\rm \pi}}} \,{\rm e}^{-{r^2}/{2 \sigma^2}}, \end{equation}

where $\sigma =\varDelta _b / 2 \sqrt {2 \ln 2}$ is the kernel width and $r$ is the distance between the particle and centre of the base cell. The resolved fluid variables are then coarsened to the base grid using the same mollification kernel (2.20) as

(2.21)\begin{equation} (1-\phi) {\boldsymbol{b}_{\mathcal{M}}}( {\boldsymbol{x}} ) = \iiint {\boldsymbol{b} {g_{\mathcal{M}}}} ( {{\boldsymbol{x}} - {\boldsymbol{x'}}} )\,{\rm{d}}\kern0.07em\boldsymbol{x'}. \end{equation}

Here $\boldsymbol {b}$ is any quantity associated with the fluid phase at the DNS grid, while ${\boldsymbol {b}_{\mathcal {M}}}( {\boldsymbol {x}} )$ is the corresponding fluid variable coarsened to the base Eulerian grid.

Then, the Eulerian field at the base grid obtained from the first step is further filtered to obtain the filtered quantity at any given filtering length scale. The filtered solid volume fraction is computed by

(2.22)\begin{equation} \overline{{\phi}({\boldsymbol{x}},t)} = \iiint {\phi}( {{\boldsymbol{x'}},t} )G ( {{\boldsymbol{x}} - {\boldsymbol{x'}}} )\,{\rm{d}}\kern0.07em\boldsymbol{x'}, \end{equation}

where $G ( {{\boldsymbol {x}} - {\boldsymbol {x'}}} )$ is a filter kernel that satisfies $\iiint G ( {{\boldsymbol {x}}} )\,\textrm {{d}}V = 1$ and a top-hat filter kernel is employed in this study. Due to the limitation of domain thickness in this work, we use a pseudo-two-dimensional (pseudo-2-D) filtering cell: the length of which along the limited direction ($y$ direction) is kept to be $L_y$ and the width of the filter cell in both remaining directions is $\varDelta _f$.

The filtered fluid and particle velocities are defined as

(2.23)$$\begin{gather} \widetilde{{{\boldsymbol{u}}_f}( {{\boldsymbol{x}},t} )} = \frac{1}{1-{\overline{{\phi}( {{\boldsymbol{x}},t} )} }}\iiint (1-{\phi} ( {{\boldsymbol{x'}},t} )){{\boldsymbol{u}}_f}( {{\boldsymbol{x'}},t} )G( {{\boldsymbol{x}} - {\boldsymbol{x'}}} )\,{\rm{d}}\kern0.07em{\boldsymbol{x'}}, \end{gather}$$

(2.24)$$\begin{gather}\widetilde{{{\boldsymbol{u}_s}}( {{\boldsymbol{x}},t} )} = \frac{1}{{\overline{{\phi}( {{\boldsymbol{x}},t} )} }}\iiint {\phi}( {{\boldsymbol{x'}},t} ){{\boldsymbol{u}}_s}( {{\boldsymbol{x'}},t} )G( {{\boldsymbol{x}} - {\boldsymbol{x'}}} )\,{\rm{d}}\kern0.07em{\boldsymbol{x'}}. \end{gather}$$

The filtered fluid–particle force is defined as

(2.25)\begin{equation} \overline{{\boldsymbol{f}}({\boldsymbol{x}},t)} = \iiint {\boldsymbol{f}_p}( {{\boldsymbol{x'}},t} )G ( {{\boldsymbol{x}} - {\boldsymbol{x'}}} )\,{\rm{d}}\kern0.07em\boldsymbol{x'}. \end{equation}

The pressure gradient force is excluded from the total hydrodynamic force, and the filtered fluid–particle force $\bar {\boldsymbol {f}}$ includes the drag force $\overline {\boldsymbol {f}_d}$, lift force $\overline {\boldsymbol {f}_l}$, added-mass force $\overline {\boldsymbol {f}_{am}}$ and history force $\overline {\boldsymbol {f}_{bs}}$. The lift force can be neglected for the sedimentation system. The added-mass and history forces depend on the relative acceleration between the particle and the undisturbed flow. From a scaling analysis, Ling, Parmar & Balachandar (2013) demonstrated that the added-mass and history forces are not important compared with the drag force for a large particle-to-fluid density ratio, and therefore, the added-mass and history forces are neglected in this work. For the homogeneous case when the density ratio is 2, the unsteady forces may be important for the instantaneous hydrodynamic force, however, only the average drag is studied for the homogeneous case and the contribution of the unsteady forces to the average drag should be unimportant, since the acceleration and deceleration of particles largely cancel each other after averaging.

We note that it is volume averaging (denoted by ‘$\bar {X}$’) for the solid volume fraction and the drag force, whereas it is phase averaging(weighted with the volume fraction and denoted by ‘$\tilde {X}$’) for the fluid and solid velocities. The phase averaging was also called Favre averaging in some previous works (Rubinstein et al. 2017; Liu et al. 2020; Zhao et al. 2021).

We define the force parallel to the filtered slip velocity as drag force and a dimensionless filtered drag force (actually drag correction coefficient) as

(2.26)\begin{equation} F = \frac{{( {{V_p}/\bar \phi } )\bar{\boldsymbol{f}} \boldsymbol{\cdot} \widetilde{{{\boldsymbol{u}}_{slip}}}}}{{3{\rm \pi} \mu {d_p} \widetilde{{{\boldsymbol{u}}_{slip}}} \boldsymbol{\cdot} \widetilde{{{\boldsymbol{u}}_{slip}}}}}, \end{equation}

where $\widetilde {\boldsymbol {u}_{slip}}=\widetilde {\boldsymbol {u}_{f}}-\widetilde {\boldsymbol {u}_{s}}$ is the filtered slip velocity and a filtered Reynolds number based on the magnitude of the filtered slip velocity between the two phases is defined as

(2.27)\begin{equation} \widetilde{Re_m} = \frac{{(1 - \bar \phi )| { \widetilde{{{\boldsymbol{u}}_{slip}}}} |{d_p}}}{\nu }. \end{equation}

We note that the filter procedure is almost equivalent to the one employed in the post-process of both interface-resolved simulations (Rubinstein et al. 2017; Liu et al. 2020; Zhao et al. 2021) and DEM-CFD simulations (Ozel et al. 2017). Figure 2 shows the instantaneous filtered fields of case 3 with $\varDelta _f/d_p=8$ in a 2-D plane.

Figure 2. Snapshots of the filtered variable field along a slice obtained for case 3 and $\varDelta _f/d_p = 8$ for (a) the particle volume fraction $\bar \phi$; (b) the non-dimensionalized vertical fluid velocity $\widetilde {{{u}_{f,1}}}/{u_t}$; (c) the non-dimensionalized vertical particle velocity $\widetilde {{{u}_{s,1}}}/{u_t}$; (d) the non-dimensionalized vertical slip velocity $\widetilde {{{u}_{slip,1}}}/{u_t}$.

To put domain-averaged analysis into context, we employ an averaging operator $\langle {( {\cdot } )} \rangle$ throughout this work to denote a time averaging and spatial filtering over the entire time and space, which is not a function of either the spatial coordinate or time.

3. Results and discussion

The drag models obtained from small-scale systems with an homogeneous and isotropic microstructure are usually referred to as microscopic drag models. Some microscopic models are summarized in table 2, which can typically be expressed as functions of the microscopic solid volume fraction $\phi$ and microscopic Reynolds number ${Re_m}$ that is defined as $Re_m = (1 - \phi )u_{slip}{d_p}/\nu$, where $u_{slip}$ represents the phase-averaged velocity difference (i.e. slip velocity).

Table 2. Some previous drag correlations for the homogeneous system.

In figure 3 the conditionally averaged dimensionless drag force is plotted as a function of $\varDelta _f/d_p$ at $(\bar \phi, \widetilde {Re_m})= (0.1,3)$ and $(0.1,15)$ for both inhomogeneous (cases 1 to 5) and homogeneous (cases H1 and H2) cases in table 1. The drag force $F$ is computed using a range of filter sizes ($\varDelta _f/d_p$ from 4 to 40) and bin averaged with both $\bar \phi$ and $\widetilde {Re_m}$. From figure 3, the averaged dimensionless drags for inhomogeneous flows ($\rho _r = 100$) and homogeneous flows ($\rho _r = 2$) are close to each other at a small filter size of $\varDelta _f = 4d_p$. There is a roughly 30 % reduction in averaged dimensionless drag at both $\widetilde {Re_m}=3$ and 15 when the filter size increases to $20d_p$ for the inhomogeneous flows, while the averaged $F$ is independent of the filter size for the homogeneous flows. The standard deviations of $F$ shown by error bars indicate that the fluctuations in $F$ at the same $\bar \phi$ and $\widetilde {Re_m}$ are much larger in the inhomogeneous systems. Thus, the homogeneous drag correlations are insufficient for capturing the actual drag in the inhomogeneous particle-laden flows and further drag correction is necessary. By comparing the microscopic drag models in table 2, the model of Xia et al. (2022a) agrees best with our results of the homogeneous cases, therefore, we take the microscopic model of Xia et al. (2022a) as the microscopic drag model in our work for further modelling of the filtered drag force.

Figure 3. Averaged dimensionless drag force $F$ as a function of normalized filter sizes $\varDelta _f/d_p$, at $\widetilde {Re_m} = 3$, 15 and $\bar \phi = 0.1$. The error bars are given by measuring the standard deviation from the bin-averaged dimensionless drag force.

The variation of average drag force with filter size may be related to the size of the particle cluster, which can be quantified by the pair distribution function $P(\boldsymbol {r})$. Due to the pseudo-2D particle cluster, we calculate the pair distribution function in a 2-D plane. In the cylindrical coordinate frame, $P(\boldsymbol {r})$ is defined as

(3.1)\begin{equation} P(\boldsymbol{r})=P(r, \theta)=\frac{H(r, \theta)}{n t_s \Delta S}, \end{equation}

where $r$ is the distance from the centre of a reference particle to the centre of a neighbouring particle, $\theta$ is the polar angle (measured from the positive $x$ axis), $H(r,\theta )$ is the histogram of particle pairs, $\Delta S=r\Delta r \Delta \theta$ is the area of the sampling bin, $n$ is the average particle number density in the plane and $t_s$ is the total number of sampling points. By averaging $P(r,\theta )$ over $\theta$, one obtains the radial distribution function $g(r)$.

The specific goal here is to identify the cluster structure. Figure 4(a) shows the radial distribution functions $g(r)$ as a function of normalized interparticle distance $r/d_p$. It is observed that $g(r)$ is greater than unity for $r/d_p < 20$ for clustered flows, whereas it becomes much more uniform for the cases of $\rho _r = 2$. In figure 4(b,c) the pair distribution function in the $x$ and $z$ directions are plotted respectively, and figure 4(d) shows the particle pair distributions for both the inhomogeneous and homogeneous systems. The relatively strong anisotropy in $P(r,\theta )$ is observed at a relatively small $Ga$ of 10 and $\rho _r=100$, and the anisotropy weakens as $Ga$ increases. The results in figure 4 show that the size of the clusters in the horizontal direction characterized by $P(r,0)>1$ is around $20d_p$ for all cases of clustered flows, corresponding to the filter size beyond which the average drag force does not change much, as shown in figure 3. From our results, the saturation in the average drag is correlated with the size of the clusters in the horizontal direction; however, further studies are needed to confirm this correlation.

Figure 4. Radial distribution function $g(r)$ and pair distribution function $P(r,\theta )$ for cases of $\langle \phi \rangle = 0.1$.

The filtered drag force in coarse cells can be corrected from the microscopic drag force $F_{mic}$ for the homogeneous system at the same filtered solid volume fraction and filtered Reynolds number by introducing a correction function or the mesoscale index $H$ (Igci & Sundaresan 2011; Milioli et al. 2013; Sarkar et al. 2016; Schneiderbauer 2017), i.e.

(3.2)\begin{equation} F = ( {1 + H} ){F_{mic}}, \end{equation}

namely,

(3.3)\begin{equation} H=\frac{F-F_{mic}(\bar{\phi}, \widetilde{R e_m})}{F_{mic}(\bar{\phi}, \widetilde{R e_m})}. \end{equation}

The mesoscale correction $H$ is the ultimate closure that is required for the filtered drag force in the present study.

3.1. Effects of subgrid quantities on filtered drag

The drift velocity and variance of the solid volume fraction were found to be key subgrid quantities due to their ability to capture subgrid-scale structures (Parmentier et al. 2012; Ozel et al. 2013). The drift velocity is defined as

(3.4)\begin{equation} {{\boldsymbol{u}}_{drift}} = {( {\widetilde{{{\boldsymbol{u}}_f}}})_p} - { {\widetilde{{{\boldsymbol{u}}_f}}}}, \end{equation}

where $(\widetilde {\boldsymbol {u}_f})_p = \overline {\phi {\boldsymbol {u}_f}} /\bar {\phi }$ denotes the filtered fluid velocity seen by the particles. The variance of the solid volume fraction is defined as

(3.5)\begin{equation} \overline {{{( {\phi '} )}^2}} = \overline {{\phi^2}} - {\bar \phi^2}. \end{equation}

The effect of subgrid quantities on the filtered drag can be derived via a series expansion to the microscopic drag force (Schneiderbauer 2017; Chen et al. 2020). The drift velocity and variance of the solids volume fraction were identified as crucial factors for the drag force modification in filtered two-fluid simulations.

A scaled drift flux is defined as the normalized component of the drift flux in the direction of the slip velocity (Rubinstein et al. 2017):

(3.6)\begin{equation} {v_d} = \frac{{( {\bar \phi {u_{drift,\beta }}} ) \widetilde{{u_{slip,\beta }}}}}{{\widetilde{{u_{slip,\beta }}} \widetilde{{u_{slip,\beta }}}}}. \end{equation}

Following the work of Ozel et al. (2017), a scaled variance of the solid volume fraction $\overline {{{( {\phi '} )}^2}} /\bar \phi ( {1 - \bar \phi } )$ is employed as a second marker in the drag closure.

The conditionally averaged mesoscale correction $H$ obtained from the IR-DNS data as a function of the scaled drift flux and variance of the solid volume fraction for $\varDelta _f/d_p = 24$ is shown in figure 5. The results in figure 5(a) show that $H$ has a strong positive correlation with the normalized drift velocity and is affected by the scaled variance of the solid volume fraction. From figure 5(b), $H$ increases with increasing scaled variance of the solid volume fraction, being consistent with the previous observations of Ozel et al. (2017). For a large drift velocity, $H$ is largely independent of the fluctuation of the solid volume fraction. For a typical coarse grid, the ratio between the filtered drag force and microscale drag force (i.e. $1+H$) is in the range of 0 to 2. Such a result suggests that the combination of scaled drift flux and variance of the solid volume fraction provides a significant mesoscale drag correction.

Figure 5. Variation of the averaged mesoscale correction $H$ with (a) the scaled drift flux $v_d/\bar \phi$, (b) the scaled variance of the solid volume fraction $\overline {{{( {\phi '} )}^2}} /\bar \phi ( {1 - \bar \phi } )$ at $\varDelta _f/d_p = 24$ for cases 1–5. The symbols represent the results obtained from our simulations and the lines denote the fitting function in (3.7).

3.2. Modelling the mesoscale correction coefficient

According to the observation in figure 5, we employ a function given as

(3.7)\begin{align} H &= f\left( {\frac{{{v_d}}}{{\bar \phi }},\frac{{\overline {{{( {\phi '} )}^2}} }}{{\bar \phi ( {1 - \bar \phi } )}}} \right)\nonumber\\ & = \frac{v_d}{\bar{\phi}}+A \frac{\overline{(\phi^{\prime})^2}}{\bar{\phi}(1-\bar{\phi})}\left(\frac{v_d}{\bar{\phi}}+1\right)^2\nonumber\\ &\quad +B\left(\frac{\overline{(\phi^{\prime})^2}}{ \bar{\phi}(1-\bar{\phi})}\right)^2\left(\frac{v_d}{\bar{\phi}}+1\right)+C \frac{\overline{(\phi^{\prime})^2}}{\bar{\phi}(1- \bar{\phi})}\left(\frac{v_d}{\bar{\phi}}+1\right). \end{align}

If we set $A=B=C=0$, the simplest approximation of $H$ in terms of the scaled drift flux (Parmentier et al. 2012) is given by

(3.8)\begin{equation} H_{{Par}}=\frac{v_d}{\bar{\phi}}. \end{equation}

The model of Ozel et al. (2017) can be obtained if $A=B=0$ and $C=2.25\pm 0.5$ as

(3.9)\begin{equation} H_{{Ozel}}=\frac{v_d}{\bar{\phi}}+C\frac{\overline{(\phi^{\prime})^2}}{\bar{\phi}(1-\bar{\phi})} \left(\frac{v_d}{\bar{\phi}}+1\right). \end{equation}

In the present study, the constants are obtained by fitting the data from cases 1 to 5 in table 1 at various filter sizes. Figure 6 shows that all constants are independent of filter size for $\varDelta _f/d_p > 12$. We note that $H$ is insensitive to the coefficients at small filter sizes, with the dependence of drag force on the filter size being hidden in the subgrid quantities. We just employ $A = 7.3$, $B= 25.4$, $C = -3.3$ for all filter sizes and the corresponding mesoscale correction coefficient is denoted as $H_{{Present}}$.

Figure 6. Model coefficients $A$, $B$ and $C$ in (3.7) as a function of filter size.

We first assess the accuracy of the proposed filtered drag model by parity plots between the exact drag force from the IR-DNS data and the predicted drag force using (3.2) and (3.7). The results in figure 7 show that the drag model with actual drift velocity and variance of the solid volume fraction computed from the IR-DNS data has a good performance in drag prediction.

Figure 7. The parity plots between the exact and predicted drag force $F$ for (a) $\varDelta _f/d_p = 12$, (b) $\varDelta _f/d_p = 20$, (c) $\varDelta _f/d_p = 28$, (d) $\varDelta _f/d_p = 36$.

To further quantify the statistical error of the models, we define the relative error between the exact values and the predicted values as

(3.10)\begin{equation} {e}(y;\hat y) = \frac{{{{\hat y}_i} - {y_i}}}{{{y_i}}}, \end{equation}

where $y$ and $\hat y$ are the exact and predicted data, respectively. The Pearson correlation coefficient $r$ between the exact value obtained from IR-DNS and that predicted by the proposed model is calculated to quantify the performance of the proposed model. The Pearson correlation coefficient is defined as

(3.11)\begin{equation} r(x ; y)=\frac{\displaystyle\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\displaystyle\sqrt{\sum_{i=1}^n(x_i-\bar{x})^2} \sqrt{\sum_{i=1}^n(y_i-\bar{y})^2}}, \end{equation}

with $x$ and $y$ being two sets of data, and $\bar {x}$ and $\bar {y}$ being their mean values. Higher $r(x;y)$ indicates a better linear correlation between two datasets. In Appendix B we evaluate the drag models by computing the probability distribution functions (p.d.f.s) of the relative error and Pearson correlation coefficients between the drag force obtained from mesoscale models and the homogeneous model of Xia et al. (2022b) and the exact forces obtained from the IR-DNS. The results indicate that our new correction yields the highest correlation coefficients.

The subgrid-quantity-dependent filtered drag correlation (3.7) may be affected by the computational domain size. In Appendix C we have studied the effect of domain size on the drag prediction and demonstrated that the present filtered drag model with actual subgrid quantities is applicable to different domain sizes.

3.3. Priori analysis for subgrid-scale quantities

3.3.1. Existing subgrid closures

In the two-fluid model the actual subgrid drift velocity $\boldsymbol {u}_{drift}$ and variance of the solid volume fraction ${\overline {{{( {\phi '} )}^2}}}$ are not available and need to be modelled. One classical closure approach for the subgrid quantities is the dynamic scale-similarity model (Parmentier et al. 2012; Rubinstein et al. 2017), which approximates the subgrid quantities via available information at or above the filter grid length scale. Figure 8 shows the conditionally averaged drift flux $v_d$ and variance of the solid volume fraction $\overline {{{( {\phi '} )}^2}}$ as a function of ${\bar \phi /{\phi _{max}}}$ with ${\phi _{max}}=0.64$. The results show that the profiles of $v_d$ and $\overline {{{( {\phi '} )}^2}}$ are similar, and they can be modelled as a function of the filter solid volume fraction and filter size.

Figure 8. Averaged (a) drift flux $v_d$ and (b) variance of the solid volume fraction $\overline {{{( {\phi '} )}^2}}$ as a function of the normalized filtered solid volume fraction ${\bar \phi /{\phi _{max}}}$ with ${\phi _{max}}=0.64$ for various filter sizes. The symbols correspond to results from the interface-resolved simulations and the lines correspond to the predictions of (3.12)–(3.15) with $\kappa _1=\kappa _2=1$.

Following the work of Rubinstein et al. (2017), the scale-similar forms for the drift flux and variance of the solid volume fraction are assumed to be

(3.12)$$\begin{gather} \bar \phi {u^{{Rub}}_{drift,\beta }} = {\kappa_{1}}{h_1}( {\bar \phi /{\phi_{max}} } ){f_1}( {{\varDelta_f/d_p}} )\widetilde{{u_{slip,\beta }}}, \end{gather}$$

(3.13)$$\begin{gather}{\overline {{{( {\phi '} )}^2}}}^{{Rub}} = {\kappa_2}{h_2}( {\bar \phi /{\phi_{max}}} ){f_2}( {{\varDelta_f/d_p}} ). \end{gather}$$

The functions ${h}_i( {\bar \phi /{\phi _{max}} } )$ and $f_i( {{\varDelta _f/d_p}} )$ with $i=1,2$ have the following forms:

(3.14)$$\begin{gather} h_i( {\bar \phi /{\phi_{max}}} ) = {( {\bar \phi /{\phi_{max}}} )^{n_i}}{( {1 - \bar \phi /{\phi_{max}}})^{m_i}}, \end{gather}$$

(3.15)$$\begin{gather}f_i( {{\varDelta_f/d_p}} ) = \frac{{{(\varDelta_f/d_p})^2}}{{k_i + {(\varDelta_f/d_p)^2}}}. \end{gather}$$

Based on our data from the simulations in table 1, we obtain $m_1 = 4.84$, $n_1 = 2.59$, $m_2 = 4.88$, $n_2 = 2.32$, $k_1 = 235.5$ and $k_2 = 269.7$. We note that the model constants $\kappa _1$ and $\kappa _2$ are computed dynamically through the surrounding resolved flow field (see Rubinstein et al. 2017 for details). From figure 8, the above model can predict very well the averaged drift velocity $v_d$ and variance of the solid volume fraction $\overline {{{( {\phi '} )}^2}}$ for a given particle volume fraction and a given filter size, simply with $\kappa _1=\kappa _2=1$.

Another widely used closure for the subgrid quantities is the gradient approach based on a Taylor expansion (Ozel et al. 2013) as

(3.16)$$\begin{gather} {\bar \phi u^{{Grad}}_{drift,\beta }} = {D}\varDelta_f^2\frac{{\partial \bar \phi }}{{\partial {x_j}}}\frac{{\partial \widetilde{{u_{f,\beta }}}}}{{\partial {x_j}}}, \end{gather}$$

(3.17)$$\begin{gather}{\overline {{{( {\phi '} )}^2}}}^{{Grad}} = E{\varDelta_f^2}\frac{{\partial \bar \phi }}{{\partial {x_j}}}\frac{{\partial \bar \phi }}{{\partial {x_j}}}, \end{gather}$$

where $D$ and $E$ are the model coefficients that can be determined from a linear regression of the IR-DNS data. It should be noted that the local gradient in the filtered particle volume fraction and fluid velocity used in the gradient-based models can be computed at different scales using the IR-DNS data. In this study, all gradient values are computed at both fine-grid (i.e. particle-size) and coarse-grid (i.e. filter-size) scales using the second-order central difference scheme, as shown in figure 9. We use the superscripts ‘fine’ and ‘coarse’ to denote fine-grid (particle-size) and coarse-grid (filter-size) scale discretization on the computation of gradient quantities, respectively.

Figure 9. Schemes for the computation of the gradient of filtered quantities by the central difference scheme. (a) Fine-grid (i.e. particle-size) scale discretization, (b) coarse-grid (i.e. filter-size) scale discretization.

Figure 10 shows the average drift flux $v_d$ and variance of the solid volume fraction ${\overline {{{( {\phi '} )}^2}}}$ as a function of $G_{v_d}$ and $G_{{\overline {{{( {\phi '} )}^2}}}}$, which are defined by

(3.18a,b)\begin{equation} G_{v_d} = { \varDelta_f^2}\frac{{\partial \bar \phi }}{{\partial {x_j}}}\frac{{\partial \widetilde{u_{f,\beta }}}}{{\partial {x_j}}}\frac{{ \widetilde{u_{slip,\beta }}}}{{\widetilde{u_{slip,\beta }} \widetilde{u_{slip,\beta }}}}\quad {\rm{and}}\quad G_{{\overline {{{( {\phi '} )}^2}}}} = \varDelta_f^2 \frac{{\partial \bar \phi }}{{\partial {x_j}}}\frac{{\partial \bar \phi }}{{\partial {x_j}}}. \end{equation}

The gradient model ((3.16)–(3.17)) can be written as $v^{{Grad}}_d = D G_{v_d}$ and ${\overline {{{( {\phi '} )}^2}}}^{{Grad}} = E G_{{\overline {{{( {\phi '} )}^2}}}}$. The results in figure 10 indicate a robust positive correlation between the averaged drift flux $v_d$ (or variance of the solid volume fraction ${\overline {{{( {\phi '} )}^2}}}$) with $G_{v_d}$ (or $G_{{\overline {{{( {\phi '} )}^2}}}}$). The subgrid-scale quantities can be effectively modelled using the closures (3.16)–(3.17). The model coefficients mainly depend on the filter size. Explicit functions are obtained by data fitting as

(3.19)\begin{equation} \left.\begin{array}{@{}c@{}} D_{{fine}} ={-} 0.094 + 0.23\exp ( { 0.023{\varDelta_f}/{d_p}} ), \\ D_{{coarse}} ={-} 0.094 + 0.04\exp ( { 0.14{\varDelta_f}/{d_p}} ), \\ E_{{fine}} = 0.27 + 0.07\exp ( { 0.048{\varDelta_f}/{d_p}} ),\\ E_{{coarse}} ={-}0.75 + 0.93\exp ( { 0.055{\varDelta_f}/{d_p}} ). \end{array}\right\} \end{equation}

The model coefficients increase with increasing filter size and are larger than the theoretical value of 1/12 derived by Taylor expansion (Parmentier et al. 2012), especially for the coarse-grid scale.

Figure 10. Averaged (a,b) drift flux $v_d$ and (c,d) variance of the solid volume fraction $\overline {{{( {\phi '} )}^2}}$ as a function of $G_{v_d}$ or $G_{{\overline {{{( {\phi '} )}^2}}}}$ for various filter sizes, using (a,c) the fine grid (particle size) and (b,d) the coarse grid (filter size) for the gradient computation.

The drift velocity can also be expressed as a function of the covariance of the solid volume fraction and the fluid-phase velocity (Cloete et al. 2019), which can be further modelled as a function of the filtered fluid-phase velocity fluctuation and the variance of the solid volume fraction (Schneiderbauer 2017; Rauchenzauner & Schneiderbauer 2020) as

(3.20)\begin{equation} {\bar \phi u^{{Sch}}_{drift,\beta }} = {\zeta_\phi }\frac{{\sqrt {2{k_{f,\beta }}\overline {{({\phi '}^2)}} } }}{{ {1 - \bar \phi }}}. \end{equation}

In the above equation, $2{k_{f,\beta }} = \widetilde {{u'_\beta }{u'_\beta }}$ is the Reynolds normal stress (here no summation for ${\beta }$). The correlation coefficient ${\zeta _\phi }$ can be obtained via data fitting or dynamically computing (Rauchenzauner & Schneiderbauer 2020). In the work of Schneiderbauer (2017), an algebraic expression for $\overline {({{\phi '}^2})}$ is derived as

(3.21)\begin{equation} {\overline {{{( {\phi '} )}^2}}}^{{Sch}} = \frac{{8\xi_\phi^2{k_f}}}{{{{\left( {\dfrac{{\partial \widetilde{{u_{f,k}}}}}{{\partial {x_k}}} + {C_\phi }{C_\varepsilon }\dfrac{{k_f^{1/2}}}{{{l_{mf}}}}} \right)}^2}}}\frac{{\partial \bar \phi }}{{\partial {x_i}}}\frac{{\partial \bar \phi }}{{\partial {x_i}}}, \end{equation}

where $l_{mf}=C_\nu \varDelta _f$ denotes the mixing length of the fluid phase; $C_\phi = 0.4$, $C_\nu = 0.4$ and $C_\varepsilon = 1$ are model constants suggested by Schneiderbauer (2017). We note that (3.21) can be regarded as a gradient-based model for the variance of the solid volume fraction. Explicit functions for the coefficients $\xi _\phi$ and $\zeta _\phi$ appearing in (3.20) and (3.21) are obtained by data fitting:

(3.22)\begin{equation} \left.\begin{array}{@{}c@{}} \zeta_\phi ={-} 0.43 + 0.26\exp ( { - 0.09{\varDelta_f}/{d_p}}), \\ \xi_{\phi,{fine}} = 0.48+0.1 \exp ( { 0.024{\varDelta_f}/{d_p}}), \\ \xi_{\phi,{coarse}} ={-}1.14+1.83 \exp ( { 0.02{\varDelta_f}/{d_p}}). \end{array}\right\} \end{equation}

As shown in figure 11, the function (3.20) is a good approximation of the drift flux. The coefficient $\zeta _\phi$ is found to be negative, with its magnitude reaching up to 0.43 as the filter size increases.

Figure 11. (a) Conditionally averaged drift flux $v_d$. The lines correspond to the predictions of (3.20). (b) Model constant ${\zeta _\phi }$ in (3.20) as a function of the dimensionless filter size.

The performance of existing models for subgrid quantities is evaluated in Appendix B, which shows that the model proposed by Schneiderbauer (2017) for the drift velocity is desirable at relatively large filter sizes, while the gradient model yields the lowest correlation coefficient at large filter sizes with coarse-grid discretization.

One drawback of the gradient model is that when the gradient of the particle volume fraction is close to zero, it predicts a very small $| {v_d} |$ and ${\overline {{{( {\phi '} )}^2}}}$, but the actual $| {v_d} |$ and ${\overline {{{( {\phi '} )}^2}}}$ may not be so small, particularly at large filter sizes (see figure 10). In figure 12 we depict the conditionally averaged variance of the solid volume fraction $\overline {{{( {\phi '} )}^2}}$ as a function of ${\bar \phi /{\phi _{max}}}$ in the limit of $( {\partial \bar \phi /\partial {x_j}} )( {\partial \bar \phi /\partial {x_j}} ) \to 0$. Interestingly, the profiles of $\overline {{{( {\phi '} )}^2}}$ are similar to those for all gradients of the solid volume fraction shown in figure 8(b). Based on present simulation data, we propose the following model for the variance of the solid volume fraction $\overline {{{( {\phi '} )}^2}}$:

(3.23)\begin{equation} {\overline {{{( {\phi '})}^2}}}^{{Present}} = E'{\varDelta_f^2}\frac{{\partial \bar \phi }}{{\partial {x_j}}}\frac{{\partial \bar \phi }}{{\partial {x_j}}}+{( {\bar \phi /{\phi_{max}}} )^{{n_2}}}{( {1 - \bar \phi /{\phi_{max}}} )^{{m_2}}}\frac{{{{({\varDelta_f}/{d_p})}^2}}}{{{k_2} + {{({\varDelta_f}/{d_p})}^2}}}. \end{equation}

Here the values of $n_2$, $m_2$, $k_2$ are the same as those in the scale-similarity model of Rubinstein et al. (2017) presented earlier (3.13), and $E' = 0.05+0.17 \exp ( { -0.05{\varDelta _f}/{d_p}} )$. The value of $E'$ is weakly dependent of filter size and close to the theoretical value 1/12 as derived from Taylor expansion (Parmentier et al. 2012) for $\varDelta _f/d_p$ greater than 24. Note that the gradient of the solid volume fraction in our model above is calculated with the filter size (i.e. coarse grid). Our model for $\overline {{{( {\phi '} )}^2}}$ is essentially a combination of the gradient model and the scale-similarity model.

Figure 12. Averaged variance of the solid volume fraction $\overline {{{( {\phi '} )}^2}}$ as a function of ${\bar \phi /{\phi _{max}}}$ with ${\phi _{max}}=0.64$ for various filter sizes when $( {\partial \bar \phi /\partial {x_j}} )( {\partial \bar \phi /\partial {x_j}} ) \to 0$. The symbols correspond to results from the interface-resolved simulations and the lines correspond to the predictions of $( {\bar \phi /{\phi _{max}}} )^{{n_2}}( {1 - \bar \phi /{\phi _{max}}} )^{{m_2}}({\varDelta _f}/{d_p})^2/[ {k_2} + {{({\varDelta _f}/{d_p})}^2} ]$.

Due to the great performance of the model of Schneiderbauer (2017) for the prediction of the drift velocity, we choose this closure to compute the drift velocity:

(3.24)\begin{equation} {\bar \phi u^{{Present}}_{drift,\beta }} = {\zeta_\phi }\sqrt {2{k_{f,\beta }}{{\overline {({{\phi '}^2})} }^{{{Present}}}}} /( {1 - \bar \phi } ). \end{equation}

Here $\overline {{{( {\phi '} )}^2}}$ is computed with (3.23).

Figure 13(a,b) presents the predictions of the drift velocity and the variance of the solid volume fraction using the present model (3.23)–(3.24) with a coarse-grid scale discretization at $\varDelta _f/d_p = 20$. Figure 13(c,d) shows the relative error at different filter sizes. It can be seen that the predictions are in good agreement with the actual values and the performance of the model on the relative error is even better for a larger filter size.

Figure 13. Parity plots between the exact and predicted (a) drift flux $v_d$ and (b) variance of the solid volume fraction $\overline {{{( {\phi '} )}^2}}$ for $\varDelta _f/d_p = 20$. Probability density functions of the relative error between the actual and predicted (c) drift flux $v_d$ and (d) variance of the solid volume fraction $\overline {{{( {\phi '} )}^2}}$ using model (3.23)–(3.24) at different filter sizes.

The predicted $v_d$ and $\overline {{{( {\phi '} )}^2}}$ and the actual quantities conditionally averaged for $\bar \phi =0.1$ and 0.2 are compared in figure 14. The superscripts ‘fine’ and ‘coarse’ refer to fine-grid (particle-size) and coarse-grid (filter-size) scale discretization on the computation of gradient quantities, respectively, while ‘actual’ in the model of Schneiderbauer (2017) indicates that the actual ${\overline {{{( {\phi '} )}^2}}}$ obtained from the interface-resolved simulations is employed. The results show that all mesoscale drag models can predict well the average $v_d$ and $\overline {{{( {\phi '} )}^2}}$ for different filter sizes and different solid volume fractions. The present subgrid closure agrees best with our DNS results. The error analysis in Appendix B shows that the present model exhibits a better performance on the correlation coefficient and relative error than the other models using the resolved quantities at a coarse-grid scale. In addition, our model can be utilized without dynamic adjustment of the model constant, unlike the dynamic scale-similarity model.

Figure 14. Comparison of $v_d$ and $\overline {{{( {\phi '} )}^2}}$ averaged for the solid volume fraction (a,c) $\bar \phi = 0.1$ and (b,d) $\bar \phi = 0.2$. The lines correspond to the actual $v_d$ and $\overline {{{( {\phi '} )}^2}}$ from the interface-resolved simulations, and the symbols correspond to the predictions.

It is difficult to predict the subgrid quantities using the resolved filtered quantities at a coarse-grid scale accurately because the resolved filtered quantities are insufficient for intrinsically characterizing the feature of subgrid-scale heterogeneous structures. It is possible to improve the prediction of the subgrid quantities by using additional transport equations or developing neural-network-based models for the subgrid quantities (Zhang et al. 2020; Zhu et al. 2021), which is outside the scope of the present study.

3.3.2. Evaluation of filtered drag models with models for subgrid quantities

We have evaluated the performance of filtered drag models with the actual subgrid quantities as well as the closure models for the subgrid quantities. In this part, further assessment of the filtered drag models with the closures of subgrid quantities is conducted. A coarse grid is normally used in the simulations of industrial-scale particle-laden flows with the two-fluid model, therefore, we only consider the models using the filter size to compute the gradient, i.e. the ‘coarse’ gradient model.

We plot the conditionally averaged dimensionless drag force (both the actual drag and predicted drag) as a function of $\varDelta _f/d_p$ at $(\bar \phi, \widetilde {Re_m}) = (0.1,15)$ and $(\bar \phi, \widetilde {Re_m}) = (0.2,2.0)$ in figure 15(a,b). The results show that all mesoscale drag models using the resolved filtered quantities at a coarse-grid scale can predict well the averaged drag and are superior to the homogeneous drag model. Figure 15(c,d) presents the comparisons between the predicted and actual drag forces at $\varDelta _f/d_p=16$ and $\varDelta _f/d_p=28$. Our model agrees best with our DNS results, although the improvement over the previous mesoscale drag models is not so significant.

Figure 15. Averaged dimensionless drag force $F$ as a function of normalized filter sizes $\varDelta _f/d_p$ at (a) $(\bar \phi, \widetilde {Re_m}) = (0.1,15)$ and (b) $(\bar \phi, \widetilde {Re_m}) = (0.2,2.0)$. Predictions of the drag force $F$ from subgrid-quantity-dependent drag models at (c) $\varDelta _f/d_p = 16$ and (d) $\varDelta _f/d_p = 28$. The error bars are given by measuring the standard deviation from the bin-averaged dimensionless drag force. In (a) and (b) the green dashed lines represent the prediction of the homogeneous drag model of Xia et al. (2022a) and the black solid lines correspond to the exact $F$ from the interface-resolved simulations.

In Appendix D we further assess the filtered drag models for small fluidized beds and the results show that all mesoscale models are applicable to small-scale systems.

4. Concluding remarks

In the present work, we have investigated the mesoscale drag models for clustered particle-laden flows based on the data from large-scale interfaced-resolved simulations. The main conclusions are drawn as follows.

1. (i) For the homogeneous system, the filtered drag is independent of the filter size, whereas for the clustered particle-laden flows, the averaged drag becomes smaller than the homogeneous drag at the filter size above 4 particle diameters. The drag reduction saturates at the filter size being comparable to the cluster size in the horizontal direction in our simulations. The dimensionless averaged drag force decreases as the drift velocity decreases at the same variance of the solid volume fraction, or as the variance of the solid volume fraction decreases at the same drift velocity.

2. (ii) A new correlation (3.7) for the mesoscale drag correction is proposed as a function of the drift flux and the variance of the solid volume fraction, based on the modification of the existing correlations. The new correlation with the actual drift flux and variance of the solid volume fraction is in excellent agreement with the DNS results.

3. (iii) The existing models of the subgrid drift velocity and variance of the solid volume fraction have been assessed in an a priori manner. A new model for the variance of the solid volume fraction (3.23) is proposed by combining the gradient model and the scale-similarity model. The good performances of this model for the variance of the solid volume fraction (3.23) and the model of Schneiderbauer (2017) for the drift velocity (3.24) with the variance of the solid volume fraction computed with our model have been demonstrated.

4. (iv) All mesoscale models considered (i.e. scale-similarity model, gradient model, Schneiderbauer model and our model) can predict well the filtered drag with comparable accuracy (our model is slightly better), and are superior to the homogeneous drag model for the inhomogeneous system.

In addition, in the appendices, we show that our correlation for the mesoscale drag correction using the actual subgrid quantities with the same parameter values is applicable to the different computational domain sizes. All mesoscale drag models can predict well the filtered drag and are superior to the homogeneous drag model for a small-scale system. Nevertheless, since the cluster structure in this small-scale fluidization system is significantly different from the large system studied, the values of the constants in the scale-similarity model need to be adjusted for better agreement with the DNS data. For the other models, the same parameter values obtained from the large-scale system are used for the small-scale system, and thus, this may be deemed as a posteriori evaluation. Nevertheless, more a posteriori evaluations are certainly necessary to fully evaluate the models. We plan to conduct further a posteriori tests by comparing to the experiments on the clustered particle-laden flows, after the complete filtered two-fluid models (including the turbulence models and particle stress models) are established.

We believe that the present study provides important insights on the modelling of the drag force of the clustered particle-laden flows and the new model is a useful complement to existing mesoscale drag models.

In the present work we were only concerned with the effects of particle distribution heterogeneity on the drag, but the effects of heterogeneity on particle collision stress and Reynolds stresses are also important for the accurate simulation of clustered particle-laden flows, which are good subjects of future studies.