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)6velocity_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)6velocity_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)6velocity_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)6velocity_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];
}
}
}