Hello to everyone,
I am trying to understand how the multiple relaxation model works, but I feel overdosed of info, which I don’t really understand.
Does anyone have a plain code of multiple relaxation where the collision routine is well explained to just try to assimilate what’s going on with that huge matrix?
Thanks in advance.
Hi,
This is just a quick comment about MRT models.
They are indeed a source of much confusion, but although a full analysis is very difficult, the implementation and practical usage is actually not as bad as it first seems. In particular, the matrix isn’t huge. For D2Q9 it’s just 9x9. It’s this size for every grid point, but it’s fully local. Moreover, we do not need to do any matrix inversions in the code.
The idea behind MRT is to do the streaming step on the f_i and the collisions on the moments of f_i. That is, streaming is done in the usual way, but instead of doing collisions on f_i, we do them on its moments.
The D2Q9 LBE has 9 independent moments. Three of these are the conserved quantities - density and the two components of momentum. Three more are given by the components of the momentum flux, Pi_{ab}=sum_i(f_i*c_{ia}*c_{ib}). The remaining 3 do not appear in the Navier-Stokes equations.
Let’s write the moments of f_i as m. Note that m is a vector with 9 components - one for each moment. We can write the relationship between the f_i and its moments as
m=M*f
where f is a 9-dimesnional vector f=<f0,f1,…f8> and M I will call the transformation matrix. M is a 9*9 matrix and each row corresponds to a different moment. For example, the first row is made of of just 1s, because this is how we define density (sum(f_i)). The next row is for the x-component of momentum, ie <0,1,0,-1,0,1,-1,-1,1>, and the third for the y-component.
There are different ways of constructing the other rows, but they should form a basis. The most common approach is to use Gram-Schmidt orthogonlisation. This is how d’Humieres did it, and Lallemand and Luo, for example. Dellar, on the other hand, constructs the matrix using the first 6 (hydrodynamic) moments mentioned above, and then extends the Hermite polynomial analogy to find the remaining three. Please see those papers for further details.
The MRT operator is sometimes written as C=M^-1SM(f_i-f0_i), where f0_i are the equilibria and S is a diagonal matrix of collision times. Note that by the transformation above, if the 9 entries of S are all 1/tau then we get back to BGK.
So, one way to implement the MRT would be to find C by finding the inverse of M, multiplying the matrices, etc. But this is expensive and pointless.
Now, since the relationship between m and f is a linear, invertible, transformation/mapping, we can easily convert between the two. So, instead of colliding f_i, we can relax it’s moments:
m’=m-1/tau_j*(m-m0)
where m’ means the post collision m, tau_j is the relaxation time of the jth moment, and m0 is the equilibrium of the moment. For example, if you were using Dellar’s MRT which uses the xx component of Pi, you’d say
Pi’_xx=Pi_xx-1/tau*(Pi_xx-rho/3-u^2)
Once we have found the post-collision m, we can find the post-collision f_i by our linear transformation:
f’=M^-1*m’
and then do the streaming. Note that we do not need need to invert the matrix M at every timestep. We just need to know what f_i looks like in terms of its moments. Then we simply write this into our code, using m’.
The simplest MRT is probably the one which sets the ‘ghost’ moments to be zero and chooses their collision times so that they decay straight to this equilibrium. In Dellar’s MRT, for example, this means setting tauJ=tauN=1/2. This means that the “ghost” moments never actually appear in the code. This then reduces to the model used by Tony Ladd in the 90s and was more recently given the name “regularized LBE”.
I hope this is of some help, even if only to say that we do not do any matrix inversions in the code. The work by the author’s I mentioned will help fill in the gaps.
Good luck!
Hi guys,
I wrote a FORTRAN code for Poiseuille flow past a cylinder SRT-LBM D2Q2 and it works properly. Then I removed the cylinder and implemented MRT-LBM instead of SRT and I’ve got good results. But once I tried to insert the cylinder again with MRT, it doesn’t work. If anybody is willing to help me I can send him my code.
Thank you in advance