Modeling on Effect of Particle Sediment on Fluid Flow and Heat Transfer of Solid–Fluid Suspension

Abstract

A two-way coupling particle flux model is proposed for studying the multi-component solid–fluid suspension. The suspension mixture is treated as a non-linear single-phase fluid and the migration of the solid particles is modeled by a particle flux equation. The proposed particle flux model takes the effects of the particle migration on the transport of the suspension’s momentum and internal energy into account. Two benchmark problems are calculated to study the performance of the proposed particle flux model, i.e., flow in a sudden expansion straight channel and flow between two rotating cylinders. It is found that the particle flux model converges without numerical stability issue with the commonly used PISO-SIMPLE transient solver, and the effect of the particle migration is evident on both velocity profile and temperature distribution.

1. Introduction

Solid–fluid suspensions are common in both engineering applications and geographic process. The transport of the particle sediments is fundamental for fluvial, estuarine and coastal morphology [1], sequestration processes in porous media [2], and formation of coal particle slurries with oil or water for coal feeders of power plants [3,4]. The rheological properties of these solid–fluid suspensions, such as the solid particles migration under flow, the viscosity and the heat transfer characteristics of the suspension mixture, are main issues affecting the successful operation of numerous engineering facilities [4,5,6].

For the flow of the solid–fluid suspension, when the density of the solid particles is remarkably larger than the density of the carrying fluid, particles sedimentation under gravity occurs. Such process could be studied numerically by the Computation Fluid Dynamics-Discrete Element Method (CFD-DEM) [7,8] and the two-fluid method [3,9,10], at the expense of added computational cost from the introduction of Lagrangian discrete phase or another fluid phase. A lower-cost variant takes the advantage of particle flux equation to model the effect of solid particles, so that both phases of the suspension could be calculated with single-phase transport equations. Compared with CFD-DEM and two-fluid method, the particle flux model typically has higher numerical stability, since only a scalar transport equation is added into the original Naiver–Stokes equations of the single-phase fluid. The advantage on numerical stability could be important for problems with complex geometry or extreme flow condition.

The most widely used particle flux method was proposed by Phillips et al. [11,12,13,14]. However, as a one-way coupling model, it considers only the effect of fluid on the motion of particles, the effect of the migration/transport of the solid particles back on the conservation of the linear momentum and internal energy are neglected. To address such effect, we derived a new particle flux model based on Phillips’ model [14], considering the effect of the particles migration/transport on the conservation of the linear momentum and internal energy of the fluid phase. The particle migration is completed using Phillips’ model with Tetlow’s empirical correlation [15]. Thermal diffusivity is formulated following Jeffery’s model, which considers the second order effects of thermal conductivity based on the solid volume fraction [4,16].

As the solid–fluid suspension typically shows complex non-linear rheological properties, it is important to properly model the constitutive relations within the stress tensor and heat flux. Although in some cases the solid–fluid suspension could be approximated by the Bingham or the Power law [17], under more common conditions the stress tensor strongly depends on the volume fraction of the solid component and the shear rate of the fluid flow. Krieger et al. [18] provided a general expression for the dependency of the viscosity on both solid volume fraction and shear rate. In addition, in applications such as the high speed bearing flow and deep sea drilling, the effect of by pressure and temperature on viscosity [4,5] must also be considered.

The paper is organized as follows: In the Methods section, the effect of the particle migration on the transport of the linear momentum and internal energy are investigated, based on Phillips’s model [14]. Then we introduce the constitutive equations for shear stress, heat flux and particles transport flux implemented in the proposed two-way coupling particle flux model. In the Results section, two benchmark problems are calculated by the proposed particle flux model and the results are discussed.

2. Methods

2.1. Governing Equations

For the conservation of mass, for both solid and fluid components based on the mixture theory

$$\frac{\partial {\rho }_f}{\partial t}+div\left({\rho }_f\mathit{v}\right)=div\left(-\mathit{N}{\rho }_{f0}\right)$$

(1)

$$\frac{\partial {\rho }_s}{\partial t}+div\left({\rho }_s\mathit{v}\right)=div\left(\mathit{N}{\rho }_{s0}\right)$$

(2)

where ${\rho }_f=\left(1-\varphi \right){\rho }_{f0}=\alpha {\rho }_{f0}$, ${\rho }_s=\varphi {\rho }_{s0}$, ${\rho }_{f0}$ and ${\rho }_{s0}$ are the pure density of the fluid and the solid particle, in the reference configuration, therefore the mixture density is defined as${\rho }_m=\alpha {\rho }_{f0}+\varphi {\rho }_{s0}$; $\varphi$ is the volume fraction of the solid particles, where $0\le \varphi <{\varphi }_{max}<1$. $\mathit{v}$ is the velocity of the suspension mixture. $\mathit{N}$ is the solid particle flux, which measures the transport/migration of the volume fraction of the solid component (for more details please refer to Equations (9) and (18)). Thus, $\mathit{N}{\rho }_{s0}$ and $-\mathit{N}{\rho }_{f0}$ represent the mass transport of each component caused by the particle flux. Combining Equations (1) and (2), we can acquire

$$\frac{\partial \left(\frac{{\rho }_f}{{\rho }_{f0}}\right)}{\partial t}+div\left(\frac{{\rho }_f}{{\rho }_{f0}}\mathit{v}\right)+\frac{\partial \left(\frac{{\rho }_s}{{\rho }_{s0}}\right)}{\partial t}+div\left(\frac{{\rho }_s}{{\rho }_{s0}}\mathit{v}\right)=div\left(-\mathit{N}\right)+div\left(\mathit{N}\right)$$

(3)

$$div \mathit{v}=0$$

(4)

While the conservation of linear momentum can be written as

$$\begin{array}{l}\frac{\partial \left(\alpha {\rho }_{f0}+\varphi {\rho }_{s0}\right)\mathit{v}}{\partial t}+div\left(\left(\alpha {\rho }_{f0}+\varphi {\rho }_{s0}\right)\mathit{v}\mathit{v}\right)\\ =div\left(\mathit{T}\right)+div\left(\mathit{N}\left({\rho }_{s0}-{\rho }_{f0}\right)\mathit{v}\right)+\left(\alpha {\rho }_{f0}+\varphi {\rho }_{s0}\right)\mathit{b}\end{array}$$

(5)

substituting Equations (1) and (2) into (5), we can acquire

$$\begin{array}{c}\left(\alpha {\rho }_{f0}+\varphi {\rho }_{s0}\right)\left[\frac{\partial \mathit{v}}{\partial t}+\nabla \left(\mathit{v}\right)\mathit{v}\right]\\ =div\left(\mathit{T}\right)-div\left(\mathit{N}\left({\rho }_{s0}-{\rho }_{f0}\right)\right)\mathit{v}+div\left(\mathit{N}\left({\rho }_{s0}-{\rho }_{f0}\right)\mathit{v}\right)+\left(\alpha {\rho }_{f0}+\varphi {\rho }_{s0}\right)\mathit{b}\end{array}$$

(6)

In the current paper, we consider gravity $\mathit{g}$ as the only body force. Term $div\left(\mathit{N}\left({\rho }_{s0}-{\rho }_{f0}\right)\mathit{v}\right)-div\left(\mathit{N}\left({\rho }_{s0}-{\rho }_{f0}\right)\right)\mathit{v}$ is the linear momentum change caused by the particle flux $\mathit{N}$.

As for the conservation of energy, following the previous derivation we have

$$\begin{array}{l}\frac{\partial \left(\alpha {\rho }_{f0}C{p}_{f0}+\varphi {\rho }_{s0}C{p}_{s0}\right)\theta }{\partial t}+div\left(\left(\alpha {\rho }_{f0}C{p}_{f0}+\varphi {\rho }_{s0}C{p}_{s0}\right)\theta \mathit{v}\right)\\ \hspace{1em}=\mathit{T}:\mathit{L}+div\left(q\right)+div\left(\mathit{N}\left({\rho }_{s0}C{p}_{s0}-{\rho }_{f0}C{p}_{f0}\right)\theta \right)\end{array}$$

(7)

$$\begin{array}{l}\left(\alpha {\rho }_{f0}C{p}_{f0}+\varphi {\rho }_{s0}C{p}_{s0}\right)\left[\frac{\partial \theta }{\partial t}+\nabla \left(\theta \right)·\mathit{v}\right]\\ \hspace{1em}=\mathit{T}:\mathit{L}+div\left(\mathit{q}\right)+div\left(\mathit{N}\left({\rho }_{s0}C{p}_{s0}-{\rho }_{f0}C{p}_{f0}\right)\theta \right)-div\left(\mathit{N}\left({\rho }_{s0}C{p}_{s0}-{\rho }_{f0}C{p}_{f0}\right)\right)\theta \end{array}$$

(8)

where L is the gradient of velocity, q is the heat flux vector, $\alpha {\rho }_{f0}C{p}_{f0}+\varphi {\rho }_{s0}C{p}_{s0}={\rho }_{m}C{p}_m$ is the mixture heat capacity and $div\left(\mathit{N}\left({\rho }_{s0}C{p}_{s0}-{\rho }_{f0}C{p}_{f0}\right)\theta \right)-div\left(\mathit{N}\left({\rho }_{s0}C{p}_{s0}-{\rho }_{f0}C{p}_{f0}\right)\right)\theta$ is the internal energy change due to the transport of the particle flux $\mathit{N}$.

The conservation of solid concentration is simply

$$\frac{\partial \varphi }{\partial t}+\mathit{v}\frac{\partial \varphi }{\partial x}=-div\mathit{N}$$

(9)

where $\mathit{N}$ is the particle diffusive flux.

2.2. Constitutive Equations

Roscoe [19] provided a thorough review of modeling and characterizing gel-like materials. In general, gels behave as non-Newtonian fluids (NNF). (Non-linear) NNF exhibit some unusual characteristics. For example, some NNF show normal stress effects, as manifested in phenomena such as die-swell or rod-climbing (see the book by Schowalter [20]). Some NNF have yield stress or show viscoelastic effects, while some have viscosities depending on shear rate and temperature or even pressure (see the books by Macosko [21], Carreau, De Kee and Chhabra [22]). Perhaps the simplest model, from an engineering perspective, which can capture the shear-thinning (or shear-thickening) effect, is the power-law model. In current work, for the shear stress tensor, we assume that the viscous stress of the solid–fluid suspension can be represented by a generalized power-law fluid model,

$${\mathit{T}}_v=-p\mathit{I}+\left({\mu }_t+{\mu }_{fr}{\left(1-\frac{\varphi }{{\varphi }_{max}}\right)}^{-\beta }{\Pi }^{\frac{m}{2}}\right)\mathit{D}$$

(10)

$$\Pi =2tr{\mathit{D}}^2, \mathit{D}=\frac{1}{2}\left(\mathit{L}+{\mathit{L}}^T\right), \mathit{L}=gradv$$

(11)

where $\varphi$ is the volume fraction, ${\mu }_t$ is the turbulent eddy viscosity, ${\mu }_{fr}$ is the reference suspending fluid viscosity, and $m$ is the shear dependency constant. In the above equations, ${\left(1-\frac{\varphi }{{\varphi }_m}\right)}^{-\beta }$ depends on Krieger’s work which provides a correlation between the viscosity and particle concentration, formulated by both theoretical analysis and experimental data ($\mu ={\mu }_{fr}{\left(1-\varphi /0.68\right)}^{-1.82}$) [18]. ${\Pi }^{\frac{m}{2}}$ represents the effect of volume fraction and shear-dependency. In current work, we assume that ${\mu }_{fr}$ is a constant and predefined by the user.

The heat flux vector is modeled as follows:

$$\mathit{q}=-k_m\nabla \theta$$

(12)

where$k_m$ is the effective thermal conductivity. Various models are available to evaluate the effective thermal conductivity for solid–fluid suspensions [23,24,25], in current study we apply the model proposed by Jeffery [16], keeping consistent with our previous study [4]. Jeffery’s model considers the second order effects of thermal conductivity based on the solid volume fraction [26]:

$$k_m={\kappa }_f\left[1+3\xi \varphi +\widehat{\xi }{\varphi }^2\right]+O\left({\varphi }^3\right)$$

(13)

where

$$\widehat{\xi }=3{\xi }^2+\frac{3{\xi }^3}{4}+\frac{9{\xi }^3}{16}\left(\frac{\omega +2}{2\omega +3}\right)+\frac{3{\xi }^4}{26}+\dots$$

(14)

and

$$\xi =\frac{\omega -1}{\omega +2}$$

(15)

For the constitution of particle flux, we do not take the Brownian motion into account since the size of the particles of interest is sufficiently large. Then, the diffusive particle flux becomes

$$\mathit{N}={\mathit{N}}_c+{\mathit{N}}_{\mu }+{\mathit{N}}_g+{\mathit{N}}_t$$

(16)

where ${\mathit{N}}_c$, ${\mathit{N}}_{\mu }$, ${\mathit{N}}_g$ and ${\mathit{N}}_t$ are the flux due to particles collision, spatial variations in the viscosity, gravity and turbulence, respectively.

Following Phillips’ model [14], particle flux due to collisions ${\mathit{N}}_c$ and non-uniform viscosity ${\mathit{N}}_{\mu }$ are evaluated as

$${\mathit{N}}_c=-a^2\varphi K_c\nabla \left(\dot{\gamma }\varphi \right)$$

(17)

$${\mathit{N}}_{\mu }=-a^2{\varphi }^2\dot{\gamma }{K}_{\mu }\nabla \left(ln \mu \right)$$

(18)

$$\dot{\gamma }={\left(2{\mathit{D}}_{ij}{\mathit{D}}_{ij}\right)}^{1/2}={\left(\Pi \right)}^2$$

(19)

where a is the characteristic particle length (e.g., radius), $\dot{\gamma }$ is the local shear rate, $\mu$ is the effective viscosity, and $K_c,$ $K_{\mu }$ are empirically determined coefficients. According to Phillips et al. [14], ${\mathit{N}}_c$ accounts for the effect of the spatially varying interaction frequency due to the collision between particles, while ${\mathit{N}}_{\mu }$ accounts for the effect of the spatially varying viscosity. Mathematically ${\mathit{N}}_c$ and ${\mathit{N}}_{\mu }$ are responsible for the particles migrations causing non-homogenous distribution of the particles. Seifu et al. [27] noted that the particle flux could be given as $\mathit{N}=-a^2{K}_c{\varphi }^2\dot{\gamma }\nabla \left(ln\left(\dot{\gamma }\varphi {\mu }^{K_{\mu }/K_c}\right)\right)$. We could see that $ln\left(\dot{\gamma }\varphi {\mu }^{K_{\mu }/K_c}\right)$ is a field potential incorporating mechanisms that cause migration of particles. Phillips et al. and Subia et al. [6,14] assumed that $K_c,$ $K_{\mu }$, are constants. In this paper, we apply the empirical correlation proposed by Tetlow et al. [15]

$$\frac{K_c}{K_{\mu }}=1.25\varphi +0.075; K_{\mu }=0.62$$

(20)

As for the gravity-induced particle flux, based on the work of Acrivos et al. [28] and Mavromoustaki [11], ${\mathit{N}}_g$ could be modeled as

$${\mathit{N}}_g=\frac{2}{9}\varphi f\left(\varphi \right)\frac{a^2\left({\rho }_{s0}-{\rho }_{f0}\right)}{{\mu }_{fr}}\mathit{g}$$

(21)

where

$$f\left(\varphi \right)=\frac{\left(1-\varphi \right){\mu }_{fr}}{\mu \left(\varphi \right)}$$

(22)

This gravity flux model has been applied to various studies, such as falling of thin film [11]. For other mathematical models of gravity induced particle flux, please refer to [1,29].

The particle flux due to turbulence reads

$${\mathit{N}}_t=\frac{{\nu }_t}{Sc}\nabla \left(\varphi \right)$$

(23)

where ${\nu }_t$ is the eddy viscosity and $Sc$ is the Schmidt number which is assumed as 0.9 for current study [30,31].

By substituting the Equations (10) and (11) into (6), (12)–(15) into (8), and (16)–(23) into (9), we obtain a set of expanded governing equations for the particle flux model. For the expanded form of governing equations based on the original Phillips’s model [7], please refer to Appendix A.

Conservation of mass:

$$div\left(\mathit{v}\right)=0$$

(24)

Conservation of linear momentum:

$$\begin{array}{l}\left(\alpha {\rho }_{f0}+\varphi {\rho }_{s0}\right)\left[\frac{\partial \mathit{v}}{\partial t}+\nabla \left(\mathit{v}\right)\mathit{v}\right]\\ \hspace{1em}=-grad p+div \left(\left({\mu }_t+{\mu }_{fr}{\left(1-\frac{\varphi }{{\varphi }_{max}}\right)}^{-\beta }{\left(2tr{\mathit{D}}^2\right)}^{\frac{m}{2}}\right)\mathit{D}\right)\\ \hspace{1em}-div\left(\mathit{N}\left({\rho }_{s0}-{\rho }_{f0}\right)\right)\mathit{v}+div\left(\mathit{N}\left({\rho }_{s0}-{\rho }_{f0}\right)\mathit{v}\right)+\left(\alpha {\rho }_{f0}+\varphi {\rho }_{s0}\right)\mathit{g}\end{array}$$

(25)

Conservation of energy:

$$\begin{array}{l}\left(\alpha {\rho }_{f0}C{p}_{f0}+\varphi {\rho }_{s0}C{p}_{s0}\right)\left[\frac{\partial \theta }{\partial t}+\nabla \left(\theta \right)·\mathit{v}\right]\\ \hspace{1em}=\left(\left({\mu }_t+{\mu }_{fr}{\left(1-\frac{\varphi }{{\varphi }_{max}}\right)}^{-\beta }{\left(2tr{\mathit{D}}^2\right)}^{\frac{m}{2}}\right)\mathit{D}\right):\mathit{L}\\ \hspace{1em}+div \left(\left(\frac{{\nu }_t{\rho }_f}{Pr}+{\kappa }_f\left(1+3\xi \varphi +\widehat{\xi }{\varphi }^2\right)\right)\nabla \theta \right)\\ \hspace{1em}+div\left(\mathit{N}\left({\rho }_{s0}C{p}_{s0}-{\rho }_{f0}C{p}_{f0}\right)\theta \right)-div\left(\mathit{N}\left({\rho }_{s0}C{p}_{s0}-{\rho }_{f0}C{p}_{f0}\right)\right)\theta \end{array}$$

(26)

Conservation of particles concentration:

$$\begin{array}{l}\frac{\partial \varphi }{\partial t}+\mathit{v}\frac{\partial \varphi }{\partial x}=div\left(a^2\varphi K_c\dot{\gamma }\nabla \left(\varphi \right)\right)+div\left(a^2\varphi K_c\nabla \left(\dot{\gamma }\right)\varphi +a^2\dot{\gamma }{K}_{\mu }\nabla \left(\mathit{ln}\mu \right){\varphi }^2\right)\\ \hspace{1em}\hspace{1em}+div\left(\frac{2}{9}\varphi f\left(\varphi \right)\frac{a^2\left({\rho }_{s0}-{\rho }_{f0}\right)}{{\mu }_{fr}}\mathit{g}\right)+div\left(\frac{{\nu }_t}{Sc}\nabla \left(\varphi \right)\right)\end{array}$$

(27)

3. Results and Discussion

To evaluate the performance of the two-way coupling particle flux model, two benchmark problems were calculated and analyzed in this chapter—namely, the flow in straight sudden expansion channel and the flow between two rotating cylinders (circles). The simulation results would also be compared with the original Phillips’ model [14]. Numerical simulation of both benchmark problems was finished by a user-defined solver developed based on the open sourced CFD toolbox OpenFOAM [16]. The NS equation was solved by the standard PISO-SIMPLE (so-called PIMPLE) algorism, a transient segregated solver iterating between pressure and velocity to enforce mass conservation over momentum equations. The pressure-velocity iteration would be renewed to the next time step should initial (normalized) residual of pressure was lower than $1\times {10}^{-5}$ and that of the velocity was lower than $1\times {10}^{-6}$. During the transient calculations, the maximum Courant (maxCo) number was kept below 0.9 everywhere in the flow filed, to resolve the time-dependent information [32].

For spatial discretization, central difference (CD) interpolation was used for gradient and Laplacian terms, while convection terms were interpolated by a so-called “bounded Gauss linear” scheme, which improve the boundedness of the simulation by taking advantage of the incompressible part of velocity divergence term (non-zero before convergence was reached). For temporal discretization, implicit Euler scheme was applied for time derivative terms. As both cases were wall bounded flow, Wilcox’s k-omega model was used for the treatment of turbulence effect, and coefficients were chosen according to [33]. The underlying assumption of unsteady Reynolds-averaged Navier–Stokes (URANS) turbulence method was that the local turbulence scale estimated by RANS is sufficiently smaller than the dominating large-scale fluctuation [34,35,36,37,38,39]. Standard wall functions were applied for turbulence variables such as k (turbulent kinetic energy) and omega (specific dissipation of k) on all solid walls, where no-slip and zero gradient boundary conditions were used for velocity and pressure, respectively. The value of the gravity $\mathit{g}$, was set as a constant value of $\left(0,-10,0\right)\mathrm{m}·{\mathrm{s}}^{-2}$. To reduce computational cost of transient simulation, both calculations were performed on two-dimensional (2D) grid(s). It was acknowledged that the three-dimensional dynamics of vortices stretching could not be captured by 2D grid with URANS turbulence treatment. The purpose of the benchmark calculation served essentially as a test on the model performance and its numerical stability, by comparison with the original Phillips’s model. When the local small transient scale was not dominant or not the major interest of the study, the computationally cheaper 2D calculation could also provide valuable information, especially for the testing of the new model and numerical method [40,41].

3.1. Flow in Straight Sudden Expansion Channel

Figure 1 and panel (a) of Figure 2 showed the computational grid and geometry of the straight sudden expansion channel. The inlet velocity was 0.5 $\mathrm{m}·{\mathrm{s}}^{-1}$, and the inlet volume fraction of the solid component was 0.05. The radius of particle was 1 $\mathrm{mm}$. Temperature of the inlet flow was set as 393 $\mathrm{K}$, as the suspending fluid was water, laminar viscosity was set as a constant of 1 $cP$ in accordance with the temperature. The density of the suspending fluid and the solid particles were 1000 and 3000 $\mathrm{kg}·{\mathrm{m}}^{-3}$. As a result, if we used half the step height at the sudden expansion as the reference length scale, corresponding hydraulic Reynolds would be $1\times {10}^5$, verifying the application of turbulence model. The mass heat capacity of the water was set as 4023 $\mathrm{J}/\left(\mathrm{kg}·\mathrm{K}\right)$. To enhance the effect of the particle flux on the fluid temperature, we assumed a high heat capacity for the particles, which was equal to that of the carrying water. On all solid walls, temperature was kept constant at 293 $\mathrm{K}$. For the initial condition of the solid phase, velocity, temperature and solid volume fraction were 0 $\mathrm{m}·{\mathrm{s}}^{-1}$, 293 $\mathrm{K}$ and 0, respectively. It should be pointed out that to demonstrate the path of the particles sediment, the gravity body force would not take effect until the sudden expansion, in other words, there was no gravity in the inlet channel before the sudden expansion.

Governing equations of the suspension was discretized by a 2D structured grid, a constant spatial resolution of 1 $\times$ 1 $\mathrm{mm}$ was applied throughout the entire flow domain, as shown in Figure 1, and from then on it would be referred to as the “regular” grid or simply the “grid used”. To confirm grid independence of solution, simulation was conducted with different levels of spatial resolutions. Two of them could be found in Figure 1, namely a “coarse” and “fine” grid, as the corresponding step size was 2 and 0.5 $\mathrm{mm}$, respectively. Calculations on all three grids converged reasonably without numerical stability issues. Additionally, the velocity in the X direction (${\mathit{v}}_x$) was sampled at the geometric center of the expanded channel, namely at the location (x = 0.125$\mathrm{mm}$, y = 0). As shown in panel (d) of Figure 1. In fact, for the “fine” and “regular” grid, both velocity profile and solid volume fraction distribution converged into quasi steady state after flow time (t) = around 20 s, due to the strong diffusion effect of RANS turbulent model. However, velocity profile of the “coarse” grid calculation kept oscillating, even the ensemble mean value was quite different from the predictions of finer grid, as shown in Figure 1d. As refining the “regular” grid brought little difference on the simulation results, we assumed that the evenly distributed 1 mm step size was sufficient for current study, at least in terms of grid independency. The discussions of analysis of the simulation results were based on the “regular” grid. For a more detailed verification study of the numerical setup, please check Appendix B.

Simulation was performed for 60 s of flow time for all the cases with different parameters, for both the new particle flux model and the original Phillips model, total wall-clock times were around $1\times {10}^4$ s (~2.8 h) on a dual-channel machine with DDR4 DRAM clocked at 3200 MHz. The reason was that in typical CFD calculations, bottleneck of computation workload was the iteration between p-U coupling, especially the pressure equation, so the extra equations introduced by the new particle flux model brought little CPU overhead.

Panel (b), (c) and (d) of Figure 2 showed the velocity, solid volume fraction and temperature fields after t = 2 s. Inside each panel, the top contour showed the simulation results using the original Phillips’ model [14] but without considering the gravity; the middle contour represents results by Phillips’ model [14] considering the gravity effect; while simulation results of the proposed particle flux model considering the gravity effect were displayed in the bottom contour. According to the top contours of all the three panels, we could see that when Phillips’s model [14] was applied and gravity was neglected, all fields were symmetric along the channel symmetric line. Whereas when gravity body force was considered by Phillips’ model [14], in other words, when there is gravity but the effect of the particle flux on the fluid was ignored, particles sediment near the bottom wall was obvious after the sudden expansion, which also slightly modified the velocity and temperature profiles. However, as shown in the bottom contours of panel (b), (c) and (d), when the effect of the particle flux on the fluid was also considered (namely applying the proposed particle flux model), the high velocity and high temperature region shifted substantially towards the particles sediment direction. One explanation could be that the settling particle carried the momentum and the temperature from center region towards the bottom wall, so that the particle flux enhanced the convection of momentum and heat transfer. More interestingly, if one looked at the volume fraction field shown in panel (c), it could be found that the proposed two-way coupling particle flux model predicted a streamlined particles sediment shape, the solid component moves more like a projectile motion. Such streamlined shape of the particles sediment should be attributed to the different velocity field predicted by the proposed model, which estimated higher convection rate of momentum transfer from the high-speed center region to the bottom part of the channel.

Figure 3, Figure 4 and Figure 5 showed the evolution of the velocity, solid volume fraction and temperature fields, predicted by the proposed two-way coupling particle flux model with gravity considered. From the first 1 to 3 s, the high-velocity region gradually moved towards the direction of the gravity. At 20 s, the velocity field was almost fully developed. Along the flow direction, the high velocity region firstly leaned towards the bottom wall after the expansion, and then climbed towards the opposite direction, induced by the evolution of volume fraction distribution. It could be found that almost all the solid particles settled down rapidly near the bottom wall, as shown in Figure 2c and Figure 4c. This should be related to the large particle size and density difference between the particles and suspending fluid. Furthermore, Figure 4 showed that as the flow evolved, the particles sediment region kept growing along the flow direction. Within the first 3 s, the sediment region maintained a streamlined shape, while as the front of the sediment was pushed further away from the inlet, the tailing slope of the interface between the sediment and suspending fluid becomes steeper. This could be explained as that the settled particles far away from the sudden expansion were transported from the surface of the upstream sediment rather than settling down from the main stream of the flow (by gravity effect). The temperature profiles of the suspension, as displayed in Figure 5, evolved towards a more uniform distribution during the process of sedimentation. The temperature distribution in the bottom corner, in particular, is even more uniform compared with the rest of flow domain, due to the effect of the particle flux.

3.2. Flow between Two Rotating Cylinders (Circles)

In this section, the performance of the proposed particle flux model was evaluated by another benchmark problem: flow between two rotating cylinders. This problem was quite common in engineering applications, such as the bearing flow and the well drilling flow in the oil and geothermal industry. Figure 6 showed the geometry of the two cylinders. The radius of the inner and the outer cylinders were 0.64 and 2.54 $\mathrm{cm}$, respectively. The inner cylinder was rotating at 180 RPM and the outer cylinder is static. The kinetic viscosity of the suspending fluid was set as $1\times {10}^5$, leading to a hydraulic Reynolds number around $1.5\times {10}^5$, so the entire flow region could be regarded as fully developed turbulence [42,43,44]. The temperature of the outer cylinder and the inner cylinder were 393 K and 293 K, respectively. To speed up the simulation, firstly a steady state calculation was carried out without considering the motion of either fluid or solid, then the temperature profile of the steady state solution was used as the initial condition of transient calculation. The initial velocity of transient calculation is 0 $\mathrm{m}·{\mathrm{s}}^{-1}$, and the initial solid volume fraction is 0.1. The gravity was along negative Z-direction. The other parameters and flow conditions were kept unchanged with the previous case.

A fully structured O-grid was used for calculation. As no special refinement was applied to the boundary layer, bunching of grid layers in both radial and circumferential direction was equally divided. To contain computational cost, grid dependence study similar to the previous section was performed, and spatial resolution was chosen so that further refinement would not improve the simulation results. As the results, typical grid step size was set to be 1 $\mathrm{mm}$ in the near wall region of outer cylinder, hence near the inner cylinder wall it was around 0.4 $\mathrm{mm}$, the total cell count was around 12,000. The transient simulation converged to a quasi-steady state after around 1 s of flow time (t = 1 s), and the calculation was continued until t = 12 s. For both proposed particle flux model and original Phillips model, total wall-clock times were around 4000 s on the same dual-channel machine used for simulation of sudden expansion flow, due to the same reason that the extra equations introduced by the new particle flux model were not the bottleneck of computation workload.

Figure 7 showed the velocity, solid volume fraction and temperature distribution after t = 0.2 s, the top panels represented the results predicted by the Phillips’ model [14], while the bottom panels represented the results by the proposed two-way coupling model. In fact, for the solid volume fraction distribution, the results by two different models were quite similar, although the effect of the particle flux could still be observed for the temperature and velocity fields, particularly around the inner cylinder.

Figure 8, Figure 9 and Figure 10 show the evolution of velocity, solid volume fraction and temperature fields at different time steps. From Figure 8, we can see that as flow time increases, the solid particles migrated towards the bottom wall of the outer cylinder. Due to the effect of the rotation of the inner cylinder, the particles settling on the top wall of the inner cylinder was not symmetric along the y direction and those particles moved gradually towards the bottom of the outer cylinder. We might infer that if the inner cylinder was static, the time required to complete the particle sediment towards the bottom of the outer cylinder would be much longer. Figure 9 shows an obvious effect of the particle sediment on the velocity field. Near the inner cylinder, the high velocity region always corresponded to the high concentration of the solid particles where the viscosity also took high values. Figure 10 shows the temperature evolution with the increasing of flow time. Near the top wall of the outer cylinder, the high temperature region was slightly shifted towards the direction of gravity, as shown in Figure 10 for t = 0 s and t = 1 s. This could be attributed to the enhancement of heat transfer termed from the particle flux, which was strongly controlled by the force of gravity due to density difference between solid and fluid.

4. Conclusions

The particle flux model by Phillips et al. [14] was further developed so that a two-way coupling particle flux equations—namely, the particle flux model could be formulated by considering the effect of the particle migration on the transport of the momentum and internal energy of the fluid phase. The proposed two-way coupling particle flux model manifested reasonable numerical stability during the calculation of two benchmark problems, whose results indicated that the effect of the particle migration on the flow and heat transfer was evident; hence, in a situation where the density of a particle is sufficiently higher than the suspending fluid, it is expected that the proposed model could provide more reasonable prediction than one-way coupling models. Further extension of the application to three-dimensional process with more complex geometry is expected though, as the proposed particle flux model brought almost no additional computational cost. Besides, although solution verification showed that the numerical methods behaved as expected in terms of estimated error and observed order, carefully designed experimental validation for more flow conditions should be carried out for further calibration and development of the proposed particle flux model.