2.1 Fictitious domain method

A parallel direct-forcing fictitious domain method (DF/FD) is employed for the interface-resolved simulations of particle-laden turbulent flows (Yu et al. Reference Yu, Lin, Shao and Wang2016b). The fictitious domain (FD) method for particulate flows was originally proposed by Glowinski et al. (Reference Glowinski, Pan, Hesla and Joseph1999). The key idea of this method is that the interior of the particles is filled with the fluids and the inner fictitious fluids are enforced to satisfy the rigid body motion constraint through a pseudo-body force, which is introduced as a distributed Lagrange multiplier in the FD formulation (Glowinski et al. Reference Glowinski, Pan, Hesla and Joseph1999). In the following, we describe the DF/FD method briefly, and the reader is referred to Yu & Shao (Reference Yu and Shao2007) for the details.

For simplicity of the description, only one particle is considered in the following exposition. Let $P(t)$ represent the solid domain and $\unicode[STIX]{x1D6FA}$ represent the entire domain, including the interior and exterior of the solid body. The bulk velocity of the suspension (i.e. the fluid–solid mixture) $u_{b}$ is kept constant in our simulations. We take the channel half-width $H$ as the characteristic length and the bulk velocity $u_{b}$ as the characteristic velocity for the non-dimensionalization scheme. The bulk Reynolds number $Re_{b}$ is defined by $Re_{b}=u_{b}H/\unicode[STIX]{x1D708}$, where $\unicode[STIX]{x1D708}$ is the fluid kinematic viscosity. The dimensionless FD formulation for the incompressible fluid can be written as follows:

(2.1)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=\frac{\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}}{Re_{b}}-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D740}\quad \text{in }\unicode[STIX]{x1D6FA}, & \displaystyle\end{eqnarray}$$

(2.2)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}=\boldsymbol{U}+\unicode[STIX]{x1D74E}_{p}\times \boldsymbol{r}\quad \text{in }P(t), & \displaystyle\end{eqnarray}$$

(2.3)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0\quad \text{in }\unicode[STIX]{x1D6FA}, & \displaystyle\end{eqnarray}$$

(2.4)$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D70C}_{r}-1)V_{p}^{\ast }\left(\frac{\text{d}\boldsymbol{U}}{\text{d}t}-Fr\frac{\boldsymbol{g}}{g}\right)=-\int _{P}\unicode[STIX]{x1D740}\,\text{d}\boldsymbol{x}, & \displaystyle\end{eqnarray}$$

(2.5)$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D70C}_{r}-1)\frac{\text{d}(\boldsymbol{J}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D74E}_{p})}{\text{d}t}=-\int _{P}\boldsymbol{r}\times \unicode[STIX]{x1D740}\,\text{d}\boldsymbol{x}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$, $p$, $\unicode[STIX]{x1D740}$, $\boldsymbol{U}$, $\unicode[STIX]{x1D74E}_{p}$, $\boldsymbol{r}$ represent the fluid velocity, fluid pressure, distributed Lagrange multiplier, particle translational velocity, particle rotational velocity and position vector with respect to the mass centre of the particle, respectively; $\unicode[STIX]{x1D70C}_{r}$ is the ratio of the particle density ratio $\unicode[STIX]{x1D70C}_{s}$ to the fluid density ratio $\unicode[STIX]{x1D70C}_{f}$, i.e. $\unicode[STIX]{x1D70C}_{r}=\unicode[STIX]{x1D70C}_{s}/\unicode[STIX]{x1D70C}_{f}$; $Fr$ represents the Froude number defined by $Fr=gH/u_{b}^{2}$, with $g$ being the gravitational acceleration; $V_{p}^{\ast }$ and $\boldsymbol{J}^{\ast }$ are the dimensionless particle volume and moment of inertia tensor, defined by $V_{p}^{\ast }=V_{p}/H^{3}$ and $\boldsymbol{J}^{\ast }=\boldsymbol{J}/\unicode[STIX]{x1D70C}_{s}H^{5}$, respectively.

A fractional-step time scheme is used to decouple the system (2.1)–(2.5) into the following two sub-problems.

Fluid sub-problem for $\boldsymbol{u}^{\ast }$ and $p$

(2.6)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\boldsymbol{u}^{\ast }-\boldsymbol{u}^{n}}{\unicode[STIX]{x0394}t}-\frac{1}{2}\frac{\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}^{\ast }}{Re_{b}}=-\unicode[STIX]{x1D735}p-\frac{1}{2}[3(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u})^{n}-(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u})^{n-1}]+\frac{1}{2}\frac{\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}^{n}}{Re_{b}}+\unicode[STIX]{x1D740}^{n}, & \displaystyle\end{eqnarray}$$

(2.7)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }=0. & \displaystyle\end{eqnarray}$$

A finite-difference-based projection method on a homogeneous half-staggered grid is used for the solution of the above fluid sub-problem. All spatial derivatives are discretized with the second-order central difference scheme.

Particle sub-problem for $\boldsymbol{U}^{n+1}$, $\unicode[STIX]{x1D74E}_{p}^{n+1}$, $\unicode[STIX]{x1D740}^{n+1}$ and $\boldsymbol{u}^{n+1}$

(2.8)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{r}V_{p}^{\ast }\frac{\boldsymbol{U}^{n+1}}{\unicode[STIX]{x0394}t}=(\unicode[STIX]{x1D70C}_{r}-1)V_{p}^{\ast }\left(\frac{\boldsymbol{U}^{n}}{\unicode[STIX]{x0394}t}+Fr\frac{\boldsymbol{g}}{g}\right)+\int _{P}\left(\frac{\boldsymbol{u}^{\ast }}{\unicode[STIX]{x0394}t}-\unicode[STIX]{x1D740}^{n}\right)\,\text{d}\boldsymbol{x}, & \displaystyle\end{eqnarray}$$

(2.9)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{r}\frac{\boldsymbol{J}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D74E}_{p}^{n+1}}{\unicode[STIX]{x0394}t}=(\unicode[STIX]{x1D70C}_{r}-1)\left[\frac{\boldsymbol{J}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D74E}_{p}^{n}}{\unicode[STIX]{x0394}t}-\unicode[STIX]{x1D74E}_{p}^{n}\times (\boldsymbol{J}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D74E}_{p}^{n})\right]+\int _{P}\boldsymbol{r}\times \left(\frac{\boldsymbol{u}^{\ast }}{\unicode[STIX]{x0394}t}-\unicode[STIX]{x1D740}^{n}\right)\,\text{d}\boldsymbol{x}. & \displaystyle\end{eqnarray}$$

Note that the above equations have been reformulated so that all the right-hand side terms are known quantities and consequently the particle velocities $\boldsymbol{U}^{n+1}$ and $\unicode[STIX]{x1D74E}_{p}^{n+1}$ are obtained without iteration. Then, the pseudo-body forces $\unicode[STIX]{x1D740}$ defined at the Lagrangian nodes are determined from

(2.10)$$\begin{eqnarray}\unicode[STIX]{x1D740}^{n+1}=\frac{\boldsymbol{U}^{n+1}+\unicode[STIX]{x1D74E}_{p}^{n+1}\times \boldsymbol{r}-\boldsymbol{u}^{\ast }}{\unicode[STIX]{x0394}t}+\unicode[STIX]{x1D740}^{n}.\end{eqnarray}$$

Finally, the fluid velocities $\boldsymbol{u}^{n+1}$ at the Eulerian nodes are corrected from

(2.11)$$\begin{eqnarray}\boldsymbol{u}^{n+1}=\boldsymbol{u}^{\ast }+\unicode[STIX]{x0394}t(\unicode[STIX]{x1D740}^{n+1}-\unicode[STIX]{x1D740}^{n}).\end{eqnarray}$$

In the above manipulations, a tri-linear function is used to transfer the fluid velocity from the Eulerian nodes to the Lagrangian nodes, and the pseudo-body force from the Lagrangian nodes to the Eulerian nodes.

2.2 Collision model

A particle–particle collision model is required to prevent the mutual penetration of particles. We assume that a repulsive force is activated when the distance between the surfaces of two spheroids is smaller than a critical value. This force acts on the two surface points of shortest distance between two particles along the connecter of these two points, and is transferred to the force on the particle mass centre plus a torque. We adopt an efficient iterative method proposed by Lin & Han (Reference Lin and Han2002) to determine the positions of two points of shortest distance between two spheroids. The force has the form

(2.12)$$\begin{eqnarray}\boldsymbol{F}_{ij}=F_{0}(1-d_{ij}/d_{c})\boldsymbol{n}_{ij},\end{eqnarray}$$

where $\boldsymbol{F}_{ij}$, $d_{ij}$ and $\boldsymbol{n}_{ij}$ are the repulsive force acting on particle $j$ from particle $i$, the gap distance and the unit vector parallel to the connecter of two shortest distance points; $d_{c}$ represents a cutoff distance and the repulsive force is activated when $d_{ij}<d_{c}$; $F_{0}$ is the magnitude of the force at contact. We set $d_{c}=h$ ($h$ being the fluid mesh size), and $F_{0}=10^{3}$. The motions of the particles due to the collision force (2.12) and due to the hydrodynamic force (2.8)–(2.11) are handled separately with a fractional-step scheme. The time step for the collision model is set to be one tenth of the latter (i.e. $\unicode[STIX]{x0394}t/10$) to circumvent the stiffness problem arising from the explicit integration scheme with a large value of $F_{0}$, as suggested by Glowinski et al. (Reference Glowinski, Pan, Hesla and Joseph1999). A collision between a particle and a wall is treated similarly with the coefficient $F_{0}$ in (2.12) doubled.

The repulsive force in our collision model is similar to the elastic (spring) force in the discrete element model (referred to as DEM) developed originally for the simulation of granular materials. In the DEM collision model, besides the spring-like repulsive force, the viscous damping force in the normal direction and the tangential friction force are also employed. In our previous work on the horizontal turbulent duct flow laden with heavy spherical particles (Lin et al. Reference Lin, Shao, Yu and Wang2017b), it was observed that the relative difference in the bulk velocity under the same pressure gradient obtained with the simple repulsive force model and the DEM model was around 4 % at a particle volume fraction of 2.36 % for the case where most particles settled down on the bottom wall and the difference became smaller for a lower particle settling effect. For the vertical channel flow here, the particles tend to migrate away from the wall and there are few particles in the near-wall region, therefore, the effect of the particle collision model at the same particle volume fraction of 2.36 % is expected to be less significant than that for the horizontal duct flow, since the particles located closer to the wall collide more frequently due to a higher velocity gradient.

Note that the repulsive force at contact $F_{0}$ is normalized with $\unicode[STIX]{x1D70C}_{f}u_{b}^{2}H^{2}$, and thus $F_{0}=10^{3}$ is much larger than the hydrodynamic force on the particles and contact without a gap between the particles is prevented. In one sense, one may think that there exist physically short-range repulsive forces between the particles (and between the particles and the wall) in our fluid–solid system.

2.3 Flow configuration and simulation settings

Figure 1.

Schematic diagram for the upward channel flow, with $x$, $y$ and $z$ representing the streamwise, transverse and spanwise coordinates, respectively.

The schematic diagram for the upward channel flow is shown in figure 1. The no-slip velocity boundary condition is imposed on the channel walls, and the periodic boundary condition is imposed in the streamwise and spanwise directions, respectively. The fluid is assumed to move upwards, and gravity points downwards. In the present study, we only consider heavy particles and refer to their gravity-driven downward relative motion with respect to the fluid as ‘sedimentation’ or ‘settling’. We denote the half-width of the channel as $H$. The computational domain spans $[0,8H]\times [-H,H]\times [0,4H]$. The bulk velocity of the suspension (i.e. the fluid–solid mixture) $u_{b}$ is kept constant in our simulations. We take $H$ as the characteristic length and the bulk velocity $u_{b}$ as the characteristic velocity, as mentioned earlier. The bulk Reynolds number defined by $Re_{b}=u_{b}H/\unicode[STIX]{x1D708}$ is 2873, which is determined from the bulk velocity of our single-phase flow with $h=H/64$ ($h$ being the mesh size) at the friction Reynolds number $Re_{\unicode[STIX]{x1D70F}}=u_{\unicode[STIX]{x1D70F}}H/\unicode[STIX]{x1D708}=180$ ($u_{\unicode[STIX]{x1D70F}}$ being the wall friction velocity). The friction velocity is defined as $u_{\unicode[STIX]{x1D70F}}=\sqrt{\unicode[STIX]{x1D70F}_{w}/\unicode[STIX]{x1D70C}_{f}}$, here $\unicode[STIX]{x1D70F}_{w}$ is the mean shear stress on the walls, and $\unicode[STIX]{x1D70C}_{f}$ denotes the fluid density. For our single-phase flow with $h=H/128$, which is used in the present study, $Re_{\unicode[STIX]{x1D70F}}=180.8$.

We consider spheroidal particles, which are characterized by three semi-axes $(a=b\neq c)$, with $c$ being the semi-symmetry axis. The aspect ratio $\unicode[STIX]{x1D706}=c/a$ distinguishes between prolate spheroids $(\unicode[STIX]{x1D706}>1)$, spheres $(\unicode[STIX]{x1D706}=1)$ and oblate spheroids $(\unicode[STIX]{x1D706}<1)$.

Table 1.

Physical parameter settings for the simulations of particulate upward flows in a vertical channel. In all cases, the particle equivalent radius is $a_{e}=H/20$.

The settings of the physical parameters for our numerical simulations are presented in table 1. For all particle-laden cases, the particle equivalent radius is $a_{e}=H/20$. We choose $\unicode[STIX]{x1D706}=2$ for the prolate particle and $\unicode[STIX]{x1D706}=1/3$ for the oblate particle. Throughout this study, the particle–fluid density ratio is $\unicode[STIX]{x1D70C}_{r}=2.0$. There are 2880 particles, corresponding to a particle volume fraction $\unicode[STIX]{x1D711}_{0}=2.36\,\%$. The length of the semi-axis in wall units computed with the friction velocity for each case is presented in table 1. We use the ratio of the Stokes free-fall velocity of a sphere $u_{s}$ to the bulk velocity $u_{b}$ to measure the settling effect of the particles in the vertical channel flow. The Stokes velocity $u_{s}$ is defined by $u_{s}=2(\unicode[STIX]{x1D70C}_{s}-\unicode[STIX]{x1D70C}_{f})ga^{2}/(9\unicode[STIX]{x1D707})$. The Froude number for the non-neutrally buoyant case is not an independent control parameter, and can be determined from the other control parameters

(2.13)$$\begin{eqnarray}Fr=\frac{gH}{u_{b}^{2}}=\frac{9(u_{s}/u_{b})}{(\unicode[STIX]{x1D70C}_{r}-1)(a/H)^{2}Re_{b}}.\end{eqnarray}$$

We refer to $u_{s}/u_{b}$ as the settling coefficient, and in the following, we often simplify it as $u_{s}$, since the dimensionless value of $u_{b}$ in our computations is unity. Two settling coefficients are considered $u_{s}/u_{b}=0.1$ and 0.3, which turn out to represent weak and moderate settling effects, respectively. The terminal settling velocity of a spherical particle in an unbounded domain is smaller than the Stokes velocity due to the fluid inertial effect, and can be determined from the Stokes velocity, the particle size and the bulk Reynolds number with the following equation:

(2.14)$$\begin{eqnarray}C_{D}{\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}_{f}u_{T}^{2}\unicode[STIX]{x03C0}a_{e}^{2}=6\unicode[STIX]{x03C0}a_{e}\unicode[STIX]{x1D707}u_{s},\end{eqnarray}$$

where the standard drag coefficient $C_{D}$ has the form

(2.15)$$\begin{eqnarray}C_{D}=\frac{24}{Re_{p}}(1+0.15Re_{p}^{0.687}),\end{eqnarray}$$

in which $Re_{p}$ represents the particle Reynolds number defined as $Re_{p}=2a_{e}u_{T}/\unicode[STIX]{x1D708}$. For the $a_{e}/H=0.05$ and $Re_{b}=2873$ values considered, $u_{T}/u_{b}=0.05125$ for $u_{s}/u_{b}=0.1$, and $u_{T}/u_{b}=0.1135$ for $u_{s}/u_{b}=0.3$. Then, $Re_{p}=14.7$ for $u_{s}/u_{b}=0.1$, and $Re_{p}=32.6$ for $u_{s}/u_{b}=0.3$. The settling coefficient, the terminal settling velocity and the difference between the fluid and solid mean velocities in the channel centre region for the spherical particle case are compared in table 2. The inter-phase mean velocity difference is smaller than the particle terminal velocity in the presence of fluid inertia, which can be easily understood, considering that both the lateral and streamwise velocity fluctuations can lead to the enhancement in the total drag compared to the no fluctuation case since the drag on the individual particle is proportional to $u^{\unicode[STIX]{x1D6FC}}$ (with $u$ being the slip velocity and $\unicode[STIX]{x1D6FC}>1$ due to fluid inertia), and the hydrodynamic interactions between particles can also hinder the particle settling. The inter-phase mean velocity difference being smaller than the particle terminal velocity was observed experimentally in an upward vertical turbulent channel flow (Kameyama et al. Reference Kameyama, Kanai, Kawashima and Ishima2014) and an isotropic turbulent flow (Doroodchi et al. Reference Doroodchi, Evans, Schwarz, Lane, Shah and Nguyen2008).

Table 2.

Comparison of the settling coefficient, the terminal settling velocity and the difference between the fluid and solid mean velocities in the channel centre region for the spherical particle case.

The grid number is $1024\times 256\times 512$, corresponding to a spatial resolution of $2a_{e}/h=12.8$, and $h^{+}=1.40$ in wall units. The dimensionless time step is 0.0008. The fluid statistics are obtained from the averaging of the data in the real fluid domain over a period of typically one thousand non-dimensional time units (i.e. approximately 63 large-eddy turnover times) after the statistically steady state is reached. The statistics for the solid-phase velocity and volume fraction are obtained from the averaging of the data in the solid domain (i.e. inside the particle boundary), whereas the orientation and rotational velocity statistics are calculated with the particle data.

The accuracy of our DF/FD code for the single-phase turbulence and good mesh convergence of the turbulence statistics for the particle-laden flows (even for the coarse mesh of 6.4 points per one particle diameter) have been demonstrated in Yu et al. (Reference Yu, Lin, Shao and Wang2016b). Our results on the turbulent channel flow laden with finite-size neutrally buoyant particles have been compared to lattice Boltzmann simulations using interpolated bounce back at the fluid–solid interfaces, and the two completely different numerical approaches yielded quantitatively similar results in general (Wang et al. Reference Wang, Peng, Guo and Yu2016). Our DF/FD code has been applied to various types of particle-laden turbulent flows including pipe flow (Wu et al. Reference Wu, Shao and Yu2011), duct flow (Lin et al. Reference Lin, Yu, Shao and Wang2017a,Reference Lin, Shao, Yu and Wangb) and channel flow laden with spherical particles (Shao et al. Reference Shao, Wu and Yu2012; Yu et al. Reference Yu, Lin, Shao and Wang2017) and spheroidal particles (Zhu et al. Reference Zhu, Yu and Shao2018), respectively. Thus, the validation of our code is not conducted here.

3.1 Fluid velocities

Figure 2.

Fluid mean velocity profiles of particle-laden turbulent channel flows. The inset shows the close up at the wall.

The fluid-phase mean velocity profiles for all seven cases studied are plotted in figure 2, with the inset showing the close up of the velocities in the immediate vicinity of the wall. From figure 2, in the bulk region (say $y>0.2H$, with $y$ being the distance away from the wall) the velocity profiles for $u_{s}=0.3$ are flatter than that for the single-phase flow, consistent with the observation that the velocity profile was largely flat in the bulk region when the particle mean settling velocity was comparable to the bulk velocity (Uhlmann Reference Uhlmann2008; Santarelli & Fröhlich Reference Santarelli and Fröhlich2015). Normally, a flatter velocity profile means a larger velocity gradient on the wall (i.e. larger flow friction). However, it is interesting that the velocity gradients for the spherical and prolate particle cases were smaller than that of the single-phase flow, as shown in the inset of figure 2. The inset of figure 2 also shows that the flow friction is largest for the spherical particles and smallest for the oblate particles in the case of $u_{s}=0.1$, as observed for the neutrally buoyant case (Zhu et al. Reference Zhu, Yu and Shao2018), whereas it is largest for the oblate particles and smallest for the spherical particles in the case of $u_{s}=0.3$, exactly opposite to the case of $u_{s}=0.1$. The flow friction coefficients are listed in table 3. This is the most important finding and the focus of the present study, and we will analyse the reason later.

Figure 3.

Profiles of the root-mean-square velocity components and the Reynolds shear stress for the fluid phase: (a) streamwise $u_{frms}$, (b) wall normal $v_{frms}$, (c) spanwise $w_{frms}$, (d) the Reynolds stress $-\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$.

Figure 3 shows the profiles of the root-mean-square (r.m.s.) velocity components and the Reynolds shear stress. The particle addition attenuates the peak streamwise r.m.s. velocity and enhances the transverse and spanwise r.m.s. velocities near the wall for $u_{s}=0.1$, as for the neutrally buoyant case (Shao et al. Reference Shao, Wu and Yu2012; Picano et al. Reference Picano, Breugem and Brandt2015; Wang et al. Reference Wang, Peng, Guo and Yu2016), whereas all r.m.s. velocity components and the Reynolds shear stress are attenuated at any transverse position for $u_{s}=0.3$, irrespective of particle shape. The effects of the particle shape on the fluid r.m.s. velocities and Reynolds stress are relatively small for $u_{s}=0.1$, and become significant for $u_{s}=0.3$. For the neutrally buoyant case, the oblate particles of $\unicode[STIX]{x1D706}=1/3$ caused a lower fluid Reynolds stress than the spherical particles and the prolate particles of $\unicode[STIX]{x1D706}=2$ (Zhu et al. Reference Zhu, Yu and Shao2018), which is one factor related to the drag reduction by the oblate particles for $\unicode[STIX]{x1D706}=1/3$. By contrast, for the non-neutrally buoyant case at $u_{s}=0.3$, all r.m.s. velocity components and the Reynolds stress are largest in the oblate case and smallest in the sphere case. For $u_{s}=0.1$, the r.m.s. velocities and Reynolds stress in the channel centre region become largest in the oblate case, due to the particle settling effect.

3.2 Solid-phase velocities

Figure 4.

(a) Solid-phase mean velocity profiles, and (b) difference between the fluid and solid mean velocities for upward vertical turbulent channel flows.

Figure 5.

Snapshots of spherical particles and low- and high-speed streaks in the plane of $y/H=-0.9$ for (a) $u_{s}=0.1$ and (b) $u_{s}=0.3$.

The mean particle-phase velocity and fluid–particle velocity difference profiles are plotted in figures 4(a) and 4(b), respectively. For all particle-laden cases, the solid mean velocities are smaller than the fluid mean velocities in the bulk region, but larger than the fluid mean velocities in the near-wall region, as observed experimentally in particle-laden upward vertical pipe flow (Tsuji et al. Reference Tsuji, Morikawa and Shiomi1984; Shokri et al. Reference Shokri, Ghaemi, Nobes and Sanders2017). The smaller particle velocity in the bulk region is due to the particle settling effect. The reason for the larger particle mean velocity near the wall was attributed to the facts that the particles can slip on the wall and are preferentially located in the fluid high-speed streaks (Tsuji et al. Reference Tsuji, Morikawa and Shiomi1984; Shokri et al. Reference Shokri, Ghaemi, Nobes and Sanders2017; Do-Quang et al. Reference Do-Quang, Amberg, Brethouwer and Johansson2014). From figure 4(a), the particle mean velocity profiles become flatter for a higher settling coefficient, with the spherical particle case being most pronounced. In the bulk region, the settling (or slip) velocity of the spherical particles is larger than the non-spherical particles, which can be explained by the observation that non-spherical settling particles tend to align their long axes perpendicular to the settling direction in the bulk region due to fluid inertia (see figure 11). The mean settling velocity of the oblate particles is smallest, since the projected (upwinding) area of the oblate particle with $\unicode[STIX]{x1D706}=1/3$ in the streamwise direction is larger than that of the prolate particle with $\unicode[STIX]{x1D706}=2$ and that of the spherical particle. Figure 4(b) shows that the particles lead the fluid more significantly on average in the near-wall region for $u_{s}=0.3$, compared to the case of $u_{s}=0.1$, although the settling velocity of a particle in quiescent fluid is larger for $u_{s}=0.3$. To seek the reason, snapshots of spherical particles and velocity streaks in the plane of $y/H=-0.9$ for $u_{s}=0.1$ and $u_{s}=0.3$ are plotted in figure 5. Figure 5 presents an evidence that the particles are preferentially located in the fluid high-speed streaks, and the comparison between figure 5(a,b) shows that such a preference is stronger for a higher settling coefficient; almost all particles are located in high-speed streaks for $u_{s}=0.3$ at the time randomly chosen.

Figure 6.

Relative mean velocities of the fluid in an annular band (cut shell) of thickness $a/4$ with respect to the spherical particle as a function of the distance between the shell and the particle for different $y$-bands of width $0.05H$, and (a) $u_{s}=0.1$ and (b) $u_{s}=0.3$. (c) Comparison between the relative mean velocity at $\unicode[STIX]{x1D6FF}r=9.875a$ and the inter-phase mean velocity difference, and (d) schematic diagram of the computational model for the relative mean velocity of the fluid with respect to the particle. The averaging is performed for the particles located in each $y$-band during statistically steady state; the legend ‘0.05–0.1’ in (a) means the $y$-band of $0.05H<y<0.1H$.

The time for the flow field in figure 5 is randomly chosen, however, it is more rigorous to inspect the relative mean velocity of the surrounding fluid with respect to the particles. Cisse, Homann & Bec (Reference Cisse, Homann and Bec2013) and Fornari et al. (Reference Fornari, Picano, Sardina and Brandt2016b) calculated the relative mean velocity by averaging the fluid velocity in the region of a spherical shell surrounding each particle. For turbulent channel flow with an inhomogeneous velocity distribution (i.e. non-constant shear rate distribution), it is not appropriate to take an entire spherical shell as the averaging region. Since the particle velocity is the velocity of the particle centre, it is better to compare the particle velocity with the velocity of the fluid in the region with the $y$-position close to that of the particle centre. Figure 6(d) shows the schematic diagram of our computational model for the relative mean velocity. The averaging region is a part of a spherical shell of thickness $a/4$ ($a$ being the particle radius) between two planes located at $y_{0}\pm d_{y}$, where $y_{0}$ is the $y$-position of the particle centre and $d_{y}$ is the half-height of the averaging region (i.e. an annular band). Zhu et al. (Reference Zhu, Yu, Shao and Deng2020) showed that $d_{y}=0.25a$ is a good choice, and therefore this value is chosen here. The mean relative velocities of the fluid with respect to the spherical particles as a function of the distance between the cut shell and the particle $\unicode[STIX]{x1D6FF}r$ for $u_{s}=0.1$ and $u_{s}=0.3$ are plotted in figures 6(a) and 6(b), respectively. Figure 6(c) shows the comparison between the relative mean velocity at $\unicode[STIX]{x1D6FF}r=9.875a$ and the inter-phase mean velocity difference. The relative mean velocity agrees well with the inter-phase mean velocity difference, except for $y\leqslant 0.125H$ at $u_{s}=0.3$. The discrepancy might be caused by the insufficient particle sample number, since there are few particles near the wall for $u_{s}=0.3$. However, the averaging period is around 1000 time units and the results are observed to change little when the averaging time is further increased.

From figure 6(a,b), the relative mean velocity of the fluid first increases with increasing distance away from the particle roughly along the $x$–$z$ plane for all $y$-positions in the channel, due to the particle sedimentation effect, and then deceases with increasing distance roughly for $y<0.5H$, which indicates that the settling particles near the wall are preferentially located in the high-speed regions. The decrease in the mean relative velocity with the distance is more pronounced for $u_{s}=0.3$ particularly in the near-wall region, and this is the statistical evidence that the particles have a stronger preference to be located in the high-speed region for a higher settling coefficient. In the near-wall region of $y<0.15H$, the particle phase-averaged velocity is larger than the fluid phase-averaged velocity, as shown in figure 4(b), and also larger than the mean velocity of the far-field fluid, as shown in figure 6(a,b), however, the particles still lag the surrounding near-field fluid due to the particle sedimentation effect (figure 6a,b). It has been recognized that the Saffman lift force can cause the particles which lag the fluid to migrate towards the channel centre in upward vertical channel flow (Yu, Phan-Thien & Tanner Reference Yu, Phan-Thien and Tanner2004; Marchioli et al. Reference Marchioli, Picciotto and Soldati2007). Thus, the Saffman effect can explain the lateral migration of the particles towards the channel centre in the entire channel.

The results from the point-particle simulations without consideration of the Saffman lift force showed that gravity had important effects on the concentration distribution of the particles with large inertia. Both the results of Marchioli et al. (Reference Marchioli, Picciotto and Soldati2007) and Nilsen et al. (Reference Nilsen, Andersson and Zhao2013) indicated that the inertial particles migrated towards the channel centre in a downward turbulent channel flow, opposite to the prediction by the Saffman effect. For the upward turbulent channel flow, migration of the inertial particles towards the wall was observed by Nilsen et al. (Reference Nilsen, Andersson and Zhao2013), and an almost constant particle volume fraction in the bulk region was observed by Marchioli et al. (Reference Marchioli, Picciotto and Soldati2007) for Stokes numbers up to 25. Their Stokes number was defined as $St=2\unicode[STIX]{x1D70C}_{r}(au_{\unicode[STIX]{x1D70F}}/\unicode[STIX]{x1D708})^{2}/9$, which equals 36 for our simulations at $\unicode[STIX]{x1D70C}_{r}=2.0$, $a/H=0.05$ and $Re_{\unicode[STIX]{x1D70F}}=180$. It was observed in the interface-resolved simulations that the finite-size heavy particles migrated towards the walls in a downward flow, as predicted by the Saffman effect (Zhu et al. Reference Zhu, Yu, Shao and Deng2020). Fornari et al. (Reference Fornari, Formenti, Picano and Brandt2016a) showed that the neutrally buoyant particles with $\unicode[STIX]{x1D70C}_{r}=2.0$ did not migrate laterally, as for the case of $\unicode[STIX]{x1D70C}_{r}=1.0$, although pronounced migration towards the channel centre was observed for $\unicode[STIX]{x1D70C}_{r}=10.0$ (Fornari et al. Reference Fornari, Formenti, Picano and Brandt2016a; Yu et al. Reference Yu, Lin, Shao and Wang2017). We believe that the Saffman effect is responsible for the particle migration towards the channel centre here, since the effects of the turbulence (via the drag force) and particle inertia cannot cause the particle migration behaviour observed here.

Figure 7.

Schematic diagram explaining why particles are distributed preferentially in high-speed streaks in a vertical channel. The arrows represent the flow direction (i.e. streamwise direction).

Now, we attempt to explain why heavy particles are distributed preferentially in high-speed streaks in an upward channel flow. It was observed that the particles in the bottom half-channel moved more slowly than the fluid on average, while they were faster in the upper half-channel for the turbulent channel flow in a horizontal channel (Shao et al. Reference Shao, Wu and Yu2012). Shao et al. (Reference Shao, Wu and Yu2012) argued that if the particle volume fraction decreases with increasing distance away from the wall in a half-channel, the particles tend to be located in the low-speed streaks, and if it increases with increasing distance away from the wall, the particles tend to be located in the high-speed streaks, due to the entrainment effect on the particles by the streamwise vortex pairs. Here, in the upflow vertical channel, the particle volume fraction in the near-wall region is lower than in the bulk region (see figure 8) and consequently the particles are more likely located in the high-speed region, since the number of the particles entrained into the high-speed region by the large-scale vortex pair is higher than that of the particles entrained out of the high-speed region, as sketched in figure 7. The particle concentration distribution for $u_{s}=0.3$ is more inhomogeneous than for $u_{s}=0.1$ (figure 8), and this explains why the particles have a stronger preferential distribution in the high-speed steaks for $u_{s}=0.3$, which leads to a higher particle mean velocity for $u_{s}=0.3$, compared to $u_{s}=0.1$, as observed in figure 4(a).

Figure 8.

Distribution of the particle volume fraction across the channel.

From figure 8, the migration towards the channel centre is most pronounced for the spherical particles and weakest for the oblate particles for both settling coefficients. However, for $u_{s}=0.1$, the particle volume fraction in the near-wall region of $y<0.05H$ is largest for the sphere case and smallest for the oblate case, which is opposite to the observation in the bulk region, but similar to the neutrally buoyant case (Zhu et al. Reference Zhu, Yu and Shao2018). The larger concentration near the wall for the neutrally buoyant spherical particles was regarded as a primary reason for the larger flow friction in this case, compared to the non-spherical particle case (Ardekani et al. Reference Ardekani, Costa, Breugem, Picano and Brandt2017). The stronger migration of the spherical particles settling in an upflow may be related to two factors. One is the direct effect of the particle shape: the settling non-spherical particles tend to align their long axes with the streamwise direction in the near-wall region, and perpendicular to the streamwise direction in the bulk region, as shown later in figures 10 and 11, and their spanwise rotation rates are reduced, as shown by Ardekani et al. (Reference Ardekani, Costa, Breugem, Picano and Brandt2017) for the neutrally buoyant case. The effect of the particle rotation on the Saffman force is complicated. Yu et al. (Reference Yu, Phan-Thien and Tanner2004) observed that the sphere rotation produced enhanced lift forces pointing towards the pipe axis in the neutrally buoyant and downflow cases for laminar pipe flow, whereas it had a negligible effect on the migration of the sphere in laminar upward pipe flow. The effect of the particle shape on the lift force for laminar flow has not been investigated, to our knowledge. The other factor is the indirect effect of the particle shape: the large-scale vortices are weakened more significantly by the spherical particles, particularly for $u_{s}=0.3$ (figures 3 and 13), and consequently, the diffusion effect of the turbulence on the particle distribution that tends to flatten the particle concentration distribution profile is weaker for spherical particles. In figure 8, the particle concentration distribution profiles for the non-spherical particles in the bulk region are much flatter than those for the spherical particles. Considering that the turbulence intensities are comparable for three types of particles for $u_{s}=0.1$ (figure 3) but the difference in the particle volume fraction profile in this case is still pronounced, both direct and indirect effects of the particle shape are expected to be important to the particle migration.

Figure 9.

Profiles of the r.m.s. velocity and Reynolds shear stress for the particle phase: (a) streamwise $u_{prms}$, (b) wall normal $v_{prms}$, (c) spanwise $w_{prms}$, (d) Reynolds stress $-\langle u_{p}^{\prime }v_{p}^{\prime }\rangle$.

The solid-phase r.m.s. velocities and kinematic Reynolds shear stresses (here meaning $-\langle u_{p}^{\prime }v_{p}^{\prime }\rangle$ without the density) for all cases are plotted in figure 9. The intensity of the solid-phase velocity fluctuation generally decreases with increasing particle settling coefficient, similar to the fluid velocity fluctuation. The solid-phase r.m.s. velocities are larger than those of the fluid in the near-wall region due to the collision between the particles and the wall. For the same settling coefficient, the particle-phase r.m.s. velocities in the bulk region for the oblate case is largest, consistent with the results of the fluid r.m.s. velocities in figure 3.

3.3 Particle orientation and rotation

Figure 10.

Probability density functions (PDFs) of the particle orientation in the near-wall region (a,c,e) and the central region (b,d,f) of the channel: (a) $\cos \unicode[STIX]{x1D6FC}$, (b) $|\text{cos}\unicode[STIX]{x1D6FC}|$, (c) $|\text{cos}\unicode[STIX]{x1D6FD}|$, (d) $|\text{cos}\unicode[STIX]{x1D6FD}|$, (e) $|\text{cos}\unicode[STIX]{x1D6FE}|$ and (f) $|\text{cos}\unicode[STIX]{x1D6FE}|$.

The orientation of a spheroidal particle is measured with the direction cosines ($|\text{cos}\,\unicode[STIX]{x1D6FC}|$, $|\text{cos}\,\unicode[STIX]{x1D6FD}|$, $|\text{cos}\,\unicode[STIX]{x1D6FE}|$), where $\unicode[STIX]{x1D6FC}$, $\unicode[STIX]{x1D6FD}$, $\unicode[STIX]{x1D6FE}$ represent the angles between the particle symmetry axis and the streamwise, wall-normal and spanwise directions, respectively. Figure 10 shows the probability density functions (PDFs) of the direction cosines ($\cos \unicode[STIX]{x1D6FC}$, $|\text{cos}\,\unicode[STIX]{x1D6FD}|$, $|\text{cos}\,\unicode[STIX]{x1D6FE}|$) in the near-wall region of $y<2a_{e}$ (i.e. $0.1H$) and ($|\text{cos}\,\unicode[STIX]{x1D6FC}|$, $|\text{cos}\,\unicode[STIX]{x1D6FD}|$, $|\text{cos}\,\unicode[STIX]{x1D6FE}|$) in the central region of the channel (the distance away from the centre plane being smaller than $0.1H$). From figure 10(a,c,e), for the near-wall region, the orientation behaviour of the spheroidal particles settling in an upward vertical channel flow is similar to the neutrally buoyant case. Figure 10(a) shows that the most probable $\cos \unicode[STIX]{x1D6FC}$ for the oblate particles in the near-wall region is approximately $-0.35$, corresponding to $\unicode[STIX]{x1D6FC}\approx 110^{\circ }$. Note that $\unicode[STIX]{x1D6FC}$ is defined in the interval of $(0,180^{\circ })$. In other words, the symmetry axis is defined as always pointing away from the channel wall. Our result implies that the most probable orientation of the long (major) axis of the oblate particles is not exactly the streamwise direction, but has an inclination angle of approximately $20^{\circ }$ with the streamwise direction, as observed in the neutrally buoyant case (Zhu et al. Reference Zhu, Yu and Shao2018). From figure 10(a), the most probable orientation of the prolate particles is close to the streamwise direction, with the PDF of a positive inclination angle being larger than that of a negative inclination angle. Figure 10(b) shows that the oblate particles tend to align their symmetry axes with the wall-normal direction (with a deviation of approximately $20^{\circ }$). The sedimentation effect attenuates the preference of the orientation and makes the orientation more isotropic for particles near the wall, as shown in figure 10(a,c,e).

For the channel centre region, the orientation of the spheroidal particles is almost isotropic for $u_{s}=0.1$, as shown in figure 10(b,d,f) and observed in the neutrally buoyant case (Mortensen et al. Reference Mortensen, Andersson, Gillissen and Boersma2008b; Zhu et al. Reference Zhu, Yu and Shao2018). However, for $u_{s}=0.3$, the sedimentation effect results in a preferential orientation: the spheroids tend to align their long axes perpendicular to the streamwise direction, with almost isotropic projected orientation in the transverse–spanwise (i.e. $y$–$z$) plane, as indicated by the peak of the p.d.f. of $|\text{cos}\,\unicode[STIX]{x1D6FC}|$ occurring at $|\text{cos}\,\unicode[STIX]{x1D6FC}|=1$ in figure 10(b) and largely the same p.d.f. for $|\text{cos}\,\unicode[STIX]{x1D6FD}|$ and $|\text{cos}\,\unicode[STIX]{x1D6FE}|$ in figure 10(d,f). It has been recognized that a non-spherical particle settling in a Newtonian fluid turns its broad side perpendicular to the settling direction due to the fluid inertial effect (Huang et al. Reference Huang, Feng and Joseph1994). As mentioned earlier, the particle Reynolds number based on the slip velocity is $Re_{p}=14.7$ for $u_{s}=0.1$, and $Re_{p}=32.6$ for $u_{s}=0.3$. Our results indicate that the fluid inertial effect on the particle orientation is overshadowed by the turbulence effect at $Re_{p}=14.7$ and prevails over the turbulence effect at $Re_{p}=32.6$ in the bulk region.

Figure 11.

Side view of the snapshots of particles in the left half-channel: (a) spheres ($u_{s}=0.1$), (b) spheres ($u_{s}=0.3$), (c) oblate particles ($u_{s}=0.1$), (d) oblate particles ($u_{s}=0.3$), (e) prolate particles ($u_{s}=0.1$) and (f) prolate particles ($u_{s}=0.3$).

The snapshots of particles in the lower half-channel are shown in figure 11. Figure 11 provides the direct evidence for the results presented above. It can be seen that the particle settling effect drives the particles away from the wall at $u_{s}=0.3$, particularly for the sphere case. The spheroids tend to align their long axes with the streamwise direction, and the oblate particles tend to align their symmetry axes with the wall-normal direction, with a positive preferential inclination angle between the long axis and the streamwise direction in the vicinity of the wall. In the bulk region, the settling effect causes the particles to turn their long axes perpendicular to the streamwise direction, which is more pronounced for a higher settling coefficient.

Figure 12.

Particle mean rotation rates: (a) tumbling components and (b) spinning components. The insets show close ups for $y>0.2H$.

We decompose the particle rotation rate ($\langle \unicode[STIX]{x1D714}_{i^{\prime }}\unicode[STIX]{x1D714}_{i^{\prime }}\rangle$) into the spinning and tumbling components. The spinning component represents the rotation about the symmetry axis, defined as $\langle \unicode[STIX]{x1D714}_{z^{\prime }}\unicode[STIX]{x1D714}_{z^{\prime }}\rangle$, and the tumbling component represents the rotation about the other two axes, defined as $0.5(\langle \unicode[STIX]{x1D714}_{x^{\prime }}\unicode[STIX]{x1D714}_{x^{\prime }}\rangle +\langle \unicode[STIX]{x1D714}_{y^{\prime }}\unicode[STIX]{x1D714}_{y^{\prime }}\rangle )$, where $z^{\prime }$ denotes the symmetry axis, and $x^{\prime }$ and $y^{\prime }$ the other two axes in the body-fixed frame (Zhu et al. Reference Zhu, Yu and Shao2018). Figure 12 shows the particle mean tumbling and spinning velocities. The orientation and rotation rates of the spherical particles are computed with the same method for the non-spherical particles, after the body-fixed frame is established arbitrarily at the initial time. Due to the isotropic shape, the mean tumbling velocity of the spheres is nearly the same as the spinning velocity in figure 12. The results in figure 12 show that the prolate particles have higher spinning velocities and lower tumbling velocities in the bulk region, compared to the oblate particles, as in the neutrally buoyant case (Zhu et al. Reference Zhu, Yu and Shao2018). This was explained by the preferential alignment of the particle long axis with the local fluid vorticity direction for the point-particle case (Zhao et al. Reference Zhao, Challabotla, Andersson and Variano2015b). It is observed in the insets of figure 12 that the particle settling effect reduces both the spinning and tumbling rates of all three types of particles.

3.4 Vortex structure

Figure 13.

Vortex structures of the particle-laden vertical channel flows for: (a) spheres ($u_{s}=0.1$), (b) spheres ($u_{s}=0.3$), (c) oblate particles ($u_{s}=0.1$), (d) oblate particles ($u_{s}=0.3$), (e) prolate particles ($u_{s}=0.1$) and (f) prolate particles ($u_{s}=0.3$). The colour of the vortices represents the fluid streamwise velocity and the brighter means a larger streamwise velocity. The particles in a half-domain are shown.

Figure 13 shows the typical vortex structures in the lower half-channel, identified with the $Q$ criterion ($Q=1$) defined by Hunt, Wray & Moin (Reference Hunt, Wray and Moin1988),

(3.1)$$\begin{eqnarray}Q={\textstyle \frac{1}{2}}(\unicode[STIX]{x1D6FA}_{ij}\unicode[STIX]{x1D6FA}_{ij}-\unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{ij}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}_{ij}=(u_{i,j}-u_{j,i})/2$ and $\unicode[STIX]{x1D61A}_{ij}=(u_{i,j}+u_{j,i})/2$ are respectively the rotation rate tensor and the strain rate tensor. From figure 13, the particles induce vortex rings around them at $u_{s}=0.3$. It is known that flow separation occurs at the critical Reynolds number of approximately 20 for flow over a sphere (Chang & Maxey Reference Chang and Maxey1994). For $u_{s}=0.1$, the particle Reynolds number is 14.7 and there is no flow separation. This may explain why no particle-induced vortex structures are identified with $Q=1$ in the bulk region for $u_{s}=0.1$ in figure 13 and the effect of the particles on the turbulent channel flow is relatively weak, as observed earlier for the fluid statistics (figures 2 and 3). The large-scale vortices are not significantly weakened by the particles for $u_{s}=0.1$, as shown in figure 13(a,c,e). By contrast, the large-scale vortices are suppressed significantly in the spherical and prolate cases for $u_{s}=0.3$, from the comparison between figure 13(a,b), and comparison between figure 13(e,f). The comparison of figure 13(b,d,f) shows that the particle-induced vortices are strongest for the spherical particles, and weakest for the oblate particles in the bulk region at $u_{s}=0.3$. The reason should be that the preferential orientation of non-spherical particles with their long axes perpendicular to the settling direction reduces their settling velocity, as shown in figure 4(b). The oblate particles for $\unicode[STIX]{x1D706}=1/3$ have a largest upwinding area when settling, and thus the slip velocity is smallest, compared to the spherical particles and the prolate particles of $\unicode[STIX]{x1D706}=2$. The particle-induced vortices should be responsible for the suppression of the large-scale vortices.

3.5 Flow friction and total drag

In this subsection, we will explore the reasons for the effects of the particles on the friction drag. Based on the spatial averaging theorem, the mean momentum equation of turbulent channel flow can be derived as below (Picano et al. Reference Picano, Breugem and Brandt2015; Yu et al. Reference Yu, Lin, Shao and Wang2017)

(3.2)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{d}}{\text{d}y}(\unicode[STIX]{x1D711}_{f}\langle \unicode[STIX]{x1D70E}_{f}\rangle _{xy}+\unicode[STIX]{x1D711}_{s}\langle \unicode[STIX]{x1D70E}_{s}\rangle _{xy})+\left(-\frac{\text{d}p_{e}}{\text{d}x}\right)\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{\text{d}}{\text{d}y}(\unicode[STIX]{x1D711}_{f}\unicode[STIX]{x1D70C}_{f}\langle -u_{f}^{\prime }v_{f}^{\prime }\rangle +\unicode[STIX]{x1D711}_{s}\unicode[STIX]{x1D70C}_{s}\langle -u_{s}^{\prime }v_{s}^{\prime }\rangle )-\unicode[STIX]{x1D711}_{s}(\unicode[STIX]{x1D70C}_{r}-\unicode[STIX]{x1D70C}_{f})g=0,\end{eqnarray}$$

where $u$ and $v$ are the local velocity components in the $x$ and $y$ directions, respectively, with the subscript ‘$f$’ meaning the fluid phase and ‘$s$’ meaning the solid phase; $\unicode[STIX]{x1D70E}_{s}$ is the solid inner stress, $g$ is the gravity, $\unicode[STIX]{x1D711}_{f}$ and $\unicode[STIX]{x1D711}_{s}$ are the fluid and solid volume fractions at a given $y$ position and $-\text{d}p_{e}/\text{d}x$ denotes the pressure gradient. Note that, here, the particle Reynolds stress means the solid-phase Reynolds stress; namely, the velocities at the Eulerian grids inside the particle boundaries rather than the particle translational velocities are used for the computation of the particle Reynolds stress. For the two-phase upward vertical channel flow, the relationship between $p_{e}$ and $\unicode[STIX]{x1D70F}_{w}$ is $-\text{d}p_{e}/\text{d}x=\unicode[STIX]{x1D70F}_{w}/H+\unicode[STIX]{x1D711}_{0}(\unicode[STIX]{x1D70C}_{s}-\unicode[STIX]{x1D70C}_{f})g$. Equation (3.2) can be rewritten as follows:

(3.3)$$\begin{eqnarray}\displaystyle -\frac{\unicode[STIX]{x1D70F}_{w}}{H} & = & \displaystyle \frac{\text{d}}{\text{d}y}(\unicode[STIX]{x1D711}_{f}\langle \unicode[STIX]{x1D70E}_{f}\rangle _{xy}+\unicode[STIX]{x1D711}_{s}\langle \unicode[STIX]{x1D70E}_{s}\rangle _{xy})+(\unicode[STIX]{x1D711}_{0}-\unicode[STIX]{x1D711}_{s})(\unicode[STIX]{x1D70C}_{r}-\unicode[STIX]{x1D70C}_{f})g\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\text{d}}{\text{d}y}(\unicode[STIX]{x1D711}_{f}\unicode[STIX]{x1D70C}_{f}\langle -u_{f}^{\prime }v_{f}^{\prime }\rangle +\unicode[STIX]{x1D711}_{s}\unicode[STIX]{x1D70C}_{s}\langle -u_{s}^{\prime }v_{s}^{\prime }\rangle ).\end{eqnarray}$$

The contributions of different stresses to the friction coefficient can be obtained by applying triple integration (Fukagata, Iwamoto & Kasagi Reference Fukagata, Iwamoto and Kasagi2002), i.e. $\int _{0}^{1}\text{d}y\int _{0}^{y}\text{d}y\int _{0}^{y}\text{d}y$, to (3.3):

(3.4)$$\begin{eqnarray}\displaystyle C_{f} & = & \displaystyle \underbrace{\frac{6}{Re_{b}}\int _{0}^{1}(1-y)\unicode[STIX]{x1D711}_{f}\left(\frac{\text{d}u^{\ast }}{\text{d}y}\right)\text{d}y}_{C_{fv}}+\underbrace{6\int _{0}^{1}(1-y)(\unicode[STIX]{x1D711}_{f}\langle -u_{f}^{\ast \prime }u_{f}^{\ast \prime }\rangle )\,\text{d}y}_{C_{fR}}\nonumber\\ \displaystyle & & \displaystyle +\,\underbrace{6\int _{0}^{1}(1-y)\frac{\unicode[STIX]{x1D711}_{s}}{\unicode[STIX]{x1D70C}_{f}u_{b}^{2}}\langle \unicode[STIX]{x1D70E}_{s}\rangle _{xy}\,\text{d}y}_{C_{pI}}+\underbrace{6\int _{0}^{1}(1-y)(\unicode[STIX]{x1D711}_{s}\unicode[STIX]{x1D70C}_{r}\langle -u_{s}^{\ast \prime }u_{s}^{\ast \prime }\rangle )\,\text{d}y}_{C_{pR}}\nonumber\\ \displaystyle & & \displaystyle +\,\underbrace{\left(-3\int _{0}^{1}(1-y)^{2}(\unicode[STIX]{x1D711}_{s}-\unicode[STIX]{x1D711}_{0})(\unicode[STIX]{x1D70C}_{r}-1)Fr\,\text{d}y\right)}_{C_{g}},\end{eqnarray}$$

where $C_{f}$ is the friction coefficient defined by $C_{f}=2\unicode[STIX]{x1D70F}_{w}/\unicode[STIX]{x1D70C}_{f}u_{b}^{2}$, and $y$ is the dimensionless coordinate normalized with $H$, with the wall position shifted to $y=0$; $u^{\ast }$ and $v^{\ast }$ denote the dimensionless velocities normalized by $u_{b}$, and $Fr$ represents the Froude number defined by $Fr=gH/u_{b}^{2}$. We define the terms of the fluid viscous stress, the fluid Reynolds stress, the particle inner stress, the particle Reynolds stress and gravity in (3.4) as $C_{fV}$, $C_{fR}$, $C_{pi}$, $C_{pR}$ and $C_{g}$, respectively. To derive (3.4), the following relations are used:

(3.5)$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D719}_{f}\langle \unicode[STIX]{x1D70E}_{f}\rangle _{xy}+\unicode[STIX]{x1D719}_{s}\langle \unicode[STIX]{x1D70E}_{s}\rangle _{xy})|_{y=0}=\unicode[STIX]{x1D70F}_{w} & \displaystyle\end{eqnarray}$$

(3.6)$$\begin{eqnarray}\displaystyle & \displaystyle \int _{0}^{1}\int _{0}^{y}f\,\text{d}y\,\text{d}y=\int _{0}^{1}(1-y)f\,\text{d}y;\int _{0}^{1}\int _{0}^{y}\int _{0}^{y}f\,\text{d}y\,\text{d}y\,\text{d}y=\frac{1}{2}\int _{0}^{1}(1-y)^{2}f\,\text{d}y. & \displaystyle\end{eqnarray}$$

The contributions of the different stresses in (3.4) are presented in table 3. For the neutrally buoyant case, the smaller flow friction for the non-spherical particles is primarily caused by the significantly reduced contribution from the particle inner stress, compared to the spherical particles (Ardekani et al. Reference Ardekani, Costa, Breugem, Picano and Brandt2017; Zhu et al. Reference Zhu, Yu and Shao2018). The drag reduction by the non-spherical particles compared to the single-phase flow is due to the reductions in the contributions from both particle inner stress and total Reynolds stress (Ardekani et al. Reference Ardekani, Costa, Breugem, Picano and Brandt2017; Eshghinejadfard et al. Reference Eshghinejadfard, Zhao and Thévenin2018; Zhu et al. Reference Zhu, Yu and Shao2018; Ardekani & Brandt Reference Ardekani and Brandt2019). For $u_{s}=0.1$ here, the reason for the smaller flow friction of the non-spherical particles compared to the spherical particles is due to the reductions in the particle inner stress, total Reynolds stress and gravity term, as shown in table 3. However, drag reduction by the non-spherical particles does not occur, because the reduction in the Reynolds stress is smaller than the extra contribution from the particle inner stress and gravity term, compared to single-phase flow. The particle inner stress is expected to be increased by the particle settling effect, because the drag on the particles is enhanced due to the increase in the slip velocity. The contribution of the particle inner stress becomes significantly larger, as $u_{s}$ is increased from 0.1 to 0.3. Nevertheless, the particle sedimentation effect at $u_{s}=0.3$ attenuates greatly the turbulence (figure 13) and gives rise to the substantial reduction in the fluid Reynolds stress, resulting in the reduction in the flow friction, particularly for the spherical and prolate particles, as shown in table 3.

Table 3.

Contributions of the different stresses in (3.4) to the $C_{f}$ (the data unit is $\times 10^{-3}$).

From table 3, the flow frictions for the spherical and prolate particles at $u_{s}=0.3$ are smaller than that of the single-phase flow, but this does not mean drag reduction, since the flow drag is normally measured by the pressure drop at the same flow rate. For the upward vertical channel flow studied, the pressure drop is the sum of the flow friction and the buoyant force on the particles, and the relation is: $-\text{d}p_{e}/\text{d}x=\unicode[STIX]{x1D70F}_{w}/H+\unicode[STIX]{x1D711}_{0}(\unicode[STIX]{x1D70C}_{s}-\unicode[STIX]{x1D70C}_{f})g$. The pressure gradients for all cases are presented in table 4. Although the particle volume fraction ($2.36\,\%$) and the settling coefficient are not large, the buoyant force on the particles is comparable to the flow friction for $u_{s}=0.3$, and thus the particle sedimentation effect does not cause drag reduction. However, it is still interesting that the sedimentation effect can lead to the transition of the flow drag of the spherical particles from being larger to being smaller than that of the non-spherical particles, as the settling coefficient increases.

Table 4.

Pressure gradients and friction Reynolds numbers for the seven cases listed in table 1 ($(-\text{d}p_{e})/(\text{d}x)$, $\unicode[STIX]{x1D70F}_{w}/H$ and $\unicode[STIX]{x1D711}_{0}(\unicode[STIX]{x1D70C}_{s}-\unicode[STIX]{x1D70C}_{f})g$ are normalized by $\unicode[STIX]{x1D70C}_{f}u_{b}^{2}/H$ and the data unit is $\times 10^{-3}$).