# Coupled 2D DEM-LBM method

Good morning !!

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();

}