Full load instability is one of the most dangerous unsteady flow phenomena, and practically restricts the zone of stable operation of the whole hydraulic power plant. This paper presents a mathematical model and a numerical method for simulation of full load instability and the resulting pressure pulsations in Francis turbines. The model consists of one-dimensional hydro-acoustic equations for a long penstock domain, and three-dimensional Reynolds averaged Navier–Stokes equations of “liquid-vapour” flow for the turbine domain. Series of computations of a high head hydraulic power plant are carried out. Investigated are the sensitivities to time step and mesh size refinements as well as the effect of turbulence model. Thoma number and operating point dependencies of the computed amplitude and frequency of pressure pulsations are compared to measurements and to the predictions of the fully one-dimensional model of the power plant. The amplitude of the computed pressure and power oscillations agree well with the available experimental data, showing the potential of the presented approach to investigate and predict high load pulsations in hydraulic power plants.