Inlet/Outlet boundary conditions

Thanks for your effort to make this open-source Lattice-Boltzmann solver (tu)

I am trying to model the flow in a bio-reactor using openLB. The geometry of the boundary is non-trivial and I
think I have to use bounce-back condition there. I tried using uniform flow profile at the inlet and outlet and it runs but the average energy steadily decreases. I suspect a loss of mass. I am pretty inexperience with LB and was wondering if you could advise me on the best choice of boundary conditions at the inlet and outlet.

If you are imposing the velocity profile on the inlet and outlet, the simulation ought to be pretty well behaving. It is possible that something went wrong when you define the boundary condition. Have you had a look at the example code poiseuille.cpp? The strategy adopted in this example to impose a velocity profile should also work in your case. If the problem persists, consider posting a piece of code to show how you implement your inlet/outlet.

I have just tried imposing parabolic velocity profile at the inlet/outlet but still I get a constant decrease of the energy indicating a mass loss to me!

In the Poiseuille case, proper velocity boundary conditions are imposed on the wall whereas I have to impose bounce-back. Should that make anydifference?
Thanks

If I understand correctly what you are doing you are defining the inlet/outlet velocity on bounce-back nodes. Is that correct?

If it is then this will have no real effect on the dynamics (of course it depends on the rest of the initialization of your system) of your simulation since the only thing done is to store the values of the velocity and the density on these nodes, but it cannot influence of your system, since fullway bounce-back can only impose zero velocity at the boundaries (you can test this case in the poiseuille example by replacing lines 71 and 72 by

lattice.defineDynamics(0,0, 1,ny-2, &instances::getBounceBack<T,DESCRIPTOR>());
lattice.defineDynamics(nx-1,nx-1, 1,ny-2, &instances::getBounceBack<T,DESCRIPTOR>());

You should rather use one of the available velocity/pressure boundary conditions (LocalBoundaryCondition, InterpBoundaryCondition, ZouHeBoundaryCondition, or any other). If you have a too complicated geometry then maybe the simplest would be to define a pressure gradient between inlet and outlet.

If not I will also recommend you to post the piece of code where you define your inlet and outlet.

That your energy is steadily decreasing does not necessarily indicate a loss of mass, but rather that your inlet/outlet condition does not work. It is OK to use bounce-back on the side walls. But bounce-back doesn’t make sense on the inlet/outlet, of course (at least not the full-way bounce-back implemented in OpenLB ). For these, use one of the standard velocity boundary conditions in OpenLB, and your program should work without problem. Actually, the simulation setup you are describing is straightforward.and should be fixed easily.

Hello,
Here’s my function which defines the geometry!
Sorry it’s long but it’s essentially geomtry definition.
Any help would be greatly appreciated.
Mathieu

const T x1 = 0.;
const T y1l = -4.11e-4;
const T y1u = 4.11e-4;
const T x2 = 1.41e-3;
const T y2l = -4.11e-4;
const T y2u = 4.11e-4;
const T x3 = 4.e-3;
const T y3l = -2.6e-3;
const T y3u = 2.6e-3;
const T x4 = 2e-2;
const T y4l = -4.11e-4;
const T y4u = 4.11e-4;
const T L = 1.91e-2;
const T L0 = 8.22e-4; // hydraulic diameter

T poiseuilleVelocity(int iY, int N, T u) {
T y = (T)iY / (T)N;
return 4.u * (y-yy);
}

void iniGeometry( BlockStructure2D<T,DESCRIPTOR>& lattice,
LBunits const& converter,
Dynamics<T, DESCRIPTOR>& bulkDynamics,
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )
{
const T omega = converter.getOmega();

const T x0 = 0.;
const T y0 = -3.2e-3;

const T y_or = -y0/L0;
const T x1_s = (x1-x0)/L0;
const T y1l_s = (y1l-y0)/L0;
const T y1u_s = (y1u-y0)/L0;
const T x2_s = (x2-x0)/L0;
const T y2l_s = (y2l-y0)/L0;
const T y2u_s = (y2u-y0)/L0;
const T x3_s = (x3-x0)/L0;
const T y3l_s = (y3l-y0)/L0;
const T y3u_s = (y3u-y0)/L0;
const T x4_s = (x4-x0)/L0;
const T y4l_s = (y4l-y0)/L0;
const T y4u_s = (y4u-y0)/L0;
const T L_s = L/L0;

const int nx = converter.getNx();
const int ny = converter.getNy();
const int n_x0 = converter.nCell(x0);
const int n_y0 = converter.nCell(y_or);
const int n_x1 = converter.nCell(x1_s);
const int n_y1l = converter.nCell(y1l_s);
const int n_y1u = converter.nCell(y1u_s);
const int n_x2 = converter.nCell(x2_s);
const int n_y2l = converter.nCell(y2l_s);
const int n_y2u = converter.nCell(y2u_s);
const int n_x3 = converter.nCell(x3_s);
const int n_y3l = converter.nCell(y3l_s);
const int n_y3u = converter.nCell(y3u_s);
const int n_x4 = converter.nCell(x4_s);
const int n_y4l = converter.nCell(y4l_s);
const int n_y4u = converter.nCell(y4u_s);
const int n_L = converter.nCell(L_s);

T *js, *jn;
js = new T[nx];
jn = new T[nx];

// Bulk dynamics
lattice.defineDynamics(0,nx-1,0,ny-1,    &bulkDynamics);

// Definition of the obstacle (bounce-back nodes)

int iY, nYmin, nYmax, iterate, iswitch=0, isw, inw, ise, ine;

T u[2] = {converter.getLatticeU(),0.};
T rho = (T)1;

cout << "n_y4l = " << n_y4l << "\n";
cout << "n_y4u = " << n_y4u << "\n";
cout << "n_y0 = " << n_y0 << "\n";
cout << "n_L = " << n_L << "\n";
cout << "n_x4 = " << n_x4 << "\n";

// Find the indices of upper and lower nodes on the wall
for (int iX=0; iX<nx; ++iX) {

	iY = 0;
iterate = 1;

while (iterate==1)
{
	if (iX<n_x2)
	{
		nYmin = n_y1l;
		nYmax = n_y1u;
		if ((iY>nYmin)&&(iY<nYmax))
		{
			if (iswitch==0)
			{
				iswitch=1;
				js[iX] = iY; 
				if (iX==0)
				{
					isw = iY;
				}
			}
		}
		else
		{	
			if (iswitch==1)
			{
				iswitch=0;
				jn[iX] = iY-1;
				if (iX==0)
				{
					inw = iY-1;
				}
			}
		}
	}
	else if ((iX>=n_x2)&&(iX<n_x3))
	{
		nYmin = (n_y3l-n_y2l)*(iX-n_x2)/(n_x3-n_x2)+n_y2l;
		nYmax = (n_y3u-n_y2u)*(iX-n_x2)/(n_x3-n_x2)+n_y2u;
		if ((iY>nYmin)&&(iY<nYmax))
		{
			if (iswitch==0)
			{
				iswitch=1;
				js[iX] = iY; 
			}
		}
		else
		{
			if (iswitch==1)
			{
				iswitch=0;
				jn[iX] = iY-1;
			}
		}
	}
	else
	{
		nYmin = (n_y4l-n_y0)*n_L/(n_L-(n_x4-iX))+n_y0;
		nYmax = (n_y4u-n_y0)*n_L/(n_L-(n_x4-iX))+n_y0;
		if ((iY>nYmin)&&(iY<nYmax))
		{
			if (iswitch==0)
			{
				iswitch=1;
				js[iX] = iY; 
				if (iX==(nx-1))
				{
					ise = iY;
				}
			}
		}
		else
		{
			
			if (iswitch==1)
			{
				iswitch=0;
				jn[iX] = iY-1;
				if (iX==(nx-1))
				{
					ine = iY-1;
				}
			}
		}
	}
	if (iY>=(ny-1)) 
		iterate = 0;
	else
		++iY; 
}
}

// Define parabolic profile in the fluid
// Define bounce-back at the boundary
for (int iX=0; iX<nx; ++iX) {

	iY = 0;
iterate = 1;

while (iterate==1)
{
	if ((iY>=js[iX])&&(iY<=jn[iX]))
	{
		T vel[] = {
            	poiseuilleVelocity(iY-js[iX], jn[iX]-js[iX], converter.getLatticeU()),
            	T()
        		};
        		lattice.get(iX,iY).defineRhoU((T)1, vel);
        		lattice.get(iX,iY).iniEquilibrium((T)1, vel);;
	}
	else
	{
		lattice.defineDynamics( iX,iX,iY,iY,
                    &instances::getBounceBack<T,DESCRIPTOR>() );
	}
	if (iY>=(ny-1)) 
		iterate = 0;
	else
		++iY; 
}
}

cout << "isw = " << isw << "\n";
cout << "inw = " << inw << "\n"; 
cout << "ise = " << ise << "\n";
cout << "ine = " << ine << "\n";

boundaryCondition.addVelocityBoundary0N(0,0,       isw,inw, omega);
boundaryCondition.addVelocityBoundary0P(nx-1,nx-1,  ise, ine, omega);

// Make the lattice ready for simulation
lattice.initialize();

delete [] js;
delete [] jn;

}

The lines

boundaryCondition.addVelocityBoundary0N(0,0, isw,inw, omega);
boundaryCondition.addVelocityBoundary0P(nx-1,nx-1, ise, ine, omega);

should be written before the imposition of the parabolic profile, because in a normal bulk node, you do not store the velocity and/or density. Therefore they are by default initialized (and therefore kept) as being zero (for the velocity) and one (for the density).

Thank you very much. The results make a lot more sense now!
Regards,

hi
i want to use pressure outlet boundary condition
can u help me?