The present thesis focuses on the CFD analysis of natural circulation flows that are related to the passive safety loops of the forthcoming generation of water nuclear reactors. These flows involve inherent unsteadiness, laminar-to-turbulent transition and thermal stratification phenomena. As an initial validation study, a range of models that solve the Reynolds-Averaged Navier-Stokes equations (RANS), within the open-source CFD code Code_Saturne, are assessed in a forced convection 1-D channel flow case. A proposed numerical form of the Analytical Wall Function (AWF) for handling near-wall regions has been implemented within the code, and is further extended to account for buoyancy effects. A range of eddy-viscosity and second-moment closures with different near-wall treatments and different turbulent heat flux approximations are then tested using 2-D steady computations of the flow and heat transfer in a high aspect ratio differentially heated cavity. In this buoyancy driven flow, the k-epsilon model with AWF shows very promising results of the heat transfer and the turbulence mixing effects. Further comparisons of RANS models are performed in a differentially heated square cavity at Ra=1.58e9 and 10e11 with air using 2-D and 3-D transient simulations. It is shown that some models tend to laminarise the flow whilst others are deficient due to the log-law wall function. Significant improvements are shown when the effects of the 3-D structures are considered. The proposed strategy, AWF with k-epsilon model, outperforms the log-law based wall function, exhibiting robustness and improved Nusselt number and turbulence level predictions. The greatest resolution of this unsteady flow was achieved through a computationally demanding Large Eddy Simulation (LES), which was shown to agree closely with existing DNS data. The final series of 2-D and 3-D calculations concern unsteady RANS and conjugate heat transfer analysis for three vertical heater vertical cooler loop configurations. A range of flow regimes are studied, with modified Rayleigh numbers within 10e9-10e13, for water. The analysis mostly concerns the steady-state though some thermal transients and secondary motion in the cooler side are examined, and a way of appropriately non-dimensionalizing the predicted quantities is proposed. The developed strategy, k-epsilon with AWF, returned satisfactory predictions of the thermal and dynamic fields whereas the EBRSM and the Launder and Sharma form of the k-epsilon required a specific numerical treatment to maintain the flow turbulent. All computations were validated against either experimental or DNS data, or 1-D correlations where appropriate.