The multi-size group (MUSIG) model is employed in this paper to simulate the gas–liquid two-phase flow in pump as turbine (PAT) since the traditional Eulerian–Eulerian two-fluid model is unable to take into account the phenomena of breakup and coalescence of bubbles. First, the simulation of gas–liquid two-phase flow in a square column is compared with the experiment to verify the accuracy of the MUSIG model. Then, the results of gas–liquid two-phase flow in PAT simulated by the MUSIG model are compared with those by the conventional uniform bubble (UB) model and find that the MUSIG model is more favorable to capture the flow pattern at high gas content compared to the UB model. Based on the MUSIG model, the internal flow characteristics, pressure fluctuation, and bubble size distribution of the PAT are analyzed. The rotation of the blades breaks a part of big bubbles into small bubbles in the volute, resulting in a smaller diameter of the bubbles entering the impeller. As the gas content increases, the number and size of vortices in the impeller flow channel increase. The vortex is formed at locations where the gas phase distribution in the impeller flow channel is concentrated. The outlet of the impeller is more prone to bubble consolidation under high gas content conditions. In conclusion, the MUSIG model can well predict the complex flow characteristics of gas–liquid two-phase inside the PAT and identify the key influencing factors of energy acquisition, which can provide support for improving the performance of the PAT design.