Resolved CFD-DEM simulations of the hydraulic conveying of coarse grains through a very-narrow elbow

This paper investigates numerically the hydraulic conveying of solids through a 90$^{\circ}$ elbow that changes the flow direction from horizontal to vertical, in the very-narrow case where the ratio of pipe to particle diameters is less than 5. We performed resolved CFD-DEM (computational fluid dynamics - discrete element method) computations, in which we made use of the IB (immersed boundary) method of the open-source code CFDEM. We investigate the effects of the water flow and particle injection rate on the transport rate and sedimentation by tracking the granular structures appearing in the pipe, the motion of individual particles, and the contact network of settled particles. We found the saturated transport rate for each water velocity and that a large number of particles settle in the elbow region for smaller velocities, forming a crystal-like lattice that persists in time, and we propose a procedure to mitigate the problem.


Introduction
The hydraulic conveying of solid particles has been used for decades in industry as an effective way to extract and/or transport solids continuously within different facilities, being employed, for example, in mining, oil, chemical and food industries, and in wastewater treatment. Basically, it consists of pumping a liquid, usually water, that entrains solid particles as it flows through a tube or channel [1]. It is a viable and efficient alternative to trucks, rails and belts for transporting continuously high amounts of commodities, such as coal, bitumen, ores and other grains, and also for conveying organic matter to and from bioreactors and through sewer systems. Because these materials appear in a broad range of sizes, shapes and physical properties, it is not uncommon to find particles with sizes comparable to that of pipes, with even the occurrence of blockages and clogging [2,3]. In the specific case of organic matter, particles usually grow along time and develop an organic film with adhesion properties [4], increasing the probabilities of clustering, and thus blockages. In addition, piping systems frequently contain horizontal and vertical portions and changes of direction via elbows and tees, making the design of conveying systems complex. Therefore, the prediction of flow patterns, pressure drops, sedimentation and erosion rates, transport rates, clustering, clogging, and jamming in liquid-solid piping systems re-mains challenging [5,6].
In addition to pattern formation and instabilities, curved pipes can enhance particle-wall and particle-particle collisions, leading to pipe erosion and degradation of the product being transported [1,7,8]. Therefore, a great part of previous studies on pneumatic and hydraulic conveying through curved pipes investigated pipe erosion and pressure drop. For example, Bourgoyne Jr [9] investigated experimentally how sand particles being conveyed by gas and liquid in diverter systems erode different geometries of tees and bends, finding that erosion rates are two orders of magnitude higher when the fluid is a gas. Later, Shirazi et al. [10] proposed an estimation model to predict wear in tees and elbows for small particle concentrations (2-3% in weight) that agreed well with experiments. Using CFD-DEM (computational fluid dynamics -discrete element method), Zhang et al. [11] investigated the erosion of a slurry flow in a 90 • elbow by changing the flow velocity and elbow orientation, and found that the force (impulse) related with the impact of the particle upon the wall depends on the fluid velocity, while the puncture position depends on the elbow orientation. Similarly, Peng and Cao [12] studied numerically the erosion in a 90 • elbow, but for a gas-solid flow, and found that the profile of the particle concentration is directed linked to the erosion profile.
Other aspect influencing the hydraulic conveying of grains is the ratio between the pipe and particle diameters (D and d, respectively). When under high confinement, such as in narrow (10 D/d 100) and very-narrow tubes (D/d 10), different granular structures are observed [2,13,14].
Ravelet et al. [15] investigated experimentally the hydraulic conveying of coarse grains in horizontal pipes for 6 [20], which couples the open-source codes OpenFOAM (www.openfoam.com) and LIGGGHTS [21,22] for the CFD and DEM computations, respectively. We investigate the effects of the water flow and particle injection rate on the transport rates and sedimentation of grains. For that, we vary the water velocities within 0.06 and and 0.10 m/s and the particle insertion rate between one-quartersaturated and fully-saturated conditions, and identify and track the granular structures appearing in the pipe, the motion of individual particles, and the contact network of settled particles. In addition, we show how the mean water flow varies along the elbow. We find the saturated transport rate for each water velocity, and show that the mean particle velocity is an increasing function of the particle insertion rate, while particle settling in the elbow region decreases with the insertion rate. We also find that a large number of particles settles in the elbow region for smaller fluid velocities, forming a crystal-like lattice that persists in time, and we propose a simple procedure to mitigate the problem. These results can be used to avoid settling and grain accumulation in elbows and to improve the hydraulic conveying of coarse grains in industrial facilities.
The next sections present the governing equations, numerical setup, results, and the conclusions.

Governing equations
Our numerical simulations are of the Lagrangian-Eulerian type, where the motion of grains is computed in a Lagrangian framework and that of fluid in an Eulerian framework. We perform CFD-DEM computations, whose main equations are described next.
The solid particles are followed in a Lagrangian way by solving, for each particle, the linear and angular momentum equations, given by Eqs. 1 and 2, respectively, where, for each particle i, u i is the linear velocity, ω i the angular velocity, m i the mass, and I i the moment of inertia. F c,ij and F c,ik are, respectively, the contact force between particles i and j and the contact force between particle i and the wall k. T c,ij and T c,ik are the torque at the particle-particle and particle-wall contacts, respectively, F f,i is the force that the fluid exerts on the particle, and g is the acceleration of gravity. We do not consider momentum variations caused directly by the fluid in Eq. 2 because they are negligible with respect to contacts [23][24][25]. Contact forces between particles i and j are decomposed in normal and tangential components [26], and are given by Eqs. 3 and 4, respectively: where n and t are subscripts representing the normal and tangential directions, k and η are the spring and dashpot coefficients, respectively, δ represents the overlap between the particles, n ij is the vector that links particle centers, u ij is the relative velocity, and u slip,ij is the slip velocity at the contact point. t ij is the tangential vector, defined as u slip,ij /| u slip,ij |. The above relations also hold for particle-wall collisions. The contact forces are modeled following the Hertz-Mindlin and Deresiewicz model [27].
The liquid motion is computed in an Eulerian frame by using the mass and momentum equations given by Eqs. 5 and 6, respectively, with initial conditions given by Eq. 7, boundary conditions by Eq. 8, and the conditions at the solid-fluid interface by Eqs. 9 and 10, where Ω f is the fluid domain, Γ is the CFD boundary, Γ s represents the interface between the solid and the fluid, u f is the fluid velocity, u s is the velocity of the solid particle, u Γ is the boundary condition for the fluid velocity, u 0 is the initial condition for the fluid velocity, σ is the stress tensor, n is a unit vector normal to the solid surface, t is the traction vector of the fluid acting on the solid surface, and ν s is the kinematic viscosity of the fluid.
Equations 9 and 10 are responsible for the coupling between both phases: Eq. 9 matches the velocity between phases (no-slip condition) and Eq. 10 represents the stress that the fluid exerts on the solid particle. Therefore, the force acting on each particle is obtained by integrating the boundary condition given by Eq. 10 over the body's boundary Γ s , where the velocity used for the force calculation is weighted by the void fraction distribution on the interface. The same weighting average can be applied to the buoyancy force on each particle, when it is considered. Finally, the force contribution caused by the fluid on particle i is given by: where Ω s is the solid domain, V Ω S represents all cells covered or partially covered by the solid, V is the cell volume, and the index c stands for "cell".

Algorithm
The resolved four-way computations are implemented in CFDEM according with the following algorithm [26]: (i) DEM outputs the particle positions and velocities at a certain time step, and those values are passed to the CFD code.
(ii) A void fraction model identifies the regions covered by the particles and their surfaces, and dynamically refines the mesh on those regions.
(iii) The fluid flow is computed without considering the presence of particles.
(iv) Particle velocities are corrected in cells where they are present.
(v) The force that the fluid exerts on the particle is computed using the velocity and the pressure fields (Eq. 11) and passed to DEM.
(vi) The flow field is corrected to satisfy mass conservation.
(vii) The pressure is once more corrected and the routine restarts from step (i).

Numerical setup
We make use of the IB method of the open-source code CFDEM [20], which couples the open-source codes OpenFOAM (www.openfoam.com) and LIGGGHTS [21,22]. In the CFD part, the fluid flow is computed using the PISO (pressure-implicit with splitting of operators) algorithm with two main pressure corrections and two non-orthogonal flux corrections. Base meshes consist of hexahedral elements; however, because we use a dynamic grid refinement utility from OpenFOAM that splits some cells, creating tetrahedral elements, the non-orthogonal correctors are necessary. We use linear interpolation (central-difference scheme) to interpolate quantities from the center to the face of cells, and, for the transient term, we use the Euler discretization scheme, which is a first-order, bounded, implicit scheme. Figure 1 shows an example of mesh refinement used in the simulations.  Fluid enters the domain with a fixed value for the velocity and zerogradient condition for pressure, exits with a zero-gradient condition for the velocity and fixed value for pressure, and has a no-slip condition at the wall.
The time steps used for the CFD and DEM were, respectively, 1 × 10 −3 and 2 × 10 −5 s, leading to a coupling time of 50 DEM time steps. Those time steps were chosen in order to keep the DEM time step less than 10 % of the Rayleigh time [28] and ensure a coupling within 10 and 50 time steps between the DEM and CFD. Particle properties were chosen to represent Nylon 6-6 (our model for organic particles) and are listed in Tab. 1. The DEM time step was small enough to capture the collision span, which is related to the particle Young's modulus [23,24]. For this reason, the Young's modulus used in the simulations is two orders of magnitude lower than the real value, enhancing the DEM time step while still capturing the particle dynamics [29]. An example of simulation is available in Ref. [30], containing CFDEM input and output files for a case with particle rate equal to 50 particles/s and fluid velocity of 0.06 m/s (scripts for post-processing the outputs are also available).  Fig. 3a. Table 2 presents the parameters of the tests of Ref. [31], which are also used here. Because in this particular case we are not interested in the collision dynamics, but rather in that of settling, the DEM time step was set to 1 ×10 −5 s and that of CFD to 1 ×10 −4 s. The domain shown in Fig. 3a was discretized in three different ways in order to produce either 3, 4.5 or 6 cells per particle diameter, labeled, respectively, h 3 , h 2 and h 1 , and representing a total of 12800, 43200 and 102400 hexahedral elements each. In addition, we made use of a dynamic mesh refinement on the interface of each particle, as shown in Fig. 1. Table 3  show that convergence, indeed, is reached by increasing the mesh refinement, and they are quite satisfactory and close to the experimental ones [31] for an initial mesh refinement of 6 cells per diameter. Case For the tested cases, we measured the particle trajectories and the flow field around the spheres, and the agreement with the experiments of Ten   and Lawler [32], specially for 5 ≤ Re s ≤ 700.

Sedimentation velocity of a suspension of spheres
Finally, we validate our IB computations against the sedimentation velocity of suspensions of particles as measured by Richardson and Zaki [33]. U s is then computed as the final velocity of particles with respect to water.
The mesh and time steps used are basically the same as those described in  Figure 6 shows the sedimentation velocity U s normalized by the terminal velocity U t as a function of the void fraction (1 -ε p ), where ε p is the particle fraction. The symbols correspond to the numerical outputs and the continuous line to a logarithmic fitting. In the present case, the terminal velocity U t = 0.12 m/s was obtained by letting one single particle fall freely in still water, which corresponds to a terminal Reynolds number of Re t = U t d/ν = 720. For this value, the Richardson-Zaki correlation [33], has n = 2.4. From the fitting of our numerical results (Fig. 6), a coefficient n = 2.33 is obtained. By considering that the Richardson-Zaki correlation was obtained empirically and some dispersion has been shown to exist [35,36], our results agree reasonably well with that correlation.

Results
We investigate in this paper the effects of the water velocity and particle insertion rate in the transport rate and particle settling. Therefore, before starting the water velocity and insertion rate analyses, we performed two simulations to evaluate the terminal velocity in a 25.4-mm-ID vertical tube of both a single particle and a particle pack with particle fraction of 0.46, and found them to be U t = 0.12 m/s and U s = 0.075 m/s, respectively. We then investigate flow conditions in which water velocities are smaller, equal and higher than U s . A very high particle rate was imposed together with the three distinct fluid velocities, and we identified the maximum particle rate that each fluid velocity is able to drag. Once this saturation limit defined, we performed other two simulations for each fluid velocity, with particle rates corresponding to approximately half and a quarter of the respective saturation values. This procedure is illustrated in Fig. 7, which shows how the water velocities were chosen (above and below U s ) and the maximum particle rate for each water velocity determined. The parameters of the simulated cases are presented in Tab. 4, where U is the cross-sectional mean velocity, Γ is the particle rate, and Re = U D/ν f is the pipe Reynolds number.
Clearly, the maximum particle rate increases as the fluid velocity increases, The magnitudes of average particle velocities were evaluated at two different regions in the pipe, one in the horizontal and the other in the vertical section, both with a length of 1D. The horizontal region is located at a distance between 1.5D and 2.5D from the inlet, and the vertical region at a distance between 3D and 4D from the outlet (red regions in Fig. 2a). The averages were computed over all grains crossing those regions during the entire simulation time. We observe that the magnitudes of the particle average velocity in the horizontal section for cases I, IV and VII have values close  to those of the fluid at the inlet, and in the vertical section the magnitudes are slightly smaller (around 95% of those in the horizontal section). Interestingly, when the particle rate is decreased (by half in cases II, V and VIII, and by three-quarters in cases III, VI and IX) the magnitudes in the horizontal section decrease much faster than in the vertical section. Therefore, for Γ < Γ max , average velocities are smaller in the horizontal with respect to the vertical section (approximately 70-80% for Γ max /2 and 50-70% for Γ max /4). Larger velocities in the vertical than in the horizontal section when the insertion rate decreases indicates that grains are settling in the elbow region.
In order to understand that, we followed both the granular structures and individual grains along time.  and elbow, indicating that these regions are populated with grains forming an organized structure, and that they become less populated as U increases. In  Fig. 12c, no tracked particle accumulated in the elbow, given the higher particle mobility and lower settling times in this case. Interestingly, in all cases some particles did not reach the end of the vertical section. Instead, they fluctuated in that region rather than going straight to the outlet, in a behavior similar to those of fluidized beds [13,14]. We inquire into the liquid flow next. Figures 13a, 13b and 13c show the time-averaged velocities of the liquid in a lateral section of the pipe (plane y = 0) for cases II, V and VIII, respectively, in which the averages were besides the curvature effects in fluid velocities, we can notice that particles are organized in radial structures (Fig. 14a).
Finally, based on the observation that for smaller insertion rates particles settle and tend to accumulate on the bottom of the elbow, we investigate next a simple procedure to remove the settled structure once it has formed.
Since decreasing Γ worsens the problem, a solution that remains is to impose a transient by increasing U , even for saturated transport, in order to put settled grains back into motion. In principle, this solution cannot be taken for granted without further investigation because the formation of granular plugs has been reported for heavier particles [2,3]. However, the present case concerns much lighter particles that may follow the liquid more closely and not crystallize or jam. Therefore, we simulated Case II (0.06 m/s) for 10 s and, after a settled lattice has formed, increased U by steps of 0.01 m/s at 1 s intervals, keeping Γ constant. Figure 15  within the range of fluid, particles and geometry investigated, an increase of roughly 30% in the liquid velocity is capable of entraining gradually grains from the settled structure, an increase of 70% keeps that structure slipping over the wall, and an increase of 170% destroys the structure. Therefore, for piping systems conveying coarse grains of organic matter, this simple procedure can be used to break structures of settled particles in elbows.

Conclusions
This paper investigated numerically the hydraulic conveying of coarse grains through a very-narrow pipe (D/d = 4.23) curved by a 90 • elbow that changed the flow direction from horizontal to vertical. The CFD-DEM computations were made for light (organic) materials (S = 1.14), and made use of the IB method of the open-source code CFDEM. We varied the water velocity and particle insertion rate, and identified and tracked the granular structures appearing in the pipe, the motion of individual particles, and the contact network of settled particles. We found the maximum amount of particles transported by a given flow condition, that for either smaller velocities or particle insertion rates more particles settle in the elbow region and form a crystal-like lattice that persists in time, and that particle velocities in the horizontal section decrease much faster than in the vertical section as the particle rate is decreased. We observed also that part of the grains fluctuated in the vertical region rather than going straight to the outlet, in a behavior similar to those of fluidized beds. Finally, once the granular lattice is formed at the bottom of the elbow, we found that increasing the fluid velocity by 30% causes grains from the settled structure to be entrained by the flow, by 70% the lattice to slip over the wall, and by 170% the lattice to be destroyed. Because low insertion rates lead to particle settling, the increase of the liquid velocity can be used to break structures of settled particles in elbows, mitigating the accumulation of coarse grains of organic matter. Our results shed light on how to avoid settling and grain accumulation in elbows and improve the hydraulic conveying of solid particles in industrial facilities.

Declaration of Competing Interest
The authors declare no conflict of interest.