Phase contrast MRI is a powerful tool for the assessment of blood flow. However, especially in the highly complex and turbulent flow that accompanies many cardiovascular diseases, phase contrast MRI may suffer from artifacts. Simulation of phase contrast MRI of turbulent flow could increase our understanding of phase contrast MRI artifacts in turbulent flows and facilitate the development of phase contrast MRI methods for the assessment of turbulent blood flow. We present a method for the simulation of phase contrast MRI measurements of turbulent flow. The method uses an Eulerian-Lagrangian approach, in which spin particle trajectories are computed from time-resolved large eddy simulations. The Bloch equations are solved for each spin for a frame of reference moving along the spins trajectory. The method was validated by comparison with phase contrast MRI measurements of velocity and intravoxel velocity standard deviation (IVSD) on a flow phantom consisting of a straight rigid pipe with a stenosis. Turbulence related artifacts, such as signal drop and ghosting, could be recognized in the measurements as well as in the simulations. The velocity and the IVSD obtained from the magnitude of the phase contrast MRI simulations agreed well with the measurements.