Diffusion with TRT collision operator
OpenLB – Open Source Lattice Boltzmann Code › Forums › on Lattice Boltzmann Methods › General Topics › Diffusion with TRT collision operator
 This topic has 2 replies, 2 voices, and was last updated 5 years, 3 months ago by tobit.

AuthorPosts

November 11, 2018 at 5:29 pm #2005tobitParticipant
Hey everybody!
After having programmed a D2Q9 and D3Q19 model with BGK, TRT collision operators and the corresponding Smagorinsky model with C++ and OpenMP (I validated the code with several benchmark scenarios), I wanted to add a basic advectiondiffusion equation that accounts for the passive scalar transport of a mass fraction Y.
Everything worked fine with the BGK operator but I ran into issues with the two relaxation time TRT collision operator. Weirdly enough my simulated physics change as I adjust the second relaxation time. If I enforce omega = omega+ the result is correct but as soon as I deviate from this, for example by enforcing the magical parameter, the model behaves in a way I can’t yet understand: If I set omega higher than omega+ the diffusion happens slower and if I set it lower the diffusion happens faster than it should. The contribution of the antisymmetric part becomes comparably large but I can’t see my mistake.I then took out the advection part, now the code is basically diffusion only but still exhibits the same behaviour. Here the relevant code segment:
const double tau_d = (6.0*diffusivity + 1.0) / 2.0;
const double magic_param = 1.0 / 6.0;const double omegap = 1/tau_d;
const double omegam = (1.0/omegap – 1.0/2.0) / (magic_param + 1.0/2.0*( 1.0/omegap – 1.0/2.0));for(unsigned int y = 0; y < NY; ++y)
{
for(unsigned int x = 0; x < NX; ++x)
{
/// temporary populations (load populations from main memory)
double gt0 = g0[D2Q9_ScalarIndex(x,y)];
double gt1 = g1[D2Q9_FieldIndex(x,y,1)];
double gt2 = g1[D2Q9_FieldIndex(x,y,2)];
double gt3 = g1[D2Q9_FieldIndex(x,y,3)];
double gt4 = g1[D2Q9_FieldIndex(x,y,4)];
double gt5 = g1[D2Q9_FieldIndex(x,y,5)];
double gt6 = g1[D2Q9_FieldIndex(x,y,6)];
double gt7 = g1[D2Q9_FieldIndex(x,y,7)];
double gt8 = g1[D2Q9_FieldIndex(x,y,8)];/// microscopic to macroscopic: mass fractions
double Ymt = gt0 + gt1 + gt2 + gt3 + gt4 + gt5 + gt6 + gt7 + gt8;
Ym[D2Q9_ScalarIndex(x,y)] = Ymt;/// equilibrium distribution: only diffusion (no velocity!)
double geq1 = Ymt*wc;
double geq2 = Ymt*wc;
double geq3 = Ymt*wc;
double geq4 = Ymt*wc;
double geq5 = Ymt*wd;
double geq6 = Ymt*wd;
double geq7 = Ymt*wd;
double geq8 = Ymt*wd;/// symmetric and antisymmetric populations
double g13p = (gt1 + gt3 – (geq1 + geq3))/2;
double g13m = (gt1 – gt3 – (geq1 – geq3))/2;
double g24p = (gt2 + gt4 – (geq2 + geq4))/2;
double g24m = (gt2 – gt4 – (geq2 – geq4))/2;
double g57p = (gt5 + gt7 – (geq5 + geq7))/2;
double g57m = (gt5 – gt7 – (geq5 – geq7))/2;
double g68p = (gt6 + gt8 – (geq6 + geq8))/2;
double g68m = (gt6 – gt8 – (geq6 – geq8))/2;/// neighbouring cells coordinates (assume periodic boundaries
unsigned int xp = (x + 1) % NX;
unsigned int yp = (y + 1) % NY;
unsigned int xm = (NX + x – 1) % NX;
unsigned int ym = (NY + y – 1) % NY;/// collision & streaming
g0[D2Q9_ScalarIndex(x,y)] = gt0 + omegap*(Ymt*w0 – gt0);
g2[D2Q9_FieldIndex(xp, y, 1)] = gt1 – omegap*g13p – omegam*g13m;
g2[D2Q9_FieldIndex(xm, y, 3)] = gt3 – omegap*g13p + omegam*g13m;
g2[D2Q9_FieldIndex(x, yp, 2)] = gt2 – omegap*g24p – omegam*g24m;
g2[D2Q9_FieldIndex(x, ym, 4)] = gt4 – omegap*g24p + omegam*g24m;
g2[D2Q9_FieldIndex(xp, yp, 5)] = gt5 – omegap*g57p – omegam*g57m;
g2[D2Q9_FieldIndex(xm, ym, 7)] = gt7 – omegap*g57p + omegam*g57m;
g2[D2Q9_FieldIndex(xm, yp, 6)] = gt6 – omegap*g68p – omegam*g68m;
g2[D2Q9_FieldIndex(xp, ym, 8)] = gt8 – omegap*g68p + omegam*g68m;
}
}I am thankful for every hint!
RegardsNovember 12, 2018 at 2:03 pm #2960mgaedtkeKeymasterHello tobit,
when you do TRT for the advection diffusion equation, make sure that omegam is the omega, which is calculated from the diffusion coefficient, and omegap is calculated with the magic parameter. For more information see http://alexkuzmin.com/wpcontent/uploads/2011/05/stat_stability.pdf
I hope this will help you fix the problem.
With best regards,
MaxNovember 12, 2018 at 3:00 pm #2961tobitParticipantThanks a lot, that solved my problem!
Funny enough I had already tried something similar: I simply swapped omegap and omegam but doing so I also changed the rest node relaxation rate and therefore I still had a similar problem.
So for the advectiondiffusion equation the rest node is still relaxed with omegap but in this case omegap is calculated with the magic parameter whereas omegam is calculated with the LBM diffusion coefficient.
Thanks again! 
AuthorPosts
 You must be logged in to reply to this topic.