Three-phase flows involving a liquid, a gas, and a solid material are common in various industrial processes, including, but not limited to, slurry flow, particle encapsulation, and multiphase bubble column reactors. Numerical simulation in this context is a crucial tool for complementing experiments and theoretical analysis. Nevertheless, many phenomena happening at the microscopic scale are often overlooked when dealing with industrial-scale simulations, leading to approximate solutions that fail to capture correctly the underlying complex physics of these systems. The reason is the absence of correct closure models for three-phase flows, which need to be developed on a smaller scale before they can be applied to real-sized systems. In this work, we present our recently developed volume of fluid-immersed boundary method (VOF-IBM) solver for three-phase simulation. The IBM is a part of the CFDEMcoupling framework, which employs OpenFOAM for the fluid phase and LIGGGHTS for the solid phase. We propose an extension of VOF technique in OpenFOAM, which takes into account the particle presence when transporting the interface. We demonstrate the accuracy of the proposed tool by comparing numerical simulation results with our experiments [1] and with existing literature. The analysed setup includes heads-on collision of a single rising bubble with a settling rigid particle in different viscous liquids. The regime transition has been quantified and the simulations are extended toward a small-scale simulation setup with multiple bubbles and particles. Statistical analysis of the bubble-particle interactions paves the way to the development of accurate closure models for the Euler-Euler simulations of industrial three-phase flows.