Numerical simulations of cylindrical bubble column operating in the churn-turbulent regime have been simulated using Euler Euler approach incorporated with the RNG k-epsilon model for liquid turbulence. Single-size bubble model, double-size bubble model, and the multiple size group model (MUSIG), including the homogeneous and inhomogeneous discrete methods, are employed in the simulations. Mass conserved formulations of breakup and coalescence rates were used in the computation of bubble size distributions. The Schiller-Naumann drag force was used in the single-size model, and the Ishii-Zuber drag force was used for the MUSIG simulations. For the double-size bubble model, an empirical drag formulation was adapted. The simulation results of time-averaged axial velocity and gas holdup obtained with the three models were compared with reported experimental data in the literature. The results showed that only MUSIG models with lift force can reproduce the measured radial distribution of gas holdup in the fully developed flow regime and that the inhomogeneous MUSIC model performs a little better than other models in the prediction of axial liquid velocity. The RNG k-epsilon model was used in all simulations, and the results confirmed that this version of k-epsilon model did yield relatively high turbulence dissipation rates and high bubble breakup rates and, thus, resulted in a rational bubble size distribution. The ad hoc manipulation of the breakup rates was avoided. The simulation results indicated profound mutual effects of drag force, mean bubble sizes, and turbulence characteristics. An increase in drag force yielded a decrease in the relative velocity between phases, the later could result in decreases in k and epsilon. A large Sauter diameter results from a low bubble breakup rate which was directly connected to the dissipation rates of turbulence. The change of Suter diameter, in turn, influenced the drag force.