Lift-off of multiple particles in a narrow channel

Abstract In this work, we perform simulations of particle laden flow in a wide and long narrow channel in a Newtonian fluid. Simulations are performed for mono-sized and equal density spheres with varying Archimedes and Reynolds number. In the simulations, different phases of particle transport - rolling, saltation and suspension are observed. During simulations the average bed height is monitored and its steady state value is used for proposing a correlation between solids volume fraction ( ϕ ), shear Reynolds number ( Re s ) and Archimedes number (Ar). This correlation is used to predict the critical shear Reynolds number for particle lift-off and a condition at which particles would occupy the whole channel at a given Archimedes number. The value of this critical Reynolds number is compared with the critical Reynolds number for single particle lift-off.


Introduction
Particle laden flows are encountered in many industrial as well as in natural processes such as biological and environmental flows, e.g. proppant transport in hydraulic fracturing. In these flows, the motion of particles is influenced by hydrodynamic effects, interparticle and particle-wall interactions, and gravity. The distribution of particles in a flow channel results from a combined effect of all above-mentioned phenomena and is therefore difficult to predict. An efficient numerical model that takes into account all these phenomena can be very helpful to understand and optimize processes involving particle transport.
Fluid-solid multiphase flow can be simulated using either an Euler-Euler (E-E) or an Euler-Lagrange (E-L) approach. In an E-E approach, both the fluid and solid phase are treated as interpenetrating continua coupled via momentum transfer (Kuipers et al., 1992;Gidaspow, 1994). This approach is commonly referred as the Two Fluid Model (TFM). In the E-L approach the motion of individual particles or parcels (i.e. group of particles) are simulated in a Lagrangian manner. The challenge of TFM is an accurate rheological and hydrodynamic description of the solid phase. Such a description is often obtained using the kinetic theory of granular flows. An advantage of the E-E compared to the E-L approach is a considerable saving of computational cost, which facilitates simulating large scale phenomena. Basu et al. (2016) performed simulations of particle transport in a channel using an E-E approach for mimicking proppant flow in a representative fracture. Although a satisfactory qualitative match with the E-L results for particle concentration was reported, large quantitative deviations were observed.
One of the E-L approaches is the Discrete Particle Model (DPM) where the motion of individual particles is tracked and the fluid velocity field is solved on a unresolved grid (Deen et al., 2007). DPM accounts for fluid-particle interactions using closure relations for momentum exchange and for particle-particle/particle-wall interactions using a collision model. This approach has been used for proppant transport studies by Blyton (2016), Zhang et al. (2017), Baldini et al. (2018). Tomac and Gutierrez (2014) also include a model for lubrication forces to study its effect on particle transport in a channel. DPM gives a more accurate representation of the system compared to E-E method, however at the expense of higher computation costs. Snider (2001) proposed a Multi Phase-Particle in Cell (MP-PIC) approach which is more accurate than the E-E approach, while computationally less expensive than DPM. In MP-PIC, a parcel is tracked as opposed to individual particles in DPM. Tsai et al. (2012) used this approach to study proppant settling in shale gas production.
The above-mentioned E-L approaches show distinctive advantages compared to the E-E approach. However, the momentum exchange between fluid and solid phases has to be modelled and still contains uncertainty due to the use of a closure relation for the drag force. There is an additional error due to the unavailability of a closure relation for the lift force in particle suspensions. This could be particularly important for particle transport in narrow channels as the effect of lift can be significant. Moreover, rotation of particles also plays an important role during their transport, but due to the unavailability of reliable correlations for the hydrodynamic torque its effect on a particle is usually ignored. These challenges can be mitigated by resolving the flow field around a particle and integrating the hydrodynamic forces and torque over the particle surface. This approach is known as Direct Numerical Simulation (DNS) and has been used widely to study transport of buoyant particles. A major advantage of this method is that only a minimal amount of physics needs to be modeled. An important unresolved force is the lubrication force. Completely resolving it is computationally too intensive because its range of influence continues to small scales that are three to four orders of magnitude smaller than the particle size.
Based on the flow conditions, particles either roll or get suspended by lift forces. This phenomenon is referred as fluidization by lift in the literature. DNS studies of fluidization by lift are performed mainly in the context of free surface flows (Kidanemariam and Uhlmann, 2014;Kempe et al., 2014;Vowinckel et al., 2014;Jain et al., Feb 2020). There are very few DNS studies focusing on fluidization by lift in narrow channels (Patankar et al., 2001;Xiaoqi, 2018). Although similar hydrodynamic and inter-particle interactions influence these systems, some differences will exist due to the different velocity profiles and range of Reynolds numbers. In this paper, we investigate a suspension of particles for different Archimedes and Reynolds numbers in a planar flow. The flow profiles and particle distributions are observed to characterize the different stages in particle suspension by channel flow. The average bed height is measured as a function of time and the equilibrium height at steady state is used to propose a correlation between the shear Reynolds number, Archimedes number and solids volume fraction.

Numerical setup
In this paper, simulations of transport of hundreds of particles in a narrow channel are performed. These simulations are particle resolved, with a coupling between particle and fluid motion implemented by means of a ghost-cell immersed boundary method. A detailed description of the simulation methodology is provided by Maitri et al. (2018).
To mimic a long and wide channel, periodic boundary conditions are implemented in the x-and z-direction whereas no-slip boundary conditions are applied at the upper and lower walls at y ¼ 0 and y ¼ L y . Flow is driven by applying a pressure gradient in x-direction. Initially, particles are placed in a simple cubic configuration as shown in Fig. 1. A parabolic flow profile is imposed in the channel as an initial condition. Simulations are first conducted with static particles such that the flow profile develops. Once a steady flow profile is obtained, particles are released to freely rotate and translate according to hydrodynamic and interparticle forces and torques. Inter-particle interactions are accounted for using lubrication and collision models. Gravity acts in the negative-y direction.
The channel dimensions are L x Â L y Â L z ¼ 20 d p Â 8 d p Â 6 d p . A uniform grid of D ¼ 5 Â 10 À5 m is used and 1 mm particles are resolved with 20 grid cells across their diameter. The Newtonian fluid density is 1000 kg=m 3 and the viscosity is 0:001Pa s. In the simulations 480 particles are arranged in 4 layers, where each layer consists of 120 (20 Â 6) particles. Unresolved lubrication force for particle-article and particle-wall interaction is modeled using an approach of Breugem (2010). Particles in current simulations have a normal coefficient of restitution e n ¼ 0:97 and a tangential coefficient of restitution e t ¼ 0:33. The friction coefficient is 0.1 for Fig. 1. Schematic of the simulation setup and initial velocity on a cross section for Res ¼ 50.

Fig. 2.
Particle distribution and velocity profiles in a cross-section at different times for Ar = 9.81 and Res = 20. Saltation phenomena is observed for few particles initially, but particles remain sedimented after the quasi-steady state is achieved.
particle-particle collisions and 0.2 for particle-wall collisions. These collisional properties correspond to a system of glass beads and alumina wall. A shear Reynolds number is defined based on the wall shear rate for an empty channel, Simulations are performed for multiple shear Reynolds and Archimedes numbers, The Reynolds number is varied by changing the applied pressure gradient whereas varying the solid density gives different Archimedes numbers (9.81, 49.05, 98.10). The values of Archimedes number are chosen to represent the typical Archimedes number in the hydraulic fracturing application. To keep the simulation setup the same, the particle diameter is kept constant and the Archimedes number is varied by varying the particle density. In a real application the situation would probably be the opposite. Particle diameters might vary at fixed material properties. Keeping fluid properties same and using the typical particle density of 2500 kg/ m 3 , chosen Archimedes numbers in this paper (9.81, 49.05, 98.10) give particle diameters of (87, 149, 188) lm, respectively. Note, however, that when varying particle diameter also varies other dimensionless quantities such as the ratio between channel dimensions and particle size, which in our simulations are kept constant. For the considered Archimedes numbers, multiple Reynolds numbers are chosen such that at the maximum Reynolds number, particles are fully suspended and occupy the whole channel.
During the simulations, where particles are freely rotating and translating, the average half 'height' of the particle bed is approximated as: In this equation, y i is the distance of the centre-of-mass of particle i from the lower wall. This definition of h is consistent with that used by Patankar et al. (2001). It is the centre-of-mass y-position of the particle bed. If a particle bed has a uniform particle distribution the height of the surface of the bed is 2 h. Also in case of a nonuniform distribution we will use 2 h as an approximation of the bed height and use it to compute the average solids volume fraction as In the current simulations, the maximum solids volume fraction is achieved when particles do not suspend and pack on the bottom of the channel. The minimum solids volume fraction is obtained when particles occupy the whole channel. Hence, the results obtained with the current simulations are valid for the range of solids fraction between p=12 (0.2618) and a maximum packing density of approximately 0.64. In this paper, time is made dimensionless using the characteristic time equal to the inverse of the shear rate at the wall (in absence of particles) at a steady state (see Eq. 1). The same shear rate at the wall was used to compute the shear Reynolds number in Eq. 1.

Results
In this section, we will present the results of transport of particles in a narrow channel. Transport of particles by channel flow occurs in three states: rolling, saltation and suspension (see . Rolling is observed when the Reynolds number is small and the corresponding lift force is insufficient to lift any of the particles. In this case, particles roll over other particles or along the wall. In saltation, the lift force on a particle overcomes the immersed weight of the particle and hence the particle is lifted off. However, after lift-off the magnitude of the lift force is not sufficient to bal- ance gravity and the particle falls back on the bed. When the Reynolds number is sufficiently high, a suspension state is observed where multiple particles remain suspended for longer times. The objective of the simulations is to understand different stages of suspension of particles by lift forces and subsequently find a correlation between solids volume fraction, shear Reynolds number and Archimedes number. An accurate correlation can also be used to predict an the average bed height as a function of the shear Reynolds and Archimedes numbers. Fig. 2 shows the particle distribution and velocity profile on a cross-section plane at various times for Ar ¼ 9:81; Re s ¼ 20. Initially one or two isolated particles are lifted. However, at longer times these particles do not stay suspended due to an insufficient lift force and they fall back on the particle bed. At steady state, the particle bed moves with the flow without suspension of any particle. The average height of the bed at the steady state is less than the initial average height. For this case, two phenomena are observed: rolling and saltation. The pressure profile, relative to the imposed macroscopic linear pressure field, before releasing the particles shows high pressures in front and low pressures behind all particles in the top layer (see Fig. 5). Due to very small velocities in the interior of the particle bed at the start of the simulation, significant pressure variations are not observed below the top layer. The pressure differences over the top-layer particles provides the form drag that sets the layer in motion. At later times also the particles in lower layers move. In Fig. 6a profiles of the axial liquid velocities are shown. At t Ã ¼ 0 there is a parabolic profile in the region between the bed surface and the upper wall, with (nearly) zero velocity inside the particle bed. At later times there is sizable axial velocity inside the particle bed. This is due to the particle motion that drags along the liquid. In the later stages particles remain layered and roll over each other. The liquid velocity profile in the bed increases with increasing distance from the lower wall. The particles stay arranged in four layers and sharper steps occur at the more porous interstitial layers between particle layers. In Fig. 6a it can be observed that within a particle-layer the axial velocity increases linearly. This linear increase within the particle layer is due to the rotation of particles with no-slip boundary conditions with the fluid.
When the shear Reynolds number is increased to 50 for Ar = 9.81 (see Fig. 3), initial observations are similar to that of the Re s ¼ 20 case. The particles start in a saltation type of motion, where single particles are lifted, disturb the flow, collide with other particles and thereby lift these particles. However, in the final state the average lift force acting on particles is higher than for Re s ¼ 20 and particles remain suspended, with an average bed height higher than the initial height. A notable observation for this case is that the bottom layer of particles does not get suspended and particles in this layer roll over the bottom wall. This can also be observed from the velocity profile shown in Fig. 6b where a slowly translating and rotating particle layer is observed at the bottom wall. Another important observation is that the initial parabolic profile in the particle-free region is adjusted and the velocity profile becomes flatter in the centre region indicating lower shear rates.
When the shear Reynolds number is further increased to 200, the particle bed suspends much quicker and particles from the top layer reach the top wall. However, later these particles are pushed downwards again. The reason for this is the interaction of particles with the top wall as well as the insufficient hydrodynamic force to balance the gravitational force. Simultaneously, the lower layer of particles is pushed upwards with higher magnitude and therefore, a monotonic increase in the average bed expansion is observed. Unlike the previous case of Re s ¼ 50, the bottom layer of particles is also suspended and particles get more concentrated at the centre of the channel. Fig. 6c shows the velocity profile at different times and a more symmetric velocity profile is observed at the steady state. A flatter velocity profile is observed in the centre, which suggests shear-thinning behavior fluids. We have also performed this simulation for a 3 times longer domain and a similar bed evolution is observed. Fig. 7 shows an example for a larger Archimedes number, namely, Ar ¼ 49:05. Compared to the smaller Archimedes number, heavier particles need a larger Reynolds number to get suspended because they possess more inertia. This means that once the particles pick up momentum it takes a longer time to loose it again. Initially the particles pick up speed in a saltation motion. Particles with larger momentum then move to the top of the channel without loosing speed. This results in a quick bed expansion. Due to  6. X-velocity profiles as a function of height for Ar ¼ 9:81. Initially, a small magnitude of velocity is observed in the particle bed and a parabolic profile is observed above the particle bed. The particle bed is pushed by the flow, where particles roll over each other which can be observed as the linear part in the velocity profile. For the highest Res of 200, a flat velocity profile is observed in the centre region which is similar to the velocity profile observed in shear-thinning fluids. frictional interaction with the top wall kinetic energy is partly lost. This causes the bed height to decrease again. Fig. 8a-c show the average bed height as a function of time. For Ar ¼ 9:81, a monotonic increase or decrease of the bed height is observed. Similar observations are also made for Ar ¼ 49:05 and Re s = 100, 200, 300. For Re s > 300, the bed height exhibits an overshoot as a function of time. It initially increases and then slowly decreases until reaching a steady state value. This trend is similar for Ar = 98.10. Note that for Ar ¼ 98:10 with Re s = 100 and 200 one would expect a gradual increase of bed height with time. The fact that an overshoot is observed is because the starting point for these simulations is based on the velocity, pressure and particle information from Re s = 300 simulations at around t = 0.4 s.

Lift-off correlation
At stationary state the bed height 2 h is computed using Eq. 3. This height depends on Re s and Ar. From the height an average solids volume fraction in the bed, /, is computed using Eq. 4. In this section, the approach to develop a correlation between /; Re s and Ar is described. There are other dimensionless variables that should enter such a correlation, namely, the aspect ratio (d p =L y ) and the global solids volume fraction (/ min ) inside the channel. However, these parameters were kept constant for all the simulations and hence their dependence is not investigated here. Our assumption is that the dependence on d p =L y and / min is weak. Fig. 7. Particle distributions and velocity profiles on a cross-section at different times for Ar = 49.05, Res = 300. In this case, all particle layers are suspended. Top layer of particles are fully suspended as they rise till the top wall but eventually descend to lower heights.
First, a fit between steady state values of the average solids volume fraction and shear Reynolds number, of the functional form: is performed separately for each Archimedes number. In case of the rolling state, bias arises if the considered shear Reynolds number is lower than the critical Re s and for the case of suspension state. Bias also arises when the Re s is higher than the Reynolds number at which the bed is fully suspended. Therefore, data points where the particle bed is fully sedimented or fully suspended are not considered.
Values of the coefficient a and the power, m, for different Archimedes number are tabulated in Table 1. A match between correlation and simulation values is shown in Fig. 9a-c. The value of m does not vary significantly for different Archimedes number, however the value of a is sensitive to the Archimedes number. Therefore we combine all data with the aim of constructing single correlation that is valid for the considered range of /; Ar and Re s . A refit with also a power-law dependence on the Archimedes number gives: The curve-fitting is performed using non-linear least squares using Levenberg-Marquardt algorithm from scipy python package. A parity plot of simulated versus predicted values of / is shown in Fig. 10. It shows a good match except for a few outliers. The range of validity for the proposed correlation is 9:81 < Ar < 98:10; 20 < Re s < 2000 and 0:2618 < / < 0:64.
The correlation Eq. 5 can have several usage. For example, for a certain bed-weight it can be used to estimate the pressure drop needed to suspend particles to a certain average solids volume fraction. By substituting the values of maximum and minimum solids volume fraction for a given Ar, we can obtain the corresponding critical shear Reynolds number at which the particle bed starts to suspend and when the particle bed is fully suspended, respectively. In Table 2 these values are compared with the values for the critical Reynolds number required for the lift-off of a single particle. This value is obtained by simulating a single particle in the channel. This particle was constrained to translate and rotate in a horizontal plane near the wall. In this simulation the average lift force on the particle was determined. Lift-off is the state where lift-force equals the difference between weight of the particle and buoyancy force. The comparison shows that the critical shear Reynolds number for the lift-off of many particles is in the same range as that of single particle lift-off condition. This finding is different for circular particles where the critical Re s for multiparticle systems is lower by an order of magnitude compared to the critical Re s for the single particle (Patankar et al. (2001)). The Table 2 suggests a trend where, with increasing Archimedes number, particles in the upper layer require a lower lift-off velocities than single particles at a wall. However, to claim such a trend more than three data points are needed. Patankar et al. (2001) followed a similar approach for 2D simulations of circular particles. In their work the fluid volume fraction was used instead of solids volume fraction. We have tried their approach with fluid volume fraction, however, the power law correlation fitted better when using solids volume fraction. Similarities of our results compared to the circular particle findings are: 1. Abrupt bed expansion and then gradual decrease of the average bed height is observed at higher Reynolds number. 2. Three states of particle transport in the fluid are observed: rolling, saltation and suspension.  Differences compared to the circular particle findings include: 1. For spherical particles, critical Reynolds number for the expansion of the bed is of the same order of magnitude compared to the critical Reynolds number for the single particle, whereas for circular particle these two values differ significantly. 2. At steady state, the average bed height is lower than the height of the channel. 3. The lift force has a stronger dependence on shear Reynolds number in spherical particles compared to circular particles.

Conclusions
In this paper, we performed direct numerical simulations of transport of particles in a Newtonian fluid in a long and wide narrow channel. The simulations are performed for multiple Reynolds numbers for three Archimedes numbers that were varied by changing the solid density. For low shear Reynolds number, a rolling motion is observed where particles remain sedimented and roll over other particles or the bottom wall. For intermediate Reynolds numbers, a saltation phenomenon is observed where few particles from the top layer are lifted-off but due to insufficient lift force they get sedimented. At higher Reynolds number, sufficiently large numbers of particles lift off and stay suspended as the hydrodynamic lift force balances the particles' immersed weight. In the suspended state at lower Reynolds numbers it is found that although first three layers of particle are lifted-off, particles in the bottom roll on the wall. Suspension of all particles is observed when Reynolds number is increased further. Monotonic increase in the average height is observed up to a certain Reynolds number after which the particle bed is found to rise abruptly and afterwards slowly descends.
Based on the numerical simulations, a correlation is obtained between solids volume fraction (/), shear Reynolds number (Re s ) and Archimedes number (Ar). This correlation is based on a power law dependence between these parameters and it can be used to find the critical value of shear Reynolds number required for the lift-off and also the condition at which particles occupy the whole channel. The critical Re s required for lift-off is found to be comparable to the critical Re s for lifting a single particle.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements
This work is supported by the programme 'Computational Sciences for Energy Research (CSER)' of the Foundation for Fundamental Research on Matter (FOM) which is now part of the Netherlands Organisation for Scientific Research Institutes (NWO-I). This research is also co-financed by Shell Global Solutions International