In this paper the three-dimensional model of gas-particle mixture turbulent flow in horizontal tubes and channels, as well as in complex pipes with bends was developed, in order to make a tool for analyzing processes in pneumatic conveying systems. Gas turbulence was modelled using k-ε model of turbulence. In getting the numerical solution for the gas phase a finite volume discretization scheme was used. Iterative procedure was based on SIMPLE algorithm. Numerical code CAST for single phase flow with colocated grid represented the basis for it. The presence of dispersed phase and its influence on gas phase was modelled by adding one additional source term in the equations of the gas phase (PSI-CELL method). Dispersed phase was treated by the Lagrangian approach. For particle motion LSD model was used [13]. The concept of parcel, the computational particle which represents the ensemble of real particles with the same performances was adopted. Particles were treated as ideal spheres. Beside the drag and the gravitational forces, the lift forces due to the particle rotation and the gradient of the gas velocity field were included in calculation. Special attention was devoted to the models of particle collision with the rough wall and mutual collisions of particles. On this way the use of the model on flows in which the particle-gas mass ratio is high was enabled. The stochastical approach was adopted, by which all the parameters with the stochastic nature in reality retain it in the model also. The modelling of roughness was performed by changing the wall surface with the 'virtual plane', whose position is detemined with the angle of inclination around horizontal axes (modelling of roughness height), and the angle of rotation around vertical axes (modelling of orientation of roughness in space). The first one was obeyed to the normal, and the second one to the uniform distribution. In the model of collisions between the particles the collision probability was calculated, and on this basis it was stochastically determined whether the collision of the particular particle with some other would happen or not. The calculation were made for particles of the mean diameter in the range 40-500 μm and for the particle-gas mass ratio 0-5. For the whole ensemble of particles in the flow: field the log-normal distribution of their diameters was adopted. The results of calculation were compared with the experimental results available in literature. For the profiles of gas velocity and pressure drop excellent agreement was achieved. The agreement of the concentration profiles of the dispersed phase was good, except for the regions of high particle concentration near the wall.