Liquid-solid two-phase flows are found in numerous operations in the chemical, petroleum, pharmaceutical and many other industries. In numerous cases, the mixture or slurry that flows is composed by a suspension of solid particles (dispersed phase) transported by a liquid (continuum phase). However, the large number and range of variables encountered in slurry flows, in the case of pipelines, cause the flow behavior of these slurry systems to vary over a wide range of conditions, and consequently, different approaches have been used to describe the behavior of different flow regimes. Therefore, there are numerous studies of particular cases that cover limited ranges of conditions. In consequence, the experimental approach is necessarily limited by geometric and physical scale factors. For these reasons, Computational Fluid Dynamics, CFD, constitutes an ideal technique for predicting the general flow behavior of these systems. CFD models in this area can be divided in two different classes: Eulerian-Eulerian and Lagrangian-Eulerian models. Differences between these models are related to the way the solid phase flow is represented. Lagrangian-Eulerian models calculate the path and motion of each particle, while Eulerian-Eulerian models treat the particle phase as a continuum and average out motion on the scale of individual particles. This work focuses on the Eulerian-Eulerian approach for modeling the flow of a mixture of sand particles and water in a horizontal pipe. Homogeneous and heterogeneous flow regimes are considered. The k-ε model was used for modeling turbulent effects. Additionally, closure of solid-phase momentum equations requires a description for the solid-phase stress. Constitutive relations for the solid-phase stress considering the inelastic nature of particle collisions based on the Gas Kinetic Theory concepts have been used. Governing equations are solved numerically using the control volume-based finite element method. An unstructured non-uniform grid was chosen to discretize the entire computational domain. A second-order scheme in space and time was used. Numerical solutions in fully developed turbulent flow were found. Results show that flow predictions are very sensitive to the restitution coefficient and pseudo-viscosity of the solid phase. The mean pressure gradients from numerical solutions were compared with results obtained using the correlations of Einstein, Thomas and Krieger for homogeneous cases and with experimental data found in the open literature for heterogeneous cases. The solutions were found to be in good agreement with both correlations and experimental data. In addition, these numerical results were closer to experimental data than results obtained using other numerical models.