Two new control volume solvers multiFluidInterFoam and rheoMultiFluidInterFoam are presented for the simulation of Newtonian and non-Newtonian n-phase flows (n≥2), respectively, fully accounting for interfacial tension and contact-angle effects for each phase. The multiFluidInterFoam solver modifies certain crucial aspects of the regular multiphaseInterFoam solver provided by OpenFOAM for Newtonian flows, improving its efficiency and robustness, but most importantly improving considerably its accuracy for surface tension driven flows. The new solver uses the volume-of-fluid (VOF) method based on the MULES method and artificial interface compression for the interface capturing, in combination with a suitable continuum surface force model; i.e. the interfacial tension coefficient decomposition method is employed to treat pairings of interfacial tension between the phases. VOF smoothing is also implemented by applying a Laplacian filter on the distribution of volume fractions, and the SIMPLEC algorithm is employed for the velocity–pressure coupling. The above algorithm has also been extended with the development of the rheoMultiFluidInterFoam solver which is capable of fully taking into account complex non-Newtonian effects (e.g. yield stress, viscoelasticity, etc.) of the involved liquids. To this end, the developed solver incorporates the RheoTool toolbox (Pimenta and Alves, 2017), utilizing a wealth of constitutive equations suitable for modeling different types of fluids with complex rheology. Our solvers are tested for both the implementation of the surface tension and of the constitutive equations in a wide range of different flow setups, considering typical benchmark two-phase and three-phase flows, such as the cases of dam break with an obstacle, spreading of a floating lens, a levitating drop, and a bubble rising in an ambient fluid; flows involving Newtonian, viscoplastic and viscoelastic liquids have been considered. Comparisons against analytical solutions or previously published numerical data are performed clearly demonstrating the capability of the new solvers to efficiently and accurately simulate interfacial tension dominated three-phase flows (and n-phase flows in general) for both simple and complex liquids.