CFD Investigation into the Effects of Surrounding Particle Location on the Drag Coefﬁcient

: In the simulation of dilute gas-solid ﬂows such as those seen in many industrial applications, the Lagrangian Particle Tracking method is used to track packets of individual particles through a converged ﬂuid ﬁeld. In the tracking of these particles, the most dominant forces acting upon the particles are those of gravity and drag. In order to accurately predict particle motion, the de-termination of the aforementioned forces become of the upmost importance, and hence an improved drag force formula was developed to incorporate the effects of particle concentration and particle Reynolds number. The present CFD study examines the individual effects of particles located both perpendicular and parallel to the ﬂow direction, as well as the effect of a particle entrain within an inﬁnite matrix of evenly distributed particles. Results show that neighbouring particles perpendicular to the ﬂow (Model 2) have an effect of increasing the drag force at close separation distances, but this becomes negligible between 5–10 particle diameters depending on particle Reynolds number (Re p ). When entrained in an inﬁnite line of particles co-aligned with the ﬂow (Model 1), the drag force is remarkably reduced at close separation distances and increases as the distance increases. The results of the inﬁnite matrix of particles (Model 3) show that, although not apparent in the individual model, the effect of side particles is experienced many particle diameters downstream.


Introduction
Dilute gas-solid flows are of considerable importance in many technical and industrial processes to efficiently transport solid particles that have a wide range of sizes (few µm to few mm) and density [1,2]. The application includes, but not limited to pneumatic conveying, fluidised beds, vertical risers, classifiers, cyclones, and flow mixing devices [3,4]. Pneumatic conveying, due to its many advantages such as simplicity and flexibility in operation, environmental compliance, and inherent safety, is widely used in the chemical, food processing, and cement industries and also in order to transport pulverised coal in thermal power plants [5]. This wide application has led to extensive research on pneumatic conveying of solids through the different pipe elements [6,7]. Despite its wide application, the design of the pipe network for pneumatic conveying is still a challenge because ostensibly small changes in the system or product frequently cause significant changes in the system performance or design [8]. Furthermore, during the flowing, these fine particles have not only the interactions with gas phase but also the collisions with each other and the wall of the pipes [9,10]. Therefore, the flow phenomena of gas-solid two-phase flows are very complex. Computational Fluid Dynamics (CFD) is usually used for modelling the flow in these systems for a number reasons ranging from the design stage to monitoring flows where experimental measurements are unavailable.
The CFD modelling relies heavily on a strong knowledge base of the fluid flow under consideration, and as such, the behaviour of particle motion is of the utmost importance when trying to simulate gas particle flows [11]. Over the years, many researchers have looked into the CFD modelling of gas-particle flows and the degree of accuracy has greatly improved with advances in turbulence modelling and computing power [12][13][14]. In many applications, the importance of the drag force in accurately predicting particle motion cannot be understated [15,16]. When looking at gas driven flows, the drag force is responsible for the acceleration of the particles. In CFD simulations, especially using the Lagrangian framework [17], the most common approach to determine the drag force on a particle is using the standard drag coefficient curve, which is based on experimental studies of a sphere in unbounded fluid flow. The most commonly accepted approximation of this curve is given by the following equation [18]: As it can be seen, the above equation is solely a function of the local particle Reynolds number, and as such, discounts other effects that may affect the drag on a particle. In a situation whereby particles are relatively spread out, this assumption can be correctly employed to accurately predict particle motion, but in many industrial flows, it is impossible to assume that the distribution will be such that particles will not interact and in turn will affect each other's motion.
A number of researchers have tried to measure the influence on the drag force of a particle in the presence of other particles. Liang, Hong and Fan [19] proved that altering the position of surrounding particles experienced drastic changes in drag coefficient. A configuration of three-coaligned particles led to a reduction of drag experienced by the centre particle compared with that of the leading particle at a separation distance of 2-3 d p . Cheng and Papanicolaou [20] calculated the analytical force on an array of particles at low Reynolds numbers and volume fraction. The analytical results showed good comparison with the phenomenological results available at the time of this work. Kim, Elghobashi and Sirignano [21] solved the full Navier-Stokes equations for spherical particle motion at a range of Reynolds numbers and particle-to-fluid density ratios. The full Navier-Stokes solution showed considerable differences to some of the more commonly used particle motion formulas and resulted in the authors proposing their own new particle motion formula. Zhang and Fan [16] proposed a new semi analytical expression for the drag force of an interactive particle due to wake effect. This work was based on the experimental findings of Liang, Hong and Fan [19] looking at separation distances up to 7 particle diameters and particle Reynolds number ranging from 54-154. Zhang and Fan [22] used the above work as a basis to predict the rise of interactive bubbles in liquids. The work showed that the new drag model and the inclusion of both the added mass and basset force provided the best agreement with the available experimental data [23].
The work of Liang, Hong and Fan [19] looked at three particles co-aligned and at separation distances up to 7 particle diameters, which corresponds to a particle volume fraction of approximately 10 −3 . In order to extend this work to investigate the effects at lower concentration values, 10 −4 , which is common in many dilute phase industrial flows, a full CFD study of this work was undertaken. The aim of this study was to produce a new drag force relation that considered the influence of particle concentration as well as particle Reynolds number. The new relation needed to be a generic equation which relied on the particle concentration or particle volume fraction, α, because the Lagrangian particle tracking method tracks individual representative isolated particles without considering the position of surrounding particles. An assumption was made that because no particle position is actually known, the particles in any given cell are evenly distributed in a cubic formation, see Figure 1. Utilising this assumption, the volume fraction can be transformed into particles separated by a uniform distance in all directions. 3 (2) where L is the distance between particle centres and l is the distance between particle centres as a ratio of particle diameter.
cubic formation, see Figure 1. Utilising this assumption, the volume fraction can be transformed into particles separated by a uniform distance in all directions.  π π π ε = = = = = (2) where L is the distance between particle centres and l is the distance between particle centres as a ratio of particle diameter.

Mathematical Models
This section covers the numerical methods used to solve the motion of both the continuous fluid and the dispersed particle motion. The continuous fluid phase is calculated by solving the governing equations. The dispersed phase is solved in the Lagrangian framework using Newton's law to describe individual representative particle motion. In present work, CFD technique is used as a tool to develop a new particle drag force equation that considered the influence of particle concentration (α) as well as particle Reynolds number. This paper used the CFX-5 package to replicate experimental work of Liang, Hong and Fan [19].

Governing Equations
A general form of the transport equation used in CFD model is shown in Equation (3). This equation describes the instantaneous transport of any variable quantity φ, which may be mass, momentum, heat, or a scalar quantity such as turbulent kinetic energy.
where is the fluid density, x is the distance in the ith, is the velocity in the ith,  is the diffusion coefficient of the variable φ and is a source or sink term for the variable φ.
For a fluid with constant density, the Navier-Stokes system of equations consist of the continuity Equation (4), derived from the principle of conservation of mass, and three equations derived from the principle of the conservation of momentum (5).

Mathematical Models
This section covers the numerical methods used to solve the motion of both the continuous fluid and the dispersed particle motion. The continuous fluid phase is calculated by solving the governing equations. The dispersed phase is solved in the Lagrangian framework using Newton's law to describe individual representative particle motion. In present work, CFD technique is used as a tool to develop a new particle drag force equation that considered the influence of particle concentration (α) as well as particle Reynolds number. This paper used the CFX-5 package to replicate experimental work of Liang, Hong and Fan [19].

Governing Equations
A general form of the transport equation used in CFD model is shown in Equation (3). This equation describes the instantaneous transport of any variable quantity φ, which may be mass, momentum, heat, or a scalar quantity such as turbulent kinetic energy.
where ρ is the fluid density, x is the distance in the ith, u is the velocity in the ith, Γ is the diffusion coefficient of the variable φ and S is a source or sink term for the variable φ. For a fluid with constant density, the Navier-Stokes system of equations consist of the continuity Equation (4), derived from the principle of conservation of mass, and three equations derived from the principle of the conservation of momentum (5).
∂ρ ∂t where S ij is the strain tensor, defined as:

Turbulence Model
In the present study, the Reynolds Averaged Navier-Stokes (RANS) equation was used to describe turbulent flows. The RANS can be used with approximations based on knowledge of the properties of flow turbulence to give approximate time-averaged solutions to the Navier-Stokes equations. The RANS equation can be written in the following form: which only differs from the laminar form by the term ρu i u j . The quantity is known as the Reynolds stress tensor and has six independent components. As the Reynolds stress tensor contains six unknown quantities combined with the unknown velocity and pressure terms, we now have 10 unknowns and only 4 equations (3 momentum and the continuity) to solve. In order to analytically solve these equations, a form of turbulence closure method is required. The turbulence model used to simulate the experimental results was the standard k-ε model [13]. A Reynolds stress turbulence model was also used in the simulation and compared for accuracy purposes only.

Particle Phase Model
The calculation of the particle motion comes directly from Newton's second law, which states: where F is the force on the particle, and m is the mass of the particle. Although this F is made up of a few different forces, the major contributing factor during dilute phase flow is the drag force between the continuum and the dispersed phases.

Drag Force Model
The drag force generally takes the form: where A CS is the cross-sectional area perpendicular to the velocity direction, ρ is the density of the continuum medium, V Rel is the relative velocity between the particle and the continuum phase and C D is the coefficient of drag of the particle. The standard model by [24] employed by CFX-5 is given by: were Re = ρ|V Rel |d µ (11) here d is the diameter of the particle, ρ is the fluid density and µ is the viscosity. This relation for the drag coefficient only relates the amount of drag to the relative velocity. For extremely dilute flows whereby the particles are well dispersed, this relation is adequate, but when the volume fraction of particles increases, the effects of neighbouring particles cannot be excluded. The following section covers the development of a new particle drag coefficient formula considering the effect of neighbouring particles on the drag force particles entrained in particle streams.

Turbulent Dispersion
The effect of the fluid turbulence on the particle motion is considered via the turbulent dispersion method. As the fluid flow is calculated as a time averaged steady-state flow, the velocity at any point can be quantified as a mean and fluctuating component of veloc- ities. The fluctuating components of velocity are calculated using Gaussian probability distribution of zero mean and standard deviation of: where k is the kinetic energy of the fluid in a particular cell. This calculated fluctuating component combined with the time-averaged mean value at that point are used to predict particle motion in that immediate area.

Particle-Wall Collision Model
The standard particle-wall collision based of the work of Shuen, Solomon, Zhang and Faeth [25] included in CFX-5 consists of particle being reflected off the wall surface based on the coefficient of restitution. This coefficient is a direct measure of the amount of energy lost during the collision process and effectively reduces the normal velocity of the rebounding particle. The figure shows the relationship between the pre and post collision velocities, both parallel and normal to wall. To model particle wall collision a rough-wall collision model similar to that of Sommerfeld [26] was implemented. This model considers the effects of slight undulations in the wall surface. These slight undulations result in slightly modified impact angles by introducing a virtual wall as seen in Figure 2. The calculation for post collision velocities is approached the same as normal wall, but the incoming velocities are adjusted to incorporate the slight impact angle change. The formulas to calculate the post-collision velocities for both sliding and non-sliding collisions are as follows: Non-sliding collision ( ) Sliding collision The calculation for post collision velocities is approached the same as normal wall, but the incoming velocities are adjusted to incorporate the slight impact angle change. The formulas to calculate the post-collision velocities for both sliding and non-sliding collisions are as follows: Non-sliding collision Sliding collision where u 1 and u 2 represents the corresponding velocities pre and post collision in the tangential direction with respect to the collision, v are the velocities normal to the wall, e is the coefficient of restitution and µ is the coefficient of friction. A non-sliding collision occurs when the following condition is valid: The equations for the standard model and improved are identical, but the incoming velocities for the improved model have been modified according the roughness angle based on the work of [26].

Particle-Particle Collision Model
The standard CFX-5 software [27] package does not account for particle-particle collision as part of the Lagrangian framework. As a consequence, the software needed to be modified to consider the effect of particle collisions through different flow systems. The modelling method implemented into the program is based on the model developed and refined by Sommerfeld [28]. The basis of this method is the creation of a fictitious particle that is used for the calculation of collision probability and if required the collision process.
As the Lagrangian particles are tracked individually through the flow domain, sampling the average particle velocities through a particular computational cell creates a fictitious particle. With the fictitious particle created, the collision probability of collisions occurring between the tracked and fictitious particles is calculated using a similar method for that of the kinetic theory of gases. The formula for the collision probability is as follows: where D is the diameter of the relevant particle, u the velocity of the relevant particle, n p is the number of particles per unit volume and ∆t is the time step. After the collision probability is calculated, this probability is compared with a random number (RN) that is generated. A collision is to be simulated if the generated RN is less than the P coll. When a collision is to be simulated, the location of collision point on the real particle is determined randomly and the resulting post collision velocities calculated using the follows formulas: where u 1 and u 1 are the pre and post collision velocities, respectively, in the tangential direction relative to the collision vector between particle center at the time of impact, e is the coefficient of restitution between the two particles, and m 1 and m 2 are the respective masses. In the normal direction, during the collision process, a particle may slide against the wall surface or depending on the friction between the faces, the particle may rotate in a non-sliding collision. The effects of the two sliding conditions are accounted for by: non-sliding collision sliding collision where e is the coefficient of restitution, µ is the coefficient of friction and m 1 and m 2 are the respective masses of the considered particles.

Description of the Simulated Cases
The experimental work of Liang, Hong and Fan [19] was used in the present work for numerical method. In their work, four different orientations were investigated. The first considered the effect of varying angle with respect to the flow direction at given separation distance. The second configuration looked a three co-aligned particle in the flow direction. The third test case considered a central test particle surrounded by 6 equally spaced particles in hexagonal configuration. The final case looked at a central test particle surround by 8 particles in cubic orientation. As previously mentioned the three co-aligned particles work forms the basis of the following study whereby separation distances from 0.25 to 7 x/d were investigated where the results for 0-3 x/d were compared with the experimental data of [19]. The published results considered the flow conditions where the particle Reynolds number equalled 54. The set up consists for the CFD model replicates those dimensions seen in the experimental work, with a pipe diameter of 15.24 cm and a particle diameter of 1.58 cm. The particles are co-aligned along the centreline of the pipe. The present CFD geometry based on the [19] experiments and inlet velocity profile is given in Figure 3. Figure 4 shows a section through the mid-plane of the pipe highlighting the mesh densities used. The mesh for these simulations was unstructured for the bulk of the domain, but regions close to the pipe wall and particle surfaces were structured in nature. These are the areas of greater gradients of flow variables. Beyond the structured region, heavy unstructured refinement was enforced around the particle surfaces and along the centreline of the pipe as this are the regions of most interest and the areas where the greater velocity gradients occur. To obtain an acceptable accuracy with acceptable computational time, various grid independency tests were conducted with different mesh resolutions. The predicted results were verified to have been no significant impact of the grid resolution on the results. A time step of 0.01 s was used in the present model. spaced particles in hexagonal configuration. The final case looked at a central test particle surround by 8 particles in cubic orientation. As previously mentioned the three co-aligned particles work forms the basis of the following study whereby separation distances from 0.25 to 7 x/d were investigated where the results for 0-3 x/d were compared with the experimental data of [19]. The published results considered the flow conditions where the particle Reynolds number equalled 54. The set up consists for the CFD model replicates those dimensions seen in the experimental work, with a pipe diameter of 15.24 cm and a particle diameter of 1.58 cm. The particles are co-aligned along the centreline of the pipe. The present CFD geometry based on the [19] experiments and inlet velocity profile is given in Figure 3. Figure 4 shows a section through the mid-plane of the pipe highlighting the mesh densities used. The mesh for these simulations was unstructured for the bulk of the domain, but regions close to the pipe wall and particle surfaces were structured in nature. These are the areas of greater gradients of flow variables. Beyond the structured region, heavy unstructured refinement was enforced around the particle surfaces and along the centreline of the pipe as this are the regions of most interest and the areas where the greater velocity gradients occur. To obtain an acceptable accuracy with acceptable computational time, various grid independency tests were conducted with different mesh resolutions. The predicted results were verified to have been no significant impact of the grid resolution on the results. A time step of 0.01 s was used in the present model.

Boundary Conditions
To accurately replicate the experimental conditions of Liang, Hong and Fan [19], a fully developed laminar pipe flow condition is used at the inlet boundary:  spaced particles in hexagonal configuration. The final case looked at a central test particle surround by 8 particles in cubic orientation. As previously mentioned the three co-aligned particles work forms the basis of the following study whereby separation distances from 0.25 to 7 x/d were investigated where the results for 0-3 x/d were compared with the ex perimental data of [19]. The published results considered the flow conditions where the particle Reynolds number equalled 54. The set up consists for the CFD model replicates those dimensions seen in the experimental work, with a pipe diameter of 15.24 cm and a particle diameter of 1.58 cm. The particles are co-aligned along the centreline of the pipe The present CFD geometry based on the [19] experiments and inlet velocity profile is given in Figure 3. Figure 4 shows a section through the mid-plane of the pipe highlighting the mesh densities used. The mesh for these simulations was unstructured for the bulk o the domain, but regions close to the pipe wall and particle surfaces were structured in nature. These are the areas of greater gradients of flow variables. Beyond the structured region, heavy unstructured refinement was enforced around the particle surfaces and along the centreline of the pipe as this are the regions of most interest and the areas where the greater velocity gradients occur. To obtain an acceptable accuracy with acceptable computational time, various grid independency tests were conducted with different mesh resolutions. The predicted results were verified to have been no significant impact of the grid resolution on the results. A time step of 0.01 s was used in the present model.

Boundary Conditions
To accurately replicate the experimental conditions of Liang, Hong and Fan [19], a fully developed laminar pipe flow condition is used at the inlet boundary:

Boundary Conditions
To accurately replicate the experimental conditions of Liang, Hong and Fan [19], a fully developed laminar pipe flow condition is used at the inlet boundary: Both the pipe wall and particle surfaces were set as no-slip boundary conditions with the velocity equal to zero. Outlet of the pipe was set to an opening type with a static pressure value equal to atmospheric pressure since it gives better convergence results. In present work the convergence criteria for all numerical simulations of all the residual were taken as 10 −4 for each scaled residual component. A glycerine/water solution of approximately 82 wt% glycerine was used in this study. Fluid properties were set according to the values given in Table 1.  Figure 5 shows the good agreement between the experimental study of Liang, Hong and Fan [19] co aligned particle work and the present CFD validation work, where cd is the drag force coefficient on the test particle and cd0 is the drag force coefficient on an isolated particle. The CFD model was able to predict the drag force accurately on all three particles.

  R
Both the pipe wall and particle surfaces were set as no-slip boundary conditio the velocity equal to zero. Outlet of the pipe was set to an opening type with a stat sure value equal to atmospheric pressure since it gives better convergence results sent work the convergence criteria for all numerical simulations of all the residu taken as 10 for each scaled residual component. A glycerine/water solution of a imately 82 wt% glycerine was used in this study. Fluid properties were set accor the values given in Table 1.  Figure 5 shows the good agreement between the experimental study of Liang and Fan [19] co aligned particle work and the present CFD validation work, whe the drag force coefficient on the test particle and cd0 is the drag force coefficien isolated particle. The CFD model was able to predict the drag force accurately on a particles. The results also accurately predicted the transition of drag force whereby the the middle particle surpasses that of the trailing particle. Figure 6 shows the veloc tour plot of the present CFD study at a separation distance of 5 particle diameter Reynolds number of 54. The leading particle is subjected to higher velocities than t ing particles which in turn leads to the higher drag force. The results also accurately predicted the transition of drag force whereby the drag on the middle particle surpasses that of the trailing particle. Figure 6 shows the velocity contour plot of the present CFD study at a separation distance of 5 particle diameters and  Figure 7 shows a close-up view of the recirculation area behind the three partic The leading particle, due to the higher incoming velocity, has a larger recirculation a and results in higher drag forces. Due to the close agreement of the predicting ability the CFD model, it is valid to assume an accurate extension of this work could be und taken.

Effect of separation distance on the coefficient of drag
To fully help understand the phenomenon at work in the particle drag of clusters particles, it is necessary to consider various cases. In the present work, three different ca were investigated: • Model 1 investigates the drag force of an average particle within an infinite line particles co aligned with the flow direction. This model will highlight the effects inline wake characteristics. • Model 2 considers a particle surrounded by an infinite plane of particles aligned p pendicular to the flow to investigate the influence of neighbouring particles. • Model 3 investigates an infinite matrix of particles, which is a combination and tension of the previous two models.
The results of cases 1 and two are for comparative purposes only and necessary fully understand the effect of position and location of particles on surrounding particl  Figure 7 shows a close-up view of the recirculation area behind the three particles. The leading particle, due to the higher incoming velocity, has a larger recirculation area and results in higher drag forces. Due to the close agreement of the predicting ability of the CFD model, it is valid to assume an accurate extension of this work could be undertaken.
To fully help understand the phenomenon at work in the particle drag of clusters of particles, it is necessary to consider various cases. In the present work, three different cases were investigated: • Model 1 investigates the drag force of an average particle within an infinite line of particles co aligned with the flow direction. This model will highlight the effects of inline wake characteristics. • Model 2 considers a particle surrounded by an infinite plane of particles aligned perpendicular to the flow to investigate the influence of neighbouring particles. • Model 3 investigates an infinite matrix of particles, which is a combination and extension of the previous two models.
The results of cases 1 and two are for comparative purposes only and necessary to fully understand the effect of position and location of particles on surrounding particles.

Particles Aligned with Flow (Model 1)
Model 1 looks at the case of infinite number particles co-aligned with the direction of flow (see, Figure 8a). This Model highlights the influence of particle wake on the drag force when aligned with the flow direction. The ranges of the simulations completed for this Model 1 include particle separation distances of 1 to 20 particle diameters and particle Reynolds numbers varying from 1 to 50.
In order to accurately model an infinite number of particles in flow direction, the inlet and outlet regions (see, Figure 8b), require a periodic boundary condition. The periodic boundary condition forces the quantities leaving the domain to be reintroduced into the domain through the inlet patch. The drag force created by the fluid moving passed the particles results in an energy loss due to friction between the shearing layers of fluid adjacent to the particle surface. The periodic boundary condition forms an endless loop of fluid that unless corrected will eventually stall. To overcome this phenomenon an equal amount of momentum is introduce evenly throughout the domain to ensure continuous flow.
• Model 2 considers a particle surrounded by an infinite plane of particles aligned perpendicular to the flow to investigate the influence of neighbouring particles. • Model 3 investigates an infinite matrix of particles, which is a combination and extension of the previous two models.
The results of cases 1 and two are for comparative purposes only and necessary to fully understand the effect of position and location of particles on surrounding particles. . Velocity vector plot of (a). leading particle, (b). middle particle, (c). trailing particle at x/d = 5d and Re = 54.

Particles Aligned with Flow (Model 1)
Model 1 looks at the case of infinite number particles co-aligned with the direction of flow (see, Figure 8a). This Model highlights the influence of particle wake on the drag force when aligned with the flow direction. The ranges of the simulations completed for this Model 1 include particle separation distances of 1 to 20 particle diameters and particle Reynolds numbers varying from 1 to 50.
In order to accurately model an infinite number of particles in flow direction, the inlet Figure 7. Velocity vector plot of (a). leading particle, (b). middle particle, (c). trailing particle at x/d = 5d and Re = 54. To test the introduction of firstly the momentum introduction method and the periodic boundary condition, a pseudo-isolated particle was simulated (separation distance of 80 pd) at various Reynolds numbers and compared with commonly accepted relation for drag force of isolated particles ( shows the close agreement between the theoretical value based on an isolated particle and the CFD simulations. In the cross-stream direction, the effects of the wall on the flow characteristics needed to be removed so a series of geometries were tested at varying distances to determine the minimum distance from the wall for the simulation. It was found the use of 20 d wall distance removed all wall effects from the flow. At a distance of 20 d, a symmetry boundary or free slip (shear stress = 0) was used. This kind of boundary condition removed the drag due to wall friction as is seen in pipe flows. To test the introduction of firstly the momentum introduction method and the periodic boundary condition, a pseudo-isolated particle was simulated (separation distance of 80 pd) at various Reynolds numbers and compared with commonly accepted relation for drag force of isolated particles (F D = 1 2 C D ρAV 2 where C D = 24 Re 1 + 0.15 Re 0.687 ). Figure 9 shows the close agreement between the theoretical value based on an isolated particle and the CFD simulations. In the cross-stream direction, the effects of the wall on the flow characteristics needed to be removed so a series of geometries were tested at varying distances to determine the minimum distance from the wall for the simulation. It was found the use of 20 d wall distance removed all wall effects from the flow. At a distance of 20 d, a symmetry boundary or free slip (shear stress = 0) was used. This kind of boundary condition removed the drag due to wall friction as is seen in pipe flows. To test the introduction of firstly the momentum introduction method and the periodic boundary condition, a pseudo-isolated particle was simulated (separation distance of 80 pd) at various Reynolds numbers and compared with commonly accepted relation for drag force of isolated particles ( shows the close agreement between the theoretical value based on an isolated particle and the CFD simulations. In the cross-stream direction, the effects of the wall on the flow characteristics needed to be removed so a series of geometries were tested at varying distances to determine the minimum distance from the wall for the simulation. It was found the use of 20 d wall distance removed all wall effects from the flow. At a distance of 20 d, a symmetry boundary or free slip (shear stress = 0) was used. This kind of boundary condition removed the drag due to wall friction as is seen in pipe flows.  Figure 10 shows the normalised drag force plot for Model 1, where f is the predicted drag force on the test particle and f0 is the corresponding drag force for an isolated particle. It can be noted that at very small separation distances, which represents high particle volume fractions, the reduction of drag is large for all Reynolds numbers. The larger the Reynolds numbers the greater the reduction in drag force. At a separation distance of 1 particle diameter and a Reynolds number of 50, the average particle drag force is approximately 18% of that of an isolated particle at the same Reynolds number. This is the phenomenon that is used by birds to fly great distance by sharing the drag force between the birds. The experimental work of Katz and Meneveau [29] experimentally studied the rise velocity of interactive bubbles and considered the relative between two identical bubble and several separation distances. The results suggest the relative velocity between the particles increased at smaller separation distances. This is consistent with the CFD predictions shown whereby a reduction in drag force would correspond to an increase velocity of the trailing bubble leading to greater relative velocity between to the two bubbles.  Figure 10 shows the normalised drag force plot for Model 1, where f is the predicted drag force on the test particle and f0 is the corresponding drag force for an isolated particle. It can be noted that at very small separation distances, which represents high particle volume fractions, the reduction of drag is large for all Reynolds numbers. The larger the Reynolds numbers the greater the reduction in drag force. At a separation distance of 1 particle diameter and a Reynolds number of 50, the average particle drag force is approximately 18% of that of an isolated particle at the same Reynolds number. This is the phenomenon that is used by birds to fly great distance by sharing the drag force between the birds. The experimental work of Katz and Meneveau [29] experimentally studied the rise velocity of interactive bubbles and considered the relative between two identical bubble and several separation distances. The results suggest the relative velocity between the particles increased at smaller separation distances. This is consistent with the CFD predictions shown whereby a reduction in drag force would correspond to an increase velocity of the trailing bubble leading to greater relative velocity between to the two bubbles. The velocity contours for some limited cases of Model 1 are shown in Figure 11. These plots highlight the reason the drag force at smaller separation distances. At 1 particle diameter separation distance, the wake from the previous particle engulfs the trailing particle leading to lower oncoming fluid flow velocity and consequentially lower drag force. The velocity contours for some limited cases of Model 1 are shown in Figure 11. These plots highlight the reason the drag force at smaller separation distances. At 1 particle diameter separation distance, the wake from the previous particle engulfs the trailing particle leading to lower oncoming fluid flow velocity and consequentially lower drag force.

Particles Perpendicular to Flow (Model 2)
Model 2 investigates the effects of particles located perpendicular to the flow (see, Figure 12a). The range of separation distance for this case includes 1 to 20 particle diameters. Model 2 focussed on the influence of neighbouring particle perpendicular to the flow direction and as such the cross-sectional dimensions were important. The length of the geometry was set to double the separation distance of the particle to the side boundaries. The geometrical design yields little effect on the drag force with the outlet set at this ratio to the side boundaries. This case required an infinite plane of particles perpendicular to the flow to be simulated. This plane is the product of using symmetry boundaries again on the wall regions parallel with the flow, but with this case, the distance to the symmetry boundary will varying accordingly to the separation distance being modelled. (See Figure 12b).

Particles Perpendicular to Flow (Model 2)
Model 2 investigates the effects of particles located perpendicular to the flow (see, Figure 12a). The range of separation distance for this case includes 1 to 20 particle diameters. Model 2 focussed on the influence of neighbouring particle perpendicular to the flow direction and as such the cross-sectional dimensions were important. The length of the geometry was set to double the separation distance of the particle to the side boundaries. The geometrical design yields little effect on the drag force with the outlet set at this ratio to the side boundaries. This case required an infinite plane of particles perpendicular to the flow to be simulated. This plane is the product of using symmetry boundaries again on the wall regions parallel with the flow, but with this case, the distance to the symmetry boundary will varying accordingly to the separation distance being modelled. (See Figure  12b).  Looking at the drag force results for Model 2 (see, Figure 13), the drag force on the test particle sharply increased compared to that of an isolated particle at the smallest separation distance whereas all the results showed a normalised drag force greater than unity with the greater separation distances showing very little difference from unity. The main Looking at the drag force results for Model 2 (see, Figure 13), the drag force on the test particle sharply increased compared to that of an isolated particle at the smallest separation distance whereas all the results showed a normalised drag force greater than unity with the greater separation distances showing very little difference from unity. The main reason the increase in drag force is due to the higher velocity experienced by the particle due to the fact that the cross section for the flow to move through decreases past the particle causing the velocity to increase to conserve continuity. This would explain why the results for the smallest geometry are so extreme because this is the case where the flow reduces the most. For example, the case of a separation distance of 1 particle diameter reduces the cross-sectional area past the particle by approximately 20%. Following along the same lines of thinking the smaller Reynolds number results provide greater force increases because the rate of change of drag coefficient and hence drag force with respect to Reynolds number is greater at lower numbers thereby leading to higher experienced drag force as the increase in seen velocity results in higher Reynolds numbers. The increased drag force attributed to the squeezing of the flow or the 'nozzle' effect was also noted by Chen and Lu [30]. In their experimental work, the drag force on 2 and 3 particles positioned side-by-side in the oncoming flow were measured. Their results suggest that at smaller separation distances, the drag force increases, which they primary conclude is based on the nozzle effect, which causes a rise in local Reynolds numbers. Although there is some difference in the magnitude of force increase between the predicted CFD results and those of the experimental work (see, Figure 14), it is noticeable from their work, that the case of three particle yielded slightly higher drag force increase compared with those of the two-particle case. With this in mind, the infinite array of particles perpendicular to flow in these CFD predictions would naturally expect to produce higher drag forces due to introduction of additional particles. The increased drag force attributed to the squeezing of the flow or the 'nozzle' effect was also noted by Chen and Lu [30]. In their experimental work, the drag force on 2 and 3 particles positioned side-by-side in the oncoming flow were measured. Their results suggest that at smaller separation distances, the drag force increases, which they primary conclude is based on the nozzle effect, which causes a rise in local Reynolds numbers. Although there is some difference in the magnitude of force increase between the predicted CFD results and those of the experimental work (see, Figure 14), it is noticeable from their work, that the case of three particle yielded slightly higher drag force increase compared with those of the two-particle case. With this in mind, the infinite array of particles perpendicular to flow in these CFD predictions would naturally expect to produce higher drag forces due to introduction of additional particles.
The velocity contour plots of Figure 15 show how the smaller separation distance greatly increases the maximum normalized velocity. The case of separation distance of 10 x/d shows almost a constant normalized velocity of close to unity resulting in a drag force ratio of approximately the same magnitude. This highlights the importance of the "nozzle" effect on the drag force for the smaller geometries caused by close neighbouring particles.
Although there is some difference in the magnitude of force increase between the predicted CFD results and those of the experimental work (see, Figure 14), it is noticeable from their work, that the case of three particle yielded slightly higher drag force increase compared with those of the two-particle case. With this in mind, the infinite array of particles perpendicular to flow in these CFD predictions would naturally expect to produce higher drag forces due to introduction of additional particles. The velocity contour plots of Figure 15 show how the smaller separation distance greatly increases the maximum normalized velocity. The case of separation distance of 10 x/d shows almost a constant normalized velocity of close to unity resulting in a drag force ratio of approximately the same magnitude. This highlights the importance of the "nozzle" effect on the drag force for the smaller geometries caused by close neighbouring particles.

Infinite Matrix of Particles (Model 3)
The final model consisted of a combination of Models 1 and 2 to form an infinite matrix of particles from which a new drag force was developed (see, Figure 16a). The geometry for Model 3 consisted of a cubic shape of dimensions based on separation distance. As would be expected being a combination of Models 1 and 2, the combination of boundary conditions were also applied to this model. As in Model 1, the periodic boundary condition and introduced momentum were used in Model 3. From Model 2, the symmetry boundary conditions were used (see, Figure 16).

Infinite Matrix of Particles (Model 3)
The final model consisted of a combination of Models 1 and 2 to form an infinite matrix of particles from which a new drag force was developed (see, Figure 16a). The geometry for Model 3 consisted of a cubic shape of dimensions based on separation distance. As would be expected being a combination of Models 1 and 2, the combination of boundary conditions were also applied to this model. As in Model 1, the periodic boundary condition and introduced momentum were used in Model 3. From Model 2, the symmetry boundary conditions were used (see, Figure 16). Results for the drag force of Model 3 are seen in Figure 17. For all the cases, there was a sharp increase in drag force at the smaller separation distances similar to that noted for Model 2. This increase again is caused due the increase of perceived velocity seen at the particle due to the constriction in the flow due to the physical presence of the particle. This would also explain why the effect is less noticeable at the greater separation distances. Also of interest is that as the Reynolds number increases the constriction effect is dramatically reduced mainly in part to fact that as rate of change of drag coefficient with respect to Reynolds decreases dramatically at higher Re.  Results for the drag force of Model 3 are seen in Figure 17. For all the cases, there was a sharp increase in drag force at the smaller separation distances similar to that noted for Model 2. This increase again is caused due the increase of perceived velocity seen at the particle due to the constriction in the flow due to the physical presence of the particle. This would also explain why the effect is less noticeable at the greater separation distances. Also of interest is that as the Reynolds number increases the constriction effect is dramatically reduced mainly in part to fact that as rate of change of drag coefficient with respect to Reynolds decreases dramatically at higher Re. Results for the drag force of Model 3 are seen in Figure 17. For all the cases, there w a sharp increase in drag force at the smaller separation distances similar to that noted f Model 2. This increase again is caused due the increase of perceived velocity seen at th particle due to the constriction in the flow due to the physical presence of the particle. Th would also explain why the effect is less noticeable at the greater separation distance Also of interest is that as the Reynolds number increases the constriction effect is drama ically reduced mainly in part to fact that as rate of change of drag coefficient with respe to Reynolds decreases dramatically at higher Re.   Figure 18 shows the velocity contour plots at various separation distances ranging from 10 to 1 x/d. At the smallest separation distance, it can be easily seen that the wake of the previous particle continues to the next, generally, this result in a reduction of drag force but due to the constriction effect previously discussed, the drag is higher. The constriction effect is highlighted in the contour's plots, as the normalised velocity of the smaller separation cases is higher than those of the 10 x/d case.
Fluids 2022, 7, x FOR PEER REVIEW Figure 18 shows the velocity contour plots at various separation distances from 10 to 1 x/d. At the smallest separation distance, it can be easily seen that the the previous particle continues to the next, generally, this result in a reduction force but due to the constriction effect previously discussed, the drag is higher. T striction effect is highlighted in the contour's plots, as the normalised velocity smaller separation cases is higher than those of the 10 x/d case. A series of comparison graphs of the three models are shown in Figure 19. T graph shows that at low Reynolds numbers, the results for Model 3 and Model most identical. This is due mainly to the constriction effect discussed early in the At a Reynolds number of 5 and above it is clear to see distinction between the Again, all the results show the sudden increase of drag at smaller separation d due to the change of flow area, something that is not seen in Model 1 because boundaries were set far enough away to minimise any area change to negligible tions. For the higher Reynolds number cases it can be seen that apart from the increase at smaller separation distance, as the particles spread, the results repli trends seen in Model 1 with the magnitudes of the drag force reduction slightly as it approaches unity. A series of comparison graphs of the three models are shown in Figure 19. The first graph shows that at low Reynolds numbers, the results for Model 3 and Model 2 are almost identical. This is due mainly to the constriction effect discussed early in the chapter. At a Reynolds number of 5 and above it is clear to see distinction between the models. Again, all the results show the sudden increase of drag at smaller separation distances due to the change of flow area, something that is not seen in Model 1 because the wall boundaries were set far enough away to minimise any area change to negligible proportions. For the higher Reynolds number cases it can be seen that apart from the sudden increase at smaller separation distance, as the particles spread, the results replicate the trends seen in Model 1 with the magnitudes of the drag force reduction slightly smaller as it approaches unity. Figure 20 shows the differences between Model 1 and 3 under the same flow conditions. It can be seen that the wake length produce in Model 1 is greater than that produced by Model 3. The difference in the models is the effect produced by the introduction of neighbouring particles upstream. It is interesting to note that effect of Model 2 beyond 5 particle diameters of separation distance was seen to be negligible, further downstream these particles impact on the wake length shortening and consequential drag force increase. The normalised drag coefficients shown in Figure 21 reinforce the influence of neighbouring particles from Model 2 as there is a clear differential between Model 1 and 3.  Figure 20 shows the differences between Model 1 and 3 under the same flow conditions. It can be seen that the wake length produce in Model 1 is greater than that produced by Model 3. The difference in the models is the effect produced by the introduction of neighbouring particles upstream. It is interesting to note that effect of Model 2 beyond 5 particle diameters of separation distance was seen to be negligible, further downstream these particles impact on the wake length shortening and consequential drag force increase. The normalised drag coefficients shown in Figure 21 reinforce the influence of neighbouring particles from Model 2 as there is a clear differential between Model 1 and 3.   The fitment of an equation to predict the drag coefficient of a spherical considering both Reynolds number and particle concentration (α) is based on the results seen in Figure  22. The equation development is based on Reynolds number between 1 and 50 and separation distance of 5 to 20 x/d. The separation distance corresponds particle volume fractions of 10 −2 and 5 × 10 −5 . To generate an equation for the data obtained from the present CFD prediction in the previous sub sections, a least squares method was utilised. A generic form of the equation was generated as seen in Equation (9). The least squares regression program solved for the constant A, B, C and D. The fitment of an equation to predict the drag coefficient of a spherical considering both Reynolds number and particle concentration (α) is based on the results seen in Figure 22. The equation development is based on Reynolds number between 1 and 50 and separation distance of 5 to 20 x/d. The separation distance corresponds particle volume fractions of 10 −2 and 5 × 10 −5 . To generate an equation for the data obtained from the present CFD prediction in the previous sub sections, a least squares method was utilised. A generic form of the equation was generated as seen in Equation (9). The least squares regression program solved for the constant A, B, C and D.  Substituting in the values from The results were more accurately predicted for those particles whose sizes were closest to the average diameter. Of the three different particle size ranges, the medium particle with an average of 226 µm produced the best CFD predictions with almost perfect predictions. The cases of the larger and smaller particle sizes (342 µm and 136 µm respectively)  Substituting in the values from Table 2, the generic equation takes the form of Equation (9).  The results were more accurately predicted for those particles whose sizes were closest to the average diameter. Of the three different particle size ranges, the medium particle with an average of 226 µm produced the best CFD predictions with almost perfect predictions. The cases of the larger and smaller particle sizes (342 µm and 136 µm respectively) while still showing improved predictions; they were not as prominent as the medium case. Several limitations to this newly developed drag model need to be acknowledged. Firstly, the particle concentration (α) range for this formula is considered to be adequate for the purpose of modelling dilute flows using the Lagrangian framework. In simulated flows where regions exceed the 0.01 volume fraction, limitations on the newly developed drag force, require the particle drag to be calculated using the standard drag model. This is not considered to be such a problem because at high particle concentrations (α) the particle motion is governed predominately by inter-particle collisions and less by drag force. Secondly, the limitation set on the particle Reynolds number is also considered to be adequate for the use in industrial type flows where the transport of small particles (less than 1 mm) is to be considered whereby generally the local particle Reynolds number would rarely exceed 50. In the case of larger particles producing the particle Reynolds numbers greater than 50 the new drag force will not be used, and the standard drag model would then be used. Finally, the present developed drag force formula has been calculated on a study of the relationship between the drag forces of mono-sized particles. Particles of differing diameter will obviously affect the resulting wake size and consequently the amount of drag reduction seen by trailing particles. A smaller particle (i.e., d p ≥ 136 µm) will not reduce the drag force as significantly on a larger particle (i.e., d p ≤ 342 µm) as would be expected if the situation were reversed. With this in mind, narrow particle size distributions are best suited for use with this formula, as the particles would be considered almost mono-sized. Any particles outside of this narrow band would be subject to incorrect reductions in drag force, particularly smaller particles, which may lead to variations from any validating experimental data. Although this limitation would exclude many applications if the new formula were to be used although the reductions in drag force may not be entirely accurate, it would still yield improvement over the standard drag force model.

Summary and Conclusions
The present work accurately reproduced the drag force results seen in the Liang, Hong and Fan [19] experimental study of three co-aligned particle at various separation distances. An extension of this work resulted in the CFD study of three models of varying particle configurations. Model 1 highlighted the important influence that particle wake has on the drag force of trailing particles. From this study, it was found that a trailing particle at small separation distances and a Reynolds number of 50 experience a reduction of up to 80% of drag force. Although not as prominent as this case all scenarios resulted in reduction of the seen drag force with this model. Model 2 considered the effect of neighbouring particle perpendicular to the flow. In all cases, the drag force increased due to a squeezing of the flow at the constriction, as expected this was much more prominent at the smaller separation distance. Most notably at a separation distance of 1 and a Reynolds number also 1 the drag force increased by over 250%. Beyond about 5 particle diameters this effect was almost negligible with drag approximating those of isolated particles at the same Reynolds numbers. For Model 3 for which the new drag equation was developed, at smaller separation distances its behaviour mimics that of Model 2 again because of the squeezing of flow. As the separation distance increased, the results followed more of the results seen for Model1. For both extremes, the increase or reduction of drag force were never as high as the results seen the models whose behaviour it was following. As the application of a new drag force model is to be used primarily for Lagrangian particle tracking where dilute flows are of most interest the range of separation distance was reduced to 5 to 20 particle diameters. This firstly made the fitment of a curve to the data easier and more accurate and secondly the importance of the drag force in dense phase flow becomes secondary to those of particle-particle collisions above and beyond a volume fraction of 0.01, which represents 5 particle diameters of separation distance.

Conflicts of Interest:
The authors have no conflict of interest to declare that are relevant to the content of this article.

Nomenclature
A surface area of the particle A CS cross-sectional area perpendicular to the velocity direction C D coefficient of drag of the particle D p particle diameter e Restitution coefficient F force on the particle F D drag force k kinetic energy of the fluid in a particular cell L distance between particle centres m p particle mass n p number of particles per unit volume P coll collision probability Re p Particle Reynolds number S φ, S Pφ source terms S ij strain tensor u velocity in the i th U p,g velocity of particle and gas respectively V Rel relative velocity x distance in the i th Greek letters ρ density of the continuum medium ε volume fraction φ variable quantity Γ φ diffusion quantity µ Viscosity