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.;
// Fluid density per link
#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,“Sortie domaine LB+LINK t %d\n”,istep);
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 !!!