The constant bubble size modeling approach (CBSM) and variable bubble size modeling approach (VBSM) are frequently employed in Eulerian–Eulerian simulation of bubble columns. However, the accuracy of CBSM is limited while the computational efficiency of VBSM needs to be improved. This work aims to develop method for bubble size modeling which has high computational efficiency and accuracy in the simulation of bubble columns. The distribution of bubble sizes is represented by a series of discrete points, and the percentage of bubbles with various sizes at gas inlet is determined by the results of computational fluid dynamics (CFD)–population balance model (PBM) simulations, whereas the influence of bubble coalescence and breakup is neglected. The simulated results of a 0.15 m diameter bubble column suggest that the developed method has high computational speed and can achieve similar accuracy as CFD–PBM modeling. Furthermore, the convergence issues caused by solving population balance equations are addressed.