The one-dimensional two-fluid model approach has been traditionally used in thermal-hydraulics codes for the analysis of transients and accidents in water–cooled nuclear power plants. This paper investigates the performance of RELAP5/MOD3 predicting vertical upward bubbly flow at low velocity conditions. For bubbly flow and vertical pipes, this code applies the drift-velocity approach, showing important discrepancies with the experiments compared. Then, we use a classical formulation of the drag coefficient approach to evaluate the performance of both approaches. This is based on the critical Weber criteria and includes several assumptions for the calculation of the interfacial area and bubble size that are evaluated in this work. A more accurate drag coefficient approach is proposed and implemented in RELAP5/MOD3. Instead of using the Weber criteria, the bubble size distribution is directly considered. This allows the calculation of the interfacial area directly from the definition of Sauter mean diameter of a distribution. The results show that only the proposed approach was able to predict all the flow characteristics, in particular the bubble size and interfacial area concentration. Finally, the computational results are analyzed and validated with cross-section area average measurements of void fraction, dispersed phase velocity, bubble size, and interfacial area concentration.