Hydrodynamic forces on assemblies of non-spherical particles: orientation and voidage effects

TL;DR

This paper develops a comprehensive model for hydrodynamic forces on assemblies of non-spherical particles, highlighting the key role of particle incident angle and orientation-independent forces, applicable across viscous and inertial regimes.

Contribution

It introduces new parameters to characterize flow and particle configuration, and proposes closure models for drag, lift, and torque in non-spherical particle assemblies.

Findings

Hydrodynamic forces are independent of mutual particle orientations.

Drag follows sine-squared scaling, and lift follows sine-cosine scaling, even in multiparticle systems.

Pressure drop can be estimated using orientation distribution and force scaling laws.

Abstract

This work provides a recipe for creating drag, lift and torque closures for static assemblies of axisymmetric, non-spherical particles. Apart from Reynolds number $Re$ and solids volume fraction $\epsilon_s$, we propose four additional parameters to characterize the flow through non-spherical particle assemblies. Two parameters consider the mutual orientations of particles (the orientation tensor eigenvalues $S_1$ and $S_2$) and two angles represent the flow direction (polar and azimuthal angles $\alpha$ and $\beta$). Interestingly, we observe that the hydrodynamic forces on the particles are independent of the mutual particle orientations. Rather, the most important parameter representing the particle configuration itself is the incident angle $\phi$ of the individual particles with respect to the incoming flow. Moreover, we observe that our earlier finding of sine-squared scaling of drag for isolated particles (Sanjeevi & Padding 2017) holds on average even for a multiparticle system in both the viscous and inertial regimes. Similarly, we observe that the average lift for a multiparticle system follows sine-cosine scaling, as is observed for isolated particles. Such findings are very helpful since the pressure drop of a packed bed or porous media can be computed just with the knowledge of orientation distribution of particles and their drag at $\phi=0^\circ$ and $\phi=90^\circ$ for a given $Re$ and $\epsilon_s$. With the identified dependent parameters, we propose drag, lift and torque closures for multiparticle systems.

Tables

Table 1. Table 1 : Details of the simulation parameters used in our simulations in LB units. L D subscript 𝐿 𝐷 L_{D} denotes the side length of the cubic domain. The range of d e ​ q subscript 𝑑 𝑒 𝑞 d_{eq} specified is respectively for 0.1 ≤ ϵ s ≤ 0.5 0.1 subscript italic-ϵ 𝑠 0.5 0.1\leq\epsilon_{s}\leq 0.5 .

$\text{Rey}$	$L_D$	$d_{eq}$	$\nu$
$0.1\le \text{Rey}\le 10$	288	$28.36-48.5$	1.3/3
	576	$56.72-97.0$	0.1 to 0.08/3
300	576	$56.72-97.0$	0.04/3
600	576	$56.72-97.0$	0.015/3
1000	768	$75.63-129.3$	0.01/3

Table 2. Table 2 : Coefficients of the fits for F D subscript 𝐹 𝐷 F_{D} and F L subscript 𝐹 𝐿 F_{L}

	$F_D$
Coefficients	$\varphi =0^{\circ }$	$\varphi ={90}^{\circ }$	$F_{L,mag}$
$a$	2	3	0.85
$b$	11.3	17.2	5.4
$c$	0.69	0.79	0.97
$d$	0.77	3	0.75
$e$	0.42	11.12	-0.92
$f$	4.84	11.12	2.66
$g$	0	0.57	1.94

Table 3. Table 3 : Coefficients of the fits for T P , m ​ a ​ g subscript 𝑇 𝑃 𝑚 𝑎 𝑔 T_{P,mag}

Coefficients	$T_{P,mag}$
$a$	0.82
$b$	1.44
$c$	1.07
$d$	5.48
$e$	0.223

Equations

$$C_{D,\varphi }=C_{D,\varphi =0^{\circ }}+(C_{D,\varphi =90^{\circ }}-C_{D,\varphi =0^{\circ }}){\sin }^2\varphi .$$

$$C_{L,\varphi }=(C_{D,\varphi =90^{\circ }}-C_{D,\varphi =0^{\circ }})\sin \varphi \cos \varphi$$

$$|f(\mathbf{r}+{\mathbf{e}}_{\alpha }\Delta t,t+\Delta t)⟩=|f(\mathbf{r},t)⟩-\text{mathsfbi}{M}^{-1}\text{mathsfbi}\hat{S}(|m(\mathbf{r},t)⟩-|\ast ⟩m^{(eq)}(\mathbf{r},t)),$$

$${\mathbf{g}}_{new}={\mathbf{g}}_{prev}+\frac{({\mathbf{u}}_{s,ref}-{\mathbf{u}}_{s,prev})}{K_p^2}\Delta t,$$

$$\text{mathsfbi}S=⟨\mathbf{p}{\mathbf{p}}^T⟩.$$

$${\mathbf{p}}_{new}=\{\begin{array}{ll}{\mathbf{p}}_{new},& \text{if }dE<0\\ {\mathbf{p}}_{new},& \text{if }dE\ge 0\text{ and }U([0,1])<\exp (-dE)\\ {\mathbf{p}}_{curr},& \text{otherwise}\end{array}$$

$$dE=A(({\mathbf{p}}_{new}\cdot \mathbf{n}{)}^2-({\mathbf{p}}_{curr}\cdot \mathbf{n}{)}^2).$$

$${\mathbf{F}}_{f\to p}=\mathbf{F}+{\mathbf{F}}_b.$$

$$\mathbf{F}={\mathbf{F}}_{f\to p}(1-{ϵ}_s).$$

$${\mathbf{F}}_{norm}$$

$${\mathbf{T}}_{norm}$$

$${\hat{\mathbf{e}}}_1$$

$${\hat{\mathbf{e}}}_2$$

$${\hat{\mathbf{e}}}_3$$

$${\sigma }_{\bar{x}}=\sigma /\sqrt{n}.$$

$$F_D=$$

$$F_2$$

$$F_L=$$

$$T_1$$

$$T_P=$$

$$T_3$$

$$F_{D,\varphi }=F_{D,\varphi =0^{\circ }}+(F_{D,\varphi =90^{\circ }}-F_{D,\varphi =0^{\circ }}){\sin }^2\varphi .$$

$$F_D(\text{Rey},{ϵ}_s)$$

$$F_{d,isol}(\text{Rey})$$

$$F_{{ϵ}_s}({ϵ}_s)$$

$$F_{\text{Rey},{ϵ}_s}(\text{Rey},{ϵ}_s)$$

$$F_{l,\varphi }(\text{Rey},\varphi )$$

$$F_{l,isol}(\text{Rey})$$

$$S_{f,\varphi }(\text{Rey},\varphi )$$

$$F_{L,\varphi }(\text{Rey},{ϵ}_s,\varphi )$$

$$F_{L,mag}(\text{Rey},{ϵ}_s)$$

$$F_{{ϵ}_s}({ϵ}_s)$$

$$F_{\text{Rey},{ϵ}_s}(\text{Rey},{ϵ}_s)$$

$$F_{L,\varphi }=(F_{D,\varphi =90^{\circ }}-F_{D,\varphi =0^{\circ }})\sin \varphi \cos \varphi .$$

$$F_{L,\varphi }=F_{L,simple}\sin \varphi \cos \varphi .$$

$$T_{p,\varphi }(\text{Rey},\varphi )$$

$$T_{p,isol}(\text{Rey})$$

Full text

1 Introduction

Accurate fluid-particle drag, lift and torque closures are required for precise Euler-Lagrangian simulations of non-spherical particles. Historically, different drag closures have been developed for assemblies of spherical particles (Beetstra et al., 2007; Tenneti et al., 2011; Tang et al., 2015). However, practical flows often involve assemblies of non-spherical particles for which there exist no closures at the moment. Even for static, mono-disperse, non-spherical particle assemblies, creating the required closures is complicated due to the different possible mutual orientations of the particles. Furthermore, there is a lack of knowledge identifying the relevant parameters that can parametrize the drag, lift and torque, which adds to the complication. Most fluidization applications involve gas-solid flows, in which case the large density ratios ensure large Stokes numbers, i.e. the typical relaxation time of the solid particle velocity is large relative to the response time of the gas (Sanjeevi et al., 2018a). It has been shown that under such conditions, it is sufficient to assume the particle configurations to be quasi-static (Rubinstein et al., 2017).

Conventionally, fluidization simulations of non-spherical particles are performed by combining isolated particle drag correlations with correlations expressing the voidage effects as determined for sphere assemblies. There have been several works in the past focussing on the drag experienced by isolated non-spherical particles. Hölzer & Sommerfeld (2008) proposed a correlation for the drag coefficient $C_{D}$ for non-spherical particles. The proposed correlation is a function of particle sphericity and crosswise-sphericity, based on the projected area, which indirectly represents the particle orientation. Their proposed correlation is based on literature data of different non-spherical particles of various shapes and aspect ratios. More recently, drag, lift and torque closures for isolated non-spherical particles have been derived based on direct numerical simulations. Zastawny et al. (2012) developed drag, lift and torque coefficients for four different non-spherical particles as a function of Reynolds number $\Rey$ and incident angle $\phi$ with respect to the incoming flow. The investigated particles have aspect ratios ranging from 1.25 to 5 and $\Rey\leq 300$. Similarly, Richter & Nikrityuk (2013) proposed fits for drag, lift, torque coefficients for cubic and ellipsoidal particles. The above mentioned literature is primarily limited to steady flow conditions. Recently, we developed drag, lift and torque closures for three different non-spherical particles from the viscous Stokes regime upto the high $\Rey$ regime of $\Rey=2000$, involving complex, unsteady flows (Sanjeevi et al., 2018a). In an earlier work (Sanjeevi & Padding, 2017), we reported the interesting finding that the drag coefficient $C_{D}$ at different incident angles $\phi$ follows a sine-squared scaling given by

[TABLE]

Likewise, we reported another interesting finding that the lift coefficient $C_{L}$ follows sine-cosine scaling at different $\phi$ as

[TABLE]

for various elongated particles. The above mentioned scaling laws must be mathematically true in the Stokes regime due to linearity of the flow fields. However, their validity in the inertial regimes is primarily due to an interesting pattern of pressure distribution contributing to the drag and lift for different incident angles (Sanjeevi & Padding, 2017). In equations 1 and 2, the drag coefficients at incident angles of 0 and 90 degrees still depend on particle shape and Reynolds number. The Reynolds number in the present work is defined as $\Rey=|\bm{u}_{s}|d_{eq}/\nu$, where $\bm{u}_{s}$ is the superficial flow velocity, $\nu$ is the kinematic viscosity of the fluid, and $d_{eq}$ is the diameter of the volume-equivalent sphere given by $d_{eq}=(6V_{p}/\pi)^{1/3}$ with $V_{p}$ the particle volume.

For multiparticle systems, various literature is available to include the voidage effects, often developed through experiments and numerical simulations. One of the most widely used expressions is that of Ergun (1952), which has been developed based on a series of packed bed experiments of different particle shapes. The only limitation of this work is that it is applicable primarily in the dense limit. Richardson & Zaki (1954) performed various sedimentation and fluidization experiments and proposed accordingly the effect of particle volume fraction on the drag. Based on the previous literature on sedimentation and packed bed experiments, Di Felice (1994) bridged the dilute and dense particulate regimes through a unified function, which also extends from low to high $\Rey$. Though the above correlations provide a good approximation, the use of such closures in Euler-Lagrangian simulations often do not represent accurate physics. This is mainly due to the inability to construct moderate solids volume fractions in experiments.

There is a growing interest to use numerical simulations to accurately develop drag closures for different Reynolds numbers $\Rey$ and solids volume fractions $\epsilon_{s}$, albeit primarily for spheres. Initially, lattice Boltzmann method (LBM) has been the choice for simulating assemblies of spheres (Hill et al., 2001; van der Hoef et al., 2005; Beetstra et al., 2007). Recently, Tenneti et al. (2011) used an immersed boundary method (IBM) to develop drag closures for static assemblies of spheres for $0.01\leq\Rey\leq 300$ and $0.1\leq\epsilon_{s}\leq 0.5$. They observed a deviation of 30% in the $\Rey$ range from 100 to 300 with respect to the earlier work of Beetstra et al. (2007). This is possible because Beetstra et al. (2007) used LBM with the conventional stair-case boundary condition to represent the sphere boundaries, for which at high $\Rey$ thinner boundary layers result in larger deviations. In this work, we use a multi-relaxation time (MRT) LBM for high $\Rey$ flows and an interpolated bounceback scheme to much more accurately represent the particle geometry. Recently, Tang et al. (2015) used an IBM based solver to create drag closures for static assemblies of spheres upto $\Rey\leq 1000$ and $\epsilon_{s}\leq 0.6$.

There are several disadvantages with combining an isolated non-spherical particle drag with a voidage function based on spheres. First, the assumption that the voidage effects are independent of particle shape is probably incorrect, since there exist different closures even for assemblies of polydisperse spheres (Beetstra et al., 2007; Holloway et al., 2010). Second, the voidage effects on lift and torque in a multiparticle system are unknown and hence are often neglected in Euler-Lagrangian simulations (Oschmann et al., 2014; Mahajan et al., 2018). Thirdly, using the same factor for voidage effects for all incident angles $\phi$ may hold in sufficiently dilute regimes but its validity in the dense limit is unknown. At the moment, only He & Tafti (2018) have discussed the drag, lift and torque for an assembly of non-spherical particles. However, they do not propose any correlations which can be used in Euler-Lagrangian simulations. This could be due to the difficulty in identifying the dependent parameters which represent the orientation effects in non-spherical, multiparticle system adequately.

In this work, we propose and subsequently identify the important dependent parameters for static, mono-disperse assemblies of axisymmetric non-spherical particles. With the identified parameters, we create the drag, lift and torque closures accordingly. Our particle of interest is a capsule-like spherocylinder of aspect ratio 4 (total length/shaft diameter). Compared to the two parameters for sphere assemblies, i.e. Reynolds number $\Rey$ and solids volume fraction $\epsilon_{s}$, we propose four additional parameters for the assembly of axisymmetric non-spherical particles. Two parameters describe the mutual orientations of the particles, namely two eigenvalues $S_{1}$ and $S_{2}$ of the orientation tensor, and two angle parameters $\alpha$ and $\beta$ represent the polar and azimuthal angles of the average flow (in the coordinate frame determined by the principal directions of the order tensor). The resulting six dimensional parameter space is adequately explored and correlations are proposed accordingly. It should be noted that the fixed nature of the particles in our simulations imply that the proposed correlations are applicable for high Stokes number flows as typically experienced by Geldart D category particles.

2 Numerical method

2.1 Lattice Boltzmann method

In the present work, we use a D3Q19, multi-relaxation time (MRT) lattice Boltzmann method (d’Humières et al., 2002) to simulate the fluid flow. The numerical method is adequately explained and validated in our previous works (Sanjeevi & Padding, 2017; Sanjeevi et al., 2018a). The evolution of particle distribution function $\ket{f}$ is computed as

[TABLE]

for position $\bm{r}$ with discrete velocities $\bm{e}_{\alpha}$ in directions $\alpha=1,2...,19$. Equation 3 is solved in a sequence of two steps namely collision and streaming. $\mathsfbi{M}$ is a $19\times 19$ transformation matrix used to transform $\ket{f}$ from velocity space to moment space $\ket{m}$ with $\ket{m}=\mathsfbi{M}\cdot\ket{f}$. Here, the ket vector $\ket{\cdot}$ implies a column vector. The relaxation matrix $\mathsfbi{\hat{S}}=\mathsfbi{M}\cdot\mathsfbi{S}\cdot\mathsfbi{M}^{-1}$ is a $19\times 19$ diagonal matrix. $\mathsfbi{\hat{S}}$ utilizes different, optimally chosen relaxation rates for different moments, thereby providing better stability compared to the single-relaxation-time LBM scheme (d’Humières et al., 2002). The matrices $\mathsfbi{M}$ and $\mathsfbi{\hat{S}}$ are similar to Huang et al. (2012) and are given in Sanjeevi et al. (2018a). The density is computed as $\rho=\sum_{\alpha}f_{\alpha}$ and the momentum as $\rho\bm{u}=\sum_{\alpha}f_{\alpha}e_{\alpha}$. The relation between the kinematic viscosity of the fluid and the dimensionless relaxation time $\tau$ is $\nu=c_{s}^{2}(\tau-1/2)\Delta t$, and the pressure $p$ is related to the density by $p=\rho c_{s}^{2}$, where $c_{s}$ is the speed of sound. A linearly interpolated bounce back scheme (Bouzidi et al., 2001; Lallemand & Luo, 2003) is used to accurately consider the curved geometry of the particle, as opposed to the traditional stair-case bounce back boundary condition. The flow is driven by a body force $\bm{g}$ and the simulated domain is periodic in all three directions. The use of the interpolated bounce back scheme within a periodic domain results in a slow mass leakage/gain in the system. Accordingly, the mass is corrected using a case 3 type correction described in Sanjeevi et al. (2018b). The results for the multiparticle system are validated in section 3.4.

The ratio of $d_{eq}/d_{min}$ equals 1.765 for the considered spherocylinder of aspect ratio 4, where $d_{min}$ implies diameter of the cylinder. The simulation parameters used in our LBM simulations are summarized in table 1. Specifically, it can observed that a good particle resolution ($d_{eq}$) is maintained for different $\Rey$. Further with increasing $\epsilon_{s}$, the $d_{eq}$ is increased accordingly to resolve increased velocity gradients at high $\epsilon_{s}$. All LBM simulations have cubic domain, each with 200 particles unless otherwise specified. At least two independent simulations are performed for each $\Rey$ and $\epsilon_{s}$ and the details of independent number of simulations are discussed later (see figure 13).

2.2 Flow control

In order to perform a simulation for a specific $\Rey$, it is required to control the superficial flow velocity $\bm{u}_{s}$ by applying a body force $\bm{g}$. The relationship between the superficial velocity and the average interstitial flow velocity $\bm{u}_{avg}$ is given by $\bm{u}_{s}=(1-\epsilon_{s})\bm{u}_{avg}$. Due to the non-spherical nature of the particles, the sum of lift forces is often non-zero, and the resultant direction of $\bm{u}_{s}$ can be different from the direction of $\bm{g}$. This necessitates the need to control both direction and magnitude of the body force. Initially, the fluid is at rest with both $\bm{u}_{s}$ and $\bm{g}$ zero. The flow is slowly ramped up by increasing $\bm{g}$ until the desired $\bm{u}_{s}$ is achieved. For each timestep, the updated gravity $\bm{g}_{new}$ is computed as

[TABLE]

where $\bm{g}_{prev}$ is the gravity from the previous timestep, $\bm{u}_{s,ref}$ is the desired reference superficial velocity, and $\bm{u}_{s,prev}$ is the superficial velocity from the previous timestep. $K_{p}$ is a time constant which controls the system response rate. The stopping criterion for the simulations is when the system $\bm{u}_{s}$ reaches 99.9% of the reference setpoint.

3 Simulation setup

3.1 Orientation tensors

In this section, we briefly explain the characterization of mutual orientations in an assembly of axisymmetric non-spherical particles with orientation tensors. We subsequently explain the use of a Maier-Saupe potential to achieve the desired particle configurations through Monte-Carlo simulations.

To describe the orientation of a single axisymmetric particle, the azimuthal and polar angles are sufficient. For a multiparticle configuration, it is important to parametrize the mutual orientations of the particles, with the least number of parameters. For this, we propose to use the orientation tensor $\mathsfbi{S}$ which is defined as the average of the dyadic products of the particle orientation vectors. In other words,

[TABLE]

Here, $\bm{p}$ is the unit orientation vector of a particle. The 3 eigenvalues (which we order as $S_{1},S_{2},S_{3}$ from small to large) characterize the type of mutual alignment, as shown in figure 1. The corresponding 3 eigenvectors define the principal directions of mutual particle alignment.

Because the trace of S is 1, only 2 eigenvalues are sufficient to specify the amount of randomness, planar random (bi-axial), or unidirectional (nematic) order. It should be noted that the tensor $\mathsfbi{S}$ is insensitive to an orientation $\bm{p}$ or $-\bm{p}$ of particles. In other words, the tensor captures essentially the mutual alignment of particles irrespective of particles oriented in positive or negative direction. Figure 1($a$) shows a completely random configuration with $S_{1}=S_{2}=S_{3}=1/3$. Figure 1($b$) shows a planar random configuration with particles primarily confined to planes (in this example with random orientations in planes normal to the $x$-direction) resulting in $S_{1}=0,S_{2}=S_{3}=1/2$, and similarly a unidirectional (nematic, in this example in the $z$-direction) configuration in figure 1($c$) with $S_{1}=S_{2}=0,S_{3}=1$. In practical conditions, particles can exhibit more complex configurations in between these extremes but can be adequately described by 2 eigenvalues $S_{1}$ and $S_{2}$. Regarding the unidirectional case, we consider only nematic configurations but not smectic because ordering of both positions and orientations is rare in fluidization conditions.

The above metrics can be used to describe the particle configuration. However, due to the nonsphericity of the particles, the flow orientation with respect to the principal directions of the particle orientations is also important. This results in two parameters, namely the polar angle ($\alpha$) and azimuthal angle ($\beta$) of the average flow velocity vector with respect to the space spanned by the 3 eigenvectors of the orientation tensor. In summary, the parameter space to be explored for our flow problem has 6 parameters, namely Reynolds number $\Rey$, solids volume fraction $\epsilon_{s}$, two particle configuration parameters $S_{1}$, $S_{2}$ and two angles $\alpha$ and $\beta$ describing the mean flow orientation with respect to the configuration.

3.2 Generation of biased particle configurations

The generation of non-overlapping configurations of the particles in a periodic domain is required as an input for the flow simulations. Further, it is also required to generate configurations of particles with a prescribed orientation tensor, which adds further complexity. In this section, we briefly describe the Monte-Carlo simulation algorithm for generating configuration of non-overlapping particles and the use of a Maier-Saupe potential (Maier & Saupe, 1959) to bias the system to produce the required orientation tensor.

As the particles are spherocylindrical in shape, a simple way to detect overlap is to find the minimum distance between two line segments. We define the line segment as the line connecting the centres of the two spheres at the extremes of the spherocylinder. If the distance between two line segments is less than the particle diameter, then the spherocylinders overlap. A fast algorithm is used to measure the shortest distance between the line segments (Vega & Lago, 1994).

Using the above overlap detection algorithm, the particles are randomly translated in small steps compared to the particle size and rotated by a small angle around a randomly chosen axis. This procedure results in a random configuration after many iterations. If a prescribed orientation tensor is required, besides the requirement of no overlap, the following Maier-Saupe potential is applied to accept or reject a new orientation of a particle. We define a director $\bm{n}$ along which the system is biased towards or against. Inside each Monte-Carlo simulation step, a new particle orientation $\bm{p}_{new}$ can be accepted or rejected from the current orientation $\bm{p}_{curr}$ based on following criteria:

[TABLE]

Here, $dE$ is the increase in Maier-Saupe potential and $U([0,1])$ is a random number uniformly distributed between 0 and 1. The mutual particle orientations emerge from the balance between the random rotations, which tend to disorder the particle orientations, and the Maier-Saupe potential, which tend to order the particle orientations. The magnitude of $A$ determines the intensity of the configuration towards the director. A planar random configuration is achieved with the plane perpendicular to the director $\bm{n}$, if $A$ is positive. A unidirectional configuration along the direction of $\bm{n}$ is achieved, if $A$ is negative. Higher absolute $A$ values result in better perfection towards the desired configuration. With the mentioned strategy, any configuration in-between the ideal cases shown in figure 1 can be achieved. Some sample configurations generated using the above mentioned algorithm are shown in figure 2. For simplicity, the eigenvectors of the orientation tensor $\mathsfbi{S}$ are considered as aligned with the Cartesian coordinate system in figure 2. The shown configurations are respectively equivalent to figure 1. For better clarity, the shown configuration has only 50 particles and the solids volume fraction $\epsilon_{s}$ is 0.1. The actual flow simulations have 200 particles and are performed for various $\epsilon_{s}$.

A common intuition may be that a random configuration would result in particles with evenly distributed values of the incident angle $\phi$. However for a random configuration, the available number of particles at different $\phi$ are not uniform, as shown in figure 3($a$). This is due to the higher probability to find particles at an angle $\phi$ near $90^{\circ}$ because the Jacobian for a spherical coordinate system scales as $\sin\phi$. Therefore, the disadvantage for a random configuration is that there are actually few data points at $\phi=0^{\circ}$ to create angle-dependent closures. On the contrary, the planar configuration with the planes parallel to the flow direction results in even particle distributions, as shown in figure 3($b$). This information is considered while we generate configurations for the flow simulations.

3.3 Forces and torques acting on a particle

For an assembly of particles, different definitions are used to report the forces (Beetstra et al., 2007; Tenneti et al., 2011; Tang et al., 2015). To ensure consistency, it is important to know the form of the reported results. For a packed bed of particles in a flow induced by a macroscopic pressure gradient $\nabla P$, each particle of volume $V_{p}$ experiences a resulting force $\bm{F}$ due to the flow and a buoyancy force $\bm{F}_{b}=-V_{p}\nabla P$ due to the pressure gradient. For such a case, the total fluid-to-particle force $\bm{F}_{f\rightarrow p}$ acting on a particle is

[TABLE]

Given $N$ particles with each of volume $V_{p}$ and total volume of the system $V$, the solids volume fraction is given by $\epsilon_{s}=NV_{p}/V$. Further, the relationship between $\bm{F}$ and $\bm{F}_{f\rightarrow p}$ is given by (Tang et al., 2015)

[TABLE]

In this work, we report the forces $\bm{F}$ due to the flow and not $\bm{F}_{f\rightarrow p}$. The effects of buoyancy on torques are unknown and hence the reported torques $\bm{T}$ are also as they are determined from the simulations. We normalize the force and torque with the Stokes drag and torque of a volume-equivalent sphere:

[TABLE]

Here, $\mu$ is the dynamic viscosity and $R_{eq}$ is the radius of the volume equivalent sphere. Let $\bm{p}$ be the normalized orientation vector of the considered particle. The local coordinate system for each particle is defined as

[TABLE]

The above defined axes are accordingly illustrated in figure 4. The incident angle $\phi$ a particle makes with respect to the incoming flow is given by $\phi=\cos^{-1}(|\hat{\bm{e}}_{1}\cdot\bm{p}|)$. We also compute the average forces and torques for different $\phi$ intervals. Due to the finite number of measurements in these intervals, there is an error on the mean $\bar{x}$ of any property $x$. We use the standard error on the mean $\sigma_{\bar{x}}$ for the errorbars, computed as

[TABLE]

Here $\sigma$ is the standard deviation of the corresponding variable $x$ and $n$ is the number of data points within the given $\phi$ interval. The normalized drag $F_{D}$ and lift $F_{L}$ can be computed from $\bm{F}_{norm}$ as

[TABLE]

Since the reported forces are without buoyancy effects, the $(1-\epsilon_{s})$ term must be considered accordingly for both drag and lift while performing Euler-Lagrangian simulations. Due to the influence of neighbouring particles, the lateral force $F_{2}$ for each individual particle may not be equal zero, as shown in figure 5 ($\Rey=100$ and $\epsilon_{s}=0.3$). However, due to symmetry, the average $F_{2}$ does equal zero. Therefore, $F_{2}$ is not considered in our further discussion. The torques about the above defined axes are

[TABLE]

Here $T_{P}$ is the pitching torque acting on a particle. We show the three different torques for a flow through a random particle configuration at $\Rey=100$ and $\epsilon_{s}=0.3$ in figure 6. It can be observed that $T_{1}$ and $T_{3}$, though having some non-zero values, are statistically zero on average due to symmetry. The non-zero values are primarily due to hydrodynamic interactions with other particles. Only the average pitching torque $T_{P}$ (or $T_{2}$) remains non-zero for different $\phi$ and varies as $\sin\phi\cos\phi$. Though individual particles experience non-zero $T_{1}$ and $T_{3}$, they become zero at $\phi=0^{\circ}$ and $\phi=90^{\circ}$ respectively, where the axis of symmetry of the particle coincides with the measured axis for torque. This implies that the hydrodynamic interaction of particles does not induce a torque (or a spin) about the axis of symmetry of the particle.

3.4 Validation

Sufficient validation has been done for our LBM code in the past for flow around isolated particles (Sanjeevi & Padding, 2017; Sanjeevi et al., 2018a). For a multiparticle configuration, we have chosen flow around a random assembly of 100 particles at $\Rey=100$ and $\epsilon_{s}=0.3$ and measure the $F_{D}$ experienced by the individual particles. The LBM results are compared with results from COMSOL Multiphysics, a body-fitted, unstructured mesh based FEM solver. The simulated LBM domain is of size $360^{3}$. The volume equivalent sphere diameter is $d_{eq}=64.4$ lattice cells. The superficial velocity $u_{s}$ is 0.0414 and the kinematic viscosity $\nu$ is 0.08/3 in lattice units. The FEM solver domain is made of 2.1 million elements. The resulting drag forces are shown in figure 7. A good agreement between LBM and FEM results can be observed. The average $F_{D}$ experienced by all particles in LBM and FEM solvers are 26.6 and 26.4 respectively. Also a good match in $F_{D}$ values for individual particles at different $\phi$ can be observed.

4 Tests of configuration independence

Given a six-dimensional parameter space, exploring each dimension with approximately 5 simulations, results in a massive $5^{6}=15625$ simulations. Furthermore, closures must be created for drag, lift and torque as a function of this six-dimensional space. Before proceeding with these simulations, we tried to identify if there are any independent parameters specifically related to the mutual orientation of particles. In this section, we will show that the flow around a non-spherical particle assembly is independent of the mutual orientation of the particles themselves. This configuration independence removes the configuration parameters $S_{1}$, $S_{2}$ and flow angle parameters $\alpha$ and $\beta$ from the parameter space to be explored. We find that, when averaged over a number of particles, the only dependence that the particles exhibit regarding orientation is the particle’s incident angle $\phi$ as in flow around single particles. Effectively, we will show that the flow problem depends only on the Reynolds number $\Rey$, solids volume fraction $\epsilon_{s}$ and the incident angle $\phi$ of individual particles with respect to the flow direction.

In the extremely dilute regimes, i.e. $\epsilon_{s}\rightarrow 0$, it is already shown that there exists a sine-squared scaling of drag for elongated non-spherical particles (Sanjeevi & Padding, 2017; Sanjeevi et al., 2018a). In this section, we discuss the results of flow around different configurations at an intermediate solids volume fraction of $\epsilon_{s}=0.3$. Results of different configurations (in respective plot insets) at an intermediate $\Rey=100$ are shown in figure 8 such as fully random, planar random with flows parallel and perpendicular to the planes and unidirectional configurations with principal directors at different angles. Though there exists scatter in the measured $F_{D}$, it can be observed that the average $F_{D}$ at different $\phi$ interval scales similar to sine-squared scaling as in our earlier works of isolated particles. In other words, the $F_{D}$ at any $\phi$ can be computed as

[TABLE]

It is important to note that the same values for $F_{D,\phi=0^{\circ}}$ and $F_{D,\phi=90^{\circ}}$ emerge for all configurations. Likewise, we also show that the scaling phenomenon extends to both Stokes and high $\Rey$ regimes in figure 9. With the sine-squared scaling behaviour (or the configuration independence) identified at $\epsilon_{s}=0$ and $\epsilon_{s}=0.3$, it can be inferred that the scaling is safely applicable in the region $0\leq\epsilon_{s}\leq 0.3$. We have verified the same at $\epsilon_{s}=0.1$ and the results are not shown here for brevity. Though we observe the results are dependent on only 3 parameters, namely $\Rey$, $\epsilon_{s}$ and $\phi$, the simulation needs to be set up for only two parameters, namely $\Rey$ and $\epsilon_{s}$. With a sufficiently random configuration, the system involves different particle orientations covering all $\phi$. A caveat with a random configuration is that there are always very few particles near $\phi=0^{\circ}$, as shown in section 3.2. Therefore, biased random configurations with more particles at $\phi=0^{\circ}$ are created and at least two simulations are performed for better statistics.

We also observe the configuration independence phenomenon at $\epsilon_{s}=0.4$. The criterion considered to declare configuration independence phenomena is that the average drag results in a given $\phi$ range of different configurations are within 10% deviation. In almost all cases, the deviations are within $\pm$ 5%. However in a dense case with $\epsilon_{s}=0.5$, several more factors such as the mutual orientations, relative positions of particles, etc. influence the results. The $F_{D}$ distribution for such dense configurations at $\Rey=100$ and $\epsilon_{s}=0.5$ are given in figure 10. Although these results can be predominantly parametrized by $\Rey$, $\epsilon_{s}$, and $\phi$, the influence of the additional parameters cannot be ignored. Therefore, specific cases of $\epsilon_{s}=0.5$ are performed with more simulations for better statistics.

For a practical fluidization or other relevant gas-solid flow simulation, the densest configuration is most likely to occur when the particles are at bottom or at rest (e.g. before the start of fluidization). In such a dense condition, the particle configuration itself is dependent on the wall geometry. For a typical bed configuration with a flat wall at the bottom, the particles also roughly align in planes parallel to the bottom wall, i.e. a planar random configuration. Pournin et al. (2005) observed the same for particles poured freely from the top. Similarly, we also observe the same for a bed containing freely poured particles settled under gravity ($\epsilon_{s}=0.54$), as shown in figure 11. The bed contains 30000 particles and it can be observed that roughly 2/3 of all particles are in the range $\phi=70-90^{\circ}$ confirming our hypothesis. Given such criteria, the most relevant regime would be to generate an accurate fit for $F_{D,\phi=90^{\circ}}$ at high $\epsilon_{s}$.

It should also be noted that with increasing aspect ratio of elongated particles, the maximum $\epsilon_{s}$ decreases for a packed bed (Williams & Philipse, 2003). This is because the locking phenomenon is stronger with high aspect ratio particles. Unless the particles are packed with their orientations aligned, the decrease in peak $\epsilon_{s}$ for high aspect ratio elongated particles is unavoidable. Also, practical applications as shown in figure 11 do not allow such long range ordering. A decreasing peak $\epsilon_{s}$ implies that the configuration independence phenomenon will be very applicable. With the observed sine-squared drag scaling, the pressure drop across a packed bed can be determined with the knowledge of the $\phi$ distribution alone.

In the subsequent sections, we will show that in the dilute and intermediate $\epsilon_{s}$ regimes, the influence of $\epsilon_{s}$ is nearly shape independent. This implies that the drag on isolated non-spherical particles can be combined with sphere-based multiparticle correlations for the voidage effect to mimic flow around assemblies of non-spherical particles upto intermediate $\epsilon_{s}$.

4.1 Explored regimes

In this section, we briefly explain the regimes explored in the current work and also explain the number of independent simulations performed per regime tested. An example of the flow stream lines for a random configuration at $\Rey=100$ and $\epsilon_{s}=0.3$ is shown in figure 12. Until solids volume fractions of $\epsilon_{s}=0.35$, the generation of randomly orientation configurations is possible, as experienced by He & Tafti (2018) for prolate spheroids of aspect ratio 2.5. In our case, we are able to achieve random configurations upto $\epsilon_{s}=0.4$. However for denser configurations, it is difficult to generate a truly random configuration. For dense configurations of $\epsilon_{s}=0.5$, the particles have a natural tendency to orient to planar random or unidirectional orientation configurations. A truly random configuration with a finite number of particles, at such solids volume fraction, is not possible. This is due to a strong orientation bias imposed by neighbouring particles due to lack of inter-particle space. The explored regimes are indicated in figure 13. Overall, at least two simulations are performed for the explored regimes. However for specific cases of dilute and intermediate $\epsilon_{s}$, we performed 5 simulations with 2 random, 1 planar random with flow aligned to the plane and 2 unidirectional configurations with flow parallel and perpendicular to the principal director. For solids fraction $\epsilon_{s}=0.5$, 3 planar random configurations with flows aligned to the plane and 2 unidirectional configurations with flows parallel and perpendicular to the principal director are performed. For cases with more simulations, the results are accordingly weighted while making the fits.

5 Results

5.1 Drag

With sine-squared scaling valid as shown in section 4, the drag experienced by a particle in a multiparticle system can be explained by the equation 22 involving only the drag experienced at $\phi=0^{\circ}$ and $\phi=90^{\circ}$. Therefore, we propose to generate fits for $F_{D,\phi=0^{\circ}}$ and $F_{D,\phi=90^{\circ}}$ as a function of $\Rey$ and $\epsilon_{s}$ as

[TABLE]

The corresponding terms are as follows:

[TABLE]

Here, $C_{d,isol}$ is the isolated particle drag at given $\Rey$ as detailed in Sanjeevi et al. (2018a) for the considered particle (fibre) for both $\phi=0^{\circ}$ and $\phi=90^{\circ}$. The coefficients in equations 25 and 26 for both $F_{D,\phi=0^{\circ}}$ and $F_{D,\phi=90^{\circ}}$ are given in table 2. The average absolute deviation of the fits and simulation data are 3.5% and 2% for $F_{D,\phi=0^{\circ}}$ and $F_{D,\phi=90^{\circ}}$, respectively.

The simulated data and corresponding fits are shown in figure 14. The fits follow the physical limits beyond the $\Rey$ range simulated as shown in figure 15. In the Stokes flow limit, it can be observed that both $\phi=0^{\circ}$ and $\phi=90^{\circ}$ normalized drag becomes independent of $\Rey$. In the high $\Rey$ limit, the normalized drag approaches a linear dependency on $\Rey$.

The ratio of the perpendicular to parallel drag $F_{D,\phi=90^{\circ}}/F_{D,\phi=0^{\circ}}$ at different $\Rey$ and $\epsilon_{s}$ is shown in figure 16. For low $\Rey$ ($\Rey=0.1$), the ratio remains constant at a value a little larger than 1 for all $\epsilon_{s}$. The reason for this is that at low $\Rey$, the particles experience stronger viscous effects. The viscous drag reduces and pressure drag increases with increasing $\phi$ at low $\Rey$. The same has been confirmed for isolated particles (Sanjeevi & Padding, 2017) and for a multiparticle system (He & Tafti, 2018). The combined viscous and pressure drag components result in a drag ratio close to 1 for the considered spherocylinders at low $\Rey$. Due to inertial dominance at moderate and large $\Rey$ ($\Rey\geq 100$) we can observe a near constant drag ratio for solids volume fractions upto $\epsilon_{s}=0.3$ and a decrease in the ratio for $\epsilon_{s}>0.3$. Further, figure 16 gives an indication that for very dense crowding, i.e. at $\epsilon_{s}>0.5$, there is a possibility that $F_{D,\phi=90^{\circ}}/F_{D,\phi=0^{\circ}}$ tends back to approximately 1. Up to moderate crowding, although the flow is disturbed due to the presence of neighbouring particles, there is sufficient inter-particle space for flow to achieve uniformity. However with increased particle crowding, there appear pronounced fluctuations in flow velocities (see also section 5.4), resulting in a reduced drag ratio at high $\epsilon_{s}$. This is an important finding because the traditional approach of Euler-Lagrangian simulations involve combining isolated non-spherical particle drag with the voidage effects based on sphere assemblies. This would result in a constant drag ratio for all $\epsilon_{s}$. This in turn could affect Euler-Lagrangian simulation results, especially in predicting the minimum fluidization velocity as there exists a dense packing of particles. This mandates the need for the current work.

Figure 17 shows a similar interesting observation: The scaling of the voidage effect $F_{D}(\epsilon_{s})/F_{D}(\epsilon_{s}=0)$ in the inertial regime (high $\Rey$ limit) is shape and orientation independent for $\epsilon_{s}\leq 0.3$. Here, we have normalized the drag with respective isolated particle drag for different $\Rey$ and $\phi$. It can be observed that all the normalized points fall on a same trend until $\epsilon_{s}=0.3$. Similar normalized $F_{D}$ for spheres from Tang et al. (2015) at $\Rey=100$ and $\Rey=1000$ also show the same trend until $\epsilon_{s}=0.3$. Here, we use the isolated sphere drag correlation of Schiller & Naumann (1935) for the normalization. The predictions of Tenneti et al. (2011) for spheres do not follow the exact trend for the voidage effects as observed from figure 17. It should be noted that Tenneti et al. (2011) explored only until $\Rey=300$ in their work and extrapolation to high $\Rey$ may not apply. Therefore, the above discussion indicates that spherical drag correlations for the voidage effect, combined with isolated non-spherical particle drag correlations can be applied to dilute suspension simulations of non-spherical particles in the inertial regimes. For a given non-spherical particle, the effect of crowding ($\epsilon_{s}$) on $F_{D}$ is different for different $Re$ and $\phi$. Figure 18 shows the voidage effect ($F_{D}$ normalized by the corresponding isolated particle drag) as a function of $\Rey$. It can be seen at low $\Rey$, the increase in drag due to crowding is comparable for both $\phi=0^{\circ}$ and $\phi=90^{\circ}$ at different $\epsilon_{s}$. At high $\Rey$, the increase in drag due to crowding with increasing $\epsilon_{s}$ is much stronger for $\phi=0^{\circ}$ compared to $\phi=90^{\circ}$. This also explains further the reason for the observed reduction in perpendicular to parallel drag ratios with increasing $\epsilon_{s}$ in figure 16.

In the previous sections, we discussed the $F_{D}$ averaged over all particles with similar $\phi$. However, the distribution of $F_{D}$ within a $\phi$ interval is itself also a function of both $\Rey$ and $\epsilon_{s}$. The standard deviations of the distribution of drag measurements, normalized by the average $F_{D}$ in the corresponding interval, are plotted in figure 19. It is important that the standard deviations are normalized by the respective average $F_{D}$, rather than against a single value, say $F_{D,\phi}=90^{\circ}$, for a given $\Rey$ and $\epsilon_{s}$. This is because with increasing $\Rey$, the ratio $F_{D,\phi=90^{\circ}}/F_{D,\phi=0^{\circ}}$ increases, as shown in figure 16 earlier. Therefore, using $F_{D,\phi=90^{\circ}}$ for normalization will make the standard deviations at $\phi=0^{\circ}$ appear insignificant at large $\Rey$.

For dilute configurations ($\epsilon_{s}=0.1$), we clearly observe that increasing $\Rey$ results in an increased $\sigma_{F_{D}}/F_{D}$ at all $\phi$. It should be noted that the absolute magnitudes of $F_{D}$ at $\Rey=1000$ are much larger than at $\Rey=0.1$. Despite the normalization by these larger values, we observe increased standard deviations for higher $\Rey$. This is because at low $\Rey$, the viscous effects dominate, resulting in long-range flow uniformity. Conversely, at high $\Rey$, the boundary layers are thinner and flow wakes are stronger. This results in high non-uniformity in the incoming flow on each particle, and thereby large fluctuations in the hydrodynamic forces. For dense particle configurations ($\epsilon_{s}=0.5$), it can be observed that $\sigma_{F_{D}}/F_{D}$ increases relative to dilute conditions, with a higher standard deviation for higher $\Rey$. The reason for higher spread in $F_{D}$ is due to the fact the particles locally encounter highly non-uniform incoming flows when there is more crowding.

5.1.1 Comparison with other literature

Given the unavailability of multiparticle correlations for non-spherical particles, we combine the available literature results on isolated non-spherical particles with voidage effects based on spheres. For this, we normalize the multiparticle drag of spheres with the isolated sphere Schiller & Naumann (1935) correlation and multiply with the isolated non-spherical particle drag. The results are shown in figures 20 and 21 for $\epsilon_{s}=0.3$ and $\epsilon_{s}=0.5$, respectively. The isolated particles drag law used are SKP (Sanjeevi et al., 2018a) and HS (Hölzer & Sommerfeld, 2008). They are accordingly combined with the multiparticle effects of TGS (Tenneti et al., 2011) and Tang et al. (2015) for spheres. In the moderately crowded regime ($\epsilon_{s}=0.3$), our earlier suggestion of combining isolated non-spherical particle drag with multiparticle effects from spheres works well. For example, the combination of SKP with Tang et al. (2015) follows nearly the same trend as that of the current work (equation 23). This can be observed for both $\phi=0^{\circ}$ and $\phi=90^{\circ}$. However for dense regimes ($\epsilon_{s}=0.5$), it can be observed that the combination of SKP with Tang et al. (2015) does not agree well with the present work for $\phi=0^{\circ}$. At the same time, the combination with the HS (Hölzer & Sommerfeld, 2008) isolated drag law seem to be closer to the current work for $\epsilon_{s}=0.5$. Such an agreement must be considered with care. The decent agreement occurs because HS possesses high drag values for $\phi=0^{\circ}$ (for $\epsilon_{s}=0$), in combination with a weak voidage effect for spheres. On the other hand, SKP with TGS or Tang et al. (2015) show decent agreement with the present work for $\phi=90^{\circ}$.

5.2 Lift

The normalized lift $F_{l,\phi}$ on a single elongated particle from Sanjeevi et al. (2018a) is given by

[TABLE]

Here, $S_{f,\phi}$ is the scaling function dependent on $\Rey$ and $\phi$. The coefficients $b_{i}$ are accordingly listed in the mentioned literature. In particular, the coefficients $b_{6}$ to $b_{9}$ describe the amount of skewness of the lift coefficient on a single elongated particle around $\phi=45^{\circ}$. In the current work, we observe the same skewness for the multiparticle system at different $\Rey$. Therefore, we assume the term $S_{f,\phi}$ remains the same for the multiparticle system. The normalized lift $F_{L}$ for a multiparticle system takes the following form:

[TABLE]

The functional form of $F_{L,mag}(\Rey,\epsilon_{s})$ remains similar to that of the drag and is given by

[TABLE]

with

[TABLE]

The corresponding coefficients are given in table 2. The proposed lift correlation has around 5% average absolute deviation with respect to the simulation results. The comparison of the $F_{L}$ from simulations and the proposed correlation is shown in figure 22.

5.2.1 A simplified lift function

In our earlier works (Sanjeevi & Padding, 2017; Sanjeevi et al., 2018a), we have shown successfully that for elongated particles, the relation between lift and drag in the Stokes flow regime can be successfully used for higher $\Rey$ flows too. In other words, $F_{L}$ at different $\phi$ can be computed as

[TABLE]

In this section, we show that equation 34 is a reasonable approximation even for a multiparticle system. This implies that the scaling law is valid not only just for different $\Rey$ but even for different $\epsilon_{s}$. Given a measured $F_{L}$ distribution from simulations at a given $\Rey$ and $\epsilon_{s}$, the data can be fitted in a simple form as

[TABLE]

Here, $F_{L,simple}$ is a fit parameter that best describes the simulation data. An example for such a fit for $\Rey=100$ and $\epsilon_{s}=0.3$ is given in figure 23. The comparison of the Stokes regime lift law (equation 34) and our hypothesis (equation 35) is shown in figure 24 and it can be observed that there is a good agreement. The highest absolute deviation observed between the equations is still less than 20% and average absolute deviation is around 12%. Therefore in Euler-Lagrangian simulations, in the absence of explicit lift data, equation 34 can be applied to include the effects of lift with acceptable accuracy. This implies that in the often-used approach of using Hölzer & Sommerfeld (2008) type drag correlations, combined with sphere-based voidage effect correlations in Euler-Lagrangian simulations, one can also include lift effects based on equation 34. In the following section, we will show the importance of including lift, as it is often of comparable magnitude to drag at high $\Rey$.

5.2.2 Importance of lift compared to drag

In Euler-Lagrangian simulations, the effect of lift forces is often neglected. This is because there is not much literature on non-spherical particle lift correlations. In this section, we analyse the magnitudes of lift compared to the drag on individual non-spherical particles at different $\Rey$ and $\epsilon_{s}$. Figure 25 shows the distributions of the magnitude of the lift force relative to the drag force on each particle $|F_{L}|/F_{D}$. It can be observed that for Stokes flow ($\Rey=0.1$), most particles experience lift which is about one order of magnitude smaller than the drag. However for high $\Rey$ ($\Rey=1000$), the distribution is much more wider spread and there are even some particles with $|F_{L}|/F_{D}=1$. This emphasizes the need for including lift in Euler-Lagrangian simulations, especially while handling Geldart D particles, where the encountered particle $\Rey$ is high. With increasing $\epsilon_{s}$, a different interesting observation can be made. In the low $\Rey$ regime, increasing $\epsilon_{s}$ results in an increased probability of particles experiencing high lift magnitudes compared to the drag. On the contrary, at high $\Rey$ ($\Rey=1000$), increasing $\epsilon_{s}$ results in the $|F_{L}|/F_{D}$ distribution skewing to the left. It should be noted that the highest $\epsilon_{s}$ shown in figure 25 is $\epsilon_{s}=0.4$ as opposed to $\epsilon_{s}=0.5$, the highest $\epsilon_{s}$ explored. This is because random configurations are not possible for $\epsilon_{s}=0.5$. To ensure consistency, all results shown in figure 25 are based on random configurations.

5.3 Torque

For an isolated non-spherical particle, the torque correlation (Sanjeevi et al., 2018a) is given by:

[TABLE]

The $\Rey$ dependent skewness terms $c_{5},c_{6},c_{7},c_{8}$ equal zero for an isolated spherocylinder resulting in a symmetric distribution for $\phi$ around $45^{\circ}$. Likewise, we also observe a near symmetric distribution of torque at different $\Rey$ and $\epsilon_{s}$ for the multiparticle configuration (see figure 26). Unlike drag and lift, for an isolated non-spherical particle, the pitching torque vanishes for all $\phi$ in the Stokes flow regime. We observe the same for the multiparticle configuration. Therefore, the proposed correlation for the torque $T_{P}$ is applicable only in the inertial regime ($10<\Rey\leq 1000$) and is given by

[TABLE]

The corresponding terms in the scaling are as follows:

[TABLE]

The average absolute deviation between equation 39 and corresponding simulation data is 3%. It should be noted that $T_{P,mag}$ in equation 40 maps only the magnitude of the torque for different $\Rey$ and $\epsilon_{s}$. The $\phi$ dependence is included separately with the $sine$ and $cosine$ terms. The comparison of $T_{P,mag}$ and the corresponding simulation data are given in figure 27. Given a symmetric form for $T_{P,\phi}$, the $T_{P,mag}$ is equal to twice the magnitude of $T_{P,\phi=45^{\circ}}$ since $\sin\phi\cos\phi=1/2$ at $\phi=45^{\circ}$. From figure 27, it can be observed that $T_{P,mag}$ roughly follows the same power law dependence on $\Rey$ for different $\epsilon_{s}$ because the slopes are similar. This is in contrast to the drag trends in figure 15, where the trend starts from zero slope at low $\real$ to a constant slope at high $\Rey$. The reason is that the torque vanishes at low $\Rey$ for all $\epsilon_{s}$. The distributions of torque $T_{P}$ for different $\Rey$, $\epsilon_{s}$ and $\phi$ are given in figure 26.

5.4 Flow histograms

In the previous sections, we discussed the influence of the flow on the hydrodynamic forces and torques on the particles. The flow around particulate assemblies can also be viewed as flow through a porous medium. In this section, we discuss the results of the influence of the particles on the flow distribution.

The probability distributions of the normalized axial flow velocities ($u_{ax}/u_{avg}$) at different $\Rey$ and $\epsilon_{s}$ for random configurations are given in figure 28. Here, the normalization is done against the average axial velocity $u_{avg}=u_{s}/(1-\epsilon_{s})$ rather than the superficial velocity $u_{s}$ to ensure a fair comparison for different $\epsilon_{s}$. Only the velocities of fluid cells are shown here and the zero velocities in the solid cells are ignored. It can be observed that with increasing $\Rey$, the spread of the velocity distribution becomes narrower. This can be simply attributed to the increased inertial effects and thinner boundary layers for increasing $\Rey$. Interestingly, the high $\Rey$ flows also demonstrate some negative velocities corresponding to wake effects. With increasing $\epsilon_{s}$, the peaks of the distribution shift towards the left and the distribution itself spreads wider. This implies that the increased presence of particle surfaces at higher $\epsilon_{s}$ pulls the velocities of fluid cells towards zero (hence the left skewness). At the same time, the fluid accelerates in the bulk regions further removed from the particle surfaces resulting in increased velocities (and hence a wider distribution) to maintain the desired $u_{s}$.

It is also interesting to investigate the velocity distributions for different configurations for a given $\Rey$ and $\epsilon_{s}$. The distributions of $u_{ax}/u_{avg}$ at $\Rey=100$ and $\epsilon_{s}=0.3$ for different configurations are plotted in figure 29. Given sufficient randomness of particles, as in random and planar random configurations (see figure 29 ($a$) and ($b$)), the velocity distributions are nearly identical. However, velocity distributions can be different for different configurations, as can be observed for the unidirectional configurations with flow parallel and perpendicular to the principal director (see figure 29 ($c$) and ($d$)). Among the different configurations shown, the unidirectional configuration with flow parallel to principal director has the least recirculation, as is evident from the least number of fluid cells with negative velocities ($u_{ax}/u_{avg}<0$). At the same time, as expected, the unidirectional configuration with flow perpendicular to principal director has the highest amount of recirculation. Overall, we can infer that there is no dependency between the configuration independence phenomenon, explained in section 4, and the flow velocity distribution of different configurations. The variation in forces at different incident angles $\phi$ is mainly arising from the pressure forces. The same can also be confirmed from the multiparticle work of He & Tafti (2018), which is also in line with our finding for isolated non-spherical particles (Sanjeevi & Padding, 2017).

6 Conclusion

The flow around assemblies of axisymmetric, non-spherical particles has been studied extensively using the multi-relaxation-time lattice Boltzmann method. The performed simulations are from the Stokes flow regime to high $\Rey$ ($0.1\leq\Rey\leq 1000$) at different solids volume fraction $\epsilon_{s}$ ($\epsilon_{s}\leq 0.5$) and different mutual orientations of particles.

In general, forces on random assemblies of spheres are only dependent on $\Rey$ and $\epsilon_{s}$. Considering the non-spherical nature of the particles, we proposed four additional parameters to describe the flow problem: two to parametrize the mutual orientation of the non-spherical particles ($S_{1}$ and $S_{2}$) and two to represent the polar and azimuthal angles ($\alpha$ and $\beta$) of the averaged flow velocity with respect to the configuration. For this, we have developed different static particle configurations using Monte-Carlo simulations. In the process, the configurations are biased to the desired amount of nematic or bi-axial orientational order with the use of a Maier-Saupe potential. The flow simulations indicate that the particle forces are configuration independent, at least for $\epsilon_{s}\leq 0.4$, implying that the four additional parameters do not influence the results. The only important parameter representing orientation dependence is the incident angle $\phi$ of individual particles with respect to the average flow direction.

The configuration independence greatly simplifies the parameter space to be explored from 6 to 3 dimensions, namely $\Rey$, $\epsilon_{s}$, and $\phi$. Of the three, the simulations are set up for only two parameters: $\Rey$ and $\epsilon_{s}$. Given a sufficiently random particle configuration, different incident angles $\phi$ are covered automatically. Another interesting result from the current work is that our previous finding of sine-squared scaling of drag for isolated non-spherical particles (Sanjeevi & Padding, 2017) applies also to multiparticle systems. In other words, given a $\Rey$ and $\epsilon_{s}$, the drag on the subset of particles oriented at an incident angle $\phi$ with respect to the superficial flow velocity can be described with the knowledge of drag at $\phi=0^{\circ}$ and $\phi=90^{\circ}$ alone. This information can be used in a packed bed to determine the pressure drop across the bed with the knowledge of $\phi$ distribution alone. In a multiparticle configuration, also the lift on a particle at an incident angle $\phi$ can be computed with good accuracy using the drag at $\phi=0^{\circ}$ and $\phi=90^{\circ}$, as in our previous work on isolated non-spherical particles. Having identified the dependent parameters, we proposed correlations for drag, lift and torque for a multiparticle configuration. During the process, we used correlations for isolated non-spherical particles and extended them to the multiparticle systems.

We have also explored the validity of the conventional approach of combining known correlations for isolated non-spherical particle drag with correlations for voidage effects based on sphere packings. We observe that in the dilute and intermediate $\epsilon_{s}$ regimes ($\epsilon_{s}\leq 0.3$), the influence of $\epsilon_{s}$ is nearly shape independent. This implies that the above conventional approach can safely be used to mimic flow around assemblies of non-spherical particles upto intermediate $\epsilon_{s}$. However for denser regimes, there is a need for multiparticle simulations and hence the need for this work. In the inertial regimes, the ratios of drag at $\phi=90^{\circ}$ and $\phi=0^{\circ}$ ($F_{D,\phi=90^{\circ}}/F_{D,\phi=0^{\circ}}$) are nearly constant until $\epsilon\leq 0.3$ and then decrease with increasing $\epsilon_{s}$. This further proves that the conventional approach is not valid for dense regimes. In the process, we have analysed the flow velocity distribution as function of $\Rey$ and $\epsilon_{s}$. Likewise, the influence of different particle configurations on the flow velocities have also been analysed.

Overall, this work provides a recipe to parametrize the drag, lift and torque experienced by non-spherical particles in multiparticle environment. To the best of the authors’ knowledge, there exists no work which parametrizes the drag, lift and torque for non-spherical particles in a multiparticle environment. Generally, lift and torque are ignored in large scale Euler-Lagrangian simulations. The proposed accurate drag, lift and torque correlations enable future Euler-Lagrangian simulations to be performed with more realistic physics.

Acknowledgements

The authors thank the European Research Council for its financial support under its consolidator grant scheme, contract no. 615096 (NonSphereFlow). The simulations are performed in LB3D code with improved boundary conditions. We thank Prof. J. Harting for the code framework. SKPS thanks Y. El Hasadi for the fruitful discussions on fitting the data and I. Mema for the data in figure 11. The work has been made possible by a grant for computation time, project number SH-351-15, financed by the Netherlands Organisation for Scientific Research (NWO).