Gas flow through static particle arrangements with a channel Resolved simulations and analytic considerations

Fractures of particle assemblies happen frequently in dense gas-solid systems leading to a notable heterogeneity in the particle conﬁguration, especially in case of cohesive powders and non-spherical particle interlocking. In this work, we investigate the inﬂuence of such heterogeneities on the hydrodynamic drag by studying the idealized case of a random arrangement of spheres with a channel-like void region. More speciﬁcally, we introduce this heterogeneity to a homogeneous particle arrangement by shifting apart two bulk regions, such that a void channel divides particle bulk. Single-relaxation-time lattice Boltzmann simulations were performed to resolve ﬂuid ﬂow through such arrested particle conﬁgura- tions and calculate the corresponding gas-particle momentum exchange and pressure drop. The calculated drag forces acting on the solids for random sphere arrangement are in good agreement with previously reported results of Hill et al. (2001b), Tenneti et al. (2011), and Tang et al. (2015). However, the overall momentum exchange obtained for conﬁgurations containing a heterogeneity is signiﬁcantly lower. Obviously, the channel allows for a by-passing of a considerable amount of the ﬂow leading to a reduced overall pressure drop and thereby underestimating the minimum ﬂuidization velocity in a ﬂu- idized bed. Based on these direct numerical simulations, we examine the overall pressure drop dependence on the characteristic length scale (i.e. width) of the channel-like heterogeneity L c , the superﬁcial Reynolds number (30 6 Re 6 300), and the solid volume fraction in the dense (i.e. bulk) region (0.4 6 / p 6 0.55). The width of the channel is varied within the order of magnitude of particle diameter D p (1 6 L c = D p 6 4 : 36), decreasing an overall solid volume fraction (0.25 6 / 6 0.55). In addition to the numerical simulations, we derive (semi)-analytic correlations for the dense bulk region as well as for the channel. As the simulations range from laminar to transitional ﬂow, providing a single pressure drop correlation is very challenging. Therefore, we estimate the channel pressure drop with the appropriate correlations selected according to calculated superﬁcial Reynolds number. For laminar ﬂow, we achieved a good agreement between a combined (i.e. bulk and channel) analytical prediction of overall pressure drop and our resolved numerical simulation. In the transitional regime, the pressure drop values are more difﬁcult to predict, with the better agreement as we approach the turbulent regime. We believe that this work can act as a basis for the development of future drag laws accounting for channel-like sub-grid heterogeneities. Elsevier creativecommons.org/licenses/by/4.0/).


Introduction
Gas-particle flows often occur in various industrial applications containing chemical (e.g. combustion, gasification) and physical processes (e.g. granulation, segregation, coating, and drying) or natural phenomena. These operations make use of advantageous properties such as good fluid-solid contact characteristics (Gilliland and Mason, 1949;Mickley and Trilling, 1949). Fundamental and quantitative understanding of gas-particle interactions is necessary for the efficient processes and design improvement.
Computational fluid dynamics (CFD) simulations solve the averaged multiphase flow equations and can augment the experimental studies supplying the detailed, otherwise obscured, data. Numerous studies accomplished these simulations with either the continuum Eulerian-Eulerian (EE) two-fluid approach (Anderson and Jackson, 1967;Gidaspow, 1994;Schneiderbauer et al., 2013) or the Euler-Lagrange (EL) approach (Zhou et al., 2010;Tsuji et al., 2014;Lichtenegger and Pirker, 2018). The Euler-Lagrange coupling can either be resolved on unresolved. The latter combines a discrete element method (DEM) for solids with a continuum approach for the fluid phase where the Eulerian grid is larger than the particle size. In these models, the averaged equation of motion is solved for the gas phase. When the flow past the particles is not resolved, a term in the mean momentum conservation equation representing the average interphase momentum transfer between the fluid phase and the solid particles is required. This term is commonly referred to as the drag law, or the drag correlation.
Different correlations for the drag laws where proposed from theoretical, experimental and numerical investigations. Various theoretical and experimental studies exist for Stokes flow through porous media, packed bed of solids or ordered arrays of spheres (Kozeny, 1927;Hasimoto, 1959;Sangani and Acrivos, 1982). As theoretical methods are typically limited to low solid volume fractions and low Reynolds number, experimentally based models have a wide use when modeling the granular flow. Some of the commonly employed empirical correlations based on experimental pressure drop measurements are those by Ergun (1952) for dense systems and Wen and Yu (1966) for dilute systems. Gidaspow (1986) proposed an improved correlation combining these two models, considering the particle bed porosity. In contrast to that, Syamlal and O'Brien (1987) for instance derived their correlation for a single particle and modified it with a relative velocity correlation.
The study of the drag is more complicated for random arrays of particles and a wider distribution of Re. With the increase of computational power, DNS has become a powerful tool for directly resolving the flow past solid particles which allows for more accurate quantification of the gas-particle forces. Numerous studies used the lattice Boltzmann method (LBM) to investigate the drag forces and the pressure drop over a wide range of Re and particle volume fractions / (Ladd, 1994;Hill et al., 2001b;Van der Hoef et al., 2005;Bogner et al., 2015). The immersed boundary method (IBM) was likewise applied for the case of fixed or moving arrays of randomly arranged monodisperse spheres (Tenneti et al., 2011;Tang et al., 2015). Tenneti et al. (2011) reported the dimensionless average fluid-particle force for random assemblies of monodisperse spheres, extending the previous work of Hill et al. (2001b) to a wider range of solid volume fractions / and Reynolds numbers Re. Tang et al. (2016) analyzed the influence of particle mobility on the gas-solid drag force and modified the existing drag correlation they obtained from simulations of stationary particles. Kriebitzsch et al. (2013) compared the finite-resolution fully resolved simulation with the drag force values on a sphere in an ordered array (as calculated by Hasimoto (1959)). They found a deviation between the results even for low Re where it is possible to obtain the exact solution.
The majority of DNS studies on drag forces deals with homogeneous systems of random spatially fixed particle arrangements, where correlations are typically given as a function of Re and /. Experimental and numerical studies of gas-particle systems such as fluidized beds found the tendency to form channels for cohesive (Baerns, 1996;Geldart, 1973;Pacek and Nienow, 1990;Raganati et al., 2018) and non-spherical particles (Liu et al., 2008;Vollmari et al., 2015;Mahajan et al., 2018) due to inter-particle forces and interlocking, respectively. In this work, we introduce a channel-like void region to the particle configuration. Channel through a packed bed allows a by-passing of a considerable part of the flow, leading to reduced forces on the particles. In larger systems, the formation of heterogeneous structures can occur on a scale smaller than the grid resolution of a few particle diameters, typically used in CFD-DEM simulations. Therefore, in an unresolved simulation, the entire here resolved domain might be represented by a single cell where only overall solid volume fraction of the cell is known. Applying inadequate drag law would lead to a vast drag overestimation for heterogeneous structures. Traditional drag law closures might not be applicable to heterogeneous configurations without understanding the sub-grid heterogeneities. Consequently, there is a clear need for a new correlation that accounts for such heterogeneous particle arrangement.
Some of the approaches dealing with the influence of heterogeneous formations on the drag force are the energy minimization method (EMMS) (Li and Kwauk, 1994;Liu et al., 2001;Wang and Li, 2007) and coarse grid filtering of the fine grid simulations (Yang et al., 2003;Igci et al., 2008;Parmentier et al., 2012;Schneiderbauer et al., 2013;Schneiderbauer and Pirker, 2014). Ma et al. (2009) performed simulations adopting Lagrangian-Lagrangian schemes to observe and understand the mechanism of particle clustering and its formation. Zhou et al. (2014) quantified the drag force dependence for heterogeneous flow past spheres over distinct dilute-dense regions (stepwise heterogeneity). For this work, however, we examine higher solid volume fractions.
In this paper, we perform LBM simulations to study the flow past heterogeneous particle assemblies containing a channel-like void region for various flow parameters (Re and /). Heterogeneity is generated in form of a channel-like fracture dividing the particle bulk into two parts which would not be captured in an unresolved simulation where cells are larger than particles. Such structure leads to higher fluid velocities in a channel-like region and lower velocities in the surrounding dense particle phase, highly influencing the pressure drop and drag forces.
The paper is organized as follows. First, we introduce the applied numerical method (LBM) for simulation of the flow (Section 2). Then, we present several different simulation setups (Section 3) and the results providing the values of the particle-fluid forces (Section 4). With the numerical findings presented, we make a semi-analytic approach combining the existing correlations for pressure drop through particle bed and channels. (Section 5). Finally, we state a conclusion and an outlook of future work and improvements (Section 6).

Lattice Boltzmann method
The lattice Boltzmann equation (LBE) originated from Ludwig Boltzmann's kinetic theory of gases (Guo et al., 2000;Bao and Meskas, 2011). The LBM applies a mesoscopic simulation approach where instead of directly solving the macroscopic fluid properties (i.e. pressure and velocity), one models the evolution of discrete particle distribution functions. The exchange of momentum is achieved through particle collision and streaming and is modeled by the Boltzmann transport equation: where f ðx; tÞ is the particle distribution function, u is the particle velocity and X is the collision operator.
The LBE can be viewed as a particular discrete form of the continuum Boltzmann equation (Eq. (1)). The domain is discretized in uniform Cartesian cells that hold a fixed number of distribution functions. For each particle on the lattice, such a discrete probability distribution function describes the probability of streaming in one particular direction (Succi, 2001;Iglberger et al., 2008;Bao and Meskas, 2011).
In this paper, we simulate the flow using the single-relaxationtime (SRT) D3Q19 model employing the Bhatnagar-Gross-Krook (BGK) equation (Bhatnagar et al., 1954), broadly used due to its simplicity (Qian et al., 1992;Chen et al., 1992). This is a threedimensional model containing 19 velocities and particle distribution functions f a (Fig. 1). Based on the BGK model, update of the distribution function is as: Here, e a dt is a vector pointing to neighboring lattice points, dt the lattice time step, e a the discrete lattice velocity in direction aand sthe dimensionless lattice relaxation time/parameter. The right side of Eq. (1) here results in BGK collision operator: The term f ðeqÞ a is the equilibrium distribution function and is calculated as: where x a is the weighting factor related to the used LBM model, q is the lattice fluid density, c ¼ 1= ffiffiffi 3 p is the lattice speed of sound, u is the lattice fluid velocity. In the actual implementation of the model the distribution function (Eq. (2)) is solved in two steps, namely the collision step and the streaming step: f a ðx i þ e a dt; t þ dtÞ ¼f a ðx i ; t þ dtÞ ð 6Þ Eq. (5) represents the collision step that models interactions and calculates the updated values of the distribution function. In the streaming step (Eq. (6)), the distribution functions are streamed to the neighboring lattice points.

Fluid-particle interaction
Bounce-back boundary conditions where a fluid particle scatters back when reaching a boundary node are commonly used to implement a no-slip condition. We use a linearly interpolated bounce-back scheme for curved boundaries proposed by Bouzidi et al. (2001) and extended by Lallemand and Luo (2003). With the linear interpolation scheme for the no-slip boundary, the parallel code requires only one layer of ghost cells as opposed to two layers of a quadratic scheme, resulting in additional communication overhead. Given adequate resolution, the improvement in solution accuracy with quadratic interpolation is negligible (Pan et al., 2006;Kruggel-Emden et al., 2016). The particle surface can cross the connection between two nodes at arbitrary distances. This is termed as the fractional distance along the direction a and is given by: An example of a simple one-dimensional bounce-back boundary scheme is shown in Fig. 2 with three possible situations, summarized in two equations, depending on the value of q a : wheref a and f a are the distribution function before and after advection, and a represents the direction opposite of a.

Numerical setup
Random sphere packages with an artificially arranged channellike void region are studied to examine the effects of such heterogeneities on the hydrodynamic drag and the pressure drop. The simulations are set up within a framework of LB3D (Schmieschek et al., 2017), a parallel-processing implementation of LBM described in the former sections.
Prior to the LBM simulation, random configurations of spherical particles need to be generated. Spherical particles are initially placed in a face-centered-cubic (FCC) arrangement. In the FCC configuration, for every sphere there is a gap enclosed by six other spheres (octahedral) and two smaller gaps enclosed by four spheres (tetrahedral) making for twelve neighbors (coordination number 12, Fig. 3). The largest fraction of space occupied by solids that can be reached with such packing is p=ð3 ffiffiffi 2 p Þ ' 0:74048. After insertion, particles are randomly moved using a Monte-Carlo numerical method. Here, we do not aim for maximum packing density as the sphere radius must be reduced so there is enough space to move the spheres and find the random positions. Additionally, the distance between the nearest neighbor spheres must be at least one lattice space unit to resolve the flow around the sphere. With those limitations, generating the densest random packing requires a high domain resolution. To create a channellike heterogeneity, half of the particles are shifted in a positive direction on the z-axis. In this way, we can control the heterogeneity magnitude of interest, i.e. the characteristic width of the channel in dependence of the particle diameter. Here, we separate two quantities: solid volume fraction / p of the dense packed region and overall solid volume fraction of the simulation domain / that takes the void channel region into account (see Fig. 3). Approximating a large system, the code is set to obey periodic boundary conditions (PBC) in every Cartesian direction of the rectangular simulation domain. A fluid element transferred out of one domain side enters the appropriate position on the opposite side. For fluid-particle interactions, interpolated bounce-back is applied to obtain high accuracy solution. The Reynolds number is based on the magnitude of superficial fluid velocity u s and the particle diameter D p , as Re ¼ u s D p =m. The superficial velocity is a hypothetical value calculated as if the fluid was the only phase and it is related to the average fluid velocity u as u s ¼ ð1 À /Þu. The kinematic viscosity of the fluid m is related to the relaxation time parameter s by the expression m ¼ ðs À 0:5Þ=b 0 , where b 0 ¼ 3 for the 19 velocity vector model used here. The relaxation time needs to be within the stability limits 1:9 > s > 0:5, effectively establishing limits on the viscosity. The maximum lattice velocity is kept sufficiently low to avoid compressibility effects. With these constraints we limit the maximum solid volume fraction to / ¼ 0:55 and Reynolds number to Re ¼ 300. Previous studies employing the LB3D code investigated the simulation dependence on the grid resolution and domain size (Sanjeevi and Padding, 2017;Sanjeevi et al., 2018). The influence of grid resolution is stronger with increasing Re, so here we keep high resolution throughout all instances (360-576 grid cells in each direction with particle diameter D p = 63 grid cells, ensuring a precisely resolved flow).
Depending on the number of particles and packing density, the simulations are performed with different domain sizes. Through all simulations, the dense region is maintained at least twice as large as the channel in agreement with the assumption that heterogeneities occur at relatively small scale compared to the particle bulk. We increase the size of the simulation box with increasing channel width to prevent the periodic artifact from strongly influencing the flow. However, completely removing any periodicity effect would require a very large simulation box, so we approach the issue with a balance between computational costs and spatial resolution in mind (see Tables 2 and 3).

Drag force and pressure drop calculation
Fluid exerts two types of force on the particles: a force f d due to the friction between the particles and the fluid at the particle surface and a buoyancy force f b resulting from the average pressure gradient in the fluid. Therefore, the total force experienced by the particles is the sum of these forces f pÀf ¼ f d þ f b . In the literature, there is some ambiguity as to whether the contribution of the pressure gradient should be included in the drag force definition (Hill et al., 2001a;Van der Hoef et al., 2005;Tenneti et al., 2011;Tang et al., 2015).
In this work, the total force acting on the particle array f f Àp is obtained from the simulations and the mean force per particle is calculated as hf f Àp i ¼ f f Àp =n p . There is no buoyancy force in our direct numerical simulation setup, therefore f b ¼ 0. Hence, we present the results for the dimensionless drag force normalized by Stokes drag as F D ¼ hf f Àp i=f Stokes , where f Stokes ¼ 3pqmD p ju s jis the drag acting on a single isolated sphere in a fluid moving at a relative velocity equal to the superficial velocity.
The total fluid-particle interaction force is related to the pressure drop over the domain as:   ÀrP ¼ where V tot is the total volume of the system that contains n p particles. It follows that the correlation between the thus calculated dimensionless drag force F D and the overall pressure drop is:

Results
For each time step, the solver outputs the values of x-, y-and zforce and torque components on every solid particle (Schmieschek et al., 2017). We display the dimensionless fluid velocity field contours (u=u s ) obtained by resolving the flow around the spheres at three different packing densities (see Fig. 4). Velocity component values in the x-direction u are normalized by a superficial velocity u s which is known beforehand and provided as an input constant value (an example of input parameters shown in Table 2). Velocity field contours for all 12 configurations with varying heterogeneity magnitude are presented in Fig. 5 for a case of packing density inside porous region of / ¼ 0:4.
Values of the dimensionless drag forces from various performed simulations are presented in the Table 3. Simulations results are obtained with the steady state reached, performing time averaging when necessary. We plot the obtained values and compare the simulation data with predictions obtained by applying established correlation from the literature (Fig. 6, a similar approach for the comparison was used by Tenneti et al. (2011) and Tang et al. (2015)). The solid line with circles () represents our simulation data, dashed lines (Hill et al., 2001b), dotted lines (Beetstra et al., Table 3 Average dimensionless drag force FD for different configurations (i.e. channel width), where np is the number of particles, Dp the particle diameter in lattice units (particle resolution), Lz the size of computational domain in the z-direction in lattice units (domain resolution), Lc the characteristic width of the channel / overall solid volume fraction over the entire domain and / p is the solid volume fraction inside the dense packed regions.     (Tenneti et al., 2011) and dash-dot lines marked with points (Tang et al., 2015) show values predicted by the existing correlations.

Discussion
The calculated values of the drag force on the particles in a flow over random sphere packings fall into an expected range between other DNS-based drag predictions (Hill et al., 2001b;Tenneti et al., 2011;Tang et al., 2015). However, even the data reported by the existing literature is diverging; the results obtained by Tenneti et al. (2011) and Beetstra et al. (2007) differ up to 30% at higher Reynolds number. Tenneti et al. (2011) attribute the discrepancy to the constant resolution used by Beetstra et al. (2007) for the entire range of Re. However, Bogner et al. (2015) eliminated significant finite resolution effects and the results did not match precisely with the former studies. The results by Bogner et al. (2015) are not shown in the Fig. 6 because their investigation is tied to lower solid volume fractions (/ ¼ 0:1 À 0:35). This shows the need for further development of accurate simulation methods in the future even for random homogeneous packings, especially in case of high Re. We suspect a more accurate interpolation scheme is needed to capture the effects of the thinner (momentum) boundary layers.
In this work, we focus on the heterogeneous structures. The inclusion of a channel-like heterogeneity to the configuration leads to a significant drag reduction when compared to the flow over a random sphere array with an equal porosity (Fig. 6, right). The obtained forces are lower as we further increase the magnitude of the heterogeneity (Fig. 7). Fluid flows through the channel region with low flow resistance which results in significantly higher velocities relative to those observed in dense regions. The difference between dimensionless velocity through two regions can be observed in Figs. 8 and 9. We observe the same behavior for the entire investigated range of Reynolds numbers and porosity but it is especially prominent for low Re and high heterogeneity magnitude.
In a large-scale granular system, the channel formation can develop on a scale smaller than the grid resolution of the unresolved simulations. Consequently, the channels might be overlooked leading to the application of inadequate drag law, where the particle distribution is assumed as being homogeneous. This inevitably would result in an overestimation of the pressure drop. Our proposed analytical description of overall fluid-solid momentum exchange might open the path towards the development of Fig. 7. Dimensionless drag force dependence on the increase of Reynolds numbers for different normalized channel widths Lc=Dp. Overall solid volume fraction / reduces with a channel width increase.   9. Averaged dimensionless radial velocity profiles for a particle bed region (left) and channel/crack region (right). For lower Re, the flow profile is close to parabolic with the expected uavg ¼ 2=3 umax, where uavg is calculated with an average velocity in the dense part up deducted. With higher Re, one can observe more evenly distributed flow through the channel width the average value higher than parabolic 2=3 umax, indicating transient flow regime where flow is not completely dominated by viscous effect. It should be noted that the velocity scales are different for each profile plot. specific sub-grid models for unresolved CFD-DEM simulations (i.e. for the application in dense flows of cohesive powders).

Overall pressure drop
For heterogeneous particle arrangements containing a channellike void region, the simulations show a significant decrease in overall pressure drop. Considering the immense impact of the channel, the heterogeneity intensity is examined comparing the velocities across the domain. For the majority of simulated cases, the velocities through dense particle areas are significantly lower than the overall superficial velocity u s . Furthermore, we aim to accompany the numerical results with a semi-analytical approach, looking into existing correlations for pressure drop in dense particle beds and in channels (laminar and turbulent).
In our analysis, we assume the flow is under significantly different conditions and regimes through the respective dense and channel regions. One can observe that this assumption holds better for configurations with a wider channel and lower Reynolds number flows (Figs. 8 and 9). However, the central channel is very clear in the average axial velocity profile for every simulated case. The flow through particle bed is estimated as uniform with an average superficial velocity u p and the channel average velocity as u c . In this manner, total flow flux through the domain is divided into two parts as: where S is total cross-sectional surface while S p and S c are crosssectional surfaces for porous and channel regions respectively. In our configuration, therefore: where L y and L z are the domain sizes in appropriate directions orthogonal to the flow and L c is the characteristic width of the channel. Because the channel is parallel to the flow, the pressure drop through the two distinct regions must be equal: dp dx Together with Eq. (12), Eq. (13) creates the basis for a simple calculation of a pressure drop through investigated configurations. Note that in the presence of a channel, the flow velocities in the dense packed regions tend to be low. Using well-known equations, it is possible to estimate the pressure drop and fluid velocity for low Reynolds number flow through dense bed of spheres where the viscous forces of the fluid dominate. Various established works on the fluid flow through packed beds of spheres exist in the literature (Carman, 1937;Ergun, 1952;Wen and Yu, 1966). Average velocity in the dense region is extracted from the simulation data and local particle Reynolds number Re p can be calculated for these regions only. Initially, we estimate the pressure drop for the flow through a packed bed of solids using the Kozeny-Carman equation.

DP
Conversely, the flow in the channel part is considerably faster. The empirical Darcy-Weisbach equation relates the pressure drop due to friction along a given length of the channel to the average velocity of the fluid flow. It is valid for fully developed, steady state and incompressible flow. where D ¼ 2L c is a hydraulic diameter and k f the dimensionless friction factor, which depends on the characteristics of the flow. In our setup, the channel geometry is equivalent to a pair of rough infinite (considering PBC) parallel plates. For laminar flow between two smooth parallel plates the friction factor is inversely proportional to the Reynolds number as k f ¼ 96=Re, matching the Hagen-Poiseuille equation. Plane Poiseuille flow is a well-known theoretical case, which is analytically derived from the Navier-Stokes equation, that can serve as an initial point for the pressure drop calculation through the channel region in the laminar regime.
The channel size can be written as a function of overall solid volume fraction and solid volume fraction of the dense packed region: Here we introduce the heterogeneity index f as the ratio between the volume-averaged velocity through the channel u c and the volume-averaged superficial velocity through the particle bed u p , representing an indicator of the heterogeneity influence on the flow profile: Calculating the pressure drop for laminar flow through two different regions with Eq. (16) and (14), it is possible to calculate the expected ratio of these distinct velocities: The velocity ratio calculated in this fashion (Eq. (19)) is independent of Re. A comparison between the velocity ratios obtained from resolved simulations and those calculated from the pressure drop correlations for the laminar flow is shown in Fig. 10. There is a good agreement between the approximation and numerical results for Re ¼ 30. The channel influence is stronger for low Re and larger channel size. This is expected because for higher Re the channel flow is in a transitional regime, while the packed bed flow remains laminar. The pressure drops in the channel and dense region are expected to remain equal to each other; however, the friction factor of the channel increases faster, and therefore the equal pressure drop assumption leads to a change in the ratio of channel velocity to packed bed velocity.
Consequently, it is necessary to adopt a different approach for higher Re. Fully developed turbulent flow in a rough duct can be defined by the Colebrook expression. However, due to the implicit nature of this correlation, an alternative form was proposed by Swamee and Jain (1976): where =D is relative roughness and Reynolds number is calculated for the channel based on the hydraulic diameter D ¼ 2L c , i.e. as Re c ¼ 2u c L c =m. We have estimated the relative roughness with average height of surface irregularities of one particle radius ( ¼ 0:5D p ). The flow in our simulations is never fully turbulent but falls into the transient flow regime. In general, the results in the transient regime are harder to predict. Nevertheless for higher Re, the calculated channel pressure drop and the dimensionless friction factor (with Eqs. (15) and (20), respectively) are closer to the numerical results as the flow is closer to truly turbulent. As it was previously indicated by the estimation of heterogeneity index, the laminar friction factor is matching the simulations for low Re (Fig. 11, right).
The pressure drop through the packed bed is easier to predict as there we do not deal with the instabilities of the transient flow. The appropriate literature correlation for the dimensionless drag force chosen here is proposed by Tenneti et al. (2011) To avoid the heterogeneity influence, modified bed superficial velocity u 0 p is obtained from the numerical results taking into account only dense regions sufficiently distant from the channel. Therefore, the correlation (Eq. (21)) can be written as a function of the bed solid volume fraction / p and particle Reynolds number, calculated as Re 0 p ¼ u 0 p D p =m, i.e. ÀrP ¼ f ðRe 0 p ; / p Þ. The pressure drop calculated in this manner predicts the numerical results very well (Fig. 11, left). To observe the change in the pressure drop with Re and channel width (Table 4), the obtained values are normalized by the pressure drop calculated from the Eq. (9) with F D ¼ 1. For the normalization, we take the overall solid volume fraction into account, as if the particles are homogeneously distributed through the simulation domain.

Conclusion and outlook
In this work, we studied fluid flow through static particle arrangements by means of highly resolved lattice-Boltzmann simulations and evaluated fluid-solid momentum exchange and pressure drop. In case of randomly arranged (i.e. homogeneously distributed) spheres our numerical results agree very well with previous studies (Hill et al., 2001b;Tenneti et al., 2011;Tang et al., 2015).
In the case of heterogeneous particle arrangements comprising a channel-like void region in between two dense bulk regions, our simulations indicate a dramatic decrease in overall pressure drop even for channel widths of only one particle diameter. We substantiated this main finding by a series of numerical simulations varying superficial Reynolds number and the dense bulk volume fraction. Table 4 Comparison between the dimensionless pressure drop ÀrP Ã LB3D obtained from the numerical simulation with the one calculated using Darcy-Weisbach equation ÀrP Ã DW . The friction factor is calculated combining the laminar (k f ¼ 96=Rec, for Re = 30) and , for Re = 100, 300) expressions for the dimensionless friction factor. We further supplemented these purely numerical findings by a combined analytical approach, by merging existing correlations for pressure drop in homogeneous particle beds and in channels. Relevant correlations were selected according to the superficial Reynolds number, as the simulations range from laminar to transitional flow. For laminar channel flow, our analytical predictions are in good agreement with the numerical results. However, friction and pressure drop in the transitional regime are characterized by instabilities and are difficult to predict. Nonetheless, as we approach the turbulent flow the results are in closer agreement with appropriate channel flow pressure correlations.
As the main message, this work underlines the significance of channel/crack formation on fluid-solid momentum exchange and overall pressure drop in dense particle arrangements. At the same time, such sub-grid heterogeneities are commonly neglected in unresolved CFD-DEM simulations of fluid flow through dense particle assemblies. In this regard, our proposed analytical description of overall fluid-solid momentum exchange in particle arrangements with channel-like heterogeneities might pave the way towards the development of specific sub-grid drag models in unresolved CFD-DEM simulations of cohesive powders and interlocking non-spherical particles.