Good morning !!
Please,
I wanted to coupled 2D DEM and LBM method to simulate fluid and moving obstacles intercation. i want to impose LBM boundaries conditions such as : in South side (Inlet) velocity and on North side (Outlet) pressure boundary. I proceed in this way but i get errors:
- compute fluid macroscopics variables (density and velocities)
- collision step
- bounce back for moving particles (obstacles)
- bounce back for walls
- streaming step
Define functions :
void computeMacros() {
int iX, iY, iPop;
realLBM rho, u_x, u_y;
for (iX = 0; iX < lx; iX++) {
for (iY = 0; iY < ly; iY++) {
rho = 0.;
u_x = 0.;
u_y = 0.;
lattice[iX][iY].rho = 0.;
lattice[iX][iY].ux = 0.;
lattice[iX][iY].uy = 0.;
if (obst[iX][iY] == FLUID) {
for (iPop = 0; iPop < Q; iPop++) {
rho += lattice[iX][iY].f[iPop];
u_x += cx[iPop] * lattice[iX][iY].f[iPop];
u_y += cy[iPop] * lattice[iX][iY].f[iPop];
}
lattice[iX][iY].rho = rho;
u_x /= rho;
u_y /= rho;
lattice[iX][iY].ux = u_x;
lattice[iX][iY].uy = u_y;
}
}
}
}
void Inlet_Velocity() {
int iX, iY, iPop;
for (iX = 1; iX < lx-1; iX++) {
iY = 0;
if (obst[iX][iY] == FLUID) {
lattice[iX][iY].ux = 0.0;
lattice[iX][iY].uy = vyMax0;
lattice[iX][iY].rho = (lattice[iX][iY].f[0] + lattice[iX][iY].f[2] + lattice[iX][iY].f[6]
+ 2. * (lattice[iX][iY].f[3] + lattice[iX][iY].f[4] + lattice[iX][iY].f[5])) / (1. - vyMax0);
lattice[iX][iY].f[8] = lattice[iX][iY].f[4] + (2./3.) * lattice[iX][iY].uy;
lattice[iX][iY].f[1] = lattice[iX][iY].f[5] + (1./6.) * lattice[iX][iY].uy - (1./2.) * (lattice[iX][iY].f[2] - lattice[iX][iY].f[6]);
lattice[iX][iY].f[7] = lattice[iX][iY].f[3] + (1./6.) * lattice[iX][iY].uy - (1./2.) * (lattice[iX][iY].f[6] - lattice[iX][iY].f[2]);
} }
}
void Outlet_Pressure() {
int iX, iY;
for (iX = 1; iX < lx-1; iX++) {
iY = ly-1;
if (obst[iX][iY] == FLUID) {
lattice[iX][iY].rho = rho_out;
lattice[iX][iY].ux = 0.0;
lattice[iX][iY].uy = (lattice[iX][iY].f[0] + lattice[iX][iY].f[2] + lattice[iX][iY].f[6]
+ 2. * (lattice[iX][iY].f[1] + lattice[iX][iY].f[7] + lattice[iX][iY].f[8])) - rho_out;
lattice[iX][iY].f[4] = lattice[iX][iY].f[8] - (2./3.) * lattice[iX][iY].uy;
lattice[iX][iY].f[5] = lattice[iX][iY].f[1] - (1./6.) * lattice[iX][iY].uy + (1./2.) * (lattice[iX][iY].f[2] - lattice[iX][iY].f[6]);
lattice[iX][iY].f[3] = lattice[iX][iY].f[7] - (1./6.) * lattice[iX][iY].uy + (1./2.) * (lattice[iX][iY].f[6] - lattice[iX][iY].f[2]);
}
}
}
void Collision_BGK() {
int iPop, iX, iY;
double tauInv = 1./tau;
// Compute density and velocity from the f’s
computeMacros();
rho = 0.;
ux = 0.;
uy = 0.;
realLBM uSqr = 0.;
realLBM feq = 0.;
realLBM fi = 0.;
for (iX = 0; iX < lx; iX++) {
for (iY = 0; iY < ly; iY++) {
rho = lattice[iX][iY].rho;
ux = lattice[iX][iY].ux;
uy = lattice[iX][iY].uy;
if (obst[iX][iY] == FLUID) {
uSqr = ux * ux + uy * uy;
feq = computeEquilibrium(iPop, rho, ux, uy, uSqr);
for (iPop = 0; iPop < Q; iPop++) {
fi = (1. - tauInv) * lattice[iX][iY].f[iPop] + tauInv * feq;
lattice[iX][iY].f[iPop] = fi;
}
}
}
}
}
void swapAfterCollision() {
int iX, iY, iPop;
static const int half = (Q-1)/2;
for (iX = 0; iX < lx; iX++) {
for (iY = 0; iY < ly; iY++) {
for (iPop = 1; iPop <= half; iPop++) {swap(&lattice[iX][iY].f[iPop], &lattice[iX][iY].f[iPop+half]);}
}
}
}
void Streaming() {
int iX, iY, iPop, next_x, next_y;
static const int half = (Q-1)/2;
for (iX = 0; iX < lx; iX++) {
for (iY = 0; iY < ly; iY++) {
for (iPop = 1; iPop <= half; iPop++) {
next_x = iX + cx[iPop];
next_y = iY + cy[iPop];
if (next_x >= 0 && next_y >= 0 && next_x < lx && next_y < ly) {swap(&lattice[iX][iY].f[iPop+half], &lattice[next_x][next_y].f[iPop]);}
}
}
}
}
void BounceBack() {
register int iX, iY, iPop, next_x, next_y;
for (iX = 0; iX < lx; iX++) {
for (iY = 0; iY < ly; iY++) {
if (obst[iX][iY] == WALL) {
for (iPop = 1; iPop < Q; iPop++) {
next_x = iX + cx[iPop];
next_y = iY + cy[iPop];
if (obst[next_x][next_y] == FLUID) {
lattice[iX][iY].f[iPop] = lattice[next_x][next_y].f[oppositeOf[iPop]];
}
}
}
}
}
}
void CollisionAndStreaming() {
// Post-collision part computation
Collision_MRT();
#ifdef ACTIVE_GRAINS
// Apply Bounce Back on moving particles
#endif
// Apply Bounce back on walls (No-slip)
BounceBack();
swapAfterCollision();
Streaming();
// Apply Boundaries Conditions
Inlet_Velocity();
Outlet_Pressure();
}
if you need more information I will give you.
thank you for your help !!!