# LBM BGK

Hi all!

I just started with the LBM. I correctly simulated 2D Poiseuille flow with a pressure gradient between the inlet and outlet (rho_in rho_out = 1.1 and = 1.0). For this, I used two distribution functions f [lx] [ly] [q] and f_old [lx] [ly] [q].

In the public void initialize (), I initialized. And the function void collision_streaming_bgk (), in the collision, I used f [x] [y] [k] éffectuer calculation, and the streaming part and the boundary conditions (pressure, bounce-back) I used f_old [x] [y] [k]. And prior to recalculate the macroscopic variables (density rho, and velocity_x velocity_y), I made a swapping between f_old [x] [y] [k] and f [x] [y] [k].

I would like to now work with only one distribution function f [x] [y] [k], without using f_old [x] [y] [k], do you have any suggestions for me as that should proceed.

Thank you for help!!!

typedef double T;

T *** f, ***f_old; // LBM populations (old and new)
T *feq; // Equilibrium pouplation
T **rho; // Fluid density
T **velocity_x; // Fluid velocity x-component
T **velocity_y; // Fluid velocity y-component

// This function computes the equilibrium function from the fluid density and velocity.

void equilibrium(T density, T vel_x, T vel_y) {

T f1 = (T)3;
T f2 = (T)9/(T)2;
T f3 = (T)3/(T)2;

feq[0] = weights[0] * density * ((T)1 - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[1] = weights[1] * density * ((T)1 + f1 * vel_x + f2 * SQ(vel_x) - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[2] = weights[2] * density * ((T)1 + f1 * vel_y + f2 * SQ(vel_y) - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[3] = weights[3] * density * ((T)1 - f1 * vel_x + f2 * SQ(-vel_x) - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[4] = weights[4] * density * ((T)1 - f1 * vel_y + f2 * SQ(-vel_y) - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[5] = weights[5] * density * ((T)1 + f1 * ( vel_x + vel_y) + f2 * SQ( vel_x + vel_y) - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[6] = weights[6] * density * ((T)1 + f1 * (-vel_x + vel_y) + f2 * SQ(-vel_x + vel_y) - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[7] = weights[7] * density * ((T)1 - f1 * ( vel_x + vel_y) + f2 * SQ(-vel_x - vel_y) - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[8] = weights[8] * density * ((T)1 + f1 * ( vel_x - vel_y) + f2 * SQ( vel_x - vel_y) - f3 * (SQ(vel_x) + SQ(vel_y)));
}

// Allocate memory for the fluid density, velocity and Initialize the fluid density and velocity. Start with unit density and zero velocity.

void initialize() {

rho = new T*[lx];
velocity_x = new T*[lx];
velocity_y = new T*[lx];

for (int X = 0; X < lx; ++X) {

rho[X] = new T[ly];

velocity_x[X] = new T[ly];
velocity_y[X] = new T[ly];
}

for (int X = 0; X < lx; ++X) {

for (int Y = 0; Y < ly; ++Y) {

rho[X][Y] = 1.;

velocity_x[X][Y] = 0.;
velocity_y[X][Y] = 0.;
}
}

// Allocate memory for the populations and Initialise density distribution function with equilibrium fo zero density.

feq = new T[Q];
f = new T**[lx];
f_old = new T**[lx];

for (int X = 0; X < lx; ++X) {

f[X] = new T*[ly];
f_old[X] = new T*[ly];

for (int Y = 0; Y < ly; ++Y) {

f[X][Y] = new T[Q];
f_old[X][Y] = new T[Q];

for (int iD = 0; iD < Q; ++iD) {

feq[iD] = 0.;

f[X][Y][iD] = 0.;
f_old[X][Y][iD] = 0.;
}
}
}

for (int X = 0; X < lx; ++X) {

for (int Y = 0; Y < ly; ++Y) {

equilibrium(rho[X][Y], velocity_x[X][Y], velocity_y[X][Y]);

for (int iD = 0; iD < Q; ++iD) {

f[X][Y][iD] = feq[iD];
f_old[X][Y][iD] = feq[iD];
}
}
}
}

/// Princal Lattice Boltzmann Method LBM BGK: Computed Collision, streaming and boundaries conditions steps.

void collision_streaming_bgk() {

// Computation Collision, streaming and Boundaries conditions steps.

for (int X = 0; X < lx; ++X) {

int Xn = (X > 0)?(X-1):(lx-1);
int Xp = (X < lx-1)?(X+1):(0);

for (int Y = 0; Y < ly; ++Y) {

int Yn = (Y > 0)?(Y-1):(ly-1);
int Yp = (Y < ly-1)?(Y+1):(0);

// Compute equilibrium
// The pequilibrium population are computed.

equilibrium(rho[X][Y], velocity_x[X][Y], velocity_y[X][Y]);

// Compute new population
// This is the lattice boltzmann equation (combined collision and streaming) including external forcing.

f[X][Y][0] = f[X][Y][0] * ((T)1 - omega) + feq[0] * omega;
f[X][Y][1] = f[X][Y][1] * ((T)1 - omega) + feq[1] * omega;
f[X][Y][2] = f[X][Y][2] * ((T)1 - omega) + feq[2] * omega;
f[X][Y][3] = f[X][Y][3] * ((T)1 - omega) + feq[3] * omega;
f[X][Y][4] = f[X][Y][4] * ((T)1 - omega) + feq[4] * omega;
f[X][Y][5] = f[X][Y][5] * ((T)1 - omega) + feq[5] * omega;
f[X][Y][6] = f[X][Y][6] * ((T)1 - omega) + feq[6] * omega;
f[X][Y][7] = f[X][Y][7] * ((T)1 - omega) + feq[7] * omega;
f[X][Y][8] = f[X][Y][8] * ((T)1 - omega) + feq[8] * omega;

// Streaming steps

f_old[X][Y][0] = f[X][Y][0];
f_old[Xp][Y][1] = f[X][Y][1];
f_old[X][Yp][2] = f[X][Y][2];
f_old[Xn][Y][3] = f[X][Y][3];
f_old[X][Yn][4] = f[X][Y][4];
f_old[Xp][Yp][5] = f[X][Y][5];
f_old[Xn][Yp][6] = f[X][Y][6];
f_old[Xn][Yn][7] = f[X][Y][7];
f_old[Xp][Yn][8] = f[X][Y][8];

}
}

// Dirichlet Pressure Boundaries conditions
// Zou and He computed Boundaries conditions

// Inlet (South side)

for (int X = 0; X < lx; ++X) {

int Y = 0;

rho[X][Y] = rho_in;

velocity_x[X][Y] = (T)0;
velocity_y[X][Y] = rho_in - (f_old[X][Y][0] + f_old[X][Y][1] + f_old[X][Y][3] + (T)2 * (f_old[X][Y][4] + f_old[X][Y][7] + f_old[X][Y][8]));

f_old[X][Y][2] = f_old[X][Y][4] + (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][5] = f_old[X][Y][7] + (T)1/(T)6 * velocity_y[X][Y] - (T)1/(T)2 * (f_old[X][Y][1]-f_old[X][Y][3]);
f_old[X][Y][6] = f_old[X][Y][8] + (T)1/(T)6 * velocity_y[X][Y] - (T)1/(T)2 * (f_old[X][Y][3]-f_old[X][Y][1]);
}

// Outlet (East side)

for (int X = 0; X < lx; ++X) {

int Y = ly-1;

rho[X][Y] = rho_out;

velocity_x[X][Y] = (T)0;
velocity_y[X][Y] = -rho_out + (f_old[X][Y][0] + f_old[X][Y][1] + f_old[X][Y][3] + (T)2 * (f_old[X][Y][2] + f_old[X][Y][5] + f_old[X][Y][6]));

f_old[X][Y][4] = f_old[X][Y][2] - (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][7] = f_old[X][Y][5] - (T)1/(T)6 * velocity_y[X][Y] + (T)1/(T)2 * (f_old[X][Y][1]-f_old[X][Y][3]);
f_old[X][Y][8] = f_old[X][Y][6] - (T)1/(T)6 * velocity_y[X][Y] + (T)1/(T)2 * (f_old[X][Y][3]-f_old[X][Y][1]);
}

// Bounce-back
// Due to the presence of the rigid walls at x = 0 and x = lx-1, the population have to be bounced back.

// left wall (West side)

for (int Y = 1; Y < ly-1; ++Y) {

int X = 0;

velocity_x[X][Y] = (T)0;
velocity_y[X][Y] = (T)0;

rho[X][Y] = velocity_x[X][Y]+f_old[X][Y][0]+f_old[X][Y][2]+f_old[X][Y][4]+(T)2*(f_old[X][Y][3]+f_old[X][Y][6]+f_old[X][Y][7]);
f_old[X][Y][1] = f_old[X][Y][3]+(T)2/(T)3velocity_x[X][Y];
f_old[X][Y][5] = f_old[X][Y][7]+(T)1/(T)6
velocity_x[X][Y]+(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])+(T)1/(T)2velocity_y[X][Y];
f_old[X][Y][8] = f_old[X][Y][6]+(T)1/(T)6
velocity_x[X][Y]-(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])-(T)1/(T)2*velocity_y[X][Y];
}

// right wall (East side)

for (int Y = 1; Y < ly-1; ++Y) {

int X = lx-1;

velocity_x[X][Y] = (T)0;
velocity_y[X][Y] = (T)0;

rho[X][Y] = -velocity_x[X][Y]+f_old[X][Y][0]+f_old[X][Y][2]+f_old[X][Y][4]+(T)2*(f_old[X][Y][1]+f_old[X][Y][5]+f_old[X][Y][8]);
f_old[X][Y][3] = f_old[X][Y][1]-(T)2/(T)3velocity_x[X][Y];
f_old[X][Y][7] = f_old[X][Y][5]-(T)1/(T)6
velocity_x[X][Y]-(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])-(T)1/(T)2velocity_y[X][Y];
f_old[X][Y][6] = f_old[X][Y][8]-(T)1/(T)6
velocity_x[X][Y]+(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])+(T)1/(T)2*velocity_y[X][Y];
}

// Corners nodes treatment

// Bottom left corner node
{
int X = 0;
int Y = 0;

rho[X][Y] = rho[X][Y+1]; // Extrapolation 1st order
velocity_x[X][Y] = (T)0;
velocity_y[X][Y] = (T)0;

f_old[X][Y][1] = f_old[X][Y][3] + (T)2/(T)3 * velocity_x[X][Y];
f_old[X][Y][2] = f_old[X][Y][4] + (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][5] = f_old[X][Y][7] + (T)1/(T)6 * (velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][6] = (T)1/(T)12 * (-velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][8] = (T)1/(T)12 * (velocity_x[X][Y]-velocity_y[X][Y]);
f_old[X][Y][0] = (T)0;

f_old[X][Y][0] = rho[X][Y] - (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);
}

// Top left corner node
{
int X = 0;
int Y = ly-1;

rho[X][Y] = rho[X][Y-1]; // Extrapolation 1st order
velocity_x[X][Y] = (T)0;
velocity_y[X][Y] = (T)0;

f_old[X][Y][1] = f_old[X][Y][3] + (T)2/(T)3 * velocity_x[X][Y];
f_old[X][Y][4] = f_old[X][Y][2] - (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][8] = f_old[X][Y][6] + (T)1/(T)6 * (velocity_x[X][Y]-velocity_y[X][Y]);
f_old[X][Y][5] = (T)1/(T)12 * (velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][7] = -(T)1/(T)12 * (velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][0] = (T)0;

f_old[X][Y][0] = rho[X][Y] - (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);
}

// Bottom right corner node
{
int X = lx-1;
int Y = 0;

rho[X][Y] = rho[X][Y+1]; // Extrapolation 1st order
velocity_x[X][Y] = (T)0;
velocity_y[X][Y] = (T)0;

f_old[X][Y][3] = f_old[X][Y][1] - (T)2/(T)3 * velocity_x[X][Y];
f_old[X][Y][2] = f_old[X][Y][4] + (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][6] = f_old[X][Y][8] + (T)1/(T)6 * (-velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][5] = (T)1/(T)12 * (velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][7] = -(T)1/(T)12 * (velocity_x[X][Y]-velocity_y[X][Y]);
f_old[X][Y][0] = (T)0;

f_old[X][Y][0] = rho[X][Y] - (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);

}

// Top right corner node
{
int X = lx-1;
int Y = ly-1;

rho[X][Y] = rho[X][Y-1]; // Extrapolation 1st order
velocity_x[X][Y] = (T)0;
velocity_y[X][Y] = (T)0;

f_old[X][Y][3] = f_old[X][Y][1] - (T)2/(T)3 * velocity_x[X][Y];
f_old[X][Y][4] = f_old[X][Y][2] - (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][7] = f_old[X][Y][5] - (T)1/(T)6 * (velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][6] = (T)1/(T)12 * (-velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][8] = (T)1/(T)12 * (velocity_x[X][Y]-velocity_y[X][Y]);
f_old[X][Y][0] = (T)0;

f_old[X][Y][0] = rho[X][Y] - (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);

}

// Swap popualations
// The present code used old and new populations which are swapped at the beginning of each time step.
// This way, the old populations are not overwritten during propagation.
// The resulting code is easier to write and debug.

T ***swap_tmp;

swap_tmp = f_old;
f_old = f;
f = swap_tmp;

// Compute fluid density and velocity
// This function computes the fluid density and velocity from the populations.

momenta();
}

// Compute fluid density and velocity

void momenta() {

for (int X = 0; X < lx; ++X) {

for (int Y = 0; Y < ly; ++Y) {

rho[X][Y] = f[X][Y][0] + f[X][Y][1] + f[X][Y][2] + f[X][Y][3] + f[X][Y][4] + f[X][Y][5] + f[X][Y][6] + f[X][Y][7] + f[X][Y][8];
velocity_x[X][Y] = (f[X][Y][1] + f[X][Y][5] + f[X][Y][8] - (f[X][Y][3] + f[X][Y][6] + f[X][Y][7])) /rho[X][Y];
velocity_y[X][Y] = (f[X][Y][2] + f[X][Y][5] + f[X][Y][6] - (f[X][Y][4] + f[X][Y][7] + f[X][Y][8])) /rho[X][Y];
}
}
}

You could define an interface function to access your array elements, calculating your array index’s based on an offset. something like:

index_x = (x-t*ex)%lx

where x is the coordinate of the node you are currently inspecting, t is the current iteration number, ex the lattice vector and lx the length. The modulus operation enforces periodicity.

An alternative approach would be a multi-dimensional implementation of STD::Rotate

edit: depending on the implementation of the modulo operator the expression may need to be ABS(index_x = (x-t*ex)%lx). It turns out mod in matlab is subtly different from % in c.

Can you tell me which type of Buonce back have you applied on the left and right wall, please?

Best regards,

Bolt,

Hi Bolt

Look this :

http://lbmworkshop.com/wp-content/.../Corners.pdf