General Drag Correlations for Particle-Fluid System

Abstract

Particle-fluid flows are commonly encountered in industrial applications. It is of great importance to understand the fundamentals governing the behavior of such a flow system for better process design, control, and optimization. Generally, the particle-fluid flow behavior is strongly influenced by the interaction forces between fluid and particles. Among the various kinds of particle-fluid interaction forces, the drag force is the most essential. This chapter reviews the modeling of drag force for particle-fluid systems: from single particle to multiple particles, monosize to multisize, spherical to nonspherical, and Newtonian fluid to non-Newtonian fluid. Typical drag correlations in the literature are compared and assessed in terms of physical meaning, consistency, and generality.

1. Introduction

Particle-fluid flows are commonly encountered in industrial applications. It is of great importance to understand the fundamentals governing the behavior of such a flow system for better process design, control, and optimization. The flow behaviors of particles and fluid are strongly influenced by their interaction forces. Thus, it is critical to model the particle-fluid interaction forces accurately when simulating particle-fluid flows. The particle-fluid interaction forces include the pressure gradient force (or buoyancy force), drag force, virtual mass force, Basset force, and lift forces. The drag force, which is usually the dominant force in many particle-fluid flow systems, is undoubtedly the most critical and most studied.

Numerous efforts have been made to quantify the drag force using experimental or numerical methods in the past decades. Early studies in this area were mainly conducted by experiments [1, 2, 3, 4]. Due to the limitations of techniques, the experimental conditions are difficult to control. Therefore, the drag force models proposed in these studies somewhat lack generality and consistency. Nevertheless, these pioneer studies [1, 2, 3, 4] provide a solid foundation for subsequent research. With the rapid development of computational technology, various numerical methods have become attractive for studying the fluid-particle interaction on a subparticle scale, such as direct numerical simulation (DNS) and Lattice-Boltzmann (LB) model. With these numerical approaches, various fluid-particle systems can be studied under well-controlled conditions, considering more important and complicated factors that affect fluid-particle drag force. Nowadays, the studies on the fluid-particle interaction associated with the development of the drag force model have been extended extensively: from single particle to multiple particles, monosize to multisize, spherical particle to nonspherical particle system, and Newtonian fluid to non-Newtonian fluid. This chapter will review the modeling of particle-fluid drag force from these perspectives.

2. Mathematical model

The flow of a fluid is governed by the momentum and mass conservation equations when ignoring the compressible and viscous heat dissipation effect:

$$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathbf{u}\right)=0$$E1

$$\frac{\partial \left(\rho \mathbf{u}\right)}{\partial t}+\nabla ·\left(\rho \mathbf{uu}\right)=\nabla p+\nabla ·\left(\mathit{\rho \nu }\nabla \mathbf{u}\right)$$E2

where $\rho$, $p$, $\mathbf{u}$, and ${\nu }_f$ are the fluid density, pressure, velocity, and fluid kinematic viscosity, respectively.

In a traditional computational fluid dynamics (CFD) method, Eqs. (1) and (2) with initial and boundary conditions are often solved by the finite volume method (FVM) or finite element method (FEM). However, due to the difficulty of grid generation and boundary treatment in the DNS, the FVM and FEM are not widely used in studying the interaction between fluid and particles on a subparticle scale. Instead of FVM and FEM, the LB model is more widely used to study the drag force in particle-fluid systems. Its detailed description can be found elsewhere [5]. For brevity, the following only briefly introduces the LB model used.

The equations of the LB model are solved in two steps. First is the collision step:

$$f_i^{+}\left(\mathbf{x},t\right)=f_i\left(\mathbf{x},t\right)-{\left[{\mathbf{M}}^{-1}{\mathbf{S}}_f·\left(\mathbf{m}-{\mathbf{m}}^{\mathit{eq}}\right)\left(\mathbf{x},t\right)\right]}_i$$E3

which is followed by a propagation step:

$$f_i\left(\mathit{x}+{\mathit{c}}_{\mathit{i}}\mathit{\delta t},t+\mathit{\delta t}\right)=f_i^{+}\left(\mathit{x},t\right)$$E4

where $\mathbf{m}=\mathbf{Mf}$ and $\mathbf{f}$ is the column vector of $f_i$ is the velocity distribution function at the lattice node $\mathbf{x}$ and the time $t$ with the discrete velocity vector ${\mathbf{c}}_{\mathit{i}}$ in the $i$th direction. The number of directions of discrete velocity vectors ${\mathbf{c}}_{\mathit{i}}$ depends on the velocity model used. For example, provided that the fluid particle can move in the 19 directions in a three-dimensional case, as shown in Figure 1, the velocity model is referred to as a D3Q19 scheme. For the D3Q19 scheme, ${\mathbf{c}}_{\mathit{i}}$ is given as

Figure 1.

Schematic illustrations of (a) D3Q19 LB method, and (b) bounce-back boundary treatment.

$${\mathbf{c}}_{\mathit{i}}=c\left[\begin{array}{ccccccccccccccccccc}0& 1& -1& 0& 0& 0& 0& 1& -1& 1& -1& 1& -1& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 1& -1& 0& 0& 1& 1& -1& -1& 0& 0& 0& 0& 1& -1& 1& -1\\ 0& 0& 0& 0& 0& 1& -1& 0& 0& 0& 0& 1& 1& -1& -1& 1& 1& -1& -1\end{array}\right]$$E5

where $c=\frac{\mathit{\delta x}}{\mathit{\delta t}}$ is the lattice speed, $\mathit{\delta x}$ is the lattice length, and $\mathit{\delta t}$ is the time step. Other velocity models are also available, such as D3Q15 and D3Q27, and their details can be found elsewhere [5, 6].

The transfer matrix $\mathbf{M}$ in Eq. (3) defines the transformation of the distribution function to the moment space, which can be chosen to be the same as that of D’Humieres et al. [6], as given in Eq. (6). Note that ${\mathbf{S}}_f$ in Eq. (7) is a diagonal matrix. Provided that all the elements of ${\mathbf{S}}_f$ equal the same value as $\tau$, Eq. (3) can be reduced to the well-known LBGK scheme, $f_i^{+}\left(\mathbf{x},t\right)=f_i\left(\mathbf{x},t\right)-\frac{f_i\left(\mathbf{x},t\right)-f_i^{\left(\mathit{eq}\right)}\left(\mathbf{x},t\right)}{\tau }$. D’Humieres et al. [6] gives an optimized value of ${\mathbf{S}}_f$, where $s_0^f=s_3^f=s_5^f=s_7^f=0$, $s_1^f=1.19$, $s_2^f=s_{10}^f=s_{12}^f=1.4$, $s_{16}^f=s_{17}^f=s_{18}^f=1.98$, $s_9^f=s_{11}^f=s_{13}^f=s_{14}^f=s_{15}^f=s_{\nu }$, $s_4^f=s_6^f=s_8^f=s_q$, and $s_q=\frac{8\left(2-s_{\nu }\right)}{\left(8-s_{\nu }\right)}$. $s_q$ are set to $s_q=\frac{8\left(2-s_{\nu }\right)}{\left(8-s_{\nu }\right)}$ to satisfy the nonslip boundary condition. It should be noted that $s_{\nu }$ has a relationship with the kinematic viscosity $\nu$, as given by Eq. (8).

$$\mathbf{M}=\left[\begin{array}{ccccccccccccccccccc}1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1\\ -30& -11& -11& -11& -11& -11& -11& 8& 8& 8& 8& 8& 8& 8& 8& 8& 8& 8& 8\\ 12& -4& -4& -4& -4& -4& -4& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1\\ 0& 1& -1& 0& 0& 0& 0& 1& -1& 1& -1& 1& -1& 1& -1& 0& 0& 0& 0\\ 0& -4& 4& 0& 0& 0& 0& 1& -1& 1& -1& 1& -1& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 1& -1& 0& 0& 1& 1& -1& -1& 0& 0& 0& 0& 1& -1& 1& -1\\ 0& 0& 0& -4& 4& 0& 0& 1& 1& -1& -1& 0& 0& 0& 0& 1& -1& 1& -1\\ 0& 0& 0& 0& 0& 1& -1& 0& 0& 0& 0& 1& 1& -1& -1& 1& 1& -1& -1\\ 0& 0& 0& 0& 0& -4& 4& 0& 0& 0& 0& 1& 1& -1& -1& 1& 1& -1& -1\\ 0& 2& 2& -1& -1& -1& -1& 1& 1& 1& 1& 1& 1& 1& 1& -2& -2& -2& -2\\ 0& -4& -4& 2& 2& 2& 2& 1& 1& 1& 1& 1& 1& 1& 1& -2& -2& -2& -2\\ 0& 0& 0& 1& 1& -1& -1& 1& 1& 1& 1& -1& -1& -1& -1& 0& 0& 0& 0\\ 0& 0& 0& -2& -2& 2& 2& 1& 1& 1& 1& -1& -1& -1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& -1& -1& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& -1& -1& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& -1& -1& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& -1& 1& -1& -1& 1& -1& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& -1& -1& 1& 1& 0& 0& 0& 0& 1& -1& 1& -1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 1& -1& -1& -1& -1& 1& 1\end{array}\right]$$E6

$${\mathbf{S}}_f=\mathit{diag}\left(s_0^f,s_1^f,s_2^f,s_3^f,s_4^f,s_5^f,s_6^f,s_7^f,s_8^f,s_9^f,s_{10}^f,s_{11}^f,s_{12}^f,s_{13}^f,s_{14}^f,s_{15}^f,s_{16}^f,s_{17}^f,s_{18}^f\right)$$E7

$$\nu =\frac{1}{3}\left(\frac{1}{s_{\nu }}-0.5\right)$$E8

${\mathbf{m}}^{\mathit{eq}}$ in Eq. (3) are the velocity equilibrium moments defined by Eq. (9). The $j_x$, $j_y,$ and $j_z$ are defined as $j_x=\rho u_x$, $j_y=\rho u_{y,}$ and $j_z=\rho u_z$. The constant ${\rho }_0$ is the mean fluid density of the system, which is set to unity for an incompressible fluid.

$${\mathbf{m}}_f^{\mathit{eq}}={\left(\begin{array}{c}\rho ,-11\rho +\frac{19\left(j_x^2+j_y^2+j_z^2\right)}{{\rho }_0},3\rho -\frac{11\left(j_x^2+j_y^2+j_z^2\right)}{2{\rho }_0},j_x,-\frac{2j_x}{3},j_y,-\frac{2j_y}{3},j_z,-\frac{2j_z}{3},\\ \frac{\left(2j_x^2-j_y^2-j_z^2\right)}{{\rho }_0},-\frac{\left(j_x^2-j_y^2-j_z^2\right)}{2{\rho }_0},\frac{\left(j_y^2-j_z^2\right)}{{\rho }_0},\frac{\left(j_y^2-j_z^2\right)}{2{\rho }_0},\frac{j_x{j}_y}{{\rho }_0},\frac{j_y{j}_z}{{\rho }_0},\frac{x{j}_z}{{\rho }_0},0,0,0\end{array}\right)}^T$$E9

When the effect of fluid rheology is considered, ν is dependent on the shear rate for a non-Newtonian fluid other than a constant for a Newtonian fluid. It is described here by the power-law model, which is suitable for a wide range of non-Newtonian fluids [7]:

$$\nu ={\nu }_0{\dot{e}}^{n-1}={\nu }_0{\left(2e_{\mathit{\alpha \beta }}{e}_{\mathit{\alpha \beta }}\right)}^{\left(n-1\right)/2}$$E10

where ${\nu }_0$ is the flow consistency index, $\dot{e}$ denotes the shear rate, and n is the power-law index. The fluid shows shear-thickening behavior when n > 1, shear-thinning behavior when n < 1, and Newtonian behavior when n = 1. In Eq. (10), Einstein’s summation convention is applied and $e_{\mathit{\alpha \beta }}$ represents an element of the tensor $\dot{e}$ at the position αβ, calculated by [8]

$$e_{\mathit{\alpha \beta }}=\frac{3}{2\mathit{\rho \tau }}\sum _{i=0}{f}_i^{\left(1\right)}{c}_{\mathit{i\alpha }}{c}_{\mathit{i\beta }}$$E11

where α and β denote spatial indices, $c_{\mathit{i\alpha }}$ and $c_{\mathit{i\beta }}$ are the elements of the vector ${\mathbf{c}}_i$ at positions α and β, and $f_i^{\left(1\right)}=f_i-f_i^{\mathrm{eq}}$ is the non-equilibrium part of the distribution function.

Through the Chapman-Enskog multiple-scale expansion, Eqs. (3) and (4) can be recovered to the Navier-Stokes equations in a low Mach number limit. The macroscopic variables such as density $\rho$, momentum density $\rho \mathbf{u}$, and stress $\mathit{\sigma }\mathit{,}$can be obtained through the moments of the velocity and temperature distributions:

$$\rho =\sum _i{f}_i,\rho \mathit{u}=\sum _i{f}_i{\mathit{c}}_{\mathit{i}},\mathit{\sigma }=\sum _i{f}_i{\mathit{c}}_i{\mathit{c}}_i$$E12

The force exerted on a solid particle is calculated by the momentum exchange method proposed by Ladd [9], given by

$$f_i\left({\mathbf{x}}_b,t+0.5\right)=2\left(f_i^{+}\left({\mathbf{x}}_f,t\right)-f_{\overline{i}}^{+}\left({\mathbf{x}}_s,t\right)\right){\mathit{c}}_i$$E13

where x_f, x_s, and x_b denote the nodes in the fluid region and solid region and at the boundary for a stationary boundary on a uniform lattice, and the subscript $\overline{i}$ indicates the opposite direction of i. The curved wall of spherical solid particles is treated by the halfway boundary scheme (Figure 1b), which sets the wall node at the mid-point between fluid and solid nodes. The selection of the halfway boundary treatment has various advantages. First, it is the simplest and most efficient boundary treatment. Secondly, compared with interpolation schemes, the halfway boundary scheme is very space-saving in modeling particles that are placed closely (e.g., in a packed bed) [10]. Thirdly, the high numerical stability of the halfway boundary scheme ensures simulation reliability under various complicated flow conditions. Lastly, the halfway boundary scheme is a second-order scheme; its accuracy in calculating forces can be as good as other interpolation schemes and better than immersed boundary schemes once a calibration step is applied [11].

3. Model validation

It is necessary to verify the validity of the mathematical model before its application for numerical experiments. The LB model has proven to be valid for studying the fluid-particle interactions under different conditions involving a single particle, two interactive particles, and packed beds of spherical and nonspherical particles over a wide range of conditions [12, 13, 14, 15]. For example, the LB prediction of drag coefficient $C_{D0}=\frac{8f_d}{\mathit{\pi \rho }{d_e}^2{U_0}^2}$ for an isolated particle in the Newtonian fluid has been compared with the experimental measurements of Schlichting and Gersten [16], and with those calculated by three correlations given by Stokes [17], Dallavalle [4], and Clift et al. [18], respectively. Figure 2a shows good agreement between the simulations and experiments, including the three correlations. Specifically, in the intermediate flow region of Re = 100–1000, the drag coefficient of the sphere for the velocity inlet and outflow conditions under open BCs falls closer to the experimental data than the one under periodic boundary conditions. This may be due to the nonuniformity in the flow far away from the sphere for high values of Re, the state of which can be different from that of the experiments. Drag coefficient in the flow of power-law fluid has also been compared against the measurements by Chhabra [19] and Peden and Luo [20] over 5 ≤ Re ≤ 100 and 0.6 ≤ n ≤ 1, as shown in Figure 2b. In the comparison, the numerical data obtained by Dhole et al. [21] using the FVM are also considered. Figure 2 shows that both the LB model and FVM can well predict the drag coefficient on particles in the power-law fluid flows at different Re and n. The prediction errors are less than 4.2% for the present LB model with a calibration step and 7.6% for the FVM.

Figure 2.

Comparison between simulated and measured drag coefficients in (a) Newtonian, and (b) non-Newtonian fluids.

4. Drag force for single-particle system

The fluid flow past an isolated particle is one of the most basic flow phenomena. It can be seen as a special case of a particle-fluid system when the particle volume fraction is close to zero. Also, it is very natural to start from a single-particle system to investigate the drag force in a complicated particle-fluid system. In the past decades, extensive efforts have been made to study various physical characteristics (e.g., drag, wake, shear, heat transfer, and vortex shedding) of the flow past an isolated particle. Tiwari et al. [22, 23] have given a comprehensive review on this topic, and interested readers can refer to the review for details.

For an isolated particle in fluid flows, the particle-fluid drag force has been well established, which is expressed as [4]

$${\mathbf{f}}_d=\frac{1}{8}{C}_{D0}\pi {\rho }_f{d}_p^2\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|\left({\mathbf{u}}_f-{\mathbf{u}}_p\right)$$E14

where ${\mathbf{f}}_d$ is the drag force, $C_{D0}$ is the drag coefficient for a single particle in fluid flows, which is determined by the Reynolds number $\mathit{Re}=\frac{{\rho }_f{d}_p\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|}{{\mu }_f}$, ${\rho }_f$ is the fluid density, $d_p$ is the particle diameter, ${\mu }_f$ is the fluid viscosity, and ${\mathbf{u}}_f$ and ${\mathbf{u}}_p$ are the fluid and particle velocities, respectively.

Figure 3 shows the relationship between $C_{D0}$ and $\mathit{Re}$ in different flow regimes established based on experimental and numerical data. It can be seen that the effect of $\mathit{Re}$ on $C_{D0}$ varies in different flow regimes. Generally, the flow regimes can be classified into nine groups: laminar regime or axisymmetric wake regime (regime I), planar symmetric wake regime (regime II), vortex shedding regime (regime III), separating vortex tubes regime (regime IV), subcritical regime (regime V), critical regime (regime VI), supercritical regime (regime VII), and transcritical regime (regime VIII). It should be noted that the regimes I–IV may also be called the transition regime. In the laminar regime, the drag coefficient $C_{D0}$ decreases significantly as $\mathit{Re}$ increases. In contrast, in the transition regime (regimes I–-IV), the drop of $C_{D0}$ decreases as $\mathit{Re}$ increases and $C_{D0}$ approaches a constant in the subcritical regime. Then, $C_{D0}$ experiences a rapid drop in the critical regime but increases again and reaches a near constant in the subsequent supercritical and transcritical regimes.

Figure 3.

Drag coefficient for the flows past a stationary sphere in different flow regimes based on experimental and numerical data [23].

Due to the complicated effect of $\mathit{Re}$ on $C_{D0}$, it is still challenging to formulate a drag coefficient correlation that can accurately describe the complicated variation of $C_{D0}$ in different flow regimes. However, a general drag coefficient can be formulated to meet the needs of engineering applications if the following requirements are satisfied:

1. $C_{D0}$ decreases in the laminar region.

2. The drop of $C_{D0}$ decreases in the transition regime.

3. $C_{D0}$ approaches a positive constant value when $\mathit{Re}$ is relatively large.

Table 1 summarizes some of the well-known drag coefficient correlations. Notably, all these drag coefficient correlations can meet the above requirements.

References	Drag correlation $C_{D0}$
Schiller and Nauman (1935) [24]	$C_{D0}=\left\{\begin{array}{c}\frac{24}{\mathit{Re}}\left(1+0.15{\mathit{Re}}^{0.687}\right)\mathit{Re}<1000\\ 0.44\mathit{Re}>1000\end{array}\right.$
Dallavalle (1948) [4]	$C_{D0}={\left(0.63+\frac{4.8}{\sqrt{\mathit{Re}}}\right)}^2$
Turton and Levenspiel (1986) [25]	$C_{D0}=\frac{24}{\mathit{Re}}\left(1+0.173{\mathit{Re}}^{0.657}\right)+\frac{0.413}{1+16300{\mathit{Re}}^{-1.09}}$

Table 1.

Drag correlations for an isolated particle in fluid flows.

5. From single-particle system to multiparticle system

Extension from an isolated particle system to a multi-particle system is a big challenge due to the complexity of the particulate system. The presence of other particles reduces the space for fluid and generates a sharp fluid velocity gradient, yielding increased shear stress on particle surfaces. The drag force enhancement is closely associated with particle configuration, particle-fluid slip velocity, and particle and fluid properties.

Generally, the methods to formulate drag correlations for multi-particle systems can be grouped into Ergun and Wen-Yu types. The Ergun-type model, also called the capillary model, focuses on the flow resistance of the whole system, exemplified by the Ergun Eq. (1). This kind of model originates from the idea of treating the void space in a porous medium as a bunch of tortuous conduits, as illustrated in Figure 4. Using equivalent tubes instead of the void space of a packed bed, one can model the total pressure drop across the whole system by considering the contributions of viscous force and inertial force. Correspondingly, its friction factor $f=\frac{d\left(∆P\right)}{\mathit{L\rho }{\left|{\mathbf{U}}_{\mathbf{s}}\right|}^2}\left(\frac{{\epsilon }^3}{1-\epsilon }\right)$ falls into the form of $f=\frac{A}{\mathit{Re}}+B$. The Ergun-type model is an easier way to formulate correlations by fitting the data in the creeping flow regime or the extremely high Reynolds flow regime. However, as pointed out previously [26], because of the original idea of treating a multiparticle system as tortuous conduits, the Ergun-type correlations are more suitable for densely packed beds but usually have a discontinuity at high porosity. For example, when particles are immersed in a highly dilute system, the drag force determined by the Ergun-type model is usually inconsistent with the value given by the drag correlation for a single particle, which is apparently against the real physical world. This defect may cause problems in some situations, such as CFD-DEM simulations, and limit the application range of the Ergun-type model. The Wen-Yu-type model, also called the submerge object model, focuses on the enhanced drag force on a particle due to neighboring particles, exemplified by the Wen-Yu Eq. (2). Such a model is more straightforward as it directly considers the effect of neighboring particles on the drag force, as illustrated in Figure 4. In the model, the enhancement of drag force on particles in a multi-particle system can be expressed as a ratio to the drag force on a single particle in an unhindered environment. Correspondingly, the drag force in the Wen-Yu type can be written in ${\mathbf{f}}_d=\frac{1}{8}{C}_D\pi {\rho }_f{d}_p^2\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|\left({\mathbf{u}}_f-{\mathbf{u}}_p\right)$, where $C_D=C_{D0}{{\epsilon }_f}^{2-\chi }$. Unlike Ergun-type models, Wen-Yu-type models naturally ensure the consistency of drag force between the predictions by the correlation for a multi-particle system under dilute conditions and the correlation for a single-particle system, leading to its popularity in CFD-DEM simulations.

Figure 4.

Schematic representation of idealization of Ergun-type model (a), and Wen-Yu-type model (b).

Irrespective of the different origins of Ergun-type and Wen-Yu-type models, the drag correlations proposed by either way can be written in the following equation:

$${\mathbf{f}}_d=\frac{1}{8}{C}_D\pi {\rho }_f{d}_p^2\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|\left({\mathbf{u}}_f-{\mathbf{u}}_p\right)$$E15

where $C_D$ is a coefficient relying on $\mathit{Re}=\frac{{\rho }_f{d}_p{\epsilon }_f\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|}{{\mu }_f}$ and ${\epsilon }_f$ is the porosity. $C_D$ is difficult to derive theoretically for a multiparticle system. Formulating an empirical correlation of $C_D$ should meet the following requirements for physical consistency and generality:

1. $C_D$ should be able to reduce to $C_{D0}$ when ${\epsilon }_f$ is approaching 1;

2. $C_D$ should be continuous over the whole application range.

3. $C_D$ should always be larger than $C_{D0}$ when ${\epsilon }_f$ is smaller than 1;

Note that the requirement (2-1) must be satisfied in a CFD-DEM simulation. Otherwise, inconsistency may occur in the dilute regime. The requirement (2-1) can be written in $C_{D0}=\underset{\epsilon \to 1}{lim}{C}_D$, which should satisfy the requirements (1-1) to (1-3). This is important for evaluating different drag correlations when new drag correlations are considered for a single particle. The requirement (2-3) ensures that $C_D$ is physically meaningful because neighboring particles reduce the space for fluid and generate a higher particle-fluid velocity in the surrounding region under the same $\mathit{Re}$, increasing the drag force.

Table 2 summarizes the correlations of drag coefficient $C_D$ reported for multi-particle systems in recent years. For easy comparison, the value of $C_D$ when ${\epsilon }_f=1$ is given to assess if $\underset{\epsilon \to 1}{lim}{C}_D$ satisfies the requirements from (1-1) to (1-3). Table 2 also indicates whether the correlations meet the requirements from (2-1) to (2-3). Note that the requirement (2-3) is not easy to be assessed directly. The values of void function $\chi =2-{log}_{{\epsilon }_f}\frac{C_D}{C_{D0}}$ (Figure 5) are calculated from different drag correlations under different conditions to solve this problem. The value of $\chi$ in Figure 5 must be positive to satisfy the requirement (2-3). Otherwise, $C_D$ would be smaller than $C_{D0}$.

References	Drag coefficient $C_D$	$C_D$ when ${\epsilon }_f\to 1$	Req (2-1)	Req (2-2)	Req (2-3)
Ergun (1952) [1]	$C_D=\frac{4}{3}\left(150\frac{1-{\epsilon }_f}{\mathit{Re}}+1.75\right)$	$C_{D,{\epsilon }_f\to 1}=2.33$	$\times$	$\surd$	$\times$
Wen and Yu (1966) [2]	$C_D=C_{D0}{{\epsilon }_f}^{2-\chi }$, $C_{D0}=\left\{\begin{array}{c}\frac{24}{\mathit{Re}}\left(1+0.15{\mathit{Re}}^{0.687}\right)\mathit{Re}<1000\\ 0.44\mathit{Re}>1000\end{array}\right.$, $\chi =3.65$	$C_{D,{\epsilon }_f\to 1}=\left\{\begin{array}{c}\frac{24}{\mathit{Re}}\left(1+0.15{\mathit{Re}}^{0.687}\right)\mathit{Re}<1000\\ 0.44\mathit{Re}>1000\end{array}\right.$	$\surd$	$\surd$	$\surd$
Gibilaro et al. (1985) [27]	$C_D=C_{D0}{{\epsilon }_f}^{2-\chi }$, $C_{D0}=\frac{23}{\mathit{Re}}+0.448$, $\chi =3.8$	$C_{D,{\epsilon }_f\to 1}=\frac{23}{\mathit{Re}}+0.448$	$\surd$	$\surd$	$\surd$
Gidaspow (1994) [28]	$C_D=\left\{\begin{array}{c}{C}_{D0}{{\epsilon }_f}^{2-\chi }{\epsilon }_f\ge 0.8\\ \frac{4}{3}\left(150\frac{1-{\epsilon }_f}{\mathit{Re}}+1.75\right){\epsilon }_f<0.8\end{array}\right.$, $C_{D0}=\left\{\begin{array}{c}\frac{24}{\mathit{Re}}\left(1+0.15{\mathit{Re}}^{0.687}\right)\mathit{Re}<1000\\ 0.44\mathit{Re}>1000\end{array}\right.$, $\chi =3.65$	$C_{D,{\epsilon }_f\to 1}=\left\{\begin{array}{c}\frac{24}{\mathit{Re}}\left(1+0.15{\mathit{Re}}^{0.687}\right)\mathit{Re}<1000\\ 0.44\mathit{Re}>1000\end{array}\right.$	$\surd$	$\times$	$\surd$
Di Felice (1994) [29]	$C_D=C_{D0}{{\epsilon }_f}^{2-\chi }$, $C_{D0}={\left(0.63+\frac{4.8}{\sqrt{\mathit{Re}}}\right)}^2$, $\chi =3.7-0.65\mathit{exp}\left[-\frac{{\left(1.5-log\mathit{Re}\right)}^2}{2}\right]$	$C_{D,{\epsilon }_f\to 1}={\left(0.63+\frac{4.8}{\sqrt{\mathit{Re}}}\right)}^2$	$\surd$	$\surd$	$\surd$
Van der Hoef et al. (2005) [30]	$C_D=\frac{4}{3}\left(180\frac{1-{\epsilon }_f}{\mathit{Re}}+18\frac{{\epsilon }_f^4\left(1+1.5\sqrt{1-{\epsilon }_f}\right)}{\mathit{Re}}\right)$	$C_{D,{\epsilon }_f\to 1}=\frac{24}{\mathit{Re}}$	$\times$	$\surd$	$\times$
Benyahia et al. (2006) [31]	$C_D=\frac{12{\epsilon }_f^{2}F}{{\mathit{Re}}_{\ast }}$, where $F=\left\{\begin{array}{c}1+\frac{3}{8}{\mathit{Re}}_{\ast }{\epsilon }_f\ge 0.99\mathit{\text{and}}{\mathit{Re}}_{\ast }\le \frac{F_2-1}{3/8-F_3}\\ F_0+F_1{{\mathit{Re}}_{\ast }}^2{\epsilon }_f<0.99\mathit{\text{and}}{\mathit{Re}}_{\ast }\le \frac{F_3+\sqrt{F_3^2-4F_1\left(F_0-F_2\right)}}{2F_1}\\ F_2+F_3{\mathit{Re}}_{\ast }\mathit{\text{others}}\end{array}\right.$,; $F_0=\left\{\begin{array}{c}\left(1-w\right)\left[\frac{1+3\sqrt{\left(1-{\epsilon }_f\right)/2}+\left(135/64\right)\left(1-{\epsilon }_f\right)ln\left(1-{\epsilon }_f\right)+17.4\left(1-{\epsilon }_f\right)}{1+0.681\left(1-{\epsilon }_f\right)-8.48{\left(1-{\epsilon }_f\right)}^2+8.16{\left(1-{\epsilon }_f\right)}^3}\right]+w\left[10\frac{1-{\epsilon }_f}{{\epsilon }_f^3}\right]0.6<{\epsilon }_f<0.99\\ 10\frac{1-{\epsilon }_f}{{\epsilon }_f^3}{\epsilon }_f\le 0.6\end{array}\right.$,; $F_1=\left\{\begin{array}{c}\sqrt{2/\left(1-{\epsilon }_f\right)}/400.9\le {\epsilon }_f<0.99\\ 0.11+0.00051exp\left[11.6\left(1-{\epsilon }_f\right)\right]{\epsilon }_f<0.9\end{array}\right.$,; $F_2=\left\{\begin{array}{c}\left(1-w\right)\left[\frac{1+3\sqrt{\left(1-{\epsilon }_f\right)/2}+\left(135/64\right)\left(1-{\epsilon }_f\right)ln\left(1-{\epsilon }_f\right)+17.89\left(1-{\epsilon }_f\right)}{1+0.681\left(1-{\epsilon }_f\right)-11.03{\left(1-{\epsilon }_f\right)}^2+15.41{\left(1-{\epsilon }_f\right)}^3}\right]+w\left[10\frac{1-{\epsilon }_f}{{\epsilon }_f^3}\right]{\epsilon }_f>0.6\\ 10\frac{1-{\epsilon }_f}{{\epsilon }_f^3}{\epsilon }_f\le 0.6\end{array}\right.$,; $F_3=\left\{\begin{array}{c}0.9351\left(1-{\epsilon }_f\right)+0.03667{\epsilon }_f>0.9047\\ 0.0673+0.212\left(1-{\epsilon }_f\right)+0.0232{\epsilon }_f^{-5}{\epsilon }_f\le 0.9047\end{array}\right.$,; ${\mathit{Re}}_{\ast }=\frac{{\rho }_f{\epsilon }_f{d}_p\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|}{2{\mu }_f}$,; $w=e^{\left(-10\left({\epsilon }_f-0.6\right)/\left(1-{\epsilon }_f\right)\right)}$	$C_{C_{D,{\epsilon }_f\to 1}}=\left\{\begin{array}{c}\frac{12}{{\mathit{Re}}_{\ast }}\left(1+0.375{\mathit{Re}}_{\ast }\right){\mathit{Re}}_{\ast }\le \frac{F_2-1}{3/8-F_3}\\ \frac{12}{{\mathit{Re}}_{\ast }}\left(1+0.03667{\mathit{Re}}_{\ast }\right)\mathit{Re}>\frac{F_2-1}{3/8-F_3}\end{array}\right.$	$\times$	$\times$	$\surd$
Mazzei and Lettieri (2007) [32]	$C_D=C_{D0}{{\epsilon }_f}^{2-\chi }$, $C_{D0}={\left(0.63+\frac{4.8}{\sqrt{\mathit{Re}}}\right)}^2$, $\chi =1-2\frac{4.8+0.42{\mathit{Re}}_{\ast }^{3/4}}{1+0.175{\mathit{Re}}_{\ast }^{3/4}}$ and ${\mathit{Re}}_{\ast }=\frac{{\rho }_f{d}_p\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|}{{\mu }_f{\epsilon }_f^{n-1}}$	$C_{D,{\epsilon }_f\to 1}={\left(0.63+\frac{4.8}{\sqrt{\mathit{Re}}}\right)}^2$	$\surd$	$\surd$	$\surd$
Beetstra et al. (2007) [26]	$C_D=\frac{4}{3}\left(180\frac{1-{\epsilon }_f}{\mathit{Re}}+18\frac{{\epsilon }_f^4\left(1+1.5\sqrt{1-{\epsilon }_f}\right)}{\mathit{Re}}+0.31\frac{{\epsilon }_f^{-1}+3\left(1-{\epsilon }_f\right){\epsilon }_f+8.4{\mathit{Re}}^{-0.343}}{1+{10}^{3\left(1-{\epsilon }_f\right)}{\mathit{Re}}^{-0.5-2\left(1-{\epsilon }_f\right)}}\right)$	$C_{C_{D,{\epsilon }_f\to 1}}=\frac{24}{\mathit{Re}}+0.413\frac{1+8.4{\mathit{Re}}^{-0.343}}{1+{\mathit{Re}}^{-0.5}}$	$\surd$	$\surd$	$\surd$
Cello et al. (2010) [33]	$C_D=\frac{24{{\epsilon }_f}^2}{\mathit{Re}}\left[K_1+K_2{\epsilon }_f^4+K_3\left(1-{\epsilon }_f^4\right)\right]$, $K_0=\frac{1-{\epsilon }_f}{1+3{\epsilon }_f}$, $K_1=\frac{1+128K_0+715K_0^2}{{\epsilon }_f\left(1+49.5K_0\right)}$, $K_2=\frac{1+0.13\mathit{Re}+0.000666{\mathit{Re}}^2}{1+0.0342\mathit{Re}+0.00000692{\mathit{Re}}^2}-1$,; $K_3=\left(\frac{2{\mathit{Re}}^2}{1+\mathit{Re}}\right)\left(\frac{-410{\epsilon }_f+9.2\times {10}^7{K}_0^{20}\mathit{Re}+1900{\epsilon }_f^2-0.066\mathit{Re}}{6600{\epsilon }_f+0.000429\mathit{Re}-43000{\epsilon }_f^2-0.000131{\mathit{Re}}^2+73800{\epsilon }_f^3}\right)$	$C_{D,{\epsilon }_f\to 1}=\frac{24}{\mathit{Re}}\left(1+K_2\right)$	$\times$	$\surd$	$\times$
Tenneti et al. (2011) [34]	$C_D=C_{D0}+\frac{4}{3}\left(104.58\frac{1-{\epsilon }_f}{\mathit{Re}}+8.64\frac{{\left(1-{\epsilon }_f\right)}^{1/3}}{\mathit{Re}{\epsilon }_f}+18{\epsilon }_f^3{\left(1-{\epsilon }_f\right)}^3\left(0.95+\frac{0.61{\left(1-{\epsilon }_f\right)}^3}{{\epsilon }_f^2}\right)\right)$,; $C_{D0}=\left\{\begin{array}{c}\frac{24}{\mathit{Re}}\left(1+0.15{\mathit{Re}}^{0.687}\right)\mathit{Re}<1000\\ 0.44\mathit{Re}>1000\end{array}\right.$	$C_{D,{\epsilon }_f\to 1}=\left\{\begin{array}{c}\frac{24}{\mathit{Re}}\left(1+0.15{\mathit{Re}}^{0.687}\right)\mathit{Re}<1000\\ 0.44\mathit{Re}>1000\end{array}\right.$	$\surd$	$\times$	$\surd$
Rong et al. (2013) [14]	$C_D=C_{D0}{{\epsilon }_f}^{2-\chi }$, $C_{D0}={\left(0.63+\frac{4.8}{\sqrt{\mathit{Re}}}\right)}^2$, $\chi =2.65\left(1+{\epsilon }_f\right)-\left(5.3-3.5{\epsilon }_f\right){\epsilon }_f^{2}exp\left[-\frac{{\left(1.5-log\mathit{Re}\right)}^2}{2}\right]$	$C_{D,{\epsilon }_f\to 1}={\left(0.63+\frac{4.8}{\sqrt{\mathit{Re}}}\right)}^2$	$\surd$	$\surd$	$\surd$
Zaidi et al. (2014) [35]	$C_D=\left\{\begin{array}{c}\frac{4}{3}\left(180\frac{1-{\epsilon }_f}{\mathit{Re}{\epsilon }_f}+18\frac{{\epsilon }_f^4\left(1+1.5\sqrt{1-{\epsilon }_f}\right)}{\mathit{Re}}+\frac{0.612}{{\epsilon }_f^{0.7}}\right)\mathit{Re}\le 200\\ \frac{4}{3}\left(\frac{196.2{\left(1-{\epsilon }_f\right)}^{0.4}}{\mathit{Re}{\epsilon }_f^{-0.3}}+\frac{0.432}{{\epsilon }_f^{0.86}}\right)\mathit{Re}\ge 200\end{array}\right.$	$C_{D,{\epsilon }_f\to 1}=\left\{\begin{array}{c}\frac{24}{\mathit{Re}}+0.816\mathit{Re}\le 200\\ 0.576\mathit{Re}\ge 200\end{array}\right.$	$\times$	$\times$	$\times$
Bogner et al. (2015) [36]	$C_D=\frac{24}{\mathit{Re}{\epsilon }_f^{2.726}}\left[1.751+0.151{\mathit{Re}}^{0.684}-0.445{\left(1+\mathit{Re}\right)}^{1.04\left(1-{\epsilon }_f\right)}-0.16{\left(1+\mathit{Re}\right)}^{0.0003\left(1-{\epsilon }_f\right)}\right]$	$C_{D,{\epsilon }_f\to 1}=\frac{24}{\mathit{Re}}\left(1.146+0.151{\mathit{Re}}^{0.684}\right)$	$\times$	$\surd$	$\surd$
Tang et al. (2015) [37]	$C_D=\frac{4}{3}\left(180\frac{1-{\epsilon }_f}{\mathit{Re}}+18\frac{{\epsilon }_f^4\left(1+1.5\sqrt{1-{\epsilon }_f}\right)}{\mathit{Re}}+18{\epsilon }_f^2\left[0.11\left(1-{\epsilon }_f\right)\left(2-{\epsilon }_f\right)-\frac{0.00456}{{\epsilon }_f^4}+\left(0.169{\epsilon }_f+\frac{0.0644}{{\epsilon }_f^4}\right){\mathit{Re}}^{-0.343}\right]\right)$	$C_{D,{\epsilon }_f\to 1}=\frac{24}{\mathit{Re}}\left(1-0.0046\mathit{Re}+0.233{\mathit{Re}}^{0.657}\right)$	$\times$	$\surd$	$\times$
Zhou and Fan (2015) [38]	$C_D=F_0\frac{24{\epsilon }_f^2}{\mathit{Re}}+\left(0.256+1.41\left(1-{\epsilon }_f\right)-5.61{\left(1-{\epsilon }_f\right)}^2+6.04{\left(1-{\epsilon }_f\right)}^3\right)\frac{\left[{\epsilon }_f^{-1}+3{\epsilon }_f\left(1-{\epsilon }_f\right)+8.4{\mathit{Re}}^{-0.343}\right]}{{\epsilon }_f\left[1+{10}^{3\left(1-{\epsilon }_f\right)}{\mathit{Re}}^{-2.5+2{\epsilon }_f}\right]}$; where $F_0=\left\{\begin{array}{c}\frac{9.9\left(1-{\epsilon }_f\right)}{{\epsilon }_f^2}+{\epsilon }_f^3\left[1+3{\left(1-{\epsilon }_f\right)}^{0.6}\right]{\epsilon }_f>0.45\\ \frac{5.87sin\left[\frac{\pi }{2}{\left(\frac{1-{\epsilon }_f}{0.637}\right)}^{1.75}\right]}{{\epsilon }_f^2}0.363\le {\epsilon }_f\le 0.45\\ \frac{5.87\left(1-{\epsilon }_f\right)}{0.637{\epsilon }_f^2}{\epsilon }_f<0.363\end{array}\right.$	$C_{D,{\epsilon }_f\to 1}=\frac{24}{\mathit{Re}}+0.256\frac{1+8.4{\mathit{Re}}^{-0.343}}{1+{\mathit{Re}}^{-0.5}}$	$\surd$	$\times$	$\surd$
Sheikh and Qiu (2018) [39]	$C_D=\frac{24}{\mathit{Re}{\epsilon }_f^{0.5}}\left[3^{1.25\left(1-{\epsilon }_f^3\right)}+\mathit{Re}\left(0.2{\mathit{Re}}^{0.5{\epsilon }_f-1}-0.053{\epsilon }_f+0.073\right)\right]$	$C_{D,{\epsilon }_f\to 1}=\frac{24}{\mathit{Re}}\left[1+\mathit{Re}\left(0.2{\mathit{Re}}^{-0.5}+0.02\right)\right]$	$\surd$	$\surd$	$\surd$
Kravets et al. (2019) [40]	$C_D=240\frac{1-{\epsilon }_f}{\mathit{Re}}+24\frac{{\epsilon }_f^4\left(1+1.5\sqrt{1-{\epsilon }_f}\right)}{\mathit{Re}}+24{{\epsilon }_f}^2\left[0.1695\left(1-{\epsilon }_f\right)\left(2-{\epsilon }_f\right)-\frac{0.004321}{{\epsilon }_f^4}+\left(0.0719{\epsilon }_f+\frac{0.02169}{{\epsilon }_f^4}\right){\mathit{Re}}^{-0.2017}\right]$	$C_{D,{\epsilon }_f\to 1}=\frac{24}{\mathit{Re}}\left[1-0.00432\mathit{Re}+0.0936{\mathit{Re}}^{0.7983}\right]$	$\times$	$\surd$	$\times$

Table 2.

Summary of the drag correlations for multiparticle systems.

Figure 5.

Exponent $\chi$ as a function of porosity ${\epsilon }_f$ and Reynolds number $\mathit{Re}$ calculated from different drag correlations.

Table 2 shows that most of the drag correlations in Ergun-type cannot satisfy the requirement (2-1), thereby generating nonphysical results in a dilute regime. Combining different correlations was employed in the past to overcome this problem. For example, the Gidaspow correlation [28], adopted in many commercial CFD software such as ANSYS Fluent, combined the Ergun correlation under dense conditions and the Wen-Yu correlation under dilute conditions. A similar treatment is also taken by Benyahia et al. [31] and Zhou and Fan [38]. However, this treatment may cause apparent discontinuity at the switching point, as shown in Figure 5, which is against the requirement (2-2). Some other investigators attempted to formulate an Ergun-type correlation by adding additional items to guarantee continuity based on DNS or LB simulation data for randomly distributed particle systems generated with the Monte Carlo method. Typical work can be exemplified by Van der Hoef et al. [30], Beetstra et al. [26], Cello et al. [33], Tang et al. [37], and Kravets et al. [40]. Among all these Ergun-type correlations, only the correlations proposed by Beetstra et al. [26] and Sheikh and Qiu [39] can meet all the requirements. The correlations of Tenneti et al. [34] and Zhou and Fan [38] can fulfill the requirements (2-1) and (2-3). However, with $\underset{\epsilon \to 1}{lim}{C}_D$, the correlation of Zhou and Fan [38] cannot be reduced to a single particle. Specifically, the value of the $\underset{\epsilon \to 1}{lim}{C}_D$ is much smaller than the experimental value at a high $\mathit{Re}$. The void function value $\chi$ predicted by the correlation of Tenneti et al. [34] has an opposite trend of other correlations. Also, a discontinuous point exists at a relatively high $\mathit{Re}$, as shown in Figure 5.

Compared with Ergun-type correlations, Wen-Yu-type drag correlations inherently have the advantage of satisfying both requirements (2-1) and (2-2) due to the introduction of the voidage function $\chi$. In fact, all the Wen-Yu-type correlations in Table 2 can satisfy the requirements from (2-1) to (2-3). One of the most widely used Wen-Yu-type correlations is the one proposed by Di Felice [29], which has good performance in both dilute and dense regimes. However, its voidage function is too simple and ignores the effect of $\mathit{Re}$ and ${\epsilon }_f$. Later, Rong et al. [14] used the LB method to study the drag force in different packed beds with a wide range of porosity generated by the DEM simulations, and on this basis, modified the voidage function of Di Felice [29] to incorporate the effect of porosity. Compared with other formulations, for example, the correlation of Beetstra et al. [26], that of Rong et al. [14] has a simpler form and better performance in a broader range of applications. Note that the existing drag correlations were often established on experimental or numerical data obtained in low-to-intermediate flow regimes (Re < 1000). The same expression has been adopted to deal with the high Re regime in applications. Nevertheless, drag force correlations for mono-size particle systems have been firmly established with the extensive efforts from different investigators in the last decade.

6. From monosize to multisize particle systems

The correlations in Table 2 are all formulated for systems with monosize spheres. However, particle systems composed of multisize spheres are encountered in most applications. The drag correlation formulated for monosize spheres cannot be directly applied to multisize particle systems. Therefore, some investigators have studied the drag forces using LB or DNS simulations for different components in multisize particle systems generated with DEM simulations or Monte Carlo methods, formulating new correlations based on numerical data.

Generally, two main treatments are used to estimate the fluid-particle drag force in a mixture of particles. The first one is to estimate the drag force on each particle directly, calculate the mean drag forces on different components, and finally obtain the total fluid-particle drag force by summing the forces of all particles. Due to the lack of reliable drag correlations for mixtures, different investigators used this treatment in early studies, such as Feng and Yu [41] and Bokkers et al. [42]. It has proven to have a poor performance by Rong et al. [13]. Thus, investigators turned their sight to the second approach, which uses an opposite calculation path. It first estimates the total fluid-particle drag force and then distributes the total force among different components to obtain their mean drag forces according to a specific rule. In this approach, the total drag force is estimated by treating the particle mixture as a monosize particle system with a representative average diameter, that is, the Sauter mean diameter $⟨d_p⟩={\left[\sum _i{x}_i/d_i\right]}^{-1}$, where $x_i$ and $d_i$ are the volume fraction and diameter of component $i$. The distributing function can be formulated in different ways. However, it should meet some basic requirements for physical consistency and generality, which have been well-established by Rong et al. [13], given by:

1. $\frac{{\stackrel{\sim }{f}}_{d,i}}{{\stackrel{\sim }{f}}_d}=1$ when all $d_i$ are equal, where ${\stackrel{\sim }{f}}_{d,i}={\mathbf{f}}_{d,i}/3\pi {\mu }_f{d}_i{\epsilon }_f\left({\mathbf{u}}_f-{\mathbf{u}}_p\right)$ is the drag force for component $i$ normalized by the Stokes drag and ${\stackrel{\sim }{f}}_d$ is normalized drag force for the equivalent mono-size particle.

2. ${\stackrel{\sim }{f}}_d={\sum }_{i=1}^N\frac{x_i{\stackrel{\sim }{f}}_{d,i}}{{y_i}^2}$, where $y_i=d_i/⟨d_p⟩$.

3. ${\stackrel{\sim }{f}}_{d,i}>0$.

Note that the requirement (3-1) is required for the drag force model to reproduce the drag law for mono-size particle systems. The requirement (3-2) ensures the sum of the drag forces of all components equals the total fluid-particle drag force. The requirement (3-3) ensures positive mean drag forces.

Table 3 summarizes the distribution functions for multisize particle systems proposed by different investigators. The predictions by different correlations for comparison are given in Figure 6. These correlations were extended from monosize particle systems through distributing functions. Therefore, the performance of the drag correlations of multisize particle systems needs to be evaluated based on the requirements for both monosize and multi-size particle systems. Only the correlations of Rong et al. [13] and Mehrabadi et al. [45] can satisfy all the requirements. Others suffer various discrepancies to different extents. For example, the correlation of Sarkar et al. [43] fails to meet the criteria (3-1) and (3-2), whereas the correlations by Yin and Sundaresan [44] and Cello et al. [33] may generate negative drag forces, which is obviously against the requirement (3-3). It should also be noted that Rong et al. [13] and Mehrabadi et al. [45] used the correlations proposed in their previous studies to estimate the total fluid-particle drag force. The discussion about these two correlations for monosize systems can be found in Section 5.

References	Treatments	Drag correlation	Req. (3-1)	Req. (3-2)	Req. (3-3)
Sarkar et al. (2009) [43]	$\frac{{\stackrel{\sim }{f}}_{d,i}}{{\stackrel{\sim }{f}}_d}=y_i+0.064y_i^3$	Van der Hoef et al. (2005) [30]	$\times$	$\times$	$\surd$
Yin and Sundaresan (2009) [44]	${\stackrel{\sim }{f}}_{d,i}=\frac{1}{{\epsilon }_f}+\left({\stackrel{\sim }{f}}_d-\frac{1}{{\epsilon }_f}\right)\left(a{y}_i+\left(1-a\right)y_i^2\right)$; $a=1-2.66\left(1-{\epsilon }_f\right)+9.096{\left(1-{\epsilon }_f\right)}^2-11.338{\left(1-{\epsilon }_f\right)}^3$	Van der Hoef et al. (2005) [30]	$\surd$	$\times$	$\times$
Cello et al. (2010) [33]	$\frac{{\stackrel{\sim }{f}}_{d,i}}{{\stackrel{\sim }{f}}_d}=y_i+\frac{1-{\epsilon }_f}{{\epsilon }_f}\left(\frac{1-{\epsilon }_f-0.27}{1-0.27}\right)\frac{\left(y_i^2-y_i\right)}{{\sum }_{j=1}^N{x}_j{y}_j}$	Cello et al. (2010) [33]	$\surd$	$\times$	$\times$
Rong et al. (2014) [13]	$\frac{{\stackrel{\sim }{f}}_{d,i}}{{\stackrel{\sim }{f}}_d}=\frac{0.5{\epsilon }_f}{{\sum }_{i=1}^N\left(x_i/y_i^2\right)}+0.5\left(1-{\epsilon }_f\right)y_i^2+0.5y_i$	Rong et al. (2013) [14]	$\surd$	$\surd$	$\surd$
Mehrabadi et al. (2016) [45]	$\frac{{\stackrel{\sim }{f}}_{d,i}}{{\stackrel{\sim }{f}}_d}={\epsilon }_f{y}_i+\left(1-{\epsilon }_f\right)y_i^2$	Tenneti et al. (2011) [34]	$\surd$	$\surd$	$\surd$
Duan et al. (2020) [46]	${\stackrel{\sim }{f}}_{d,i}=\frac{1}{{\epsilon }_f}+\left({\stackrel{\sim }{f}}_d-\frac{1}{{\epsilon }_f}{\sum }_{i=1}^N\left(\frac{x_i}{y_i^2}\right)\right)\left(0.37y_i+0.63y_i^{2-\eta }\right)$; $\eta =2.169ln\left[0.466+{\left(\frac{1-{\epsilon }_f}{{\epsilon }_{\mathit{ij},\mathit{max}}}\right)}^{0.2025}\right]$; ${\epsilon }_{\mathit{ij},\mathit{max}}=\left\{\begin{array}{c}0.313\left(1-{\chi }_{\mathit{ij}}\right)x_i+0.64,x_i<0.414\\ 0.870\left(1-{\chi }_{\mathit{ij}}\right)x_j+0.64,x_i>0.414\end{array}\right.$; ${\chi }_{\mathit{ij}}=\sqrt{d_j/d_i}\left(d_j>d_i\right)$	Zhou and Fan (2015) [38]	$\surd$	$\times$	$\times$

Table 3.

Summary of the drag coefficient for multisize particle system.

Figure 6.

${\stackrel{\sim }{f}}_{d,i}/{\stackrel{\sim }{f}}_d$ as a function of porosity ${\epsilon }_f$ and volume fraction of small particle $x_s$ calculated from different drag correlations.

7. From spherical to nonspherical particle system

Particle shape could strongly affect drag forces. Rather than perfectly spherical particles, real particles often show diverse morphology. For example, they can be cubes, cylinders, and ellipsoids, or more generally, particles of irregular shapes. Such morphological diversity adds further complexity to the modeling of particle-fluid interaction forces. In earlier years, nonspherical particles are usually treated as volume-equivalent spheres. However, this simple treatment gives inaccurate results, even for an isolated particle. Several studies have been conducted to overcome this problem, particularly based on ellipsoidal particles. Such efforts are discussed in this section.

First, the definition of the geometry of an ellipsoidal particle is given to understand the factors that should be considered in drag correlations. For a standard axis-aligned ellipsoid particle, the Cartesian coordinates are given by $\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}=1$, where a, b, and c are the principal semi-axes. In most works for ellipsoidal particles, the degenerate cases are considered, where the ellipsoid particles have two equal axes, say $a=b$, generally referred to as spheroid or ellipsoid of revolution. Variation of $a$ results in different shapes of particles, which can be represented by aspect ratio, defined as $\mathit{Ar}=\frac{c}{a}=\frac{c}{b}$. Obviously, for an oblate spheroid, $\mathit{Ar}<1$; for a sphere, $\mathit{Ar}=1$; and for a prolate spheroid, $\mathit{Ar}>1$. The shapes of oblate and prolate spheroids are schematically shown in Figure 7. Note that the angle between the direction of the main flow and the symmetric axis of the ellipsoid body is defined as incident angle $\theta$.

Figure 7.

Characteristics of spheroidal particles prolate spheroid (left) and oblate spheroid (right).

Based on a large number of experimental data, Holzer and Sommerfeld [47] established a drag correlation for an isolated nonspherical particle, given as:

$$C_D=\frac{8}{\mathit{Re}}\frac{1}{\sqrt{{\psi }_{\perp }}}+\frac{16}{\mathit{Re}}\frac{1}{\sqrt{\psi }}+\frac{3}{\sqrt{\mathit{Re}}}\frac{1}{{\psi }^{3/4}}+{0.4210}^{0.4{\left(-log\psi \right)}^{0.2}}\frac{1}{{\psi }_{\perp }}$$E16

where $\psi$ is the sphericity, which is defined as the ratio between the surface area of the equivalent-volume sphere and that of the considered particle. ${\psi }_{\perp }$ is the mean crosswise sphericity, ${\psi }_{\perp }=⟨{\psi }_{\perp ,i}⟩$, when applying to a multi-particle system. Here, the individual crosswise sphericity ${\psi }_{\perp ,i}$ represents the ratio between the cross-sectional area of the equivalent-volume sphere and the projected cross-sectional area of the considered ith particle. For spheroidal particles, these two parameters can be calculated by $\psi =\frac{{\mathit{Ar}}^{2/9}}{{\left(\frac{1+2{\mathit{Ar}}^{1.61}}{3}\right)}^{1/1.61}}$ and ${\psi }_{\perp ,i}=\frac{{\mathit{Ar}}^{2/9}}{{\left({\mathit{Ar}}^2{sin}^2{\theta }_i+{\mathit{Ar}}^2{cos}^2{\theta }_i\right)}^{1/2}}$. Later, Zastawny et al. [48] and Ouchene et al. [49] established correlations based on LB and DNS simulation data.

The procedure of extension from single sphere to multi-spheres particle systems is also used in nonspherical particle systems. In this direction, various drag correlations have been established in recent years based on numerical data. Similar to the requirements as discussed for multiparticle systems, the drag correlation $C_D\left(\mathit{Re},{\epsilon }_f,\psi \right)$ for nonspherical particles should also meet the following requirements for physical consistency and generality:

1. $C_D\left(\mathit{Re},{\epsilon }_f,\psi \right)$ should be reduced to the drag correlation for spherical particles when $\psi =1$.

2. $C_D\left(\mathit{Re},{\epsilon }_f,\psi \right)$ should be reduced to the drag correlation $C_{D0}\left(\mathit{Re},\psi \right)$ for single nonspherical particle when porosity ${\epsilon }_f$ approaches unity.

3. $C_D\left(\mathit{Re},{\epsilon }_f,\psi \right)$ should be continuous over the entire application range of $\mathit{Re}$, ${\epsilon }_f$, $\psi$.

4. $C_D\left(\mathit{Re},{\epsilon }_f,\psi \right)$ should always be larger than $C_{D0}\left(\mathit{Re},\psi \right)$ when ${\epsilon }_f$ is smaller than 1;

Table 4 summarizes the drag correlations proposed by different investigators. Figure 8 shows the values of $\chi$ calculated from different drag correlations. All the correlations can meet the requirements (4-1) & (4-2) and can be reduced to the drag correlation for a single spherical particle. These correlations are established based on those for a nonspherical particle and are extensions of earlier studies from different investigators. Specifically, Rong et al. [12] used the drag correlation of Holzer and Sommerfeld [47], and their correlation is the extension of their work for monosize and multisize particle systems [13, 14]. Li et al. [50, 51] extended the work of Zhou and Fan [38]. Cao et al. [52] extended the work of Tenneti et al. [34] and adopted the drag correlation for ellipsoids proposed by Ouchene et al. [49]. Because Li et al. [50, 51] conducted their works for oblate and prolate particles separately, the drag correlation discontinues at the switching point as $\mathit{Ar}$ approaches 1. Cao et al. [52] only examined the drag correlation for prolate particles and their drag correlation may generate the negative value of the void function at relatively high $\mathit{Ar}$ or $\mathit{Re}$. So far, only the drag correlation proposed by Rong et al. [12] can satisfy all the requirements (4-1) to (4-4). However, it should be noted that all these studies only focused on particular particle shapes and a unified drag correlation with better generality for particles of irregular shapes is still lacking.

References	Drag coefficient $C_D$	Req (4-1)	Req (4-2)	Req (4-3)	Req (4-4)
Rong et al. (2015) [13]	$C_D=C_{D0}{{\epsilon }_f}^{2-\chi }$, $\chi =\beta +\lambda$, $\beta =2.65\left(1+{\epsilon }_f\right)-\left(5.3-3.5{\epsilon }_f\right){\epsilon }_f^{2}exp\left[-\frac{{\left(1.5-log\mathit{Re}\right)}^2}{2}\right]$; $\lambda =\left(1-\psi \right)\left\{101.8{\left(\psi -0.81\right)}^2+2.4-\left(39\psi -20.6\right)exp\left[-\frac{{\left(3.5-log\mathit{Re}\right)}^2}{2}\right]\right\}$, $\psi =\frac{{\mathit{Ar}}^{2/9}}{{\left(\frac{1+2{\mathit{Ar}}^{1.61}}{3}\right)}^{1/1.61}}$	$\surd$	$\surd$	$\surd$	$\surd$
Li et al. (2019, 2020) [50, 51]	$C_D=C_{D,0^{°}}{cos}^2\theta +C_{D,{90}^{°}}{sin}^2\theta$; $C_{D,0^{°}}=C_{D,\mathit{\text{sphere}}}{F}_{D0,0^{°}}{\psi }_{0^{°}}$,; $F_{D0,0^{°}}=\left\{\begin{array}{c}\frac{8}{3}{\left[-\frac{2\mathit{Ar}}{{\mathit{Ar}}^2-1}+\frac{2{\mathit{Ar}}^2-1}{{\left({\mathit{Ar}}^2-1\right)}^{3/2}}ln\left(\frac{\mathit{Ar}+\sqrt{{\mathit{Ar}}^2-1}}{\mathit{Ar}-\sqrt{{\mathit{Ar}}^2-1}}\right)\right]}^{-1},\mathit{Ar}>1\\ 1,\mathit{Ar}=1\\ \frac{8}{3{\mathit{Ar}}^{1/3}}{\left[\frac{2\mathit{Ar}}{1-{\mathit{Ar}}^2}+\frac{2\left(1-2{\mathit{Ar}}^2\right)}{{\left(1-{\mathit{Ar}}^2\right)}^{3/2}}{tan}^{-1}\left(\frac{\sqrt{1-{\mathit{Ar}}^2}}{\mathit{Ar}}\right)\right]}^{-1},\mathit{Ar}<1\end{array}\right.$; ${\psi }_{0^{°}}=\left\{\begin{array}{c}{\left(1-{\epsilon }_f\right)}^{-0.103S+0.167}\left(\mathit{Ar}-1\right)\left(0.326-0.374{\mathit{Ar}}^{-0.34}{\left(S+0.5\right)}^{1.05}\right)+1,\mathit{Ar}>1\\ 1,\mathit{Ar}=1\\ {\left(1-{\epsilon }_f\right)}^{-0.2S+0.498}{\mathit{Ar}}^{1.65}{\left(\frac{1}{\mathit{Ar}}-1\right)}^{0.706}\left(-0.550+0.302{\mathit{Ar}}^{-2.3}{\left(S+0.5\right)}^{1.08}\right)+1,\mathit{Ar}<1\end{array}\right.$; $C_{D,{90}^{°}}=C_{D,\mathit{\text{sphere}}}{F}_{D0,{90}^{°}}{\psi }_{{90}^{°}}$; $F_{D0,{90}^{°}}=\left\{\begin{array}{c}\frac{8}{3}{\left[\frac{\mathit{Ar}}{{\mathit{Ar}}^2-1}+\frac{2{\mathit{Ar}}^2-3}{{\left({\mathit{Ar}}^2-1\right)}^{3/2}}ln\left(\mathit{Ar}+\sqrt{{\mathit{Ar}}^2-1}\right)\right]}^{-1},\mathit{Ar}>1\\ 1,\mathit{Ar}=1\\ \frac{8}{3{\mathit{Ar}}^{1/3}}{\left[-\frac{\mathit{Ar}}{1-{\mathit{Ar}}^2}-\frac{2{\mathit{Ar}}^2-3}{{\left(1-{\mathit{Ar}}^2\right)}^{3/2}}{sin}^{-1}\left(\sqrt{1-{\mathit{Ar}}^2}\right)\right]}^{-1},\mathit{Ar}<1\end{array}\right.$; ${\psi }_{{90}^{°}}=\left\{\begin{array}{c}{\left(1-{\epsilon }_f\right)}^{-0.344S+0.716}{\left(\mathit{Ar}-1\right)}^{0.281}\left(0.232{\left(S+\mathit{Ar}\right)}^{0.851}-0.401\right)+1,\mathit{Ar}>1\\ 1,\mathit{Ar}=1\\ {\left(1-{\epsilon }_f\right)}^{-0.160S+0.293}{\mathit{Ar}}^{2.18}{\left(\frac{1}{\mathit{Ar}}-1\right)}^{2.65}\left(1.54-1.30{\mathit{Ar}}^{0.266}{\left(S+0.5\right)}^{1.13}\right)+1,\mathit{Ar}<1\end{array}\right.$; $S=\frac{1}{2}\left(3\overline{{cos}^2\varphi }-1\right)$	$\surd$	$\surd$	$\times$	$\surd$
Cao et al. (2020) [52]	$C_D=C_{D0}{{\epsilon }_f}^{2-\chi }$, $\chi =\beta +\lambda$, $\beta =4.988\left(0.5139+{\epsilon }_f\right)-\left(3.175-1.493{\epsilon }_f\right){\epsilon }_f^{2}exp\left[-\frac{{\left(1.5-log\mathit{Re}\right)}^2}{2}\right]-0.5884log\mathit{Re}$; $\lambda =\left(0.5201{\mathit{Ar}}^2+0.6094\mathit{Ar}-11.74\right){{\epsilon }_f}^3+\left(-1.713{\mathit{Ar}}^2-3.467\mathit{Ar}+12.65\right){\epsilon }_f+\left(1.056{\mathit{Ar}}^2+0.4316\mathit{Ar}-4.168\right)$,	$\surd$	$\surd$	$\surd$	$\times$

Table 4.

Summary of the drag correlations for nonspherical particle system.

Figure 8.

Exponent $\chi$ as a function of porosity ${\epsilon }_f$, Reynolds numbers $\mathit{Re}$ and $\mathit{Ar}$ calculated from different drag correlations.

8. From Newtonian fluid to non-Newtonian fluid system

Another important factor affecting the fluid-particle interaction force is fluid rheology. In some applications, the fluid does not necessarily follow Newton’s law of viscosity. Fluid may present shear-thinning/shear-thickening behaviors. That is, the viscosity of fluid $\mu$ increases/decreases as the shear stress $\tau$ or shear rate $\dot{\gamma }$ increases. Table 5 lists some widely used viscosity models for non-Newtonian fluid. Among these models, the power-law model is the most popular one in engineering applications because it provides a unified and simple way of describing the rheological characteristics of shear-thinning fluid (n < 1), Newtonian fluid (n = 1), and shear-thickening fluid (n > 1). However, even with the simplest power-law fluid model, predicting the behaviors of particles and fluid or quantifying the fluid-particle interaction force in non-Newtonian fluid is extremely challenging due to the variable viscosity.

Viscosity model	Correlation
Power-law fluid model	$\mu =K{\dot{\gamma }}^{n-1}$, where $K$ is the flow consistency index and n is the power index.
Spriggs fluid model	$\mu =\left\{\begin{array}{c}{\mu }_0,\dot{\gamma }\le \dot{{\gamma }_0}\\ {\mu }_0{\left(\frac{\dot{\gamma }}{\dot{{\gamma }_0}}\right)}^{n-1},\dot{\gamma }>\dot{{\gamma }_0}\end{array}\right.$, where $\dot{{\gamma }_0}$ is the reference shear rate.
Carreau fluid model	$\mu ={\mu }_{\infty }+\left({\mu }_0-{\mu }_{\infty }\right){\left[1+{\left(\lambda \dot{\gamma }\right)}^2\right]}^{\frac{n-1}{2}}$
Bingham fluid model	$\mu =\left\{\begin{array}{c}\mu +\frac{\tau }{\dot{\gamma }},{\tau }_{\mathit{max}}\ge \tau \\ \infty ,{\tau }_{\mathit{max}}<\tau \end{array}\right.$, where ${\tau }_{\mathit{max}}$ is the maximum shear stress.
Casson fluid model	$\mu =\left\{\begin{array}{c}{\left(\sqrt{\mu }+\sqrt{\frac{\tau }{\dot{\gamma }}}\right)}^2,{\tau }_{\mathit{max}}\ge \tau \\ \infty ,{\tau }_{\mathit{max}}<\tau \end{array}\right.$

Table 5.

Summary of the viscosity models for non-Newtonian fluids.

To date, attention paid to the effect of fluid rheology on drag force is much less than other factors. Table 6 summarizes the drag correlations proposed for a non-Newtonian fluid-particle system. Note that the rheology of all the considered non-Newtonian fluids obeys the power law. As Newtonian fluid is a special case when $n=1$ for the power-law viscosity model, the drag correlation for a non-Newtonian fluid-particle system should not only satisfy requirements (2-1) to (2-3) but also the following requirements:

References	Drag coefficient $C_D$	Req (2-1)	Req (2-2)	Req (2-3)	Req (5-1)
Srinivas and Chhabra (1995) [53]	$C_D=\frac{4}{3}\left(\frac{150}{{\mathit{Re}}_p}+1.75\right)$, ${\mathit{Re}}_p=\frac{{\left(15\sqrt{2}\right)}^{1-n}{\rho }_f{{\epsilon }_f}^n{\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|}^{2-n}{d_p}^n}{K{\left(1-{\epsilon }_f\right)}^n}{\left(\frac{4n}{3n+1}\right)}^n$	$\times$	$\surd$	$\times$	$\surd$
Sabiri and Comiti (2000) [54]	$C_D=4{\tau }^3\left(\frac{16}{{\mathit{Re}}_p}+0.194\right)$, ${\mathit{Re}}_p=\frac{{\rho }_f{{\epsilon }_f}^n{\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|}^{2-n}{d_p}^n{\tau }^{2-n}}{2^{2n-3}K{\left[3\left(3n+1\right)/4n\right]}^n{\left(1-{\epsilon }_f\right)}^n}$, $\tau =1-0.41ln{\epsilon }_f$ S	$\times$	$\surd$	$\times$	$\times$
Dhole et al. (2004) [55]	$C_D=\frac{4}{3}\left(\frac{230.4}{{\mathit{Re}}_k}+0.33\right)$, ${\mathit{Re}}_k=\frac{{\rho }_f{{\epsilon }_f}^{2-n}{\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|}^{2-n}{d}_p}{\phi \left(1-{\epsilon }_f\right)}$,; $\phi =\frac{K}{12}{\left(9+\frac{3}{n}\right)}^n{\left(150S{\epsilon }_f\right)}^{\left(1-n\right)/2}$, $S=\frac{{d_p}^2{{\epsilon }_f}^3}{150{\left(1-{\epsilon }_f\right)}^2}$	$\times$	$\surd$	$\surd$	$\times$
Qi et al. [57]	$C_D=C_{D0}{{\epsilon }_f}^{2-\chi }$, $C_{D0}={\left(\frac{2.654}{n^{-1.9}+3.213}+\frac{4.8}{n^{0.08}\sqrt{\mathit{Re}}}\right)}^2$, $\mathit{Re}=\frac{{\rho }_f{{\epsilon }_f}^{2-n}{\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|}^{2-n}{d_p}^n}{K}$; $\chi =2.65n{\left(\epsilon +1\right)}^n-\left(5.3-3.5{\epsilon }^n\right){\epsilon }^{2}exp\left[-\frac{{\left(log\mathit{Re}-1.5\right)}^2}{2}\right]+\frac{2.65{\left(\epsilon +1\right)}^{2-n}-2.65n{\left(\epsilon +1\right)}^n}{1+exp\left[-2\left(log\mathit{Re}-1.5\right)\right]}$	$\surd$	$\surd$	$\surd$	$\surd$

Table 6.

Summary of the drag correlations for non-Newtonian fluid-particle system.

1. $C_D\left(\mathit{Re},{\epsilon }_f,n\right)$ should be reduced to the drag correlation $C_D$ for a Newtonian fluid-particle system in Table 2 when $n=1$.

The correlation proposed by Srinivas and Chhabra [53] is established based on the Ergun-drag correlation with a different definition of Reynolds number for non-Newtonian fluids and can satisfy requirements (2-2) and (5-2). However, it also has the same problem faced in the Ergun correlation. Sabiri and Comiti [54] and Dhole et al. [55] further revised the definition of the Reynolds number and introduced the tortuosity of packed beds to consider the effect of fluid rheology. However, their correlations cannot overcome the discrepancies of the Ergun drag correlation. Besides, the correlations of Srinivas and Chhabra [53] and Sabiri and Comiti [54] may generate a negative value of void function $\chi$, as shown in Figure 9, which is against the requirement (2-3). To overcome this problem, Qi et al. [15, 56, 57] conducted a series of studies by extending the studies of Rong et al. [13, 14] from Newtonian fluid to non-Newtonian fluid. Qi et al. [56, 57] also considered both monosize and multisize particle systems. Overall, the resulting drag correlation is consistent with the drag correlation proposed for Newtonian fluid-particle systems by Rong et al. [13, 14] and satisfies all the requirements (Table 6).

Figure 9.

Exponent $\chi$ as a function of porosity ${\epsilon }_f$, Reynolds number $\mathit{Re}$ and power-law index n calculated from different drag correlations.

9. Conclusions

Drag force correlations for particle-fluid systems are reviewed, covering from simple to complicated systems, including from single particle to multiple particles, monosize to multisize, spherical to nonspherical, and Newtonian fluid to non-Newtonian fluid. The drag correlations for mono-size and multi-size spherical particle systems are more mature. However, the physical consistency and generality of different drag correlations could be more carefully considered, which are discussed in this review. Several drag correlations show superiority in generality. The understanding of fluid-particle interactions in complicated systems involving factors such as particle shape or fluid rheology is still lacking. For example, the studies on the effect of fluid rheology are largely limited to non-Newtonian fluids obeying the power-law fluid model. Similarly, the studies on nonspherical particles are mostly limited to ellipsoids. Further studies should be conducted to generally consider these important factors to meet various engineering needs.