Hi all,
As we know the eigenvalues of S in MRT Lattice Boltzmann model are all between 0 and 2 to maintain linear stability. but I don’t know how can I choose or calculate every part of S matrix to keep the program stable. I will appreciate if you help me.
thanks in advance,
Arman
Hi Arman,
I don’t think there is a general rule for setting all the relaxation times in the collision matrix. Moreover, there is not one MRT model (at least not in my opinion) - that is, the basis is not unique.
If you are using an operator designed in the spirit of Lallemand, Luo, D’humieres and others then perhaps there are some things to remember. The relaxation time associated with the “energy” mode can be used to satisfy the no-slip boundary condition, in the same way as in the TRT model. Although, to the best of my knowledge, there is not a “magic parameter” in general, perhaps you can use the relationship between the relaxation times for Poiseuille flow in the TRT model as a guide (see work by Ginzburg, and also recent work by Luo et al).
As for the other times…well, there is some guidance in Lallemand and Luo (1998, I think) but no definite optimal choice. It may be a case of trial and error to find the most stable set. I would normally argue that since these times aren’t assocaited with the hydrodynamics then we want them to relax to their equilibrium straight away (ie, which probably means setting them to dt, or 1 in lattice units, but it depends on your model).
It is perhaps worth remembering that even if your choice of the kinetic relaxation times is not optimal (or you can’t tell if it is optimal or not), it still should be a lot more stable than BGK.
Good luck!
ps
Li-Shi gave a talk at our workshop in 2010 and his slides, which say a little bit about relaxation times in collision operators, can be found here (click on his name): http://www.maths.ox.ac.uk/groups/occam/forthcoming-events/occam-lattice-boltzmann-workshop
Dear pleb01, Really thanks for your advice, actually, my method is MRT LBE . I have simulated a fluid in 2D lid driven cavity (D2Q9) easily and my findings was verified,as well as, by changing the eigenvalues of S the result would be better but I ignored to get the maximum accurate achievement. In 3D (D3Q19), when I change the collision’s eigenvalues the results are entirely different from previous results even sometimes the program completely crash. I use the equations which is given by Lallemand, Luo, D’humieres (2002) and Du Rui (2010) unfortunately none of them have explanation how they choose and why. I have found some papers which are mentioned how to choose the collision parameters in 2D but they are not helpful for my problem. Do you have any experiment on MRT LBE D3Q19?
Hi,
I have never used an MRT D3Q19 model so I have no experience with the choice of relaxation times. However, the general guidance for D2Q9 lattice should be transferable.
This may sound obvious but make sure
a) if you use the same parameters as in the papers you quoted, you get the same results;
b) if you use “safe” values for the relaxation times (for example, set the kinetic ones equal to one) that you get results in agreement with benchmark simulations;
c) if you set ALL relaxations to be the same then you get the BGK results.
If these tests fail then I would guess you have a bug in your code. I don’t know for sure, of course, but since the MRT scheme should either be better than BGK (at least if you stay away from the stability boundaries), having all equal relaxation times is equivalent to the BGK case, and I can’t think of a reason why you can’t do cavity simulations with 3D LBE (someone please correct me if I am wrong) then this seems the most likely option.
My guess is you have a problem with the boundary conditions. Make sure you double check this. If you are using bounce back then make sure this is implemented correctly. If you are use Zou and He then remember this is not straight forward to use in 3D, and if you are using other methods then again double check everything. Note also that the cavity flow has a singularity in the stress at the top corners. This can be particularly important if you are using a method that implements boundary conditions by looking at derivatives of the velocity field.
Good luck!
ps
In my opinion, “MRT” doesn’t exclusively refer to the collision operators of lallemand and the other people we had already mentioned. This is because you can construct different matrices (with different eigenvectors, not just eigenvalues).To me, MRT is anything that uses more then one relaxation time…but I seem to be alone in thinking this! 
Hi, I think u are right I must find my problem on the other part. I will simulate 3D Lid-driven cavity base on Lallemand, Luo, D’humieres paper to get the same result, it makes me to be sure there is no problem in B.C. today I set wrong boundary condition for 2D lid-driven cavity and I obtained the same behavior of my main program. as well as, I mixed different B.Cs for different walls and I compared my result with Ghia’s result, my finding shows using Zou and He scheme for top wall and bounce back scheme for other walls are more accurate for this case (plz, notify me if I am wrong) however it requires more careful attention and academic research. I have just uploaded my matlab codes (2D lid driven cavity) in matlabcentral website. but I need at least one day to get a consent (I will share download address in next post).
you mentioned Zou and He in 3D is not simple, I know that we should modify the distribution functions for different walls but is there other difficulty?
Boundary conditions are one of the most common sources of error in LB, in my experience, and they can get complicated in 3D (because you have 19 distribution functions to be careful of.
The LB method gives results in good agreement with Ghia’s, even when the Reynolds number is 7500, so this is always a good way to test your code (in 2D). You may find that you can not reach high Reynolds numbers with the Zou and He method using the same grid resolution (257*257 nodes, lid velocity=0.1, in lattice units, is the standard). I don’t know for sure because I’ve not done it, but Zou and He could have a problem at high Re because of the corners, and because they can produce unwanted sound waves coming from the walls (they don’t guarantee mass conservation). Like I say, I may well be wrong, but be aware.
I’ve never implemented Zou and He in 3D but I know you have to be careful. You need to have enough information (or perhaps not too much) to find all the components of the velocity and also the density from the known distributions (those going into the wall), and make sure the equations for finding the f_i are linearly independent. I think this is possible, but might get a bit messy. People have studied this so a quick google search will give you more information than I can (there may even be an article on the LBwiki page). Also, if you are using a combination of bounce-back and Zou and He then you may have to think about which to use at the top corners and whether or not you can use half-way or full-way bound back on the non-moving walls (remember that Zou and H is applied exactly on the nodes, whereas to get second order results with bound-back you need to have the wall half way between grid points…this may not be an issue at all, and I’m sorry if I’m making things sounds more complicated than they are, but I thought I’d mention it just in case).
Good luck!
your points are very helpful, again thank you so much.
hi again,
my download source code is available in this post (2D Lid Driven Cavity Matlab).
I hope, it is useful for someone who would like to learn LBM (special MRT LBM).
Arman
Dear pleb01,
I have written 3D lid driven cavity program in Matlab. My program works well when I run it with BGK collision term however it dose not work in MRT. I set everything according to the papers that I mentioned above post.That’s weird, I thought it should be easy but it isn’t. Unfortunately, I don’t know which part I did wrong and I can not find my mistake or mistakes, could you check my codes?
If you let me know your email address I will send it to you.
my email address is sarman2@live.utm.my
thanks in advance,
Arman
Hi,
I’m afraid I don’t understand Matlab because I have never used it. However, I’ll reiterate some points that were raised before.
Make sure you get the BGK results if you set all relaxation parameters to be equal in your MRT code. If this is not the case then you can be quite confident that you’ve made a mistake somewhere in the collision operator.
If your MRT code works in 2D but not in 3D, and if you are 100% sure that your collision matrix is correct, then I would guess that the boundary conditions are causing the problem. This is still a little weird if the BGK case works but I wouldn’t be surprised if using Zou and He boundaries are a bit awkward for 3D MRT models (in general, Zou and He is a bit inconsistent). You might want to think about using just bounce-back, and setting the top wall moving by setting u=Ulid in the equilibrium function
The “benchmark” LB cavity paper is this one http://arxiv.org/abs/comp-gas/9401003. It’s 2d, but it may be of interest if you haven’t read it.
Good luck!
for testing the program not only I set all relaxation parameters equal but also I plot Images of the cavity in different faces so I am sure there is a problem on my codes. honestly, at first I tried to use Zou and He BC but it is quit complicated special when I added the correction terms. then I used the standard equilibrium distribution function at previous time for moving wall.Also bounce back BC is used for the rest walls. my program’s algorithm is :
1-Calculation of macroscopic ?ow quantities and equilibrium distribution
rho = sum(f) , velocity = sum ( e*f ) / rho , feq(for moving wall) = … and so on
2- Collision term
3- Propagation step
4- full Bounce back
5- go to step 1
For MRT I just change the Collision term. I checked collision matrix for several time. if you think, I should change the major components of this algorithm please tell me. when I can’t find the bug it makes me very angry.
Anyway thanks for your helping and prompt respond.
I spend most of my life being angry with my codes! Sometimes I am only a few bugs away from throwing my computer out of the window.
Are you saying that if you set all relaxation times to be the same (all of them equal to 1/tau, say) then your code does’t work? But it does if you use standard BGK? Well, this probably means you’ve found where your problem is because if all the relaxation times are equal then you have (or should have) nothing more that he BGK operator (albeit expressed in a more complicated way).
If you can’t see the bug in your code then do what us mathematicians do and embrace the pen and paper! Work out the collision operator by hand and make sure it is the same as in the paper and in the code (there is a chance there is a mistake in the paper, remember). Why not multiply out your collision matrix to see what you get? This will be very messy in 3D so it might be an idea to use something like Maple to do the matrix arithmetic. If all relaxation times are then all set equal (1/tau) you should end up with nothing more than the BGK operator.
After a lot of anger and frustration, it should work…
hehe ohh yeeeh I hope so.
I worked on my codes entire night and I checked every part of program and values. fortunately, At last I found a wrong negative sign in the collision term which is make the program works terrible :D.
for someone who wants to know how the polynomial eigenvectors in MRT and TRT come I offer to study ‘‘Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation’’ by Irina Ginzburg in 2005. Also there is a little explanation for how we choose eigenvalues of collision matrix for D3Q19 and D3Q15 in “three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase flow” by Kannan, Premnath, Abraham in 2007.
please share if u know something more.
For 3D MRT model, you can also refer to the following one.
d’Humières, D., Ginzburg, I., Krafczyk, M., Lallemand, P., and Luo, L.S. Multiple-relaxation-time lattice Boltzmann models in three dimensions. Philosophical Transactions of the Royal Society A: Mathematical, Physical & Engineering Sciences 360, 437-451, 2002.