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!

  • $\begingroup$ One should be able to reproduce basic vector properties in vector fields if the calculation is correct, e.g., generate a vector field from a scalar potential $v=grad(\phi)$, and make sure $curl(v)=0$; or make $v=curl(A)$ and see if $div(v)=0$. It is nontrivial to satisfy those properties identically, but there is a class of methods called "mimetic" that is designed for those things. $\endgroup$ Commented Dec 6, 2022 at 3:04

1 Answer 1


One way to make your field divergence free is via the Leray Projection.

The underlying idea is to identify that part of your vector field which is non-divergence free, and then subtract that from your initial field.

Another approach is to choose the x and y components of your vector field, and choose the third component via the icompressibility condition.

  • $\begingroup$ For the second case. ux and uy are chosen to be certain random fields. Now one can take the gradients off those fields by using central difference. How would one go ahead and calculate the integral over these 3D matrixes added together? $\endgroup$
    – DvB
    Commented Dec 7, 2022 at 16:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.