I am trying to simulate decaying homogeneous isotropic turbulence. As initial condition I want a divergence-free vector-field, i.e, $\mathrm{div} U = 0$.
How do I initialize random velocity field in an uniform grid such that above condition is satisfied?
A purely random initial condition will not resemble realistic turbulence and will take quite some time to reach a realistic decay state. If you want good results it's more complicated than just ensuring zero divergence. Usually the velocity field is generated in spectral space to give the proper spectral energy distribution and is then transformed to physical space to initialize the simulation.
Section 4.1 of the following paper outlines one way to do this:
Huidan Yu, Sharath S. Girimaji, Li-Shi Luo, DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method, Journal of Computational Physics, Volume 209, Issue 2, 1 November 2005, Pages 599-616, ISSN 0021-9991, http://dx.doi.org/10.1016/j.jcp.2005.03.022.
This section begins:
We conduct the simulations in a 3D periodic cube with various resolutions $N^3$. The initial incompressible homogeneous isotropic velocity field $\vec{u}_0$ (${\nabla} \cdot \vec{u}_0 = 0$) is generated in spectral space $\kappa$ with the following energy spectrum in a prescribed range $\kappa_{\min} \le \kappa \le \kappa_{\max}$ and a random phase (cf. details in [47]):
\begin{align} \hat{E}(\kappa,0) = \left\{ \begin{array} 0.038\kappa^m \exp(-0.14\kappa^2) & \kappa\in[\kappa_{\min},\kappa_{\max}],\\ 0 & \kappa\notin[\kappa_{\min},\kappa_{\max}]. \end{array} \right.~~~~~~~~~~~~(16) \end{align}Then the velocity field is transferred to physical space. In what follows, we use $m = 4$ or $2$ in Eq. (16) to investigate the effect of $m$ on the energy spectrum and other quantities.
...
[47] T. Miyauchi, T. Ishizu, Direct numerical simulation of homogeneous isotropic turbulence – decay of passive scalar fluctuation, in Preprint vol. JSME No. 914-2, 1991, pp. 166–168.
You might also want to look at:
Consistent initial conditions for the DNS of compressible turbulence
Ristorcelli, J. R. and Blaisdell, G. A., Physics of Fluids, 9, 4-6 (1997), DOI:http://dx.doi.org/10.1063/1.869152
which could also be used for incompressible turbulence by using only the first step in their procedure which involves generating an appropriate divergence-free velocity field.
That paper chooses a solenoidal (zero divergence) velocity field with spectrum $E(\kappa) = A\kappa^4\exp(-2\kappa^2/\kappa_p^2)$ and random phase components. $A$ determines the velocity amplitude. $\kappa_p$ determines the peak in the spectral energy and is set equal to $12$ in that paper.
Note that these two methods are equivalent with appropriate choices of parameters.
Edit:
As for the specifics of generating a zero-divergence velocity with a desired energy spectrum, that question has been asked here but has not been answered. Perhaps someone can go back and answer that old question now.
Create two random scalar fields $f$ and $g$ and set the velocity to: $$ u=\nabla f \times \nabla g $$ which is guaranteed to be divergence free.
If you are in 2D and if you want more physical setups, I suggest you consider potential flows.
There are various ways to construct these potential flows which are always divergence free and which satisfy certain boundary conditions.
In the wikipedia article, you find the power laws, that describes flows around plates, edges, or in corners.
Also, there are other generating functions, that may give you, e.g., the flow around an idealized airfoil.
