A computational algorithm for polydisperse bubbly flow is developed by combining quadrature-based moment methods (QBMM) with an existing two-fluid solver for gas–liquid flows. Care is taken to ensure that the two-fluid model equations are hyperbolic by generalizing the kinetic model for the bubble phase proposed by Bieseuvel and Gorissen (1990). The kinetic formulation for the bubble phase includes the full suite of interphase momentum exchange terms for bubbly flow, as well as ad hoc bubble–bubble interaction terms to model the transition from isolated bubbles to regions of pure air at very high bubble-phase volume fractions. A robust numerical algorithm to couple the QBMM approach with a gas–liquid two-fluid solver is proposed. The resulting algorithm is tested to show hyperbolicity, verified against the two-fluid model currently implemented into OpenFOAM, and validated against two sets of experiments on bubbly flows from the literature. In both cases, the computational method shows good agreement with experimental data, and improved accuracy in comparison to a two-fluid model considered for comparison purposes. The robustness of the algorithm is demonstrated on an unstructured mesh with a high superficial gas inlet velocity and source terms for coalescence and breakup. The resulting computational approach is implemented in the open-source CFD code OpenFOAM as part of the OpenQBMM project.