This work aimed to study the bubble distribution in a multiphase pump. A Euler-Euler inhomogeneous two-phase flow model coupled with a discrete particle population balance model (PBM) was used to simulate the whole flow channel of a three-stage gas-liquid two-phase centrifugal pump. Comparison of the computational fluid dynamic (CFD) simulation results with experimental data shows that the model can accurately predict the performance of the pump under various operating conditions. In addition, the liquid phase velocity distribution, gas-phase distribution, and pressure distribution of the second stage impeller at a 0.5 span of blade height under three typical working conditions were compared. Results show that the region with high local gas volume fraction (LGVF) mainly appears on the suction surface (SS) of the blade. With the increase in inlet gas volume fraction (IGVF), vortices and low velocity recirculation regions are generated at the impeller outlet and SS of the blade, the area with high LGVF increases, and gas–liquid separation occurs at the SS of the blade. The liquid phase flows out of the impeller at high velocity along the pressure surface of the blade, and the limited pressurization of fluid mainly happens at the impeller outlet. The average bubble size at the impeller outlet is the smallest while that at the impeller inlet is the largest. Under low IGVF conditions, bubbles tend to break into smaller ones, and the broken bubbles mainly concentrate at the blade pressure surface (PS) and the impeller outlet. Bubbles tend to coalesce into larger ones under high IGVF conditions. With the increase in IGVF, the bubble aggregation zone diffuses from the blade SS to the PS.