An iterative procedure for the solution of multi-fluid equations for gas-particle flows, with kinetic theory closures and frictional stress models for the description of the particle phase is presented and was implemented in the fully unstructured open source code OpenFOAM (R), capable of using arbitrarily-shaped polyhedral cells. Improved interpolation practices have been adopted to ensure the smoothness of the solution in situations where the phase effective density presents steep gradients, aiming at reproducing the discretisation typical of the staggered grid arrangement. This was achieved including the strongly varying source terms of the momentum equation in the momentum interpolation procedure. Phase coupling was implemented using the partial elimination algorithm, which was incorporated in the improved interpolation practice. An implicit treatment for the particle pressure was implemented to ensure the robustness of the algorithm whenever the particle phase approaches the packed condition, where strong and rapidly changing stresses as a function of the dispersed phase volume fraction are present. The implemented code was tested against a number of test cases, including a bubbling fluidised bed, and results were verified against the solution provided by another solver (MFIX), or against published numerical solutions. (C) 2011 Elsevier B.V. All rights reserved.