Hydrodynamic forces on monodisperse assemblies of axisymmetric elongated particles: Orientation and voidage effects

Abstract We investigate the average drag, lift, and torque on static assemblies of capsule‐like particles of aspect ratio 4. The performed simulations are from Stokes flow to high Reynolds numbers (0.1 ≤ Re ≤ 1,000) at different solids volume fraction (0.1 ≤ ɛ s ≤ 0.5). Individual particle forces as a function of the incident angle ϕ with respect to the average flow are scattered. However, the average particle force as a function of ϕ is found to be independent of mutual particle orientations for all but the highest volume fractions. On average, a sine‐squared scaling of drag and sine‐cosine scaling of lift holds for static multiparticle systems of elongated particles. For a packed bed, our findings can be utilized to compute the pressure drop with knowledge of the particle‐orientation distribution, and the average particle drag at ϕ = 0° and 90°. We propose closures for average forces to be used in Euler–Lagrange simulations of particles of aspect ratio 4.


| INTRODUCTION
Accurate fluid-particle drag, lift, and torque closures are required for precise Euler-Lagrangian simulations of nonspherical particles. Historically, different drag closures have been developed for assemblies of spherical particles. [1][2][3] However, practical flows often involve assemblies of nonspherical particles for which there exist no closures at the moment. Even for static, monodisperse, nonspherical 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, that is, the typical relaxation time of the solid particle velocity is large relative to the response time of the gas. 4 It has been shown that under such conditions, it is sufficient to assume the particle configurations to be quasistatic. 5 Conventionally, fluidization simulations of nonspherical 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 nonspherical particles. Hölzer and ing complex, unsteady flows. 4 In an earlier work, 9 we reported the interesting finding that the drag coefficient C D at different incident angles ϕ follows a sine-squared scaling given by C D,ϕ = C D,ϕ = 0 + C D,ϕ = 90 − C D,ϕ = 0 sin 2 ϕ: Likewise, we reported another interesting finding that the lift coefficient C L follows sine-cosine scaling at different ϕ as C L,ϕ = C D,ϕ = 90 − C D,ϕ = 0 sinϕcosϕ ð2Þ for various elongated particles. The abovementioned 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. 9 In Equations (1) and (2), the drag coefficients at incident angles of 0 and 90 still depend on particle shape and Reynolds number. The Reynolds number in the present work is defined as Re = |u s |d eq /ν, where u s is the superficial flow velocity, ν is the kinematic viscosity of the fluid, and d eq is the diameter of the volume-equivalent sphere given by d eq = (6V p /π) 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, 10 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  scheme to much more accurately represent the particle geometry.
Recently, Tang et al 3 used an IBM based solver to create drag closures for static assemblies of spheres upto Re ≤ 1,000 and ϵ s ≤ 0.6. We note that all mentioned works report their drag closures as the average drag on a collection of particles (typically a hundred to a few hundred) as a function of the average solids volume fraction. In reality, variations in local volume fraction and precise placement of neighboring particles will lead to a scatter in the drag per particle. However, these closures are developed for use in unresolved Euler-Lagrange (CFD-DEM) simulations, where a typical CFD cell will be as large as the entire resolved simulation box (i.e., with a cell size typically equal to 3-6 particle lengths). It is true that in reality individual drag forces can be higher or lower than the average drag, but such detail is generally not taken into account in Euler-Lagrange simulations for computational efficiency. In general, it is assumed that the particle-particle interactions (collisions) will lead to a rapid redistribution of particle velocities within a cell, making the average drag the most relevant factor.
There are also several disadvantages with combining an isolated nonspherical 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. 1,15 Second, the voidage effects on lift and torque in a multiparticle system are unknown and hence are often neglected in Euler-Lagrangian simulations. 16,17 Third, using the same factor for voidage effects for all incident angles ϕ may hold in sufficiently dilute regimes but its validity in the dense limit is unknown. At the moment, only Li et al 18

| Lattice Boltzmann method
In the present work, we use a D3Q19, MRT lattice Boltzmann method 21 to simulate the fluid flow. The numerical method is adequately explained and validated in our previous works. 4,9 The evolution of particle distribution function jf〉 is computed as for position r with discrete velocities e α in directions α = 1, 2…, 19.
Equation (3) 4 The density is computed as ρ = P α f α and the momentum as ρu = P α f α e α . The relation between the kinematic viscosity of the fluid and the dimensionless relaxation time τ is ν = c 2 s τ − 1=2 ð Þ Δt, and the pressure p is related to the density by p = ρc 2 s , where c s is the speed of sound. A linearly interpolated bounce back scheme 23,24 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 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. 25 The results for the multiparticle system are validated in subsequent sections.
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 Re. Further with increasing ϵ s , the d eq is increased accordingly to resolve increased velocity gradients at high ϵ s .
All LBM simulations have cubic domain, each with 200 particles unless otherwise specified. At least two independent simulations are performed for each Re and ϵ s and the details of independent number of simulations are discussed later (see Figure 13).

| Flow control
In order to perform a simulation for a specific Re, it is required to control the superficial flow velocity u s by applying a body force g. The relationship between the superficial velocity and the average interstitial flow velocity u avg is given by u s = (1 − ϵ s )u avg . Due to the nonspherical nature of the particles, the sum of lift forces is often non-zero, and the resultant direction of u s can be different from the direction of g. This necessitates the need to control both direction and magnitude of the body force. Initially, the fluid is at rest with both u s and g zero. The flow is slowly ramped up by increasing g until the desired u s is achieved. For each timestep, the updated gravity g new is computed as where g prev is the gravity from the previous timestep, u s, ref is the desired superficial velocity, and 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 u s reaches 99.9% of the reference setpoint.

| Orientation tensor
In this section, we briefly explain the characterization of mutual orientations in an assembly of axisymmetric nonspherical particles with an Note: L D denotes the side length of the cubic domain. The range of d eq specified is respectively for 0.1 ≤ ϵ s ≤ 0.5. orientation tensor. 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 S, also known in literature as the nematic order tensor, 18,26,27 defined as Here, p is the unit orientation vector of a particle along the axis of symmetry. The three 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 three eigenvectors define the principal directions of mutual particle alignment.
Because the trace of S is 1, only two eigenvalues are sufficient to specify the amount of randomness, planar random (biaxial), or unidirectional (nematic) order. Note that the tensor S is insensitive to an orientation p or − 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 1a shows a completely random configuration with S 1 = S 2 = S 3 = 1/3. Figure 1b 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 1c with S 1 = S 2 = 0, S 3 = 1. In practical conditions, particles can exhibit complex configurations in between these extremes but can be adequately described by two 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 (α) and azimuthal angle (β) of the average flow velocity vector with respect to the space spanned by the three eigenvectors of the orientation tensor. In summary, the parameter space to be explored for our flow problem has six parameters, namely Reynolds number Re, solids volume fraction ϵ s , two particle configuration parameters S 1 , S 2 and two angles α and β describing the mean flow orientation with respect to the configuration.

| Generation of biased particle configurations
The generation of nonoverlapping configurations of the particles in a periodic domain is required as an input for the flow simulations. 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 nonoverlapping particles and the use of a Maier-Saupe potential 28 to bias the system to 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 (i.e., sum of the radii of two interacting particles), then the spherocylinders overlap. A fast algorithm is used to measure the shortest distance between the line segments. 29 Using the above overlap detection algorithm, randomly picked particles are randomly translated in small (compared to the particle diameter) steps and rotated by a small angle around a randomly chosen axis. Because our system is always below the threshold for a spontaneous nematic order transition, this procedure results in a random configuration after many iterations. If a prescribed amount of mutual orientation is required, besides the requirement of no overlap, a Monte-Carlo procedure is applied to decide whether to accept or reject a new orientation of a particle. In detail, we choose a principal director n, which is a reference vector to which the particles are biased to align with or against (depending on the sign of the magnitude A of the Maier-Saupe potential). In the Monte-Carlo approach a new orientation p new of a randomly picked particle, having current orientation p curr , is accepted or rejected based on the following criteria: where Here, ΔE is the change in Maier-Saupe potential and U([0, 1]) is a random number uniformly distributed between 0 and 1. Of the three conditions in Equation (6), it is clear that the first condition accepts the new orientation if it leads to a lower Maier-Saupe potential. Without the second condition, the system would approach toward an ideal mutual orientation (such as perfect parallel or perfectly perpendicular particles w.r.t. the principal director, depending on the sign of A) when the Monte-Carlo simulation is run for a sufficiently long time.
With the second condition, however, increases in the Maier-Saupe potential are also accepted with a certain probability less than 1 (the larger the increase the potential, the smaller the probability of acceptance). After sufficiently long time, a balance is found between the random particle reorientations and particle orientation ordering by the Maier-Saupe potential, leading to a degree of randomness that can be controlled by the magnitude of the user specified constant A. A bias toward planar random configuration is achieved when A > 0, with more particles oriented in planes perpendicular to the director n. A bias toward unidirectional (nematic) configuration is achieved when A < 0, with more particles oriented along the direction of n.
With the mentioned strategy, any configuration in-between the ideal cases shown in Figure 1 can be achieved. Some sample configurations generated using the abovementioned algorithm are shown in A common intuition may be that a random configuration would result in particles with evenly distributed values of the incident angle ϕ. However, for a random configuration, the available number of particles at different ϕ are not uniform, as shown in Figure 3a. This is due to the higher probability to find particles at an angle ϕ near 90 because the Jacobian for a spherical coordinate system scales as sinϕ.
Therefore, the disadvantage for a random configuration is that there are actually few data points at ϕ = 0 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 3b. This information is considered while we generate configurations for the flow simulations.

| Forces and torques acting on a particle
For an assembly of particles, different definitions are used to report the forces. [1][2][3] 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 rP, each particle of volume V p experiences a resulting force F due to the flow alone, and a buoyancy force F b = −V p rP due to the pressure gradient.
For such a case, the total fluid-to-particle force F f ! p acting on a particle is Given N particles with each of volume V p and total volume of the system V, the solids volume fraction is given by ϵ s = NV p /V. Further, the relationship between F and F f ! p is given by 3 In this work, we report the forces F due to the flow and not F f ! p .
Note that in some simulation packages F f ! p is needed, in which case the correlations we report in this work should be divided by (1 − ϵ s ).
The effects of buoyancy on torques are unknown and hence the reported torques 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: Here, μ is the dynamic viscosity and R eq is the radius of the volume equivalent sphere. The Stokes torque that we used is based on the torque experienced by a rotating sphere in still fluid. 30 Let p be the normalized orientation vector of the considered particle. The local coordinate system for each particle is defined aŝ The above defined axes are accordingly illustrated in Figure 4.
The incident angle ϕ a particle makes with respect to the incoming flow is given by ϕ = cos − 1ê 1 Áp j j ð Þ. We also compute the average forces and torques for different ϕ intervals. Due to the finite number of measurements in these intervals, there is an error on the mean x of any property x. We use the standard error on the mean σ x for the errorbars, computed as Here, σ is the standard deviation of the corresponding variable x and n is the number of data points within the given ϕ interval.
Throughout this work, we use overbar ( -) to denote arithmetic averages and boldface to denote vectors.
The normalized drag F D and lift F L can be computed from F norm as Since the reported forces are without buoyancy effects, the (1 − ϵ s ) term must be considered accordingly for both drag and lift while performing Euler-Lagrangian simulations. Due to the influence of neighboring particles, the lateral force F 2 for each individual particle may not be equal zero, as shown in Figure 5 (Re = 100 and ϵ 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 T P = T 2 = T norm Áê 2 , and ð20Þ Here, T P is the pitching torque acting on a particle. We show the three different torques for a flow through a random particle F I G U R E 4 The local coordinate system of a particle. u s and F D act alongê 1 , F L alongê 3 and T P about theê 2 axis  stances, the relative density variations were observed to be at most 2%, which is why our results can be considered to be in the incompressible limit.

| Tests of configuration independence
Given a six-dimensional parameter space, exploring each dimension with approximately five simulations, results in a massive 5 6 = 15,625 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 average hydrodynamic force on a nonspherical particle 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 α and β 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 ϕ as in flow around single particles. Effectively, we will show that the average force depends only on the Reynolds number Re, In other words, the average drag F D at any ϕ can be computed as It is important to note that the same values for average F D,ϕ = 0 and F D,ϕ = 90 emerge for all configurations: the solid lines in Figure 8 are obtained as a single fit to the data from all configurations investigated at a certain ϵ s . Likewise, we also show that the scaling phenomenon extends to both Stokes and high Re regimes in Figure 9.
With the sine-squared scaling behavior (or the configuration independence) identified at ϵ s = 0 and ϵ s = 0.3, it can be inferred that the scaling is safely applicable in the region 0 ≤ ϵ s ≤ 0.3. We have verified the same at ϵ s = 0.1 and the results are not shown here for brevity.
Though we observe the results are dependent on only three parameters, namely Re, ϵ s , and ϕ, the simulation needs to be set up for only two parameters, namely Re and ϵ s . With a sufficiently random configuration, the system involves different particle orientations covering all ϕ. A caveat with a random configuration is that there are always very few particles near ϕ = 0 , as shown in Figure 3. Therefore, biased random configurations with more particles at ϕ = 0 are created and at least two simulations are performed for better statistics.
We also observe the configuration independence phenomenon at In such a dense condition, the particle configuration itself is dependent on the wall geometry. we also observe the same for a bed containing freely poured particles settled under gravity (ϵ s = 0.54), as shown in Figure 11. The bed contains 30,000 particles and it can be observed that roughly 2/3 of all particles are in the range ϕ = 70-90 confirming our hypothesis.
Given such criteria, the most relevant regime would be to generate an accurate fit for average F D,ϕ = 90 at high ϵ s , which would help to predict minimum fluidization velocity of the bed accurately.
It should also be noted that with increasing aspect ratio of elongated particles, the maximum ϵ s decreases for a packed bed. 32 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 ϵ 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 ϵ s implies that the configuration independence phenomenon will be very applicable.
With the observed sine-squared drag scaling, the pressure drop

| 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 Re = 100 and ϵ s = 0.3 is shown in Figure 12.

| Drag
With sine-squared scaling valid for all particle mutual orientations, as shown in the previous section, the average drag experienced by a particle in a multiparticle system can be explained by the Equation (22) involving only the average drag experienced at ϕ = 0 and ϕ = 90 .
Therefore, we propose to generate fits for average F D,ϕ = 0 and F D,ϕ = 90 as a function of Re and ϵ s as The corresponding terms are as follows: Here, C d, isol is the isolated particle drag at given Re as detailed in Reference 4 for the considered particle (fibre) for both ϕ = 0 and ϕ = 90 . The coefficients in Equations (25) and (26) for both average F D,ϕ = 0 and F D,ϕ = 90 are given in Table 2. The average absolute deviation of the fits and simulation data are 3.5 and 2% for F D,ϕ = 0 and F D,ϕ = 90 , respectively.
The simulated data and corresponding fits are shown in Figure 14. The fits follow the physical limits beyond the Re range simulated as shown in Figure 15. In the Stokes flow limit, it can be observed that both ϕ = 0 and ϕ = 90 normalized drag becomes independent of Re. In the high Re limit, the normalized drag approaches a linear dependency on Re.
The ratio of average perpendicular to average parallel drag F D,ϕ = 90 = F D,ϕ = 0 at different Re and ϵ s is shown in Figure 16. For low Re (Re = 0.1), the ratio remains constant at a value a little larger than 1 for all ϵ s . The reason for this is that at low Re, the particles experience stronger viscous effects. The viscous drag reduces and pressure drag increases with increasing ϕ at low Re. The same has been confirmed for isolated particles 9 and for a multiparticle system. 19 The combined viscous and pressure drag components result in a drag ratio close to 1 for the considered spherocylinders at low Re. Due to inertial dominance at moderate and large Re (Re ≥ 100) we can observe a near constant drag ratio for solids volume fractions upto ϵ s = 0.3 and a decrease in the ratio for ϵ s > 0.3. Further, Figure 16 gives an indication that for very dense crowding, that is, at ϵ s > 0.5, there is a possibility that  Re may not apply. Therefore, the above discussion indicates that spherical drag correlations for the voidage effect, combined with isolated nonspherical particle drag correlations can be applied to dilute suspension simulations of nonspherical particles in the inertial regimes. For a given nonspherical particle, the effect of crowding (ϵ s ) on F D is different for different Re and ϕ. Figure 18 shows the voidage effect (average F D normalized by the corresponding isolated particle drag) as a function of Re. It can be seen at low Re, the increase in drag due to crowding is comparable for both ϕ = 0 and ϕ = 90 at different ϵ s . At high Re, the increase in drag due to crowding with increasing ϵ s is much stronger for ϕ = 0 compared to ϕ = 90 . This also explains further the reason for the observed reduction in average perpendicular to average parallel drag ratios with increasing ϵ s in In the previous sections, we discussed the F D averaged over all particles with similar ϕ. However, the distribution of F D within a ϕ interval is itself also a function of both Re and ϵ 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 average F D at respective ϕ, rather than against a single value, say F D,ϕ = 90 , for a given Re and ϵ s . This is because with increasing Re, the ratio F D,ϕ = 90 = F D,ϕ = 0 increases, as shown in Figure 16. Therefore, using average F D,ϕ = 90 for normalization will make the standard deviations at ϕ = 0 appear insignificant at large Re.

| Lift
The normalized lift F l, ϕ on a single elongated particle from Sanjeevi et al 4 is given by Re b5 Re 24 , and ð28Þ Here, S f, ϕ is the scaling function dependent on Re and ϕ. 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 ϕ = 45 . In the current work, we observe the same skewness for the multiparticle system at different Re. Therefore, we assume the term S f, ϕ remains the same for the multiparticle system. The average lift F L for a multiparticle system takes the following form: The functional form of F L, mag (Re, ϵ s ) remains similar to that of the drag and is given by The corresponding coefficients are given in Table 2

| A simplified lift function
In our earlier works, 4,9 we have shown successfully that for isolated elongated particles, the relation between lift and drag in the Stokes flow regime can be successfully used for higher Re flows too. In other words, F L at different ϕ can be computed as 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 Re but even for different ϵ s . Given a measured average F L distribution from simulations at a given Re and ϵ s , the data can be fitted in a simple form as Here, F L, simple is a fit parameter that best describes the simulation data. An example for such a fit for Re = 100 and ϵ 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 and Sommerfeld 6 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 Re.

| 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 nonspherical particle lift correlations. In this section, we analyse the magnitudes of lift compared to the drag on individual nonspherical particles at different Re and ϵ 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 (Re = 0.1), most particles experience lift which is about one order of magnitude smaller than the drag.
However, for high Re (Re = 1,000), the distribution is much more

| Torque
For an isolated nonspherical particle, the torque correlation 4 is given by: All coefficients can be found in our previous work. 4 We note that for our particle geometry the isolated particle torque strictly increases with increasing Re (at least in the range of Re studied). It may be possible that at higher Re the torque will decrease again, as predicted by Khayat and Cox 34 for slender bodies.
The Re dependent skewness terms c 5 , c 6 , c 7 , c 8 equal zero for an isolated spherocylinder resulting in a symmetric distribution for ϕ around 45 . Likewise, we also observe a near symmetric distribution of average torque at different Re and ϵ s for the multiparticle configuration (see Figure 26). Unlike drag and lift, for an isolated nonspherical particle, the pitching torque vanishes for all ϕ in the Stokes flow regime. We observe the same for the multiparticle configuration.
Therefore, the proposed correlation for the average torque T P is applicable only in the inertial regime (10 < Re ≤ 1,000) and is given by T P,ϕ Re,ϵ s , ϕ ð Þ= T P,mag Re,ϵ s ð ÞÁsinϕcosϕ, with ð39Þ The corresponding terms in the scaling are as follows (coefficients for the fit are given in Table 3): 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 Re and ϵ s . The ϕ 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,ϕ , the T P, mag is equal to twice the magnitude of T P,ϕ = 45 since sinϕcosϕ = 1/2 at ϕ = 45 . From Figure 27, it can be observed that T P, mag roughly follows the same power law dependence on Re for different ϵ s because the slopes are similar. This is in contrast to the drag trends in

| 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 Re and ϵ s for random configurations are given in Figure 28.  The configuration independence greatly simplifies the parameter space to be explored from 6 to 3 dimensions, namely Re, ϵ s , and ϕ. Of the three, the simulations are set up for only two parameters: Re and ϵ s . Given a sufficiently random particle configuration, different incident angles ϕ are covered automatically. Another interesting result from the current work is that our previous finding of sine-squared scaling of drag for isolated nonspherical particles 9 applies also to static monodisperse assemblies containing axisymmetric, elongated particles. In other words, given a Re and ϵ s , the average drag on the subset of particles oriented at an incident angle ϕ with respect to the superficial flow velocity can be described with the knowledge of average drag at ϕ = 0 and ϕ = 90 alone. This information can be used in a packed bed to determine the pressure drop across the bed with the knowledge of ϕ distribution alone. In a multiparticle configuration, also the average lift on a particle at an incident angle ϕ can be computed with good accuracy using the average drag at ϕ = 0 and ϕ = 90 , as in our previous work on isolated nonspherical particles. Having identified the dependent parameters, we proposed correlations for average drag, lift, and torque for a multiparticle configuration of aspect ratio 4 spherocylinders. During the process, we used correlations for isolated nonspherical particles and extended them to the multiparticle systems.
We have also explored the validity of the conventional approach of combining known correlations for isolated nonspherical particle drag with correlations for voidage effects based on sphere packings.
We observe that in the dilute and intermediate ϵ s regimes (ϵ s ≤ 0.3), the influence of ϵ s is nearly shape independent. This implies that the above conventional approach can safely be used to mimic flow around assemblies of nonspherical particles upto intermediate ϵ 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 average drag at ϕ = 90 and ϕ = 0 ( F D,ϕ = 90 = F D,ϕ = 0 ) are nearly constant until ϵ ≤ 0.3 and then decrease with increasing ϵ 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 Re and ϵ s . Likewise, the influence of different particle con-