The present study proposes a novel approach to efficiently simulate 4D Flow MRI acquisitions in realistic complex flow conditions. Navier-Stokes and Bloch equations are simultaneously solved with Eulerian-Lagrangian coupling. A semi-analytic solution for the Bloch equation as well as a periodic particle re-injection strategy are implemented to reduce the computational cost. The Bloch solver and the velocity reconstruction pipeline were first validated in a steady flow configuration. The coupled 4D Flow MRI simulation procedure was validated in a complex pulsatile flow phantom cardiovascular-typical experiment. Besides, we compared simulated MR velocity data with experimental 4D Flow MRI measurements.