An OpenFOAM computational fluid dynamics model setup is proposed for simulating thermally driven winds in mountain–valley systems. As a first step, the choice of Reynolds Averaged Navier–Stokes k−ε turbulence model is validated on a 3D geometry by comparing its results vs. large-eddy simulations reported in the literature. Then, a numerical model of an idealised 2D mountain–valley system with mountain slope angle of 20° is developed to simulate thermally driven winds. A couple of top surface boundary conditions (BC) and various combinations of temperature initial conditions (IC) are tested. A transient solver for buoyant, turbulent flow of incompressible fluids is used. Contrary to classical approaches where buoyancy is set as a variable of the problem, here temperature linearly dependent with altitude is imposed as BC on the slope and successfully leads to thermally driven wind generation. The minimum fluid domain height needed to properly simulate the thermally driven winds and the effects of the different setups on the results are discussed. Slip wall BC on the top surface of the fluid domain and uniform temperature IC are found to be the most adequate choices. Finally, valleys with different widths are simulated to see how the mountain–valley geometry affects the flow behaviour, both for anabatic (daytime, up-slope) and katabatic (nighttime, down-slope) winds. The simulations correctly reproduce the acceleration and deceleration of the flow along the slope. Increasing the valley width does not significantly affect the magnitude of the thermally driven wind but does produce a displacement of the generated convective cell.