Abstract: Bubbly flows are found in a large number of chemical engineering applications. For the computational fluid dynamics (CFD) simulations of such multi-phase flows, both physical models and numerical treatment are crucial to obtain robust and accurate results. In this numerical study, we investigate the two-fluid model (TFM) under challenging conditions such as phase segregation and inversion. For the phase segregation, a singular problem arises in the phase momentum and the two-phase k – ε equations when one phase fraction approaches zero. Another numerical issue is the accurate calculation of the drag coefficient, e.g., during the phase inversion. To address the singular problem, previous studies used a non-conservative formulation after dividing by the phase fraction; in our approach, we present a robust methodology for semi-conservative and fully conservative formulations. A special numerical treatment is introduced to the phase momentum equations and the turbulence equations, which avoids the singular problem in case of phase segregation. Concerning the drag force, two novel methods, the linear and the hyperbolic blending method, respectively, are presented to obtain accurate results. For testing the new numerical treatment, the analytical solution of a two-dimensional test case is first compared with the results predicted using a semi-conservative and a fully conservative formulation. The second test case investigated is a bubble column with different superficial velocities. The results from three-dimensional simulations using the novel formulations show good agreement with the literature data. Especially when phase segregation occurs, the semi-conservative and the fully conservative formulations using the two-phase k – ε model formulation converge.