Multiple size group (MUSIG) model combined with a three-dimensional two-fluid model were employed to predict subcooled boiling flow of liquid nitrogen in a vertical upward tube. Based on the mechanism of boiling heat transfer, some important bubble model parameters were amended to be applicable to the modeling of liquid nitrogen. The distribution of different discrete bubble classes was demonstrated numerically and the distribution patterns of void fraction in the wall-heated tube were analyzed. It was found that the average void fraction increases nonlinearly along the axial direction with wall heat flux and it decreases with inlet mass flow rate and subcooled temperature. The local void fraction exhibited a U-shape distribution in the radial direction. The partition of the wall heat flux along the tube was obtained. The results showed that heat flux consumed on evaporation is the leading part of surface heat transfer at the rear region of subcooled boiling. The turning point in the pressure drop curve reflects the instability of bubbly flow. Good agreement was achieved on the local heat transfer coefficient against experimental measurements, which demonstrated the accuracy of the numerical model.