This work focuses on the development of a comprehensive numerical methodology for the study of natural circulation systems using molten salts as working fluid, in particular for its application in the design of novel nuclear reactors designs. The methodology is conceived to be able to evaluate the feasibility of different natural circulation configurations through a detailed description of the thermal-hydraulic system and a stability analysis of its dynamic behavior. For the high amount of cases to be handled in a design stage, Reduced Order Models are used to reduce the corresponding computational cost. The numerical implementation consists in a CFD-based MATLAB-OpenFOAM-Python coupling and has been tested against literature results, showing good agreement. For a first assessment of the required capabilities a 2D case of a simple natural circulation configuration known as Rayleigh-Bénard convection has been used. The resulting dynamic states of the system have been studied for a range of both geometric and physical parameters, describing the different obtained solutions. The methodology is under further development for the treatment of more complex cases and is currently being used for the design of a natural circulation experiment.