Two-phase cooling solutions employing subcooled nucleate boiling flows e.g. thermosyphons, have gained a special interest during the last few decades. This interest stems from their enhanced ability to remove extremely high heat fluxes, while keeping a uniform surface temperature. Consequently, modelling and predicting boiling flows is very important, in order to optimise the two-phase cooling operation and to increase the involved heat transfer coefficients.  In this work, a subcooled boiling model is implemented in the open-source code OpenFOAM to improve and extend its existing solver reactingTwoPhaseEulerFoam dedicated to model boiling flows. These flows are modelled using Computational Fluid Dynamics (CFD) following the Eulerian two-fluid approach. The simulations are used to evaluate and analyse the existing Active Nucleation Site Density (ANSD) models in the literature. Based on this evaluation, the accuracy of the CFD simulations using existing boiling sub-models is determined, and features leading to improve this accuracy are highlighted. In addition, the CFD simulations are used to perform a sensitivity analysis of the interfacial forces acting on bubbles during boiling flows. Finally, CFD simulation data is employed to study the Onset of Nucleate Boiling (ONB) and to propose a new model for this boiling sub-model, with an improved prediction accuracy and extended validity range. It is shown in this work that predictions associated with existing boiling sub-models are not accurate, and such sub-models need to take into account several convective boiling quantities to improve their accuracy. These quantities are the thermophysical properties of the involved materials, liquid and vapour thermodynamic properties and the heated surface micro-structure properties. Regarding the interfacial momentum transfer, it is shown that all the interfacial forces have considerable effects on boiling, except the lift force, which can be neglected without influencing the simulations' output. The new proposed ONB model takes into account convective boiling features, and it able to predict the ONB with a very good accuracy with a standard deviation of 2.7% or 0.1 K. This new ONB model is valid for a wide range of inlet Reynolds numbers, covering both regimes, laminar and turbulent, and a wide range of inlet subcoolings and applied heat fluxes.