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] = omegaE*ro[p] + omegaEj*j_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] = omegaXX*meq[p][9];
meq[p][11] = (pow(jy, 2) - pow(jz, 2)) / rho0;
meq[p][12] = omegaXX*meq[p][11];

meq[p][13] = jx

*jy / rho0;*

meq[p][14] = jyjz / rho0;

meq[p][14] = jy

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.