I'm trying to build a DNS so simulate non-Newtonian turbulent channel flow using the Lattice Boltzmann Method. Most of the papers I use initialise the system using the analytical solution of the laminar flow (which describes a parabolic flow profile in the flow direction) and ad divergent free random fluctuations to it.
I tried making these divergent free random fluctuations firstly by using the formula proposed in this post $$ u_{random} = \nabla f \times \nabla g $$ where f and g are random fluctuating scalar fields. Sadly this implementation using numpy did not generate a divergence free field
f = np.random.rand(Nx, Ny, Nz)
g = np.random.rand(Nx, Ny, Nz)
gradf = np.gradient(f)
gradg = np.gradient(g)
u = np.cross(gradf, gradg, axis=0)
since the maximum divergence here is around 1. My first question is if this is the correct method and if so how to make it more precise.
Secondly lots of posts refer to different methods to initialise turbulent DNS. This post proposes an energy spectrum function that I believe has to be implemented in the Rogallo's procedure which is further discussed here. For me simulation time is not such a restriction and since the Ragallo method seems way more difficult I was wondering if just using a simple divergence free u field is enough.
Hoop someone can help me further!