Numerical Simulation of Pneumatic Conveying in a Horizontal Pipe

A numerical simulation was attempted for pneumatic conveying of solids in a horizontal pipe. Trajectories of individual particles were calculated using equations of motion. In this simulation, the fluid drag, lift force due to particle rotation and torque on the rotating particles were taken into consideration. The friction loss due to collision of particles with a pipe wall was also calculated using impulsive equations. The pressure drop due to the presence of particles was obtained from the fluid drag acting on the particles through the momentum theorem. To avoid sliding motion of particles on a pipe bottom wall, a model of abnormal bouncing was newly proposed. Several parameters concerning the abnormal bouncing were determined empirically. It was found that the particle flow predicted by the present simulation agreed with measurements regarding particle distribution, pressure drop and particle velocities including angular velocities. In addition, this method was applied to find the effects of particle size, pipe diameter, particle density and so forth. Since particle diffusion due to the air turbulence was neglected in this analysis, the case of fine particles was not investigated.


Introduction
Technology of pneumatic conveying was developed in the middle of the 19th century for transportation of grains. Although a number of investigations have been carried out by many workers to improve facilities and refine design procedure, appropriate prediction of pressure drops is not still easy. This is because development of various industries has been increasing the number of solid materials to be conveyed and their properties of friction against a pipe wall depend largely upon each material. Therefore, it is very difficult to establish a unified theory quantitatively predicting the phenomena.
The methods of analysis for pneumatic conveying system may be classified into the follow-ing three stages.
( 1) One-dimensional analysis (highly empirical) (2) Numerical analysis based on individual particle trajectories, by neglecting the interaction between the particles and fluid. (3) Numerical analysis based on individual particle trajectories, by taking the interaction into consideration.

The first method
The analysis based upon the first method has been conducted by experiment depending upon the material and the condition. Although many monographs and papers describing this approach have been presented and various empirical factors have been introduced into the equations for pressure drop, the mechanism of friction between particle and wall is not directly considered.
The second method A trajectory of each particle is calculated numerically by computer simulation, where the values of mean particle velocity and pressure drop can be obtained as an ensemble average value.
It should be noticed that this method assumes that fluid motion is not affected by the presence of particles except for the loss of pressure energy; that is, the relationship between particle and fluid is a "one way path". This assumption is contradictory to the measurements by Laser Doppler anemometer even in a dilute-phase conveying 2 • 3 ). Fortunately it is known, however, that the assumption of the one way path causes no serious errors in pressure loss.

The third method
The analysis based upon the third method can also be carried out by computer simulations. As one of the simulations for single phase flows, the Navier-Stokes equation is usually used, by taking momentum exchange between particles and fluid into consideration. Although this method has been already used for several analyses of two-phase flows, more research will be necessary before it is applicable to pneumatic conveying. This is because there are several barriers to be overcome.
First, the situation of pneumatic conveying which differs from other types of two-phase flows is that collision of the particles against the wall should be dealt with. The collision reduces kinetic energy of the particles. The particles which lose the kinetic energy are accelerated by fluid through the drag force which results in the loss of pressure energy of the fluid. If the problem of the collision remains unsolved, the numerical solution of the Navier-Stokes equation pursued in detail will be meaningless. Therefore it is a pre-requisite to clear up the collision problem in the second method.
Secondly, the pipe flow used in pneumatic conveying is usually turbulent, and turbulence is also greatly influenced by the presence of particles. The problem how to deal with the interaction between turbulence and particle motion is beyond the scope of the present knowledge of pneumatic conveying, though it is interesting from the viewpoint of fluid mechanics. The difficulty due to turbulence could be avoided, if some convenient assumptions and models would be adopted in analyses. However, the third method based on such KONA No.3 ( 1 985) simple assumptions is not regarded superior to the first or second one.
The object of the present paper is to apply the second method to particle flow in pneumatic conveying. For convenience of computer simulation, the discussion is limited within the condition; particles used are assumed to be spherical because calculation becomes simplified, and large because they are little affected by fluid turbulence.

Outline of the present simulation
The behavior of particles transported in a pipeline is simple in a sence that particles are discharged from a feeder to a pipeline and travel under the influence of fluid-dynamic and gravitational forces with occasional collision against the wall. This is all that should be taken into consideration in the simulation. However, each particle shows a different trajectory, because the initial position and velocity have different values and the collision also causes a kind of randomness even in the case of spherical particles. This problem will be discussed in the later section. In the present simulation, the initial conditions of particle motion are given by random numbers. When equations of motion are given as well as equations of energy loss due to collision, the whole motion of each particle can be calculated. Practical data required for design works such as mean particle velocity, distribution of particles and pressure drop are obtainable by averaging the results of the calculation for many particles.
Fluid-dynamic forces acting on the particles are drag, lift and viscous dissipation of particle rotation. Existing knowledge of fluid mechanics gives sufficient information on these forces and thus one has no difficulty in estimating them. On the other hand, the effects of collision are difficult to estimate. If the particle size is much smaller than the pipe diameter, it can be assumed that the particle collides with a flat plate instead of the curved plane of a pipe wall. Change in particle velocity due to such collision can be estimated by using impulsive equations with the coefficients of restitution e and kinetic friction f. The coefficients of restitution and kinetic friction are dependent upon both the materials of particles and pipe wall. These values are obtained in an experiment where a particle falls on a plate if the material of the test plate and the pipe wall is the same. On the other hand, friction factors used in an one-dimensional analysis must be determined by the actual data of experiments.
The similar calculation was already applied to the simulation of pneumatic conveying including pipe bends 5 ). The results were satisfactory with respect to the effects of bends. However, this simulation neglected a factor which was necessary to enable the particles to continue bouncing motion. Therefore when a particle repeats the bouncing motion many times in a horizontal and long pipe with the coefficient of restitution e less than unity, the incidence angle of the particle against the wall decreases till the particle begins to slide on the wall. In the present work, artificial irregularity is assumed in the bouncing motion to prevent such sliding so that the simulation is applicable to a long pipe.

Collision between a single particle and wall
Consider a three dimensional collision of a spherical particle with a flat plate as shown in Fig. I. When the coefficients of restitution and kinetic friction are given, particle motions before and after the collision can be estimated by solving impulsive equations. This may be based upon classical dynamics. To solve the collision problem, the following assumptions were used. (1) Plastic deformation and fracture are neglected. (2) There exists a period when the particle slides on the plate. (3) Once the particle stops sliding motion, it slides no longer. ( 4) The friction between the particle and plate obeys Coulomb's law. The distance between the center of particle mass and contact point is kept constant a throughout the process of collision. The process of collision is divided into two periods: the period during which the material in collision is compressed and the period during which this compression is released. The former and the latter are called the compression and recovery periods, respectively. The coefficient of restitution e is defined by (1) where Jn(t) is the normal component of the impulse of the force which acts on a particle in the compression period and Jn ( 2 ) is the one in the recovery period. Coulomb's law says that the impulse of the friction force is the product of the impulse of the normal force and the coefficient of kinetic friction f. Thus when the particle slides on the X-Z plane, the impulse Jt of the friction force is expressed as (2) where i and k are unit vectors corresponding to X and Z directions, satisfying the following relation (3) Signs of Ex and Ez are defined to be equal to those of v,(o)+aw<o) and v,<o)_aw<o) respectively X X X X' · In general, the particle comes into collision incorporated with particle rotation. Thus u<o), the velocity of the point, on which the particle is in contact with the wall is the sum of the velocities of translation and rotation, that is, The impulsive equations have different forms according to the following three cases.
( 1) Case I: the particle stops sliding in the compression period (2) Case II: the particle stops sliding in therecovery period. (3) Case ill: the particle continues to slide throughout the collision. The procedure to obtain the solution of Case I alone is shown in the following.   Table I shows the velocities, angular velocities and impulses for each moment and period which must be considered in Case I. The known variables are only v<o)and w<o) All the velocities have three components corresponding to X, Y and Z directions, and thus there are 18 unknowns regarding these velocities.
Furthermore, the impulses J(s), J(r) and f 2 ) are expressed as The above equations means that there are 6 unknowns regarding the impulses. Therefore, unknown variables of Case I total 26 at the beginning, since ex and e 2 are added as unknowns.
Considering the momentum exchange at each moment, we have the following equations.
The boundary conditions are (16) It is found that there are 26 equations (Eq. (3) and the three components of Eqs. (8) to (16)) are given for 26 unknowns, and then the problem can be solved if v<o) and w<o) are known.
The solutions of Case I are shown in the column ( 1) of Table 2. The condition that the particle stops sliding in the compression period is also shown in terms of the inequality in the same column which is deduced from the following condition, The solutions of Case ll and Case Ill, which are the same to each other, are shown in the column (2) of Table 2. 4. Determination of the coefficients of restitution and kinetic friction Table 2 also provides a method of estimating the coefficients of restitution and kinetic friction if the particle velocities before and after collision are known. Consider a particle falling from a certain height to a plate inclined at the angle of a to the vertical line as shown in Fig. 2. That is, a two-dimensional collision with zero w<o) is dealt with. Velocity ratios Vx IV l) and Vy I V ) 0 ) can be measured by using photographic method such as a camera or others. The coefficient of restitution is given by (18) The ratio VxiV} 0 ) depends on a, e and[, i.e.
The velocity ratios calculated by Eqs. (18)

1 Coordinate system and initial condition
For practical pneumatic lines, particles are usually supplied downward from a hopper to a pipeline. Thus, the coordinate system shown in Fig. 4 is chosen in the present paper. In this system, x is the axial distance along the pipe, y, the horizontal axis normal to x, and z, the vertical axis in which the upward direction is taken positive. Initial conditions of each particle were specified as follows. Although both actual falling velocities and positions are considered to be random, it is sufficient to take only the falling position as a random variable for the simulation. The initial position of a particle is given by Yo and 8 as shown in Fig 4 (b).
Values of Yo and 8 are given by the random numbers distributed uniformly within allowable ranges. The initial velocity of the particle i = (dz/dt) depends on the method of feeding.
Fortunately, some empirical formulas are avail-able for such falling velocities when the particles are discharged from the hopper 8 ).

2 Fluid-dynamic forces and equations of motion
The equation of particle motion is written by (25) where X is the radius vector expressed as xi + yj + zk. The values of Fv, FL and FG are forces due to drag, lift and gravity, respectively. Viscous dissipation due to rotation of the particle must be also taken into account. The assumptions used in the calculation of particle trajectories are as follows. ( 1) Existing results obtained in a uniform and stationary flow field are applied to the fluid-dynamic forces acting on the particle. (2) The velocity distribution of the air obeys the 1/7 power law, and it is not affected by the presence of particles. Also, particle motion is not influenced by air turbulence. (3) Collision between particles is neglected.
The magnitude of the drag force is expressed as (26) where CD is the drag coefficient, A, cross sectional area of the particle, P, fluid density and U is the relative velocity defined by by many workers. Next, let us consider the lift force due to the particle rotation. Collision of a particle against the wall causes the particle to rotate at a considerably high rate of rotation, for instance, 10 3 ~ 1 0 4 rad/sec in the case of a coarse particle, which was confirmed in an experiment by Matsumoto,et al. 9 >. Therefore calculated trajectories are different whether the particle rotation is taken into account or not. When a particle rotates with translational velocity U as shown in Fig. 5, the magnitude of the lift force is expressed as where eL is the lift coefficient given as a function of the velocity ratio a, a=aw/U The following equations are relation between eL and a, assumed for the (30) by referring to the results of Maccol 10 >. The vector U of the relative velocity as well as w has three components of x, y and z directions, but taking all the components is so complicated that the following approximation is used. For the relative velocity, only one component is considered. The vector of angular velocity w is given as (32) This suggests that the lift force is assumed to be applied only in the (y, z) plane. As a result of the above approximation, the magnitude of the lift and velocity ratio a can The gravity force Fe is given by Finally, one obtains the equations of motion in the following form.

3 Pressure drop
If the purpose of the present simulation were to obtain only the particle trajectories, it would not be of value from the viewpoint of industry, because what design workers really want to know is the pressure drop rather than the particle trajectories. The pressure drop can The total pressure drop which occurs in the region dx is assumed to be the sum of the pressure drop due to air, dpa, and the additional one, dPs· We can use some empirical equations such as Blasius' or Plandtl's formula to estimate the pressure drop dpa. The pressure drop dps is caused by the drag force exerted by fluid on particles. Hence, the equation of momentum balance in the region dx becomes (46) where S is the cross sectional area of the pipe and N is the number of particles contained in the region dx. Definition of the mass flow ratio of particle to air gives Trajectories and pressure distributions are obtained by integrating Eqs. (43) and (48) simultaneously. In order to make the averaged dps KONA No. 3 (1985) converge to a certain value, calculations must be made for many particles with different initial conditions.

4 Change in particle rotation due to viscous dissipation
First, the authors explain the reason to consider the viscous dissipation for particle rotation. As was mentioned before, particles conveyed in the pipe rotate at high angular velocities. Obviously, the increase in the angular velocity is limited to a certain level in actual pneumatic conveying; however when a factor which suppresses the rotation is not considered in the simulation, the angular velocity continuously increases along the flow owing to repeated collision against the wall. This tendency is marked particularly when the particlewall collision takes place in the lower half of the pipe section. As shown in Section 3, the longitudinal velocity after the collision is influenced by the rotation before the collision. In the simulation where a factor of the dissipation is ignored, the longitudinal velocity increases with increasing rotation so that the relative velocity between the air and particle decreases. As a result, calculated pressure drops become smaller than measurements at a point far from the feeder.

(49)
where Tis the torque exserted on the rotating particle and R is the Reynolds number defined by a 2 w /v. When the dissipation is taken into consideration, the angular velocity settles asymptotically to a certain value as shown later.
5. 5 Abnormal bouncing of particle against the wall What are necessary for the simulation are almost explained in the foregoing section except for one thing that is "abnormal bouncing". If a spherical particle repeats bouncing motion according to the relationship shown in Table 2, the height of the particle rebounded from the wall gradually decreases and it can not rebound  Fig. 7 Bouncing motion in the pipe at last even at a high velocity of air, as shown in Fig. 7 (b). For a practical pneumatic conveying, however, one finds that particles repeat bouncing motion steadily as shown in Fig. 7(a). At least for large particles, the Magnus force, air turbulence and lift force due to velocity shear can not be the causes which prevent the particles from settling on the bottom wall. The reason why the particles continue the bouncing motion in a long horizontal pipe is that the normal component of the velocity of the particle rebounded from the wall does not diminish. Such bouncing motion is defined as abnormal bouncing in this paper.
There have been several models describing the abnormal bouncing proposed, for example, the pipe wall has a slight roughness 12 ) or particle shape deviates from that of a sphere. In the present work, the authors propose newly a virtual wall model. This model states that when the incidence angle becomes smaller than a certain value, the wall is replaced by a virtual one with the angle a against the true wall, as shown in Fig. 8. The velocity after collision with this virtual wall must be determined by the impulsive equations, the results of which are shown in Table 2. In this model, the energy conservation principle holds. Although the particle gains more energy in the normal direction, it loses much more energy in the longitudinal direction due to the collision against the inclined wall.
The relation between the angle a and incidence angle 8 is assumed as This model was qualitatively justified by the experiment of the particle bouncing on a rotat-46 -------- Fig. 8 Virtual wall model of abnormal bouncing ing disc 13 ). It is further assumed that the virtual wall makes an angle against the tangential plane of the pipe wall. The angle 'Y, which is called "yaw angle" in this paper, varies in the range and the value of 'Y is given by random numbers. The reason why the yaw angle is needed is that without a randomness due to the yaw angle particles tend to concentrate on the same vertical line in a long pipe and so particle distribution becomes unrealistic.
The next step is to determine the parameters {3, 8 and 'Y 0 • These parameters were determined in a purely empirical way at the present stage of analysis. What the authors aim at in the present simulation is to explain satisfactorily all the quantities such as translational particle velocity, angular velocity, particle distribution and pressure drop. Good agreements between simulation and experiment can be obtained by choosing proper values of {3, 8 and 'Y 0 • However, it is almost impossible to adjust the parameters to every possible condition of conveying. Hence the values of parameters were determined empirically by comparing simulated results with test data. The same parameter were applied to other cases under different conditions. The experiment to determine the parameters was made under the following conditions: Pipe where Fr is defined as Fr = fl / v'i/5. } (52) 6. Results

1 Particle trajectory
Particle trajectories calculated for polyethylene pellets with 1 .1 mm and 50 )..!ill (Ps = 923 kg/m 3 ) are presented in Fig. 9. It is seen in Fig. 9 that a particle constantly repeats bouncing on the wall and that a free path increases with increasing air velocity. On the other hand, the particle with 50 )..!ill does not show the bouncing motion or if any, it bounces with very short steps. At a high air velocity it does not touch the wall within the distance shown in Fig. 9. This means that abnormal bouncing is not enough to make a small particle suspended in the horizontal flow. Diffusion due to air turbulence must be taken into account in order to avoid such contradiction. The minimum size of particle with the density of the order 10 3 kg/m 3 to which the present simulation is applicable is about 0.5 mm.
6. 2 Particle distribution and angular velocity Figure 10 presents particle distributions at various longitudinal sections. The Magnus lift force is taken into account in Fig. 10 (b), while it is not in Fig. 10(a). The results of Fig. 10 (b) are more realistic than those of Fig. 10 (a). However, one should notice that the Magnus force is not only the factor to give a good agreement with experimental distributions, because the calculated distribution is easily changed by the parameter ~, o and 'Y 0 . What is meant in Fig. 10 is that the Magnus force clearly affects the particle distribution even for large particle although it is not enough to make the particles being suspended.
To show the state of distribution quantitatively, the pipe section is divided into 24 subsections as shown in Fig. 11, and the height of the distribution center is defined by where fi is the number of particles contained in the ith divided section and z~; is the height of gravity center of the ith sub-section measured from the bottom line y'. Variations of Zg along the pipeline are shown in Fig. 12 in which plotted points represent the results of KONA No.3 (1985) the simulation and three curves represent measurements. It is found that agreements between the simulations and experiments are satisfactory, although there are some differences observed if one compares the results in detail. Figure 13 indicates that the angular velocity settles at an asymptotic value owing to viscous dissipation of particle rotation.

3 Pressure drop
As mentioned before, the authors place a priority of this work on predicting the pressure drop, because the pressure drop has been always investigated as a central subject for pneumatic conveying. Many workers have been seeking a similarity law of the pressure drop like that of the single phase flow. Generally, the situation in the two-phase flows is so complicated that it is hard to establish the similarity law to the same preciseness as the single phase flow. However it has been found 15 ) that the following coefficient As plotted against the Froude number Fr * = c/..Jil) provides a comparatively good similarity in the pressure drop data of coarse particles.
where Pds is the dispersed density of the particles. Hence the results of pressure drops are shown in the form of As vs F/ in this paper.
The values of As are indicated in Fig. 14 where plotted points represent the simulated result and the dashed curves show the range of corresponding measurements. It is seen that agreement between them is fairly good.
All the foregoing results are limited to the condition that the particle diameter is 1.1 mm and the pipe diameter is 52 or 42 mm. The results predicted by the simulation are presented in the following section, where the effects of various factors such as particle and pipe diameters and particle density are shown. The following results are mainly described about the coefficient of the additional pressure drop As.
6. 4 Effects of particle and pipe diameters Figure 15 shows the relation between As and F/ calculated for various particle sizes. The dashed curves given for comparison are the same ones that are shown in Fig. 14. The results of the simulation are within the range between the curves. As the particle size becomes larger, the calculated results tend to deviate from the range of the curves. The effects of the pipe diameter are shown in Fig. 16, under the condition that pipe diameters are considerably different. Although the difference of the diameters of the two pipes with 27 and 200 mm is about one figure, the simulation indicates that there is little influence of the pipe diameter on the relation of As and F/. As mentioned before, the pressure drop data plotted in the form of As and F/ showed a comparatively good similarity. The present results support such an experimental fact to some extent.

5 Effects of particle density and kinetic friction
Results of the cases of Ps = 500 kg/m 3 and 2500 kg/m 3 are compared to see the effects of particle density. Figure 17 shows the variation of particle to air velocity ratios. The particle and pipe diameters used are 1.1 and 52 mm,  Figure 19 shows the effects of friction coefficient f on the relation of As and F/. It is found that the higher value off (for instance f = 0.5) gives slightly higher values of As for the same F/.

Conclusion
Trajectories and pressure drop coefficient in pneumatic conveying for large spherical particles could be predicted successfully by the present simulation.