Streaming and Collision (SWAP algorithm)

Lattice Boltzmann Method LBM: Algorithms
#1

Hello,

I am trying to simulate Poisuille flow by LBM. I coded Streaming and Collision by SWAP algorithm but getting the results (X-velocity, Y-velocity, and Velocity Magnitude) like that (attached as picture).

Velocity-Y
Velocity-X
Velocity-Magnitude

[code=“cpp”]
int fStreamingSwap() {
const int half = (Q - 1) / 2;

for (int ix = 1; ix <= Nx - 1; ++ix) {
	for (int iy = 1; iy <= Ny - 1; ++iy) {
		for (int iq = 1; iq <= half; ++iq) {
			int nextX = ix - (int)c[iq][0];
			int nextY = iy - (int)c[iq][1];
			if (nextX >= 1 && nextY >= 1 && nextX <= Nx-1 && nextY <= Ny-1 ) {
				std::swap(fpop[ix][iy][iq+half], fpop[nextX][nextY][iq]);
			}
		}
	}
}

return 0;

}



[code="cpp"]
int fBgkCollisionSwap() {
	const int half = (Q - 1) / 2;
	fMacroscopicVariables();
	fComputeEquilibrium();
	for (int ix = 1; ix <= Nx - 1; ++ix) {
		for (int iy = 1; iy <= Ny - 1; ++iy) {
			for (int iq = 0; iq < Q; ++iq) {
				fpop[ix][iy][iq] *= (1.0 - omega);
				fpop[ix][iy][iq] += omega*fEq[ix][iy][iq];
			}
			for (int iq = 1; iq <= half; ++iq) {
				std::swap(fpop[ix][iy][iq], fpop[ix][iy][iq+half]);
			}
		}
	}

	return 0;
}

Kindly suggest what happened.
Best Regards

Ali

#2

Hi Ali,

If the numbering of your links is correct (conjugates are (Q-1)/2 links apart, the first four links point in negative x and y directions), the only thing I can think of is that you need to make sure the boundary points at ix = 0, ix = Nx, iy = 0 and iy = Ny are all dealt with correctly.

My suspicion is that while you may have dealt with boundary conditions at iy = 0 and iy = Ny correctly, I’m not sure you have been copying distribution functions from ix = Nx-1 to ix = 0 and from ix = 1 to ix = Nx to give the correct periodic boundary conditions. If these boundary points are essentially empty, you will get that striping effect using this form of the swap propagation algorithm. (You might want to check what the total mass of the fluid is at each time step: it should be constant but I think it might be decreasing gradually over time.)

Regards,

Michael

#3

Thanks Michael, can u share ur email so I will sen you my code to look into it?

Best Regards
Ali

#4

Total Mass is 0.0000 in my case. Will it be Zero or One?

This is how I am Computing total Mass:

[code=“cpp”]
int MassConservation() {

double mass = 0.0;

for (int ix = 1; ix <= Nx-1; ++ix) {
	for (int iy=1; iy <= Ny-1; iy++) {
		mass += rho[ix][iy];
	}
}

return 0;

}



How to apply Boundary conditions, my solution domain is:

[code="cpp"]
int fInitializeSolutionDomain() {
// Virtual Nodes

	for (int ix = 0; ix <= Nx; ++ix) {
		nodes[ix][0]  = -1.0;
		nodes[ix][Ny] = -1.0;
	}

	for (int iy = 0; iy <= Ny; ++iy) {
		nodes[0][iy]  = -1.0;
		nodes[Nx][iy] = -1.0;
	}
	std::cout << "Virtual nodes initialized..." << std::endl;
	
	// for inlet/outlet  // periodic b/c
	for (int iy = 2; iy <= Ny-2; ++iy) {
		nodes[1][iy]    = 0.1;
		nodes[Nx-1][iy] = 0.5;
	}
	std::cout << "Inlet nodes initialized..." << std::endl;
// fluid domain
	for (int ix = 1; ix <= Nx-1; ++ix) {
		for (int iy = 1; iy <= Ny-1; ++iy) {
			nodes[ix][iy] = 1.0;
		}
	}
	std::cout << "Fluid domains initialized..." << std::endl;

// for no slip
	for (int ix = 1; ix <= Nx - 1; ++ix) {
		nodes[ix][1] = 0.0;
		nodes[ix][Ny-1] = 0.0;
	}
	std::cout << "No-slip nodes initialized..." << std::endl;
	return 0;
}

and the lattice constants are defined as:

[code=“cpp”]
static const double lbw[9] = { 4./9., 1./36., 1./9., 1./36., 1./9.,
1./36., 1./9., 1./36., 1./9. };
static const double c[9][2] = {
{ 0., 0. },
{ -1., 1. }, { -1., 0. }, { -1., -1. }, { 0., -1. },
{ 1., -1. }, { 1., 0. }, { 1., 1. }, { 0., 1. }
};



Upper and Lower have no-Slip boundary condition whereas West =Inlet and East=Outlet.

Kindly comment and suggest boundary condition implementation.

Regards Ali
#5

Hi Ali,

The total system mass should not be zero! If you started with a density of 1 at each grid point, the total mass of the system should (in theory) be equal to 1 * the number of fluid grid points in the system and remain at that value for the kind of simulation you are carrying out.

The routines for calculating system mass and the lattice constants seem sensible, but I’m not sure what is happening with the boundary conditions. (What does the array “nodes” do in the code?)

Regards,

Michael

#6

The array “node” is defining the solution domain, actually. node = 1.0, is for fluid, 0.0 for no-slip BC. Did you see my “mass conservation function” in my last post, Is there any modification needed?

Regarding BCs"

West Wall = Inlet (Velocity profile or uniform flow)
East Wall = Outlet (Gradient BC)
Top and bottom walls = (either no-slip or free-slip BCs)

Best Regards
Ali