LBM 3D MRT

Lattice Boltzmann Method LBM: Algorithms
#1

Hello everyone!
I am new to LBM and this forum, so I apologize if I lack general understanding of things.
I wrote a code SRT (BGK) D3Q19 LBM with Zou He Boundary conditions for simulating Duct flow and it works fine.
In order to increase the accuracy of my output I’ve decided to implement MRT (Multi-Relaxation Time scheme). The problem I am having is that the symmetry is no longer there for a square duct. My understanding of MRT (from papers) is that we need to transfer all of the distribution fuctions to a momentum space (by multiplying it by collision matrix, that I have derived uniquely for my judiciously chosen arrangment of e vectors) and then apply the collision step with different relaxation (1/tau or omega) for each of the distribution functions (in this momentum space) and then transfer back to continue with streaming step (multiplying by inverse matrix M-1).

What I did to find the error but wasn’t able to:
I made all of the omegas equal in the S matrix but couldn’t recover the BGK results, however when I used f_eq to find m_eq(m_eq=M*f_eq) instead of their own formulas I was able to get the same results as BGK.
Here is that piece of my code that I believe should be working but doenst seem like it:

for (int p = 0; p < num_cells; p++) {
jx = ro[p] * ux[p]; jy = ro[p] * uy[p]; jz = ro[p] * uz[p]; j_squ = pow(jx, 2) + pow(jy, 2) + pow(jz, 2);
meq[p][0] = ro[p];
meq[p][1] = -11 * ro[p] + 19 * j_squ / rho0;
meq[p][2] = omegaEro[p] + omegaEjj_squ / rho0;
meq[p][3] = jx;
meq[p][4] = -2 * jx / 3;
meq[p][5] = jy;
meq[p][6] = -2 * jy / 3;
meq[p][7] = jz;
meq[p][8] = -2 * jz / 3;
meq[p][9] = (2 * pow(jx, 2) - (pow(jy, 2) + pow(jz, 2))) / rho0;
meq[p][10] = omegaXXmeq[p][9];
meq[p][11] = (pow(jy, 2) - pow(jz, 2)) / rho0;
meq[p][12] = omegaXXmeq[p][11];
meq[p][13] = jxjy / rho0;
meq[p][14] = jyjz / rho0;
meq[p][15] = jx*jz / rho0;
meq[p][16] = 0;
meq[p][17] = 0;
meq[p][18] = 0;

	}

I am not sure if this is what’s causing the error.
But I am 100% sure that I got the right collision matrix, as I re-checked 5 times.

f0 f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 f12 f13 f14 f15 f16 f17 f18
0 1 -1 0 0 0 0 1 -1 -1 1 1 -1 -1 1 0 0 0 0
0 0 0 1 -1 0 0 1 1 -1 -1 0 0 0 0 1 -1 -1 1
0 0 0 0 0 1 -1 0 0 0 0 1 1 -1 -1 1 1 -1 -1

double M[19][19] = {
{ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ -30, -11, -11, -11, -11, -11, -11, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8 },
{ 12, -4, -4, -4, -4, -4, -4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1, 0, 0, 0, 0 },
{ 0, -4, 4, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1, 0, 0, 0, 0 },
{ 0, 0, 0, 1, -1, 0, 0, 1, 1, -1, -1, 0, 0, 0, 0, 1, -1, -1, 1 },
{ 0, 0, 0, -4, 4, 0, 0, 1, 1, -1, -1, 0, 0, 0, 0, 1, -1, -1, 1 },
{ 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, 1, -1, -1, 1, 1, -1, -1 },
{ 0, 0, 0, 0, 0, -4, 4, 0, 0, 0, 0, 1, 1, -1, -1, 1, 1, -1, -1 },
{ 0, 2, 2, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -2, -2, -2, -2 },
{ 0, -4, -4, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, -2, -2, -2, -2 },
{ 0, 0, 0, 1, 1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0 },
{ 0, 0, 0, -2, -2, 2, 2, 1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0 },
{ 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0 },
{ 0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 1, -1, 1, 1, -1, 0, 0, 0, 0 },
{ 0, 0, 0, 0, 0, 0, 0, -1, -1, 1, 1, 0, 0, 0, 0, 1, -1, -1, 1 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, -1, -1, -1, -1, 1, 1 },
};
double S[19] = { 0, 1.19, 1.4, 0, 1.2, 0, 1.2, 0, 1.2, omega, 1.4, omega, 1.4, omega, omega, omega, 1.98, 1.98, 1.98 };

I’d really appreciate your help.

#2

I fixed my code, the problem was in the collision matrix, on the sixth time recheking it I saw a mistake, I was not using the right arrangement of vectors, took me 20 hours to fix((

Now that my MRT code can give BGK results, I have changed my S matrix, and now I have a problem, when I set omega = 0 for conserved parameters, it goes to infinity (, not sure why though
double S[19] = { 0, 1.19, 1.4, 0, 1.2, 0, 1.2, 0, 1.2, omega, 1.4, omega, 1.4, omega, omega, omega, 1.98, 1.98, 1.98 };
I mean these 0s

#3

I will admit, I don’t know why you get infinity when you set the relaxation frequencies for conserved properties (rho, jx, jy and jz) to 0, but setting these to 1 should definitely still work. Are you implementing any forces in your MRT scheme?

Regards,

Michael

#4

I think you have a mistake in meq[p][1] = -11 * ro[p] + 19 * j_squ / rho0;
you should multiply by rho0 instead of division.

#5

e
int e[Q][3] = {
{ 0,0,0},//0
{ 1,0,0},//1
{-1,0,0},//2
{0, 1,0},//3
{0,-1,0},//4
{0,0, 1},//5
{0,0,-1},//6
{ 1, 1,0},//7
{-1, 1,0},//9
{ 1,-1,0},//8
{-1,-1,0},//10
{ 1, 0,1},//11
{-1, 0,1},//12
{ 1,0,-1},//13
{-1,0,-1},//14
{0, 1, 1},//15
{0,-1, 1},//17
{0, 1,-1},//16
{0,-1,-1}};//18
[\code]
for
M=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ;
-30 -11 -11 -11 -11 -11 -11 8 8 8 8 8 8 8 8 8 8 8 8;
12 -4 -4 -4 -4 -4 -4 1 1 1 1 1 1 1 1 1 1 1 1;
0 1 -1 0 0 0 0 1 -1 1 -1 1 -1 1 -1 0 0 0 0;
0 -4 4 0 0 0 0 1 -1 1 -1 1 -1 1 -1 0 0 0 0;
0 0 0 1 -1 0 0 1 1 -1 -1 0 0 0 0 1 -1 1 -1;
0 0 0 -4 4 0 0 1 1 -1 -1 0 0 0 0 1 -1 1 -1;
0 0 0 0 0 1 -1 0 0 0 0 1 1 -1 -1 1 1 -1 -1;
0 0 0 0 0 -4 4 0 0 0 0 1 1 -1 -1 1 1 -1 -1;
0 2 2 -1 -1 -1 -1 1 1 1 1 1 1 1 1 -2 -2 -2 -2;
0 -4 -4 2 2 2 2 1 1 1 1 1 1 1 1 -2 -2 -2 -2;
0 0 0 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 0 0 0 0;
0 0 0 -2 -2 2 2 1 1 1 1 -1 -1 -1 -1 0 0 0 0;
0 0 0 0 0 0 0 1 -1 -1 1 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 -1 1;
0 0 0 0 0 0 0 0 0 0 0 0 1 -1 -1 1 0 0 0;
0 0 0 0 0 0 0 1 -1 1 -1 -1 1 -1 1 0 0 0 0;
0 0 0 0 0 0 0 -1 -1 1 1 0 0 0 0 1 -1 1 -1;
0 0 0 0 0 0 0 0 0 0 0 1 1 -1 -1 -1 -1 1 1];

#6

I fixed all of that(problem was in boundary conditions), now it is independent of those omegas (for conserved parameters). And my MRT version of the code seems working ok (it gives adequate results)
The only thing is that I am simulating flow through porous media and as you know permeability is dependent on relaxation parameter for SRT.
In contrast, I have read in the literature that permeability should be independent of free parameter for MRT, but only if you carefully choose your omegas in the S matrix. With parameters that I wrote earlier S[19] = { 0, 1.19, 1.4, 0, 1.2, 0, 1.2, 0, 1.2, omega, 1.4, omega, 1.4, omega, omega, omega, 1.98, 1.98, 1.98 }; it is still depends on omega somewhat. (however slightly less that for SRT)

So my questions is what set of omegas I should use?

P.S. I tried TRT as well and it gives consistent results at any omega, but I want MRT to work not TRT, as I am planning to implement it into multiphase Shan-Chen model. I know I might be wrong in some places, but your advice can be extremely helpful. Thanks

#7

Thanks for the reply

But it is not the case, as rho0 is set as unity so it doesnt matter if I divide or multiply by it