# LBM D2Q9 Poiseuille

Hello all

I’m a beginner with the Lattice Boltzmann method. And I like to have your help, because I have any problems with the calculation of the pressure field in my rectangular duct for boundary conditions of Von Neumann (Velocity inlet/Outlet flow in side East and West):

CBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBC
I O
I O
I O
CBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBC

I: velocity Inlet flow
O: velocity Outlet flow
C: corner nodes
B: Bounce back

I would also like to know how to deal specifically with the corner nodes. Then I’d like to change my boundary conditions. At first, I’d imposed a velocity in inlet and pressure in outlet. Then I’d like an imposed pressure gradient between the inlet and outlet.

Here is my code, I know if I respect much LBM algorithm, but as I have said I am a novice with this method, I come through the forum, solicit your support, your comments, corrections and your explanation.

My code:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <unistd.h>

// Grid size in x and y dimensions
#define lx 301
#define ly 81

// Lattice Boltzmann D2Q9
#define Q 9

// D2Q9 weighting factors
double t0 = 4./9.;
double t1 = 1./9.;
double t2 = 1./36.;

#define density 0.1

// BGK collision time
double tau;

// Collision parameter
double omega;

// Maximum number of iteration
#define MaxStep 100000001

// Output frequency
#define NUMAX 1000

// Length of rectangular cylinder
#define a 21
#define b 21

// Distance from inlet boundary
#define x1 80

// Coordonate of rectangular point
int x2 = x1 + a;

// Computing control
int NUMB;

// Iteration counter
int istep;

// Velocity inlet
double Uw = 0.02;

// Reynolds number Re
double Re = 100.;

// Reference length for Reynolds number
double XYMAX;

// Density function of distribution
double f[lx][ly][Q];

// ftemp… array of temporare storage
double ftemp[lx][ly][Q];

// Solid wall nodes array
int wall[lx][ly];

//
/* Define solid nodes wall boundary */
/
/

void Wall_coordinate(void)

{
int x,y,obst_x,obst_y,obst_r;

obst_x = (lx-1)/5;
obst_y = (ly-1)/2;
obst_r = (ly-1)/10;

for (y = 0; y < ly; y++) {

for (x = 0; x < lx; x++) {

wall[x][y] = 0;
}
}

for (x = 0; x < lx; x++) {

wall[x][0] = 1;
wall[x][ly-1] = 1;
}

for (x = 0; x < lx; x++) {

for (y = 0; y < ly; y++) {

if (((x - obst_x) * (x - obst_x) + (y - obst_y) * (y - obst_y)) <= (obst_r * obst_r)) wall[x][y] = 1;
}
}

/*for (x = x1; x < x2; x++) {

wall[x][y1] = 1;
wall[x][y2] = 1;
}*/

/*for (y = y1; y < y2; y++) {

wall[x1][y] = 1;
wall[x2][y] = 1;
}*/

}

/********************************/
/
Initialize density function /
/
********************************/

void Init_density(void)

{
int x,y;
double w0,w1,w2;

// Compute weighting factors

w0 = t0 * density;
w1 = t1 * density;
w2 = t2 * density;

// Loop over computational domain

for (x = 0; x < lx; x++) {

for (y = 0; y < ly; y++) {

f[x][y][0] = w0;
f[x][y][1] = w1;
f[x][y][2] = w1;
f[x][y][3] = w1;
f[x][y][4] = w1;
f[x][y][5] = w2;
f[x][y][6] = w2;
f[x][y][7] = w2;
f[x][y][8] = w2;
}
}

}

//
/* Inlet & outlet velocity bounday condition */
/
/

void Inlet_outlet_BC(void)

{

//local variables
int i,x,y;
double C1,C2,C3,U,V,uv[Q],Dxy;

//square speed of sound
C1 = 1./3.;
C2 = 2. * C1 * C1;
C3 = 2. * C1;

// Inlet velocity Bondary Condition

for (y = 0; y < ly; y++) {

x = 0;

Dxy = density;

// give velocity on the boundary
U = Uw;
V = 0.0;

// square velocity
uv[0] = U * U + V * V;
uv[1] = U;
uv[2] = V;
uv[3] = - U;
uv[4] = - V;
uv[5] = U + V;
uv[6] = - U + V;
uv[7] = - U - V;
uv[8] = U - V;

f[x][y][0] = t0 * Dxy * (1. - uv[0]/C3);
f[x][y][1] = t1 * Dxy * (1. + uv[1]/C1 + uv[1] * uv[1]/C2 - uv[0]/C3);
f[x][y][2] = t1 * Dxy * (1. + uv[2]/C1 + uv[2] * uv[2]/C2 - uv[0]/C3);
f[x][y][3] = t1 * Dxy * (1. + uv[3]/C1 + uv[3] * uv[3]/C2 - uv[0]/C3);
f[x][y][4] = t1 * Dxy * (1. + uv[4]/C1 + uv[4] * uv[4]/C2 - uv[0]/C3);
f[x][y][5] = t2 * Dxy * (1. + uv[5]/C1 + uv[5] * uv[5]/C2 - uv[0]/C3);
f[x][y][6] = t2 * Dxy * (1. + uv[6]/C1 + uv[6] * uv[6]/C2 - uv[0]/C3);
f[x][y][7] = t2 * Dxy * (1. + uv[7]/C1 + uv[7] * uv[7]/C2 - uv[0]/C3);
f[x][y][8] = t2 * Dxy * (1. + uv[8]/C1 + uv[8] * uv[8]/C2 - uv[0]/C3);
}

// Outlet velocity Bondary Condition

for (y = 0; y < ly; y++) {

x = lx-1;

Dxy = density;

// give velocity on the boundary
U = Uw;
V = 0.0;

// square velocity
uv[0] = U * U + V * V;
uv[1] = U;
uv[2] = V;
uv[3] = - U;
uv[4] = - V;
uv[5] = U + V;
uv[6] = - U + V;
uv[7] = - U - V;
uv[8] = U - V;

f[x][y][0] = t0 * Dxy * (1. - uv[0]/C3);
f[x][y][1] = t1 * Dxy * (1. + uv[1]/C1 + uv[1] * uv[1]/C2 - uv[0]/C3);
f[x][y][2] = t1 * Dxy * (1. + uv[2]/C1 + uv[2] * uv[2]/C2 - uv[0]/C3);
f[x][y][3] = t1 * Dxy * (1. + uv[3]/C1 + uv[3] * uv[3]/C2 - uv[0]/C3);
f[x][y][4] = t1 * Dxy * (1. + uv[4]/C1 + uv[4] * uv[4]/C2 - uv[0]/C3);
f[x][y][5] = t2 * Dxy * (1. + uv[5]/C1 + uv[5] * uv[5]/C2 - uv[0]/C3);
f[x][y][6] = t2 * Dxy * (1. + uv[6]/C1 + uv[6] * uv[6]/C2 - uv[0]/C3);
f[x][y][7] = t2 * Dxy * (1. + uv[7]/C1 + uv[7] * uv[7]/C2 - uv[0]/C3);
f[x][y][8] = t2 * Dxy * (1. + uv[8]/C1 + uv[8] * uv[8]/C2 - uv[0]/C3);
}

}

//
/
Propagate fluid densities to their next neighbour nodes /
/
/

void Stream(void)

{

//local variables
int x,y,xn,xp,yn,yp;

//loop over all f

for (x = 0; x < lx; x++) {

xn = (x > 0)?(x-1):(lx-1);
xp = (x < lx-1)?(x+1):(0);

for (y = 0; y < ly; y++) {

yn = (y > 0)?(y-1):(ly-1);
yp = (y < ly-1)?(y+1):(0);

ftemp[x][y][0] = f[x][y][0];
ftemp[xp][y][1] = f[x][y][1];
ftemp[x][yp][2] = f[x][y][2];
ftemp[xn][y][3] = f[x][y][3];
ftemp[x][yn][4] = f[x][y][4];
ftemp[xp][yp][5] = f[x][y][5];
ftemp[xn][yp][6] = f[x][y][6];
ftemp[xn][yn][7] = f[x][y][7];
ftemp[xp][yn][8] = f[x][y][8];
}
}

}

//
/
Fluid densities are rotated. By the next propagation step /
/
this results in a bounce back from obstacle nodes. /
/
/

void Bounceback(void)

{

//local variables
int x,y;

//loop over all f

for (x = 0; x < lx; x++) {

for (y = 0; y < ly; y++) {

if (wall[x][y] == 1) {

f[x][y][1] = ftemp[x][y][3];
f[x][y][2] = ftemp[x][y][4];
f[x][y][3] = ftemp[x][y][1];
f[x][y][4] = ftemp[x][y][2];
f[x][y][5] = ftemp[x][y][7];
f[x][y][6] = ftemp[x][y][8];
f[x][y][7] = ftemp[x][y][5];
f[x][y][8] = ftemp[x][y][6];
}
}
}
}

/**************************************/
/
One-step density Collision process /
/
**************************************/

void Collision(void)

{
//local variables
int i,x,y;
double C1,C2,C3,U,V,uv[Q],Dxy,feq[Q];

//square speed of sound
C1 = 1./3.;
C2 = 2. * C1 * C1;
C3 = 2. * C1;

for (x = 0; x < lx; x++) {

for (y = 0; y < ly; y++) {

Dxy = 0.;

//only fluid nodes are considered here

if (wall[x][y] == 0) {

// calculate local density Dxy
for (i = 0; i < Q; i++) { Dxy += ftemp[x][y]; }

// calculate velocity components
U = (ftemp[x][y][1] + ftemp[x][y][5] + ftemp[x][y][8] - (ftemp[x][y][3] + ftemp[x][y][6] + ftemp[x][y][7])) / Dxy;
V = (ftemp[x][y][2] + ftemp[x][y][5] + ftemp[x][y][6] - (ftemp[x][y][4] + ftemp[x][y][7] + ftemp[x][y][8])) / Dxy;

// square velocity
uv[0] = U * U + V * V;
uv[1] = U;
uv[2] = V;
uv[3] = - U;
uv[4] = - V;
uv[5] = U + V;
uv[6] = - U + V;
uv[7] = - U - V;
uv[8] = U - V;

feq[0] = t0 * Dxy * (1. - uv[0]/C3);
feq[1] = t1 * Dxy * (1. + uv[1]/C1 + uv[1] * uv[1]/C2 - uv[0]/C3);
feq[2] = t1 * Dxy * (1. + uv[2]/C1 + uv[2] * uv[2]/C2 - uv[0]/C3);
feq[3] = t1 * Dxy * (1. + uv[3]/C1 + uv[3] * uv[3]/C2 - uv[0]/C3);
feq[4] = t1 * Dxy * (1. + uv[4]/C1 + uv[4] * uv[4]/C2 - uv[0]/C3);
feq[5] = t2 * Dxy * (1. + uv[5]/C1 + uv[5] * uv[5]/C2 - uv[0]/C3);
feq[6] = t2 * Dxy * (1. + uv[6]/C1 + uv[6] * uv[6]/C2 - uv[0]/C3);
feq[7] = t2 * Dxy * (1. + uv[7]/C1 + uv[7] * uv[7]/C2 - uv[0]/C3);
feq[8] = t2 * Dxy * (1. + uv[8]/C1 + uv[8] * uv[8]/C2 - uv[0]/C3);

for (i = 0; i < Q; i++) {f[x][y] = (1. - omega) * ftemp[x][y] + omega * feq;}

}
}
}

}

/********************************************************************************************/
/
Compute total density to verify of the system that no divergence occurs when iterating /
/
********************************************************************************************/

void check_density(void)

{
//local variables
int x,y,iLB;
double sum = 0.;

for (x = 0; x < lx; x++)

{

for (y = 0; y < ly; y++)

{
for (iLB = 0; iLB < Q; iLB++){ sum +=f[x][y][iLB];}
}
}

printf("–> Total density in the system %lf\n",sum);
}

//
/* Output of results to file ’ densities000000.vtyk ’ */
/
/

void Write_densities(void)

{
//local variables
int i,x,y;
double U,V,P,W,Dxy,Dxy0,cs,pasxyz;

//square speed of sound
cs = 1./3.;

pasxyz = 1./lx;

printf(“istep = %d\n”,istep);

char nomfic[25];

FILE * sortie;
sprintf(nomfic,“Visual/densities%.6i.vtk”,istep);
sortie = fopen(nomfic, “w”);

fprintf(sortie,"# vtk DataFile Version 2.0\n");
fprintf(sortie,“ASCII\n”);
fprintf(sortie,“DATASET RECTILINEAR_GRID\n”);
fprintf(sortie,“DIMENSIONS %d %d 1\n”,lx,ly);
fprintf(sortie,“X_COORDINATES %d float\n”,lx);

for (i = 0; i < lx; i++) { fprintf(sortie,"%e ",(float)i * pasxyz);}

fprintf(sortie,"\n");
fprintf(sortie,“Y_COORDINATES %d float\n”,ly);

for (i = 0; i < ly; i++) { fprintf(sortie,"%e ",(float)i * pasxyz);}

fprintf(sortie,"\n");
fprintf(sortie,“Z_COORDINATES 1 float\n”);
fprintf(sortie,“0\n”);

fprintf(sortie,“POINT_DATA %d\n”,lx*ly);
fprintf(sortie,“SCALARS Pression float 1\n”);
fprintf(sortie,“LOOKUP_TABLE default\n”);

for (y = 0; y < ly; y++) {

for (x = 0; x < lx; x++) {

Dxy = 0.;

// pressure = average pressure
if (wall[x][y] == 1) { P = 0.;}

else {

// calculate local Density Dxy
for (i = 0; i < Q; i++) { Dxy += f[x][y];}

// calculate pressure
P = Dxy * cs;
}

fprintf(sortie,"%.4lf\n",P);
}
}

fprintf(sortie,“VECTORS VecVelocity float\n”);

for (y = 0; y < ly; y++)

{
for (x = 0; x < lx; x++)

{
Dxy = 0.;

if (wall[x][y] == 1)

{
U = 0.;
V = 0.;
}

else

{
// calculate local Density Dxy
for (i = 0; i < Q; i++) { Dxy += f[x][y];}

// calculate velocity components

U = (f[x][y][1] + f[x][y][5] + f[x][y][8] - (f[x][y][3] + f[x][y][6] + f[x][y][7]))/Dxy;
V = (f[x][y][2] + f[x][y][5] + f[x][y][6] - (f[x][y][4] + f[x][y][7] + f[x][y][8]))/Dxy;
}

fprintf(sortie,"%.4lf %.4lf 0.\n",U,V);
}
}

fclose(sortie);

}

// Principal function main()

int main(int argc, char** argv)

{

system(“mkdir Visual”);

// Begin Initialisation

// Call wall coordinate function
Wall_coordinate();

// Call init density function
Init_density();

XYMAX = x2 - x1;

tau = Uw * XYMAX * 3./Re + 1./2.;
omega = 1./tau;
Re = Uw * XYMAX * 3./(1./omega - 1./2.);

printf("\n");
printf(“BGK collision time\n”);
printf(“tau = %.4lf\n”,tau);
printf("\n");
printf(“Reynolds number\n”);
printf(“Re = %.2lf\n”,Re);
printf("\n");
printf(“Grid size dimensions\n”);
printf(“L = %d , H = %d\n”,lx,ly);
printf("\n");

// Initial call check densities
check_density();

// End Initialisation !!!

// Begin Iterations

NUMB = 0;

for (istep = 0; istep < MaxStep; istep++) {

if (NUMB <= NUMAX) NUMB += 1;

// Call Inlet-outlet boundary condition
Inlet_outlet_BC();

// Call Stream density propagation
Stream();

// Call bounc back from solid wall
Bounceback();

// Call density Collision
Collision();

/* each MaxStep/10 iteration … */

if (NUMB > NUMAX) {

/* Output the data velocity, vorticity and pressure */

// Call write output files vtk
Write_densities();

// Call check densities
check_density();

NUMB = 1;
}

}

// End Iterations !!!

return 0;

}

thank you very much !!!