Gas-liquid-solid three-phase flows (GLSTPF) have important applications in nature and industries, but they are complex and numerical simulation is challenging. In this study, a new solver coupling Eulerian-Eulerian-Lagrangian approach with population balance model (PBM) was developed in OpenFOAM to simulate GLSTPF. In the new solver, the effects of solid particles on gas-phase or liquid-phase holdups were included. Moreover, interphase forces of gas-solid and liquid-solid were calculated separately, and all the forces exerted on a solid particle were added to calculate the particle movement. Furthermore, the new solver was successfully compared with the experimental results of bubble size distribution and phase velocities in a three-phase bubble column. The introduction of PBM significantly improved the predictions of bubble rise and solid velocities (by up to 20%) and phase holdups (by up to 30%). The gas phase oscillated in the bubble column due to flow turbulence, and both the liquid and solid phases tended to circulate in the bubble column. The liquid vorticity exhibited irregular distributions, and different scales of vortex structure were observed with higher maximum radial vorticity (10 s−1) than axial vorticity (4 s−1), demonstrating the existence of swirling and oscillating flow in the column.