Hi Adithya,

the best way to impose an initial velocity field that is divergence-free is to impose that every single synthetic eddy satisfies du/dx+dv/dy+dw/dz = 0. You can then simply sum their individual velocity fields at every point of your domain, when you generate multiple anisotropic eddies with random location. There are many examples of such model eddies such as Townsend's eddy (you can find the precise form in many turbulence books).

Hope the answer helps!

There is a fairly simple method to construct all possible divergence-free eddies.

Define a vector field.  This can be absolutely any smooth function you desire (as long as it is differentiable).  This means pick a function (analytic or discrete) for the x,y, and z components of a vector.

Define your eddy velocity field to be the curl of this vector.

These eddies are always incompressible.

FFT can do this job in Cartesian coordinates.  I have a piece of Matlab code from some paper in JPO. It can decompose any flow field into the sum of a non-divergent part (curl of something) and a divergent part (gradient of something).

nft = 2^nextpow2(max(Ny,Nx)*1.2);% min 20% zero padding, Nx and Ny are grid points of the 2D flow field.

nf2 = nft/2;

kv = (-nf2:(nf2-1))/nft;

thetamn= fftshift(atan2(kv(ones(nft,1),:)',kv(ones(nft,1),:)));

% PS decomposition

%VXm_geostr and VYm_geostr are U and V components of the velocity vector field

vkdiv = fft2(VXm_geostr,nft,nft).*cos(thetamn)+fft2(VYm_geostr,nft,nft).*sin(thetamn);

vav = real(ifft2(vkdiv.*cos(thetamn)));

Pu = vav(1:Ny,1:Nx);

vav = real(ifft2(vkdiv.*sin(thetamn)));

Pv = vav(1:Ny,1:Nx);

%(Pu, Pv) is the divergent part of original (U,V) field

vkdiv = fft2(VYm_geostr,nft,nft).*cos(thetamn)-fft2(VXm_geostr,nft,nft).*sin(thetamn);

vav = real(ifft2(-vkdiv.*sin(thetamn)));

Su = vav(1:Ny,1:Nx);

vav = real(ifft2(vkdiv.*cos(thetamn)));

Sv = vav(1:Ny,1:Nx);

clear VXm_geostr VYm_geostr vkdiv vav;

%(Su, Sv) is the non-divergent part of original (U,V) field.

Hope it helps..