Unresolved mass motions are frequently detected in flares from extreme ultraviolet (EUV) observations, which are often regarded as turbulence. Non-thermal broadening of EUV emission lines caused by turbulence can be found at the entire flare region including flare loop top, legs, footpoints and the region above the looptop. Peaks in non-thermal velocity values tend to show up above the high density flare loops, reaching 100-200 km s^−1, while footpoints have lower non-thermal velocities of a few tens km s^−1. However, how this turbulence forms during the flare is still largely a mystery. Using a three-dimensional (3D) magnetohydrodynamic (MHD) simulation, we demonstrate how turbulent motions widely distribute throughout a flaring region, and can originate from a single source. The turbulence forms as a result of an intricate non-linear interaction between the reconnection outflows and the magnetic arcades below the reconnection site, in which the shear- flow driven Kelvin-Helmholtz Instability (KHI) plays a key role for generating turbulent vortices. The turbulence is produced above high density flare loops, and then propagates to chromospheric footpoints along the magnetic field as Alfvénic perturbations. The simulated strength and spatial distribution of the volume-filling turbulent motions show excellent agreement with observational results as revealed by synthetic views in EUV and by fitted Hinode-EIS spectra.