Hello,
This is my first post on this forum and first of all I want to express my appreciation for this forum, since it has helped me a number of times with several problems related to LBM. For my master’s thesis I am working on the free-surface LBM, which has been popularized by Nils Thurey. Currently I am working on Regularized LBM (RLBM), as explained by Jonas Latt in his thesis [1]. I have two questions about RLBM, one concerns the derivation of RLBM and the other one concerns RLBM and adaptive time stepping.
First the derivation question. For the derivation of RLBM Jonas Latt combines two equations to get to the end result [2]. Since this was not entirely straight forward to me I decided to write out this derivation. This led me to the conclusion that f^{(1)} approx t_i / (2cs^4) * Pi^{neq}_ab and not f^{(1)} = t_i / (2cs^4) * Pi^{neq}_ab as Jonas Latt in his paper states. I have wrote the derivation in Latex and placed it online at http://home.kpn.nl/riji99ou/derivation.pdf (tex-file can be found at http://home.kpn.nl/riji99ou/derivation.tex if you want to add/change anything). Am I making a mistake in my derivation/reasoning?
Secondly, RLBM and adaptive time stepping. According to Jonas Latt’s thesis, we should scale the variables rho, u, and Pi^(1) appropriately when the time step changes. Then we can recompute the distribution functions with the RLBM procedure. As a checking mechanism we can recompute the velocity and density from the rescaled distribution functions, these should be equal to the rescaled velocity and density. However, they’re not in my case, the difference is of order 10^-3. I have posted my code below and perhaps/probably I am doing something wrong. Help would greatly be appreciated.
N.B.: In his thesis Jonas Latt converts u to a dimensionless formulation differently in equation 4.1 and 4.11. I think 4.1 is correct and 4.11 not. Correct me if I am wrong.
[1] J. Latt, Hydrodynamic limit of lattice Boltzmann equations, http://archive-ouverte.unige.ch/vital/access/services/Download/unige:464/THESIS
[2] J. Latt and B. Chopard, Lattice Boltzmann method with regularized pre-collision distribution functions, Mathematics and Computers in Simulation, Volume 72, Issues 2-6.
void LBMsim::regularizedChangeTimeStep(float dtNew_p) {
//Determine average density
double rhoTotal = 0.0;
int numFluid = 0;
for (int y = 1; y < SIZE - 1; ++y) {
for (int x = 1; x < SIZE - 1; ++x) {
if (flags[x][y] == IS_FLUID) {
rhoTotal += rho[x][y];
++numFluid;
}
}
}
double rhoAvg = rhoTotal / numFluid;
//Scale global parameters
double st = dtNew_p / dt_p;
dt_p = dtNew_p;
gravity_lb *= (st * st);
double nuNew_lb = nu_p * dtNew_p/(dx_p*dx_p); //New viscosity
tau = nuNew_lb / cs2 + 0.5;
//Scale cell parameters
for (int y = 1; y < SIZE - 1; ++y) {
for (int x = 1; x < SIZE - 1; ++x) {
if (flags[x][y] == IS_FLUID || flags[x][y] == IS_INTERFACE) {
double eq[9];
double neq[9];
double piNeq_xx = 0.0, piNeq_yy = 0.0, piNeq_xy = 0.0;
for (int i = 0; i < 9; ++i) {
double fi = f[next][x][y][i];
//Compute equilibrium distribution function
eq[i] = fEq(i, rho[x][y], ux[x][y], uy[x][y]);
//Compute off-equilibrium distribution function
neq[i] = fi - eq[i];
//Compute off-equilibrium momentum flux tensor
piNeq_xx += neq[i] * ex[i] * ex[i];
piNeq_yy += neq[i] * ey[i] * ey[i];
piNeq_xy += neq[i] * ex[i] * ey[i];
}
//Scale off-equilibrium momentum flux tensor
piNeq_xx *= st;
piNeq_yy *= st;
piNeq_xy *= st;
double f1[9];
for (int i = 0; i < 9; ++i) {
//Compute first-order term of distribution function
double Qixx = ex[i] * ex[i] - cs2;
double Qiyy = ey[i] * ey[i] - cs2;
double Qixy = ex[i] * ey[i];
f1[i] = w[i] / (2.0 * cs2 * cs2) * (Qixx * piNeq_xx + Qiyy * piNeq_yy + 2.0 * Qixy * piNeq_xy);
//Replace f by regularized counterpart
f[next][x][y][i] = eq[i] + f1[i];
}
//Scale velocity and density
double rhoOld = rho[x][y];
ux[x][y] *= st;
uy[x][y] *= st;
rho[x][y] = st * (rhoOld - rhoAvg) + rhoAvg;
}
}
}
//Check
for (int y = 1; y < SIZE - 1; ++y) {
for (int x = 1; x < SIZE - 1; ++x) {
if (flags[x][y] == IS_FLUID) {
double rhoCell = 0.0, uxCell = 0.0, uyCell = 0.0;
for (int i = 0; i < 9; ++i) {
double fi = f[next][x][y][i];
rhoCell += fi; uxCell += ex[i] * fi; uyCell += ey[i] * fi;
}
uxCell /= rhoCell;
uyCell /= rhoCell;
//The following output shows difference of order 10^-3
cout << uxCell - ux[x][y] << " ";
cout << uyCell - uy[x][y] << " ";
cout << rhoCell - rho[x][y] << "\n";
}
}
}
exit(0);
}
